LifeV
AssemblyElemental.cpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  @file
29  @brief File containing the procedures for the local assembly of the differential operators
30 
31  @contributor Samuel Quinodoz <samuel.quinodoz@epfl.ch>
32  @mantainer Samuel Quinodoz <samuel.quinodoz@epfl.ch>
33  */
34 
35 #ifndef ELEMOPER_CPP
36 #define ELEMOPER_CPP 1
37 
38 #include <lifev/core/fem/AssemblyElemental.hpp>
39 #include <boost/multi_array.hpp>
40 
41 namespace LifeV
42 {
43 
44 namespace AssemblyElemental
45 {
46 void mass (MatrixElemental& localMass,
47  const CurrentFE& massCFE,
48  const Real& coefficient,
49  const UInt& fieldDim)
50 {
51  const UInt nbFEDof (massCFE.nbFEDof() );
52  const UInt nbQuadPt (massCFE.nbQuadPt() );
53  Real localValue (0);
54 
55  MatrixElemental::matrix_type mat_tmp (nbFEDof, nbFEDof);
56 
57  // Loop over the basis functions (diagonal part)
58  for (UInt iDof (0); iDof < nbFEDof; ++iDof)
59  {
60  localValue = 0.0;
61  for (UInt iQuadPt (0); iQuadPt < nbQuadPt; ++iQuadPt)
62  {
63  localValue += massCFE.phi (iDof, iQuadPt)
64  * massCFE.phi (iDof, iQuadPt)
65  * massCFE.wDetJacobian (iQuadPt);
66  }
67  localValue *= coefficient;
68 
69  // Add on the local matrix
70  mat_tmp (iDof, iDof) = localValue;
71  }
72 
73  // Loop over the basis functions (extradiagonal part)
74  for (UInt iDof (0); iDof < nbFEDof; ++iDof)
75  {
76  // Build the local matrix only where needed:
77  // Lower triangular + diagonal parts
78  for (UInt jDof (0); jDof < iDof; ++jDof)
79  {
80  localValue = 0.0;
81 
82  //Loop on the quadrature nodes
83  for (UInt iQuadPt (0); iQuadPt < nbQuadPt; ++iQuadPt)
84  {
85  localValue += massCFE.phi (iDof, iQuadPt)
86  * massCFE.phi (jDof, iQuadPt)
87  * massCFE.wDetJacobian (iQuadPt);
88  }
89 
90  localValue *= coefficient;
91 
92  // Add on the local matrix
93  mat_tmp (iDof, jDof) = localValue;
94  mat_tmp (jDof, iDof) = localValue;
95  }
96  }
97 
98  // Copying the mass in all the diagonal blocks (just one for scalar problem)
99  for (UInt iDim (0); iDim < fieldDim; ++iDim)
100  {
101  MatrixElemental::matrix_view mat = localMass.block (iDim, iDim);
102  mat += mat_tmp;
103  }
104 }
105 
106 void stiffness (MatrixElemental& localStiff,
107  const CurrentFE& stiffCFE,
108  const Real& coefficient,
109  const UInt& fieldDim)
110 {
111  const UInt nbFEDof (stiffCFE.nbFEDof() );
112  const UInt nbQuadPt (stiffCFE.nbQuadPt() );
113  Real localValue (0);
114 
115  MatrixElemental::matrix_type mat_tmp (nbFEDof, nbFEDof);
116 
117  // Loop over the basis functions (diagonal part)
118  for (UInt iDof (0); iDof < nbFEDof; ++iDof)
119  {
120  localValue = 0.0;
121 
122  // Loop on the quadrature noodes
123  for (UInt iQuadPt (0); iQuadPt < nbQuadPt; ++iQuadPt)
124  {
125  for (UInt iDim (0); iDim < stiffCFE.nbLocalCoor(); ++iDim)
126  {
127  localValue += stiffCFE.dphi (iDof, iDim, iQuadPt)
128  * stiffCFE.dphi (iDof, iDim, iQuadPt)
129  * stiffCFE.wDetJacobian (iQuadPt);
130  }
131  }
132  localValue *= coefficient;
133 
134  // Add on the local matrix
135  mat_tmp (iDof, iDof) = localValue;
136  }
137 
138  // Loop over the basis functions (extradiagonal part)
139  for (UInt iDof (0); iDof < nbFEDof ; ++iDof)
140  {
141  for (UInt jDof (0); jDof < iDof; ++jDof)
142  {
143  localValue = 0.0;
144 
145  // Loop on the quadrature nodes
146  for (UInt iQuadPt (0); iQuadPt < nbQuadPt; ++iQuadPt)
147  {
148  for (UInt iDim (0); iDim < stiffCFE.nbLocalCoor(); ++iDim)
149  {
150  localValue += stiffCFE.dphi (iDof, iDim, iQuadPt)
151  * stiffCFE.dphi (jDof, iDim, iQuadPt)
152  * stiffCFE.wDetJacobian (iQuadPt);
153  }
154  }
155 
156 
157  localValue *= coefficient;
158 
159  // Put in the local matrix
160  mat_tmp (iDof, jDof) = localValue;
161  mat_tmp (jDof, iDof) = localValue;
162  }
163  }
164 
165  // Copying the diffusion in all the diagonal blocks
166  for ( UInt iDim (0); iDim < fieldDim; ++iDim)
167  {
168  MatrixElemental::matrix_view mat = localStiff.block (iDim, iDim);
169  mat += mat_tmp;
170  }
171 }
172 
173 void advectionNewton ( Real coef, VectorElemental& vel,
174  MatrixElemental& elmat, const CurrentFE& fe,
175  int iblock, int jblock )
176 {
177  Real sum, sumDerivative;
178  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock, jblock );
179 
180  //Assemble the local matrix
181  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
182  {
183  for ( UInt j = 0; j < fe.nbFEDof(); j++ )
184  {
185  sum = 0.;
186  for ( UInt iq = 0; iq < fe.nbQuadPt(); iq++ )
187  {
188  //evaluate derivative
189  VectorElemental::vector_view velicoor = vel.block ( iblock );
190  sumDerivative = 0.;
191  for ( UInt k = 0; k < fe.nbFEDof(); k++ )
192  {
193  sumDerivative += velicoor ( k ) * fe.dphi ( k, jblock, iq );
194  }
195  //end evaluate derivative
196 
197  sum += fe.phi (j, iq) * sumDerivative * fe.phi (i, iq) * fe.weightDet ( iq );
198  }
199  mat_icomp ( i, j ) += sum * coef;
200  }
201  }
202 }
203 
204 void grad ( MatrixElemental& elmat,
205  const CurrentFE& uCFE,
206  const CurrentFE& pCFE,
207  const UInt& fieldDim)
208 {
209  const UInt nbPFEDof (pCFE.nbFEDof() );
210  const UInt nbUFEDof (uCFE.nbFEDof() );
211  const UInt nbQuadPt (pCFE.nbQuadPt() );
212  Real localValue (0);
213 
214  for (UInt iFDim (0); iFDim < fieldDim; ++iFDim)
215  {
216  MatrixElemental::matrix_view localView = elmat.block (iFDim, 0);
217 
218  for (UInt iDofP (0); iDofP < nbPFEDof; ++iDofP)
219  {
220  for (UInt iDofU (0); iDofU < nbUFEDof; ++iDofU)
221  {
222  localValue = 0.0;
223  for (UInt iQuadPt (0); iQuadPt < nbQuadPt; ++iQuadPt)
224  {
225  localValue += uCFE.dphi (iDofU, iFDim, iQuadPt)
226  * pCFE.phi (iDofP, iQuadPt)
227  * uCFE.wDetJacobian (iQuadPt);
228  }
229  localView (iDofU, iDofP) -= localValue;
230  }
231  }
232  }
233 }
234 
235 void divergence ( MatrixElemental& elmat,
236  const CurrentFE& uCFE,
237  const CurrentFE& pCFE,
238  const UInt& fieldDim,
239  const Real& coefficient)
240 {
241  const UInt nbPFEDof (pCFE.nbFEDof() );
242  const UInt nbUFEDof (uCFE.nbFEDof() );
243  const UInt nbQuadPt (pCFE.nbQuadPt() );
244  Real localValue (0);
245 
246  for (UInt iFDim (0); iFDim < fieldDim; ++iFDim)
247  {
248  MatrixElemental::matrix_view localView = elmat.block (0, iFDim);
249 
250  for (UInt iDofP (0); iDofP < nbPFEDof; ++iDofP)
251  {
252  for (UInt iDofU (0); iDofU < nbUFEDof; ++iDofU)
253  {
254  localValue = 0.0;
255  for (UInt iQuadPt (0); iQuadPt < nbQuadPt; ++iQuadPt)
256  {
257  localValue += uCFE.dphi (iDofU, iFDim, iQuadPt)
258  * pCFE.phi (iDofP, iQuadPt)
259  * uCFE.wDetJacobian (iQuadPt);
260  }
261  localView (iDofP, iDofU) -= coefficient * localValue;
262  }
263  }
264  }
265 }
266 
267 void stiffStrain (MatrixElemental& localStiff,
268  const CurrentFE& stiffCFE,
269  const Real& coefficient,
270  const UInt& fieldDim)
271 {
272  const UInt nbFEDof (stiffCFE.nbFEDof() );
273  const UInt nbQuadPt (stiffCFE.nbQuadPt() );
274  Real localValue (0);
275  const Real newCoefficient (coefficient * 0.5);
276 
277  stiffness (localStiff, stiffCFE, newCoefficient, fieldDim); // for the stiff part, we exploit the existing routine
278 
279  for ( UInt iFDim (0); iFDim < fieldDim; ++iFDim )
280  {
281  for ( UInt jFDim (0); jFDim < fieldDim; ++jFDim )
282  {
283  MatrixElemental::matrix_view localView = localStiff.block ( iFDim, jFDim );
284 
285  for ( UInt iDof (0); iDof < nbFEDof; ++iDof )
286  {
287  for ( UInt jDof (0); jDof < nbFEDof; ++jDof )
288  {
289  localValue = 0.0;
290  for ( UInt iQuadPt (0); iQuadPt < nbQuadPt; ++iQuadPt )
291  {
292  localValue += stiffCFE.dphi ( iDof, jFDim, iQuadPt )
293  * stiffCFE.dphi ( jDof, iFDim, iQuadPt )
294  * stiffCFE.wDetJacobian ( iQuadPt );
295  }
296  localView ( iDof, jDof ) += newCoefficient * localValue;
297  }
298  }
299  }
300  }
301 }
302 
303 void bodyForces (VectorElemental& localForce,
304  const CurrentFE& massRhsCFE,
305  const function_Type& fun,
306  const Real& t,
307  const UInt& fieldDim)
308 {
309 
310  const UInt nbFEDof (massRhsCFE.nbFEDof() );
311  const UInt nbQuadPt (massRhsCFE.nbQuadPt() );
312  std::vector<Real> fValues (nbQuadPt, 0.0);
313  Real localValue (0.0);
314 
315  // Assemble the local diffusion
316  for (UInt iterFDim (0); iterFDim < fieldDim; ++iterFDim)
317  {
318  VectorElemental::vector_view localView = localForce.block (iterFDim);
319 
320  // Compute the value of f in the quadrature nodes
321  for (UInt iQuadPt (0); iQuadPt < nbQuadPt; ++iQuadPt)
322  {
323  fValues[iQuadPt] = fun (t,
324  massRhsCFE.quadNode (iQuadPt, 0),
325  massRhsCFE.quadNode (iQuadPt, 1),
326  massRhsCFE.quadNode (iQuadPt, 2),
327  iterFDim);
328  }
329 
330  // Loop over the basis functions
331  for (UInt iDof (0); iDof < nbFEDof ; ++iDof)
332  {
333  localValue = 0.0;
334 
335  //Loop on the quadrature nodes
336  for (UInt iQuadPt (0); iQuadPt < nbQuadPt; ++iQuadPt)
337  {
338  localValue += fValues[iQuadPt]
339  * massRhsCFE.phi (iDof, iQuadPt)
340  * massRhsCFE.wDetJacobian (iQuadPt);
341  }
342 
343  // Add on the local matrix
344  localView (iDof) = localValue;
345  }
346  }
347 
348 }
349 
350 //
351 //----------------------------------------------------------------------
352 // Element matrix operator
353 //----------------------------------------------------------------------
354 //Real coef(t,x,y,z,u)
355 /*
356  Mass matrix: \int coef(t,x,y,z,u) v_i v_j
357 */
358 
359 //
360 // coeff*Mass
361 //
362 void mass ( Real coef, MatrixElemental& elmat, const CurrentFE& fe,
363  int iblock, int jblock )
364 /*
365  Mass matrix: \int v_i v_j
366 */
367 {
368  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
369  UInt i, ig;
370  int iloc, jloc;
371  Real s, coef_s;
372  //
373  // diagonal
374  //
375  for ( i = 0; i < fe.nbDiag(); i++ )
376  {
377  iloc = fe.patternFirst ( i );
378  s = 0;
379  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
380  {
381  s += fe.phi ( iloc, ig ) * fe.phi ( iloc, ig ) * fe.weightDet ( ig );
382  }
383  mat ( iloc, iloc ) += coef * s;
384  }
385  //
386  // extra diagonal
387  //
388  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
389  {
390  iloc = fe.patternFirst ( i );
391  jloc = fe.patternSecond ( i );
392  s = 0;
393  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
394  {
395  s += fe.phi ( iloc, ig ) * fe.phi ( jloc, ig ) * fe.weightDet ( ig );
396  }
397  coef_s = coef * s;
398  mat ( iloc, jloc ) += coef_s;
399  mat ( jloc, iloc ) += coef_s;
400  }
401 }
402 //
403 // coeff*Mass
404 //
405 void mass ( Real coef, MatrixElemental& elmat, const CurrentFE& fe,
406  int iblock, int jblock, UInt nb )
407 /*
408  Mass matrix: \int v_i v_j (nb blocks on the diagonal, nb>1)
409 */
410 {
411  Matrix mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
412  UInt i, ig;
413  int iloc, jloc;
414  Real s, coef_s;
415  mat_tmp = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
416  //
417  // diagonal
418  //
419  for ( i = 0; i < fe.nbDiag(); i++ )
420  {
421  iloc = fe.patternFirst ( i );
422  s = 0;
423  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
424  {
425  s += fe.phi ( iloc, ig ) * fe.phi ( iloc, ig ) * fe.weightDet ( ig );
426  }
427  mat_tmp ( iloc, iloc ) += coef * s;
428  }
429  //
430  // extra diagonal
431  //
432  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
433  {
434  iloc = fe.patternFirst ( i );
435  jloc = fe.patternSecond ( i );
436  s = 0;
437  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
438  {
439  s += fe.phi ( iloc, ig ) * fe.phi ( jloc, ig ) * fe.weightDet ( ig );
440  }
441  coef_s = coef * s;
442  mat_tmp ( iloc, jloc ) += coef_s;
443  mat_tmp ( jloc, iloc ) += coef_s;
444  }
445  // copy on the components
446  for ( UInt icomp = 0; icomp < nb; icomp++ )
447  {
448  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icomp, jblock + icomp );
449  for ( i = 0; i < fe.nbDiag(); i++ )
450  {
451  iloc = fe.patternFirst ( i );
452  mat_icomp ( iloc, iloc ) += mat_tmp ( iloc, iloc );
453  }
454  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
455  {
456  iloc = fe.patternFirst ( i );
457  jloc = fe.patternSecond ( i );
458  mat_icomp ( iloc, jloc ) += mat_tmp ( iloc, jloc );
459  mat_icomp ( jloc, iloc ) += mat_tmp ( jloc, iloc );
460  }
461  }
462 }
463 
464 //
465 // coeff[q]*Mass
466 //
467 void mass ( const std::vector<Real>& coef, MatrixElemental& elmat, const CurrentFE& fe,
468  int iblock, int jblock, UInt nb )
469 /*
470  Mass matrix: \int v_i v_j (nb blocks on the diagonal, nb>1)
471 */
472 {
473  Matrix mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
474  UInt i, ig;
475  int iloc, jloc;
476  Real s;//, coef_s;
477  mat_tmp = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
478  //
479  // diagonal
480  //
481  for ( i = 0; i < fe.nbDiag(); i++ )
482  {
483  iloc = fe.patternFirst ( i );
484  s = 0;
485  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
486  {
487  s += coef[ig] * fe.phi ( iloc, ig ) * fe.phi ( iloc, ig ) * fe.weightDet ( ig );
488  }
489  mat_tmp ( iloc, iloc ) = s;
490  }
491  //
492  // extra diagonal
493  //
494  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
495  {
496  iloc = fe.patternFirst ( i );
497  jloc = fe.patternSecond ( i );
498  s = 0;
499  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
500  {
501  s += coef[ig] * fe.phi ( iloc, ig ) * fe.phi ( jloc, ig ) * fe.weightDet ( ig );
502  }
503  mat_tmp ( iloc, jloc ) += s;
504  mat_tmp ( jloc, iloc ) += s;
505  }
506  // copy on the components
507  for ( UInt icomp = 0; icomp < nb; icomp++ )
508  {
509  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icomp, jblock + icomp );
510  for ( i = 0; i < fe.nbDiag(); i++ )
511  {
512  iloc = fe.patternFirst ( i );
513  mat_icomp ( iloc, iloc ) += mat_tmp ( iloc, iloc );
514  }
515  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
516  {
517  iloc = fe.patternFirst ( i );
518  jloc = fe.patternSecond ( i );
519  mat_icomp ( iloc, jloc ) += mat_tmp ( iloc, jloc );
520  mat_icomp ( jloc, iloc ) += mat_tmp ( jloc, iloc );
521  }
522  }
523 }
524 
525 
526 
527 
528 
529 
530 void stiff_divgrad ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
531 {
532 
533  double s;
534  // div u^k at quad pts
535  std::vector<Real> duk (fe.nbQuadPt() );
536 
537  // loop on quadrature points
538  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
539  {
540  s = 0;
541  // loop on space coordinates
542  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
543  {
544  // s = 0.0; Alessandro
545  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
546  {
547  s += fe.phiDer ( i, icoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // costruzione di \div u^k at a quadrature point
548  }
549 
550  // duk[ ig ] = s; Alessandro
551  }// chiude il ciclo su icoor
552  duk[ ig ] = s;
553  }// chiude il ciclo su ig
554 
555  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
556 
557  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
558  {
559  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
560  {
561  s = 0.0;
562  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
563  {
564  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
565  {
566  s += duk[ ig ] * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
567  }
568  }
569  mat_tmp ( i, j ) = coef * s;
570  }
571  }
572 
573  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
574  {
575  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
576  mat += mat_tmp;
577  }
578 
579 }
580 
581 // Stiffness matrix: coef * ( (\div u) \grad u_k : \grad v ) controllato!!!
582 void stiff_divgrad_2 ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
583 {
584 
585  double s;
586  boost::multi_array<Real, 3> guk (boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
587 
588  // loop on quadrature points
589  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
590  {
591  // loop on space coordinates
592  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
593  {
594  // loop on space coordinates
595  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
596  {
597  s = 0.0;
598  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
599  {
600  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
601  }
602  guk[ icoor ][ jcoor ][ ig ] = s;
603  }
604  }
605  }
606 
607  //
608  // blocks (icoor,jcoor) of elmat
609  //
610 
611  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
612  {
613  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
614  {
615  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
616 
617  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
618  {
619  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
620  {
621  s = 0;
622  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
623  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
624  {
625  s += fe.phiDer ( j, jcoor, ig ) * guk[ icoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
626  }
627  mat ( i, j ) += coef * s;
628  }
629  }
630  }
631  }
632 }
633 
634 // Stiffness matrix: coef * ( \grad u_k : \grad u_k) *( \grad u : \grad v ) controllato!!!
635 void stiff_gradgrad ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
636 {
637 
638  double s, s1;
639  std::vector<Real> gguk (fe.nbQuadPt() );
640 
641  // loop on quadrature points
642  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
643  {
644  s = 0;
645  // loop on space coordinates
646  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
647  {
648  for ( UInt l = 0; l < fe.nbLocalCoor(); l++ )
649  {
650  s1 = 0;
651  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
652  {
653  s1 += fe.phiDer ( i, l, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ];
654  }
655  s += s1 * s1;
656  }
657  }// chiude il ciclo su icoor
658  gguk[ ig ] = s;
659  }// chiude il ciclo su ig
660 
661  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
662 
663  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
664  {
665  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
666  {
667  s = 0.0;
668  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
669  {
670  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
671  {
672  s += gguk[ ig ] * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
673  }
674  }
675  mat_tmp ( i, j ) = coef * s;
676  }
677  }
678 
679  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
680  {
681  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
682  mat += mat_tmp;
683  }
684 }
685 
686 // Stiffness matrix: coef * ( \grad u_k : \grad u) *( \grad u_k : \grad v ) controllato!!!
687 void stiff_gradgrad_2 ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
688 {
689 
690  double s;
691  boost::multi_array<Real, 3> guk (
692  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
693 
694  // loop on quadrature points
695  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
696  {
697  // loop on space coordinates
698  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
699  {
700  // loop on space coordinates
701  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
702  {
703  s = 0.0;
704  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
705  {
706  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
707  }
708  guk[ icoor ][ jcoor ][ ig ] = s;
709  }
710  }
711  }
712 
713  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
714  {
715  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
716  {
717  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
718 
719  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
720  {
721  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
722  {
723  s = 0.0;
724  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
725  {
726  for ( UInt l = 0; l < fe.nbLocalCoor(); ++l )
727  {
728  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
729  {
730  s += guk[ jcoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) * guk[ icoor ][ k ][ ig ] * fe.phiDer (i, k, ig ) * fe.weightDet ( ig ); // il ciclo sui nodi di quadratura
731  }
732  } // fa l'integrale
733  }
734  mat ( i, j ) += coef * s;
735  }
736  }
737  }
738  }
739 }
740 
741 // Stiffness matrix: coef * ( \grad d^k \grad d : \grad v )controllato!!!
742 void stiff_dergrad_gradbis ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
743 {
744 
745  double s;
746  boost::multi_array<Real, 3> guk (
747  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
748 
749  // loop on quadrature points
750  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
751  {
752 
753  // loop on space coordinates
754  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
755  {
756 
757  // loop on space coordinates
758  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
759  {
760  s = 0.0;
761  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
762  {
763  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
764  }
765  guk[ icoor ][ jcoor ][ ig ] = s;
766  }
767  }
768  }
769  //
770  // blocks (icoor,jcoor) of elmat
771  //
772 
773  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
774  {
775  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
776  {
777 
778  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
779 
780  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
781  {
782  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
783  {
784  s = 0.0;
785  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
786  {
787  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
788  {
789  s += guk[ icoor ][ jcoor ][ ig ] * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
790  }
791  }
792  mat ( i, j ) += coef * s;
793  }
794  }
795  }
796  }
797 }
798 
799 // Stiffness matrix: coef * ( \grad u^k [\grad d]^T : \grad v ) controllato!!!
800 void stiff_dergrad_gradbis_Tr ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
801 {
802 
803  double s;
804  boost::multi_array<Real, 3> guk (
805  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
806 
807  // loop on quadrature points
808  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
809  {
810 
811  // loop on space coordinates
812  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
813  {
814 
815  // loop on space coordinates
816  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
817  {
818  s = 0.0;
819  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
820  {
821  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
822  }
823  guk[ icoor ][ jcoor ][ ig ] = s;
824  }
825  }
826  }
827  //
828  // blocks (icoor,jcoor) of elmat
829  //
830 
831  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
832  {
833  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
834  {
835 
836  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
837 
838  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
839  {
840  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
841  {
842  s = 0.0;
843  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
844  {
845  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
846  {
847  s += guk[ icoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.phiDer ( i, jcoor, ig ) * fe.weightDet ( ig );
848  }
849  }
850  mat ( i, j ) += coef * s;
851  }
852  }
853  }
854  }
855 }
856 
857 // Stiffness matrix: coef * ( \grad \delta d \grad d^k : \grad v ) controllato!!!
858 void stiff_dergrad_gradbis_2 ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
859 {
860 
861  double s;
862  boost::multi_array<Real, 3> guk (
863  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
864 
865  // loop on quadrature points
866  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
867  {
868 
869  // loop on space coordinates
870  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
871  {
872 
873  // loop on space coordinates
874  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
875  {
876  s = 0.0;
877  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
878  {
879  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
880  }
881  guk[ icoor ][ jcoor ][ ig ] = s;
882  }
883  }
884  }
885  //
886  // blocks (icoor,jcoor) of elmat
887  //
888  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
889 
890  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
891  {
892  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
893  {
894  s = 0.0;
895  for ( UInt l = 0; l < fe.nbLocalCoor(); ++l )
896  {
897  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
898  {
899  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
900  {
901  s += guk[ l ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
902  }
903  }
904  }
905  mat_tmp ( i, j ) = coef * s;
906  }
907  }
908 
909  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
910  {
911  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
912  mat += mat_tmp;
913  }
914 }
915 
916 // Stiffness matrix: coef * ( \grad \delta d [\grad d^k]^T : \grad v )
917 void stiff_dergrad_gradbis_Tr_2 ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
918 {
919 
920  double s;
921  boost::multi_array<Real, 3> guk (
922  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
923 
924  // loop on quadrature points
925  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
926  {
927 
928  // loop on space coordinates
929  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
930  {
931 
932  // loop on space coordinates
933  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
934  {
935  s = 0.0;
936  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
937  {
938  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
939  }
940  guk[ icoor ][ jcoor ][ ig ] = s;
941  }
942  }
943  }
944  //
945  // blocks (icoor,jcoor) of elmat
946  //
947  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
948 
949  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
950  {
951  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
952  {
953  s = 0.0;
954  for ( UInt l = 0; l < fe.nbLocalCoor(); ++l )
955  {
956  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
957  {
958  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
959  {
960  s += guk[ k ][ l ][ ig ] * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
961  }
962  }
963  }
964  mat_tmp ( i, j ) = coef * s;
965  }
966  }
967 
968  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor ) // copia del blocco sulla diagonale
969  {
970  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
971  mat += mat_tmp;
972  }
973 }
974 
975 // Stiffness matrix: coef * ( \grad u^k [\grad u^k]^T \grad u : \grad v ) controllato!!!
976 void stiff_gradgradTr_gradbis ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
977 {
978 
979  double s;
980  boost::multi_array<Real, 3> guk_gukT (
981  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
982 
983  // loop on quadrature points
984  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
985  {
986 
987  // loop on space coordinates
988  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
989  {
990 
991  // loop on space coordinates
992  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
993  {
994  s = 0.0;
995  for ( UInt n = 0; n < fe.nbLocalCoor(); n++ )
996  {
997  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
998  {
999  for ( UInt j = 0; j < fe.nbFEDof(); j++ )
1000  {
1001  s += fe.phiDer ( i, n, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ] * fe.phiDer ( j, n, ig ) * uk_loc.vec() [ j + jcoor * fe.nbFEDof() ] ; // \grad u^k [\grad u^k]^T at each quadrature point
1002  }
1003  }
1004  }
1005  guk_gukT[ icoor ][ jcoor ][ ig ] = s;
1006  }
1007  }
1008  }
1009  //
1010  // blocks (icoor,jcoor) of elmat
1011  //
1012 
1013  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
1014  {
1015  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
1016  {
1017 
1018  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor ); // estrae il blocco (icoor, jcoor)
1019 
1020  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1021  {
1022  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1023  {
1024  s = 0;
1025  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
1026  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1027  {
1028  s += fe.phiDer ( i, k, ig ) * guk_gukT[ icoor ][ jcoor ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1029  }
1030  mat ( i, j ) += coef * s;
1031  }
1032  }
1033  }
1034  }
1035 }
1036 
1037 // Stiffness matrix: coef * ( \grad u^k [\grad u]^T \grad u^k : \grad v )controllato!!!
1038 void stiff_gradgradTr_gradbis_2 ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
1039 {
1040 
1041  double s;
1042  boost::multi_array<Real, 3> guk (
1043  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
1044 
1045  // loop on quadrature points
1046  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
1047  {
1048 
1049  // loop on space coordinates
1050  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
1051  {
1052 
1053  // loop on space coordinates
1054  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
1055  {
1056  s = 0.0;
1057  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
1058  {
1059  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
1060  }
1061  guk[ icoor ][ jcoor ][ ig ] = s;
1062  }
1063  }
1064  }
1065 
1066  //
1067  // blocks (icoor,jcoor) of elmat
1068  //
1069  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
1070  {
1071  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
1072  {
1073 
1074  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor ); // estrae il blocco (icoor, jcoor)
1075 
1076  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1077  {
1078  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1079  {
1080  s = 0;
1081  for ( UInt l = 0; l < fe.nbLocalCoor(); ++l )
1082  {
1083  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
1084  {
1085  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1086  {
1087  s += guk[ icoor ][ l ][ ig ] * guk[ jcoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
1088  }
1089  }
1090  }
1091  mat ( i, j ) += coef * s;
1092  }
1093  }
1094  }
1095  }
1096 }
1097 
1098 // Stiffness matrix: coef * ( \grad u [\grad u^k]^T \grad u^k : \grad v )controllato!!!
1099 void stiff_gradgradTr_gradbis_3 ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
1100 {
1101 
1102  double s;
1103  boost::multi_array<Real, 3> guk_gukT (
1104  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
1105 
1106  // attenzione in questa funzione si deve usare il trasposto
1107  // loop on quadrature points // (\grad u^k [\grad u^k]^T )^T
1108  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
1109  {
1110 
1111  // loop on space coordinates
1112  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
1113  {
1114 
1115  // loop on space coordinates
1116  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
1117  {
1118  s = 0.0;
1119  for ( UInt n = 0; n < fe.nbLocalCoor(); n++ )
1120  {
1121  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
1122  {
1123  for ( UInt j = 0; j < fe.nbFEDof(); j++ )
1124  {
1125  s += fe.phiDer ( i, n, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ] * fe.phiDer ( j, n, ig ) * uk_loc.vec() [ j + jcoor * fe.nbFEDof() ] ; // \grad u^k [\grad u^k]^T at each quadrature point
1126  }
1127  }
1128  }
1129  guk_gukT[ icoor ][ jcoor ][ ig ] = s;
1130  }
1131  }
1132  }
1133 
1134  //
1135  // blocks (icoor,jcoor) of elmat
1136  //
1137  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
1138 
1139  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1140  {
1141  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1142  {
1143  s = 0.0;
1144  for ( UInt l = 0; l < fe.nbLocalCoor(); ++l )
1145  {
1146  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
1147  {
1148  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1149  {
1150  s += guk_gukT[ k ][ l ][ ig ] * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
1151  }
1152  }
1153  }
1154  mat_tmp ( i, j ) = coef * s;
1155  }
1156  }
1157 
1158  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor ) // copia del blocco sulla diagonale
1159  {
1160  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
1161  mat += mat_tmp;
1162  }
1163 }
1164 
1165 
1166 
1167 
1168 
1169 
1170 
1171 void ipstab_grad ( const Real coef,
1172  MatrixElemental& elmat,
1173  const CurrentFE& fe1,
1174  const CurrentFE& fe2,
1175  const CurrentFEManifold& bdfe,
1176  int iblock, int jblock )
1177 {
1178  /*
1179  Interior penalty stabilization: coef*\int_{face} grad u1_i . grad v1_j
1180  */
1181 
1182  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
1183 
1184 
1185 
1186  Real sum, sum1, sum2;
1187  UInt icoor, jcoor, i, j;
1188  UInt ig;
1189  Real x[ 3 ], rx1[ 3 ], drp1[ 3 ], rx2[ 3 ], drp2[ 3 ];
1190 
1191  boost::multi_array<Real, 3> phid1 (
1192  boost::extents[fe1.nbFEDof()][fe1.nbLocalCoor()][bdfe.nbQuadPt()]);
1193  boost::multi_array<Real, 3> phid2 (
1194  boost::extents[fe2.nbFEDof()][fe2.nbLocalCoor()][bdfe.nbQuadPt()]);
1195 
1196  Real b1[ 3 ], b2[ 3 ];
1197 
1198  fe1.coorMap ( b1[ 0 ], b1[ 1 ], b1[ 2 ], 0, 0, 0 ); // translation fe1
1199  fe2.coorMap ( b2[ 0 ], b2[ 1 ], b2[ 2 ], 0, 0, 0 ); // translation fe2
1200 
1201  for ( UInt ig = 0; ig < bdfe.nbQuadPt(); ++ig )
1202  {
1203  // first derivatives on quadrature points
1204  bdfe.coorQuadPt ( x[ 0 ], x[ 1 ], x[ 2 ], ig ); // quadrature points coordinates
1205 
1206  // local coordinates of the quadrature point
1207  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1208  {
1209  sum1 = 0;
1210  sum2 = 0;
1211  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1212  {
1213  sum1 += fe1.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b1[ jcoor ] );
1214  sum2 += fe2.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b2[ jcoor ] );
1215  }
1216  rx1[ icoor ] = sum1;
1217  rx2[ icoor ] = sum2;
1218  }
1219 
1220  for ( i = 0; i < fe1.nbFEDof(); ++i )
1221  {
1222 
1223  // first derivative on the reference element
1224  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1225  {
1226  drp1[ icoor ] = fe1.refFE().dPhi ( i, icoor, rx1[ 0 ], rx1[ 1 ], rx1[ 2 ] );
1227  drp2[ icoor ] = fe2.refFE().dPhi ( i, icoor, rx2[ 0 ], rx2[ 1 ], rx2[ 2 ] );
1228  }
1229 
1230  // first derivative on the current element
1231  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1232  {
1233  sum1 = 0;
1234  sum2 = 0;
1235  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1236  {
1237  sum1 += fe1.tInvJac ( icoor, jcoor, 0 ) * drp1[ jcoor ];
1238  sum2 += fe2.tInvJac ( icoor, jcoor, 0 ) * drp2[ jcoor ];
1239  }
1240  phid1[ i ][ icoor ][ ig ] = sum1;
1241  phid2[ i ][ icoor ][ ig ] = sum2;
1242  }
1243  }
1244  }
1245 
1246 
1247 
1248  // Loop on rows
1249  for ( i = 0; i < fe1.nbFEDof(); ++i )
1250  {
1251  // Loop on columns
1252  for ( j = 0; j < fe2.nbFEDof(); ++j )
1253  {
1254  sum = 0.0;
1255  // Loop on coordinates
1256  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1257  for ( ig = 0; ig < bdfe.nbQuadPt() ; ++ig )
1258  {
1259  sum += phid1[ i ][ icoor ][ ig ] * phid2[ j ][ icoor ][ ig ] * bdfe.wRootDetMetric ( ig );
1260  }
1261  mat ( i, j ) += coef * sum;
1262  }
1263  }
1264 
1265 }
1266 
1267 
1268 
1269 
1270 
1271 
1272 void ipstab_grad ( const Real coef,
1273  MatrixElemental& elmat,
1274  const CurrentFE& fe1,
1275  const CurrentFE& fe2,
1276  const CurrentFEManifold& bdfe,
1277  int iblock, int jblock,
1278  int nb )
1279 {
1280  /*
1281  Interior penalty stabilization: coef*\int_{face} grad u1_i . grad v1_j
1282  */
1283 
1284 
1285  MatrixElemental::matrix_type mat_tmp ( fe1.nbFEDof(), fe2.nbFEDof() );
1286 
1287 
1288  Real sum, sum1, sum2;
1289  UInt icoor, jcoor, i, j;
1290  UInt ig;
1291  Real x[ 3 ], rx1[ 3 ], drp1[ 3 ], rx2[ 3 ], drp2[ 3 ];
1292 
1293  boost::multi_array<Real, 3> phid1 (
1294  boost::extents[fe1.nbFEDof()][fe1.nbLocalCoor()][bdfe.nbQuadPt()]);
1295  boost::multi_array<Real, 3> phid2 (
1296  boost::extents[fe2.nbFEDof()][fe2.nbLocalCoor()][bdfe.nbQuadPt()]);
1297 
1298  Real b1[ 3 ], b2[ 3 ];
1299 
1300  fe1.coorMap ( b1[ 0 ], b1[ 1 ], b1[ 2 ], 0, 0, 0 ); // translation fe1
1301  fe2.coorMap ( b2[ 0 ], b2[ 1 ], b2[ 2 ], 0, 0, 0 ); // translation fe2
1302 
1303  for ( UInt ig = 0; ig < bdfe.nbQuadPt(); ++ig )
1304  {
1305  // first derivatives on quadrature points
1306  bdfe.coorQuadPt ( x[ 0 ], x[ 1 ], x[ 2 ], ig ); // quadrature points coordinates
1307 
1308  // local coordinates of the quadrature point
1309  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1310  {
1311  sum1 = 0;
1312  sum2 = 0;
1313  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1314  {
1315  sum1 += fe1.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b1[ jcoor ] );
1316  sum2 += fe2.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b2[ jcoor ] );
1317  }
1318  rx1[ icoor ] = sum1;
1319  rx2[ icoor ] = sum2;
1320  }
1321 
1322  for ( i = 0; i < fe1.nbFEDof(); ++i )
1323  {
1324 
1325  // first derivative on the reference element
1326  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1327  {
1328  drp1[ icoor ] = fe1.refFE().dPhi ( i, icoor, rx1[ 0 ], rx1[ 1 ], rx1[ 2 ] );
1329  drp2[ icoor ] = fe2.refFE().dPhi ( i, icoor, rx2[ 0 ], rx2[ 1 ], rx2[ 2 ] );
1330  }
1331 
1332  // first derivative on the current element
1333  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1334  {
1335  sum1 = 0;
1336  sum2 = 0;
1337  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1338  {
1339  sum1 += fe1.tInvJac ( icoor, jcoor, 0 ) * drp1[ jcoor ];
1340  sum2 += fe2.tInvJac ( icoor, jcoor, 0 ) * drp2[ jcoor ];
1341  }
1342  phid1[ i ][ icoor ][ ig ] = sum1;
1343  phid2[ i ][ icoor ][ ig ] = sum2;
1344  }
1345  }
1346  }
1347 
1348 
1349  // Loop on rows
1350  for ( i = 0; i < fe1.nbFEDof(); ++i )
1351  {
1352  // Loop on columns
1353  for ( j = 0; j < fe2.nbFEDof(); ++j )
1354  {
1355  sum = 0.0;
1356  // Loop on coordinates
1357  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1358  for ( ig = 0; ig < bdfe.nbQuadPt() ; ++ig )
1359  {
1360  sum += phid1[ i ][ icoor ][ ig ] * phid2[ j ][ icoor ][ ig ] * bdfe.wRootDetMetric ( ig );
1361  }
1362  mat_tmp ( i, j ) = coef * sum;
1363  }
1364  }
1365 
1366 
1367  // copy on the components
1368  for ( int icomp = 0; icomp < nb; icomp++ )
1369  {
1370  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icomp, jblock + icomp );
1371  mat_icomp += mat_tmp;
1372  }
1373 }
1374 
1375 
1376 
1377 
1378 
1379 void ipstab_bgrad ( const Real coef,
1380  MatrixElemental& elmat,
1381  const CurrentFE& fe1,
1382  const CurrentFE& fe2,
1383  const VectorElemental& beta,
1384  const CurrentFEManifold& bdfe,
1385  int iblock, int jblock,
1386  int nb )
1387 {
1388  /*
1389  Interior penalty stabilization: coef*\int_{face} (\beta1 . grad u1_i) . (\beta2 . grad v2_j)
1390  */
1391 
1392  MatrixElemental::matrix_type mat_tmp ( fe1.nbFEDof(), fe2.nbFEDof() );
1393 
1394  Real sum, sum1, sum2;
1395  UInt i, j;
1396  UInt icoor, jcoor;
1397  UInt ig;
1398 
1399  //
1400  // convection velocity \beta on the boundary quadrature points
1401  //
1402  boost::multi_array<Real, 2> b (
1403  boost::extents[fe1.nbLocalCoor()][bdfe.nbQuadPt()]);
1404 
1405  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1406  {
1407  for ( ig = 0; ig < bdfe.nbQuadPt(); ig++ )
1408  {
1409  sum = 0;
1410  for ( i = 0; i < bdfe.nbFEDof(); ++i )
1411  {
1412  sum += bdfe.phi ( i, ig ) * beta.vec() [ icoor * bdfe.nbLocalCoor() + i ];
1413  }
1414  b[ icoor ][ ig ] = sum;
1415  }
1416  }
1417 
1418 
1419  //
1420  // shape function first derivaties on the boundary quadrature points
1421  //
1422  // this should be improved!!!
1423  //
1424 
1425  Real x[ 3 ], rx1[ 3 ], drp1[ 3 ], rx2[ 3 ], drp2[ 3 ];
1426 
1427  boost::multi_array<Real, 3> phid1 (
1428  boost::extents[fe1.nbFEDof()][fe1.nbLocalCoor()][bdfe.nbQuadPt()]);
1429  boost::multi_array<Real, 3> phid2 (
1430  boost::extents[fe2.nbFEDof()][fe2.nbLocalCoor()][bdfe.nbQuadPt()]);
1431 
1432  Real b1[ 3 ], b2[ 3 ];
1433 
1434  fe1.coorMap ( b1[ 0 ], b1[ 1 ], b1[ 2 ], 0, 0, 0 ); // translation fe1
1435  fe2.coorMap ( b2[ 0 ], b2[ 1 ], b2[ 2 ], 0, 0, 0 ); // translation fe2
1436 
1437  for ( UInt ig = 0; ig < bdfe.nbQuadPt(); ++ig )
1438  {
1439  // first derivatives on quadrature points
1440  bdfe.coorQuadPt ( x[ 0 ], x[ 1 ], x[ 2 ], ig ); // quadrature points coordinates
1441 
1442  // local coordinates of the quadrature point
1443  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1444  {
1445  sum1 = 0;
1446  sum2 = 0;
1447  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1448  {
1449  sum1 += fe1.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b1[ jcoor ] );
1450  sum2 += fe2.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b2[ jcoor ] );
1451  }
1452  rx1[ icoor ] = sum1;
1453  rx2[ icoor ] = sum2;
1454  }
1455 
1456  for ( i = 0; i < fe1.nbFEDof(); ++i )
1457  {
1458 
1459  // first derivative on the reference element
1460  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1461  {
1462  drp1[ icoor ] = fe1.refFE().dPhi ( i, icoor, rx1[ 0 ], rx1[ 1 ], rx1[ 2 ] );
1463  drp2[ icoor ] = fe2.refFE().dPhi ( i, icoor, rx2[ 0 ], rx2[ 1 ], rx2[ 2 ] );
1464  }
1465 
1466  // first derivative on the current element
1467  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1468  {
1469  sum1 = 0;
1470  sum2 = 0;
1471  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1472  {
1473  sum1 += fe1.tInvJac ( icoor, jcoor, 0 ) * drp1[ jcoor ];
1474  sum2 += fe2.tInvJac ( icoor, jcoor, 0 ) * drp2[ jcoor ];
1475  }
1476  phid1[ i ][ icoor ][ ig ] = sum1;
1477  phid2[ i ][ icoor ][ ig ] = sum2;
1478  }
1479  }
1480  }
1481 
1482  // Loop on rows
1483  for ( i = 0; i < fe1.nbFEDof(); ++i )
1484  {
1485  // Loop on columns
1486  for ( j = 0; j < fe2.nbFEDof(); ++j )
1487  {
1488  sum = 0.0;
1489  // Loop on coordinates
1490  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1491  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1492  for ( ig = 0; ig < bdfe.nbQuadPt(); ig++ )
1493  sum += phid1[ i ][ icoor ][ ig ] * phid2[ j ][ jcoor ][ ig ]
1494  * b[ icoor ][ ig ] * b[ jcoor ][ ig ]
1495  * bdfe.wRootDetMetric ( ig );
1496  mat_tmp ( i, j ) = coef * sum;
1497  }
1498  }
1499 
1500  // copy on the components
1501  for ( int icomp = 0; icomp < nb; icomp++ )
1502  {
1503  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icomp, jblock + icomp );
1504  mat_icomp += mat_tmp;
1505  }
1506 
1507 }
1508 
1509 
1510 
1511 
1512 
1513 void ipstab_div ( const Real coef, MatrixElemental& elmat, const CurrentFE& fe1, const CurrentFE& fe2,
1514  const CurrentFEManifold& bdfe, int iblock, int jblock )
1515 {
1516  /*
1517  Interior penalty stabilization: coef*\int_{face} div u . div v
1518  */
1519 
1520  Real sum, sum1, sum2;
1521  UInt i, j, icoor, jcoor;
1522  UInt ig;
1523  Real x[ 3 ], rx1[ 3 ], drp1[ 3 ], rx2[ 3 ], drp2[ 3 ];
1524 
1525  boost::multi_array<Real, 3> phid1 (
1526  boost::extents[fe1.nbFEDof()][fe1.nbLocalCoor()][bdfe.nbQuadPt()]);
1527  boost::multi_array<Real, 3> phid2 (
1528  boost::extents[fe2.nbFEDof()][fe2.nbLocalCoor()][bdfe.nbQuadPt()]);
1529 
1530  Real b1[ 3 ], b2[ 3 ];
1531 
1532  fe1.coorMap ( b1[ 0 ], b1[ 1 ], b1[ 2 ], 0, 0, 0 ); // translation fe1
1533  fe2.coorMap ( b2[ 0 ], b2[ 1 ], b2[ 2 ], 0, 0, 0 ); // translation fe2
1534 
1535  for ( UInt ig = 0; ig < bdfe.nbQuadPt(); ++ig )
1536  {
1537  // first derivatives on quadrature points
1538  bdfe.coorQuadPt ( x[ 0 ], x[ 1 ], x[ 2 ], ig ); // quadrature points coordinates
1539 
1540  // local coordinates of the quadrature point
1541  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1542  {
1543  sum1 = 0;
1544  sum2 = 0;
1545  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1546  {
1547  sum1 += fe1.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b1[ jcoor ] );
1548  sum2 += fe2.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b2[ jcoor ] );
1549  }
1550  rx1[ icoor ] = sum1;
1551  rx2[ icoor ] = sum2;
1552  }
1553 
1554  for ( i = 0; i < fe1.nbFEDof(); ++i )
1555  {
1556 
1557  // first derivative on the reference element
1558  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1559  {
1560  drp1[ icoor ] = fe1.refFE().dPhi ( i, icoor, rx1[ 0 ], rx1[ 1 ], rx1[ 2 ] );
1561  drp2[ icoor ] = fe2.refFE().dPhi ( i, icoor, rx2[ 0 ], rx2[ 1 ], rx2[ 2 ] );
1562  }
1563 
1564  // first derivative on the current element
1565  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1566  {
1567  sum1 = 0;
1568  sum2 = 0;
1569  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1570  {
1571  sum1 += fe1.tInvJac ( icoor, jcoor, 0 ) * drp1[ jcoor ];
1572  sum2 += fe2.tInvJac ( icoor, jcoor, 0 ) * drp2[ jcoor ];
1573  }
1574  phid1[ i ][ icoor ][ ig ] = sum1;
1575  phid2[ i ][ icoor ][ ig ] = sum2;
1576  }
1577  }
1578  }
1579 
1580  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1581  {
1582  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1583  {
1584  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icoor, jblock + jcoor );
1585  // Loop on rows
1586  for ( i = 0; i < fe1.nbFEDof(); ++i )
1587  {
1588  // Loop on columns
1589  for ( j = 0; j < fe2.nbFEDof(); ++j )
1590  {
1591  sum = 0.0;
1592  for ( ig = 0; ig < bdfe.nbQuadPt(); ig++ )
1593  {
1594  sum += phid1[ i ][ icoor ][ ig ] * phid2[ j ][ jcoor ][ ig ] * bdfe.wRootDetMetric ( ig );
1595  }
1596  mat_icomp ( i, j ) += coef * sum;
1597  }
1598  }
1599  }
1600  }
1601 
1602 }
1603 
1604 void ipstab_bagrad ( const Real coef, MatrixElemental& elmat,
1605  const CurrentFE& fe1, const CurrentFE& fe2,
1606  const VectorElemental& beta, const CurrentFEManifold& bdfe,
1607  int iblock, int jblock )
1608 {
1609 
1610  // Interior penalty stabilization:
1611  // coef < |\beta . n|^2 / |\beta| \grad p1, \grad q2 >
1612 
1613  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
1614 
1615  Real sum, sum1, sum2;
1616  UInt icoor, jcoor;
1617  UInt i, j, ig;
1618 
1619  boost::multi_array<Real, 3> phid1 (
1620  boost::extents[fe1.nbFEDof()][fe1.nbLocalCoor()][bdfe.nbQuadPt()]);
1621  boost::multi_array<Real, 3> phid2 (
1622  boost::extents[fe2.nbFEDof()][fe2.nbLocalCoor()][bdfe.nbQuadPt()]);
1623 
1624  std::vector<Real> x (3), rx1 (3), drp1 (3), rx2 (3), drp2 (3);
1625  std::vector<Real> b1 (3), b2 (3);
1626 
1627  fe1.coorMap ( b1[ 0 ], b1[ 1 ], b1[ 2 ], 0, 0, 0 ); // translation fe1
1628  fe2.coorMap ( b2[ 0 ], b2[ 1 ], b2[ 2 ], 0, 0, 0 ); // translation fe2
1629 
1630  //
1631  // convection velocity term |\beta . n|^2 / |\beta|
1632  // on the boundary quadrature points
1633  //
1634  std::vector<Real> ba2 (bdfe.nbQuadPt() );
1635 
1636  for ( ig = 0; ig < bdfe.nbQuadPt(); ig++ )
1637  {
1638  sum1 = 0;
1639  sum2 = 0;
1640  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1641  {
1642  for ( i = 0; i < bdfe.nbLocalCoor(); ++i )
1643  {
1644  Real betaLoc = bdfe.phi ( i, ig ) *
1645  beta.vec() [ icoor * bdfe.nbLocalCoor() + i ];
1646  sum1 += betaLoc * bdfe.normal (icoor, ig);
1647  sum2 += betaLoc * betaLoc;
1648  }
1649  }
1650  ba2[ ig ] = sum2 == 0 ? 0 : sum1 * sum1 / std::pow ( sum2, 0.5 );
1651  }
1652 
1653  for ( UInt ig = 0; ig < bdfe.nbQuadPt(); ++ig )
1654  {
1655  // first derivatives on quadrature points
1656  bdfe.coorQuadPt ( x[ 0 ], x[ 1 ], x[ 2 ], ig ); // quadrature points coordinates
1657 
1658  // local coordinates of the quadrature point
1659  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1660  {
1661  sum1 = 0;
1662  sum2 = 0;
1663  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1664  {
1665  sum1 += fe1.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b1[ jcoor ] );
1666  sum2 += fe2.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b2[ jcoor ] );
1667  }
1668  rx1[ icoor ] = sum1;
1669  rx2[ icoor ] = sum2;
1670  }
1671 
1672  for ( i = 0; i < fe1.nbFEDof(); ++i )
1673  {
1674 
1675  // first derivative on the reference element
1676  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1677  {
1678  drp1[ icoor ] = fe1.refFE().dPhi ( i, icoor, rx1[ 0 ], rx1[ 1 ], rx1[ 2 ] );
1679  drp2[ icoor ] = fe2.refFE().dPhi ( i, icoor, rx2[ 0 ], rx2[ 1 ], rx2[ 2 ] );
1680  }
1681 
1682  // first derivative on the current element
1683  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1684  {
1685  sum1 = 0;
1686  sum2 = 0;
1687  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1688  {
1689  sum1 += fe1.tInvJac ( icoor, jcoor, 0 ) * drp1[ jcoor ];
1690  sum2 += fe2.tInvJac ( icoor, jcoor, 0 ) * drp2[ jcoor ];
1691  }
1692  phid1[ i ][ icoor ][ ig ] = sum1;
1693  phid2[ i ][ icoor ][ ig ] = sum2;
1694  }
1695  }
1696  }
1697 
1698 
1699 
1700  // Loop on rows
1701  for ( i = 0; i < fe1.nbFEDof(); ++i )
1702  {
1703  // Loop on columns
1704  for ( j = 0; j < fe2.nbFEDof(); ++j )
1705  {
1706  sum = 0.0;
1707  // Loop on coordinates
1708  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1709  for ( ig = 0; ig < bdfe.nbQuadPt() ; ++ig )
1710  sum += ba2[ ig ] *
1711  phid1[ i ][ icoor ][ ig ] *
1712  phid2[ j ][ icoor ][ ig ] *
1713  bdfe.wRootDetMetric ( ig );
1714  mat ( i, j ) += coef * sum;
1715  }
1716  }
1717 
1718 }
1719 
1720 
1721 // coef < |\beta . n| \grad p1, \grad q2 >
1722 // p1 lives in fe1
1723 // q2 lives in fe2
1724 // beta lives in fe3
1725 
1726 
1727 void ipstab_bagrad ( const Real coef,
1728  MatrixElemental& elmat,
1729  const CurrentFE& fe1,
1730  const CurrentFE& fe2,
1731  const CurrentFE& fe3,
1732  const VectorElemental& beta,
1733  const CurrentFEManifold& bdfe,
1734  int iblock, int jblock )
1735 {
1736 
1737  // Interior penalty stabilization:
1738  // coef < |\beta.n| \grad p1, \grad q2 >
1739 
1740  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
1741 
1742  Real sum, sum1, sum2;
1743  UInt icoor, jcoor;
1744  UInt i, j, ig;
1745 
1746  boost::multi_array<Real, 3> phid1 (
1747  boost::extents[fe1.nbFEDof()][fe1.nbLocalCoor()][bdfe.nbQuadPt()]);
1748  boost::multi_array<Real, 3> phid2 (
1749  boost::extents[fe2.nbFEDof()][fe2.nbLocalCoor()][bdfe.nbQuadPt()]);
1750 
1751  std::vector<Real> x (3), rx1 (3), drp1 (3), rx2 (3), drp2 (3);
1752  std::vector<Real> b1 (3), b2 (3);
1753 
1754  fe1.coorMap ( b1[ 0 ], b1[ 1 ], b1[ 2 ], 0, 0, 0 ); // translation fe1
1755  fe2.coorMap ( b2[ 0 ], b2[ 1 ], b2[ 2 ], 0, 0, 0 ); // translation fe2
1756 
1757  //
1758  // convection velocity term |\beta . n|
1759  // on the boundary quadrature points
1760  //
1761  std::vector<Real> bn (bdfe.nbQuadPt() );
1762 
1763  for ( ig = 0; ig < bdfe.nbQuadPt(); ig++ )
1764  {
1765  sum1 = 0;
1766  sum2 = 0;
1767 
1768  for ( icoor = 0; icoor < fe3.nbLocalCoor(); ++icoor )
1769  {
1770  for ( i = 0; i < fe3.nbFEDof(); ++i )
1771  {
1772  Real betaLoc = fe3.phi ( i, ig ) *
1773  beta.vec() [ icoor * fe3.nbFEDof() + i ];
1774  sum1 += betaLoc * bdfe.normal (icoor, ig);
1775  }
1776  }
1777  bn[ ig ] = std::fabs (sum1);
1778  }
1779 
1780  for ( UInt ig = 0; ig < bdfe.nbQuadPt(); ++ig )
1781  {
1782  // first derivatives on quadrature points
1783  bdfe.coorQuadPt ( x[ 0 ], x[ 1 ], x[ 2 ], ig ); // quadrature points coordinates
1784 
1785  // local coordinates of the quadrature point
1786  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1787  {
1788  sum1 = 0;
1789  sum2 = 0;
1790  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1791  {
1792  sum1 += fe1.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b1[ jcoor ] );
1793  sum2 += fe2.tInvJac ( jcoor, icoor, 0 ) * ( x[ jcoor ] - b2[ jcoor ] );
1794  }
1795  rx1[ icoor ] = sum1;
1796  rx2[ icoor ] = sum2;
1797  }
1798 
1799  for ( i = 0; i < fe1.nbFEDof(); ++i )
1800  {
1801 
1802  // first derivative on the reference element
1803  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1804  {
1805  drp1[ icoor ] = fe1.refFE().dPhi ( i, icoor, rx1[ 0 ], rx1[ 1 ], rx1[ 2 ] );
1806  drp2[ icoor ] = fe2.refFE().dPhi ( i, icoor, rx2[ 0 ], rx2[ 1 ], rx2[ 2 ] );
1807  }
1808 
1809  // first derivative on the current element
1810  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1811  {
1812  sum1 = 0;
1813  sum2 = 0;
1814  for ( jcoor = 0; jcoor < fe1.nbLocalCoor(); ++jcoor )
1815  {
1816  sum1 += fe1.tInvJac ( icoor, jcoor, 0 ) * drp1[ jcoor ];
1817  sum2 += fe2.tInvJac ( icoor, jcoor, 0 ) * drp2[ jcoor ];
1818  }
1819  phid1[ i ][ icoor ][ ig ] = sum1;
1820  phid2[ i ][ icoor ][ ig ] = sum2;
1821  }
1822  }
1823  }
1824 
1825 
1826 
1827  // Loop on rows
1828  for ( i = 0; i < fe1.nbFEDof(); ++i )
1829  {
1830  // Loop on columns
1831  for ( j = 0; j < fe2.nbFEDof(); ++j )
1832  {
1833  sum = 0.0;
1834  // Loop on coordinates
1835  for ( icoor = 0; icoor < fe1.nbLocalCoor(); ++icoor )
1836  for ( ig = 0; ig < bdfe.nbQuadPt() ; ++ig )
1837  sum += bn[ ig ] *
1838  phid1[ i ][ icoor ][ ig ] *
1839  phid2[ j ][ icoor ][ ig ] *
1840  bdfe.wRootDetMetric ( ig );
1841  mat ( i, j ) += coef * sum;
1842  }
1843  }
1844 
1845 }
1846 
1847 
1848 void stiff ( Real coef, MatrixElemental& elmat, const CurrentFE& fe,
1849  int iblock, int jblock )
1850 /*
1851  Stiffness matrix: coef*\int grad v_i . grad v_j
1852 */
1853 {
1854  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
1855  UInt iloc, jloc;
1856  UInt i, icoor, ig;
1857  double s, coef_s;
1858  //
1859  // diagonal
1860  //
1861  for ( i = 0; i < fe.nbDiag(); i++ )
1862  {
1863  iloc = fe.patternFirst ( i );
1864  s = 0;
1865  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
1866  {
1867  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
1868  s += fe.phiDer ( iloc, icoor, ig ) * fe.phiDer ( iloc, icoor, ig )
1869  * fe.weightDet ( ig );
1870  }
1871  mat ( iloc, iloc ) += coef * s;
1872  }
1873  //
1874  // extra diagonal
1875  //
1876  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
1877  {
1878  iloc = fe.patternFirst ( i );
1879  jloc = fe.patternSecond ( i );
1880  s = 0;
1881  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
1882  {
1883  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
1884  s += fe.phiDer ( iloc, icoor, ig ) * fe.phiDer ( jloc, icoor, ig ) *
1885  fe.weightDet ( ig );
1886  }
1887  coef_s = coef * s;
1888  mat ( iloc, jloc ) += coef_s;
1889  mat ( jloc, iloc ) += coef_s;
1890  }
1891 }
1892 
1893 
1894 void stiff ( Real coef, Real ( *fct ) ( Real, Real, Real ), MatrixElemental& elmat,
1895  const CurrentFE& fe, int iblock, int jblock )
1896 /*
1897  Stiffness matrix: coef*\int grad v_i . grad v_j
1898 */
1899 {
1900  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
1901  UInt iloc, jloc;
1902  UInt i, icoor, ig;
1903  double s, coef_s, coef_f;
1904  //
1905  // diagonal
1906  //
1907  for ( i = 0; i < fe.nbDiag(); i++ )
1908  {
1909  iloc = fe.patternFirst ( i );
1910  s = 0;
1911  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
1912  {
1913  coef_f = fct ( fe.quadPt ( ig, 0 ), fe.quadPt ( ig, 1 ), fe.quadPt ( ig, 2 ) );
1914  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
1915  s += coef_f * fe.phiDer ( iloc, icoor, ig ) * fe.phiDer ( iloc, icoor, ig )
1916  * fe.weightDet ( ig );
1917  }
1918  mat ( iloc, iloc ) += coef * s;
1919  }
1920  //
1921  // extra diagonal
1922  //
1923  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
1924  {
1925  iloc = fe.patternFirst ( i );
1926  jloc = fe.patternSecond ( i );
1927  s = 0;
1928  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
1929  {
1930  coef_f = fct ( fe.quadPt ( ig, 0 ), fe.quadPt ( ig, 1 ), fe.quadPt ( ig, 2 ) );
1931  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
1932  s += coef_f * fe.phiDer ( iloc, icoor, ig ) * fe.phiDer ( jloc, icoor, ig ) *
1933  fe.weightDet ( ig );
1934  }
1935  coef_s = coef * s;
1936  mat ( iloc, jloc ) += coef_s;
1937  mat ( jloc, iloc ) += coef_s;
1938  }
1939 }
1940 //
1941 
1942 void stiff ( Real coef, MatrixElemental& elmat, const CurrentFE& fe,
1943  int iblock, int jblock, int nb )
1944 /*
1945  Stiffness matrix: coef*\int grad v_i . grad v_j (nb blocks on the diagonal, nb>1)
1946 */
1947 {
1948 
1949 
1950  Matrix mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
1951  mat_tmp = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
1952 
1953  UInt iloc, jloc;
1954  UInt i, icoor, ig;
1955  double s, coef_s;
1956  //
1957  // diagonal
1958  //
1959  for ( i = 0; i < fe.nbDiag(); i++ )
1960  {
1961  iloc = fe.patternFirst ( i );
1962  s = 0;
1963  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
1964  {
1965  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
1966  s += fe.phiDer ( iloc, icoor, ig ) * fe.phiDer ( iloc, icoor, ig )
1967  * fe.weightDet ( ig );
1968  }
1969  mat_tmp ( iloc, iloc ) += coef * s;
1970  }
1971  //
1972  // extra diagonal
1973  //
1974  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
1975  {
1976  iloc = fe.patternFirst ( i );
1977  jloc = fe.patternSecond ( i );
1978  s = 0;
1979  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
1980  {
1981  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
1982  s += fe.phiDer ( iloc, icoor, ig ) * fe.phiDer ( jloc, icoor, ig ) *
1983  fe.weightDet ( ig );
1984  }
1985  coef_s = coef * s;
1986  mat_tmp ( iloc, jloc ) += coef_s;
1987  mat_tmp ( jloc, iloc ) += coef_s;
1988  }
1989  // copy on the components
1990  for ( int icomp = 0; icomp < nb; icomp++ )
1991  {
1992  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icomp, jblock + icomp );
1993  for ( i = 0; i < fe.nbDiag(); i++ )
1994  {
1995  iloc = fe.patternFirst ( i );
1996  mat_icomp ( iloc, iloc ) += mat_tmp ( iloc, iloc );
1997  }
1998  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
1999  {
2000  iloc = fe.patternFirst ( i );
2001  jloc = fe.patternSecond ( i );
2002  mat_icomp ( iloc, jloc ) += mat_tmp ( iloc, jloc );
2003  mat_icomp ( jloc, iloc ) += mat_tmp ( jloc, iloc );
2004  }
2005  }
2006 }
2007 
2008 
2009 void stiff ( const std::vector<Real>& coef, MatrixElemental& elmat, const CurrentFE& fe,
2010  int iblock, int jblock, int nb )
2011 /*
2012  Stiffness matrix: coef*\int grad v_i . grad v_j (nb blocks on the diagonal, nb>1)
2013  with coef given in a vector (one element per quadrature point)
2014 */
2015 {
2016 
2017  Matrix mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
2018  mat_tmp = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2019 
2020  UInt iloc, jloc;
2021  UInt i, icoor, ig;
2022  double s;//, coef_s;
2023  //
2024  // diagonal
2025  //
2026  for ( i = 0; i < fe.nbDiag(); i++ )
2027  {
2028  iloc = fe.patternFirst ( i );
2029  s = 0;
2030  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2031  {
2032  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
2033  s += coef[ig] * fe.phiDer ( iloc, icoor, ig ) * fe.phiDer ( iloc, icoor, ig )
2034  * fe.weightDet ( ig );
2035  }
2036  mat_tmp ( iloc, iloc ) += s;
2037  }
2038  //
2039  // extra diagonal
2040  //
2041  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2042  {
2043  iloc = fe.patternFirst ( i );
2044  jloc = fe.patternSecond ( i );
2045  s = 0;
2046  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2047  {
2048  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
2049  s += coef[ig] * fe.phiDer ( iloc, icoor, ig ) * fe.phiDer ( jloc, icoor, ig ) *
2050  fe.weightDet ( ig );
2051  }
2052  mat_tmp ( iloc, jloc ) += s;
2053  mat_tmp ( jloc, iloc ) += s;
2054  }
2055  // copy on the components
2056  for ( int icomp = 0; icomp < nb; icomp++ )
2057  {
2058  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icomp, jblock + icomp );
2059  for ( i = 0; i < fe.nbDiag(); i++ )
2060  {
2061  iloc = fe.patternFirst ( i );
2062  mat_icomp ( iloc, iloc ) += mat_tmp ( iloc, iloc );
2063  }
2064  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2065  {
2066  iloc = fe.patternFirst ( i );
2067  jloc = fe.patternSecond ( i );
2068  mat_icomp ( iloc, jloc ) += mat_tmp ( iloc, jloc );
2069  mat_icomp ( jloc, iloc ) += mat_tmp ( jloc, iloc );
2070  }
2071  }
2072 }
2073 
2074 
2075 // VC - December 2004
2076 //
2077 void stiff_curl ( Real coef, MatrixElemental& elmat, const CurrentFE& fe,
2078  int iblock, int jblock, int /*nb*/ )
2079 
2080 
2081 /*
2082  Stiffness matrix: coef*\int curl v_i . curl v_j (nb blocks on the diagonal, nb>1)
2083 */
2084 {
2085 
2086 
2087  Matrix mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
2088  mat_tmp = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2089  Matrix mat_tmp11 ( fe.nbFEDof(), fe.nbFEDof() );
2090  mat_tmp11 = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2091  Matrix mat_tmp12 ( fe.nbFEDof(), fe.nbFEDof() );
2092  mat_tmp12 = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2093  Matrix mat_tmp13 ( fe.nbFEDof(), fe.nbFEDof() );
2094  mat_tmp13 = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2095  Matrix mat_tmp21 ( fe.nbFEDof(), fe.nbFEDof() );
2096  mat_tmp21 = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2097  Matrix mat_tmp22 ( fe.nbFEDof(), fe.nbFEDof() );
2098  mat_tmp22 = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2099  Matrix mat_tmp23 ( fe.nbFEDof(), fe.nbFEDof() );
2100  mat_tmp23 = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2101  Matrix mat_tmp31 ( fe.nbFEDof(), fe.nbFEDof() );
2102  mat_tmp31 = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2103  Matrix mat_tmp32 ( fe.nbFEDof(), fe.nbFEDof() );
2104  mat_tmp32 = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2105  Matrix mat_tmp33 ( fe.nbFEDof(), fe.nbFEDof() );
2106  mat_tmp33 = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2107 
2108 
2109 
2110  UInt iloc, jloc;
2111  UInt i, ig;
2112  double s, coef_s;
2113 
2114  // diagonal 11
2115  //
2116  for ( i = 0; i < fe.nbDiag(); i++ )
2117  {
2118  iloc = fe.patternFirst ( i );
2119  s = 0;
2120  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2121  {
2122  s = fe.phiDer ( iloc, 1, ig ) * fe.phiDer ( iloc, 1, ig ) * fe.weightDet ( ig )
2123  + fe.phiDer ( iloc, 2, ig ) * fe.phiDer ( iloc, 2, ig ) * fe.weightDet ( ig ) ;
2124  }
2125  mat_tmp11 ( iloc, iloc ) += coef * s;
2126  }
2127  // extra diagonal 11
2128  //
2129  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2130  {
2131  iloc = fe.patternFirst ( i );
2132  jloc = fe.patternSecond ( i );
2133  s = 0;
2134  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2135  {
2136  s = fe.phiDer ( iloc, 1, ig ) * fe.phiDer ( jloc, 1, ig ) * fe.weightDet ( ig )
2137  + fe.phiDer ( iloc, 2, ig ) * fe.phiDer ( jloc, 2, ig ) * fe.weightDet ( ig ) ;
2138  }
2139  coef_s = coef * s;
2140  mat_tmp11 ( iloc, jloc ) += coef_s;
2141  mat_tmp11 ( jloc, iloc ) += coef_s;
2142  }
2143 
2144  // diagonal 12
2145  //
2146  for ( i = 0; i < fe.nbDiag(); i++ )
2147  {
2148  iloc = fe.patternFirst ( i );
2149  s = 0;
2150  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2151  {
2152  s = - fe.phiDer ( iloc, 1, ig ) * fe.phiDer ( iloc, 0, ig ) * fe.weightDet ( ig );
2153  }
2154  mat_tmp12 ( iloc, iloc ) -= coef * s;
2155  }
2156  // extra diagonal 12
2157  //
2158  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2159  {
2160  iloc = fe.patternFirst ( i );
2161  jloc = fe.patternSecond ( i );
2162  s = 0;
2163  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2164  {
2165  s = - fe.phiDer ( iloc, 1, ig ) * fe.phiDer ( jloc, 0, ig ) * fe.weightDet ( ig );
2166  }
2167  coef_s = coef * s;
2168  mat_tmp12 ( iloc, jloc ) -= coef_s;
2169  mat_tmp12 ( jloc, iloc ) -= coef_s;
2170  }
2171 
2172  // diagonal 13
2173  //
2174  for ( i = 0; i < fe.nbDiag(); i++ )
2175  {
2176  iloc = fe.patternFirst ( i );
2177  s = 0;
2178  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2179  {
2180  s = - fe.phiDer ( iloc, 2, ig ) * fe.phiDer ( iloc, 0, ig ) * fe.weightDet ( ig );
2181  }
2182  mat_tmp13 ( iloc, iloc ) -= coef * s;
2183  }
2184  // extra diagonal 13
2185  //
2186  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2187  {
2188  iloc = fe.patternFirst ( i );
2189  jloc = fe.patternSecond ( i );
2190  s = 0;
2191  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2192  {
2193  s = - fe.phiDer ( iloc, 2, ig ) * fe.phiDer ( jloc, 0, ig ) * fe.weightDet ( ig );
2194  }
2195  coef_s = coef * s;
2196  mat_tmp13 ( iloc, jloc ) -= coef_s;
2197  mat_tmp13 ( jloc, iloc ) -= coef_s;
2198  }
2199 
2200  // diagonal 21
2201  //
2202  for ( i = 0; i < fe.nbDiag(); i++ )
2203  {
2204  iloc = fe.patternFirst ( i );
2205  s = 0;
2206  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2207  {
2208  s = - fe.phiDer ( iloc, 0, ig ) * fe.phiDer ( iloc, 1, ig ) * fe.weightDet ( ig );
2209  }
2210  mat_tmp21 ( iloc, iloc ) -= coef * s;
2211  }
2212  // extra diagonal 21
2213  //
2214  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2215  {
2216  iloc = fe.patternFirst ( i );
2217  jloc = fe.patternSecond ( i );
2218  s = 0;
2219  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2220  {
2221  s = - fe.phiDer ( iloc, 0, ig ) * fe.phiDer ( jloc, 1, ig ) * fe.weightDet ( ig );
2222  }
2223  coef_s = coef * s;
2224  mat_tmp21 ( iloc, jloc ) -= coef_s;
2225  mat_tmp21 ( jloc, iloc ) -= coef_s;
2226  }
2227 
2228  // diagonal 22
2229  //
2230  for ( i = 0; i < fe.nbDiag(); i++ )
2231  {
2232  iloc = fe.patternFirst ( i );
2233  s = 0;
2234  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2235  {
2236  s = fe.phiDer ( iloc, 0, ig ) * fe.phiDer ( iloc, 0, ig ) * fe.weightDet ( ig )
2237  + fe.phiDer ( iloc, 2, ig ) * fe.phiDer ( iloc, 2, ig ) * fe.weightDet ( ig ) ;
2238  }
2239  mat_tmp22 ( iloc, iloc ) += coef * s;
2240  }
2241  // extra diagonal 22
2242  //
2243  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2244  {
2245  iloc = fe.patternFirst ( i );
2246  jloc = fe.patternSecond ( i );
2247  s = 0;
2248  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2249  {
2250  s = fe.phiDer ( iloc, 0, ig ) * fe.phiDer ( jloc, 0, ig ) * fe.weightDet ( ig )
2251  + fe.phiDer ( iloc, 2, ig ) * fe.phiDer ( jloc, 2, ig ) * fe.weightDet ( ig ) ;
2252  }
2253  coef_s = coef * s;
2254  mat_tmp22 ( iloc, jloc ) += coef_s;
2255  mat_tmp22 ( jloc, iloc ) += coef_s;
2256  }
2257 
2258  // diagonal 23
2259  //
2260  for ( i = 0; i < fe.nbDiag(); i++ )
2261  {
2262  iloc = fe.patternFirst ( i );
2263  s = 0;
2264  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2265  {
2266  s = - fe.phiDer ( iloc, 2, ig ) * fe.phiDer ( iloc, 1, ig ) * fe.weightDet ( ig );
2267  }
2268  mat_tmp23 ( iloc, iloc ) -= coef * s;
2269  }
2270  // extra diagonal 23
2271  //
2272  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2273  {
2274  iloc = fe.patternFirst ( i );
2275  jloc = fe.patternSecond ( i );
2276  s = 0;
2277  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2278  {
2279  s = - fe.phiDer ( iloc, 2, ig ) * fe.phiDer ( jloc, 1, ig ) * fe.weightDet ( ig );
2280  }
2281  coef_s = coef * s;
2282  mat_tmp23 ( iloc, jloc ) -= coef_s;
2283  mat_tmp23 ( jloc, iloc ) -= coef_s;
2284  }
2285 
2286  // diagonal 31
2287  //
2288  for ( i = 0; i < fe.nbDiag(); i++ )
2289  {
2290  iloc = fe.patternFirst ( i );
2291  s = 0;
2292  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2293  {
2294  s = - fe.phiDer ( iloc, 0, ig ) * fe.phiDer ( iloc, 2, ig ) * fe.weightDet ( ig );
2295  }
2296  mat_tmp31 ( iloc, iloc ) -= coef * s;
2297  }
2298  // extra diagonal 31
2299  //
2300  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2301  {
2302  iloc = fe.patternFirst ( i );
2303  jloc = fe.patternSecond ( i );
2304  s = 0;
2305  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2306  {
2307  s = - fe.phiDer ( iloc, 0, ig ) * fe.phiDer ( jloc, 2, ig ) * fe.weightDet ( ig );
2308  }
2309  coef_s = coef * s;
2310  mat_tmp31 ( iloc, jloc ) -= coef_s;
2311  mat_tmp31 ( jloc, iloc ) -= coef_s;
2312  }
2313 
2314  // diagonal 32
2315  //
2316  for ( i = 0; i < fe.nbDiag(); i++ )
2317  {
2318  iloc = fe.patternFirst ( i );
2319  s = 0;
2320  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2321  {
2322  s = - fe.phiDer ( iloc, 1, ig ) * fe.phiDer ( iloc, 2, ig ) * fe.weightDet ( ig );
2323  }
2324  mat_tmp32 ( iloc, iloc ) -= coef * s;
2325  }
2326  // extra diagonal 32
2327  //
2328  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2329  {
2330  iloc = fe.patternFirst ( i );
2331  jloc = fe.patternSecond ( i );
2332  s = 0;
2333  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2334  {
2335  s = - fe.phiDer ( iloc, 1, ig ) * fe.phiDer ( jloc, 2, ig ) * fe.weightDet ( ig );
2336  }
2337  coef_s = coef * s;
2338  mat_tmp32 ( iloc, jloc ) -= coef_s;
2339  mat_tmp32 ( jloc, iloc ) -= coef_s;
2340  }
2341 
2342  // diagonal 33
2343  //
2344  for ( i = 0; i < fe.nbDiag(); i++ )
2345  {
2346  iloc = fe.patternFirst ( i );
2347  s = 0;
2348  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2349  {
2350  s = fe.phiDer ( iloc, 0, ig ) * fe.phiDer ( iloc, 0, ig ) * fe.weightDet ( ig )
2351  + fe.phiDer ( iloc, 1, ig ) * fe.phiDer ( iloc, 1, ig ) * fe.weightDet ( ig ) ;
2352  }
2353  mat_tmp33 ( iloc, iloc ) += coef * s;
2354  }
2355  // extra diagonal 33
2356  //
2357  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2358  {
2359  iloc = fe.patternFirst ( i );
2360  jloc = fe.patternSecond ( i );
2361  s = 0;
2362  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2363  {
2364  s = fe.phiDer ( iloc, 1, ig ) * fe.phiDer ( jloc, 1, ig ) * fe.weightDet ( ig )
2365  + fe.phiDer ( iloc, 2, ig ) * fe.phiDer ( jloc, 2, ig ) * fe.weightDet ( ig ) ;
2366  }
2367  coef_s = coef * s;
2368  mat_tmp33 ( iloc, jloc ) += coef_s;
2369  mat_tmp33 ( jloc, iloc ) += coef_s;
2370  }
2371 
2372  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + 0, jblock + 0 );
2373  for ( i = 0; i < fe.nbDiag(); i++ )
2374  {
2375  iloc = fe.patternFirst ( i );
2376  mat_icomp ( iloc, iloc ) += mat_tmp11 ( iloc, iloc );
2377  }
2378  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2379  {
2380  iloc = fe.patternFirst ( i );
2381  jloc = fe.patternSecond ( i );
2382  mat_icomp ( iloc, jloc ) += mat_tmp11 ( iloc, jloc );
2383  mat_icomp ( jloc, iloc ) += mat_tmp11 ( jloc, iloc );
2384  }
2385 
2386  mat_icomp = elmat.block ( iblock + 0, jblock + 1 );
2387  for ( i = 0; i < fe.nbDiag(); i++ )
2388  {
2389  iloc = fe.patternFirst ( i );
2390  mat_icomp ( iloc, iloc ) -= mat_tmp12 ( iloc, iloc );
2391  }
2392  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2393  {
2394  iloc = fe.patternFirst ( i );
2395  jloc = fe.patternSecond ( i );
2396  mat_icomp ( iloc, jloc ) -= mat_tmp12 ( iloc, jloc );
2397  mat_icomp ( jloc, iloc ) -= mat_tmp12 ( jloc, iloc );
2398  }
2399 
2400  mat_icomp = elmat.block ( iblock + 0, jblock + 2 );
2401  for ( i = 0; i < fe.nbDiag(); i++ )
2402  {
2403  iloc = fe.patternFirst ( i );
2404  mat_icomp ( iloc, iloc ) -= mat_tmp13 ( iloc, iloc );
2405  }
2406  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2407  {
2408  iloc = fe.patternFirst ( i );
2409  jloc = fe.patternSecond ( i );
2410  mat_icomp ( iloc, jloc ) -= mat_tmp13 ( iloc, jloc );
2411  mat_icomp ( jloc, iloc ) -= mat_tmp13 ( jloc, iloc );
2412  }
2413 
2414  mat_icomp = elmat.block ( iblock + 1, jblock + 0 );
2415  for ( i = 0; i < fe.nbDiag(); i++ )
2416  {
2417  iloc = fe.patternFirst ( i );
2418  mat_icomp ( iloc, iloc ) -= mat_tmp21 ( iloc, iloc );
2419  }
2420  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2421  {
2422  iloc = fe.patternFirst ( i );
2423  jloc = fe.patternSecond ( i );
2424  mat_icomp ( iloc, jloc ) -= mat_tmp21 ( iloc, jloc );
2425  mat_icomp ( jloc, iloc ) -= mat_tmp21 ( jloc, iloc );
2426  }
2427 
2428  mat_icomp = elmat.block ( iblock + 1, jblock + 1 );
2429  for ( i = 0; i < fe.nbDiag(); i++ )
2430  {
2431  iloc = fe.patternFirst ( i );
2432  mat_icomp ( iloc, iloc ) += mat_tmp22 ( iloc, iloc );
2433  }
2434  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2435  {
2436  iloc = fe.patternFirst ( i );
2437  jloc = fe.patternSecond ( i );
2438  mat_icomp ( iloc, jloc ) += mat_tmp22 ( iloc, jloc );
2439  mat_icomp ( jloc, iloc ) += mat_tmp22 ( jloc, iloc );
2440  }
2441 
2442  mat_icomp = elmat.block ( iblock + 1, jblock + 2 );
2443  for ( i = 0; i < fe.nbDiag(); i++ )
2444  {
2445  iloc = fe.patternFirst ( i );
2446  mat_icomp ( iloc, iloc ) -= mat_tmp23 ( iloc, iloc );
2447  }
2448  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2449  {
2450  iloc = fe.patternFirst ( i );
2451  jloc = fe.patternSecond ( i );
2452  mat_icomp ( iloc, jloc ) -= mat_tmp23 ( iloc, jloc );
2453  mat_icomp ( jloc, iloc ) -= mat_tmp23 ( jloc, iloc );
2454  }
2455 
2456  mat_icomp = elmat.block ( iblock + 2, jblock + 0 );
2457  for ( i = 0; i < fe.nbDiag(); i++ )
2458  {
2459  iloc = fe.patternFirst ( i );
2460  mat_icomp ( iloc, iloc ) -= mat_tmp31 ( iloc, iloc );
2461  }
2462  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2463  {
2464  iloc = fe.patternFirst ( i );
2465  jloc = fe.patternSecond ( i );
2466  mat_icomp ( iloc, jloc ) -= mat_tmp31 ( iloc, jloc );
2467  mat_icomp ( jloc, iloc ) -= mat_tmp31 ( jloc, iloc );
2468  }
2469 
2470  mat_icomp = elmat.block ( iblock + 2, jblock + 1 );
2471  for ( i = 0; i < fe.nbDiag(); i++ )
2472  {
2473  iloc = fe.patternFirst ( i );
2474  mat_icomp ( iloc, iloc ) -= mat_tmp32 ( iloc, iloc );
2475  }
2476  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2477  {
2478  iloc = fe.patternFirst ( i );
2479  jloc = fe.patternSecond ( i );
2480  mat_icomp ( iloc, jloc ) -= mat_tmp32 ( iloc, jloc );
2481  mat_icomp ( jloc, iloc ) -= mat_tmp32 ( jloc, iloc );
2482  }
2483 
2484  mat_icomp = elmat.block ( iblock + 2, jblock + 2 );
2485  for ( i = 0; i < fe.nbDiag(); i++ )
2486  {
2487  iloc = fe.patternFirst ( i );
2488  mat_icomp ( iloc, iloc ) += mat_tmp33 ( iloc, iloc );
2489  }
2490  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2491  {
2492  iloc = fe.patternFirst ( i );
2493  jloc = fe.patternSecond ( i );
2494  mat_icomp ( iloc, jloc ) += mat_tmp33 ( iloc, jloc );
2495  mat_icomp ( jloc, iloc ) += mat_tmp ( jloc, iloc );
2496  }
2497 }
2498 
2499 
2500 /*
2501  Stiffness matrix: coef * ( div u , div v )
2502 */
2503 void stiff_div ( Real coef, MatrixElemental& elmat, const CurrentFE& fe )
2504 {
2505 
2506  double s;
2507 
2508  //
2509  // blocks (icoor,jcoor) of elmat
2510  //
2511  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2512  {
2513  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
2514  {
2515 
2516  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2517 
2518  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2519  {
2520  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2521  {
2522  s = 0;
2523  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2524  {
2525  s += fe.phiDer ( i, icoor, ig ) * fe.phiDer ( j, jcoor, ig ) * fe.weightDet ( ig );
2526  }
2527  mat ( i, j ) += coef * s;
2528  }
2529  }
2530  }
2531  }
2532 }
2533 
2534 
2535 
2536 /*
2537  Stiffness matrix: coef * ( [\grad u^k]^T \grad d : \grad v )
2538 */
2539 void stiff_dergradbis ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
2540 {
2541 
2542  double s;
2543  boost::multi_array<Real, 3> guk (
2544  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
2545 
2546  // loop on quadrature points
2547  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
2548  {
2549 
2550  // loop on space coordinates
2551  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
2552  {
2553 
2554  // loop on space coordinates
2555  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
2556  {
2557  s = 0.0;
2558  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
2559  {
2560  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
2561  }
2562  guk[ icoor ][ jcoor ][ ig ] = s;
2563  }
2564  }
2565  }
2566  //
2567  // blocks (icoor,jcoor) of elmat
2568  //
2569 
2570  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2571  {
2572  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
2573  {
2574 
2575  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2576 
2577  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2578  {
2579  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2580  {
2581  s = 0;
2582  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
2583  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2584  {
2585  s += fe.phiDer ( i, k, ig ) * guk[ jcoor ][ icoor ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2586  }
2587  mat ( i, j ) += coef * s;
2588  }
2589  }
2590  }
2591  }
2592 }
2593 
2594 
2595 
2596 
2597 /*
2598  Stiffness matrix: coef * ( [\grad u]^T \grad u^k [\grad u^k]^T \grad u : \grad v ) for Newton on St-Venant
2599 */
2600 void stiff_dergrad ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
2601 {
2602 
2603  double s;
2604 
2605  boost::multi_array<Real, 3> guk (
2606  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
2607 
2608  // loop on quadrature points
2609  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
2610  {
2611 
2612  // loop on space coordinates
2613  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
2614  {
2615 
2616  // loop on space coordinates
2617  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
2618  {
2619  s = 0.0;
2620  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
2621  {
2622  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
2623  }
2624  guk[ icoor ][ jcoor ][ ig ] = s;
2625  }
2626  }
2627  }
2628  //
2629  // blocks (icoor,jcoor) of elmat
2630  //
2631 
2632  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2633  {
2634  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
2635  {
2636 
2637  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2638 
2639  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2640  {
2641  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2642  {
2643  s = 0;
2644  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
2645  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2646  {
2647  s += fe.phiDer ( i, k, ig ) * ( guk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, icoor, ig )
2648  + guk[ jcoor ][ icoor ][ ig ] * fe.phiDer ( j, k, ig ) ) * fe.weightDet ( ig );
2649  }
2650  mat ( i, j ) += coef * s;
2651  }
2652  }
2653  }
2654  }
2655 }
2656 
2657 
2658 
2659 
2660 
2661 //
2662 // coef * ( \tr { [\grad u^k]^T \grad u }, \div v ) for Newton on St-Venant
2663 //
2664 //
2665 void stiff_derdiv ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
2666 {
2667  boost::multi_array<Real, 3> guk (
2668  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
2669 
2670  Real s;
2671 
2672  // loop on quadrature points
2673  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
2674  {
2675 
2676  // loop on space coordinates
2677  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
2678  {
2679 
2680  // loop on space coordinates
2681  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
2682  {
2683  s = 0.0;
2684  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
2685  {
2686  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
2687  }
2688  guk[ icoor ][ jcoor ][ ig ] = s;
2689  }
2690  }
2691  }
2692  //
2693  // blocks (icoor,jcoor) of elmat
2694  //
2695  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2696  {
2697  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
2698  {
2699 
2700  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2701 
2702  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2703  {
2704  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2705  {
2706  s = 0;
2707  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
2708  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
2709  {
2710  s += fe.phiDer ( i, icoor, ig ) * guk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2711  }
2712  mat ( i, j ) += coef * s;
2713  }
2714  }
2715  }
2716  }
2717 }
2718 
2719 
2720 
2721 
2722 void stiff_strain ( Real coef, MatrixElemental& elmat, const CurrentFE& fe )
2723 /*
2724  Stiffness matrix: coef * ( e(u) , e(v) )
2725 */
2726 {
2727  double s;
2728  double tmp = coef * 0.5;
2729 
2730  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
2731 
2732  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2733  {
2734  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2735  {
2736  s = 0;
2737  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2738  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2739  {
2740  s += fe.phiDer ( i, icoor, ig ) * fe.phiDer ( j, icoor, ig ) * fe.weightDet ( ig );
2741  }
2742  mat_tmp ( i, j ) = tmp * s;
2743  }
2744  }
2745  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2746  {
2747  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
2748  mat += mat_tmp;
2749  }
2750 
2751  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2752  {
2753  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
2754  {
2755  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2756  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2757  {
2758  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2759  {
2760  s = 0;
2761  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2762  {
2763  s += fe.phiDer ( i, jcoor, ig ) * fe.phiDer ( j, icoor, ig ) * fe.weightDet ( ig );
2764  }
2765  mat ( i, j ) += tmp * s;
2766  }
2767  }
2768  }
2769  }
2770 }
2771 
2772 void mass_divw ( Real coef, const VectorElemental& w_loc, MatrixElemental& elmat, const CurrentFE& fe,
2773  int iblock, int jblock, UInt nb )
2774 /*
2775  modified mass matrix: ( div w u,v )
2776 */
2777 {
2778 
2779  Matrix mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
2780  mat_tmp = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2781 
2782  UInt i, icomp, ig, icoor, iloc, jloc;
2783  Real s, coef_s;
2784  std::vector<Real> divw (fe.nbQuadPt() );
2785 
2786  // divw at quadrature nodes
2787  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2788  {
2789  divw[ ig ] = 0.0;
2790  for ( icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2791  for ( i = 0; i < fe.nbFEDof(); ++i )
2792  {
2793  divw[ ig ] += fe.phiDer ( i, icoor, ig ) * w_loc.vec() [ i + icoor * fe.nbFEDof() ];
2794  }
2795  }
2796 
2797  //
2798  // diagonal
2799  //
2800  for ( i = 0; i < fe.nbDiag(); i++ )
2801  {
2802  iloc = fe.patternFirst ( i );
2803  s = 0;
2804  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2805  {
2806  s += divw[ ig ] * fe.phi ( iloc, ig ) * fe.phi ( iloc, ig ) * fe.weightDet ( ig );
2807  }
2808  mat_tmp ( iloc, iloc ) += coef * s;
2809  }
2810  //
2811  // extra diagonal
2812  //
2813  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2814  {
2815  iloc = fe.patternFirst ( i );
2816  jloc = fe.patternSecond ( i );
2817  s = 0;
2818  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2819  {
2820  s += divw[ ig ] * fe.phi ( iloc, ig ) * fe.phi ( jloc, ig ) * fe.weightDet ( ig );
2821  }
2822  coef_s = coef * s;
2823  mat_tmp ( iloc, jloc ) += coef_s;
2824  mat_tmp ( jloc, iloc ) += coef_s;
2825  }
2826  // copy on the components
2827  for ( icomp = 0; icomp < nb; icomp++ )
2828  {
2829  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icomp, jblock + icomp );
2830  for ( i = 0; i < fe.nbDiag(); i++ )
2831  {
2832  iloc = fe.patternFirst ( i );
2833  mat_icomp ( iloc, iloc ) += mat_tmp ( iloc, iloc );
2834  }
2835  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2836  {
2837  iloc = fe.patternFirst ( i );
2838  jloc = fe.patternSecond ( i );
2839  mat_icomp ( iloc, jloc ) += mat_tmp ( iloc, jloc );
2840  mat_icomp ( jloc, iloc ) += mat_tmp ( jloc, iloc );
2841  }
2842  }
2843 }
2844 
2845 
2846 void mass_divw (const std::vector<Real>& coef, const VectorElemental& w_loc, MatrixElemental& elmat, const CurrentFE& fe,
2847  int iblock, int jblock, UInt nb )
2848 /*
2849  modified mass matrix: ( div w u,v )
2850 */
2851 {
2852 
2853  Matrix mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
2854  mat_tmp = ZeroMatrix ( fe.nbFEDof(), fe.nbFEDof() );
2855 
2856  UInt i, icomp, ig, icoor, iloc, jloc;
2857  Real s;
2858  std::vector<Real> divw (fe.nbQuadPt() );
2859 
2860  // divw at quadrature nodes
2861  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2862  {
2863  divw[ ig ] = 0.0;
2864  for ( icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2865  for ( i = 0; i < fe.nbFEDof(); ++i )
2866  {
2867  divw[ ig ] += fe.phiDer ( i, icoor, ig ) * w_loc.vec() [ i + icoor * fe.nbFEDof() ];
2868  }
2869  }
2870 
2871  //
2872  // diagonal
2873  //
2874  for ( i = 0; i < fe.nbDiag(); i++ )
2875  {
2876  iloc = fe.patternFirst ( i );
2877  s = 0;
2878  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2879  {
2880  s += coef[ig] * divw[ ig ] * fe.phi ( iloc, ig ) * fe.phi ( iloc, ig ) * fe.weightDet ( ig );
2881  }
2882  mat_tmp ( iloc, iloc ) += s;
2883  }
2884  //
2885  // extra diagonal
2886  //
2887  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2888  {
2889  iloc = fe.patternFirst ( i );
2890  jloc = fe.patternSecond ( i );
2891  s = 0;
2892  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2893  {
2894  s += coef[ig] * divw[ ig ] * fe.phi ( iloc, ig ) * fe.phi ( jloc, ig ) * fe.weightDet ( ig );
2895  }
2896  mat_tmp ( iloc, jloc ) += s;
2897  mat_tmp ( jloc, iloc ) += s;
2898  }
2899  // copy on the components
2900  for ( icomp = 0; icomp < nb; icomp++ )
2901  {
2902  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icomp, jblock + icomp );
2903  for ( i = 0; i < fe.nbDiag(); i++ )
2904  {
2905  iloc = fe.patternFirst ( i );
2906  mat_icomp ( iloc, iloc ) += mat_tmp ( iloc, iloc );
2907  }
2908  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
2909  {
2910  iloc = fe.patternFirst ( i );
2911  jloc = fe.patternSecond ( i );
2912  mat_icomp ( iloc, jloc ) += mat_tmp ( iloc, jloc );
2913  mat_icomp ( jloc, iloc ) += mat_tmp ( jloc, iloc );
2914  }
2915  }
2916 }
2917 
2918 
2919 
2920 
2921 void mass_gradu ( Real coef, const VectorElemental& u0_loc, MatrixElemental& elmat, const CurrentFE& fe )
2922 /*
2923  modified mass matrix: ( grad u0 u,v )
2924 */
2925 {
2926 
2927  UInt ig, icoor, jcoor, i, j;
2928  Real s;
2929 
2930  boost::multi_array<Real, 3> gu0 (
2931  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()][fe.nbLocalCoor()]);
2932 
2933  //
2934  // grad u0 at quadrature nodes
2935  //
2936  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
2937  {
2938  for ( icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2939  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
2940  {
2941  gu0[ ig ][ icoor ][ jcoor ] = 0.0;
2942  for ( i = 0; i < fe.nbFEDof(); ++i )
2943  {
2944  gu0[ ig ][ icoor ][ jcoor ] += fe.phiDer ( i, jcoor, ig ) * u0_loc.vec() [ i + icoor * fe.nbFEDof() ];
2945  }
2946  }
2947  }
2948  //
2949  // blocks (icoor,jcoor) of elmat
2950  //
2951  for ( icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
2952  {
2953  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
2954  {
2955 
2956  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2957 
2958  for ( i = 0; i < fe.nbFEDof(); ++i )
2959  {
2960  for ( j = 0; j < fe.nbFEDof(); ++j )
2961  {
2962  s = 0;
2963  for ( ig = 0; ig < fe.nbQuadPt(); ++ig )
2964  {
2965  s += gu0[ ig ][ icoor ][ jcoor ] * fe.phi ( i, ig ) * fe.phi ( j, ig ) * fe.weightDet ( ig );
2966  }
2967  mat ( i, j ) += coef * s;
2968  }
2969  }
2970  }
2971  }
2972 }
2973 
2974 
2975 //
2976 //
2977 // \! Streamline diffusion
2978 //
2979 //
2980 void stiff_sd ( Real coef, const VectorElemental& vec_loc, MatrixElemental& elmat, const CurrentFE& fe, const CurrentFE& fe2,
2981  int iblock, int jblock, int nb )
2982 /*
2983  Stiffness matrix for SD Stabilization
2984 */
2985 {
2986  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
2987  UInt iloc, jloc;
2988  UInt i, icoor, ig, jcoor;
2989  std::vector<Real> coef_v (fe.nbLocalCoor() );
2990 
2991  double s, coef_s;
2992  // int nbN1=fe.nbFEDof();
2993  UInt nbN2 = fe2.nbFEDof();
2994  //
2995  // diagonal
2996  //
2997  for ( i = 0; i < fe.nbDiag(); i++ )
2998  {
2999  iloc = fe.patternFirst ( i );
3000  s = 0;
3001  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3002  {
3003  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3004  {
3005  coef_v[ icoor ] = 0.;
3006  }
3007 
3008  // computation of the convection term in the quadrature nodes
3009  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3010  {
3011  for ( UInt iloc = 0; iloc < nbN2; iloc++ )
3012  {
3013  coef_v[ icoor ] += vec_loc.vec() [ iloc + icoor * nbN2 ] * fe2.phi ( iloc, ig );
3014  }
3015  }
3016 
3017  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3018  {
3019  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3020  {
3021  s += coef_v[ icoor ] * fe.phiDer ( iloc, icoor, ig ) * coef_v[ jcoor ] * fe.phiDer ( iloc, jcoor, ig )
3022  * fe.weightDet ( ig );
3023  }
3024  }
3025  }
3026  mat ( iloc, iloc ) += coef * s;
3027  }
3028  //
3029  // extra diagonal
3030  //
3031  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
3032  {
3033  iloc = fe.patternFirst ( i );
3034  jloc = fe.patternSecond ( i );
3035  s = 0;
3036  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3037  {
3038  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3039  {
3040  coef_v[ icoor ] = 0.;
3041  }
3042 
3043  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3044  {
3045  for ( UInt iloc = 0; iloc < nbN2; iloc++ )
3046  {
3047  coef_v[ icoor ] += vec_loc.vec() [ iloc + icoor * nbN2 ] * fe2.phi ( iloc, ig );
3048  }
3049  }
3050 
3051  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3052  {
3053  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3054  {
3055  s += coef_v[ icoor ] * fe.phiDer ( iloc, icoor, ig ) * coef_v[ jcoor ] * fe.phiDer ( jloc, jcoor, ig )
3056  * fe.weightDet ( ig );
3057  }
3058  }
3059  }
3060  coef_s = coef * s;
3061  mat ( iloc, jloc ) += coef_s;
3062  mat ( jloc, iloc ) += coef_s;
3063  }
3064  // copy on the other components (if necessary, i.e. if nb>1)
3065  for ( int icomp = 1; icomp < nb; icomp++ )
3066  {
3067  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icomp, jblock + icomp );
3068  for ( i = 0; i < fe.nbDiag(); i++ )
3069  {
3070  iloc = fe.patternFirst ( i );
3071  mat_icomp ( iloc, iloc ) += mat ( iloc, iloc );
3072  }
3073  for ( i = fe.nbDiag(); i < fe.nbDiag() + fe.nbUpper(); i++ )
3074  {
3075  iloc = fe.patternFirst ( i );
3076  jloc = fe.patternSecond ( i );
3077  mat_icomp ( iloc, jloc ) += mat ( iloc, jloc );
3078  mat_icomp ( jloc, iloc ) += mat ( jloc, iloc );
3079  }
3080  }
3081 }
3082 
3083 
3084 //
3085 void grad ( const int icoor, Real coef, MatrixElemental& elmat,
3086  const CurrentFE& fe_u, const CurrentFE& fe_p,
3087  int iblock, int jblock )
3088 /*
3089  \int q_j \frac{\partial v_i}{\partial x_icoor}
3090 */
3091 {
3092  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
3093  UInt ig;
3094  UInt i, j;
3095  double s;
3096  for ( i = 0; i < fe_u.nbFEDof(); i++ )
3097  {
3098  for ( j = 0; j < fe_p.nbFEDof(); j++ )
3099  {
3100  s = 0;
3101  for ( ig = 0; ig < fe_u.nbQuadPt(); ig++ )
3102  // Be careful the minus is here and not in coef!!!!
3103 
3104  // wrong for different quadrules of fe_u and fe_p !! Martin P.
3105  // s -= fe_p.phi(j,ig)*fe_u.phiDer(i,icoor,ig)*fe_u.weightDet(ig);
3106 
3107  s -= fe_p.refFE().phi ( j, fe_u.quadRule().quadPointCoor ( ig, 0 ), fe_u.quadRule().quadPointCoor ( ig, 1 ),
3108  fe_u.quadRule().quadPointCoor ( ig, 2 ) ) * fe_u.phiDer ( i, icoor, ig ) * fe_u.weightDet ( ig );
3109  mat ( i, j ) += coef * s;
3110  }
3111  }
3112 }
3113 
3114 void div ( const int icoor, Real coef, MatrixElemental& elmat,
3115  const CurrentFE& fe_u, const CurrentFE& fe_p,
3116  int iblock, int jblock )
3117 /*
3118  \int q_i \frac{\partial v_j}{\partial x_icoor}
3119 */
3120 {
3121  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
3122  UInt ig;
3123  UInt i, j;
3124  double s;
3125  for (i = 0; i < fe_u.nbFEDof(); i++)
3126  {
3127  for (j = 0; j < fe_p.nbFEDof(); j++)
3128  {
3129  s = 0;
3130  for ( ig = 0; ig < fe_u.nbQuadPt(); ig++ )
3131  {
3132  s -= fe_u.phi (i, ig) * fe_p.phiDer (j, icoor, ig) * fe_u.weightDet ( ig );
3133  }
3134  mat ( i, j ) += coef * s;
3135  }
3136  }
3137 }
3138 
3139 void grad_div ( Real coef_grad, Real coef_div, MatrixElemental& elmat,
3140  const CurrentFE& fe_u, const CurrentFE& fe_p,
3141  int block_pres )
3142 /*
3143  \int q_j \frac{\partial v_i}{\partial x_icoor}
3144 */
3145 {
3146  double s;
3147  int iblock = block_pres - fe_u.nbLocalCoor();
3148  for ( UInt icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
3149  {
3150  MatrixElemental::matrix_view mat_grad = elmat.block ( iblock + icoor, block_pres );
3151  MatrixElemental::matrix_view mat_div = elmat.block ( block_pres , iblock + icoor );
3152  for ( UInt i = 0; i < fe_u.nbFEDof(); i++ )
3153  {
3154  for ( UInt j = 0; j < fe_p.nbFEDof(); j++ )
3155  {
3156  s = 0;
3157  for ( UInt ig = 0; ig < fe_u.nbQuadPt(); ig++ )
3158  {
3159  s -= fe_p.phi ( j, ig ) * fe_u.phiDer ( i, icoor, ig ) * fe_u.weightDet ( ig );
3160  }
3161  mat_grad ( i, j ) += coef_grad * s;
3162  mat_div ( j, i ) += coef_div * s;
3163  }
3164  }
3165  }
3166 }
3167 //
3168 void stab_stokes ( Real visc, Real coef_stab, MatrixElemental& elmat,
3169  const CurrentFE& fe, int block_pres )
3170 {
3171  MatrixElemental::matrix_view mat = elmat.block ( block_pres, block_pres );
3172  Real s, h = fe.diameter();
3173  Real fh2 = coef_stab * h * h / ( 2 * visc );
3174  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
3175  {
3176  for ( UInt j = 0; j < fe.nbFEDof(); j++ )
3177  {
3178  s = 0;
3179  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
3180  {
3181  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3182  {
3183  s += fe.phiDer ( i, icoor, ig ) * fe.phiDer ( j, icoor, ig ) * fe.weightDet ( ig );
3184  }
3185  }
3186  mat ( i, j ) -= fh2 * s;
3187  }
3188  }
3189 }
3190 
3191 /*
3192  * Fixed by Umberto Villa, Jan 2010
3193  */
3194 void advection ( Real coef, VectorElemental& vel,
3195  MatrixElemental& elmat, const CurrentFE& fe, int iblock, int jblock, int nb )
3196 {
3197  Matrix mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
3198  Real v_grad, s;
3199  Matrix v ( fe.nbQuadPt(), fe.nbLocalCoor() );
3200 
3201  //Evaluate the advective field at the quadrature nodes
3202  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3203  {
3204  VectorElemental::vector_view velicoor = vel.block ( icoor );
3205  for ( UInt iq = 0; iq < fe.nbQuadPt(); iq++ )
3206  {
3207  s = 0.;
3208  for ( UInt k = 0; k < fe.nbFEDof(); k++ )
3209  {
3210  s += velicoor ( k ) * fe.phi ( k, iq ); // velocity on the intgt point
3211  }
3212  v (iq, icoor) = s;
3213  }
3214  }
3215 
3216  //Assemble the local matrix
3217  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
3218  {
3219  for ( UInt j = 0; j < fe.nbFEDof(); j++ )
3220  {
3221  s = 0.;
3222  for ( UInt iq = 0; iq < fe.nbQuadPt(); iq++ )
3223  {
3224  v_grad = 0.;
3225  for ( int icoor = 0; icoor < ( int ) fe.nbLocalCoor(); icoor++ )
3226  {
3227  v_grad += v (iq, icoor) * fe.phiDer ( j, icoor, iq );
3228  }
3229 
3230  s += v_grad * fe.phi ( i, iq ) * fe.weightDet ( iq );
3231  }
3232  mat_tmp ( i, j ) = s * coef;
3233  }
3234  }
3235 
3236  // copy on the components
3237  for ( int icomp = 0; icomp < nb; icomp++ )
3238  {
3239  MatrixElemental::matrix_view mat_icomp = elmat.block ( iblock + icomp, jblock + icomp );
3240  for ( UInt i = 0; i < fe.nbDiag(); i++ )
3241  {
3242  for ( UInt j = 0; j < fe.nbDiag(); j++ )
3243  {
3244  mat_icomp ( i, j ) += mat_tmp ( i, j );
3245  }
3246  }
3247  }
3248 }
3249 
3250 void grad ( const int icoor, const VectorElemental& vec_loc, MatrixElemental& elmat,
3251  const CurrentFE& fe1, const CurrentFE& fe2,
3252  int iblock, int jblock )
3253 /*
3254  \int q_j \frac{\partial v_i}{\partial x_icoor}
3255 */
3256 {
3257  //
3258  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
3259 
3260  if ( iblock == jblock )
3261  {
3262  UInt iq;
3263  UInt i, j;
3264  double s, coef;
3265  UInt nbN1 = fe1.nbFEDof();
3266  for ( i = 0; i < nbN1; i++ )
3267  {
3268  for ( j = 0; j < fe2.nbFEDof(); j++ )
3269  {
3270  s = 0;
3271  for ( iq = 0; iq < fe1.nbQuadPt(); iq++ )
3272  {
3273  coef = 0;
3274  for ( UInt iloc = 0; iloc < nbN1; iloc++ )
3275  {
3276  coef += vec_loc.vec() [ iloc + icoor * nbN1 ] * fe1.phi ( iloc, iq );
3277  }
3278 
3279  s += coef * fe2.phi ( i, iq ) * fe1.phiDer ( j, icoor, iq ) * fe1.weightDet ( iq );
3280  } // Loop on quadrature nodes
3281 
3282  mat ( i, j ) += s;
3283  } //Loop on j
3284  } // Loop on i
3285  } // if
3286 
3287 }
3288 
3289 //
3290 // \! Gradient operator in the skew-symmetric form for NS Problems
3291 // A. Veneziani - December 2002
3292 // \!
3293 void grad_ss ( const int icoor, const VectorElemental& vec_loc, MatrixElemental& elmat,
3294  const CurrentFE& fe1, const CurrentFE& fe2,
3295  int iblock, int jblock )
3296 /*
3297  \int vloc(icoor) \frac{\partial v_i}{\partial x_icoor} v_j + 1/2*\frac{\partial v_icoor}{\partial x_icoor} v_i v_j
3298 */
3299 {
3300  //
3301  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
3302 
3303  if ( iblock == jblock )
3304  {
3305  UInt iq;
3306  UInt i, j;
3307  double s, coef, coef_div;
3308  UInt nbN1 = fe1.nbFEDof();
3309  for ( i = 0; i < nbN1; i++ )
3310  {
3311  for ( j = 0; j < fe2.nbFEDof(); j++ )
3312  {
3313  s = 0;
3314  for ( iq = 0; iq < fe1.nbQuadPt(); iq++ )
3315  {
3316  coef = 0;
3317  coef_div = 0;
3318 
3319  for ( UInt iloc = 0; iloc < nbN1; iloc++ )
3320  {
3321  coef += vec_loc.vec() [ iloc + icoor * nbN1 ] * fe1.phi ( iloc, iq );
3322  coef_div += vec_loc.vec() [ iloc + icoor * nbN1 ] * fe1.phiDer ( iloc, icoor, iq );
3323  }
3324 
3325  s += ( coef * fe1.phiDer ( j, icoor, iq ) + 0.5 * coef_div * fe1.phi ( j, iq ) ) * fe2.phi ( i, iq ) * fe1.weightDet ( iq );
3326  } // Loop on quadrature nodes
3327 
3328  mat ( i, j ) += s;
3329  } //Loop on j
3330  } // Loop on i
3331  } // if
3332 
3333 }
3334 
3335 // /!
3336 // Gradient operator where the convective term is based on a local vector
3337 // living on the basis given by fe3
3338 // It is useful for advection diffusion problems driven by a NS problem
3339 // !/
3340 void grad ( const int icoor,
3341  const VectorElemental& vec_loc,
3342  MatrixElemental& elmat,
3343  const CurrentFE& fe1,
3344  const CurrentFE& fe2,
3345  const CurrentFE& fe3,
3346  int iblock, int jblock )
3347 /*
3348  \int q_j \frac{\partial v_i}{\partial x_icoor}
3349 */
3350 {
3351  //
3352  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
3353 
3354  if ( iblock == jblock )
3355  {
3356  UInt iq;
3357  UInt i, j;
3358  double s, coef;
3359  UInt nbN1 = fe1.nbFEDof();
3360  UInt nbN3 = fe3.nbFEDof();
3361 
3362  for ( i = 0; i < nbN1; i++ )
3363  {
3364  for ( j = 0; j < fe2.nbFEDof(); j++ )
3365  {
3366  s = 0;
3367  for ( iq = 0; iq < fe1.nbQuadPt(); iq++ )
3368  {
3369  coef = 0;
3370 
3371  for ( UInt iloc = 0; iloc < nbN3; iloc++ )
3372  {
3373  coef += vec_loc.vec() [iloc + icoor * nbN3] * fe3.phi (iloc, iq);
3374  }
3375 
3376  s += coef * fe2.phi (i, iq) * fe1.phiDer (j, icoor, iq) * fe1.weightDet (iq);
3377  } // Loop on quadrature nodes
3378 
3379  mat (i, j) += s;
3380  } //Loop on j
3381  } // Loop on i
3382  } // if
3383 
3384 }
3385 
3386 
3387 // Convective term with the velocity given in the quadrature nodes
3388 void grad ( const int& icoor,
3389  const std::vector<Real>& localVector,
3390  MatrixElemental& elmat,
3391  const CurrentFE& currentFE1,
3392  const CurrentFE& currentFE2,
3393  const int& iblock,
3394  const int& jblock)
3395 {
3396  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
3397 
3398  // This term concerns only the diagonal blocks (same components)
3399  if ( iblock == jblock )
3400  {
3401 
3402  for ( UInt iNode1 (0); iNode1 < currentFE1.nbFEDof(); iNode1++ )
3403  {
3404  for ( UInt jNode2 (0); jNode2 < currentFE2.nbFEDof(); jNode2++ )
3405  {
3406  Real sum (0.0);
3407  for ( UInt iq (0); iq < currentFE1.nbQuadPt(); iq++ )
3408  {
3409  // Velocity in the quadrature node, component icoor
3410  double coef (localVector[icoor * currentFE1.nbQuadPt() + iq]);
3411 
3412  sum += coef
3413  * currentFE2.phi (iNode1, iq)
3414  * currentFE1.phiDer (jNode2, icoor, iq)
3415  * currentFE1.weightDet (iq);
3416  } // Loop on quadrature nodes
3417 
3418  mat (iNode1, jNode2) += sum;
3419  } //Loop on j
3420  } // Loop on i
3421  } // if
3422 
3423 }
3424 //
3425 //
3426 /*
3427 //////////////////////
3428 // NONLINEAR TERMS
3429 /////////////////////
3430 //
3431 //
3432 //--------------------
3433 // Jacobian: 2*\Sum V_k \Int \phi_k \phi_i \phi_j
3434 // and
3435 // Vector F(V) = \Sum V_k \Sum V_j \Int \phi_k \phi_i \phi_j
3436 //-------------------
3437 void quad(std::vector<Real> coef, MatrixElemental& elmat, VectorElemental& elvec,
3438 const CurrentFE& fe,int iblock=0,int jblock=0)
3439 {
3440 MatrixElemental::matrix_view mat = elmat.block(iblock,jblock);
3441 int i,ig,iq,siz;
3442 int iloc,jloc,qloc;
3443 Real s,coef_s;
3444 siz=coef.size();
3445 ASSERT(siz==fe.nbDiag,
3446 "Error in building Local Matrix of the quadratic term");
3447 //
3448 // diagonal
3449 //
3450 for(i=0;i<fe.nbDiag();i++){
3451 iloc = fe.patternFirst(i);
3452 s = 0;
3453 for (iq=0;i<siz;++i){
3454 qloc = fe.patternFirst(iq);
3455 for(ig=0;ig<fe.nbQuadPt();ig++)
3456 s += coef(iq)*fe.phi(qloc,ig)*fe.phi(iloc,ig)*fe.phi(iloc,ig)*
3457 fe.weightDet(ig);
3458 }
3459 mat(iloc,iloc) += 2*s;
3460 elvec(iloc) += coef(iloc)*s;
3461 }
3462 //
3463 // extra diagonal
3464 //
3465 for(i=fe.nbDiag();i<fe.nbDiag()+fe.nbUpper();i++){
3466 iloc = fe.patternFirst(i);
3467 jloc = fe.patternSecond(i);
3468 s = 0;
3469 for (iq=0;i<siz;++i){
3470 qloc = fe.patternFirst(iq);
3471 for(ig=0;ig<fe.nbQuadPt();ig++)
3472 s += coef(iq)*fe.phi(qloc,ig)*
3473 fe.phi(iloc,ig)*fe.phi(jloc,ig)*fe.weightDet(ig);
3474 }
3475 mat(iloc,jloc) += 2*s;
3476 mat(jloc,iloc) += 2*s;
3477 }
3478 }
3479 */
3480 //----------------------------------------------------------------------
3481 // Element vector operator
3482 //----------------------------------------------------------------------
3483 void source ( Real constant, VectorElemental& elvec, const CurrentFE& fe, int iblock )
3484 {
3485  UInt i, ig;
3486  VectorElemental::vector_view vec = elvec.block ( iblock );
3487  Real s;
3488  for ( i = 0; i < fe.nbFEDof(); i++ )
3489  {
3490  s = 0;
3491  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3492  {
3493  s += fe.phi ( i, ig ) * fe.weightDet ( ig );
3494  }
3495  vec ( i ) += constant * s;
3496  }
3497 }
3498 
3499 void source ( Real coef, VectorElemental& f, VectorElemental& elvec, const CurrentFE& fe,
3500  int fblock, int eblock )
3501 /*
3502  compute \int f \phi_i
3503  where f is given on the dof of this element
3504  fblock is the block of the values read in f
3505  eblock is the block where the result is writen in the elvec
3506 */
3507 {
3508  UInt i, ig;
3509  VectorElemental::vector_view vec = elvec.block ( eblock );
3510  VectorElemental::vector_view vecf = f.block ( fblock );
3511  Real f_ig;
3512 
3513  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3514  {
3515  f_ig = 0.;
3516  for ( i = 0; i < fe.nbFEDof(); i++ )
3517  {
3518  f_ig += vecf ( i ) * fe.phi ( i, ig );
3519  }
3520  for ( i = 0; i < fe.nbFEDof(); i++ )
3521  {
3522  vec ( i ) += coef * f_ig * fe.phi ( i, ig ) * fe.weightDet ( ig );
3523  }
3524  }
3525 }
3526 
3527 
3528 void source_mass (const std::vector<Real>& constant, VectorElemental& elvec, const CurrentFE& currentFe, const int& iblock)
3529 {
3530  VectorElemental::vector_view vec = elvec.block ( iblock );
3531  for (UInt iterNode (0); iterNode < currentFe.nbFEDof(); ++iterNode )
3532  {
3533  for ( UInt iterQuadNode (0); iterQuadNode < currentFe.nbQuadPt(); iterQuadNode++ )
3534  {
3535  vec (iterNode) += constant[iterQuadNode]
3536  * currentFe.phi ( iterNode, iterQuadNode )
3537  * currentFe.weightDet ( iterQuadNode );
3538  }
3539  }
3540 }
3541 
3542 void source_stiff (const std::vector<Real>& constant, VectorElemental& elvec, const CurrentFE& currentFe, const int& iblock)
3543 {
3544  VectorElemental::vector_view vec = elvec.block ( iblock );
3545  const UInt nbQuadPt (currentFe.nbQuadPt() );
3546  for (UInt iterNode (0); iterNode < currentFe.nbFEDof(); ++iterNode )
3547  {
3548  for ( UInt iterQuadNode (0); iterQuadNode < nbQuadPt; iterQuadNode++ )
3549  {
3550  for (UInt iterGradComp (0); iterGradComp < currentFe.nbLocalCoor(); ++iterGradComp)
3551  {
3552  vec (iterNode) += constant[iterQuadNode + iterGradComp * nbQuadPt]
3553  * currentFe.phiDer ( iterNode, iterGradComp, iterQuadNode )
3554  * currentFe.weightDet ( iterQuadNode );
3555  }
3556  }
3557  }
3558 }
3559 
3560 
3561 
3562 void source_divuq (Real alpha, VectorElemental& uLoc, VectorElemental& elvec, const CurrentFE& fe_u, const CurrentFE& fe_p, int iblock )
3563 {
3564  UInt i, j, ic, iq;
3565  VectorElemental::vector_view vec = elvec.block ( iblock );
3566  Real s;
3567 
3568  for (i = 0; i < fe_p.nbFEDof(); i++)
3569  {
3570  s = 0;
3571  for (iq = 0; iq < fe_p.nbQuadPt(); ++iq )
3572  for (j = 0; j < fe_u.nbFEDof(); ++j)
3573  for (ic = 0; ic < fe_u.nbLocalCoor(); ++ic)
3574  {
3575  s += uLoc[ic * fe_u.nbFEDof() + j] * fe_u.phiDer (j, ic, iq) * fe_p.phi (i, iq) * fe_p.weightDet ( iq );
3576  }
3577 
3578  vec ( i ) += s * alpha;
3579  }
3580 }
3581 
3582 
3583 void source_gradpv (Real alpha, VectorElemental& pLoc, VectorElemental& elvec, const CurrentFE& fe_p, const CurrentFE& fe_u, int iblock )
3584 {
3585  UInt i, j, iq;
3586  VectorElemental::vector_view vec = elvec.block ( iblock );
3587  Real s;
3588 
3589  for ( i = 0; i < fe_u.nbFEDof(); i++ )
3590  {
3591  s = 0;
3592  for (iq = 0; iq < fe_u.nbQuadPt(); ++iq )
3593  for (j = 0; j < fe_p.nbFEDof(); ++j)
3594  {
3595  s += pLoc[j] * fe_p.phiDer (j, iblock, iq) * fe_u.phi (i, iq) * fe_u.weightDet ( iq );
3596  }
3597  vec ( i ) += s * alpha;
3598  }
3599 }
3600 
3601 
3602 
3603 
3604 void source_fhn ( Real coef_f, Real coef_a, VectorElemental& u, VectorElemental& elvec, const CurrentFE& fe,
3605  int fblock, int eblock )
3606 /*
3607  compute \int coef_f u(1-u)(u-coef_a) \phi_i
3608  (right-hand side for the Fitzhugh-Nagumo equations)
3609  where u is given on the dof of this element
3610  fblock is the block of the values read in f
3611  eblock is the block where the result is writen in the elvec
3612 */
3613 {
3614  UInt i, ig;
3615  VectorElemental::vector_view vec = elvec.block ( eblock );
3616  VectorElemental::vector_view vecu = u.block ( fblock );
3617  Real f_ig;
3618 
3619  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3620  {
3621  f_ig = 0.;
3622  for ( i = 0; i < fe.nbFEDof(); i++ )
3623  {
3624  f_ig += vecu ( i ) * fe.phi ( i, ig );
3625  }
3626  for ( i = 0; i < fe.nbFEDof(); i++ )
3627  {
3628  vec ( i ) += coef_f * f_ig * ( 1 - f_ig ) * ( f_ig - coef_a ) * fe.phi ( i, ig ) * fe.weightDet ( ig );
3629  }
3630  }
3631 
3632 }
3633 
3634 //! \f$(beta\cdot\nabla u^k, v )\f$
3635 void source_advection ( const Real& coefficient, const VectorElemental& beta_loc, const VectorElemental& uk_loc,
3636  VectorElemental& elvec, const CurrentFE& fe )
3637 {
3638  boost::multi_array<Real, 2> guk (
3639  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()]);
3640  boost::multi_array<Real, 2> conv (
3641  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()]);
3642  std::vector<Real> beta (fe.nbLocalCoor() );
3643 
3644  Real s;
3645 
3646  UInt ig, icoor, jcoor, i;
3647 
3648  // loop on quadrature points
3649  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3650  {
3651 
3652  // Interpolating beta
3653  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3654  {
3655  // Evaluate beta in the Gauss Point
3656  s = 0.0;
3657  for ( i = 0; i < fe.nbFEDof(); i++ )
3658  {
3659  s += fe.phi ( i, ig ) * beta_loc.vec() [ i + icoor * fe.nbFEDof() ];
3660  }
3661  beta[ icoor ] = s;
3662  }
3663 
3664  // Interpolation of grad u
3665  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3666  {
3667  // loop on the derivative variable
3668  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3669  {
3670  // Evaluate the derivative in the gauss point
3671  s = 0.0;
3672  for ( i = 0; i < fe.nbFEDof(); i++ )
3673  {
3674  // grad u^k at a quadrature point
3675  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ];
3676  }
3677  guk[ icoor ][ jcoor ] = s;
3678  }
3679  }
3680 
3681  // beta*(\grad u^k) at each quadrature point
3682  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3683  {
3684  s = 0.0;
3685  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3686  {
3687  s += beta[ icoor ] * guk[ jcoor ][ icoor ];
3688  }
3689  conv[ ig ][ jcoor ] = s;
3690  }
3691  }
3692 
3693  //
3694  // Numerical integration
3695  //
3696 
3697  // loop on coordinates, i.e. loop on elementary vector blocks
3698  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3699  {
3700 
3701  VectorElemental::vector_view vec = elvec.block ( icoor );
3702 
3703  // loop on nodes, i.e. loop on components of this block
3704  for ( i = 0; i < fe.nbFEDof(); i++ )
3705  {
3706 
3707  // loop on quadrature points
3708  s = 0;
3709  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3710  {
3711  s += conv[ ig ][ icoor ] * fe.phi ( i, ig ) * fe.wDetJacobian ( ig );
3712  }
3713  vec ( i ) += coefficient * s;
3714  }
3715  }
3716 }
3717 
3718 // coef * ( - \grad w^k :[I\div d - (\grad d)^T] u^k + convect^T[I\div d - (\grad d)^T] (\grad u^k)^T , v ) for Newton FSI
3719 //
3720 // Remark: convect = u^n-w^k relative vel.
3721 //
3722 void source_mass1 ( Real coef, const VectorElemental& uk_loc, const VectorElemental& wk_loc, const VectorElemental& convect_loc,
3723  const VectorElemental& d_loc, VectorElemental& elvec, const CurrentFE& fe )
3724 {
3725 
3726  boost::multi_array<Real, 2> A (
3727  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()]);
3728  boost::multi_array<Real, 2> B (
3729  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()]);
3730  boost::multi_array<Real, 2> uk (
3731  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()]);
3732  boost::multi_array<Real, 3> guk (
3733  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()][fe.nbLocalCoor()]);
3734  std::vector<Real> dw (fe.nbLocalCoor() );
3735  std::vector<Real> aux (fe.nbQuadPt() );
3736  std::vector<Real> convect (fe.nbLocalCoor() );
3737  boost::multi_array<Real, 2> convect_A (
3738  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()]);
3739 
3740  Real s, sA, sB, sG;
3741 
3742 
3743  UInt icoor, jcoor, ig, i;
3744  // loop on quadrature points
3745  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3746  {
3747 
3748  // loop on space coordindates
3749  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3750  {
3751 
3752  // each compontent of uk at each quadrature points
3753  s = 0.0;
3754  for ( i = 0; i < fe.nbFEDof(); i++ )
3755  {
3756  s += fe.phi ( i, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ];
3757  }
3758  uk[ ig ][ icoor ] = s;//uk_x(pt_ig), uk_y(pt_ig), uk_z(pt_ig)
3759 
3760  // each compontent of convect at this quadrature point
3761  s = 0.0;
3762  for ( i = 0; i < fe.nbFEDof(); i++ )
3763  {
3764  s += fe.phi ( i, ig ) * convect_loc.vec() [ i + icoor * fe.nbFEDof() ];
3765  }
3766  convect[ icoor ] = s;
3767 
3768 
3769  // loop on space coordindates
3770  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3771  {
3772  sB = 0.0;
3773  sA = 0.0;
3774  sG = 0.0;
3775  for ( i = 0; i < fe.nbFEDof(); i++ )
3776  {
3777  sG += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at each quadrature point
3778  sB -= fe.phiDer ( i, jcoor, ig ) * wk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad (- w^k) at this quadrature point
3779  sA -= fe.phiDer ( i, icoor, ig ) * d_loc.vec() [ i + jcoor * fe.nbFEDof() ]; // - (\grad d) ^T at this quadrature point
3780  }
3781  guk[ ig ][ icoor ][ jcoor ] = sG; // \grad u^k at each quadrature point
3782  B[ icoor ][ jcoor ] = sB; // \grad (convect) at this quadrature point
3783  A[ icoor ][ jcoor ] = sA; // -(\grad d) ^T at this quadrature point
3784  }
3785  }
3786 
3787  s = 0.0;
3788  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3789  {
3790  s -= A[ jcoor ][ jcoor ]; // \div d at this quadrature point ( - trace( A ) )
3791  }
3792 
3793  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3794  {
3795  A[ jcoor ][ jcoor ] += s; // I\div d - (\grad d)^T at this quadrature point (A+I(-tr(A)))
3796  }
3797 
3798  s = 0;
3799  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3800  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3801  {
3802  s += B[ icoor ][ jcoor ] * A[ icoor ][ jcoor ]; // \grad (-w^k):[I\div d - (\grad d)^T] at each quadrature point
3803  }
3804  aux[ ig ] = s;
3805 
3806  s = 0;
3807  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3808  {
3809  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3810  {
3811  s += convect[ icoor ] * A[ icoor ][ jcoor ]; // convect^T [I\div d - (\grad d)^T]
3812  }
3813  convect_A[ ig ][ jcoor ] = s;
3814  }
3815  // std::cout<<"aux "<<aux[ig]<<std::endl;
3816 
3817  // for(icoor=0; icoor<fe.nbLocalCoor(); ++icoor)
3818  // for(jcoor=0; jcoor<fe.nbLocalCoor(); ++jcoor)
3819  // {
3820  // std::cout<<" Atrue ["<<icoor<<"]["<<jcoor<<"] = "<<A[icoor][jcoor]<<std::endl;
3821  // }
3822 
3823  }
3824 
3825  // At this point we have:
3826  // v \grad u^k at each quadrature point: guk
3827  // v convect^T [I\div d - (\grad d)^T] at each quadrature point: convect_A
3828  // v \grad (-w^k):[I\div d - (\grad d)^T]: aux
3829 
3830 
3831  //
3832  // Numerical integration
3833  //
3834 
3835  // loop on coordinates, i.e. loop on elementary vector blocks
3836  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3837  {
3838 
3839  // the block iccor of the elementary vector
3840  VectorElemental::vector_view vec = elvec.block ( icoor );
3841 
3842  // loop on nodes, i.e. loop on components of this block
3843  for ( i = 0; i < fe.nbFEDof(); i++ )
3844  {
3845 
3846  // loop on quadrature points
3847  s = 0;
3848  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3849  {
3850 
3851  // \grad ( - w^k ):[I\div d - (\grad d)^T] \phi_i
3852  s += aux[ ig ] * uk[ ig ][ icoor ] * fe.phi ( i, ig ) * fe.weightDet ( ig );
3853 
3854  // convect^T [I\div d - (\grad d)^T] (\grad u^k)^T \phi_i
3855  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3856  {
3857  s += convect_A[ ig ][ jcoor ] * guk[ ig ][ icoor ][ jcoor ] * fe.phi ( i, ig ) * fe.weightDet ( ig );
3858  }
3859  }
3860  vec ( i ) += coef * s;
3861  }
3862  }
3863 }
3864 
3865 
3866 
3867 //
3868 // coef * ( \grad u^k dw, v ) for Newton FSI
3869 //
3870 //
3871 void source_mass2 ( Real coef, const VectorElemental& uk_loc, const VectorElemental& dw_loc,
3872  VectorElemental& elvec, const CurrentFE& fe )
3873 {
3874 
3875  boost::multi_array<Real, 2> guk (
3876  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()]);
3877  std::vector<Real> dw (fe.nbLocalCoor() );
3878  boost::multi_array<Real, 2> aux (
3879  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()]);
3880 
3881  Real s;
3882 
3883  UInt ig, icoor, jcoor, i;
3884 
3885  // loop on quadrature points
3886  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3887  {
3888 
3889  // loop on space coordinates
3890  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3891  {
3892 
3893  // each compontent (icoor) of dw at this quadrature point
3894  s = 0.0;
3895  for ( i = 0; i < fe.nbFEDof(); i++ )
3896  {
3897  s += fe.phi ( i, ig ) * dw_loc.vec() [ i + icoor * fe.nbFEDof() ];
3898  }
3899  dw[ icoor ] = s;
3900 
3901  // loop on space coordinates
3902  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3903  {
3904  s = 0.0;
3905  for ( i = 0; i < fe.nbFEDof(); i++ )
3906  {
3907  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at a quadrature point
3908  }
3909  guk[ icoor ][ jcoor ] = s;
3910  }
3911  }
3912 
3913  // (\grad u^k)dw at each quadrature point
3914  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3915  {
3916  s = 0.0;
3917  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3918  {
3919  s += guk[ icoor ][ jcoor ] * dw[ jcoor ];
3920  }
3921  aux[ ig ][ icoor ] = s;
3922  }
3923  }
3924 
3925  //
3926  // Numerical integration
3927  //
3928 
3929  // loop on coordinates, i.e. loop on elementary vector blocks
3930  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3931  {
3932 
3933  VectorElemental::vector_view vec = elvec.block ( icoor );
3934 
3935  // loop on nodes, i.e. loop on components of this block
3936  for ( i = 0; i < fe.nbFEDof(); i++ )
3937  {
3938 
3939  // loop on quadrature points
3940  s = 0;
3941  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3942  {
3943  s += aux[ ig ][ icoor ] * fe.phi ( i, ig ) * fe.weightDet ( ig );
3944  }
3945  vec ( i ) += coef * s;
3946  }
3947  }
3948 }
3949 
3950 
3951 
3952 //
3953 // coef * ( \grad u^n :[2 I \div d - (\grad d)^T] u^k , v ) for Newton FSI
3954 //
3955 //
3956 void source_mass3 ( Real coef, const VectorElemental& un_loc, const VectorElemental& uk_loc, const VectorElemental& d_loc,
3957  VectorElemental& elvec, const CurrentFE& fe )
3958 {
3959 
3960  boost::multi_array<Real, 2> B (
3961  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()]);
3962  boost::multi_array<Real, 2> A (
3963  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()]);
3964  boost::multi_array<Real, 2> uk (
3965  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()]);
3966  std::vector<Real> aux (fe.nbQuadPt() );
3967 
3968  Real s, sA, sB;
3969 
3970 
3971  UInt icoor, jcoor, ig, i;
3972  // loop on quadrature points
3973  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
3974  {
3975 
3976  // loop on space coordindates
3977  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
3978  {
3979 
3980  // each compontent of uk at each quadrature points
3981  s = 0.0;
3982  for ( i = 0; i < fe.nbFEDof(); i++ )
3983  {
3984  s += fe.phi ( i, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ];
3985  }
3986  uk[ ig ][ icoor ] = s;
3987 
3988 
3989  // loop on space coordindates
3990  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
3991  {
3992  sB = 0.0;
3993  sA = 0.0;
3994  for ( i = 0; i < fe.nbFEDof(); i++ )
3995  {
3996  sB += fe.phiDer ( i, jcoor, ig ) * un_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^n at this quadrature point
3997  sA -= fe.phiDer ( i, icoor, ig ) * d_loc.vec() [ i + jcoor * fe.nbFEDof() ]; // - (\grad d) ^T at this quadrature point
3998  }
3999  B[ icoor ][ jcoor ] = sB; // \grad u^n at this quadrature point
4000  A[ icoor ][ jcoor ] = sA; // -(\grad d) ^T at this quadrature point
4001  }
4002  }
4003 
4004  s = 0.0;
4005  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4006  {
4007  s -= A[ jcoor ][ jcoor ]; // \div d at this quadrature point ( - trace( A ) )
4008  }
4009 
4010  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4011  {
4012  A[ jcoor ][ jcoor ] += 2 * s; // 2 * I\div d - (\grad d)^T at this quadrature point
4013  }
4014 
4015  s = 0;
4016  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4017  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4018  {
4019  s += B[ icoor ][ jcoor ] * A[ icoor ][ jcoor ]; // \grad u^n:[2 * I\div d - (\grad d)^T] at each quadrature point
4020  }
4021  aux[ ig ] = s;
4022  }
4023 
4024  // At this point we have:
4025  // v u^k at each quadrature point: uk
4026  // v \grad u^n:[ 2 * I\div d - (\grad d)^T]: aux
4027 
4028 
4029  //
4030  // Numerical integration
4031  //
4032 
4033  // loop on coordinates, i.e. loop on elementary vector blocks
4034  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4035  {
4036 
4037  // the block iccor of the elementary vector
4038  VectorElemental::vector_view vec = elvec.block ( icoor );
4039 
4040  // loop on nodes, i.e. loop on components of this block
4041  for ( i = 0; i < fe.nbFEDof(); i++ )
4042  {
4043  // loop on quadrature points
4044  s = 0;
4045  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
4046  // \grad u^n:[2 * I\div d - (\grad d)^T] u^k \phi_i
4047  {
4048  s += aux[ ig ] * uk[ ig ][ icoor ] * fe.phi ( i, ig ) * fe.weightDet ( ig );
4049  }
4050  vec ( i ) += coef * s;
4051  }
4052  }
4053 }
4054 
4055 
4056 
4057 
4058 
4059 
4060 //
4061 // coef * ( [-p^k I + 2*mu e(u^k)] [I\div d - (\grad d)^T] , \grad v ) for Newton FSI
4062 //
4063 void source_stress ( Real coef, Real mu, const VectorElemental& uk_loc, const VectorElemental& pk_loc,
4064  const VectorElemental& d_loc, VectorElemental& elvec, const CurrentFE& fe_u,
4065  const CurrentFE& fe_p )
4066 {
4067 
4068  boost::multi_array<Real, 3> B (
4069  boost::extents[fe_u.nbQuadPt()][fe_u.nbLocalCoor()][fe_u.nbLocalCoor()]);
4070  boost::multi_array<Real, 2> A (
4071  boost::extents[fe_u.nbLocalCoor()][fe_u.nbLocalCoor()]);
4072  boost::multi_array<Real, 2> guk (
4073  boost::extents[fe_u.nbLocalCoor()][fe_u.nbLocalCoor()]);
4074  boost::multi_array<Real, 2> sigma (
4075  boost::extents[fe_u.nbLocalCoor()][fe_u.nbLocalCoor()]);
4076 
4077  Real s, sA, sG, pk;
4078 
4079  UInt icoor, jcoor, kcoor, ig, i;
4080 
4081  // loop on quadrature points
4082  for ( ig = 0; ig < fe_u.nbQuadPt(); ig++ )
4083  {
4084 
4085  // loop on space coordinates
4086  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4087  {
4088 
4089  // loop on space coordindates
4090  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4091  {
4092  sA = 0.0;
4093  sG = 0.0;
4094  for ( i = 0; i < fe_u.nbFEDof(); i++ )
4095  {
4096  sG += fe_u.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe_u.nbFEDof() ]; // \grad u^k at this quadrature point
4097  sA -= fe_u.phiDer ( i, icoor, ig ) * d_loc.vec() [ i + jcoor * fe_u.nbFEDof() ]; // - (\grad d) ^T at this quadrature point
4098  }
4099  guk[ icoor ][ jcoor ] = sG;
4100  A[ icoor ][ jcoor ] = sA;
4101  }
4102  }
4103 
4104  s = 0.0;
4105  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4106  {
4107  s -= A[ jcoor ][ jcoor ]; // \div d at a quadrature point ( - trace( A ) )
4108  }
4109 
4110  // std::cout<<"div = "<< s <<std::endl;
4111 
4112  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4113  {
4114  A[ jcoor ][ jcoor ] += s; // I\div d - (\grad d)^T
4115  }
4116 
4117  pk = 0.0;
4118  for ( i = 0; i < fe_p.nbFEDof(); i++ )
4119  {
4120  pk += fe_p.phi ( i, ig ) * pk_loc.vec() [ i ]; // p^k at this quadrature point
4121  }
4122 
4123 
4124  // sigma = [-p^k I + 2*mu e(u^k)] a quadrature point
4125  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4126  {
4127  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4128  {
4129  sigma[ icoor ][ jcoor ] = mu * ( guk[ icoor ][ jcoor ] + guk[ jcoor ][ icoor ] );
4130  }
4131  sigma[ icoor ][ icoor ] -= pk;
4132  }
4133 
4134  // [-p^k I + 2*mu e(u^k)] [I\div d - (\grad d)^T] at each quadrature point
4135  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4136  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4137  {
4138  s = 0;
4139  for ( kcoor = 0; kcoor < fe_u.nbLocalCoor(); kcoor++ )
4140  {
4141  s += sigma[ icoor ][ kcoor ] * A[ kcoor ][ jcoor ];
4142  }
4143  B[ ig ][ icoor ][ jcoor ] = s;
4144  }
4145  }
4146 
4147  // for(icoor=0; icoor<fe_u.nbLocalCoor(); ++icoor)
4148  // for(jcoor=0; jcoor<fe_u.nbLocalCoor(); ++jcoor)
4149  // {
4150  // std::cout<<" Atrue ["<<icoor<<"]["<<jcoor<<"] = "<<A[icoor][jcoor]<<std::endl;
4151  // }
4152 
4153  // for(icoor=0; icoor<fe_u.nbLocalCoor(); ++icoor)
4154  // for(jcoor=0; jcoor<fe_u.nbLocalCoor(); ++jcoor)
4155  // {
4156  // double l=0.;
4157  // for(int e=0; e<fe_u.nbQuadPt(); ++e)
4158  // l+=B[e][icoor][jcoor];
4159  // std::cout<<" Btrue ["<<icoor<<"]["<<jcoor<<"] = "<<l<<std::endl;
4160  // }
4161 
4162  //
4163  // Numerical integration
4164  //
4165 
4166  // loop on coordinates, i.e. loop on elementary vector blocks
4167  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4168  {
4169 
4170  VectorElemental::vector_view vec = elvec.block ( icoor );
4171 
4172  // loop on nodes, i.e. loop on components of this block
4173  for ( i = 0; i < fe_u.nbFEDof(); i++ )
4174  {
4175 
4176  // loop on quadrature points
4177  s = 0;
4178  for ( ig = 0; ig < fe_u.nbQuadPt(); ig++ )
4179 
4180  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4181  {
4182  s += B[ ig ][ icoor ][ jcoor ] * fe_u.phiDer ( i, jcoor, ig ) * fe_u.weightDet ( ig );
4183  }
4184  vec ( i ) += coef * s;
4185  }
4186  }
4187 }
4188 
4189 
4190 
4191 //
4192 // + \mu ( \grad u^k \grad d + [\grad d]^T[\grad u^k]^T : \grad v )
4193 //
4194 void source_stress2 ( Real coef, const VectorElemental& uk_loc, const VectorElemental& d_loc, VectorElemental& elvec, const CurrentFE& fe_u )
4195 {
4196 
4197  boost::multi_array<Real, 3> A (
4198  boost::extents[fe_u.nbQuadPt()][fe_u.nbLocalCoor()][fe_u.nbLocalCoor()]);
4199  boost::multi_array<Real, 2> guk (
4200  boost::extents[fe_u.nbLocalCoor()][fe_u.nbLocalCoor()]);
4201  boost::multi_array<Real, 2> gd (
4202  boost::extents[fe_u.nbLocalCoor()][fe_u.nbLocalCoor()]);
4203 
4204  Real su, sd, s;
4205 
4206  UInt icoor, jcoor, kcoor, ig, i;
4207 
4208  // loop on quadrature points
4209  for ( ig = 0; ig < fe_u.nbQuadPt(); ig++ )
4210  {
4211 
4212  // loop on space coordinates
4213  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4214  {
4215 
4216  // loop on space coordindates
4217  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4218  {
4219  su = 0.0;
4220  sd = 0.0;
4221  for ( i = 0; i < fe_u.nbFEDof(); i++ )
4222  {
4223  su += fe_u.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe_u.nbFEDof() ]; // \grad u^k at this quadrature point
4224  sd += fe_u.phiDer ( i, jcoor, ig ) * d_loc.vec() [ i + icoor * fe_u.nbFEDof() ]; // \grad d at this quadrature point
4225  }
4226  guk[ icoor ][ jcoor ] = su;
4227  gd[ icoor ][ jcoor ] = sd;
4228  }
4229  }
4230 
4231 
4232  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4233  {
4234  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4235  {
4236  s = 0;
4237  for ( kcoor = 0; kcoor < fe_u.nbLocalCoor(); kcoor++ )
4238  {
4239  s += guk[ icoor ][ kcoor ] * gd[ kcoor ][ jcoor ] + gd[ kcoor ][ icoor ] * guk[ jcoor ][ kcoor ];
4240  }
4241  A[ ig ][ icoor ][ jcoor ] = s;
4242  }
4243  }
4244 
4245  // for(icoor=0; icoor<fe_u.nbLocalCoor(); ++icoor)
4246  // for(jcoor=0; jcoor<fe_u.nbLocalCoor(); ++jcoor)
4247  // {
4248  // std::cout<<" gdtrue ["<<icoor<<"]["<<jcoor<<"] = "<<gd[icoor][jcoor]<<std::endl;
4249  // }
4250  // for(icoor=0; icoor<fe_u.nbLocalCoor(); ++icoor)
4251  // for(jcoor=0; jcoor<fe_u.nbLocalCoor(); ++jcoor)
4252  // {
4253  // std::cout<<" Atrue ["<<icoor<<"]["<<jcoor<<"] = "<<A[ig][icoor][jcoor]<<std::endl;
4254  // }
4255  }
4256 
4257  //
4258  // Numerical integration
4259  //
4260  // loop on coordinates, i.e. loop on elementary vector blocks
4261  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4262  {
4263 
4264  VectorElemental::vector_view vec = elvec.block ( icoor );
4265 
4266  // loop on nodes, i.e. loop on components of this block
4267  for ( i = 0; i < fe_u.nbFEDof(); i++ )
4268  {
4269 
4270  // loop on quadrature points
4271  s = 0;
4272  for ( ig = 0; ig < fe_u.nbQuadPt(); ig++ )
4273 
4274  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4275  {
4276  s += fe_u.phiDer ( i, jcoor, ig ) * A[ ig ][ icoor ][ jcoor ] * fe_u.weightDet ( ig );
4277  }
4278  vec ( i ) += coef * s;
4279  }
4280  }
4281 }
4282 
4283 
4284 
4285 
4286 
4287 //
4288 // coef * ( (\grad u^k):[I\div d - (\grad d)^T] , q ) for Newton FSI
4289 //
4290 void source_press ( Real coef, const VectorElemental& uk_loc, const VectorElemental& d_loc, VectorElemental& elvec,
4291  const CurrentFE& fe_u, const CurrentFE& fe_p, int iblock )
4292 {
4293  boost::multi_array<Real, 2> A (
4294  boost::extents[fe_u.nbLocalCoor()][fe_u.nbLocalCoor()]);
4295  boost::multi_array<Real, 2> guk (
4296  boost::extents[fe_u.nbLocalCoor()][fe_u.nbLocalCoor()]);
4297  std::vector<Real> aux (fe_u.nbQuadPt() );
4298 
4299  VectorElemental::vector_view vec = elvec.block ( iblock );
4300 
4301  Real s, sA, sG;
4302  UInt icoor, jcoor, ig, i;
4303 
4304 
4305  // loop on quadrature points
4306  for ( ig = 0; ig < fe_u.nbQuadPt(); ig++ )
4307  {
4308 
4309  // loop on space coordinates
4310  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4311  {
4312 
4313  // loop on space coordinates
4314  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4315  {
4316  sA = 0.0;
4317  sG = 0.0;
4318  for ( i = 0; i < fe_u.nbFEDof(); i++ )
4319  {
4320  sG += fe_u.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe_u.nbFEDof() ]; // \grad u^k at a quadrature point
4321  sA -= fe_u.phiDer ( i, icoor, ig ) * d_loc.vec() [ i + jcoor * fe_u.nbFEDof() ]; // - (\grad d) ^T at a quadrature point
4322  }
4323  guk[ icoor ][ jcoor ] = sG;
4324  A[ icoor ][ jcoor ] = sA;
4325  }
4326  }
4327 
4328  s = 0.0;
4329  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4330  {
4331  s -= A[ jcoor ][ jcoor ]; // \div d at this quadrature point ( - trace( A ) )
4332  }
4333 
4334  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4335  {
4336  A[ jcoor ][ jcoor ] += s; // I\div d - (\grad d)^T at this quadrature point
4337  }
4338 
4339  s = 0;
4340  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4341  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4342  {
4343  s += guk[ icoor ][ jcoor ] * A[ icoor ][ jcoor ]; // \grad u^k : [I\div d - (\grad d)^T] at each quadrature point
4344  }
4345  aux[ ig ] = s;
4346  }
4347 
4348  //
4349  // Numerical integration
4350  //
4351 
4352  // Loop on nodes, i.e. loop on elementary vector components
4353  for ( i = 0; i < fe_p.nbFEDof(); i++ )
4354  {
4355 
4356  // loop on quadrature points
4357  s = 0;
4358  for ( ig = 0; ig < fe_u.nbQuadPt(); ig++ )
4359  {
4360  //std::cout << i << " " << ig << " " << fe_p.phi(i, ig) << " " << aux[ig] << std::endl;
4361  s += aux[ ig ] * fe_p.phi ( i, ig ) * fe_u.weightDet ( ig );
4362  }
4363  vec [ i ] += coef * s;
4364  }
4365 }
4366 
4367 
4368 //
4369 // coef * ( [I\div d - (\grad d)^T - \grad d] \grap p, \grad q ) for Newton FSI
4370 //
4371 void source_press2 ( Real coef, const VectorElemental& p_loc, const VectorElemental& d_loc, VectorElemental& elvec,
4372  const CurrentFE& fe, int iblock )
4373 {
4374  boost::multi_array<Real, 2> A (
4375  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()]);
4376  boost::multi_array<Real, 2> B (
4377  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()]);
4378  boost::multi_array<Real, 2> aux (
4379  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()]);
4380  std::vector<Real> gpk (fe.nbLocalCoor() );
4381 
4382  VectorElemental::vector_view vec = elvec.block ( iblock );
4383 
4384  Real s, sA, sG;
4385  UInt icoor, jcoor, ig, i;
4386 
4387 
4388  // loop on quadrature points
4389  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
4390  {
4391 
4392 
4393  // loop on space coordinates
4394  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4395  {
4396 
4397  sG = 0.0;
4398  for ( i = 0; i < fe.nbFEDof(); i++ )
4399  {
4400  sG += fe.phiDer ( i, icoor, ig ) * p_loc.vec() [ i ]; // \grad p^k at a quadrature point
4401  }
4402  gpk[ icoor ] = sG;
4403 
4404  // loop on space coordinates
4405  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4406  {
4407  sA = 0.0;
4408  for ( i = 0; i < fe.nbFEDof(); i++ )
4409  {
4410  sA -= fe.phiDer ( i, icoor, ig ) * d_loc.vec() [ i + jcoor * fe.nbFEDof() ]; // - (\grad d) ^T at a quadrature point
4411  }
4412  A[ icoor ][ jcoor ] = sA;
4413  B[ jcoor ][ icoor ] = sA;
4414  }
4415  }
4416 
4417  s = 0.0;
4418  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4419  {
4420  s -= A[ jcoor ][ jcoor ]; // \div d at this quadrature point ( - trace( A ) )
4421  }
4422 
4423  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4424  {
4425  A[ jcoor ][ jcoor ] += s; // I\div d - (\grad d)^T at this quadrature point
4426  }
4427 
4428  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4429  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4430  {
4431  A[ icoor ][ jcoor ] += B[ icoor ][ jcoor ]; // I\div d - (\grad d)^T - \grad d at this quadrature point
4432  }
4433 
4434  s = 0;
4435  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4436  {
4437  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4438  {
4439  s += A[ icoor ][ jcoor ] * gpk[jcoor]; // [I\div d - (\grad d)^T -\grad d] \grad p^k at each quadrature point
4440  }
4441  aux[ ig ][icoor] = s;
4442  }
4443  }
4444 
4445  //
4446  // Numerical integration
4447  //
4448 
4449  // Loop on nodes, i.e. loop on elementary vector components
4450  for ( i = 0; i < fe.nbFEDof(); i++ )
4451  {
4452 
4453  // loop on quadrature points
4454  s = 0;
4455  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
4456  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4457  {
4458  s += aux[ ig ][icoor] * fe.phiDer ( i, icoor, ig ) * fe.weightDet ( ig );
4459  }
4460  vec [ i ] += coef * s;
4461  }
4462 }
4463 
4464 
4465 
4466 
4467 
4468 //
4469 //Cholesky decomposition
4470 void choldc ( KNM<Real>& a, KN<Real>& p )
4471 {
4472  int i, j, k;
4473  Real sum;
4474 
4475  int n = a.N();
4476  for ( i = 0; i < n; i++ )
4477  {
4478  for ( j = i; j < n; j++ )
4479  {
4480  for ( sum = a ( i, j ), k = i - 1; k >= 0; k-- )
4481  {
4482  sum -= a ( i, k ) * a ( j, k );
4483  }
4484  if ( i == j )
4485  {
4486  p ( i ) = std::sqrt ( sum );
4487  }
4488  else
4489  {
4490  a ( j, i ) = sum / p ( i );
4491  }
4492  }
4493  }
4494 }
4495 //
4496 //Cholesky solution
4497 void cholsl ( KNM<Real>& a, KN<Real>& p, KN<Real>& b, KN<Real>& x )
4498 {
4499  int i, k;
4500  Real sum;
4501 
4502  int n = a.N();
4503  for ( i = 0; i < n; i++ )
4504  {
4505  for ( sum = b ( i ), k = i - 1; k >= 0; k-- )
4506  {
4507  sum -= a ( i, k ) * x ( k );
4508  }
4509  x ( i ) = sum / p ( i );
4510  }
4511  for ( i = n - 1; i >= 0; i-- )
4512  {
4513  for ( sum = x ( i ), k = i + 1; k < n; k++ )
4514  {
4515  sum -= a ( k, i ) * x ( k );
4516  }
4517  x ( i ) = sum / p ( i );
4518  }
4519 }
4520 
4521 //
4522 // coef * ( (\grad u^k):[I\div d - (\grad d)^T] , q ) for Newton FSI
4523 //
4524 void source_press ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat,
4525  const CurrentFE& fe_u, const CurrentFE& fe_p, ID mmDim , int iblock )
4526 {
4527  boost::multi_array<Real, 4> A (
4528  boost::extents[fe_u.nbLocalCoor()][fe_u.nbLocalCoor()][fe_u.nbFEDof()][fe_u.nbLocalCoor()]);
4529  boost::multi_array<Real, 2> guk (
4530  boost::extents[fe_u.nbLocalCoor()][fe_u.nbLocalCoor()]);
4531  boost::multi_array<Real, 3> aux (
4532  boost::extents[fe_u.nbQuadPt()][fe_u.nbFEDof()][fe_u.nbLocalCoor()]);
4533 
4534  Real l, sA, sG;
4535  UInt icoor, jcoor, ig;
4536 
4537 
4538  // loop on quadrature points
4539  for ( ig = 0; ig < fe_u.nbQuadPt(); ig++ )
4540  {
4541 
4542 
4543  ////INIT
4544  for (UInt p = 0; p < fe_u.nbLocalCoor(); ++p)
4545  {
4546  for (UInt q = 0; q < fe_u.nbLocalCoor(); ++q)
4547  {
4548  for (UInt d = 0; d < fe_u.nbFEDof(); ++d)
4549  {
4550  for (UInt e = 0; e < fe_u.nbLocalCoor(); ++e)
4551  {
4552  A[p][q][d][e] = 0.;
4553  }
4554  }
4555  //guk[p][q]=0.;
4556  }
4557  //for(int h=0; h<fe_u.nbFEDof(); ++h)
4558  //aux[ig][h][p]=0.;
4559  }
4560  ////END INIT
4561 
4562  // loop on space coordinates
4563  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4564  {
4565 
4566  //for ( i = 0;i < fe_u.nbFEDof();i++ )
4567  //for ( short k = 0;k < fe_u.nbLocalCoor();k++ )
4568  //sA[i][k] = 0.0;
4569  // loop on space coordinates
4570  for ( UInt kcoor = 0; kcoor < fe_u.nbLocalCoor(); kcoor++ )
4571  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4572  {
4573  sG = 0.0;
4574  for ( UInt i = 0; i < fe_u.nbFEDof(); i++ )
4575  {
4576  sG += fe_u.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe_u.nbFEDof() ]; // \grad u^k at a quadrature point
4577  sA = -fe_u.phiDer ( i, icoor, ig ) /**fe_u.phi( i, ig )*/ ;/** d_loc.vec() [ i + jcoor * fe_u.nbFEDof() ];*/ // - (\grad d) ^T at a quadrature point
4578  A[ icoor ][ jcoor ][i][jcoor] = sA; // -\delta_{jcoor kcoor} \partial_{icoor}
4579  }
4580  guk[ icoor ][ jcoor ] = sG;
4581  }
4582  }
4583 
4584 
4585 
4586 
4587 
4588  std::vector<Real> z (fe_u.nbLocalCoor() );
4589  for ( UInt i = 0; i < fe_u.nbFEDof(); i++ )
4590  {
4591  // for ( int kcoor = 0;kcoor < fe_u.nbLocalCoor();kcoor++ )
4592  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4593  {
4594  //if(icoor==jcoor)
4595  z[ jcoor ] = A[ jcoor ][ jcoor ][i][jcoor]; // -\delta_{jcoor kcoor} \partial_{icoor} + \delta_{jcoor icoor}\partial_{kcoor}
4596  }
4597 
4598  for ( jcoor = 0; jcoor < fe_p.nbLocalCoor(); jcoor++ )
4599  {
4600  for (UInt kcoor = 0; kcoor < fe_p.nbLocalCoor(); kcoor++ )
4601  {
4602  //if(icoor==jcoor)
4603  A[ jcoor ][ jcoor ][i][kcoor] -= z[ kcoor ]; // -\delta_{jcoor kcoor} \partial_{icoor} + \delta_{jcoor icoor}\partial_{kcoor}
4604  }
4605  }
4606 
4607 
4608 
4609  // for ( jcoor = 0;jcoor < fe_u.nbLocalCoor();jcoor++ )
4610  // {
4611  // //if(kcoor==jcoor)
4612  // for ( icoor = 0;icoor < fe_u.nbLocalCoor();icoor++ )
4613  // {
4614  // // if(icoor==jcoor)
4615  // A[ icoor ][ jcoor ][i][kcoor] -= A[ kcoor ][ icoor ][i][jcoor]; // \delta_{jcoor icoor}\partial_{kcoor}
4616  // //else
4617  // //f[icoor][jcoor]=0.;
4618  // }
4619  // }
4620  // }
4621  // for ( UInt kcoor = 0;kcoor < fe_u.nbLocalCoor();kcoor++ )
4622  // {
4623  // for ( jcoor = 0;jcoor < fe_u.nbLocalCoor();jcoor++ )
4624  // {
4625  // for ( icoor = 0;icoor < fe_u.nbLocalCoor();icoor++ )
4626  // {
4627  // //for ( jcoor = 0;jcoor < fe_u.nbLocalCoor();jcoor++ )
4628  // if(icoor==jcoor)
4629  // A[ icoor ][ jcoor ][i][kcoor] += f[icoor][jcoor]; // -\delta_{jcoor kcoor} \partial_{icoor} + \delta_{jcoor icoor}\partial_{kcoor}
4630  // // for ( i = 0;i < fe_u.nbFEDof();i++ )
4631  // // for(jcoor=0;jcoor<fe_u.nbLocalCoor();++jcoor)
4632  // // s[i][jcoor] = 0.0;
4633  // //for ( jcoor = 0;jcoor < fe_u.nbLocalCoor();jcoor++ )
4634  // //{
4635  // }
4636  // }
4637  // }
4638  // for ( UInt kcoor = 0;kcoor < fe_u.nbLocalCoor();kcoor++ )
4639  // {
4640  for ( UInt kcoor = 0; kcoor < fe_u.nbLocalCoor(); kcoor++ )
4641  {
4642  l = 0;
4643  //for ( kcoor = 0;kcoor < fe_u.nbLocalCoor();kcoor++ )
4644  for ( icoor = 0; icoor < fe_u.nbLocalCoor(); icoor++ )
4645  for ( jcoor = 0; jcoor < fe_u.nbLocalCoor(); jcoor++ )
4646  {
4647  l += guk[ icoor ][ jcoor ] * A[ icoor ][ jcoor ][i][kcoor]; // \grad u^k : [I\div d - (\grad d)^T] at each quadrature point
4648  //}
4649  }
4650  aux[ ig ][i][kcoor] = l;
4651  //}
4652  }
4653  //std::cout<<"aux[ ig ][i][jcoor] "<< aux[ ig ][i][jcoor] <<std::endl;
4654  }
4655  }
4656  //
4657  // Numerical integration
4658  //
4659  for ( UInt kcoor = 0; kcoor < fe_u.nbLocalCoor(); kcoor++ )
4660  {
4661  // for ( short l = 0;i < fe_u.nbFEDof();i++ )
4662  //for ( jcoor = 0;jcoor < fe_u.nbLocalCoor();jcoor++ )
4663  double l = 0.;
4664 
4665  MatrixElemental::matrix_view mat = elmat.block ( iblock, kcoor );
4666  for ( UInt j = 0; j < mmDim; j++ )
4667  for ( UInt i = 0; i < fe_p.nbFEDof(); i++ )
4668  {
4669  mat (i, j) = 0.;
4670  }
4671 
4672  // Loop on nodes, i.e. loop on elementary vector components
4673  for ( UInt j = 0; j < mmDim; j++ )
4674  for (UInt i = 0; i < fe_p.nbFEDof(); i++ )
4675  {
4676  l = 0.;
4677  // loop on quadrature points
4678  for ( ig = 0; ig < fe_u.nbQuadPt(); ig++ )
4679  {
4680  // std::cout << ig << " " << fe_p.phi(i, ig) << std::endl;
4681  l += aux[ ig ][j][kcoor] * fe_p.phi ( i, ig ) * fe_u.weightDet ( ig );
4682  }
4683  mat ( i , j) += coef * l;
4684  //std::cout<<"mat [ i , j] "<<mat [ i , j]<<std::endl;
4685  //std::cout<<"l "<<l<<std::endl;
4686  }
4687  }
4688 }
4689 
4690 
4691 
4692 //shape_terms_vel:
4693 //source_mass1
4694 // coef * ( - \grad w^k :[I\div d - (\grad d)^T] u^k + convect^T[I\div d - (\grad d)^T] (\grad u^k)^T , v ) for Newton FSI
4695 //
4696 // Remark: convect = u^n-w^k relative vel.
4697 //
4698 //source_stress
4699 // coef * ( [-p^k I + 2*mu e(u^k)] [I\div d - (\grad d)^T] , \grad v ) for Newton FSI
4700 //
4701 //source_stress2
4702 // \mu ( \grad u^k \grad d + [\grad d]^T[\grad u^k]^T : \grad v )
4703 //
4704 //source_mass3
4705 // 0.5 * ( \grad u^n :[2 I \div d - (\grad d)^T] u^k , v ) for Newton FSI
4706 //
4707 //optionally (if full implicit):
4708 //source_mass2
4709 // -rho * ( \grad u^k dw, v )
4710 //
4711 //convective term
4712 // rho * ( \grad u^k du, v )
4713 //
4714 void shape_terms (
4715  //const VectorElemental& d_loc,
4716  Real rho,
4717  Real mu,
4718  const VectorElemental& un_loc,
4719  const VectorElemental& uk_loc,
4720  const VectorElemental& wk_loc,
4721  const VectorElemental& convect_loc,
4722  const VectorElemental& pk_loc,
4723  MatrixElemental& elmat,
4724  const CurrentFE& fe,
4725  const CurrentFE& fe_p,
4726  ID /*mmDim*/,
4727  MatrixElemental& /*elmatP*/,
4728  int /*iblock*/,
4729  bool wImplicit,
4730  Real alpha,
4731  std::shared_ptr<MatrixElemental> elmat_convect
4732 )
4733 {
4734  // I div d - (grad d)^T
4735  boost::multi_array<Real, 4> eta (
4736  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbFEDof()][fe.nbLocalCoor()]);
4737  // I div d - (grad d)^T
4738  boost::multi_array<Real, 4> etaMass3 (
4739  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbFEDof()][fe.nbLocalCoor()]);
4740  // u^k
4741  boost::multi_array<Real, 2> uk (
4742  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()]);
4743  // grad u^k
4744  boost::multi_array<Real, 3> guk (
4745  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()][fe.nbLocalCoor()]);
4746  // div u^k
4747  std::vector<Real> duk (fe.nbQuadPt() );
4748  // grad u^k
4749  boost::multi_array<Real, 3> gun (
4750  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()][fe.nbLocalCoor()]);
4751  // convect
4752  std::vector<Real> convect (fe.nbLocalCoor() );
4753  // convect^T [ I div d - (grad d)^T ]
4754  boost::multi_array<Real, 4> convect_eta (
4755  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()][fe.nbFEDof()][fe.nbLocalCoor()]);
4756  // - p^k I + 2 mu e(u^k)
4757  boost::multi_array<Real, 2> sigma (
4758  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()]);
4759  // ( - p^k I + 2 mu e(u^k) )( I div d - (grad d)^T )
4760  boost::multi_array<Real, 5> B (
4761  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbFEDof()][fe.nbLocalCoor()]);
4762  // gd
4763  boost::multi_array<Real, 4> gd (
4764  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbFEDof()][fe.nbLocalCoor()]);
4765  // grad u^k grad d + (grad d)^T (grad u^k)^T
4766  boost::multi_array<Real, 5> A2 (
4767  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbFEDof()][fe.nbLocalCoor()]);
4768 
4769  boost::multi_array<Real, 4> aux (
4770  boost::extents[fe.nbQuadPt()][fe.nbLocalCoor()][fe.nbFEDof()][fe.nbLocalCoor()]);
4771  boost::multi_array<Real, 2> BMass (
4772  boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()]);
4773  boost::multi_array<Real, 3> auxMass (
4774  boost::extents[fe.nbQuadPt()][fe.nbFEDof()][fe.nbLocalCoor()]);
4775  boost::multi_array<Real, 3> auxMass3 (
4776  boost::extents[fe.nbQuadPt()][fe.nbFEDof()][fe.nbLocalCoor()]);
4777 
4778  Real s, sA, sB, sGk, sGn, pk;
4779 
4780  UInt icoor, jcoor, ig, kcoor, i, j;
4781 
4782  // loop on quadrature points
4783 
4784  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
4785  {
4786  ////INIT
4787  for (UInt p = 0; p < fe.nbLocalCoor(); ++p)
4788  {
4789  uk[ig][p] = 0.;
4790  convect[p] = 0.;
4791  for (UInt q = 0; q < fe.nbLocalCoor(); ++q)
4792  {
4793  guk[ig][p][q] = 0.;
4794  gun[ig][p][q] = 0.;
4795  sigma[p][q] = 0.;
4796  BMass[p][q] = 0.;
4797  for (UInt d = 0; d < fe.nbFEDof(); ++d)
4798  {
4799  auxMass[ ig ][d][q] = 0.;
4800  auxMass3[ ig ][d][q] = 0.;
4801  aux[ig][p][d][q] = 0.;
4802  convect_eta[ig][p][d][q] = 0.;
4803  for (UInt e = 0; e < fe.nbLocalCoor(); ++e)
4804  {
4805  gd[p][q][d][e] = 0.;
4806  eta[p][q][d][e] = 0.;
4807  etaMass3[p][q][d][e] = 0.;
4808  A2[ig][p][q][d][e] = 0.;
4809  B[ig][p][q][d][e] = 0.;
4810  }
4811  }
4812  }
4813  }
4814  ////////END INIT
4815 
4816  // loop on space coordindates
4817  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4818  {
4819 
4820  // each compontent of uk at each quadrature points
4821  s = 0.0;
4822  for ( i = 0; i < fe.nbFEDof(); i++ )
4823  {
4824  s += fe.phi ( i, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ];
4825  }
4826  uk[ ig ][ icoor ] = s;//uk_x(pt_ig), uk_y(pt_ig), uk_z(pt_ig)
4827 
4828  // each compontent of convect at this quadrature point
4829  s = 0.0;
4830  for ( i = 0; i < fe.nbFEDof(); i++ )
4831  {
4832  s += fe.phi ( i, ig ) * convect_loc.vec() [ i + icoor * fe.nbFEDof() ];
4833  }
4834  convect[ icoor ] = s;
4835 
4836 
4837  // loop on space coordindates
4838  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4839  {
4840  sB = 0.0;
4841  sGk = 0.0;
4842  sGn = 0.0;
4843  for ( i = 0; i < fe.nbFEDof(); i++ )
4844  {
4845  gd[ icoor ][ jcoor ][i][jcoor] = fe.phiDer ( i, jcoor, ig ) /** d_loc.vec() [ i + jcoor * fe.nbFEDof() ]*/;
4846 
4847  sGk += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at each quadrature point
4848  sGn += fe.phiDer ( i, jcoor, ig ) * un_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad u^k at each quadrature point
4849  sB -= fe.phiDer ( i, jcoor, ig ) * wk_loc.vec() [ i + icoor * fe.nbFEDof() ]; // \grad (- w^k) at this quadrature point
4850  sA = -fe.phiDer ( i, icoor, ig ) /** d_loc.vec() [ i + jcoor * fe.nbFEDof() ]*/; // - (\grad d) ^T at this quadrature point
4851  eta[ icoor ][ jcoor ][i][ jcoor ] = sA; // -(\grad d) ^T at this quadrature point
4852  etaMass3[ icoor ][ jcoor ][i][ jcoor ] = sA; // -(\grad d) ^T at this quadrature point
4853  }
4854  guk[ ig ][ icoor ][ jcoor ] = sGk; // \grad u^k at each quadrature point
4855  gun[ ig ][ icoor ][ jcoor ] = sGn; // \grad u^n at each quadrature point
4856  BMass[ icoor ][ jcoor ] = sB;
4857  }
4858  }
4859 
4860  //!a part of source_stress
4861  pk = 0.0;
4862  for ( i = 0; i < fe_p.nbFEDof(); i++ )
4863  {
4864  pk += fe_p.phi ( i, ig ) * pk_loc.vec() [ i ]; // p^k at this quadrature point
4865  }
4866 
4867  // sigma = [-p^k I + 2*mu e(u^k)] a quadrature point
4868  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4869  {
4870  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4871  {
4872  sigma[ icoor ][ jcoor ] = mu * ( guk[ig][ icoor ][ jcoor ] + guk[ig][ jcoor ][ icoor ] );
4873  }
4874  sigma[ icoor ][ icoor ] -= pk;
4875  }
4876 
4877  //!building the tensor \f$\eta = [I\nabla\cdot d - (\nabla d)^T]\f$
4878  for ( i = 0; i < fe.nbFEDof(); i++ )
4879  {
4880 
4881  std::vector<Real> z (fe.nbLocalCoor() );
4882  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4883  {
4884  //if(icoor==jcoor)
4885  z[ jcoor ] = eta[ jcoor ][ jcoor ][i][jcoor]; // -\delta_{jcoor kcoor} \partial_{icoor} + \delta_{jcoor icoor}\partial_{kcoor}
4886  }
4887  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4888  {
4889  for ( kcoor = 0; kcoor < fe.nbLocalCoor(); kcoor++ )
4890  {
4891  eta[ jcoor ][ jcoor ][i][kcoor] -= z[kcoor]; // -\delta_{jcoor kcoor} \partial_{icoor} + \delta_{jcoor icoor}\partial_{kcoor}
4892  }
4893  }
4894 
4895  //!source_mass1
4896 
4897  for ( kcoor = 0; kcoor < fe.nbLocalCoor(); kcoor++ )
4898  {
4899  s = 0;
4900  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4901  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4902  {
4903  s += BMass[ icoor ][ jcoor ] * eta[ icoor ][ jcoor ][i][kcoor]; // \grad (-w^k):[I\div d - (\grad d)^T] at each quadrature point
4904  }
4905  auxMass[ ig ][i][kcoor] = s;
4906 
4907  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4908  {
4909  s = 0.;
4910  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4911  {
4912  s += convect[ icoor ] * eta[ icoor ][ jcoor ][i][kcoor]; // convect^T [I\div d - (\grad d)^T]
4913  }
4914  convect_eta[ ig ][ jcoor ][i][kcoor] = s;
4915  }
4916  }
4917  //! source_stress
4918 
4919  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4920  {
4921  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4922  {
4923  for ( kcoor = 0; kcoor < fe.nbLocalCoor(); kcoor++ )
4924  {
4925  s = 0.;
4926  for (UInt zcoor = 0; zcoor < fe.nbLocalCoor(); zcoor++ )
4927  {
4928  s += sigma[ icoor ][ zcoor ] * eta[ zcoor ][ jcoor ][i][kcoor];
4929  }
4930  B[ ig ][ icoor ][ jcoor ][i][kcoor] = s;
4931  }
4932  }
4933  }
4934 
4935 
4936  //! source_stress2
4937  for ( kcoor = 0; kcoor < fe.nbLocalCoor(); kcoor++ )
4938  {
4939  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4940  {
4941  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4942  {
4943  s = 0.;
4944  for ( UInt zcoor = 0; zcoor < fe.nbLocalCoor(); zcoor++ )
4945  // \grad u^k \grad d + [\grad d]^T[\grad u^k]^T at each quadrature point
4946  {
4947  s += guk[ig][ icoor ][ zcoor ] * gd[ zcoor ][ jcoor ][i][kcoor] + gd[ zcoor ][ icoor ][i][kcoor] * guk[ig][ jcoor ][ zcoor ];
4948  }
4949  A2[ ig ][ icoor ][ jcoor ][i][kcoor] = s;
4950  }
4951  }
4952  }
4953 
4954  if (wImplicit) //source_mass2 and convective term derivative
4955  {
4956  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4957  {
4958  //s = 0.0;
4959  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4960  {
4961  aux[ ig ][ icoor ][i][jcoor] = guk[ig][ icoor ][ jcoor ] * fe.phi ( i, ig ) /**alpha*/;
4962  }
4963 
4964  }
4965  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4966  {
4967  duk[ig] += guk[ig][ jcoor ][ jcoor ]; //div uk
4968  }
4969  }
4970 
4971 
4972  //source_mass3
4973  // coef * ( \grad u^n :[2 I \div d - (\grad d)^T] u^k , v ) for Newton FSI
4974  //
4975  //
4976 
4977  //!building the tensor \f$\eta = [I\nabla\cdot d - (\nabla d)^T]\f$
4978  // for ( int kcoor = 0;kcoor < fe.nbLocalCoor();kcoor++ )
4979 
4980  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4981  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4982  {
4983  for (kcoor = 0; kcoor < fe.nbLocalCoor(); kcoor++ )
4984  {
4985  //if(icoor==jcoor)
4986  etaMass3[ icoor ][ jcoor ][i][kcoor] -= 2 * z[kcoor]; //! -2*\delta_{jcoor, kcoor} \partial_{icoor} + \delta_{jcoor, icoor}\partial_{kcoor}
4987  }
4988  }
4989 
4990  for (kcoor = 0; kcoor < fe.nbLocalCoor(); ++kcoor)
4991  {
4992  s = 0;
4993  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
4994  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
4995  {
4996  s += gun[ig][ icoor ][ jcoor ] * etaMass3/**/[ icoor ][ jcoor ][i][kcoor]; // \grad u^n:[/*2*/ * I\div d - (\grad d)^T] at each quadrature point
4997  }
4998  auxMass3[ ig ][i][kcoor] = s;
4999  //if fullImplicit un=uk
5000  }
5001 
5002  ///////////////////////////////////////////////////////////////////////
5003  }
5004 
5005  }
5006 
5007 
5008 
5009  //
5010  // Numerical integration
5011  //
5012  Real g = 0.;
5013  // loop on coordinates, i.e. loop on elementary vector blocks
5014  for ( icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
5015  {
5016  for ( kcoor = 0; kcoor < fe.nbLocalCoor(); kcoor++ )
5017  {
5018 
5019  // the block iccor of the elementary vector
5020  MatrixElemental::matrix_view mat = elmat.block ( icoor, kcoor );
5021  std::shared_ptr<MatrixElemental::matrix_view> mat_convect;
5022 
5023  if (elmat_convect.get() )
5024  {
5025  mat_convect.reset (new MatrixElemental::matrix_view (elmat_convect->block ( icoor, kcoor ) ) );
5026  }
5027  // loop on nodes, i.e. loop on components of this block
5028  for ( i = 0; i < fe.nbFEDof(); i++ )
5029  {
5030  for ( j = 0; j < fe.nbFEDof(); j++ )
5031  {
5032 
5033  // loop on quadrature points
5034  s = 0.;
5035  g = 0.;
5036 
5037  for ( ig = 0; ig < fe.nbQuadPt(); ig++ )
5038  {
5039  // - convect^T [I\div d - (\grad d)^T] (u^k)^T :\grad\phi_i
5040  for ( jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
5041  {
5042  //source_mass1
5043  s += convect_eta[ ig ][ jcoor ][j][kcoor] * guk[ ig ][ icoor ][jcoor] * fe.phi ( i, ig ) * fe.weightDet ( ig ) * rho;
5044  //source_stress
5045  s += B[ ig ][ icoor ][ jcoor ][j][kcoor] * fe.phiDer ( i, jcoor, ig ) * fe.weightDet ( ig );
5046  //source_stress2
5047  s -= fe.phiDer ( i, jcoor, ig ) * A2[ ig ][ icoor ][ jcoor ][j][kcoor] * fe.weightDet ( ig ) * mu;
5048 
5049 
5050  }
5051  //source_mass1
5052  s += auxMass[ ig ][j][kcoor] * uk[ ig ][ icoor ] * fe.phi ( i, ig ) * fe.weightDet ( ig ) * rho;
5053 
5054  //source_mass3
5055  s += 0.5 * auxMass3[ ig ][j][kcoor] * uk[ ig ][ icoor ] * fe.phi ( i, ig ) * fe.weightDet ( ig ) * rho;
5056 
5057  if (wImplicit)
5058  {
5059  //source_mass2{
5060  s -= aux[ ig ][ icoor ][j][kcoor] * fe.phi ( i, ig ) * fe.weightDet ( ig ) * rho * alpha;
5061  //convective term
5062  g += aux[ ig ][ icoor ][j][kcoor] * fe.phi ( i, ig ) * fe.weightDet ( ig ) * rho;
5063  }
5064  }
5065  mat ( i , j ) += s;
5066  if (wImplicit && mat_convect.get() )
5067  {
5068  (*mat_convect) ( i , j ) += g;
5069  }
5070  }
5071  }
5072  }
5073  }
5074 }
5075 
5076 
5077 //----------------------------------------------------------------------
5078 // Compute the gradient in the Hdiv space, i.e. the opposite and transpose of the divergence matrix.
5079 void grad_Hdiv ( Real coef, MatrixElemental& elmat, const CurrentFE& dualFE, const CurrentFE& primalFE, int iblock, int jblock )
5080 {
5081  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
5082  Real sumdivphi (0.);
5083  const QuadratureRule& dualQuadRule ( dualFE.quadRule() );
5084 
5085  // Loop over all the degrees of freedom of the dual variable.
5086  for ( UInt i (0); i < dualFE.nbFEDof(); ++i )
5087  {
5088  // Loop over all the degrees of freedom of the primal variable.
5089  for ( UInt j (0); j < primalFE.nbFEDof(); ++j )
5090  {
5091  sumdivphi = 0.;
5092  // Loop over all the quadrature points.
5093  for ( UInt ig (0); ig < dualFE.nbQuadPt(); ++ig )
5094  {
5095  // There is no jacobian because of property of Piola transformation.
5096  sumdivphi -= primalFE.phi ( j, ig ) * dualFE.divPhiRef ( i, ig ) * dualQuadRule.weight ( ig );
5097  }
5098  mat ( i, j ) += coef * sumdivphi;
5099  }
5100  }
5101 }
5102 
5103 //----------------------------------------------------------------------
5104 // Compute the divergence in the Hdiv space.
5105 void div_Hdiv ( Real coef, MatrixElemental& elmat, const CurrentFE& dualFE, const CurrentFE& primalFE, int iblock, int jblock )
5106 {
5107  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
5108  Real sumdivphi (0.);
5109  const QuadratureRule& dualQuadRule ( dualFE.quadRule() );
5110 
5111  // Loop over all the degrees of freedom of the dual variable.
5112  for ( UInt i (0); i < dualFE.nbFEDof(); ++i )
5113  {
5114  // Loop over all the degrees of freedom of the primal variable.
5115  for ( UInt j (0); j < primalFE.nbFEDof(); ++j )
5116  {
5117  sumdivphi = 0.;
5118  // Loop over all the quadrature points.
5119  for ( UInt ig (0); ig < dualFE.nbQuadPt(); ++ig )
5120  {
5121  // There is no jacobian because of property of Piola transformation.
5122  sumdivphi += primalFE.phi ( j, ig ) * dualFE.divPhiRef ( i, ig ) * dualQuadRule.weight ( ig );
5123  }
5124  mat ( j, i ) += coef * sumdivphi;
5125  }
5126  }
5127 }
5128 
5129 //----------------------------------------------------------------------
5130 // Compute a Hdiv function dot product with the outwart unit normal times a hybrid function.
5131 void TP_VdotN_Hdiv ( Real coef, MatrixElemental& elmat, const ReferenceFEHybrid& hybridFE,
5132  const ReferenceFEHybrid& dualDotNFE, int iblock, int jblock )
5133 {
5134  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
5135  UInt nbnode;
5136  Real sum (0.);
5137 
5138  // Loop over all the staticBdFE.
5139  for ( UInt nf (0); nf < hybridFE.numberBoundaryFE(); ++nf )
5140  {
5141  // Take the staticBdFE of the hybrid finite element.
5142  const CurrentFEManifold& boundaryElementHybridFE ( hybridFE[ nf ] );
5143  // Take the staticBdFE of the Hdiv function dot product outward unit normal.
5144  const CurrentFEManifold& boundaryElementDualDotNFE ( dualDotNFE[ nf ] );
5145  nbnode = boundaryElementHybridFE.nbFEDof();
5146 
5147  // Loop over all the the degrees of freedom of the dual dot normal variable.
5148  for ( UInt i (0); i < nbnode; ++i )
5149  {
5150  // Loop over all the degrees of freedom of the hybrid variable.
5151  for ( UInt j (0); j < nbnode; ++j )
5152  {
5153  sum = 0.;
5154  // Loop over all the quadrature point.
5155  for ( UInt ig (0); ig < boundaryElementHybridFE.nbQuadPt(); ++ig )
5156  {
5157  // Using the Piola transform properties.
5158  sum += boundaryElementHybridFE.phi ( j , ig ) *
5159  boundaryElementDualDotNFE.phi ( i , ig ) *
5160  boundaryElementHybridFE.wRootDetMetric ( ig );
5161  }
5162  // The matrix is block diagonal, so the size of the blocks is bdfe.nbNode.
5163  mat ( nf * nbnode + i, nf * nbnode + j ) += sum * coef;
5164  }
5165  }
5166  }
5167 }
5168 
5169 //----------------------------------------------------------------------
5170 // Compute the mass matrix for hybrid variable.
5171 void TP_TP_Hdiv ( Real coef, MatrixElemental& elmat, const ReferenceFEHybrid& hybridFE, int iblock, int jblock )
5172 {
5173  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
5174  UInt nbnode;
5175  Real sum (0.);
5176 
5177  // Loop over all the staticBdFE.
5178  for ( UInt nf (0); nf < hybridFE.numberBoundaryFE(); ++nf )
5179  {
5180  // Take the staticBdFE of the hybrid finite element.
5181  const CurrentFEManifold& boundaryElementHybridFE ( hybridFE[ nf ] );
5182  nbnode = boundaryElementHybridFE.nbFEDof();
5183 
5184  // Loop over all the degrees of freedom of the first hybrid variable.
5185  for ( UInt i (0); i < nbnode; ++i )
5186  {
5187  // Loop over all the the degrees of freedom of the second hybrid variable.
5188  for ( UInt j (0); j < nbnode; ++j )
5189  {
5190  sum = 0.;
5191  // Loop over all the quadrature point.
5192  for ( UInt ig (0); ig < boundaryElementHybridFE.nbQuadPt() ; ++ig )
5193  // Using the Piola transform properties.
5194  sum += boundaryElementHybridFE.phi ( j , ig ) *
5195  boundaryElementHybridFE.phi ( i , ig ) *
5196  boundaryElementHybridFE.wRootDetMetric ( ig );
5197 
5198  // The matrix is block diagonal, so the size of the blocks is bdfe.nbNode.
5199  mat ( nf * nbnode + i, nf * nbnode + j ) += sum * coef;
5200  }
5201  }
5202  }
5203 }
5204 
5205 
5206 //----------------------------------------------------------------------
5207 // Compute the mass matrix in Hdiv with a real scalar coefficient.
5208 void mass_Hdiv ( Real coef, MatrixElemental& elmat, const CurrentFE& dualFE, int iblock, int jblock )
5209 {
5210  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
5211  Real sum (0.);
5212 
5213  // Loop over all the degrees of freedom of the first dual variable.
5214  for ( UInt j (0); j < dualFE.nbFEDof(); ++j )
5215  {
5216  // Loop over all the degrees of freedom of the second dual variable.
5217  for ( UInt i (0); i < dualFE.nbFEDof() /* by symmetry j+1 */; ++i )
5218  {
5219  sum = 0.;
5220  // Loop over all the quadrature point.
5221  for ( UInt ig (0); ig < dualFE.nbQuadPt(); ++ig )
5222  {
5223  // Loop over all the space dimension, e.g. in 3D three times, in 2D two times.
5224  for ( UInt icoor (0); icoor < dualFE.nbLocalCoor() ; ++icoor )
5225  {
5226  sum += dualFE.phi ( j, icoor, ig ) * dualFE.phi ( i, icoor, ig ) * dualFE.wDetJacobian ( ig );
5227  }
5228  }
5229  // Beware coef is the inverse of the permeability, not the permeability.
5230  mat ( i, j ) += sum * coef;
5231  }
5232  }
5233 }
5234 
5235 //----------------------------------------------------------------------
5236 // Compute the mass matrix in Hdiv with a real matrix coefficient.
5237 void mass_Hdiv ( Matrix const& Invperm, MatrixElemental& elmat, const CurrentFE& dualFE, int iblock, int jblock )
5238 {
5239  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
5240  Real partialSum (0.), sum (0.);
5241 
5242  // Loop over all the degrees of freedom of the first dual variable.
5243  for ( UInt j (0); j < dualFE.nbFEDof(); ++j )
5244  {
5245  // Loop over all the degrees of freedom of the second dual variable.
5246  for ( UInt i (0); i < dualFE.nbFEDof() /* by symmetry j+1 */; ++i )
5247  {
5248  sum = 0.;
5249  // Loop over all the quadrature point.
5250  for ( UInt ig (0); ig < dualFE.nbQuadPt(); ++ig )
5251  {
5252  // partialSum = phi[icoor]^T * K^-1 * phi[jcoor] for all icorr and jcoor.
5253  partialSum = 0.;
5254  /* Loop over all the space dimension of the first dual variable,
5255  e.g. in 3D three times, in 2D two times.*/
5256  for ( UInt icoor (0); icoor < dualFE.nbLocalCoor(); ++icoor )
5257  {
5258  /* Loop over all the space dimension of the first dual variable,
5259  e.g. in 3D three times, in 2D two times.*/
5260  for ( UInt jcoor (0); jcoor < dualFE.nbLocalCoor(); ++jcoor )
5261  {
5262  // Invperm is the inverse of the permeability.
5263  partialSum += ( Invperm ( icoor, jcoor ) *
5264  dualFE.phi ( j, jcoor, ig ) *
5265  dualFE.phi ( i, icoor, ig ) );
5266  }
5267  }
5268  sum += partialSum * dualFE.wDetJacobian ( ig );
5269  }
5270  mat ( i, j ) += sum ;
5271  }
5272  }
5273 }
5274 
5275 //----------------------------------------------------------------------
5276 // Compute the mass matrix in Hdiv with a real function coefficient.
5277 void mass_Hdiv ( Real ( *InvpermFun ) ( const Real&, const Real&, const Real& ),
5278  MatrixElemental& elmat, const CurrentFE& dualFE, int iblock, int jblock )
5279 {
5280  MatrixElemental::matrix_view mat = elmat.block ( iblock, jblock );
5281  Real sum (0.), x (0.), y (0.), z (0.);
5282 
5283  // Loop over all the degrees of freedom of the first dual variable.
5284  for ( UInt j (0); j < dualFE.nbFEDof(); ++j )
5285  {
5286  // Loop over all the degrees of freedom of the second dual variable.
5287  for ( UInt i (0); i < dualFE.nbFEDof() /* by symmetry j+1 */; ++i )
5288  {
5289  sum = 0.;
5290  // Loop over all the quadrature point.
5291  for ( UInt ig (0); ig < dualFE.nbQuadPt(); ++ig )
5292  {
5293  // Get the coordinate in the current element of the current quadrature point.
5294  x = dualFE.quadNode (ig, 0);
5295  y = dualFE.quadNode (ig, 1);
5296  z = dualFE.quadNode (ig, 2);
5297 
5298  // Loop over all the space dimension, e.g. in 3D three times, in 2D two times.
5299  for ( UInt icoor (0); icoor < dualFE.nbLocalCoor(); ++icoor )
5300  {
5301  // Caution Invperm is the inverse of the permeability.
5302  sum += InvpermFun ( x, y, z ) *
5303  dualFE.phi ( j, icoor, ig ) *
5304  dualFE.phi ( i, icoor, ig ) *
5305  dualFE.wDetJacobian ( ig );
5306  }
5307  }
5308  mat ( i, j ) += sum;
5309  }
5310  }
5311 }
5312 
5313 //----------------------------------------------------------------------
5314 // Compute the source vector in Hdiv with a real vector.
5315 void source_Hdiv ( const Vector& source, VectorElemental& elvec, const CurrentFE& dualFE, int iblock )
5316 {
5317  VectorElemental::vector_view vec = elvec.block ( iblock );
5318  Real sum (0.);
5319 
5320  // Loop over all the degrees of freedom of the dual variable.
5321  for ( UInt i (0); i < dualFE.nbFEDof(); ++i )
5322  {
5323  // Loop over all the quadrature point.
5324  for ( UInt ig (0); ig < dualFE.nbQuadPt(); ++ig )
5325  {
5326  sum = 0.;
5327  /* Loop over all the space dimension of the dual variable,
5328  e.g. in 3D three times, in 2D two times.*/
5329  for ( UInt icoor (0); icoor < dualFE.nbLocalCoor(); ++icoor )
5330  {
5331  // sum[icoor] = source[icoor] * phi[icoor] for all icorr.
5332  sum += source ( icoor ) * dualFE.phi ( i, icoor, ig );
5333 
5334  }
5335 
5336  vec ( i ) += sum * dualFE.wDetJacobian ( ig );
5337  }
5338  }
5339 }
5340 
5341 } // Namespace AssemblyElemental
5342 
5343 } // namespace LifeV
5344 
5345 #endif
LifeV::AssemblyElemental::source_gradpv
void source_gradpv(Real alpha, VectorElemental &pLoc, VectorElemental &elvec, const CurrentFE &fe_p, const CurrentFE &fe_u, int iblock)
Definition: AssemblyElemental.cpp:3583
LifeV::AssemblyElemental::grad_ss
void grad_ss(const int icoor, const VectorElemental &vec_loc, MatrixElemental &elmat, const CurrentFE &fe1, const CurrentFE &fe2, int iblock, int jblock)
Convective term with a local vector coefficient for Navier-Stokes problem in Skew-Symmetric form...
Definition: AssemblyElemental.cpp:3293
LifeV::AssemblyElemental::advection
void advection(Real coef, VectorElemental &vel, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock, int nb)
Definition: AssemblyElemental.cpp:3194
LifeV::QuadratureRule::weight
const Real & weight(const UInt &ig) const
weight(ig) is the ig-th quadrature weight
Definition: QuadratureRule.hpp:284
LifeV::CurrentFE::patternFirst
int patternFirst(int i) const
patternFirst(i): row index in the element matrix of the i-th term of the pattern
Definition: CurrentFE.hpp:716
LifeV::AssemblyElemental::stiff_curl
void stiff_curl(Real coef, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock, int)
Definition: AssemblyElemental.cpp:2077
LifeV::AssemblyElemental::mass
void mass(Real coef, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock)
coef(t,x,y,z,u)
Definition: AssemblyElemental.cpp:362
LifeV::AssemblyElemental::stiff_gradgrad
void stiff_gradgrad(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:635
LifeV::AssemblyElemental::mass
void mass(MatrixElemental &localMass, const CurrentFE &massCFE, const Real &coefficient, const UInt &fieldDim)
Elementary mass for constant mass coefficient.
Definition: AssemblyElemental.cpp:46
LifeV::AssemblyElemental::stiff
void stiff(Real coef, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock, int nb)
Definition: AssemblyElemental.cpp:1942
LifeV::AssemblyElemental::ipstab_bagrad
void ipstab_bagrad(const Real coef, MatrixElemental &elmat, const CurrentFE &fe1, const CurrentFE &fe2, const CurrentFE &fe3, const VectorElemental &beta, const CurrentFEManifold &bdfe, int iblock, int jblock)
p1 lives in fe1 q2 lives in fe2 beta lives in fe3
Definition: AssemblyElemental.cpp:1727
LifeV::AssemblyElemental::cholsl
void cholsl(KNM< Real > &a, KN< Real > &p, KN< Real > &b, KN< Real > &x)
Definition: AssemblyElemental.cpp:4497
LifeV::AssemblyElemental::source_advection
void source_advection(const Real &coefficient, const VectorElemental &beta_loc, const VectorElemental &uk_loc, VectorElemental &elvec, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:3635
LifeV::AssemblyElemental::ipstab_grad
void ipstab_grad(const Real coef, MatrixElemental &elmat, const CurrentFE &fe1, const CurrentFE &fe2, const CurrentFEManifold &bdfe, int iblock, int jblock, int nb)
Definition: AssemblyElemental.cpp:1272
LifeV::AssemblyElemental::stiff_gradgradTr_gradbis_2
void stiff_gradgradTr_gradbis_2(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:1038
LifeV::AssemblyElemental::ipstab_bagrad
void ipstab_bagrad(const Real coef, MatrixElemental &elmat, const CurrentFE &fe1, const CurrentFE &fe2, const VectorElemental &beta, const CurrentFEManifold &bdfe, int iblock, int jblock)
Definition: AssemblyElemental.cpp:1604
LifeV::CurrentFEManifold
A class for a finite element on a manifold.
Definition: CurrentFEManifold.hpp:86
LifeV::AssemblyElemental::source_fhn
void source_fhn(Real coef_f, Real coef_a, VectorElemental &u, VectorElemental &elvec, const CurrentFE &fe, int fblock, int eblock)
Definition: AssemblyElemental.cpp:3604
LifeV::AssemblyElemental::mass
void mass(const std::vector< Real > &coef, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock, UInt nb)
Mass term with coefficients given for each quadrature point.
Definition: AssemblyElemental.cpp:467
LifeV::AssemblyElemental::TP_VdotN_Hdiv
void TP_VdotN_Hdiv(Real coef, MatrixElemental &elmat, const ReferenceFEHybrid &hybridFE, const ReferenceFEHybrid &dualDotNFE, int iblock, int jblock)
Definition: AssemblyElemental.cpp:5131
LifeV::AssemblyElemental::source_mass2
void source_mass2(Real coef, const VectorElemental &uk_loc, const VectorElemental &dw_loc, VectorElemental &elvec, const CurrentFE &fe)
for Newton FSI
Definition: AssemblyElemental.cpp:3871
LifeV::AssemblyElemental::source_divuq
void source_divuq(Real alpha, VectorElemental &uLoc, VectorElemental &elvec, const CurrentFE &fe_u, const CurrentFE &fe_p, int iblock)
Definition: AssemblyElemental.cpp:3562
LifeV::AssemblyElemental::stiff_divgrad_2
void stiff_divgrad_2(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:582
LifeV::AssemblyElemental::stiff_sd
void stiff_sd(Real coef, const VectorElemental &vec_loc, MatrixElemental &elmat, const CurrentFE &fe, const CurrentFE &fe2, int iblock, int jblock, int nb)
StreamLine Diffusion.
Definition: AssemblyElemental.cpp:2980
LifeV::AssemblyElemental::stiff_dergrad
void stiff_dergrad(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
for Newton on St-Venant
Definition: AssemblyElemental.cpp:2600
LifeV::AssemblyElemental::stiff_div
void stiff_div(Real coef, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:2503
LifeV::CurrentFE::phi
Real phi(UInt node, UInt component, UInt quadNode) const
Getter for basis function values (vectorial FE case)
Definition: CurrentFE.hpp:479
LifeV::AssemblyElemental::stiff_dergradbis
void stiff_dergradbis(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:2539
LifeV::AssemblyElemental::grad
void grad(const int &icoor, const std::vector< Real > &localVector, MatrixElemental &elmat, const CurrentFE &currentFE1, const CurrentFE &currentFE2, const int &iblock, const int &jblock)
Conective term with a local vector given by quadrature node.
Definition: AssemblyElemental.cpp:3388
LifeV::AssemblyElemental::stiff_dergrad_gradbis_Tr_2
void stiff_dergrad_gradbis_Tr_2(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:917
LifeV::AssemblyElemental::advectionNewton
void advectionNewton(Real coef, VectorElemental &vel, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock)
Assemble the term .
Definition: AssemblyElemental.cpp:173
LifeV::AssemblyElemental::stiff
void stiff(const std::vector< Real > &coef, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock, int nb)
Stiff term with coefficient given for each quadrature node.
Definition: AssemblyElemental.cpp:2009
LifeV::AssemblyElemental::shape_terms
void shape_terms(Real rho, Real mu, const VectorElemental &un_loc, const VectorElemental &uk_loc, const VectorElemental &wk_loc, const VectorElemental &convect_loc, const VectorElemental &pk_loc, MatrixElemental &elmat, const CurrentFE &fe, const CurrentFE &fe_p, ID, MatrixElemental &, int, bool wImplicit, Real alpha, std::shared_ptr< MatrixElemental > elmat_convect)
Shape terms for the CE system in FSI Newton.
Definition: AssemblyElemental.cpp:4714
LifeV::AssemblyElemental::ipstab_grad
void ipstab_grad(const Real coef, MatrixElemental &elmat, const CurrentFE &fe1, const CurrentFE &fe2, const CurrentFEManifold &bdfe, int iblock, int jblock)
Definition: AssemblyElemental.cpp:1171
LifeV::AssemblyElemental::mass_Hdiv
void mass_Hdiv(Real coef, MatrixElemental &elmat, const CurrentFE &dualFE, int iblock, int jblock)
Definition: AssemblyElemental.cpp:5208
LifeV::AssemblyElemental::stiff_gradgradTr_gradbis_3
void stiff_gradgradTr_gradbis_3(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:1099
LifeV::CurrentFE::nbUpper
UInt nbUpper() const
Getter for the number of entries in the pattern.
Definition: CurrentFE.hpp:413
LifeV::CurrentFE::nbDiag
UInt nbDiag() const
Getter for the diagonal entries in the pattern.
Definition: CurrentFE.hpp:407
LifeV::AssemblyElemental::stab_stokes
void stab_stokes(Real visc, Real coef_stab, MatrixElemental &elmat, const CurrentFE &fe, int block_pres)
Definition: AssemblyElemental.cpp:3168
LifeV::AssemblyElemental::mass_divw
void mass_divw(Real coef, const VectorElemental &w_loc, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock, UInt nb)
Definition: AssemblyElemental.cpp:2772
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::AssemblyElemental::source
void source(Real coef, VectorElemental &f, VectorElemental &elvec, const CurrentFE &fe, int fblock, int eblock)
Definition: AssemblyElemental.cpp:3499
LifeV::AssemblyElemental::source_Hdiv
void source_Hdiv(const Vector &source, VectorElemental &elvec, const CurrentFE &dualFE, int iblock)
Definition: AssemblyElemental.cpp:5315
LifeV::AssemblyElemental::mass_Hdiv
void mass_Hdiv(Matrix const &Invperm, MatrixElemental &elmat, const CurrentFE &dualFE, int iblock, int jblock)
Definition: AssemblyElemental.cpp:5237
LifeV::AssemblyElemental::stiff_dergrad_gradbis_2
void stiff_dergrad_gradbis_2(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:858
LifeV::AssemblyElemental::source_mass1
void source_mass1(Real coef, const VectorElemental &uk_loc, const VectorElemental &wk_loc, const VectorElemental &convect_loc, const VectorElemental &d_loc, VectorElemental &elvec, const CurrentFE &fe)
for Newton FSI
Definition: AssemblyElemental.cpp:3722
LifeV::AssemblyElemental::grad_div
void grad_div(Real coef_grad, Real coef_div, MatrixElemental &elmat, const CurrentFE &fe_u, const CurrentFE &fe_p, int block_pres)
Definition: AssemblyElemental.cpp:3139
LifeV::AssemblyElemental::grad
void grad(const int icoor, const VectorElemental &vec_loc, MatrixElemental &elmat, const CurrentFE &fe1, const CurrentFE &fe2, int iblock, int jblock)
Convective term with a local vector coefficient (useful for Navier-Stokes problem) ...
Definition: AssemblyElemental.cpp:3250
LifeV::CurrentFE::nbLocalCoor
UInt nbLocalCoor() const
Getter for the number of local coordinates.
Definition: CurrentFE.hpp:425
LifeV::AssemblyElemental::ipstab_bgrad
void ipstab_bgrad(const Real coef, MatrixElemental &elmat, const CurrentFE &fe1, const CurrentFE &fe2, const VectorElemental &beta, const CurrentFEManifold &bdfe, int iblock, int jblock, int nb)
Definition: AssemblyElemental.cpp:1379
LifeV::AssemblyElemental::bodyForces
void bodyForces(VectorElemental &localForce, const CurrentFE &massRhsCFE, const function_Type &fun, const Real &t, const UInt &fieldDim)
Definition: AssemblyElemental.cpp:303
LifeV::CurrentFE::phiDer
Real phiDer(UInt node, UInt derivative, UInt quadNode) const
Old accessor, use dphi instead.
Definition: CurrentFE.hpp:548
LifeV::CurrentFE::dphi
Real dphi(UInt node, UInt derivative, UInt quadNode) const
Getter for the derivatives of the basis functions.
Definition: CurrentFE.hpp:486
LifeV::ID
uint32_type ID
IDs.
Definition: LifeV.hpp:194
LifeV::AssemblyElemental::mass
void mass(Real coef, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock, UInt nb)
Definition: AssemblyElemental.cpp:405
LifeV::AssemblyElemental::stiff_gradgrad_2
void stiff_gradgrad_2(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:687
LifeV::AssemblyElemental::source_stiff
void source_stiff(const std::vector< Real > &constant, VectorElemental &elvec, const CurrentFE &currentFe, const int &iblock)
Assembly for the source term where is a given by the values in the quadrature nodes.
Definition: AssemblyElemental.cpp:3542
LifeV::AssemblyElemental::grad
void grad(const int icoor, Real coef, MatrixElemental &elmat, const CurrentFE &fe_u, const CurrentFE &fe_p, int iblock, int jblock)
Definition: AssemblyElemental.cpp:3085
LifeV::CurrentFE::quadRule
const QuadratureRule & quadRule() const
Getter for the quadrature rule.
Definition: CurrentFE.hpp:395
LifeV::AssemblyElemental::stiffStrain
void stiffStrain(MatrixElemental &localStiff, const CurrentFE &stiffCFE, const Real &coefficient, const UInt &fieldDim)
Definition: AssemblyElemental.cpp:267
LifeV::AssemblyElemental::stiff_strain
void stiff_strain(Real coef, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:2722
LifeV::AssemblyElemental::stiff_derdiv
void stiff_derdiv(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
for Newton on St-Venant
Definition: AssemblyElemental.cpp:2665
LifeV::AssemblyElemental::function_Type
std::function< const Real(const Real &, const Real &, const Real &, const Real &, const ID &) > function_Type
Use the portable syntax of the boost function.
Definition: AssemblyElemental.hpp:84
LifeV::AssemblyElemental::TP_TP_Hdiv
void TP_TP_Hdiv(Real coef, MatrixElemental &elmat, const ReferenceFEHybrid &hybridFE, int iblock, int jblock)
Definition: AssemblyElemental.cpp:5171
LifeV::AssemblyElemental::mass_Hdiv
void mass_Hdiv(Real(*InvpermFun)(const Real &, const Real &, const Real &), MatrixElemental &elmat, const CurrentFE &dualFE, int iblock, int jblock)
Definition: AssemblyElemental.cpp:5277
LifeV::CurrentFE::diameter
Real diameter() const
return the diameter of the element in the 1-norm
Definition: CurrentFE.cpp:595
LifeV::CurrentFE::quadNode
Real quadNode(UInt node, UInt coordinate) const
Getter for the quadrature nodes.
Definition: CurrentFE.hpp:444
LifeV::AssemblyElemental::stiff_gradgradTr_gradbis
void stiff_gradgradTr_gradbis(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:976
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
LifeV::CurrentFE
CurrentFE - A primordial class for the assembly of the local matrices/vectors retaining the values on...
Definition: CurrentFE.hpp:243
LifeV::AssemblyElemental::ipstab_div
void ipstab_div(const Real coef, MatrixElemental &elmat, const CurrentFE &fe1, const CurrentFE &fe2, const CurrentFEManifold &bdfe, int iblock, int jblock)
Definition: AssemblyElemental.cpp:1513
LifeV::AssemblyElemental::source_stress2
void source_stress2(Real coef, const VectorElemental &uk_loc, const VectorElemental &d_loc, VectorElemental &elvec, const CurrentFE &fe_u)
Definition: AssemblyElemental.cpp:4194
LifeV::ReferenceFEHybrid::operator[]
const CurrentFEManifold & operator[](const ID &i) const
Extracting a CurrentFEManifold from the faces list.
Definition: ReferenceFEHybrid.hpp:164
LifeV::CurrentFE::nbFEDof
UInt nbFEDof() const
Getter for the number of nodes.
Definition: CurrentFE.hpp:377
LifeV::AssemblyElemental::divergence
void divergence(MatrixElemental &elmat, const CurrentFE &uCFE, const CurrentFE &pCFE, const UInt &fieldDim, const Real &coefficient)
Definition: AssemblyElemental.cpp:235
LifeV::AssemblyElemental::div
void div(const int icoor, Real coef, MatrixElemental &elmat, const CurrentFE &fe_u, const CurrentFE &fe_p, int iblock, int jblock)
Definition: AssemblyElemental.cpp:3114
LifeV::CurrentFE::wDetJacobian
Real wDetJacobian(UInt quadNode) const
Getter for the weighted jacobian determinant.
Definition: CurrentFE.hpp:458
LifeV::AssemblyElemental::choldc
void choldc(KNM< Real > &a, KN< Real > &p)
Definition: AssemblyElemental.cpp:4470
LifeV::CurrentFE::phi
Real phi(UInt node, UInt quadNode) const
Getter for basis function values (scalar FE case)
Definition: CurrentFE.hpp:472
LifeV::AssemblyElemental::source_press
void source_press(Real coef, const VectorElemental &uk_loc, const VectorElemental &d_loc, VectorElemental &elvec, const CurrentFE &fe_u, const CurrentFE &fe_p, int iblock)
for Newton FSI
Definition: AssemblyElemental.cpp:4290
LifeV::AssemblyElemental::stiff_dergrad_gradbis_Tr
void stiff_dergrad_gradbis_Tr(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:800
LifeV::AssemblyElemental::source_press2
void source_press2(Real coef, const VectorElemental &p_loc, const VectorElemental &d_loc, VectorElemental &elvec, const CurrentFE &fe, int iblock)
Definition: AssemblyElemental.cpp:4371
LifeV::AssemblyElemental::stiff
void stiff(Real coef, Real(*fct)(Real, Real, Real), MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock)
Definition: AssemblyElemental.cpp:1894
LifeV::AssemblyElemental::mass_divw
void mass_divw(const std::vector< Real > &coef, const VectorElemental &w_loc, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock, UInt nb)
Idem mass_divw , but with coefficient given by quadrature node.
Definition: AssemblyElemental.cpp:2846
LifeV::AssemblyElemental::stiffness
void stiffness(MatrixElemental &localStiff, const CurrentFE &stiffCFE, const Real &coefficient, const UInt &fieldDim)
Elementary stiffness for constant coefficient.
Definition: AssemblyElemental.cpp:106
LifeV::AssemblyElemental::stiff_divgrad
void stiff_divgrad(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:530
LifeV::CurrentFE::patternSecond
int patternSecond(int i) const
patternSecond(i): column index in the element matrix of the i-th term of the pattern ...
Definition: CurrentFE.hpp:721
LifeV::AssemblyElemental::grad
void grad(const int icoor, const VectorElemental &vec_loc, MatrixElemental &elmat, const CurrentFE &fe1, const CurrentFE &fe2, const CurrentFE &fe3, int iblock, int jblock)
Convective term with a local vector coefficient (useful for Navier-Stokes problem+adv-diff) ...
Definition: AssemblyElemental.cpp:3340
LifeV::CurrentFE::weightDet
Real weightDet(UInt quadNode) const
Old accessor, use wDetJacobian instead.
Definition: CurrentFE.hpp:527
LifeV::QuadratureRule
QuadratureRule - The basis class for storing and accessing quadrature rules.
Definition: QuadratureRule.hpp:111
LifeV::AssemblyElemental::source_mass
void source_mass(const std::vector< Real > &constant, VectorElemental &elvec, const CurrentFE &currentFe, const int &iblock)
Assembly for the source term where is a given by the values in the quadrature nodes.
Definition: AssemblyElemental.cpp:3528
LifeV::AssemblyElemental::source_stress
void source_stress(Real coef, Real mu, const VectorElemental &uk_loc, const VectorElemental &pk_loc, const VectorElemental &d_loc, VectorElemental &elvec, const CurrentFE &fe_u, const CurrentFE &fe_p)
for Newton FSI
Definition: AssemblyElemental.cpp:4063
LifeV::AssemblyElemental::grad
void grad(MatrixElemental &elmat, const CurrentFE &uCFE, const CurrentFE &pCFE, const UInt &fieldDim)
Definition: AssemblyElemental.cpp:204
LifeV::CurrentFE::divPhiRef
Real divPhiRef(UInt node, UInt quadNode) const
Getter for the divergence of a vectorial FE in the reference frame.
Definition: CurrentFE.hpp:500
LifeV::AssemblyElemental::mass_gradu
void mass_gradu(Real coef, const VectorElemental &u0_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:2921
LifeV::AssemblyElemental::grad_Hdiv
void grad_Hdiv(Real coef, MatrixElemental &elmat, const CurrentFE &dualFE, const CurrentFE &primalFE, int iblock, int jblock)
Definition: AssemblyElemental.cpp:5079
LifeV::AssemblyElemental::div_Hdiv
void div_Hdiv(Real coef, MatrixElemental &elmat, const CurrentFE &dualFE, const CurrentFE &primalFE, int iblock, int jblock)
Definition: AssemblyElemental.cpp:5105
LifeV::AssemblyElemental::source
void source(Real constant, VectorElemental &elvec, const CurrentFE &fe, int iblock)
Definition: AssemblyElemental.cpp:3483
LifeV::CurrentFE::nbQuadPt
UInt nbQuadPt() const
Getter for the number of quadrature nodes.
Definition: CurrentFE.hpp:365
LifeV::AssemblyElemental::stiff
void stiff(Real coef, MatrixElemental &elmat, const CurrentFE &fe, int iblock, int jblock)
Definition: AssemblyElemental.cpp:1848
LifeV::UInt
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191
LifeV::AssemblyElemental::source_mass3
void source_mass3(Real coef, const VectorElemental &un_loc, const VectorElemental &uk_loc, const VectorElemental &d_loc, VectorElemental &elvec, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:3956
LifeV::ReferenceFEHybrid::numberBoundaryFE
const UInt & numberBoundaryFE() const
Return the number of boundary elements of the reference element.
Definition: ReferenceFEHybrid.hpp:177
LifeV::Matrix
boost::numeric::ublas::matrix< Real > Matrix
Definition: AssemblyElemental.hpp:66
LifeV::CurrentFE::coorMap
void coorMap(Real &x, Real &y, Real &z, Real xi, Real eta, Real zeta) const
Definition: CurrentFE.cpp:648
LifeV::AssemblyElemental::source_press
void source_press(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe_u, const CurrentFE &fe_p, ID mmDim, int iblock)
for Newton FSI
Definition: AssemblyElemental.cpp:4524
LifeV::AssemblyElemental::stiff_dergrad_gradbis
void stiff_dergrad_gradbis(Real coef, const VectorElemental &uk_loc, MatrixElemental &elmat, const CurrentFE &fe)
Definition: AssemblyElemental.cpp:742
LifeV::CurrentFE::quadPt
Real quadPt(UInt node, UInt coordinate) const
Old accessor, use quadNode instead.
Definition: CurrentFE.hpp:520