LifeV
CurrentFE.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 CurrentFE class implementation
30 
31  @author Jean-Frederic Gerbeau
32  @date 00-04-2002
33 
34  @contributor V. Martin
35  Samuel Quinodoz <samuel.quinodoz@epfl.ch>
36  @mantainer Samuel Quinodoz <samuel.quinodoz@epfl.ch>
37  */
38 
39 #include <lifev/core/fem/CurrentFE.hpp>
40 
41 namespace LifeV
42 {
43 CurrentFE::CurrentFE ( const ReferenceFE& refFE, const GeometricMap& geoMap, const QuadratureRule& qr )
44  :
45  M_nbNode ( refFE.nbDof() ),
46  M_nbLocalCoor ( refFE.nbLocalCoor() ),
47  M_nbDiag ( refFE.nbDiag() ),
48  M_nbUpper ( refFE.nbUpper() ) ,
49  M_nbPattern ( refFE.nbPattern() ),
50 
51  M_nbGeoNode ( geoMap.nbDof() ),
52  M_nbQuadPt ( qr.nbQuadPt() ),
53 
54  M_refFE ( &refFE ),
55  M_geoMap ( &geoMap),
56  M_quadRule (new QuadratureRule (qr) ),
57 
58 
59  M_cellNodes (boost::extents[geoMap.nbDof()][nDimensions]),
60  M_quadNodes (boost::extents[M_nbQuadPt][nDimensions]),
61 
62  M_dphiGeometricMap (boost::extents[M_nbGeoNode][M_nbLocalCoor][M_nbQuadPt]),
63  M_jacobian (boost::extents[M_nbLocalCoor][M_nbLocalCoor][M_nbQuadPt]),
64  M_detJacobian (boost::extents[M_nbQuadPt]),
65  M_wDetJacobian (boost::extents[M_nbQuadPt]),
66  M_tInverseJacobian (boost::extents[M_nbLocalCoor][M_nbLocalCoor][M_nbQuadPt]),
67 
68  M_phi (boost::extents[M_nbNode][M_refFE->feDim()][M_nbQuadPt]),
69  M_dphi (boost::extents[M_nbNode][M_nbLocalCoor][M_nbQuadPt]),
70  M_d2phi (boost::extents[M_nbNode][M_nbLocalCoor][M_nbLocalCoor][M_nbQuadPt]),
71  M_phiVect (boost::extents[M_nbNode][M_refFE->feDim()][M_nbQuadPt]),
72 
73  M_dphiRef (boost::extents[M_nbNode][M_nbLocalCoor][M_nbQuadPt]),
74  M_d2phiRef (boost::extents[M_nbNode][M_nbLocalCoor][M_nbLocalCoor][M_nbQuadPt]),
75  M_divPhiRef (boost::extents[M_nbNode][M_nbQuadPt]),
76 
77  M_cellNodesUpdated (false),
78  M_quadNodesUpdated (false),
79 
80  M_dphiGeometricMapUpdated (false),
81  M_jacobianUpdated (false),
82  M_detJacobianUpdated (false),
83  M_wDetJacobianUpdated (false),
84  M_tInverseJacobianUpdated (false),
85 
86  M_phiUpdated (false),
87  M_dphiUpdated (false),
88  M_d2phiUpdated (false),
89  M_phiVectUpdated (false),
90 
91  M_dphiRefUpdated (false),
92  M_d2phiRefUpdated (false),
93  M_divPhiRefUpdated (false)
94 
95 {
96  for ( UInt iterQuad (0); iterQuad < M_nbQuadPt; ++iterQuad )
97  {
98  for ( UInt iterNode (0); iterNode < M_nbNode; ++iterNode )
99  {
100  // --- PHI ---
101  if (M_refFE->hasPhi() )
102  {
103  for (UInt iterFEDim (0); iterFEDim < M_refFE->feDim(); ++iterFEDim)
104  {
105  M_phi[iterNode][iterFEDim][iterQuad] = M_refFE->phi ( iterNode, iterFEDim,
106  M_quadRule->quadPointCoor (iterQuad) );
107  }
108  }
109 
110  // --- DIV ---
111  if (M_refFE->hasDivPhi() )
112  {
113  M_divPhiRef[iterNode][iterQuad] = M_refFE->divPhi (iterNode,
114  M_quadRule->quadPointCoor (iterQuad) );
115  };
116 
117 
118  for ( UInt icoor (0); icoor < M_nbLocalCoor; ++icoor )
119  {
120  // --- DPHI ---
121  if (M_refFE->hasDPhi() )
122  {
123  M_dphiRef[iterNode][icoor][iterQuad] = M_refFE->dPhi ( iterNode,
124  icoor,
125  M_quadRule->quadPointCoor (iterQuad) );
126  };
127 
128  // --- D2PHI
129  if (M_refFE->hasD2Phi() )
130  {
131  for ( UInt jcoor (0); jcoor < M_nbLocalCoor; ++jcoor )
132  {
133  M_d2phiRef[iterNode][icoor][jcoor][iterQuad] = M_refFE->d2Phi ( iterNode,
134  icoor, jcoor,
135  M_quadRule->quadPointCoor (iterQuad) );
136  }
137  }
138  }
139  }
140  for ( UInt k (0); k < M_nbGeoNode; ++k )
141  {
142  for ( UInt icoor (0); icoor < M_nbLocalCoor; ++icoor )
143  {
144  M_dphiGeometricMap[k][icoor][iterQuad] = M_geoMap->dPhi ( k, icoor,
145  M_quadRule->quadPointCoor (iterQuad) );
146  }
147  }
148  }
149 
150  M_phiUpdated = true;
151  M_dphiRefUpdated = true;
152  M_divPhiRefUpdated = true;
153  M_d2phiRefUpdated = true;
154  M_dphiGeometricMapUpdated = true;
155 }
156 
157 CurrentFE::CurrentFE ( const ReferenceFE& refFE, const GeometricMap& geoMap ) :
158  M_nbNode (refFE.nbDof() ),
159  M_nbLocalCoor (refFE.nbLocalCoor() ),
160  M_nbDiag (refFE.nbDiag() ),
161  M_nbUpper (refFE.nbUpper() ),
162  M_nbPattern (refFE.nbPattern() ),
163 
164  M_nbGeoNode (geoMap.nbDof() ),
165  M_nbQuadPt (0),
166 
167  M_refFE (&refFE),
168  M_geoMap (&geoMap),
169  M_quadRule (0),
170 
171  M_cellNodes (boost::extents[geoMap.nbDof()][nDimensions]),
172 
173 
174  M_cellNodesUpdated (false),
175  M_quadNodesUpdated (false),
176  M_dphiGeometricMapUpdated (false),
177  M_jacobianUpdated (false),
178  M_detJacobianUpdated (false),
179  M_wDetJacobianUpdated (false),
180  M_tInverseJacobianUpdated (false),
181  M_phiUpdated (false),
182  M_dphiUpdated (false),
183  M_d2phiUpdated (false),
184  M_phiVectUpdated (false),
185  M_dphiRefUpdated (false),
186  M_d2phiRefUpdated (false),
187  M_divPhiRefUpdated (false)
188 
189 {
190  // Nothing to be done here
191 }
192 
193 CurrentFE::CurrentFE (const CurrentFE& fe) :
194  M_nbNode (fe.M_nbNode),
195  M_nbLocalCoor (fe.M_nbLocalCoor),
196  M_nbDiag (fe.M_nbDiag),
197  M_nbUpper (fe.M_nbUpper) ,
198  M_nbPattern (fe.M_nbPattern),
199  M_currentId (fe.M_currentId),
200  M_currentLocalId (fe.M_currentLocalId),
201  M_nbGeoNode (fe.M_nbGeoNode),
202  M_nbQuadPt (fe.M_nbQuadPt),
203 
204  M_refFE (fe.M_refFE),
205  M_geoMap (fe.M_geoMap),
206  M_quadRule (new QuadratureRule (*fe.M_quadRule) ),
207 
208  M_cellNodes (fe.M_cellNodes),
209  M_quadNodes (fe.M_quadNodes),
210  M_dphiGeometricMap (fe.M_dphiGeometricMap),
211  M_jacobian (fe.M_jacobian),
212  M_detJacobian (fe.M_detJacobian),
213  M_wDetJacobian (fe.M_wDetJacobian),
214  M_tInverseJacobian (fe.M_tInverseJacobian),
215  M_phi (fe.M_phi),
216  M_dphi (fe.M_dphi),
217  M_d2phi (fe.M_d2phi),
218  M_phiVect (fe.M_phiVect),
219  M_dphiRef (fe.M_dphiRef),
220  M_d2phiRef (fe.M_d2phiRef),
221  M_divPhiRef (fe.M_divPhiRef),
222 
223  M_cellNodesUpdated (fe.M_cellNodesUpdated),
224  M_quadNodesUpdated (fe.M_quadNodesUpdated),
225  M_dphiGeometricMapUpdated (fe.M_dphiGeometricMapUpdated),
226  M_jacobianUpdated (fe.M_jacobianUpdated),
227  M_detJacobianUpdated (fe.M_detJacobianUpdated),
228  M_wDetJacobianUpdated (fe.M_wDetJacobianUpdated),
229  M_tInverseJacobianUpdated (fe.M_tInverseJacobianUpdated),
230  M_phiUpdated (fe.M_phiUpdated),
231  M_dphiUpdated (fe.M_dphiUpdated),
232  M_d2phiUpdated (fe.M_d2phiUpdated),
233  M_phiVectUpdated (fe.M_phiVectUpdated),
234  M_dphiRefUpdated (fe.M_dphiRefUpdated),
235  M_d2phiRefUpdated (fe.M_d2phiRefUpdated),
236  M_divPhiRefUpdated (fe.M_divPhiRefUpdated)
237 {
238  // Nothing to be done here
239 }
240 
241 //----------------------------------------------------------------------
242 
243 // Added by S. Quinodoz
244 void CurrentFE::coorBackMap (Real x, Real y, Real z,
245  Real& xi, Real& eta, Real& zeta) const
246 {
247  // If the Mapping is not P1, this simple "back mapping" does not
248  // work properly, as it suppose that the mapping is P1. If another
249  // mapping (higher order, or for different shape like Q1) is needed
250  // a Newton algorithm has to be implemented to solve the problem
251  // M(x)-p = 0
252  // with M the mapping, p the given point and x its coordinates in
253  // the reference frame. However, this is more expensive and inexact.
254 
255  // ASSERT(geoMap.name.compare("P1")==0," The 'coorBackMap' function should NOT be used in this case! ");
256 
257  // In the P1 mapping case, we want to solve the problem
258  // P = \lamda_1 V_1 + \lambda_2 V_2 + \lambda_3 V_3
259  // where V_i is the ith vertice of the cell because we have then
260  // \hat{P} = \lambda_1 \hat{V}_1 + \lambda_2 \hat{V}_2 + \lambda_3 \hat{V}_3
261  // where \hat{P} is our point in the reference frame
262  //
263  // For "simplicity", we use Cramer's formula (no need for an heavy solver)
264 
265 
266  if (M_nbLocalCoor == 1)
267  {
268 
269  // 1D Case
270 
271  const Real a11 = M_cellNodes[1][0] - M_cellNodes[0][0];
272  const Real b1 = x - M_cellNodes[0][0];
273 
274  Real lambda = b1 / a11;
275 
276  xi = lambda * 1; // 1= refFE.xi(1)-refFE.xi(0);
277  eta = 0;
278  zeta = 0;
279 
280  }
281  else if ( M_nbLocalCoor == 2)
282  {
283 
284  // 2D Case
285 
286  const Real a11 = M_cellNodes[1][0] - M_cellNodes[0][0];
287  const Real a12 = M_cellNodes[2][0] - M_cellNodes[0][0];
288  const Real a21 = M_cellNodes[1][1] - M_cellNodes[0][1];
289  const Real a22 = M_cellNodes[2][1] - M_cellNodes[0][1];
290  const Real b1 = x - M_cellNodes[0][0];
291  const Real b2 = y - M_cellNodes[0][1];
292 
293  Real D = a11 * a22 - a21 * a12;
294  Real lambda_1 = b1 * a22 - b2 * a12;
295  Real lambda_2 = a11 * b2 - a21 * b1;
296 
297  lambda_1 /= D;
298  lambda_2 /= D;
299 
300  xi = lambda_1 * 1 + lambda_2 * 0; // 1=refFE.xi(1) -refFE.xi(0) ; 0=refFE.xi(2) -refFE.xi(0);
301  eta = lambda_1 * 0 + lambda_2 * 1; // 0=refFE.eta(1) -refFE.eta(0) ; 1=refFE.eta(2) -refFE.eta(0);
302  zeta = 0;
303 
304  }
305  else if (M_nbLocalCoor == 3)
306  {
307 
308  // 3D Case
309 
310  const Real a11 = M_cellNodes[1][0] - M_cellNodes[0][0];
311  const Real a12 = M_cellNodes[2][0] - M_cellNodes[0][0];
312  const Real a13 = M_cellNodes[3][0] - M_cellNodes[0][0];
313  const Real a21 = M_cellNodes[1][1] - M_cellNodes[0][1];
314  const Real a22 = M_cellNodes[2][1] - M_cellNodes[0][1];
315  const Real a23 = M_cellNodes[3][1] - M_cellNodes[0][1];
316  const Real a31 = M_cellNodes[1][2] - M_cellNodes[0][2];
317  const Real a32 = M_cellNodes[2][2] - M_cellNodes[0][2];
318  const Real a33 = M_cellNodes[3][2] - M_cellNodes[0][2];
319  const Real b1 = x - M_cellNodes[0][0];
320  const Real b2 = y - M_cellNodes[0][1];
321  const Real b3 = z - M_cellNodes[0][2];
322 
323  const Real D = a11 * a22 * a33 + a31 * a12 * a23 + a21 * a32 * a13 - a11 * a32 * a23 - a31 * a22 * a13 - a21 * a12 * a33;
324  Real lambda_1 = b1 * a22 * a33 + b3 * a12 * a23 + b2 * a32 * a13 - b1 * a32 * a23 - b3 * a22 * a13 - b2 * a12 * a33;
325  Real lambda_2 = a11 * b2 * a33 + a31 * b1 * a23 + a21 * b3 * a13 - a11 * b3 * a23 - a31 * b2 * a13 - a21 * b1 * a33;
326  Real lambda_3 = a11 * a22 * b3 + a31 * a12 * b2 + a21 * a32 * b1 - a11 * a32 * b2 - a31 * a22 * b1 - a21 * a12 * b3;
327 
328  lambda_1 /= D;
329  lambda_2 /= D;
330  lambda_3 /= D;
331 
332  xi = lambda_1 * 1 + lambda_2 * 0 + lambda_3 * 0;
333  eta = lambda_1 * 0 + lambda_2 * 1 + lambda_3 * 0;
334  zeta = lambda_1 * 0 + lambda_2 * 0 + lambda_3 * 1;
335 
336  }
337  else
338  {
339  ERROR_MSG ("Impossible dimension to invert coordinates");
340  }
341 
342 }
343 
344 // Compute the Jacobian at the given point
345 // this means that jacobian(P1,P2,P3,i,j)
346 // computes d x_i / d zeta_j
347 // where x are the global coordinates
348 // zeta the reference coordinates
349 
350 Real CurrentFE::pointJacobian (Real hat_x, Real hat_y, Real hat_z,
351  Int comp_x, Int comp_zeta) const
352 {
353  Real jac (0);
354  GeoVector P (3);
355  P[0] = hat_x;
356  P[1] = hat_y;
357  P[2] = hat_z;
358 
359  for ( UInt i = 0; i < M_nbGeoNode; i++ )
360  {
361  jac += M_cellNodes[i][comp_x] * M_geoMap->dPhi ( i , comp_zeta , P );
362  };
363 
364  return jac;
365 }
366 
367 // Compute the determinant of the Jacobian at the given point
368 Real CurrentFE::pointDetJacobian (Real hat_x, Real hat_y, Real hat_z) const
369 {
370  Real det (0);
371 
372  switch (M_nbLocalCoor)
373  {
374  case 1:
375  {
376  Real a ( pointJacobian (hat_x, hat_y, hat_z, 0, 0) );
377 
378  det = a;
379 
380  break;
381  }
382  case 2:
383  {
384  Real a ( pointJacobian (hat_x, hat_y, hat_z, 0, 0) );
385  Real b ( pointJacobian (hat_x, hat_y, hat_z, 0, 1) );
386  Real c ( pointJacobian (hat_x, hat_y, hat_z, 1, 0) );
387  Real d ( pointJacobian (hat_x, hat_y, hat_z, 1, 1) );
388 
389  det = a * d - b * c;
390 
391  break;
392  }
393  case 3:
394  {
395  Real a ( pointJacobian (hat_x, hat_y, hat_z, 0, 0) );
396  Real b ( pointJacobian (hat_x, hat_y, hat_z, 0, 1) );
397  Real c ( pointJacobian (hat_x, hat_y, hat_z, 0, 2) );
398  Real d ( pointJacobian (hat_x, hat_y, hat_z, 1, 0) );
399  Real e ( pointJacobian (hat_x, hat_y, hat_z, 1, 1) );
400  Real f ( pointJacobian (hat_x, hat_y, hat_z, 1, 2) );
401  Real g ( pointJacobian (hat_x, hat_y, hat_z, 2, 0) );
402  Real h ( pointJacobian (hat_x, hat_y, hat_z, 2, 1) );
403  Real i ( pointJacobian (hat_x, hat_y, hat_z, 2, 2) );
404 
405  Real ei (e * i);
406  Real fh (f * h);
407  Real bi (b * i);
408  Real ch (c * h);
409  Real bf (b * f);
410  Real ce (c * e);
411 
412  det = a * (ei - fh) + d * (ch - bi) + g * ( bf - ce);
413 
414  break;
415  }
416  default:
417  ERROR_MSG ( "Dimension (M_nbLocalCoor): only 1, 2 or 3!" );
418  break;
419  }
420 
421  return det;
422 }
423 
424 Real CurrentFE::pointInverseJacobian (Real hat_x, Real hat_y, Real hat_z,
425  Int comp_x, Int comp_zeta) const
426 {
427  if ( M_nbLocalCoor == 1 )
428  {
429 
430  return 1 / pointJacobian (hat_x, hat_y, hat_z, comp_x, comp_zeta);
431  }
432  else if (M_nbLocalCoor == 2)
433  {
434 
435  Real a11 = pointJacobian (hat_x, hat_y, hat_z, 0, 0);
436  Real a12 = pointJacobian (hat_x, hat_y, hat_z, 0, 1);
437  Real a21 = pointJacobian (hat_x, hat_y, hat_z, 1, 0);
438  Real a22 = pointJacobian (hat_x, hat_y, hat_z, 1, 1);
439 
440  Int total (comp_x + comp_zeta);
441  Int mysign (1);
442  if (total % 2 == 1)
443  {
444  mysign = -1;
445  }
446 
447  return mysign * pointJacobian (hat_x, hat_y, hat_z, comp_zeta, comp_x) / (a11 * a22 - a12 * a21);
448 
449  }
450  else if (M_nbLocalCoor == 3)
451  {
452 
453  Real a11 = pointJacobian (hat_x, hat_y, hat_z, 0, 0);
454  Real a12 = pointJacobian (hat_x, hat_y, hat_z, 0, 1);
455  Real a13 = pointJacobian (hat_x, hat_y, hat_z, 0, 2);
456  Real a21 = pointJacobian (hat_x, hat_y, hat_z, 1, 0);
457  Real a22 = pointJacobian (hat_x, hat_y, hat_z, 1, 1);
458  Real a23 = pointJacobian (hat_x, hat_y, hat_z, 1, 2);
459  Real a31 = pointJacobian (hat_x, hat_y, hat_z, 2, 0);
460  Real a32 = pointJacobian (hat_x, hat_y, hat_z, 2, 1);
461  Real a33 = pointJacobian (hat_x, hat_y, hat_z, 2, 2);
462 
463  Real det = a11 * a22 * a33 + a12 * a23 * a31 + a13 * a21 * a32 - a11 * a23 * a32 - a13 * a22 * a31 - a12 * a21 * a33;
464 
465  std::vector<std::vector< Real > > cof (3, std::vector<Real> (3, 0) );
466 
467  cof[0][0] = (a22 * a33 - a23 * a32);
468  cof[0][1] = - (a21 * a33 - a31 * a23);
469  cof[0][2] = (a21 * a32 - a31 * a22);
470  cof[1][0] = - (a12 * a33 - a32 * a13);
471  cof[1][1] = (a11 * a33 - a13 * a31);
472  cof[1][2] = - (a11 * a32 - a31 * a12);
473  cof[2][0] = (a12 * a23 - a13 * a22);
474  cof[2][1] = - (a11 * a23 - a13 * a21);
475  cof[2][2] = (a11 * a22 - a12 * a21);
476 
477  // inverse need the TRANSPOSE!
478  return cof[comp_x][comp_zeta] / det;
479  }
480  else
481  {
482  ERROR_MSG ( "Dimension (M_nbLocalCoor): only 1, 2 or 3!" );
483  };
484  return 0.;
485 }
486 
487 //----------------------------------------------------------------------
488 
489 void CurrentFE::update (const std::vector<std::vector<Real> >& pts, flag_Type upFlag)
490 {
491  ASSERT (M_nbQuadPt != 0 || upFlag == UPDATE_ONLY_CELL_NODES, " No quadrature rule defined, cannot update!");
492 
493  M_cellNodesUpdated = false;
494  if ( (upFlag & UPDATE_ONLY_CELL_NODES) != 0)
495  {
496  computeCellNodes (pts);
497  };
498 
499  M_quadNodesUpdated = false;
500  if ( (upFlag & UPDATE_ONLY_QUAD_NODES) != 0)
501  {
502  computeQuadNodes();
503  };
504 
505  M_jacobianUpdated = false;
506  if ( (upFlag & UPDATE_ONLY_JACOBIAN) != 0)
507  {
508  computeJacobian();
509  };
510 
511  M_tInverseJacobianUpdated = false;
512  if ( (upFlag & UPDATE_ONLY_T_INVERSE_JACOBIAN) != 0)
513  {
514  computeTInverseJacobian();
515  };
516 
517  M_detJacobianUpdated = false;
518  if ( (upFlag & UPDATE_ONLY_DET_JACOBIAN) != 0)
519  {
520  computeDetJacobian();
521  };
522 
523  M_wDetJacobianUpdated = false;
524  if ( (upFlag & UPDATE_ONLY_W_DET_JACOBIAN) != 0)
525  {
526  computeWDetJacobian();
527  };
528 
529  M_dphiUpdated = false;
530  if ( (upFlag & UPDATE_ONLY_DPHI) != 0)
531  {
532  computeDphi();
533  };
534 
535  M_d2phiUpdated = false;
536  if ( (upFlag & UPDATE_ONLY_D2PHI) != 0)
537  {
538  computeD2phi();
539  };
540 
541  M_phiVectUpdated = false;
542  if ( (upFlag & UPDATE_ONLY_PHI_VECT) != 0)
543  {
544  computePhiVect();
545  };
546 }
547 
548 void CurrentFE::update (const std::vector<GeoVector>& pts, flag_Type upFlag)
549 {
550  std::vector< std::vector <Real > > newpts (pts.size(), std::vector<Real> (pts[0].size() ) );
551 
552  for (UInt iPt (0); iPt < pts.size(); ++iPt)
553  {
554  for (UInt iCoord (0); iCoord < pts[iPt].size(); ++iCoord)
555  {
556  newpts[iPt][iCoord] = pts[iPt][iCoord];
557  }
558  }
559 
560  update (newpts, upFlag);
561 }
562 
563 Real CurrentFE::measure() const
564 {
565  ASSERT (M_wDetJacobianUpdated, "Weighted jacobian determinant is not updated!");
566 
567  Real meas (0.0);
568  for ( UInt iterQuadNode (0); iterQuadNode < M_nbQuadPt; ++iterQuadNode )
569  {
570  meas += M_wDetJacobian[iterQuadNode];
571  }
572  return meas;
573 }
574 
575 void CurrentFE::barycenter ( Real& x, Real& y, Real& z ) const
576 {
577  ASSERT (M_cellNodesUpdated, "Cell nodes are not updated!");
578 
579  x = 0.0;
580  y = 0.0;
581  z = 0.0;
582 
583  for (UInt iNode (0); iNode < M_nbGeoNode; ++iNode)
584  {
585  x += M_cellNodes[iNode][0];
586  y += M_cellNodes[iNode][1];
587  z += M_cellNodes[iNode][2];
588  }
589 
590  x /= M_nbGeoNode;
591  y /= M_nbGeoNode;
592  z /= M_nbGeoNode;
593 }
594 
595 Real CurrentFE::diameter() const
596 {
597  ASSERT (M_cellNodesUpdated, "Cell nodes are not updated!");
598 
599  Real s (0.0);
600  Real h (0.0);
601  for ( UInt i (0); i < M_nbGeoNode - 1; ++i)
602  {
603  for ( UInt j (i + 1); j < M_nbGeoNode; ++j )
604  {
605  s = 0.0;
606  for ( Int icoor (0); icoor < (Int) nDimensions; ++icoor )
607  {
608  s += std::fabs ( M_cellNodes[i][icoor] - M_cellNodes[j][icoor] );
609  }
610  if ( s > h )
611  {
612  h = s;
613  }
614  }
615  }
616 
617  return h;
618 }
619 
620 Real CurrentFE::diameter2() const
621 {
622  ASSERT (M_cellNodesUpdated, "Cell nodes are not updated!");
623 
624  Real s (0.0);
625  Real d (0.0);
626  Real h (0.0);
627  for ( UInt i (0); i < M_nbGeoNode - 1; ++i)
628  {
629  for ( UInt j (i + 1); j < M_nbGeoNode; ++j)
630  {
631  s = 0.;
632  for ( Int icoor (0); icoor < (Int) nDimensions; ++icoor)
633  {
634  d = ( M_cellNodes[i][icoor] - M_cellNodes[j][icoor] );
635  d = d * d;
636  s += d;
637  }
638  s = std::sqrt (s);
639  if ( s > h )
640  {
641  h = s;
642  }
643  }
644  }
645  return h;
646 }
647 
648 void CurrentFE::coorMap ( Real& x, Real& y, Real& z,
649  Real xi, Real eta, Real zeta ) const
650 {
651  ASSERT (M_cellNodesUpdated, "Cell nodes are not updated!");
652 
653  x = 0.0;
654  y = 0.0;
655  z = 0.0;
656 
657  // Note: we assume (already in the constructor) that the 2nd dimension of
658  // M_cellNodes is nDimensions, that is 3. If in the future nDimensions!=3
659  // this piece of code needs to be modified.
660 
661  Real geoMap;
662  for (UInt iNode (0); iNode < M_nbGeoNode; ++iNode)
663  {
664  geoMap = M_geoMap->phi ( iNode, xi, eta, zeta );
665  x += M_cellNodes[iNode][0] * geoMap;
666  y += M_cellNodes[iNode][1] * geoMap;
667  z += M_cellNodes[iNode][2] * geoMap;
668  }
669 }
670 
671 GeoVector CurrentFE::coorMap (const GeoVector& P) const
672 {
673  ASSERT (M_cellNodesUpdated, "Cell nodes are not updated!");
674 
675  GeoVector Pcurrent ( nDimensions, 0. );
676 
677  for ( UInt i (0); i < M_nbGeoNode; ++i )
678  {
679  for (UInt j (0); j < nDimensions; ++j)
680  {
681  Pcurrent[j] += M_cellNodes[i][j] * M_geoMap->phi (i, P);
682  }
683  }
684 
685  return Pcurrent;
686 }
687 
688 
689 void CurrentFE::quadRuleVTKexport ( const std::string& filename) const
690 {
691  ASSERT (M_quadNodesUpdated, "Quad nodes are not updated! No export possible");
692 
693  std::ofstream output (filename.c_str() );
694  ASSERT (!output.fail(), " Unable to open the file for the export of the quadrature ");
695 
696  // Header
697  output << "# vtk DataFile Version 3.0" << std::endl;
698  output << "LifeV : Quadrature export" << std::endl;
699  output << "ASCII" << std::endl;
700  output << "DATASET POLYDATA" << std::endl;
701  output << "POINTS " << M_nbQuadPt << " float" << std::endl;
702 
703  for (UInt i (0); i < M_nbQuadPt; ++i)
704  {
705  output << M_quadNodes[i][0] << " " << M_quadNodes[i][1] << " " << M_quadNodes[i][2] << std::endl;
706  };
707 
708  output << "VERTICES " << M_nbQuadPt << " " << 2 * M_nbQuadPt << std::endl;
709 
710  for (UInt i (0); i < M_nbQuadPt; ++i)
711  {
712  output << 1 << " " << i << std::endl;
713  };
714 
715  output.close();
716 }
717 
718 
719 void CurrentFE::setQuadRule (const QuadratureRule& newQuadRule)
720 {
721  // Set the quadrature
722  if (M_quadRule != 0)
723  {
724  delete M_quadRule;
725  };
726  M_quadRule = new QuadratureRule (newQuadRule);
727 
728  // Adjust the constants related to the quadrature
729  M_nbQuadPt = UInt ( newQuadRule.nbQuadPt() );
730 
731  // Resize all the arrays that need it
732  M_quadNodes.resize (boost::extents[M_nbQuadPt][3]);
733 
734  M_dphiGeometricMap.resize (boost::extents[M_nbGeoNode][M_nbLocalCoor][M_nbQuadPt]);
735  M_jacobian.resize (boost::extents[M_nbLocalCoor][M_nbLocalCoor][M_nbQuadPt]);
736  M_detJacobian.resize (boost::extents[M_nbQuadPt]);
737  M_wDetJacobian.resize (boost::extents[M_nbQuadPt]);
738  M_tInverseJacobian.resize (boost::extents[M_nbLocalCoor][M_nbLocalCoor][M_nbQuadPt]);
739 
740  M_phi.resize (boost::extents[M_nbNode][M_refFE->feDim()][M_nbQuadPt]);
741  M_dphi.resize (boost::extents[M_nbNode][M_nbLocalCoor][M_nbQuadPt]);
742  M_d2phi.resize (boost::extents[M_nbNode][M_nbLocalCoor][M_nbLocalCoor][M_nbQuadPt]);
743  M_phiVect.resize (boost::extents[M_nbNode][M_refFE->feDim()][M_nbQuadPt]);
744 
745  M_dphiRef.resize (boost::extents[M_nbNode][M_nbLocalCoor][M_nbQuadPt]);
746  M_d2phiRef.resize (boost::extents[M_nbNode][M_nbLocalCoor][M_nbLocalCoor][M_nbQuadPt]);
747  M_divPhiRef.resize (boost::extents[M_nbNode][M_nbQuadPt]);
748 
749  // All the updates have to be done again
750  M_cellNodesUpdated = false;
751  M_quadNodesUpdated = false;
752 
753  M_dphiGeometricMapUpdated = false;
754  M_jacobianUpdated = false;
755  M_detJacobianUpdated = false;
756  M_wDetJacobianUpdated = false;
757  M_tInverseJacobianUpdated = false;
758 
759  M_phiUpdated = false;
760  M_dphiUpdated = false;
761  M_d2phiUpdated = false;
762  M_phiVectUpdated = false;
763 
764  M_dphiRefUpdated = false;
765  M_d2phiRefUpdated = false;
766  M_divPhiRefUpdated = false;
767 
768  // Update what can be updated
769  for ( UInt iterQuad (0); iterQuad < M_nbQuadPt; ++iterQuad )
770  {
771  for ( UInt iterNode (0); iterNode < M_nbNode; ++iterNode )
772  {
773  // --- PHI ---
774  if (M_refFE->hasPhi() )
775  {
776  for (UInt iterFEDim (0); iterFEDim < M_refFE->feDim(); ++iterFEDim)
777  {
778  M_phi[iterNode][iterFEDim][iterQuad] = M_refFE->phi ( iterNode, iterFEDim,
779  M_quadRule->quadPointCoor (iterQuad) );
780  }
781  }
782 
783  // --- DIV ---
784  if (M_refFE->hasDivPhi() )
785  {
786  M_divPhiRef[iterNode][iterQuad] = M_refFE->divPhi (iterNode,
787  M_quadRule->quadPointCoor (iterQuad) );
788  }
789 
790  for ( UInt icoor (0); icoor < M_nbLocalCoor; ++icoor )
791  {
792  // --- DPHI ---
793  if (M_refFE->hasDPhi() )
794  {
795  M_dphiRef[iterNode][icoor][iterQuad] = M_refFE->dPhi ( iterNode,
796  icoor,
797  M_quadRule->quadPointCoor (iterQuad) );
798  }
799 
800  // --- D2PHI ---
801  if (M_refFE->hasD2Phi() )
802  {
803  for ( UInt jcoor (0); jcoor < M_nbLocalCoor; ++jcoor )
804  {
805  M_d2phiRef[iterNode][icoor][jcoor][iterQuad] = M_refFE->d2Phi ( iterNode,
806  icoor, jcoor,
807  M_quadRule->quadPointCoor (iterQuad) );
808  }
809  }
810  }
811  }
812  for ( UInt k (0); k < M_nbGeoNode; ++k )
813  {
814  for ( UInt icoor (0); icoor < M_nbLocalCoor; ++icoor )
815  {
816  M_dphiGeometricMap[k][icoor][iterQuad] = M_geoMap->dPhi ( k, icoor,
817  M_quadRule->quadPointCoor (iterQuad) );
818  }
819  }
820  }
821 
822  M_phiUpdated = true;
823  M_dphiRefUpdated = true;
824  M_divPhiRefUpdated = true;
825  M_d2phiRefUpdated = true;
826  M_dphiGeometricMapUpdated = true;
827 }
828 
829 void CurrentFE::computeCellNodes ( const std::vector< std::vector<Real> >& pts)
830 {
831  for (UInt iterNode (0); iterNode < M_nbGeoNode; ++iterNode)
832  for (UInt iterCoord (0); iterCoord < nDimensions; ++iterCoord)
833  {
834  M_cellNodes[iterNode][iterCoord] = pts[iterNode][iterCoord];
835  }
836 
837  M_cellNodesUpdated = true;
838 }
839 
840 void CurrentFE::computeQuadNodes()
841 {
842  ASSERT (M_cellNodesUpdated, "Missing update: cellNodes");
843 
844  for ( UInt iterQuadNode (0); iterQuadNode < M_nbQuadPt; ++iterQuadNode)
845  {
846  GeoVector quadNode (coorMap (M_quadRule->quadPointCoor (iterQuadNode) ) );
847  for (UInt i = 0; i < quadNode.size(); i++)
848  {
849  M_quadNodes[iterQuadNode][i] = quadNode[i];
850  }
851  }
852  M_quadNodesUpdated = true;
853 }
854 
855 void CurrentFE::computeDphiGeometricMap()
856 {
857  for (UInt iterQuad (0); iterQuad < M_nbQuadPt ; ++iterQuad)
858  {
859  for (UInt iterNode (0); iterNode < M_nbGeoNode; ++iterNode)
860  {
861  for (UInt iterCoor (0); iterCoor < M_nbLocalCoor; ++iterCoor)
862  {
863  M_dphiGeometricMap[iterNode][iterCoor][iterQuad] = M_geoMap->dPhi (iterNode, iterCoor,
864  M_quadRule->quadPointCoor (iterQuad) );
865  }
866  }
867  }
868  M_dphiGeometricMapUpdated = true;
869 }
870 
871 void CurrentFE::computeJacobian()
872 {
873  ASSERT (M_cellNodesUpdated, "Missing update: cellNodes");
874  ASSERT (M_dphiGeometricMapUpdated, "Missing update: dphiGeometricMap");
875 
876  Real partialSum;
877 
878  for (UInt iterQuad (0); iterQuad < M_nbQuadPt ; ++iterQuad)
879  {
880  for (UInt icoord (0); icoord < M_nbLocalCoor; ++icoord)
881  {
882  for (UInt jcoord (0); jcoord < M_nbLocalCoor; ++jcoord)
883  {
884  partialSum = 0;
885  for (UInt iterNode (0); iterNode < M_nbGeoNode; ++iterNode)
886  {
887  partialSum += M_cellNodes[iterNode][icoord] * M_dphiGeometricMap[iterNode][jcoord][iterQuad];
888  }
889  M_jacobian[icoord][jcoord][iterQuad] = partialSum;
890  }
891  }
892  }
893  M_jacobianUpdated = true;
894 }
895 
896 void CurrentFE::computeTInverseJacobian()
897 {
898  ASSERT (M_jacobianUpdated, "Missing update: jacobian");
899 
900  for (UInt iterQuad (0); iterQuad < M_nbQuadPt ; ++iterQuad)
901  {
902  switch (M_nbLocalCoor)
903  {
904  case 1:
905  {
906  Real a ( M_jacobian[0][0][iterQuad] );
907 
908  M_tInverseJacobian[0][0][iterQuad] = 1.0 / a ;
909 
910  break;
911  }
912  case 2:
913  {
914  Real a ( M_jacobian[0][0][iterQuad] );
915  Real b ( M_jacobian[0][1][iterQuad] );
916  Real c ( M_jacobian[1][0][iterQuad] );
917  Real d ( M_jacobian[1][1][iterQuad] );
918 
919  Real det ( a * d - b * c );
920 
921  M_tInverseJacobian[0][0][iterQuad] = d / det ;
922  M_tInverseJacobian[0][1][iterQuad] = -c / det ;
923 
924  M_tInverseJacobian[1][0][iterQuad] = -b / det ;
925  M_tInverseJacobian[1][1][iterQuad] = a / det ;
926 
927  break;
928  }
929  case 3:
930  {
931  Real a ( M_jacobian[0][0][iterQuad] );
932  Real b ( M_jacobian[0][1][iterQuad] );
933  Real c ( M_jacobian[0][2][iterQuad] );
934  Real d ( M_jacobian[1][0][iterQuad] );
935  Real e ( M_jacobian[1][1][iterQuad] );
936  Real f ( M_jacobian[1][2][iterQuad] );
937  Real g ( M_jacobian[2][0][iterQuad] );
938  Real h ( M_jacobian[2][1][iterQuad] );
939  Real i ( M_jacobian[2][2][iterQuad] );
940 
941  Real ei ( e * i );
942  Real fh ( f * h );
943  Real bi ( b * i );
944  Real ch ( c * h );
945  Real bf ( b * f );
946  Real ce ( c * e );
947  Real det ( a * ( ei - fh ) + d * ( ch - bi ) + g * ( bf - ce ) );
948 
949  M_tInverseJacobian[0][0][iterQuad] = ( ei - fh ) / det ;
950  M_tInverseJacobian[0][1][iterQuad] = ( -d * i + f * g ) / det ;
951  M_tInverseJacobian[0][2][iterQuad] = ( d * h - e * g ) / det ;
952 
953  M_tInverseJacobian[1][0][iterQuad] = ( -bi + ch ) / det ;
954  M_tInverseJacobian[1][1][iterQuad] = ( a * i - c * g ) / det ;
955  M_tInverseJacobian[1][2][iterQuad] = ( -a * h + b * g ) / det ;
956 
957  M_tInverseJacobian[2][0][iterQuad] = ( bf - ce ) / det ;
958  M_tInverseJacobian[2][1][iterQuad] = ( -a * f + c * d ) / det ;
959  M_tInverseJacobian[2][2][iterQuad] = ( a * e - b * d ) / det ;
960 
961  break;
962  }
963  default:
964  ERROR_MSG ( "Dimension (M_nbLocalCoor): only 1, 2 or 3!" );
965  break;
966  }
967 
968  }
969 
970  M_tInverseJacobianUpdated = true;
971 }
972 
973 void CurrentFE::computeDetJacobian()
974 {
975  ASSERT (M_jacobianUpdated, "Missing update: jacobian");
976 
977  for (UInt iterQuad (0); iterQuad < M_nbQuadPt ; ++iterQuad)
978  {
979  switch (M_nbLocalCoor)
980  {
981  case 1:
982  {
983  Real a ( M_jacobian[0][0][iterQuad] );
984  Real det ( a);
985  M_detJacobian[iterQuad] = det;
986 
987  break;
988  }
989  case 2:
990  {
991  Real a ( M_jacobian[0][0][iterQuad] );
992  Real b ( M_jacobian[0][1][iterQuad] );
993  Real c ( M_jacobian[1][0][iterQuad] );
994  Real d ( M_jacobian[1][1][iterQuad] );
995 
996  Real det ( a * d - b * c);
997  M_detJacobian[iterQuad] = det;
998  break;
999  }
1000  case 3:
1001  {
1002  Real a ( M_jacobian[0][0][iterQuad] );
1003  Real b ( M_jacobian[0][1][iterQuad] );
1004  Real c ( M_jacobian[0][2][iterQuad] );
1005  Real d ( M_jacobian[1][0][iterQuad] );
1006  Real e ( M_jacobian[1][1][iterQuad] );
1007  Real f ( M_jacobian[1][2][iterQuad] );
1008  Real g ( M_jacobian[2][0][iterQuad] );
1009  Real h ( M_jacobian[2][1][iterQuad] );
1010  Real i ( M_jacobian[2][2][iterQuad] );
1011 
1012  Real ei (e * i);
1013  Real fh (f * h);
1014  Real bi (b * i);
1015  Real ch (c * h);
1016  Real bf (b * f);
1017  Real ce (c * e);
1018 
1019  Real det ( a * (ei - fh) + d * (ch - bi) + g * ( bf - ce) );
1020  M_detJacobian[iterQuad] = det;
1021  break;
1022  }
1023  default:
1024  ERROR_MSG ( "Dimension (M_nbLocalCoor): only 1, 2 or 3!" );
1025  break;
1026  }
1027  }
1028  M_detJacobianUpdated = true;
1029 }
1030 
1031 void CurrentFE::computeWDetJacobian()
1032 {
1033  ASSERT (M_detJacobianUpdated, "Missing update: determinant of the jacobian");
1034 
1035  for (UInt iterQuad (0); iterQuad < M_nbQuadPt ; ++iterQuad)
1036  {
1037  M_wDetJacobian[iterQuad] = M_detJacobian[iterQuad] * M_quadRule->weight (iterQuad);
1038  }
1039  M_wDetJacobianUpdated = true;
1040 }
1041 
1042 
1043 void CurrentFE::computeDphiRef()
1044 {
1045  for (UInt iterQuad (0); iterQuad < M_nbQuadPt ; ++iterQuad)
1046  {
1047  for (UInt iterNode (0); iterNode < M_nbNode; ++iterNode)
1048  {
1049  for (UInt iterCoor (0); iterCoor < M_nbLocalCoor; ++iterCoor)
1050  {
1051  M_dphiRef[iterNode][iterCoor][iterQuad] = M_refFE->dPhi (iterNode, iterCoor,
1052  M_quadRule->quadPointCoor (iterQuad) );
1053  }
1054  }
1055  }
1056  M_dphiRefUpdated = true;
1057 }
1058 
1059 void CurrentFE::computeDphi()
1060 {
1061  ASSERT (M_jacobianUpdated, "Missing update: tInverseJacobian");
1062  ASSERT (M_dphiRefUpdated, "Missing update: dphiRef");
1063 
1064  Real partialSum (0);
1065 
1066  for (UInt iterQuad (0); iterQuad < M_nbQuadPt ; ++iterQuad)
1067  {
1068  for (UInt iterNode (0); iterNode < M_nbNode; ++iterNode)
1069  {
1070  for (UInt iterCoor (0); iterCoor < M_nbLocalCoor; ++iterCoor)
1071  {
1072  partialSum = 0;
1073  for (UInt iter (0); iter < M_nbLocalCoor; ++iter)
1074  {
1075  partialSum += M_tInverseJacobian[iterCoor][iter][iterQuad] * M_dphiRef[iterNode][iter][iterQuad];
1076  }
1077  M_dphi[iterNode][iterCoor][iterQuad] = partialSum;
1078  }
1079  }
1080  }
1081  M_dphiUpdated = true;
1082 }
1083 
1084 
1085 void CurrentFE::computeD2phiRef()
1086 {
1087  for (UInt iterQuad (0); iterQuad < M_nbQuadPt ; ++iterQuad)
1088  {
1089  for (UInt iterNode (0); iterNode < M_nbNode; ++iterNode)
1090  {
1091  for (UInt iterCoor (0); iterCoor < M_nbLocalCoor; ++iterCoor)
1092  {
1093  for (UInt iterCoor2 (0); iterCoor2 < M_nbLocalCoor; ++iterCoor2)
1094  {
1095  M_d2phiRef[iterNode][iterCoor][iterCoor2][iterQuad] = M_refFE->d2Phi (iterNode, iterCoor, iterCoor2,
1096  M_quadRule->quadPointCoor (iterQuad) );
1097  }
1098  }
1099  }
1100  }
1101  M_d2phiRefUpdated = true;
1102 }
1103 
1104 
1105 void CurrentFE::computeD2phi()
1106 {
1107  ASSERT (M_jacobianUpdated, "Missing update: tInverseJacobian");
1108  ASSERT (M_d2phiRefUpdated, "Missing update: d2phiRef");
1109 
1110  Real partialSum (0);
1111  for ( UInt iterQuad (0); iterQuad < M_nbQuadPt ; ++iterQuad )
1112  {
1113  for ( UInt iterNode (0); iterNode < M_nbNode; ++iterNode )
1114  {
1115  for ( UInt icoor (0); icoor < M_nbLocalCoor; ++icoor )
1116  {
1117  for ( UInt jcoor (0); jcoor < M_nbLocalCoor; ++jcoor )
1118  {
1119  partialSum = 0.;
1120  for ( UInt k1 (0); k1 < M_nbLocalCoor; ++k1 )
1121  {
1122  for ( UInt k2 (0) ; k2 < M_nbLocalCoor; ++k2 )
1123  {
1124  partialSum += M_tInverseJacobian[icoor][k1][iterQuad]
1125  * M_d2phiRef[iterNode][k1][k2][iterQuad]
1126  * M_tInverseJacobian[jcoor][k2][iterQuad];
1127  }
1128  }
1129  M_d2phi[iterNode][icoor][jcoor][iterQuad] = partialSum;
1130  }
1131  }
1132  }
1133  }
1134  M_d2phiUpdated = true;
1135 }
1136 
1137 
1138 void CurrentFE::computePhiVect()
1139 {
1140  ASSERT (M_detJacobianUpdated, "Missing update: detJacobian");
1141 
1142  Real sum (0.);
1143  for ( UInt ig (0); ig < M_nbQuadPt; ++ig )
1144  {
1145  for ( UInt idof (0); idof < M_nbNode; ++idof )
1146  {
1147  for ( UInt icoor (0); icoor < M_nbLocalCoor; ++icoor )
1148  {
1149  sum = 0.;
1150  for ( UInt jcoor (0); jcoor < M_nbLocalCoor; ++jcoor )
1151  {
1152  sum += M_jacobian[ icoor ][ jcoor ][ ig ] * M_phi[ idof ][ jcoor ][ ig ];
1153  }
1154  M_phiVect[ idof ][ icoor ][ ig ] = sum / M_detJacobian[ ig ];
1155  }
1156  }
1157  }
1158 
1159  M_phiVectUpdated = true;
1160 }
1161 
1162 } // Namespace LifeV
LifeV::CurrentFE::computeDphiRef
void computeDphiRef()
Compute the value of the derivatives in the reference frame.
Definition: CurrentFE.cpp:1043
LifeV::GeometricMap
GeometricMap - Structure for the geometrical mapping.
Definition: GeometricMap.hpp:59
LifeV::CurrentFE::computeDetJacobian
void computeDetJacobian()
Compute the determinant of the jacobian.
Definition: CurrentFE.cpp:973
LifeV::CurrentFE::computeWDetJacobian
void computeWDetJacobian()
Compute the determinant of the jacobian in the quadrature nodes.
Definition: CurrentFE.cpp:1031
LifeV::CurrentFE::pointInverseJacobian
Real pointInverseJacobian(Real hat_x, Real hat_y, Real hat_z, int compx, int compzeta) const
Definition: CurrentFE.cpp:424
LifeV::CurrentFE::M_nbNode
const UInt M_nbNode
Definition: CurrentFE.hpp:611
LifeV::UPDATE_ONLY_PHI_VECT
const flag_Type UPDATE_ONLY_PHI_VECT(1024)
LifeV::UPDATE_ONLY_W_DET_JACOBIAN
const flag_Type UPDATE_ONLY_W_DET_JACOBIAN(32)
LifeV::flag_Type
uint32_type flag_Type
bit-flag with up to 32 different flags
Definition: LifeV.hpp:197
LifeV::CurrentFE::M_quadRule
QuadratureRule * M_quadRule
Definition: CurrentFE.hpp:627
LifeV::CurrentFE::M_wDetJacobianUpdated
bool M_wDetJacobianUpdated
Definition: CurrentFE.hpp:660
LifeV::UPDATE_ONLY_JACOBIAN
const flag_Type UPDATE_ONLY_JACOBIAN(8)
LifeV::UPDATE_ONLY_DPHI
const flag_Type UPDATE_ONLY_DPHI(128)
LifeV::CurrentFE::diameter2
Real diameter2() const
return the diameter of the element in the 2-norm
Definition: CurrentFE.cpp:620
LifeV::CurrentFE::computePhiVect
void computePhiVect()
Compute the value of the vectorial FE using Piola transform.
Definition: CurrentFE.cpp:1138
LifeV::CurrentFE::computeCellNodes
void computeCellNodes(const std::vector< std::vector< Real > > &pts)
Update only the nodes of the cells to the current one.
Definition: CurrentFE.cpp:829
LifeV::CurrentFE::coorBackMap
void coorBackMap(Real x, Real y, Real z, Real &xi, Real &eta, Real &zeta) const
Definition: CurrentFE.cpp:244
LifeV::CurrentFE::M_refFE
const ReferenceFE * M_refFE
Definition: CurrentFE.hpp:625
LifeV::CurrentFE::CurrentFE
CurrentFE(const CurrentFE &)
Definition: CurrentFE.cpp:193
LifeV::CurrentFE::computeDphiGeometricMap
void computeDphiGeometricMap()
Compute the values of the derivatives of the mapping in the quadrature nodes.
Definition: CurrentFE.cpp:855
LifeV::Int
int32_type Int
Generic integer data.
Definition: LifeV.hpp:188
LifeV::CurrentFE::computeQuadNodes
void computeQuadNodes()
Update the location of the quadrature in the current cell.
Definition: CurrentFE.cpp:840
LifeV::GeoVector
boost::numeric::ublas::vector< Real > GeoVector
Definition: QuadraturePoint.hpp:52
LifeV::CurrentFE::M_geoMap
const GeometricMap * M_geoMap
Definition: CurrentFE.hpp:626
LifeV::CurrentFE::coorMap
GeoVector coorMap(const GeoVector &P) const
Definition: CurrentFE.cpp:671
LifeV::QuadratureRule::QuadratureRule
QuadratureRule(const QuadratureRule &qr)
Copy constructor.
Definition: QuadratureRule.cpp:68
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::CurrentFE::update
void update(const std::vector< GeoVector > &pts, flag_Type upFlag)
Update method using only point coordinates. It used the flags, as defined in this page...
Definition: CurrentFE.cpp:548
LifeV::CurrentFE::computeD2phiRef
void computeD2phiRef()
Compute the value of the second derivatives in the reference frame.
Definition: CurrentFE.cpp:1085
LifeV::CurrentFE::M_nbLocalCoor
const UInt M_nbLocalCoor
Definition: CurrentFE.hpp:612
LifeV::CurrentFE::computeD2phi
void computeD2phi()
Compute the value of the derivatives in the current element.
Definition: CurrentFE.cpp:1105
ERROR_MSG
#define ERROR_MSG(A)
Definition: LifeAssert.hpp:69
LifeV::CurrentFE::M_nbPattern
const UInt M_nbPattern
Definition: CurrentFE.hpp:615
LifeV::UPDATE_ONLY_DET_JACOBIAN
const flag_Type UPDATE_ONLY_DET_JACOBIAN(4096)
LifeV::UPDATE_ONLY_QUAD_NODES
const flag_Type UPDATE_ONLY_QUAD_NODES(2)
LifeV::QuadratureRule::quadPointCoor
const GeoVector & quadPointCoor(const UInt &ig) const
quadPointCoor(ig) is the full coordinates of the quadrature point ig
Definition: QuadratureRule.hpp:298
ASSERT
#define ASSERT(X, A)
Definition: LifeAssert.hpp:90
LifeV::CurrentFE::M_detJacobianUpdated
bool M_detJacobianUpdated
Definition: CurrentFE.hpp:659
LifeV::CurrentFE::M_tInverseJacobianUpdated
bool M_tInverseJacobianUpdated
Definition: CurrentFE.hpp:661
LifeV::CurrentFE::computeTInverseJacobian
void computeTInverseJacobian()
Compute the transposed inverse of the jacobian in the quadrature nodes.
Definition: CurrentFE.cpp:896
LifeV::CurrentFE::M_nbUpper
const UInt M_nbUpper
Definition: CurrentFE.hpp:614
LifeV::CurrentFE::measure
virtual Real measure() const
Return the measure of the current element.
Definition: CurrentFE.cpp:563
LifeV::ReferenceElement::nbDof
const UInt & nbDof() const
Return the number of degrees of freedom for this reference element.
Definition: ReferenceElement.hpp:233
LifeV::CurrentFE::diameter
Real diameter() const
return the diameter of the element in the 1-norm
Definition: CurrentFE.cpp:595
LifeV::CurrentFE::computeDphi
void computeDphi()
Compute the value of the derivatives in the current element.
Definition: CurrentFE.cpp:1059
LifeV::UPDATE_ONLY_D2PHI
const flag_Type UPDATE_ONLY_D2PHI(512)
LifeV::CurrentFE::M_nbDiag
const UInt M_nbDiag
Definition: CurrentFE.hpp:613
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::computeJacobian
void computeJacobian()
Compute the jacobian in the quadrature nodes.
Definition: CurrentFE.cpp:871
LifeV::UPDATE_ONLY_T_INVERSE_JACOBIAN
const flag_Type UPDATE_ONLY_T_INVERSE_JACOBIAN(16)
LifeV::ReferenceFE
The class for a reference Lagrangian finite element.
Definition: ReferenceFE.hpp:115
LifeV::CurrentFE::barycenter
void barycenter(Real &x, Real &y, Real &z) const
return the barycenter of the element
Definition: CurrentFE.cpp:575
LifeV::nDimensions
const UInt nDimensions(NDIM)
LifeV::CurrentFE::CurrentFE
CurrentFE(const ReferenceFE &refFE, const GeometricMap &geoMap)
Constructor without quadrature specification.
Definition: CurrentFE.cpp:157
LifeV::CurrentFE::pointJacobian
Real pointJacobian(Real hat_x, Real hat_y, Real hat_z, int compx, int compzeta) const
Definition: CurrentFE.cpp:350
LifeV::CurrentFE::CurrentFE
CurrentFE(const ReferenceFE &refFE, const GeometricMap &geoMap, const QuadratureRule &qr)
Full constructor with default quadrature.
Definition: CurrentFE.cpp:43
LifeV::CurrentFE::setQuadRule
void setQuadRule(const QuadratureRule &newQuadRule)
Setter for the quadrature rule.
Definition: CurrentFE.cpp:719
LifeV::ReferenceElement::phi
Real phi(UInt i, const Real &x, const Real &y, const Real &z) const
return the value of the i-th basis function on point (x,y,z)
Definition: ReferenceElement.hpp:268
LifeV::QuadratureRule
QuadratureRule - The basis class for storing and accessing quadrature rules.
Definition: QuadratureRule.hpp:111
LifeV::CurrentFE::M_dphiGeometricMapUpdated
bool M_dphiGeometricMapUpdated
Definition: CurrentFE.hpp:657
LifeV::CurrentFE::pointDetJacobian
Real pointDetJacobian(Real hat_x, Real hat_y, Real hat_z) const
Definition: CurrentFE.cpp:368
LifeV::QuadratureRule::nbQuadPt
const UInt & nbQuadPt() const
Getter for the number of quadrature points.
Definition: QuadratureRule.hpp:305
LifeV::CurrentFE::quadRuleVTKexport
void quadRuleVTKexport(const std::string &filename) const
Export the quadrature rule on the current FE.
Definition: CurrentFE.cpp:689
LifeV::UInt
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191
LifeV::CurrentFE::coorMap
void coorMap(Real &x, Real &y, Real &z, Real xi, Real eta, Real zeta) const
Definition: CurrentFE.cpp:648
LifeV::UPDATE_ONLY_CELL_NODES
const flag_Type UPDATE_ONLY_CELL_NODES(1)