LifeV
structure/fem/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/structure/fem/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 = new double[nDimensions];
1487  double* WI = new double[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 = new double[length];
1497  double* VL = new double[length];
1498 
1499  Int LWORK = 9;
1500  Int INFO = 0;
1501 
1502  double* WORK = new double[LWORK];
1503 
1504  double* A = new double[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, WI, 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  delete[] WR;
1522  delete[] WI;
1523  delete[] VL;
1524  delete[] VR;
1525  delete[] WORK;
1526  delete[] A;
1527 }
1528 
1529 //! ***********************************************************************************************
1530 //! METHODS FOR THE ST. VENANT KIRCHHOFF PENALIZED LAW
1531 //! ***********************************************************************************************
1532 //! Stiffness vector isochoric part ---------------------------------------------------------------
1533 
1534 // Source term : Int { ( \frac{lambda}{2} * Ic_iso - \frac{3}{2}*lambda - mu ) (F : \nabla v) - 1/3 * (Ic) * (F^-T : \nabla v) ) }
1535 void source_P1iso_VKPenalized ( Real lambda,
1536  Real mu,
1537  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1538  const boost::multi_array<Real, 3 >& Fk,
1539  const std::vector<Real>& Ic_isok,
1540  const std::vector<Real>& Ic_k,
1541  const std::vector<Real>& Jack_k,
1542  VectorElemental& elvec,
1543  const CurrentFE& fe )
1544 {
1545 
1546  Real s;
1547 
1548  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1549  {
1550  VectorElemental::vector_view vec = elvec.block ( icoor );
1551  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1552  {
1553  s = 0.0;
1554  for ( UInt k = 0; k < nDimensions; ++k )
1555  {
1556  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1557  {
1558  s += std::pow ( Jack_k[ ig ], (-2.0 / 3.0) ) * ( ( lambda / 2.0 ) * Ic_isok[ ig ] - (3.0 / 2.0) * lambda - mu ) *
1559  (Fk[ icoor ][ k ][ ig ] - (1.0 / 3.0) * Ic_k[ ig ] * FkMinusTransposed[ icoor ][ k ][ ig ] ) * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
1560 
1561  }
1562  }
1563  vec ( i ) += s;
1564  }
1565  }
1566 }
1567 
1568 
1569 // Source term : Int { ( 2 * mu * Jk^(-4.0/3.0) ) * ( (F*C : \nabla v) - 1/3 * (Ic_Squared) * (F^-T : \nabla v) ) }
1570 void source_P2iso_VKPenalized ( Real mu,
1571  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1572  const boost::multi_array<Real, 3 >& FkCk,
1573  const std::vector<Real>& Ic_Squared,
1574  const std::vector<Real>& Jk,
1575  VectorElemental& elvec,
1576  const CurrentFE& fe )
1577 {
1578 
1579  Real s;
1580 
1581  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1582  {
1583  VectorElemental::vector_view vec = elvec.block ( icoor );
1584  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1585  {
1586  s = 0.0;
1587  for ( UInt k = 0; k < nDimensions; ++k )
1588  {
1589  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1590  {
1591  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 );
1592  }
1593  }
1594  vec ( i ) += s;
1595  }
1596  }
1597 }
1598 
1599 //! Jacobian matrix of the first Piola-Kirchhoff tensor for the VK Penalized law
1600 //! 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 }
1601 void stiff_Jac_P1iso_VKPenalized_0term ( Real lambda, Real mu,
1602  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1603  const boost::multi_array<Real, 3 >& Fk,
1604  const std::vector<Real>& Jk ,
1605  const std::vector<Real>& Ic_k ,
1606  const std::vector<Real>& IcIso_k ,
1607  MatrixElemental& elmat,
1608  const CurrentFE& fe )
1609 {
1610  Real s;
1611 
1612  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1613  {
1614  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1615  {
1616  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1617  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1618  {
1619  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1620  {
1621  s = 0.0;
1622  for ( UInt l = 0; l < nDimensions; ++l )
1623  {
1624  for ( UInt k = 0; k < nDimensions; ++k )
1625  {
1626  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1627  {
1628  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) ) *
1629  fe.phiDer ( i, l, ig ) * ( Fk[ icoor ][ l ][ ig ] - (1.0 / 3.0) * Ic_k[ ig ] * FkMinusTransposed[ icoor ][ l ][ ig ] ) *
1630  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1631  }
1632  }
1633  }
1634  mat ( i, j ) += s;
1635  }
1636  }
1637  }
1638  }
1639 }
1640 
1641 //! 1. Stiffness term :int{ J^(-2/3) * (lambda / 2) * ( (-2/3) * Ic_k * J^(-2/3) * F^-T:\nabla \delta ) * ( F : \nabla \v )}
1642 void stiff_Jac_P1iso_VKPenalized_1term ( Real coeff,
1643  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1644  const boost::multi_array<Real, 3 >& Fk,
1645  const std::vector<Real>& Jk ,
1646  const std::vector<Real>& Ic_k ,
1647  MatrixElemental& elmat,
1648  const CurrentFE& fe )
1649 {
1650  Real s;
1651 
1652  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1653  {
1654  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1655  {
1656  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1657  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1658  {
1659  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1660  {
1661  s = 0.0;
1662  for ( UInt l = 0; l < nDimensions; ++l )
1663  {
1664  for ( UInt k = 0; k < nDimensions; ++k )
1665  {
1666  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1667  {
1668  s += ( -2.0 / 3.0 ) * Ic_k[ ig ] * std::pow ( Jk[ig], (-4.0 / 3.0) ) *
1669  fe.phiDer ( i, l, ig ) * Fk[ icoor ][ l ][ ig ] *
1670  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1671  }
1672  }
1673  }
1674  mat ( i, j ) += coeff * s;
1675  }
1676  }
1677  }
1678  }
1679 }
1680 
1681 //! 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 ) }
1682 void stiff_Jac_P1iso_VKPenalized_2term ( Real coef,
1683  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1684  const std::vector<Real>& Jk,
1685  const std::vector<Real>& Ic_k,
1686  MatrixElemental& elmat,
1687  const CurrentFE& fe )
1688 {
1689  Real s;
1690 
1691  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1692  {
1693  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1694  {
1695  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1696  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1697  {
1698  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1699  {
1700  s = 0.0;
1701  for ( UInt l = 0; l < nDimensions; ++l )
1702  {
1703  for ( UInt k = 0; k < nDimensions; ++k )
1704  {
1705  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1706  {
1707  s += ( ( 2.0 / 9.0 ) * std::pow ( Jk[ig], (-4.0 / 3.0) ) * Ic_k[ig] * Ic_k[ig] ) *
1708  FkMinusTransposed[ icoor ][ l ][ ig ] * fe.phiDer ( i, l, ig ) *
1709  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1710  }
1711  }
1712  }
1713  mat ( i, j ) += coef * s;
1714  }
1715  }
1716  }
1717  }
1718 }
1719 
1720 //! 2. Stiffness term : Int { J^(-2/3) * coeff * ( 2 * J^(-2/3) ) * ( F : \nabla \delta ) ( F : \nabla \v )}
1721 void stiff_Jac_P1iso_VKPenalized_3term ( Real coef,
1722  const boost::multi_array<Real, 3 >& Fk,
1723  const std::vector<Real>& Jk,
1724  MatrixElemental& elmat,
1725  const CurrentFE& fe )
1726 {
1727  Real s;
1728 
1729  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1730  {
1731  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1732  {
1733  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1734  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1735  {
1736  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1737  {
1738  s = 0.0;
1739  for ( UInt l = 0; l < nDimensions; ++l )
1740  {
1741  for ( UInt k = 0; k < nDimensions; ++k )
1742  {
1743  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1744  {
1745  s += 2.0 * std::pow ( Jk[ig], (-4.0 / 3.0) ) *
1746  fe.phiDer ( i, l, ig ) * Fk[ icoor ][ l ][ ig ] *
1747  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1748  }
1749  }
1750  }
1751  mat ( i, j ) += coef * s;
1752  }
1753  }
1754  }
1755  }
1756 }
1757 
1758 //! 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 )}
1759 void stiff_Jac_P1iso_VKPenalized_4term ( Real coef,
1760  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1761  const boost::multi_array<Real, 3 >& Fk,
1762  const std::vector<Real>& Jk ,
1763  const std::vector<Real>& Ic_k,
1764  MatrixElemental& elmat,
1765  const CurrentFE& fe )
1766 {
1767  Real s;
1768 
1769  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1770  {
1771  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1772  {
1773  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1774  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1775  {
1776  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1777  {
1778  s = 0.0;
1779  for ( UInt l = 0; l < nDimensions; ++l )
1780  {
1781  for ( UInt k = 0; k < nDimensions; ++k )
1782  {
1783  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1784  {
1785  s += ( ( -2.0 / 3.0 ) * std::pow ( Jk[ig], (-4.0 / 3.0) ) * Ic_k[ig] ) *
1786  fe.phiDer ( i, l, ig ) * FkMinusTransposed[ icoor ][ l ][ ig ] *
1787  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1788  }
1789  }
1790  }
1791  mat ( i, j ) += coef * s;
1792  }
1793  }
1794  }
1795  }
1796 }
1797 
1798 //! 5. Stiffness term : int { J^(-2.0/3.0) * ( (lambda/2) * Ic_isok - ( (3/2)*lambda + mu ) ) * \nabla \delta : \nabla v }
1799 void stiff_Jac_P1iso_VKPenalized_5term ( Real coef,
1800  Real secondCoef,
1801  const std::vector<Real>& Jk,
1802  const std::vector<Real>& Ic_isok,
1803  MatrixElemental& elmat,
1804  const CurrentFE& fe )
1805 {
1806  Real s;
1807 
1808  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
1809  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1810  {
1811  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1812  {
1813  s = 0.0;
1814  for ( UInt k = 0; k < nDimensions; ++k )
1815  {
1816  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1817  {
1818  s += std::pow ( Jk[ ig ], (-2.0 / 3.0) ) * ( (coef / 2.0) * Ic_isok[ ig ] - ( (3.0 / 2.0) * coef + secondCoef ) ) *
1819  fe.phiDer ( i, k, ig ) * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1820  }
1821  }
1822  mat_tmp ( i, j ) = s;
1823  }
1824  }
1825 
1826  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1827  {
1828  //! copy of diagonal block
1829  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
1830  mat += mat_tmp;
1831  }
1832 }
1833 
1834 //! 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 ) }
1835 void stiff_Jac_P1iso_VKPenalized_6term ( Real coef, Real secondCoef,
1836  const std::vector<Real>& Jk,
1837  const std::vector<Real>& Ic_isok,
1838  const boost::multi_array<Real, 3 >& Fk,
1839  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1840  MatrixElemental& elmat,
1841  const CurrentFE& fe )
1842 {
1843  Real s;
1844 
1845  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1846  {
1847  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1848  {
1849  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1850  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1851  {
1852  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1853  {
1854  s = 0.0;
1855  for ( UInt l = 0; l < nDimensions; ++l )
1856  {
1857  for ( UInt k = 0; k < nDimensions; ++k )
1858  {
1859  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1860  {
1861  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 ) ) *
1862  fe.phiDer ( i, l, ig ) * FkMinusTransposed[ icoor ][ l ][ ig ] *
1863  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1864  }
1865  }
1866  }
1867  mat ( i, j ) += s;
1868  }
1869  }
1870  }
1871  }
1872 }
1873 
1874 //! 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 }
1875 void stiff_Jac_P1iso_VKPenalized_7term ( Real coef,
1876  Real secondCoef,
1877  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1878  const std::vector<Real>& Ic_isok,
1879  const std::vector<Real>& Ic_k,
1880  const std::vector<Real>& Jk,
1881  MatrixElemental& elmat,
1882  const CurrentFE& fe )
1883 {
1884  Real s;
1885 
1886  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1887  {
1888  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1889  {
1890  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1891  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1892  {
1893  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1894  {
1895  s = 0.0;
1896  for ( UInt l = 0; l < nDimensions; ++l )
1897  {
1898  for ( UInt k = 0; k < nDimensions; ++k )
1899  {
1900  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1901  {
1902  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 ] ) *
1903  FkMinusTransposed[ icoor ][ l ][ ig ] * fe.phiDer ( j, l, ig ) *
1904  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
1905  }
1906  }
1907  }
1908  mat ( i, j ) += s;
1909  }
1910  }
1911  }
1912  }
1913 }
1914 
1915 //! 8. Stiffness term : int { ( -4.0/3.0) * ( mu * J^(-4/3) ) * ( F^-T: \grad \delta ) * ( F C ) : \nabla v }
1916 void stiff_Jac_P1iso_VKPenalized_8term ( Real coef, const std::vector<Real>& Jack_k,
1917  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1918  const boost::multi_array<Real, 3 >& FkCk,
1919  MatrixElemental& elmat,
1920  const CurrentFE& fe )
1921 {
1922 
1923  Real s;
1924 
1925  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1926  {
1927  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1928  {
1929  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1930  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1931  {
1932  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1933  {
1934  s = 0.0;
1935  for ( UInt l = 0; l < nDimensions; ++l )
1936  {
1937  for ( UInt k = 0; k < nDimensions; ++k )
1938  {
1939  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1940  {
1941  s += (-4.0 / 3.0) * std::pow ( Jack_k[ ig ], (-4.0 / 3.0) ) *
1942  fe.phiDer ( i, l, ig ) * FkCk[ icoor ][ l ][ ig ] *
1943  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1944  }
1945  }
1946  }
1947  mat ( i, j ) += s * coef;
1948  }
1949  }
1950  }
1951  }
1952 }
1953 
1954 
1955 //! 8. Stiffness term : int { ( 4.0/9.0) * ( mu * J^(-4/3) ) Ic_kSquared * (F^-T : \grad \delta ) * F^-T : \nabla \v }
1956 void stiff_Jac_P1iso_VKPenalized_9term ( Real coef, const std::vector<Real>& Jack_k,
1957  const std::vector<Real>& Ic_kSquared,
1958  const boost::multi_array<Real, 3 >& FkMinusTransposed,
1959  MatrixElemental& elmat,
1960  const CurrentFE& fe )
1961 {
1962  Real s;
1963  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
1964  {
1965  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
1966  {
1967  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
1968  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
1969  {
1970  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
1971  {
1972  s = 0.0;
1973  for ( UInt l = 0; l < nDimensions; ++l )
1974  {
1975  for ( UInt k = 0; k < nDimensions; ++k )
1976  {
1977  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
1978  {
1979  s += ( (-4.0 / 3.0) * std::pow ( Jack_k[ ig ], (-4.0 / 3.0) ) ) *
1980  fe.phiDer ( i, l, ig ) * (-1.0 / 3.0) * Ic_kSquared[ ig ] * FkMinusTransposed[ icoor ][ l ][ ig ] *
1981  FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
1982  }
1983  }
1984  }
1985  mat ( i, j ) += s * coef;
1986  }
1987  }
1988  }
1989  }
1990 }
1991 
1992 
1993 
1994 
1995 //! 10. Stiffness term : int { ( mu * J^(-4/3) ) * ( \nabla \delta * C ) : \nabla v }
1996 void stiff_Jac_P1iso_VKPenalized_10term ( Real coef,
1997  const std::vector<Real>& Jack_k,
1998  const boost::multi_array<Real, 3 >& Ck,
1999  MatrixElemental& elmat,
2000  const CurrentFE& fe )
2001 {
2002  Real s;
2003 
2004  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2005  {
2006  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
2007  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2008  {
2009  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2010  {
2011  s = 0.0;
2012  for ( UInt l = 0; l < nDimensions; ++l )
2013  {
2014  for ( UInt k = 0; k < nDimensions; ++k )
2015  {
2016  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2017  {
2018  s += std::pow ( Jack_k[ ig ], (-4.0 / 3.0) ) * Ck[ l ][ k ][ ig ] * fe.phiDer ( i, k, ig ) *
2019  fe.phiDer ( j, l, ig ) * fe.weightDet ( ig );
2020  }
2021  }
2022  }
2023  mat ( i, j ) += coef * s;
2024  }
2025  }
2026  }
2027 }
2028 
2029 //! 6. Stiffness term : int { ( mu * J^(-4/3) ) * (F [\nabla \delta]^T F ) : \nabla \v }
2030 void stiff_Jac_P1iso_VKPenalized_11term ( Real coef,
2031  const std::vector<Real>& Jk,
2032  const boost::multi_array<Real, 3 >& Fk,
2033  MatrixElemental& elmat,
2034  const CurrentFE& fe )
2035 {
2036  Real s;
2037 
2038  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2039  {
2040  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2041  {
2042  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2043  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2044  {
2045  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2046  {
2047  s = 0.0;
2048  for ( UInt l = 0; l < nDimensions; ++l )
2049  {
2050  for ( UInt k = 0; k < nDimensions; ++k )
2051  {
2052  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2053  {
2054  s += std::pow ( Jk[ ig ], (-4.0 / 3.0) ) * fe.phiDer ( i, l, ig ) * Fk[ l ][ jcoor ][ ig ] *
2055  Fk[ icoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2056  }
2057  }
2058  }
2059  mat ( i, j ) += coef * s;
2060  }
2061  }
2062  }
2063  }
2064 }
2065 
2066 //! 6. Stiffness term : int { ( mu * J^(-4/3) ) * (F * F^T * [\nabla \delta] ) : \nabla \v }
2067 void stiff_Jac_P1iso_VKPenalized_12term ( Real coef,
2068  const std::vector<Real>& Jk,
2069  const boost::multi_array<Real, 3 >& Fk,
2070  MatrixElemental& elmat,
2071  const CurrentFE& fe )
2072 {
2073  Real s;
2074 
2075  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2076  {
2077  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2078  {
2079  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2080  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2081  {
2082  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2083  {
2084  s = 0.0;
2085  for ( UInt p = 0; p < nDimensions; p++ )
2086  {
2087  for ( UInt k = 0; k < nDimensions; k++ )
2088  {
2089  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2090  {
2091  s += std::pow ( Jk[ig], (-4.0 / 3.0) ) * Fk[ icoor ][ p ][ ig ] * fe.phiDer ( i, k, ig ) *
2092  Fk[ jcoor ][ p ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2093  }
2094  }
2095  }
2096  mat ( i, j ) += coef * s;
2097  }
2098  }
2099  }
2100  }
2101 }
2102 
2103 //! 11. Stiffness term : int { ( mu * J^(-4/3) * ( (1/3) * Ic_SquaredK * ( F^-T [\nabla \delta ]^T F^-T) : \nabla v ) }
2104 void stiff_Jac_P1iso_VKPenalized_13term ( Real coef,
2105  const std::vector<Real>& Jk,
2106  const std::vector<Real>& Ic_kSquared,
2107  const boost::multi_array<Real, 3 >& FkMinusTransposed,
2108  MatrixElemental& elmat,
2109  const CurrentFE& fe )
2110 {
2111  Real s;
2112 
2113  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2114  {
2115  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2116  {
2117  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2118  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2119  {
2120  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2121  {
2122  s = 0.0;
2123  for ( UInt l = 0; l < nDimensions; ++l )
2124  {
2125  for ( UInt k = 0; k < nDimensions; ++k )
2126  {
2127  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2128  {
2129  s += ( std::pow ( Jk[ig], (-4.0 / 3.0) ) ) * ( Ic_kSquared[ig] / 3.0 ) *
2130  fe.phiDer ( j, l, ig ) * FkMinusTransposed[ icoor ][ l ][ ig ] *
2131  fe.phiDer ( i, k, ig ) * FkMinusTransposed[ jcoor ][ k ][ ig ] * fe.weightDet ( ig );
2132 
2133  }
2134  }
2135  }
2136  mat ( i, j ) += s * coef;
2137  }
2138  }
2139  }
2140  }
2141 }
2142 
2143 //! 12. Stiffness term : int { ( mu * J^(-4/3) ) * ( (-4.0/3.0) * ( FkCk : \nabla \delta ) ) * F^-T : \nabla v ) }
2144 void stiff_Jac_P1iso_VKPenalized_14term ( Real coef,
2145  const std::vector<Real>& Jk,
2146  const boost::multi_array<Real, 3 >& FkCk,
2147  const boost::multi_array<Real, 3 >& FkMinusTransposed,
2148  MatrixElemental& elmat,
2149  const CurrentFE& fe )
2150 {
2151  Real s;
2152 
2153  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2154  {
2155  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2156  {
2157  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2158  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2159  {
2160  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2161  {
2162  s = 0.0;
2163  for ( UInt l = 0; l < nDimensions; ++l )
2164  {
2165  for ( UInt k = 0; k < nDimensions; ++k )
2166  {
2167  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2168  {
2169  s += ( std::pow ( Jk[ig], (-4.0 / 3.0) ) ) * ( -4.0 / 3.0 ) *
2170  fe.phiDer ( i, l, ig ) * FkMinusTransposed[ icoor ][ l ][ ig ] *
2171  ( FkCk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) ) * fe.weightDet ( ig );
2172  }
2173  }
2174  }
2175  mat ( i, j ) += s * coef;
2176  }
2177  }
2178  }
2179  }
2180 }
2181 
2182 //! ***********************************************************************************************
2183 //! END OF VENANT-KIRCHHOFF PENALIZED MODEL
2184 //! ***********************************************************************************************
2185 
2186 //! ***********************************************************************************************
2187 //! SECOND ORDER EXPONENTIAL MODEL
2188 //! ***********************************************************************************************
2189 //! Methods for the Second Order Exponential law
2190 //! Assemble the first Piola-Kirchhoff tensor
2191 //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) ) }
2192 void source_P1iso_SecondOrderExponential ( Real coef, Real coefExp,
2193  const boost::multi_array<Real, 3 >& CofFk,
2194  const boost::multi_array<Real, 3 >& Fk,
2195  const std::vector<Real>& Jk,
2196  const std::vector<Real>& Ic_isok,
2197  VectorElemental& elvec,
2198  const CurrentFE& fe )
2199 {
2200 
2201  Real s;
2202 
2203  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2204  {
2205  VectorElemental::vector_view vec = elvec.block ( icoor );
2206  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2207  {
2208  s = 0.0;
2209  for ( UInt k = 0; k < nDimensions; ++k )
2210  {
2211  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2212  {
2213  s += (Ic_isok[ ig ] - 3.0) * exp ( coefExp * (Ic_isok[ ig ] - 3.0) * (Ic_isok[ ig ] - 3.0) ) *
2214  (std::pow ( Jk[ ig ], (-2.0 / 3.0) ) * Fk[ icoor ][ k ][ ig ]
2215  - 1.0 / 3.0 * ( Ic_isok[ ig ] / Jk[ ig ] ) * CofFk[ icoor ][ k ][ ig ] ) * fe.phiDer ( i, k, ig ) * fe.weightDet ( ig );
2216  }
2217  }
2218  vec ( i ) += s * coef;
2219  }
2220  }
2221 }
2222 
2223 //! 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 ) }
2224 // In this method, the coef that is passed is -4/3 * alpha, therefore it can multiply the sum at the end.
2225 void stiff_Jac_P1iso_SecondOrderExp_1term ( Real coef,
2226  Real coefExp,
2227  const boost::multi_array<Real, 3 >& CofFk,
2228  const boost::multi_array<Real, 3 >& Fk,
2229  const std::vector<Real>& Jk ,
2230  const std::vector<Real>& Ic_isok,
2231  MatrixElemental& elmat,
2232  const CurrentFE& fe )
2233 {
2234  Real s;
2235 
2236  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2237  {
2238  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2239  {
2240  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2241  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2242  {
2243  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2244  {
2245  s = 0.0;
2246  for ( UInt l = 0; l < nDimensions; ++l )
2247  {
2248  for ( UInt k = 0; k < nDimensions; ++k )
2249  {
2250  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2251  {
2252  s += ( std::pow ( Jk[ig], -5.0 / 3.0) * exp (coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) ) *
2253  ( (Ic_isok[ig] - 3.0) * ( 1.0 + 2 * coefExp * (Ic_isok[ig] - 3.0) * Ic_isok[ig] ) + Ic_isok[ig] ) *
2254  fe.phiDer ( i, l, ig ) * Fk[ icoor ][ l ][ ig ] *
2255  CofFk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2256  }
2257  }
2258  }
2259  mat ( i, j ) += coef * s;
2260  }
2261  }
2262  }
2263  }
2264 }
2265 
2266 //! 2. Stiffness term : Int { 4 * coef * J^(-4/3) * exp( coefExp*( Ic_iso - 3)*( Ic_iso - 3) ) *
2267 //! ( 1.0 + 2 * coefExp * (Ic_isoK - 3)^2 ) * ( F : \nabla \delta ) ( F : \nabla \v )}
2268 // When the method is called, the coef parameter stores already 4 * alpha. This is why it is used at the end.
2269 void stiff_Jac_P1iso_SecondOrderExp_2term ( Real coef,
2270  Real coefExp,
2271  const boost::multi_array<Real, 3 >& Fk,
2272  const std::vector<Real>& Jk,
2273  const std::vector<Real>& Ic_isok,
2274  MatrixElemental& elmat,
2275  const CurrentFE& fe )
2276 {
2277  Real s;
2278 
2279  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2280  {
2281  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2282  {
2283  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2284  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2285  {
2286  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2287  {
2288  s = 0.0;
2289  for ( UInt l = 0; l < nDimensions; ++l )
2290  {
2291  for ( UInt k = 0; k < nDimensions; ++k )
2292  {
2293  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2294  {
2295  s += std::pow ( Jk[ig], -4.0 / 3.0 ) * exp ( coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) *
2296  ( 1.0 + 2 * coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) *
2297  fe.phiDer ( i, l, ig ) * Fk[ icoor ][ l ][ ig ] *
2298  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2299  }
2300  }
2301  }
2302  mat ( i, j ) += coef * s;
2303  }
2304  }
2305  }
2306  }
2307 }
2308 
2309 //! 3. Stiffness term:
2310 //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) ) ) *
2311 // ( CofF : \nabla \delta ) ( CofF : \nabla \v )}
2312 // The coef that is passed stores 4/9 * alpha.
2313 void stiff_Jac_P1iso_SecondOrderExp_3term ( Real coef, Real coefExp,
2314  const boost::multi_array<Real, 3 >& CofFk,
2315  const std::vector<Real>& Jk,
2316  const std::vector<Real>& Ic_isok,
2317  MatrixElemental& elmat,
2318  const CurrentFE& fe )
2319 {
2320  Real s;
2321 
2322  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2323  {
2324  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2325  {
2326  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2327  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2328  {
2329  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2330  {
2331  s = 0.0;
2332  for ( UInt l = 0; l < nDimensions; ++l )
2333  {
2334  for ( UInt k = 0; k < nDimensions; ++k )
2335  {
2336  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2337  {
2338  s += ( ( 1 / (Jk[ig] * Jk[ig]) ) * Ic_isok[ig] * exp ( coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) ) *
2339  ( (Ic_isok[ig] - 3.0) * ( 1.0 + 2.0 * coefExp * (Ic_isok[ig] - 3.0) * Ic_isok[ig] ) + Ic_isok[ig] ) *
2340  CofFk[ icoor ][ l ][ ig ] * fe.phiDer ( i, l, ig ) *
2341  CofFk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2342  }
2343  }
2344  }
2345  mat ( i, j ) += coef * s;
2346  }
2347  }
2348  }
2349  }
2350 }
2351 
2352 
2353 //! 4. Stiffness term:
2354 //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) ) *
2355 // ( F : \nabla \delta ) ( CofF : \nabla \v ) }
2356 // As in other methods, the coef that is passed is -4.0/3.0 * alpha.
2357 void stiff_Jac_P1iso_SecondOrderExp_4term ( Real coef, Real coefExp,
2358  const boost::multi_array<Real, 3 >& CofFk,
2359  const boost::multi_array<Real, 3 >& Fk,
2360  const std::vector<Real>& Jk ,
2361  const std::vector<Real>& Ic_isok,
2362  const std::vector<Real>& /*Ic_k*/,
2363  MatrixElemental& elmat,
2364  const CurrentFE& fe )
2365 {
2366  Real s;
2367 
2368  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2369  {
2370  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2371  {
2372  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2373  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2374  {
2375  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2376  {
2377  s = 0.0;
2378  for ( UInt l = 0; l < nDimensions; ++l )
2379  {
2380  for ( UInt k = 0; k < nDimensions; ++k )
2381  {
2382  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2383  {
2384  s += std::pow ( Jk[ig], (-5.0 / 3.0) ) * exp ( coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) *
2385  ( Ic_isok[ig] + (Ic_isok[ig] - 3.0) * ( 1.0 + 2.0 * coefExp * (Ic_isok[ig] - 3.0) * Ic_isok[ig] ) ) *
2386  fe.phiDer ( i, l, ig ) * CofFk[ icoor ][ l ][ ig ] *
2387  Fk[ jcoor ][ k ][ ig ] * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2388  }
2389  }
2390  }
2391  mat ( i, j ) += coef * s;
2392  }
2393  }
2394  }
2395  }
2396 }
2397 
2398 //! 5. Stiffness term :
2399 //Int {coef * J^(-2/3) * exp( coefExp*( Ic_iso - 3)*( Ic_iso - 3)) * ( Ic_iso - 3)* (\nabla \delta: \nabla \v)}
2400 void stiff_Jac_P1iso_SecondOrderExp_5term ( Real coef,
2401  Real coefExp,
2402  const std::vector<Real>& Jk,
2403  const std::vector<Real>& Ic_isok,
2404  MatrixElemental& elmat,
2405  const CurrentFE& fe )
2406 {
2407  Real s;
2408 
2409  MatrixElemental::matrix_type mat_tmp ( fe.nbFEDof(), fe.nbFEDof() );
2410  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2411  {
2412  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2413  {
2414  s = 0.0;
2415  for ( UInt k = 0; k < nDimensions; ++k )
2416  {
2417  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2418  {
2419  s += (Ic_isok[ ig ] - 3.0) * exp ( coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) *
2420  std::pow ( Jk[ig], -2.0 / 3.0 ) * fe.phiDer ( i, k, ig ) * fe.phiDer ( j, k, ig ) * fe.weightDet ( ig );
2421  }
2422  }
2423  mat_tmp ( i, j ) = coef * s;
2424  }
2425  }
2426 
2427  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2428  {
2429  //! copy of diagonal block
2430  MatrixElemental::matrix_view mat = elmat.block ( icoor, icoor );
2431  mat += mat_tmp;
2432  }
2433 }
2434 
2435 //! 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 }
2436 void stiff_Jac_P1iso_SecondOrderExp_6term ( Real coef,
2437  Real coefExp,
2438  const boost::multi_array<Real, 3 >& CofFk,
2439  const std::vector<Real>& Jk,
2440  const std::vector<Real>& Ic_isok,
2441  MatrixElemental& elmat,
2442  const CurrentFE& fe )
2443 {
2444  Real s;
2445 
2446  for ( UInt icoor = 0; icoor < nDimensions; ++icoor )
2447  {
2448  for ( UInt jcoor = 0; jcoor < nDimensions; ++jcoor )
2449  {
2450  MatrixElemental::matrix_view mat = elmat.block ( icoor, jcoor );
2451  for ( UInt i = 0; i < fe.nbFEDof(); ++i )
2452  {
2453  for ( UInt j = 0; j < fe.nbFEDof(); ++j )
2454  {
2455  s = 0.0;
2456  for ( UInt l = 0; l < nDimensions; ++l )
2457  {
2458  for ( UInt k = 0; k < nDimensions; ++k )
2459  {
2460  for ( UInt ig = 0; ig < fe.nbQuadPt(); ++ig )
2461  {
2462  s += ( 1 / (Jk[ig] * Jk[ig]) ) * (Ic_isok[ig] - 3.0) *
2463  exp ( coefExp * (Ic_isok[ig] - 3.0) * (Ic_isok[ig] - 3.0) ) * Ic_isok[ig] *
2464  fe.phiDer ( j, l, ig ) * CofFk[ icoor ][ l ][ ig ] *
2465  fe.phiDer ( i, k, ig ) * CofFk[ jcoor ][ k ][ ig ] * fe.weightDet ( ig );
2466  }
2467  }
2468  }
2469  mat ( i, j ) += s * coef;
2470  }
2471  }
2472  }
2473  }
2474 }
2475 
2476 //! ***********************************************************************************************
2477 //! END OF SECOND ORDER EXPONENTIAL MODEL
2478 //! ***********************************************************************************************
2479 
2480 } //! End namespace AssemblyElementalStructure
2481 
2482 } //! End namespace LifeV
2483 
2484 
2485 #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