LifeV
fsi_blocks/solver/AssemblyElementalStructure.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 for the structural problem
30 
31  @author Gianmarco Mengaldo <gianmarco.mengaldo@gmail.com>
32  @author Paolo Tricerri <gianmarco.mengaldo@gmail.com>
33  @mantainer Paolo Tricerri <paolo.tricerri@epfl.ch>
34  */
35 
36 #ifndef ELEMOPERSTRUCTURE_CPP
37 #define ELEMOPERSTRUCTURE_CPP 1
38 
39 #include <lifev/fsi_blocks/solver/AssemblyElementalStructure.hpp>
40 #include <boost/multi_array.hpp>
41 
42 namespace LifeV
43 {
44 
45 namespace AssemblyElementalStructure
46 {
47 
48 void computeGradientLocalDisplacement (boost::multi_array<Real, 3>& gradientLocalDisplacement,
49  const VectorElemental& uk_loc, const CurrentFE& fe )
50 {
51  // \grad u^k at each quadrature poInt
52  Real s;
53 
54  // loop on quadrature poInts
55  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
56  {
57 
58  // loop on space coordinates
59  for ( UInt icoor = 0; icoor < nDimensions; icoor++ )
60  {
61 
62  // loop on space coordinates
63  for ( UInt jcoor = 0; jcoor < nDimensions; jcoor++ )
64  {
65  s = 0.0;
66  for (UInt i = 0; i < fe.nbFEDof(); i++ )
67  {
68  // \grad u^k at a quadrature poInt
69  s += fe.phiDer ( i, jcoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ];
70  }
71  gradientLocalDisplacement[ icoor ][ jcoor ][ ig ] = s;
72  }
73  }
74  }
75 }
76 
77 void computeLocalDeformationGradient (const VectorElemental& uk_loc, std::vector<Epetra_SerialDenseMatrix>& tensorF, const CurrentFE& fe )
78 {
79  // \grad u^k at each quadrature poInt
80  Real s;
81 
82  for ( Int k = 0; k < static_cast<Int> (fe.nbQuadPt() ); k++)
83  {
84  // loop on space coordinates
85  for ( Int icoor = 0; icoor < static_cast<Int> (nDimensions); icoor++ )
86  {
87  // loop on space coordinates
88  for ( Int jcoor = 0; jcoor < static_cast<Int> (nDimensions); jcoor++ )
89  {
90  s = 0.0;
91  for (Int i = 0; i < static_cast<Int> (fe.nbFEDof() ); i++ )
92  {
93  // \grad u^k at a quadrature poInt
94  s += fe.phiDer ( i, jcoor, k ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ];
95  }
96  tensorF[k] ( icoor, jcoor ) = s;
97 
98  if ( icoor == jcoor )
99  {
100  tensorF[k] ( icoor, jcoor ) += 1.0;
101  }
102  }
103  }
104  }
105 }
106 
107 void computeLocalDeformationGradientWithoutIdentity (const VectorElemental& uk_loc, std::vector<Epetra_SerialDenseMatrix>& tensorF, const CurrentFE& fe )
108 {
109  // \grad u^k at each quadrature poInt
110  Real s;
111 
112  for ( Int k = 0; k < static_cast<Int> (fe.nbQuadPt() ); k++)
113  {
114  // loop on space coordinates
115  for ( Int icoor = 0; icoor < static_cast<Int> (nDimensions); icoor++ )
116  {
117  // loop on space coordinates
118  for ( Int jcoor = 0; jcoor < static_cast<Int> (nDimensions); jcoor++ )
119  {
120  s = 0.0;
121  for (Int i = 0; i < static_cast<Int> (fe.nbFEDof() ); i++ )
122  {
123  // \grad u^k at a quadrature poInt
124  s += fe.phiDer ( i, jcoor, k ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ];
125  }
126  tensorF[k] ( icoor, jcoor ) = s;
127  }
128  }
129  }
130 }
131 
132 // The methods for linear elastic model (stiff_strain and stiff_div) are implemented in AssemblyElemental.cpp
133 
134 //! Methods for St. Venant Kirchhoff model
135 //! Methods for the stiffness matrix
136 
137 //! \f$ coef \cdot ( trace { [\nabla u^k]^T \nabla u }, \nabla\cdot v ) \f$
138 void stiff_derdiv ( Real coef, const boost::multi_array<Real, 3>& gradientLocalDisplacement, MatrixElemental& elmat, const CurrentFE& fe )
139 {
140  Real s;
141 
142  //
143  // blocks (icoor,jcoor) of elmat
144  //
145  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
146  {
147  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
148  {
149 
150  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
151 
152  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
153  {
154  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
155  {
156  s = 0;
157  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
158  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
159  {
160  s += fe.phiDer ( i, icoor, ig ) * gradientLocalDisplacement[ jcoor ][ k ][ ig ]
161  * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
162  }
163  mat ( i, j ) += coef * s;
164  }
165  }
166  }
167  }
168 }
169 
170 
171 
172 //! \f$ coef \cdot ( [\nabla u^k]^T \nabla u : \nabla v )\f$
173 void stiff_dergradbis ( Real coef, const boost::multi_array<Real, 3>& gradientLocalDisplacement,
174  MatrixElemental& elmat, const CurrentFE& fe )
175 {
176 
177  Real s;
178 
179  //
180  // blocks (icoor,jcoor) of elmat
181  //
182 
183  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
184  {
185  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
186  {
187 
188  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
189 
190  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
191  {
192  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
193  {
194  s = 0;
195  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
196  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
197  s += fe.phiDer ( i, k, ig ) * gradientLocalDisplacement[ jcoor][ icoor ][ ig ] *
198  fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
199  mat ( i, j ) += coef * s;
200  }
201  }
202  }
203  }
204 }
205 
206 
207 
208 //! \f$ coef * ( (\div u_k) \grad u : \grad v )
209 void stiff_divgrad ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
210 {
211 
212  Real s;
213  // local vector for \div u^k at each quadrature poInt
214  //Real duk[ fe.nbQuadPt() ];
215  std::vector<Real > duk (fe.nbQuadPt(), 0.0);
216 
217  // loop on quadrature poInts
218  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
219  {
220  s = 0;
221  // loop on space coordinates
222  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
223  {
224  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
225  {
226  // construction of \div u^k at a quadrature poInt
227  s += fe.phiDer ( i, icoor, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ];
228  }
229  }
230  duk[ ig ] = s;
231  }
232 
233  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
234 
235  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
236  {
237  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
238  {
239  s = 0.0;
240  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
241  {
242  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
243  {
244  s += duk[ ig ] * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
245  }
246  }
247  mat_tmp ( i, j ) = coef * s;
248  }
249  }
250 
251  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
252  {
253  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
254  mat += mat_tmp;
255  }
256 
257 }
258 
259 //! \f$ coef * ( \grad u_k : \grad u_k) * ( \grad u : \grad v )
260 void stiff_gradgrad ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
261 {
262  Real s, s1;
263  // (\grad u_k : \grad u_k) at each quadrature poInt
264  //Real gguk[ fe.nbQuadPt() ];
265  std::vector<Real > gguk (fe.nbQuadPt(), 0.0);
266 
267  // loop on quadrature poInts
268  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
269  {
270  s = 0;
271  // loop on space coordinates
272  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
273  {
274  for ( UInt l = 0; l < fe.nbLocalCoor(); l++ )
275  {
276  s1 = 0;
277  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
278  {
279  s1 += fe.phiDer ( i, l, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ];
280  }
281  s += s1 * s1;
282  }
283  }
284  gguk[ ig ] = s;
285  }
286 
287  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
288 
289  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
290  {
291  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
292  {
293  s = 0.0;
294  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
295  {
296  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
297  {
298  s += gguk[ ig ] * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
299  }
300  }
301  mat_tmp ( i, j ) = coef * s;
302  }
303  }
304 
305  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
306  {
307  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
308  mat += mat_tmp;
309  }
310 }
311 
312 void stiff_dergrad_gradbis ( Real coef, const boost::multi_array<Real, 3>& gradientLocalDisplacement,
313  MatrixElemental& elmat, const CurrentFE& fe )
314 {
315 
316  Real s;
317 
318  //
319  // blocks (icoor,jcoor) of elmat
320  //
321 
322  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
323  {
324  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
325  {
326 
327  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
328 
329  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
330  {
331  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
332  {
333  s = 0.0;
334  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
335  {
336  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
337  s += gradientLocalDisplacement[ icoor ][ jcoor ][ ig ] * fe.phiDer ( i, k, ig ) *
338  fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
339  }
340  mat ( i, j ) += coef * s;
341  }
342  }
343  }
344  }
345 }
346 
347 
348 
349 // coef * ( \grad u^k [\grad u]^T : \grad v )
350 void stiff_dergrad_gradbis_Tr ( Real coef, const boost::multi_array<Real, 3>& gradientLocalDisplacement,
351  MatrixElemental& elmat, const CurrentFE& fe )
352 {
353 
354  Real s;
355  //
356  // blocks (icoor,jcoor) of elmat
357  //
358 
359  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
360  {
361  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
362  {
363 
364  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
365 
366  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
367  {
368  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
369  {
370  s = 0.0;
371  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
372  {
373  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
374  s += gradientLocalDisplacement[ icoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) *
375  fe.phiDer ( i, jcoor, ig ) * fe.weightDet ( ig );
376  }
377  mat ( i, j ) += coef * s;
378  }
379  }
380  }
381  }
382 }
383 
384 
385 
386 // coef * ( \grad u^k [\grad u^k]^T \grad u : \grad v )
387 void stiff_gradgradTr_gradbis ( Real coef, const VectorElemental& uk_loc, MatrixElemental& elmat, const CurrentFE& fe )
388 {
389 
390  Real s;
391 
392  //! \grad u^k [\grad u^k]^T at each quadrature poInt
393  //Real guk_gukT[ fe.nbLocalCoor() ][ fe.nbLocalCoor() ][ fe.nbQuadPt() ];
394  boost::multi_array<Real, 3> guk_gukT (boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
395 
396  // loop on quadrature poInts
397  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
398  {
399 
400  // loop on space coordinates
401  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
402  {
403 
404  // loop on space coordinates
405  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
406  {
407  s = 0.0;
408  for ( UInt n = 0; n < fe.nbLocalCoor(); n++ )
409  {
410  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
411  {
412  for ( UInt j = 0; j < fe.nbFEDof(); j++ )
413  {
414  //! \grad u^k [\grad u^k]^T at each quadrature poInt
415  s += fe.phiDer ( i, n, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ] *
416  fe.phiDer ( j, n, ig ) * uk_loc.vec() [ j + jcoor * fe.nbFEDof() ];
417  }
418  }
419  }
420  guk_gukT[ icoor ][ jcoor ][ ig ] = s;
421  }
422  }
423  }
424 
425  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
426  {
427  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
428  {
429 
430  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
431 
432  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
433  {
434  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
435  {
436  s = 0;
437  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
438  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
439  s += fe.phiDer ( i, k, ig ) * guk_gukT[ icoor ][ jcoor ][ ig ] *
440  fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
441  mat ( i, j ) += coef * s;
442  }
443  }
444  }
445  }
446 }
447 // End of methods for the stiffness matrix (St. Venant-Kirchhoff material)
448 
449 //! Methods for the jacobian (St. Venant-Kirchhoff material)
450 
451 //! \f$ coef \cdot ( [\nabla u]^T \nabla u^k + [\nabla u^k]^T \nabla u : \nabla v )\f$
452 void stiff_dergrad ( Real coef, const boost::multi_array<Real, 3>& gradientLocalDisplacement,
453  MatrixElemental& elmat, const CurrentFE& fe )
454 {
455 
456  Real s;
457 
458  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
459  {
460  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
461  {
462 
463  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
464 
465  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
466  {
467  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
468  {
469  s = 0;
470  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
471  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
472  {
473  s += fe.phiDer ( i, k, ig ) *
474  ( gradientLocalDisplacement[ jcoor ][ k ][ ig ] * fe.phiDer ( j, icoor, ig )
475  + gradientLocalDisplacement[ jcoor ][ icoor ][ ig ] * fe.phiDer ( j, k, ig ) ) * fe.weightDet ( ig );
476  }
477  mat ( i, j ) += coef * s;
478  }
479  }
480  }
481  }
482 }
483 
484 
485 
486 // coef * ( (\div u) \grad u_k : \grad v )
487 void stiff_divgrad_2 ( Real coef, const boost::multi_array<Real, 3>& gradientLocalDisplacement,
488  MatrixElemental& elmat, const CurrentFE& fe )
489 {
490 
491  Real s;
492 
493  //
494  // blocks (icoor,jcoor) of elmat
495  //
496 
497  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
498  {
499  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
500  {
501  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
502 
503  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
504  {
505  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
506  {
507  s = 0;
508  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
509  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
510  s += fe.phiDer ( j, jcoor, ig ) *
511  gradientLocalDisplacement[ icoor ][ k ][ ig ] *
512  fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
513  mat ( i, j ) += coef * s;
514  }
515  }
516  }
517  }
518 }
519 
520 // coef * ( \grad u_k : \grad u) *( \grad u_k : \grad v )
521 void stiff_gradgrad_2 ( Real coef, const boost::multi_array<Real, 3>& gradientLocalDisplacement, MatrixElemental& elmat, const CurrentFE& fe )
522 {
523  Real s;
524 
525  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
526  {
527  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
528  {
529  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
530 
531  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
532  {
533  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
534  {
535  s = 0.0;
536  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
537  {
538  for ( UInt l = 0; l < fe.nbLocalCoor(); ++l )
539  {
540  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
541  s += gradientLocalDisplacement[ jcoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
542  gradientLocalDisplacement[ icoor ][ k ][ ig ] * fe.phiDer (i, k, ig ) *
543  fe.weightDet ( ig );
544  }
545  }
546  mat ( i, j ) += coef * s;
547  }
548  }
549  }
550  }
551 }
552 
553 // coef * ( \grad \delta u [\grad u^k]^T : \grad v )
554 void stiff_dergrad_gradbis_2 ( Real coef, const boost::multi_array<Real, 3>& gradientLocalDisplacement,
555  MatrixElemental& elmat, const CurrentFE& fe )
556 {
557 
558  Real s;
559 
560  //
561  // blocks (icoor,jcoor) of elmat
562  //
563  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
564 
565  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
566  {
567  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
568  {
569  s = 0.0;
570  for ( UInt l = 0; l < fe.nbLocalCoor(); ++l )
571  {
572  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
573  {
574  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
575  s += gradientLocalDisplacement[ l ][ k ][ ig ] * fe.phiDer ( i, k, ig ) *
576  fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
577  }
578  }
579  mat_tmp ( i, j ) = coef * s;
580  }
581  }
582 
583  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
584  {
585  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
586  mat += mat_tmp;
587  }
588 }
589 
590 void stiff_dergrad_gradbis_Tr_2 ( Real coef, const boost::multi_array<Real, 3>& gradientLocalDisplacement,
591  MatrixElemental& elmat, const CurrentFE& fe )
592 {
593 
594  Real s;
595  //
596  // blocks (icoor,jcoor) of elmat
597  //
598  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
599 
600  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
601  {
602  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
603  {
604  s = 0.0;
605  for ( UInt l = 0; l < fe.nbLocalCoor(); ++l )
606  {
607  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
608  {
609  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
610  {
611  s += gradientLocalDisplacement[ k ][ l ][ ig ] * fe.phiDer ( i, k, ig ) *
612  fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
613  }
614  }
615  }
616  mat_tmp ( i, j ) = coef * s;
617  }
618  }
619 
620  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
621  {
622  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
623  mat += mat_tmp;
624  }
625 }
626 
627 
628 
629 // coef * ( \grad u^k [\grad u]^T \grad u^k : \grad v )
630 void stiff_gradgradTr_gradbis_2 ( Real coef, const boost::multi_array<Real, 3>& gradientLocalDisplacement,
631  MatrixElemental& elmat, const CurrentFE& fe )
632 {
633 
634  Real s;
635 
636  //
637  // blocks (icoor,jcoor) of elmat
638  //
639  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
640  {
641  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); ++jcoor )
642  {
643 
644  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor ); // it extracts the (icoor, jcoor) block
645 
646  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
647  {
648  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
649  {
650  s = 0;
651  for ( UInt l = 0; l < fe.nbLocalCoor(); ++l )
652  {
653  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
654  {
655  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
656  {
657  s += gradientLocalDisplacement[ icoor ][ l ][ ig ] *
658  gradientLocalDisplacement[ jcoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) *
659  fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
660  }
661  }
662  }
663  mat ( i, j ) += coef * s;
664  }
665  }
666  }
667  }
668 }
669 // coef * ( \grad u [\grad u^k]^T \grad u^k : \grad v )
670 void stiff_gradgradTr_gradbis_3 ( Real coef, const VectorElemental& uk_loc,
671  MatrixElemental& elmat, const CurrentFE& fe )
672 {
673 
674  Real s;
675  // \grad u^k [\grad u^k]^T at each quadrature poInt
676  //Real guk_gukT[ fe.nbLocalCoor() ][ fe.nbLocalCoor() ][ fe.nbQuadPt() ];
677  boost::multi_array<Real, 3> guk_gukT (boost::extents[fe.nbLocalCoor()][fe.nbLocalCoor()][fe.nbQuadPt()]);
678 
679  // loop on quadrature poInts // (\grad u^k [\grad u^k]^T )^T
680  for ( UInt ig = 0; ig < fe.nbQuadPt(); ig++ )
681  {
682 
683  // loop on space coordinates
684  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); icoor++ )
685  {
686 
687  // loop on space coordinates
688  for ( UInt jcoor = 0; jcoor < fe.nbLocalCoor(); jcoor++ )
689  {
690  s = 0.0;
691  for ( UInt n = 0; n < fe.nbLocalCoor(); n++ )
692  {
693  for ( UInt i = 0; i < fe.nbFEDof(); i++ )
694  {
695  for ( UInt j = 0; j < fe.nbFEDof(); j++ )
696  {
697  // \grad u^k [\grad u^k]^T at each quadrature poInt
698  s += fe.phiDer ( i, n, ig ) * uk_loc.vec() [ i + icoor * fe.nbFEDof() ] *
699  fe.phiDer ( j, n, ig ) * uk_loc.vec() [ j + jcoor * fe.nbFEDof() ] ;
700  }
701  }
702  }
703  guk_gukT[ icoor ][ jcoor ][ ig ] = s;
704  }
705  }
706  }
707 
708  //
709  // blocks (icoor,jcoor) of elmat
710  //
711  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
712 
713  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
714  {
715  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
716  {
717  s = 0.0;
718  for ( UInt l = 0; l < fe.nbLocalCoor(); ++l )
719  {
720  for ( UInt k = 0; k < fe.nbLocalCoor(); ++k )
721  {
722  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
723  {
724  s += guk_gukT[ k ][ l ][ ig ] * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
725  }
726  }
727  }
728  mat_tmp ( i, j ) = coef * s;
729  }
730  }
731 
732  for ( UInt icoor = 0; icoor < fe.nbLocalCoor(); ++icoor )
733  {
734  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
735  mat += mat_tmp;
736  }
737 }
738 // End of St. Venant Kirchhoff model
739 
740 //! ***********************************************************************************************
741 //! METHODS SHARED BETWEEN NEO-HOOKEAN AND EXPONENTIAL MODELS
742 //! ***********************************************************************************************
743 //! The volumetric part is the same for Neo-Hookean and Expoential models is the same
744 
745 //! STIFFFNESS VECTOR -----------------------------------------------------------------------------
746 //! Volumetric part--------------------------------------------------------------------------------
747 //! Source term source_Pvol: Int { coef /2* (J^2 - J + log(J) ) * 1/J * (CofF : \nabla v) }
748 void source_Pvol ( Real coef,
749  const boost::multi_array<Real, 3 >& CofFk,
750  const std::vector<Real>& Jk,
751  VectorElemental& elvec,
752  const CurrentFE& fe )
753 {
754  Real s;
755 
756  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
757  {
758  // block (icoor) of elvec
759  VectorElemental::vector_view vec = elvec.block ( icoor );
760  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
761  {
762  s = 0.0;
763  for ( UInt k = 0; k < nDimensions; ++k )
764  {
765  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
766  {
767  s += ( Jk[ig] * Jk[ig] - Jk[ig] + log ( Jk[ig] ) ) * ( 1 / Jk[ig] ) *
768  CofFk[icoor][ k][ ig] * fe.phiDer (i, k, ig) * fe.weightDet (ig);
769  }
770  }
771  vec (i) += coef * s;
772  }
773  }
774 }
775 
776 //! JACOBIAN MATRIX -------------------------------------------------------------------------------
777 //! Volumetric part--------------------------------------------------------------------------------
778 
779 //! 1. Jacobian matrix: Int { 1/2 * coef * ( 2 - 1/J + 1/J^2 ) * ( CofF : \nabla \delta ) (CofF : \nabla v) }
780 void stiff_Jac_Pvol_1term ( Real coef,
781  const boost::multi_array<Real, 3 >& CofFk,
782  const std::vector<Real>& Jk,
783  MatrixElemental& elmat,
784  const CurrentFE& fe )
785 {
786  Real s;
787 
788  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
789  {
790  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
791  {
792  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
793  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
794  {
795  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
796  {
797  s = 0.0;
798  for ( UInt l = 0; l < nDimensions; ++l )
799  {
800  for ( UInt k = 0; k < nDimensions; ++k )
801  {
802  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
803  {
804  s += ( 2.0 - ( 1 / Jk[ig] ) + ( 1 / ( Jk[ig] * Jk[ig] ) ) ) *
805  CofFk[ jcoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
806  CofFk[ icoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) *
807  fe.weightDet ( ig );
808  }
809  }
810  }
811  mat ( i, j ) += s * coef;
812  }
813  }
814  }
815  }
816 }
817 
818 
819 //! 2. Stiffness matrix: Int { 1/2 * coef * ( 1/J - 1 - log(J)/J^2 ) * ( CofF [\nabla \delta]^t CofF ) : \nabla v }
820 void stiff_Jac_Pvol_2term ( Real coef,
821  const boost::multi_array<Real, 3 >& CofFk,
822  const std::vector<Real>& Jk,
823  MatrixElemental& elmat,
824  const CurrentFE& fe )
825 {
826  Real s;
827 
828  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
829  {
830  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
831  {
832  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
833  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
834  {
835  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
836  {
837  s = 0.0;
838  for ( UInt l = 0; l < nDimensions; ++l )
839  {
840  for ( UInt k = 0; k < nDimensions; ++k )
841  {
842  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
843  {
844  s += ( ( 1 / Jk[ig] ) - 1. - ( 1 / ( Jk[ig] * Jk[ig] ) ) * log ( Jk[ig] ) ) *
845  CofFk[ icoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
846  CofFk[ jcoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
847  }
848  }
849  }
850  mat ( i, j ) += s * coef;
851  }
852  }
853  }
854  }
855 }
856 
857 
858 //! ***********************************************************************************************
859 //! METHODS FOR NEO-HOOKEAN MODEL
860 //! ***********************************************************************************************
861 //! Stiffness vector isochoric part ---------------------------------------------------------------
862 
863 //! Source term source_P1iso_NH: Int { coef * ( J^(-2/3) * (F : \nabla v) - 1/3 * (Ic_iso / J) (CofF : \nabla v) ) }
864 void source_P1iso_NH ( Real coef,
865  const boost::multi_array<Real, 3 >& CofFk,
866  const boost::multi_array<Real, 3 >& Fk,
867  const std::vector<Real>& Jk,
868  const std::vector<Real>& Ic_isok ,
869  VectorElemental& elvec,
870  const CurrentFE& fe )
871 {
872  Real s1, s2;
873 
874  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
875  {
876  VectorElemental::vector_view vec = elvec.block ( icoor );
877  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
878  {
879  s1 = 0.0;
880  s2 = 0.0;
881  for ( UInt k = 0; k < nDimensions; ++k )
882  {
883  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
884  {
885  s1 += pow ( Jk[ ig ], (-2.0 / 3.0) ) * Fk[ icoor ][ k ][ ig ] *
886  fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
887 
888  s2 += 1.0 / 3.0 * ( Ic_isok[ ig ] * ( 1 / Jk[ig] ) ) *
889  CofFk[ icoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
890  }
891  }
892  vec ( i ) += (s1 - s2) * coef;
893  }
894  }
895 }
896 //! -----------------------------------------------------------------------------------------------
897 
898 
899 
900 
901 
902 
903 //! Jacobian matrix isochoric part ----------------------------------------------------------------
904 
905 //! 1. Jacobian matrix : Int { -2/3 * coef * J^(-5/3) *( CofF : \nabla \delta ) ( F : \nabla \v ) }
906 void stiff_Jac_P1iso_NH_1term ( Real coef,
907  const boost::multi_array<Real, 3 >& CofFk,
908  const boost::multi_array<Real, 3 >& Fk,
909  const std::vector<Real>& Jk ,
910  MatrixElemental& elmat,
911  const CurrentFE& fe )
912 {
913  Real s;
914 
915  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
916  {
917  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
918  {
919  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
920  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
921  {
922  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
923  {
924  s = 0.0;
925  for ( UInt l = 0; l < nDimensions; ++l )
926  {
927  for ( UInt k = 0; k < nDimensions; ++k )
928  {
929  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
930  {
931  s += pow ( Jk[ig], -5. / 3. ) *
932  Fk[ jcoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
933  CofFk[ icoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
934  }
935  }
936  }
937  mat ( i, j ) += coef * s;
938  }
939  }
940  }
941  }
942 }
943 
944 
945 
946 //! 2. Stiffness matrix: Int { 2/9 * coef * ( Ic_iso / J^2 )( CofF : \nabla \delta ) ( CofF : \nabla \v ) }
947 void stiff_Jac_P1iso_NH_2term ( Real coef,
948  const boost::multi_array<Real, 3 >& CofFk,
949  const std::vector<Real>& Jk ,
950  const std::vector<Real>& Ic_isok,
951  MatrixElemental& elmat,
952  const CurrentFE& fe )
953 {
954  Real s;
955 
956  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
957  {
958  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
959  {
960  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
961  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
962  {
963  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
964  {
965  s = 0.0;
966  for ( UInt l = 0; l < nDimensions; ++l )
967  {
968  for ( UInt k = 0; k < nDimensions; ++k )
969  {
970  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
971  {
972  s += ( 1 / (Jk[ig] * Jk[ig] ) ) * Ic_isok[ig] *
973  CofFk[ jcoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
974  CofFk[ icoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
975  }
976  }
977  }
978  mat ( i, j ) += coef * s;
979  }
980  }
981  }
982  }
983 }
984 
985 
986 
987 //! 3. Stiffness matrix : Int { coef * J^(-2/3) (\nabla \delta : \nabla \v)}
988 void stiff_Jac_P1iso_NH_3term ( Real coef,
989  const std::vector<Real>& Jk,
990  MatrixElemental& elmat,
991  const CurrentFE& fe )
992 {
993  Real s;
994 
995  //! assembling diagonal block
996  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
997 
998  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
999  {
1000  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1001  {
1002  s = 0.0;
1003  for ( UInt k = 0; k < nDimensions; ++k )
1004  {
1005  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1006  {
1007  s += pow ( Jk[ig], -2. / 3.) * fe.phiDer ( i, k, ig ) *
1008  fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1009  }
1010  }
1011  mat_tmp ( i, j ) = coef * s;
1012  }
1013  }
1014 
1015  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1016  {
1017  //! copy of diagonal block
1018  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
1019  mat += mat_tmp;
1020  }
1021 }
1022 
1023 //! 4. Stiffness matrix : Int { -2/3 * coef * J^(-5/3) ( F : \nabla \delta ) ( CofF : \nabla \v ) }
1024 void stiff_Jac_P1iso_NH_4term ( Real coef,
1025  const boost::multi_array<Real, 3 >& CofFk,
1026  const boost::multi_array<Real, 3 >& Fk,
1027  const std::vector<Real>& Jk ,
1028  MatrixElemental& elmat,
1029  const CurrentFE& fe )
1030 {
1031  Real s;
1032 
1033  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1034  {
1035  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1036  {
1037  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1038  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1039  {
1040  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1041  {
1042  s = 0.0;
1043  for ( UInt l = 0; l < nDimensions; ++l )
1044  {
1045  for ( UInt k = 0; k < nDimensions; ++k )
1046  {
1047  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1048  {
1049  s += pow ( Jk[ig], -5. / 3. ) *
1050  Fk[ icoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) *
1051  CofFk[ jcoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
1052  }
1053  }
1054  }
1055  mat ( i, j ) += coef * s;
1056  }
1057  }
1058  }
1059  }
1060 }
1061 
1062 
1063 
1064 //! 5. Stiffness matrix : Int { 1/3 * coef * J^(-2) * Ic_iso * (CofF [\nabla \delta]^t CofF ) : \nabla \v }
1065 void stiff_Jac_P1iso_NH_5term ( Real coef,
1066  const boost::multi_array<Real, 3 >& CofFk,
1067  const std::vector<Real>& Jk ,
1068  const std::vector<Real>& Ic_isok,
1069  MatrixElemental& elmat,
1070  const CurrentFE& fe )
1071 {
1072  Real s;
1073 
1074  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1075  {
1076  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1077  {
1078  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1079  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1080  {
1081  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1082  {
1083  s = 0.0;
1084  for ( UInt l = 0; l < nDimensions; ++l )
1085  {
1086  for ( UInt k = 0; k < nDimensions; ++k )
1087  {
1088  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1089  {
1090  s += ( 1 / ( Jk[ig] * Jk[ig] ) ) * Ic_isok[ig] *
1091  CofFk[ icoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
1092  CofFk[ jcoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
1093  }
1094  }
1095  }
1096  mat ( i, j ) += s * coef;
1097  }
1098  }
1099  }
1100  }
1101 }
1102 //! ***********************************************************************************************
1103 //! END OF NEO-HOOKEAN MODEL
1104 //! ***********************************************************************************************
1105 
1106 //! ***********************************************************************************************
1107 //! METHODS FOR EXPONENTIAL MODEL
1108 //! ***********************************************************************************************
1109 
1110 //! Stiffness vector isochoric part ---------------------------------------------------------------
1111 
1112 // Source term : Int { coef * exp(coefExp *( Ic_iso -3 )) * ( J^(-2/3)* (F : \nabla v) - 1/3 * (Ic_iso / J) * (CofF : \nabla v) ) }
1113 void source_P1iso_Exp ( Real coef,
1114  Real coefExp,
1115  const boost::multi_array<Real, 3 >& CofFk,
1116  const boost::multi_array<Real, 3 >& Fk,
1117  const std::vector<Real>& Jk,
1118  const std::vector<Real>& Ic_isok,
1119  VectorElemental& elvec,
1120  const CurrentFE& fe )
1121 {
1122 
1123  Real s;
1124 
1125  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1126  {
1127  VectorElemental::vector_view vec = elvec.block ( icoor );
1128  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1129  {
1130  s = 0.0;
1131  for ( UInt k = 0; k < nDimensions; ++k )
1132  {
1133  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1134  {
1135  s += exp ( coefExp * ( Ic_isok[ ig ] - 3.0 ) ) *
1136  (pow ( Jk[ ig ], (-2.0 / 3.0) ) * Fk[ icoor ][ k ][ ig ] -
1137  1.0 / 3.0 * ( 1 / Jk[ ig ] ) * Ic_isok[ ig ] *
1138  CofFk[ icoor ][ k ][ ig ] ) * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
1139 
1140  }
1141  }
1142  vec ( i ) += s * coef;
1143  }
1144  }
1145 }
1146 
1147 
1148 
1149 
1150 //! Jacobian matrix isochoric part ----------------------------------------------------------------
1151 
1152 //! 1. Stiffness term : Int { - 2/3 *coef * J^(-5/3) * exp( coefExp*( Ic_iso - 3) )* ( 1. + coefExp * Ic_iso ) * ( CofF : \nabla \delta ) ( F : \nabla \v ) }
1153 void stiff_Jac_P1iso_Exp_1term ( Real coef,
1154  Real coefExp,
1155  const boost::multi_array<Real, 3 >& CofFk,
1156  const boost::multi_array<Real, 3 >& Fk,
1157  const std::vector<Real>& Jk ,
1158  const std::vector<Real>& Ic_isok,
1159  MatrixElemental& elmat,
1160  const CurrentFE& fe )
1161 {
1162  Real s;
1163 
1164  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1165  {
1166  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1167  {
1168  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1169  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1170  {
1171  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1172  {
1173  s = 0.0;
1174  for ( UInt l = 0; l < nDimensions; ++l )
1175  {
1176  for ( UInt k = 0; k < nDimensions; ++k )
1177  {
1178  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1179  {
1180  s += std::pow ( Jk[ig], -5. / 3. ) * exp ( coefExp * ( Ic_isok[ig] - 3 ) ) *
1181  ( 1. + coefExp * Ic_isok[ig] ) *
1182  CofFk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) *
1183  Fk[ icoor ][ l ][ ig ] * fe.phiDer ( i, l, ig ) * fe.weightDet ( ig );
1184  }
1185  }
1186  }
1187  mat ( i, j ) += coef * s;
1188  }
1189  }
1190 
1191  }
1192  }
1193 }
1194 
1195 //! 2. Stiffness term : Int { 2 * coef * coefExp * J^(-4/3) * exp( coefExp*( Ic_iso - 3) ) * ( F : \nabla \delta ) ( F : \nabla \v )}
1196 void stiff_Jac_P1iso_Exp_2term ( Real coef,
1197  Real coefExp,
1198  const boost::multi_array<Real, 3 >& Fk,
1199  const std::vector<Real>& Jk,
1200  const std::vector<Real>& Ic_isok,
1201  MatrixElemental& elmat,
1202  const CurrentFE& fe )
1203 {
1204  Real s;
1205 
1206  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1207  {
1208  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1209  {
1210  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1211  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1212  {
1213  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1214  {
1215  s = 0.0;
1216  for ( UInt l = 0; l < nDimensions; ++l )
1217  {
1218  for ( UInt k = 0; k < nDimensions; ++k )
1219  {
1220  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1221  {
1222  s += std::pow ( Jk[ig], -4.0 / 3.0 ) * exp ( coefExp * ( Ic_isok[ig] - 3 ) ) *
1223  Fk[ jcoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
1224  Fk[ icoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
1225  }
1226  }
1227  }
1228  mat ( i, j ) += coef * s;
1229  }
1230  }
1231  }
1232  }
1233 }
1234 
1235 //! 3. Stiffness term: Int { 2.0/9.0 * coef * J^-2 * Ic_iso * exp( coefExp*( Ic_iso - 3) ) * ( 1. + coefExp * Ic_iso )( CofF : \nabla \delta ) ( CofF : \nabla \v )}
1236 void stiff_Jac_P1iso_Exp_3term ( Real coef, Real coefExp,
1237  const boost::multi_array<Real, 3 >& CofFk,
1238  const std::vector<Real>& Jk,
1239  const std::vector<Real>& Ic_isok,
1240  MatrixElemental& elmat,
1241  const CurrentFE& fe )
1242 {
1243  Real s;
1244 
1245  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1246  {
1247  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1248  {
1249  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1250  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1251  {
1252  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1253  {
1254  s = 0.0;
1255  for ( UInt l = 0; l < nDimensions; ++l )
1256  {
1257  for ( UInt k = 0; k < nDimensions; ++k )
1258  {
1259  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1260  {
1261  s += ( 1 / ( Jk[ig] * Jk[ig] ) ) * exp ( coefExp * ( Ic_isok[ig] - 3 ) ) *
1262  ( 1. + coefExp * Ic_isok[ig] ) * Ic_isok[ig] *
1263  CofFk[ jcoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
1264  CofFk[ icoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
1265 
1266  }
1267  }
1268  }
1269  mat ( i, j ) += coef * s;
1270  }
1271  }
1272  }
1273  }
1274 }
1275 
1276 //! 4. Stiffness term: Int { -2.0/3.0 * coef * J^(-5/3) * exp( coefExp*( Ic_iso - 3) ) * ( 1. + coefExp * Ic_iso )( F : \nabla \delta ) ( CofF : \nabla \v ) }
1277 void stiff_Jac_P1iso_Exp_4term ( Real coef, Real coefExp,
1278  const boost::multi_array<Real, 3 >& CofFk,
1279  const boost::multi_array<Real, 3 >& Fk,
1280  const std::vector<Real>& Jk ,
1281  const std::vector<Real>& Ic_isok,
1282  MatrixElemental& elmat,
1283  const CurrentFE& fe )
1284 {
1285  Real s;
1286 
1287  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1288  {
1289  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1290  {
1291  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1292  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1293  {
1294  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1295  {
1296  s = 0.0;
1297  for ( UInt l = 0; l < nDimensions; ++l )
1298  {
1299  for ( UInt k = 0; k < nDimensions; ++k )
1300  {
1301  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1302  {
1303  s += std::pow ( Jk[ig], -5. / 3. ) * exp ( coefExp * ( Ic_isok[ig] - 3 ) ) *
1304  ( 1. + coefExp * Ic_isok[ig] ) *
1305  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) *
1306  CofFk[ icoor ][ l ][ ig ] * fe.phiDer ( i, l, ig ) * fe.weightDet ( ig );
1307  }
1308  }
1309  }
1310  mat ( i, j ) += coef * s;
1311 
1312  }
1313  }
1314  }
1315  }
1316 }
1317 
1318 //! 5. Stiffness term : Int {coef * J^(-2/3) * exp( coefExp*( Ic_iso - 3)) (\nabla \delta: \nabla \v)}
1319 void stiff_Jac_P1iso_Exp_5term ( Real coef,
1320  Real coefExp,
1321  const std::vector<Real>& Jk,
1322  const std::vector<Real>& Ic_isok,
1323  MatrixElemental& elmat,
1324  const CurrentFE& fe )
1325 {
1326  Real s;
1327 
1328  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
1329  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1330  {
1331  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1332  {
1333  s = 0.0;
1334  for ( UInt k = 0; k < nDimensions; ++k )
1335  {
1336  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1337  {
1338  s += std::pow (Jk[ig], -2.0 / 3.0) * exp ( coefExp * ( Ic_isok[ig] - 3 ) ) *
1339  fe.phiDer ( i, k, ig ) * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1340  }
1341  }
1342  mat_tmp ( i, j ) = coef * s;
1343  }
1344  }
1345 
1346  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1347  {
1348  //! copy of diagonal block
1349  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
1350  mat += mat_tmp;
1351  }
1352 }
1353 
1354 //! 6. Stiffness term : Int { 1.0/3.0 * coef * J^(-2) * Ic_iso * exp(coefExp( Ic_iso - 3)) * (CofF [\nabla \delta]^t CofF ) : \nabla \v }
1355 void stiff_Jac_P1iso_Exp_6term ( Real coef,
1356  Real coefExp,
1357  const boost::multi_array<Real, 3 >& CofFk,
1358  const std::vector<Real>& Jk,
1359  const std::vector<Real>& Ic_isok,
1360  MatrixElemental& elmat,
1361  const CurrentFE& fe )
1362 {
1363  Real s;
1364 
1365  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1366  {
1367  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1368  {
1369  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1370  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1371  {
1372  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1373  {
1374  s = 0.0;
1375  for ( UInt l = 0; l < nDimensions; ++l )
1376  {
1377  for ( UInt k = 0; k < nDimensions; ++k )
1378  {
1379  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1380  {
1381  s += ( 1 / ( Jk[ig] * Jk[ig] ) ) * Ic_isok[ig] *
1382  exp ( coefExp * ( Ic_isok[ig] - 3 ) ) *
1383  CofFk[ icoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
1384  CofFk[ jcoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
1385  }
1386  }
1387  }
1388  mat ( i, j ) += s * coef;
1389  }
1390  }
1391  }
1392  }
1393 }
1394 
1395 //! ***********************************************************************************************
1396 //! END OF EXPONENTIAL MODEL
1397 //! ***********************************************************************************************
1398 
1399 //! ***********************************************************************************************
1400 //! METHODS FOR TENSORIAL CALCULUS
1401 //! ***********************************************************************************************
1402 
1403 //! Computation of the Right Cauchy Green tensor given the tensor F.
1404 void computeInvariantsRightCauchyGreenTensor (std::vector<LifeV::Real>& invariants,
1405  const Epetra_SerialDenseMatrix& tensorF,
1406  Epetra_SerialDenseMatrix& cofactorF)
1407 {
1408 
1409  //Computation of the invariants
1410  //At the moment, only the first one is really computed.
1411  //The others are not still used in the constitutive laws
1412 
1413  Real C11 (0);
1414  Real C22 (0);
1415  Real C33 (0);
1416 
1417  //It is not rescaled by the determinant. It is done inside the method to compute the local Piola
1418  //cofactorF.Scale(invariants[3]);
1419 
1420  C11 = tensorF (0, 0) * tensorF (0, 0) + tensorF (1, 0) * tensorF (1, 0) + tensorF (2, 0) * tensorF (2, 0);
1421  C22 = tensorF (0, 1) * tensorF (0, 1) + tensorF (1, 1) * tensorF (1, 1) + tensorF (2, 1) * tensorF (2, 1);
1422  C33 = tensorF (0, 2) * tensorF (0, 2) + tensorF (1, 2) * tensorF (1, 2) + tensorF (2, 2) * tensorF (2, 2);
1423 
1424  invariants[0] = C11 + C22 + C33; //First invariant C
1425  invariants[1] = 0.0; //Second invariant C
1426  invariants[2] = 0.0; //Third invariant C
1427  invariants[3] = tensorF (0, 0) * ( tensorF (1, 1) * tensorF (2, 2) - tensorF (1, 2) * tensorF (2, 1) ) - tensorF (0, 1) * ( tensorF (1, 0) * tensorF (2, 2) - tensorF (1, 2) * tensorF (2, 0) ) + tensorF (0, 2) * ( tensorF (1, 0) * tensorF (2, 1) - tensorF (1, 1) * tensorF (2, 0) ); //Determinant F
1428 
1429  //Computation of the Cofactor of F
1430  cofactorF ( 0 , 0 ) = ( tensorF (1, 1) * tensorF (2, 2) - tensorF (1, 2) * tensorF (2, 1) );
1431  cofactorF ( 0 , 1 ) = - ( tensorF (1, 0) * tensorF (2, 2) - tensorF (2, 0) * tensorF (1, 2) );
1432  cofactorF ( 0 , 2 ) = ( tensorF (1, 0) * tensorF (2, 1) - tensorF (1, 1) * tensorF (2, 0) );
1433  cofactorF ( 1 , 0 ) = - ( tensorF (0, 1) * tensorF (2, 2) - tensorF (0, 2) * tensorF (2, 1) );
1434  cofactorF ( 1 , 1 ) = ( tensorF (0, 0) * tensorF (2, 2) - tensorF (0, 2) * tensorF (2, 0) );
1435  cofactorF ( 1 , 2 ) = - ( tensorF (0, 0) * tensorF (2, 1) - tensorF (2, 0) * tensorF (0, 1) );
1436  cofactorF ( 2 , 0 ) = ( tensorF (0, 1) * tensorF (1, 2) - tensorF (0, 2) * tensorF (1, 1) );
1437  cofactorF ( 2 , 1 ) = - ( tensorF (0, 0) * tensorF (1, 2) - tensorF (0, 2) * tensorF (1, 0) );
1438  cofactorF ( 2 , 2 ) = ( tensorF (0, 0) * tensorF (1, 1) - tensorF (1, 0) * tensorF (0, 1) );
1439 
1440  cofactorF.Scale (1 / invariants[3]);
1441 }
1442 
1443 //! Computation of the Right Cauchy Green tensor given the tensor F.
1444 void computeInvariantsRightCauchyGreenTensor (std::vector<LifeV::Real>& invariants,
1445  const Epetra_SerialDenseMatrix& tensorF)
1446 {
1447 
1448  invariants[0] = 0.0; //First invariant C
1449  invariants[1] = 0.0; //Second invariant C
1450  invariants[2] = 0.0; //Third invariant C
1451  invariants[3] = tensorF (0, 0) * ( tensorF (1, 1) * tensorF (2, 2) - tensorF (1, 2) * tensorF (2, 1) ) - tensorF (0, 1) * ( tensorF (1, 0) * tensorF (2, 2) - tensorF (1, 2) * tensorF (2, 0) ) + tensorF (0, 2) * ( tensorF (1, 0) * tensorF (2, 1) - tensorF (1, 1) * tensorF (2, 0) ); //Determinant F
1452 
1453 }
1454 
1455 
1456 void computeCauchyStressTensor (Epetra_SerialDenseMatrix& cauchy,
1457  Epetra_SerialDenseMatrix& firstPiola,
1458  LifeV::Real det,
1459  Epetra_SerialDenseMatrix& tensorF)
1460 {
1461 
1462  firstPiola.Scale ( 1 / det );
1463  cauchy.Multiply ('N', 'T', 1.0, firstPiola, tensorF, 0.0);
1464 
1465 }
1466 
1467 void computeEigenvalues (const Epetra_SerialDenseMatrix& cauchy,
1468  std::vector<LifeV::Real>& eigenvaluesR,
1469  std::vector<LifeV::Real>& eigenvaluesI)
1470 
1471 {
1472 
1473  // LAPACK wrapper of Epetra
1474  Epetra_LAPACK lapack;
1475 
1476  //List of flags for Lapack Function
1477  //For documentation, have a look at http://www.netlib.org/lapack/double/dgeev.f
1478 
1479  char JOBVL = 'N';
1480  char JOBVR = 'N';
1481 
1482  //Size of the matrix
1483  Int Dim = cauchy.RowDim();
1484 
1485  //Arrays to store eigenvalues (their number = nDimensions)
1486  double WR[nDimensions];
1487  double WI[nDimensions];
1488 
1489  //Number of eigenvectors
1490  Int LDVR = nDimensions;
1491  Int LDVL = nDimensions;
1492 
1493  //Arrays to store eigenvectors
1494  Int length = nDimensions * 3;
1495 
1496  double VR[length];
1497  double VL[length];
1498 
1499  Int LWORK = 9;
1500  Int INFO = 0;
1501 
1502  double WORK[LWORK];
1503 
1504  double A[length];
1505 
1506  for (UInt i (0); i < nDimensions; i++)
1507  for (UInt j (0); j < nDimensions; j++)
1508  {
1509  A[nDimensions * i + j] = cauchy (i, j);
1510  }
1511 
1512  lapack.GEEV (JOBVL, JOBVR, Dim, A /*cauchy*/, Dim, &WR[0], &WI[0], VL, LDVL, VR, LDVR, WORK, LWORK, &INFO);
1513  ASSERT_PRE ( !INFO, "Calculation of the Eigenvalues failed!!!" );
1514 
1515  for ( UInt i (0); i < nDimensions; i++ )
1516  {
1517  eigenvaluesR[i] = WR[i];
1518  eigenvaluesI[i] = WI[i];
1519  }
1520 
1521 }
1522 
1523 //! ***********************************************************************************************
1524 //! METHODS FOR THE ST. VENANT KIRCHHOFF PENALIZED LAW
1525 //! ***********************************************************************************************
1526 //! Stiffness vector isochoric part ---------------------------------------------------------------
1527 
1528 // Source term : Int { ( \frac{lambda}{2} * Ic_iso - \frac{3}{2}*lambda - mu ) (F : \nabla v) - 1/3 * (Ic) * (F^-T : \nabla v) ) }
1529 void source_P1iso_VKPenalized ( Real lambda,
1530  Real mu,
1531  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1532  const boost::multi_array<Real, 3 >& Fk,
1533  const std::vector<Real>& Ic_isok,
1534  const std::vector<Real>& Ic_k,
1535  const std::vector<Real>& Jack_k,
1536  VectorElemental& elvec,
1537  const CurrentFE& fe )
1538 {
1539 
1540  Real s;
1541 
1542  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1543  {
1544  VectorElemental::vector_view vec = elvec.block ( icoor );
1545  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1546  {
1547  s = 0.0;
1548  for ( UInt k = 0; k < nDimensions; ++k )
1549  {
1550  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1551  {
1552  s += std::pow ( Jack_k[ ig ], (-2.0 / 3.0) ) * ( ( lambda / 2.0 ) * Ic_isok[ ig ] - (3.0 / 2.0) * lambda - mu ) *
1553  (Fk[ icoor ][ k ][ ig ] - (1.0 / 3.0) * Ic_k[ ig ] * FkMinusTransposed[ icoor ][ k ][ ig ] ) * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
1554 
1555  }
1556  }
1557  vec ( i ) += s;
1558  }
1559  }
1560 }
1561 
1562 
1563 // Source term : Int { ( 2 * mu * Jk^(-4.0/3.0) ) * ( (F*C : \nabla v) - 1/3 * (Ic_Squared) * (F^-T : \nabla v) ) }
1564 void source_P2iso_VKPenalized ( Real mu,
1565  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1566  const boost::multi_array<Real, 3 >& FkCk,
1567  const std::vector<Real>& Ic_Squared,
1568  const std::vector<Real>& Jk,
1569  VectorElemental& elvec,
1570  const CurrentFE& fe )
1571 {
1572 
1573  Real s;
1574 
1575  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1576  {
1577  VectorElemental::vector_view vec = elvec.block ( icoor );
1578  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1579  {
1580  s = 0.0;
1581  for ( UInt k = 0; k < nDimensions; ++k )
1582  {
1583  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1584  {
1585  s += ( mu * std::pow ( Jk[ ig ], (-4.0 / 3.0) ) ) * (FkCk[ icoor ][ k ][ ig ] - (1.0 / 3.0) * Ic_Squared[ ig ] * FkMinusTransposed[ icoor ][ k ][ ig ] ) * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
1586  }
1587  }
1588  vec ( i ) += s;
1589  }
1590  }
1591 }
1592 
1593 //! Jacobian matrix of the first Piola-Kirchhoff tensor for the VK Penalized law
1594 //! 1. Stiffness term : int { -(2.0/3.0) * Jk^(-2.0/3.0) * ( (lambda/2) * Ic_isok - ( (3/2)*lambda + mu ) ) * F^-T:\nabla \delta ) * ( F - (1.0/3.0) * Ic_k * F^-T ): \nabla \v }
1595 void stiff_Jac_P1iso_VKPenalized_0term ( Real lambda, Real mu,
1596  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1597  const boost::multi_array<Real, 3 >& Fk,
1598  const std::vector<Real>& Jk ,
1599  const std::vector<Real>& Ic_k ,
1600  const std::vector<Real>& IcIso_k ,
1601  MatrixElemental& elmat,
1602  const CurrentFE& fe )
1603 {
1604  Real s;
1605 
1606  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1607  {
1608  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1609  {
1610  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1611  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1612  {
1613  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1614  {
1615  s = 0.0;
1616  for ( UInt l = 0; l < nDimensions; ++l )
1617  {
1618  for ( UInt k = 0; k < nDimensions; ++k )
1619  {
1620  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1621  {
1622  s += (-2.0 / 3.0) * std::pow (Jk[ ig ], (-2.0 / 3.0) ) * ( ( (lambda / 2.0) * IcIso_k[ ig ] ) - ( (3.0 / 2.0) * lambda + mu) ) *
1623  fe.phiDer ( i, l, ig ) * ( Fk[ icoor ][ l ][ ig ] - (1.0 / 3.0) * Ic_k[ ig ] * FkMinusTransposed[ icoor ][ l ][ ig ] ) *
1624  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1625  }
1626  }
1627  }
1628  mat ( i, j ) += s;
1629  }
1630  }
1631  }
1632  }
1633 }
1634 
1635 //! 1. Stiffness term :int{ J^(-2/3) * (lambda / 2) * ( (-2/3) * Ic_k * J^(-2/3) * F^-T:\nabla \delta ) * ( F : \nabla \v )}
1636 void stiff_Jac_P1iso_VKPenalized_1term ( Real coeff,
1637  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1638  const boost::multi_array<Real, 3 >& Fk,
1639  const std::vector<Real>& Jk ,
1640  const std::vector<Real>& Ic_k ,
1641  MatrixElemental& elmat,
1642  const CurrentFE& fe )
1643 {
1644  Real s;
1645 
1646  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1647  {
1648  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1649  {
1650  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1651  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1652  {
1653  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1654  {
1655  s = 0.0;
1656  for ( UInt l = 0; l < nDimensions; ++l )
1657  {
1658  for ( UInt k = 0; k < nDimensions; ++k )
1659  {
1660  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1661  {
1662  s += ( -2.0 / 3.0 ) * Ic_k[ ig ] * std::pow ( Jk[ig], (-4.0 / 3.0) ) *
1663  fe.phiDer ( i, l, ig ) * Fk[ icoor ][ l ][ ig ] *
1664  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1665  }
1666  }
1667  }
1668  mat ( i, j ) += coeff * s;
1669  }
1670  }
1671  }
1672  }
1673 }
1674 
1675 //! 3. Stiffness term:int { J^(-2/3) * (lambda / 2) * ( ( 2/9 ) * J^(-2/3) * Ic_k^2 ) * ( F^-T : \nabla \delta ) ( F^-T : \nabla \v ) }
1676 void stiff_Jac_P1iso_VKPenalized_2term ( Real coef,
1677  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1678  const std::vector<Real>& Jk,
1679  const std::vector<Real>& Ic_k,
1680  MatrixElemental& elmat,
1681  const CurrentFE& fe )
1682 {
1683  Real s;
1684 
1685  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1686  {
1687  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1688  {
1689  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1690  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1691  {
1692  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1693  {
1694  s = 0.0;
1695  for ( UInt l = 0; l < nDimensions; ++l )
1696  {
1697  for ( UInt k = 0; k < nDimensions; ++k )
1698  {
1699  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1700  {
1701  s += ( ( 2.0 / 9.0 ) * std::pow ( Jk[ig], (-4.0 / 3.0) ) * Ic_k[ig] * Ic_k[ig] ) *
1702  FkMinusTransposed[ icoor ][ l ][ ig ] * fe.phiDer ( i, l, ig ) *
1703  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1704  }
1705  }
1706  }
1707  mat ( i, j ) += coef * s;
1708  }
1709  }
1710  }
1711  }
1712 }
1713 
1714 //! 2. Stiffness term : Int { J^(-2/3) * coeff * ( 2 * J^(-2/3) ) * ( F : \nabla \delta ) ( F : \nabla \v )}
1715 void stiff_Jac_P1iso_VKPenalized_3term ( Real coef,
1716  const boost::multi_array<Real, 3 >& Fk,
1717  const std::vector<Real>& Jk,
1718  MatrixElemental& elmat,
1719  const CurrentFE& fe )
1720 {
1721  Real s;
1722 
1723  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1724  {
1725  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1726  {
1727  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1728  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1729  {
1730  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1731  {
1732  s = 0.0;
1733  for ( UInt l = 0; l < nDimensions; ++l )
1734  {
1735  for ( UInt k = 0; k < nDimensions; ++k )
1736  {
1737  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1738  {
1739  s += 2.0 * std::pow ( Jk[ig], (-4.0 / 3.0) ) *
1740  fe.phiDer ( i, l, ig ) * Fk[ icoor ][ l ][ ig ] *
1741  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1742  }
1743  }
1744  }
1745  mat ( i, j ) += coef * s;
1746  }
1747  }
1748  }
1749  }
1750 }
1751 
1752 //! 4. Stiffness term:int{J^(-2/3) * (lambda/2) * ( -2.0/3.0 * J^(-2/3) * Ic_k ) * ( F : \nabla \delta)*(F^-T : \nabla \v )}
1753 void stiff_Jac_P1iso_VKPenalized_4term ( Real coef,
1754  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1755  const boost::multi_array<Real, 3 >& Fk,
1756  const std::vector<Real>& Jk ,
1757  const std::vector<Real>& Ic_k,
1758  MatrixElemental& elmat,
1759  const CurrentFE& fe )
1760 {
1761  Real s;
1762 
1763  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1764  {
1765  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1766  {
1767  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1768  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1769  {
1770  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1771  {
1772  s = 0.0;
1773  for ( UInt l = 0; l < nDimensions; ++l )
1774  {
1775  for ( UInt k = 0; k < nDimensions; ++k )
1776  {
1777  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1778  {
1779  s += ( ( -2.0 / 3.0 ) * std::pow ( Jk[ig], (-4.0 / 3.0) ) * Ic_k[ig] ) *
1780  fe.phiDer ( i, l, ig ) * FkMinusTransposed[ icoor ][ l ][ ig ] *
1781  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1782  }
1783  }
1784  }
1785  mat ( i, j ) += coef * s;
1786  }
1787  }
1788  }
1789  }
1790 }
1791 
1792 //! 5. Stiffness term : int { J^(-2.0/3.0) * ( (lambda/2) * Ic_isok - ( (3/2)*lambda + mu ) ) * \nabla \delta : \nabla v }
1793 void stiff_Jac_P1iso_VKPenalized_5term ( Real coef,
1794  Real secondCoef,
1795  const std::vector<Real>& Jk,
1796  const std::vector<Real>& Ic_isok,
1797  MatrixElemental& elmat,
1798  const CurrentFE& fe )
1799 {
1800  Real s;
1801 
1802  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
1803  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1804  {
1805  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1806  {
1807  s = 0.0;
1808  for ( UInt k = 0; k < nDimensions; ++k )
1809  {
1810  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1811  {
1812  s += std::pow ( Jk[ ig ], (-2.0 / 3.0) ) * ( (coef / 2.0) * Ic_isok[ ig ] - ( (3.0 / 2.0) * coef + secondCoef ) ) *
1813  fe.phiDer ( i, k, ig ) * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1814  }
1815  }
1816  mat_tmp ( i, j ) = s;
1817  }
1818  }
1819 
1820  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1821  {
1822  //! copy of diagonal block
1823  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
1824  mat += mat_tmp;
1825  }
1826 }
1827 
1828 //! 6. Stiffness term: int { J^(-2.0/3.0) * ( (lambda/2) * Ic_isok - ( (3/2)*lambda + mu ) ) * ( (-2/3) * ( F :\nabla \delta ) ) * ( F^-T : \nabla v ) }
1829 void stiff_Jac_P1iso_VKPenalized_6term ( Real coef, Real secondCoef,
1830  const std::vector<Real>& Jk,
1831  const std::vector<Real>& Ic_isok,
1832  const boost::multi_array<Real, 3 >& Fk,
1833  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1834  MatrixElemental& elmat,
1835  const CurrentFE& fe )
1836 {
1837  Real s;
1838 
1839  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1840  {
1841  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1842  {
1843  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1844  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1845  {
1846  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1847  {
1848  s = 0.0;
1849  for ( UInt l = 0; l < nDimensions; ++l )
1850  {
1851  for ( UInt k = 0; k < nDimensions; ++k )
1852  {
1853  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1854  {
1855  s += std::pow ( Jk[ ig ], (-2.0 / 3.0) ) * (-2.0 / 3.0) * ( (coef / 2.0) * Ic_isok[ ig ] - ( (3.0 / 2.0) * coef + secondCoef ) ) *
1856  fe.phiDer ( i, l, ig ) * FkMinusTransposed[ icoor ][ l ][ ig ] *
1857  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1858  }
1859  }
1860  }
1861  mat ( i, j ) += s;
1862  }
1863  }
1864  }
1865  }
1866 }
1867 
1868 //! 7. Stiffness term : int { ( J^(-2/3) * (lambda/2) * Ic_isok - ( (3/2)*lambda + mu ) ) * ( (1/3) * Ic_k * ( F^-T \nabla \delta^T F-T ) : \nabla v }
1869 void stiff_Jac_P1iso_VKPenalized_7term ( Real coef,
1870  Real secondCoef,
1871  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1872  const std::vector<Real>& Ic_isok,
1873  const std::vector<Real>& Ic_k,
1874  const std::vector<Real>& Jk,
1875  MatrixElemental& elmat,
1876  const CurrentFE& fe )
1877 {
1878  Real s;
1879 
1880  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1881  {
1882  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1883  {
1884  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1885  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1886  {
1887  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1888  {
1889  s = 0.0;
1890  for ( UInt l = 0; l < nDimensions; ++l )
1891  {
1892  for ( UInt k = 0; k < nDimensions; ++k )
1893  {
1894  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1895  {
1896  s += std::pow ( Jk[ ig ], (-2.0 / 3.0) ) * ( (coef / 2.0) * Ic_isok[ ig ] - ( (3.0 / 2.0) * coef + secondCoef ) ) * ( (1.0 / 3.0) * Ic_k[ ig ] ) *
1897  FkMinusTransposed[ icoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
1898  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
1899  }
1900  }
1901  }
1902  mat ( i, j ) += s;
1903  }
1904  }
1905  }
1906  }
1907 }
1908 
1909 //! 8. Stiffness term : int { ( -4.0/3.0) * ( mu * J^(-4/3) ) * ( F^-T: \grad \delta ) * ( F C ) : \nabla v }
1910 void stiff_Jac_P1iso_VKPenalized_8term ( Real coef, const std::vector<Real>& Jack_k,
1911  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1912  const boost::multi_array<Real, 3 >& FkCk,
1913  MatrixElemental& elmat,
1914  const CurrentFE& fe )
1915 {
1916 
1917  Real s;
1918 
1919  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1920  {
1921  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1922  {
1923  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1924  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1925  {
1926  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1927  {
1928  s = 0.0;
1929  for ( UInt l = 0; l < nDimensions; ++l )
1930  {
1931  for ( UInt k = 0; k < nDimensions; ++k )
1932  {
1933  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1934  {
1935  s += (-4.0 / 3.0) * std::pow ( Jack_k[ ig ], (-4.0 / 3.0) ) *
1936  fe.phiDer ( i, l, ig ) * FkCk[ icoor ][ l ][ ig ] *
1937  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1938  }
1939  }
1940  }
1941  mat ( i, j ) += s * coef;
1942  }
1943  }
1944  }
1945  }
1946 }
1947 
1948 
1949 //! 8. Stiffness term : int { ( 4.0/9.0) * ( mu * J^(-4/3) ) Ic_kSquared * (F^-T : \grad \delta ) * F^-T : \nabla \v }
1950 void stiff_Jac_P1iso_VKPenalized_9term ( Real coef, const std::vector<Real>& Jack_k,
1951  const std::vector<Real>& Ic_kSquared,
1952  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1953  MatrixElemental& elmat,
1954  const CurrentFE& fe )
1955 {
1956  Real s;
1957  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1958  {
1959  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1960  {
1961  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1962  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1963  {
1964  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1965  {
1966  s = 0.0;
1967  for ( UInt l = 0; l < nDimensions; ++l )
1968  {
1969  for ( UInt k = 0; k < nDimensions; ++k )
1970  {
1971  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1972  {
1973  s += ( (-4.0 / 3.0) * std::pow ( Jack_k[ ig ], (-4.0 / 3.0) ) ) *
1974  fe.phiDer ( i, l, ig ) * (-1.0 / 3.0) * Ic_kSquared[ ig ] * FkMinusTransposed[ icoor ][ l ][ ig ] *
1975  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1976  }
1977  }
1978  }
1979  mat ( i, j ) += s * coef;
1980  }
1981  }
1982  }
1983  }
1984 }
1985 
1986 
1987 
1988 
1989 //! 10. Stiffness term : int { ( mu * J^(-4/3) ) * ( \nabla \delta * C ) : \nabla v }
1990 void stiff_Jac_P1iso_VKPenalized_10term ( Real coef,
1991  const std::vector<Real>& Jack_k,
1992  const boost::multi_array<Real, 3 >& Ck,
1993  MatrixElemental& elmat,
1994  const CurrentFE& fe )
1995 {
1996  Real s;
1997 
1998  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1999  {
2000  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
2001  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2002  {
2003  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2004  {
2005  s = 0.0;
2006  for ( UInt l = 0; l < nDimensions; ++l )
2007  {
2008  for ( UInt k = 0; k < nDimensions; ++k )
2009  {
2010  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2011  {
2012  s += std::pow ( Jack_k[ ig ], (-4.0 / 3.0) ) * Ck[ l ][ k ][ ig ] * fe.phiDer ( i, k, ig ) *
2013  fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
2014  }
2015  }
2016  }
2017  mat ( i, j ) += coef * s;
2018  }
2019  }
2020  }
2021 }
2022 
2023 //! 6. Stiffness term : int { ( mu * J^(-4/3) ) * (F [\nabla \delta]^T F ) : \nabla \v }
2024 void stiff_Jac_P1iso_VKPenalized_11term ( Real coef,
2025  const std::vector<Real>& Jk,
2026  const boost::multi_array<Real, 3 >& Fk,
2027  MatrixElemental& elmat,
2028  const CurrentFE& fe )
2029 {
2030  Real s;
2031 
2032  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2033  {
2034  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2035  {
2036  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2037  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2038  {
2039  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2040  {
2041  s = 0.0;
2042  for ( UInt l = 0; l < nDimensions; ++l )
2043  {
2044  for ( UInt k = 0; k < nDimensions; ++k )
2045  {
2046  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2047  {
2048  s += std::pow ( Jk[ ig ], (-4.0 / 3.0) ) * fe.phiDer ( i, l, ig ) * Fk[ l ][ jcoor ][ ig ] *
2049  Fk[ icoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2050  }
2051  }
2052  }
2053  mat ( i, j ) += coef * s;
2054  }
2055  }
2056  }
2057  }
2058 }
2059 
2060 //! 6. Stiffness term : int { ( mu * J^(-4/3) ) * (F * F^T * [\nabla \delta] ) : \nabla \v }
2061 void stiff_Jac_P1iso_VKPenalized_12term ( Real coef,
2062  const std::vector<Real>& Jk,
2063  const boost::multi_array<Real, 3 >& Fk,
2064  MatrixElemental& elmat,
2065  const CurrentFE& fe )
2066 {
2067  Real s;
2068 
2069  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2070  {
2071  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2072  {
2073  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2074  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2075  {
2076  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2077  {
2078  s = 0.0;
2079  for ( UInt p = 0; p < nDimensions; p++ )
2080  {
2081  for ( UInt k = 0; k < nDimensions; k++ )
2082  {
2083  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2084  {
2085  s += std::pow ( Jk[ig], (-4.0 / 3.0) ) * Fk[ icoor ][ p ][ ig ] * fe.phiDer ( i, k, ig ) *
2086  Fk[ jcoor ][ p ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2087  }
2088  }
2089  }
2090  mat ( i, j ) += coef * s;
2091  }
2092  }
2093  }
2094  }
2095 }
2096 
2097 //! 11. Stiffness term : int { ( mu * J^(-4/3) * ( (1/3) * Ic_SquaredK * ( F^-T [\nabla \delta ]^T F^-T) : \nabla v ) }
2098 void stiff_Jac_P1iso_VKPenalized_13term ( Real coef,
2099  const std::vector<Real>& Jk,
2100  const std::vector<Real>& Ic_kSquared,
2101  const boost::multi_array<Real, 3 >& FkMinusTransposed,
2102  MatrixElemental& elmat,
2103  const CurrentFE& fe )
2104 {
2105  Real s;
2106 
2107  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2108  {
2109  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2110  {
2111  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2112  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2113  {
2114  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2115  {
2116  s = 0.0;
2117  for ( UInt l = 0; l < nDimensions; ++l )
2118  {
2119  for ( UInt k = 0; k < nDimensions; ++k )
2120  {
2121  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2122  {
2123  s += ( std::pow ( Jk[ig], (-4.0 / 3.0) ) ) * ( Ic_kSquared[ig] / 3.0 ) *
2124  fe.phiDer ( j, l, ig ) * FkMinusTransposed[ icoor ][ l ][ ig ] *
2125  fe.phiDer ( i, k, ig ) * FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.weightDet ( ig );
2126 
2127  }
2128  }
2129  }
2130  mat ( i, j ) += s * coef;
2131  }
2132  }
2133  }
2134  }
2135 }
2136 
2137 //! 12. Stiffness term : int { ( mu * J^(-4/3) ) * ( (-4.0/3.0) * ( FkCk : \nabla \delta ) ) * F^-T : \nabla v ) }
2138 void stiff_Jac_P1iso_VKPenalized_14term ( Real coef,
2139  const std::vector<Real>& Jk,
2140  const boost::multi_array<Real, 3 >& FkCk,
2141  const boost::multi_array<Real, 3 >& FkMinusTransposed,
2142  MatrixElemental& elmat,
2143  const CurrentFE& fe )
2144 {
2145  Real s;
2146 
2147  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2148  {
2149  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2150  {
2151  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2152  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2153  {
2154  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2155  {
2156  s = 0.0;
2157  for ( UInt l = 0; l < nDimensions; ++l )
2158  {
2159  for ( UInt k = 0; k < nDimensions; ++k )
2160  {
2161  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2162  {
2163  s += ( std::pow ( Jk[ig], (-4.0 / 3.0) ) ) * ( -4.0 / 3.0 ) *
2164  fe.phiDer ( i, l, ig ) * FkMinusTransposed[ icoor ][ l ][ ig ] *
2165  ( FkCk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) ) * fe.weightDet ( ig );
2166  }
2167  }
2168  }
2169  mat ( i, j ) += s * coef;
2170  }
2171  }
2172  }
2173  }
2174 }
2175 
2176 //! ***********************************************************************************************
2177 //! END OF VENANT-KIRCHHOFF PENALIZED MODEL
2178 //! ***********************************************************************************************
2179 
2180 //! ***********************************************************************************************
2181 //! SECOND ORDER EXPONENTIAL MODEL
2182 //! ***********************************************************************************************
2183 //! Methods for the Second Order Exponential law
2184 //! Assemble the first Piola-Kirchhoff tensor
2185 //int { 2 * alpha * ( Ic1_iso - 3 ) * exp(gamma *( Ic1_iso -3 )^2) * ( J1^(-2/3)* (F1 : \nabla v) - 1/3 * (Ic1_iso / J1) * (CofF1 : \nabla v) ) }
2186 void source_P1iso_SecondOrderExponential ( Real coef, Real coefExp,
2187  const boost::multi_array<Real, 3 >& CofFk,
2188  const boost::multi_array<Real, 3 >& Fk,
2189  const std::vector<Real>& Jk,
2190  const std::vector<Real>& Ic_isok,
2191  VectorElemental& elvec,
2192  const CurrentFE& fe )
2193 {
2194 
2195  Real s;
2196 
2197  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2198  {
2199  VectorElemental::vector_view vec = elvec.block ( icoor );
2200  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2201  {
2202  s = 0.0;
2203  for ( UInt k = 0; k < nDimensions; ++k )
2204  {
2205  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2206  {
2207  s += (Ic_isok[ ig ] - 3.0) * exp ( coefExp * (Ic_isok[ ig ] - 3.0) * (Ic_isok[ ig ] - 3.0) ) *
2208  (std::pow ( Jk[ ig ], (-2.0 / 3.0) ) * Fk[ icoor ][ k ][ ig ]
2209  - 1.0 / 3.0 * ( Ic_isok[ ig ] / Jk[ ig ] ) * CofFk[ icoor ][ k ][ ig ] ) * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
2210  }
2211  }
2212  vec ( i ) += s * coef;
2213  }
2214  }
2215 }
2216 
2217 //! 1. Stiffness term : Int { ( - 4/3 *coef * J^(-5/3) * exp( coefExp*( Ic_iso - 3)^2 ) * ( tr_isoCk + ( tr_isoCk - 3 ) * ( 2 *coefExp * ( tr_isoCk - 3 ) + 1 ) ) ) * ( CofF : \nabla \delta ) ( F : \nabla \v ) }
2218 // In this method, the coef that is passed is -4/3 * alpha, therefore it can multiply the sum at the end.
2219 void stiff_Jac_P1iso_SecondOrderExp_1term ( Real coef,
2220  Real coefExp,
2221  const boost::multi_array<Real, 3 >& CofFk,
2222  const boost::multi_array<Real, 3 >& Fk,
2223  const std::vector<Real>& Jk ,
2224  const std::vector<Real>& Ic_isok,
2225  MatrixElemental& elmat,
2226  const CurrentFE& fe )
2227 {
2228  Real s;
2229 
2230  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2231  {
2232  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2233  {
2234  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2235  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2236  {
2237  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2238  {
2239  s = 0.0;
2240  for ( UInt l = 0; l < nDimensions; ++l )
2241  {
2242  for ( UInt k = 0; k < nDimensions; ++k )
2243  {
2244  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2245  {
2246  s += ( std::pow ( Jk[ig], -5.0 / 3.0) * exp (coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) ) *
2247  ( (Ic_isok[ig] - 3.0) * ( 1.0 + 2 * coefExp * (Ic_isok[ig] - 3.0) * Ic_isok[ig] ) + Ic_isok[ig] ) *
2248  fe.phiDer ( i, l, ig ) * Fk[ icoor ][ l ][ ig ] *
2249  CofFk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2250  }
2251  }
2252  }
2253  mat ( i, j ) += coef * s;
2254  }
2255  }
2256  }
2257  }
2258 }
2259 
2260 //! 2. Stiffness term : Int { 4 * coef * J^(-4/3) * exp( coefExp*( Ic_iso - 3)*( Ic_iso - 3) ) *
2261 //! ( 1.0 + 2 * coefExp * (Ic_isoK - 3)^2 ) * ( F : \nabla \delta ) ( F : \nabla \v )}
2262 // When the method is called, the coef parameter stores already 4 * alpha. This is why it is used at the end.
2263 void stiff_Jac_P1iso_SecondOrderExp_2term ( Real coef,
2264  Real coefExp,
2265  const boost::multi_array<Real, 3 >& Fk,
2266  const std::vector<Real>& Jk,
2267  const std::vector<Real>& Ic_isok,
2268  MatrixElemental& elmat,
2269  const CurrentFE& fe )
2270 {
2271  Real s;
2272 
2273  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2274  {
2275  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2276  {
2277  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2278  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2279  {
2280  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2281  {
2282  s = 0.0;
2283  for ( UInt l = 0; l < nDimensions; ++l )
2284  {
2285  for ( UInt k = 0; k < nDimensions; ++k )
2286  {
2287  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2288  {
2289  s += std::pow ( Jk[ig], -4.0 / 3.0 ) * exp ( coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) *
2290  ( 1.0 + 2 * coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) *
2291  fe.phiDer ( i, l, ig ) * Fk[ icoor ][ l ][ ig ] *
2292  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2293  }
2294  }
2295  }
2296  mat ( i, j ) += coef * s;
2297  }
2298  }
2299  }
2300  }
2301 }
2302 
2303 //! 3. Stiffness term:
2304 //Int { ( 4.0/9.0 * alpha * J^-2 * exp( gamma*( Ic_iso - 3)^2 ) * ( Ic_isoK + 2 * gamma * (Ic_isoK - 3)^2 * Ic_isoK + (Ic_isoK - 3) ) ) *
2305 // ( CofF : \nabla \delta ) ( CofF : \nabla \v )}
2306 // The coef that is passed stores 4/9 * alpha.
2307 void stiff_Jac_P1iso_SecondOrderExp_3term ( Real coef, Real coefExp,
2308  const boost::multi_array<Real, 3 >& CofFk,
2309  const std::vector<Real>& Jk,
2310  const std::vector<Real>& Ic_isok,
2311  MatrixElemental& elmat,
2312  const CurrentFE& fe )
2313 {
2314  Real s;
2315 
2316  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2317  {
2318  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2319  {
2320  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2321  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2322  {
2323  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2324  {
2325  s = 0.0;
2326  for ( UInt l = 0; l < nDimensions; ++l )
2327  {
2328  for ( UInt k = 0; k < nDimensions; ++k )
2329  {
2330  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2331  {
2332  s += ( ( 1 / (Jk[ig] * Jk[ig]) ) * Ic_isok[ig] * exp ( coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) ) *
2333  ( (Ic_isok[ig] - 3.0) * ( 1.0 + 2.0 * coefExp * (Ic_isok[ig] - 3.0) * Ic_isok[ig] ) + Ic_isok[ig] ) *
2334  CofFk[ icoor ][ l ][ ig ] * fe.phiDer ( i, l, ig ) *
2335  CofFk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2336  }
2337  }
2338  }
2339  mat ( i, j ) += coef * s;
2340  }
2341  }
2342  }
2343  }
2344 }
2345 
2346 
2347 //! 4. Stiffness term:
2348 //Int { (-4.0/3.0 * alpha * J^(-5/3) * exp( gamma*( Ic_iso - 3)*( Ic_iso - 3) ) * ( Ic_isoK + ( Ic_isok - 3.0 ) * (2*gamma*(Ic_isok - 3)Ic + 1) ) *
2349 // ( F : \nabla \delta ) ( CofF : \nabla \v ) }
2350 // As in other methods, the coef that is passed is -4.0/3.0 * alpha.
2351 void stiff_Jac_P1iso_SecondOrderExp_4term ( Real coef, Real coefExp,
2352  const boost::multi_array<Real, 3 >& CofFk,
2353  const boost::multi_array<Real, 3 >& Fk,
2354  const std::vector<Real>& Jk ,
2355  const std::vector<Real>& Ic_isok,
2356  const std::vector<Real>& /*Ic_k*/,
2357  MatrixElemental& elmat,
2358  const CurrentFE& fe )
2359 {
2360  Real s;
2361 
2362  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2363  {
2364  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2365  {
2366  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2367  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2368  {
2369  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2370  {
2371  s = 0.0;
2372  for ( UInt l = 0; l < nDimensions; ++l )
2373  {
2374  for ( UInt k = 0; k < nDimensions; ++k )
2375  {
2376  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2377  {
2378  s += std::pow ( Jk[ig], (-5.0 / 3.0) ) * exp ( coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) *
2379  ( Ic_isok[ig] + (Ic_isok[ig] - 3.0) * ( 1.0 + 2.0 * coefExp * (Ic_isok[ig] - 3.0) * Ic_isok[ig] ) ) *
2380  fe.phiDer ( i, l, ig ) * CofFk[ icoor ][ l ][ ig ] *
2381  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2382  }
2383  }
2384  }
2385  mat ( i, j ) += coef * s;
2386  }
2387  }
2388  }
2389  }
2390 }
2391 
2392 //! 5. Stiffness term :
2393 //Int {coef * J^(-2/3) * exp( coefExp*( Ic_iso - 3)*( Ic_iso - 3)) * ( Ic_iso - 3)* (\nabla \delta: \nabla \v)}
2394 void stiff_Jac_P1iso_SecondOrderExp_5term ( Real coef,
2395  Real coefExp,
2396  const std::vector<Real>& Jk,
2397  const std::vector<Real>& Ic_isok,
2398  MatrixElemental& elmat,
2399  const CurrentFE& fe )
2400 {
2401  Real s;
2402 
2403  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
2404  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2405  {
2406  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2407  {
2408  s = 0.0;
2409  for ( UInt k = 0; k < nDimensions; ++k )
2410  {
2411  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2412  {
2413  s += (Ic_isok[ ig ] - 3.0) * exp ( coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) *
2414  std::pow ( Jk[ig], -2.0 / 3.0 ) * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2415  }
2416  }
2417  mat_tmp ( i, j ) = coef * s;
2418  }
2419  }
2420 
2421  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2422  {
2423  //! copy of diagonal block
2424  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
2425  mat += mat_tmp;
2426  }
2427 }
2428 
2429 //! 6. Stiffness term : Int { 2.0/3.0 * coef * J^(-2) * Ic_iso * exp(coefExp( Ic_iso - 3)) * (CofF [\nabla \delta]^t CofF ) : \nabla \v }
2430 void stiff_Jac_P1iso_SecondOrderExp_6term ( Real coef,
2431  Real coefExp,
2432  const boost::multi_array<Real, 3 >& CofFk,
2433  const std::vector<Real>& Jk,
2434  const std::vector<Real>& Ic_isok,
2435  MatrixElemental& elmat,
2436  const CurrentFE& fe )
2437 {
2438  Real s;
2439 
2440  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2441  {
2442  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2443  {
2444  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2445  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2446  {
2447  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2448  {
2449  s = 0.0;
2450  for ( UInt l = 0; l < nDimensions; ++l )
2451  {
2452  for ( UInt k = 0; k < nDimensions; ++k )
2453  {
2454  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2455  {
2456  s += ( 1 / (Jk[ig] * Jk[ig]) ) * (Ic_isok[ig] - 3.0) *
2457  exp ( coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) * Ic_isok[ig] *
2458  fe.phiDer ( j, l, ig ) * CofFk[ icoor ][ l ][ ig ] *
2459  fe.phiDer ( i, k, ig ) * CofFk[ jcoor ][ k ][ ig ] * fe.weightDet ( ig );
2460  }
2461  }
2462  }
2463  mat ( i, j ) += s * coef;
2464  }
2465  }
2466  }
2467  }
2468 }
2469 
2470 //! ***********************************************************************************************
2471 //! END OF SECOND ORDER EXPONENTIAL MODEL
2472 //! ***********************************************************************************************
2473 
2474 } //! End namespace AssemblyElementalStructure
2475 
2476 } //! End namespace LifeV
2477 
2478 
2479 #endif
LifeV::Int
int32_type Int
Generic integer data.
Definition: LifeV.hpp:188
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::CurrentFE::nbLocalCoor
UInt nbLocalCoor() const
Getter for the number of local coordinates.
Definition: CurrentFE.hpp:425
ASSERT_PRE
#define ASSERT_PRE(X, A)
Definition: LifeAssert.hpp:96
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::CurrentFE::nbFEDof
UInt nbFEDof() const
Getter for the number of nodes.
Definition: CurrentFE.hpp:377
LifeV::nDimensions
const UInt nDimensions(NDIM)
LifeV::CurrentFE::nbQuadPt
UInt nbQuadPt() const
Getter for the number of quadrature nodes.
Definition: CurrentFE.hpp:365
LifeV::UInt
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191