LifeV
EvaluateNodalExpressionVectorElement.hpp
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 the LifeV library
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
21  License along with this library; if not, see <http://www.gnu.org/licenses/>
22 
23 
24 *******************************************************************************
25 */
26 //@HEADER
27 
28 /*!
29  * @file
30  @brief This file contains the definition of the IntegrateVectorElement class.
31 
32  @date 06/2011
33  @author Samuel Quinodoz <samuel.quinodoz@epfl.ch>
34  */
35 
36 #ifndef EVALUATENODALEXPRESSION_VECTOR_ELEMENT_HPP
37 #define EVALUATENODALEXPRESSION_VECTOR_ELEMENT_HPP
38 
39 #include <lifev/core/LifeV.hpp>
40 
41 #include <lifev/core/fem/QuadratureRule.hpp>
42 #include <lifev/eta/fem/ETCurrentFE.hpp>
43 #include <lifev/eta/fem/MeshGeometricMap.hpp>
44 #include <lifev/eta/fem/QRAdapterBase.hpp>
45 
46 #include <lifev/eta/expression/ExpressionToEvaluation.hpp>
47 
48 #include <lifev/eta/array/ETVectorElemental.hpp>
49 
50 
51 
52 namespace LifeV
53 {
54 
55 namespace ExpressionAssembly
56 {
57 
58 //! The class to actually perform the loop over the elements to assemble a vector
59 /*!
60  @author Samuel Quinodoz <samuel.quinodoz@epfl.ch>
61 
62  This class is used to store the data required for the assembly of a vector and
63  perform that assembly with a loop over the elements, and then, for each elements,
64  using the Evaluation corresponding to the Expression (this convertion is done
65  within a typedef).
66  */
67 template < typename MeshType, typename SolutionSpaceType, typename ExpressionType, typename QRAdapterType>
68 class EvaluateNodalExpressionVectorElement
69 {
70 public:
71 
72  //! @name Public Types
73  //@{
74 
75  //! Type of the Evaluation
76  typedef typename ExpressionToEvaluation < ExpressionType,
77  SolutionSpaceType::field_dim,
78  0,
79  MeshType::S_geoDimensions >::evaluation_Type evaluation_Type;
80 
81  //@}
82 
83 
84  //! @name Constructors, destructor
85  //@{
86 
87  //! Full data constructor
88  EvaluateNodalExpressionVectorElement (const std::shared_ptr<MeshType>& mesh,
89  const QRAdapterType& qrAdapter,
90  const std::shared_ptr<SolutionSpaceType>& testSpace,
91  const ExpressionType& expression);
92 
93  //! Copy constructor
94  EvaluateNodalExpressionVectorElement ( const EvaluateNodalExpressionVectorElement < MeshType, SolutionSpaceType, ExpressionType, QRAdapterType>& integrator);
95 
96  //! Destructor
97  ~EvaluateNodalExpressionVectorElement();
98 
99  //@}
100 
101 
102  //! @name Operator
103  //@{
104 
105  //! Operator wrapping the addTo method
106  template <typename VectorType>
107  inline void operator>> (VectorType& vector)
108  {
109  addTo (vector);
110  }
111 
112  //! Operator wrapping the addTo method (for shared_ptr)
113  template <typename VectorType>
114  inline void operator>> (std::shared_ptr<VectorType> vector)
115  {
116  addTo (*vector);
117  }
118 
119 
120  //@}
121 
122  //! @name Methods
123  //@{
124 
125  //! Ouput method
126  void check (std::ostream& out = std::cout);
127 
128  //! Method that performs the assembly
129  /*!
130  The loop over the elements is located right
131  in this method. Everything for the assembly is then
132  performed: update the values, update the local vector,
133  sum over the quadrature nodes, assemble in the global
134  vector.
135  */
136  template <typename VectorType>
137  void addTo (VectorType& vec);
138 
139  //@}
140 
141 private:
142 
143  //! @name Private Methods
144  //@{
145 
146  // No default constructor
147  EvaluateNodalExpressionVectorElement();
148 
149  //@}
150 
151  // Pointer on the mesh
152  std::shared_ptr<MeshType> M_mesh;
153 
154  // Quadrature to be used
155  QRAdapterType M_qrAdapter;
156 
157  // Shared pointer on the Space
158  std::shared_ptr<SolutionSpaceType> M_solutionSpace;
159 
160  // Tree to compute the values for the assembly
161  evaluation_Type M_evaluation;
162 
163  ETCurrentFE<MeshType::S_geoDimensions, 1>* M_globalCFE_std;
164  ETCurrentFE<MeshType::S_geoDimensions, 1>* M_globalCFE_adapted;
165 
166  ETCurrentFE<SolutionSpaceType::space_dim, SolutionSpaceType::field_dim>* M_testCFE_std;
167  ETCurrentFE<SolutionSpaceType::space_dim, SolutionSpaceType::field_dim>* M_testCFE_adapted;
168 
169  ETVectorElemental M_elementalVector;
170 };
171 
172 
173 // ===================================================
174 // IMPLEMENTATION
175 // ===================================================
176 
177 // ===================================================
178 // Constructors & Destructor
179 // ===================================================
180 
181 template < typename MeshType, typename SolutionSpaceType, typename ExpressionType, typename QRAdapterType>
182 EvaluateNodalExpressionVectorElement < MeshType, SolutionSpaceType, ExpressionType, QRAdapterType>::
183 EvaluateNodalExpressionVectorElement (const std::shared_ptr<MeshType>& mesh,
184  const QRAdapterType& qrAdapter,
185  const std::shared_ptr<SolutionSpaceType>& solutionSpace,
186  const ExpressionType& expression)
187  : M_mesh (mesh),
188  M_qrAdapter (qrAdapter),
189  M_solutionSpace (solutionSpace),
190  M_evaluation (expression),
191 
192  M_testCFE_std (new ETCurrentFE<SolutionSpaceType::space_dim, SolutionSpaceType::field_dim> (solutionSpace->refFE(), solutionSpace->geoMap(), qrAdapter.standardQR() ) ),
193  M_testCFE_adapted (new ETCurrentFE<SolutionSpaceType::space_dim, SolutionSpaceType::field_dim> (solutionSpace->refFE(), solutionSpace->geoMap(), qrAdapter.standardQR() ) ),
194 
195  M_elementalVector (SolutionSpaceType::field_dim * solutionSpace->refFE().nbDof() )
196 
197 {
198  switch (MeshType::geoShape_Type::BasRefSha::S_shape)
199  {
200  case LINE:
201  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feSegP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
202  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feSegP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
203  break;
204  case TRIANGLE:
205  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTriaP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR());
206  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTriaP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR());
207  break;
208  case QUAD:
209  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feQuadQ0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR());
210  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feQuadQ0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR());
211  break;
212  case TETRA:
213  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR());
214  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR());
215  break;
216  case HEXA:
217  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feHexaQ0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR());
218  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feHexaQ0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR());
219  break;
220  default:
221  ERROR_MSG ("Unrecognized element shape");
222  }
223 
224  M_evaluation.setQuadrature (qrAdapter.standardQR() );
225  M_evaluation.setGlobalCFE (M_globalCFE_std);
226  M_evaluation.setTestCFE (M_testCFE_std);
227 }
228 
229 
230 template < typename MeshType, typename SolutionSpaceType, typename ExpressionType, typename QRAdapterType>
231 EvaluateNodalExpressionVectorElement < MeshType, SolutionSpaceType, ExpressionType, QRAdapterType>::
232 EvaluateNodalExpressionVectorElement ( const EvaluateNodalExpressionVectorElement < MeshType, SolutionSpaceType, ExpressionType, QRAdapterType>& integrator)
233  : M_mesh (integrator.M_mesh),
234  M_qrAdapter (integrator.M_qrAdapter),
235  M_solutionSpace (integrator.M_solutionSpace),
236  M_evaluation (integrator.M_evaluation),
237 
238  M_testCFE_std (new ETCurrentFE<SolutionSpaceType::space_dim, SolutionSpaceType::field_dim> (M_solutionSpace->refFE(), M_solutionSpace->geoMap(), integrator.M_qrAdapter.standardQR() ) ),
239  M_testCFE_adapted (new ETCurrentFE<SolutionSpaceType::space_dim, SolutionSpaceType::field_dim> (M_solutionSpace->refFE(), M_solutionSpace->geoMap(), integrator.M_qrAdapter.standardQR() ) ),
240 
241  M_elementalVector (integrator.M_elementalVector)
242 {
243  switch (MeshType::geoShape_Type::BasRefSha::S_shape)
244  {
245  case LINE:
246  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feSegP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
247  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feSegP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
248  break;
249  case TRIANGLE:
250  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTriaP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR());
251  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTriaP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR());
252  break;
253  case QUAD:
254  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feQuadQ0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR());
255  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feQuadQ0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR());
256  break;
257  case TETRA:
258  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR());
259  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR());
260  break;
261  case HEXA:
262  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feHexaQ0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR());
263  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feHexaQ0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR());
264  break;
265  default:
266  ERROR_MSG ("Unrecognized element shape");
267  }
268  M_evaluation.setQuadrature (integrator.M_qrAdapter.standardQR() );
269  M_evaluation.setGlobalCFE (M_globalCFE_std);
270  M_evaluation.setTestCFE (M_testCFE_std);
271 }
272 
273 
274 template < typename MeshType, typename SolutionSpaceType, typename ExpressionType, typename QRAdapterType>
275 EvaluateNodalExpressionVectorElement < MeshType, SolutionSpaceType, ExpressionType, QRAdapterType>::
276 ~EvaluateNodalExpressionVectorElement()
277 {
278  delete M_globalCFE_std;
279  delete M_globalCFE_adapted;
280  delete M_testCFE_std;
281  delete M_testCFE_adapted;
282 }
283 
284 // ===================================================
285 // Methods
286 // ===================================================
287 
288 template < typename MeshType, typename SolutionSpaceType, typename ExpressionType, typename QRAdapterType>
289 void
290 EvaluateNodalExpressionVectorElement < MeshType, SolutionSpaceType, ExpressionType, QRAdapterType>::
291 check (std::ostream& out)
292 {
293  out << " Checking the integration : " << std::endl;
294  M_evaluation.display (out);
295  out << std::endl;
296  out << " Elemental vector : " << std::endl;
297  M_elementalVector.showMe (out);
298  out << std::endl;
299 }
300 
301 template < typename MeshType, typename SolutionSpaceType, typename ExpressionType, typename QRAdapterType>
302 template <typename VectorType>
303 void
304 EvaluateNodalExpressionVectorElement < MeshType, SolutionSpaceType, ExpressionType, QRAdapterType>::
305 addTo (VectorType& vec)
306 {
307  UInt nbElements (M_mesh->numElements() );
308  UInt nbQuadPt_std (M_qrAdapter.standardQR().nbQuadPt() );
309  UInt nbSolDof (M_solutionSpace->refFE().nbDof() );
310 
311  // Defaulted to true for security
312  bool isPreviousAdapted (true);
313 
314  for (UInt iElement (0); iElement < nbElements; ++iElement)
315  {
316  // Zeros out the elemental vector
317  M_elementalVector.zero();
318 
319  // Update the quadrature rule adapter
320  M_qrAdapter.update (iElement);
321 
322 
323  if (M_qrAdapter.isAdaptedElement() )
324  {
325  // Reset the quadrature in the different structures
326  M_evaluation.setQuadrature ( M_qrAdapter.adaptedQR() );
327  M_globalCFE_adapted -> setQuadratureRule ( M_qrAdapter.adaptedQR() );
328  M_testCFE_adapted -> setQuadratureRule ( M_qrAdapter.adaptedQR() );
329 
330  // Reset the CurrentFEs in the evaluation
331  M_evaluation.setGlobalCFE ( M_globalCFE_adapted );
332  M_evaluation.setTestCFE ( M_testCFE_adapted );
333 
334  // Update with the correct element
335  M_evaluation.update (iElement);
336 
337  // Update the currentFEs
338  M_globalCFE_adapted->update (M_mesh->element (iElement), evaluation_Type::S_globalUpdateFlag | ET_UPDATE_WDET);
339  M_testCFE_adapted->update (M_mesh->element (iElement), evaluation_Type::S_testUpdateFlag);
340 
341 
342  // Assembly
343  for (UInt iblock (0); iblock < SolutionSpaceType::field_dim; ++iblock)
344  {
345  // Set the row global indices in the local vector
346  for (UInt i (0); i < nbSolDof; ++i)
347  {
348  M_elementalVector.setRowIndex
349  (i + iblock * nbSolDof,
350  M_solutionSpace->dof().localToGlobalMap (iElement, i) + iblock * M_solutionSpace->dof().numTotalDof() );
351  }
352 
353  // Make the assembly
354  for (UInt iQuadPt (0); iQuadPt < M_qrAdapter.adaptedQR().nbQuadPt(); ++iQuadPt)
355  {
356  for (UInt i (0); i < nbSolDof; ++i)
357  {
358  M_elementalVector.element (i + iblock * nbSolDof) +=
359  M_evaluation.value_qi (iQuadPt, i + iblock * nbSolDof);
360 
361  }
362  }
363  }
364 
365  // Finally, set the flag
366  isPreviousAdapted = true;
367  }
368  else
369  {
370  // Check if the last one was adapted
371  if (isPreviousAdapted)
372  {
373  M_evaluation.setQuadrature ( M_qrAdapter.standardQR() );
374  M_evaluation.setGlobalCFE ( M_globalCFE_std );
375  M_evaluation.setTestCFE ( M_testCFE_std );
376 
377  isPreviousAdapted = false;
378  }
379 
380 
381  // Update the currentFEs
382  M_globalCFE_std->update (M_mesh->element (iElement), evaluation_Type::S_globalUpdateFlag | ET_UPDATE_WDET);
383  M_testCFE_std->update (M_mesh->element (iElement), evaluation_Type::S_testUpdateFlag );
384 
385  // Update the evaluation
386  M_evaluation.update (iElement);
387 
388  // Loop on the blocks
389  for (UInt iblock (0); iblock < SolutionSpaceType::field_dim; ++iblock)
390  {
391  // Set the row global indices in the local vector
392  for (UInt i (0); i < nbSolDof; ++i)
393  {
394  M_elementalVector.setRowIndex
395  (i + iblock * nbSolDof,
396  M_solutionSpace->dof().localToGlobalMap (iElement, i) + iblock * M_solutionSpace->dof().numTotalDof() );
397  }
398 
399  // Make the assembly
400  for (UInt iQuadPt (0); iQuadPt < nbQuadPt_std; ++iQuadPt)
401  {
402  for (UInt i (0); i < nbSolDof; ++i)
403  {
404  M_elementalVector.element (i + iblock * nbSolDof) +=
405  M_evaluation.value_qi (iQuadPt, i + iblock * nbSolDof);
406 
407  }
408  }
409  }
410 
411  }
412 
413  M_elementalVector.pushToGlobal (vec);
414  }
415 }
416 
417 
418 } // Namespace ExpressionAssembly
419 
420 } // Namespace LifeV
421 
422 #endif
LifeV::feHexaQ0
const ReferenceFEScalar feHexaQ0("Lagrange Q0 on a hexaedra", FE_Q0_3D, HEXA, 0, 0, 0, 1, 1, 3, fct_Q0_3D, derfct_Q0_3D, der2fct_Q0_3D, refcoor_Q0_3D, STANDARD_PATTERN, &feQuadQ0, &lagrangianTransform)
LifeV::feSegP0
const ReferenceFEScalar feSegP0("Lagrange P0 on a segment", FE_P0_1D, LINE, 0, 1, 0, 0, 1, 1, fct_P0_1D, derfct_P0_1D, der2fct_P0_1D, refcoor_P0_1D, STANDARD_PATTERN, &fePointP0, &lagrangianTransform)
LifeV::feTetraP0
const ReferenceFEScalar feTetraP0("Lagrange P0 on a tetraedra", FE_P0_3D, TETRA, 0, 0, 0, 1, 1, 3, fct_P0_3D, derfct_P0_3D, der2fct_P0_3D, refcoor_P0_3D, STANDARD_PATTERN, &feTriaP0, &lagrangianTransform)
LifeV::feTriaP0
const ReferenceFEScalar feTriaP0("Lagrange P0 on a triangle", FE_P0_2D, TRIANGLE, 0, 0, 1, 0, 1, 2, fct_P0_2D, derfct_P0_2D, der2fct_P0_2D, refcoor_P0_2D, STANDARD_PATTERN, &feSegP0, &lagrangianTransform)
ERROR_MSG
#define ERROR_MSG(A)
Definition: LifeAssert.hpp:69
LifeV::feQuadQ0
const ReferenceFEScalar feQuadQ0("Lagrange Q0 on a quadrangle", FE_Q0_2D, QUAD, 0, 0, 1, 0, 1, 2, fct_Q0_2D, derfct_Q0_2D, der2fct_Q0_2D, refcoor_Q0_2D, STANDARD_PATTERN, &feSegP0, &lagrangianTransform)
LifeV::UInt
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191