LifeV
IntegrateMatrixElementLSAdapted.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 IntegrateMatrixElementLSAdapted class.
31 
32  @date 06/2011
33  @author Samuel Quinodoz <samuel.quinodoz@epfl.ch>
34  */
35 
36 #ifndef INTEGRATE_MATRIX_ELEMENT_LS_ADAPTED_HPP
37 #define INTEGRATE_MATRIX_ELEMENT_LS_ADAPTED_HPP
38 
39 #include <lifev/core/LifeV.hpp>
40 
41 #include <lifev/level_set/fem/LevelSetQRAdapter.hpp>
42 #include <lifev/eta/fem/ETCurrentFE.hpp>
43 #include <lifev/eta/fem/MeshGeometricMap.hpp>
44 
45 #include <lifev/eta/expression/ExpressionToEvaluation.hpp>
46 
47 #include <lifev/eta/array/ETMatrixElemental.hpp>
48 
49 #include <boost/shared_ptr.hpp>
50 
51 
52 
53 namespace LifeV
54 {
55 
56 namespace ExpressionAssembly
57 {
58 
59 
60 //! The class to actually perform the loop over the elements to assemble a matrix
61 /*!
62  @author Samuel Quinodoz <samuel.quinodoz@epfl.ch>
63 
64  This class is used to store the data required for the assembly of a matrix and
65  perform that assembly with a loop over the elements, and then, for each elements,
66  using the Evaluation corresponding to the Expression (This convertion is done
67  within a typedef).
68 
69  The speciality of this class with respect to the LifeV::ExpressionAssembly::IntegrateMatrixElement class
70  is that the quadrature can be adapted.
71  */
72 template < typename MeshType,
73  typename TestSpaceType,
74  typename SolutionSpaceType,
75  typename ExpressionType,
76  typename LSFESpaceType,
77  typename LSVectorType >
78 class IntegrateMatrixElementLSAdapted
79 {
80 public:
81 
82  //! @name Public Types
83  //@{
84 
85  //! Type of the Evaluation
86  typedef typename ExpressionToEvaluation < ExpressionType,
87  TestSpaceType::S_fieldDim,
88  SolutionSpaceType::S_fieldDim,
89  3 >::evaluation_Type evaluation_Type;
90 
91  //! Type of the adapter
92  typedef LevelSetQRAdapter<LSFESpaceType, LSVectorType> QRAdapter_Type;
93 
94  //@}
95 
96 
97  //! @name Constructors, destructor
98  //@{
99 
100  //! Full data constructor
101  IntegrateMatrixElementLSAdapted (const std::shared_ptr<MeshType>& mesh,
102  const QRAdapter_Type& QRAdapter,
103  const std::shared_ptr<TestSpaceType>& testSpace,
104  const std::shared_ptr<SolutionSpaceType>& solutionSpace,
105  const ExpressionType& expression);
106 
107  //! Copy constructor
108  IntegrateMatrixElementLSAdapted (const IntegrateMatrixElementLSAdapted
109  < MeshType, TestSpaceType, SolutionSpaceType,
110  ExpressionType, LSFESpaceType, LSVectorType > & integrator);
111 
112  //! Destructor
113  ~IntegrateMatrixElementLSAdapted();
114 
115  //@}
116 
117 
118  //! @name Operators
119  //@{
120 
121  //! Operator wrapping the addTo method
122  template <typename MatrixType>
123  inline void operator>> (MatrixType& mat)
124  {
125  addTo (mat);
126  }
127 
128  template <typename MatrixType>
129  inline void operator>> (std::shared_ptr<MatrixType> mat)
130  {
131  addTo (*mat);
132  }
133 
134 
135  //@}
136 
137 
138  //! @name Methods
139  //@{
140 
141  //! Ouput method
142  void check (std::ostream& out = std::cout);
143 
144  //! Method that performs the assembly
145  /*!
146  The loop over the elements is located right
147  in this method. Everything for the assembly is then
148  performed: update the values, update the local matrix,
149  sum over the quadrature nodes, assemble in the global
150  matrix.
151  */
152  template <typename MatrixType>
153  void addTo (MatrixType& mat);
154 
155  //! Method that performs the assembly
156  /*!
157  The loop over the elements is located right
158  in this method. Everything for the assembly is then
159  performed: update the values, update the local matrix,
160  sum over the quadrature nodes, assemble in the global
161  matrix.
162 
163  Specialized for the case where the matrix is passed as a shared_ptr
164  */
165  template <typename MatrixType>
166  inline void addTo (std::shared_ptr<MatrixType> mat)
167  {
168  ASSERT (mat != 0, " Cannot assemble with an empty matrix");
169  addTo (*mat);
170  }
171 
172  //@}
173 
174 private:
175 
176  //! @name Private Methods
177  //@{
178 
179  //! No empty constructor
180  IntegrateMatrixElementLSAdapted();
181 
182  //@}
183 
184  // Pointer on the mesh
185  std::shared_ptr<MeshType> M_mesh;
186 
187  // Quadrature rule adapter to be used
188  QRAdapter_Type M_QRAdapter;
189 
190  // Shared pointer on the Spaces
191  std::shared_ptr<TestSpaceType> M_testSpace;
192  std::shared_ptr<SolutionSpaceType> M_solutionSpace;
193 
194  // Tree to compute the values for the assembly
195  evaluation_Type M_evaluation;
196 
197  // Duplicated CurrentFE
198  ETCurrentFE<3, 1>* M_globalCFE_unadapted;
199  ETCurrentFE<3, 1>* M_globalCFE_adapted;
200 
201  ETCurrentFE<3, TestSpaceType::S_fieldDim>* M_testCFE_unadapted;
202  ETCurrentFE<3, TestSpaceType::S_fieldDim>* M_testCFE_adapted;
203 
204  ETCurrentFE<3, SolutionSpaceType::S_fieldDim>* M_solutionCFE_unadapted;
205  ETCurrentFE<3, SolutionSpaceType::S_fieldDim>* M_solutionCFE_adapted;
206 
207  ETMatrixElemental M_elementalMatrix;
208 };
209 
210 
211 // ===================================================
212 // IMPLEMENTATION
213 // ===================================================
214 
215 // ===================================================
216 // Constructors & Destructor
217 // ===================================================
218 
219 template < typename MeshType,
220  typename TestSpaceType,
221  typename SolutionSpaceType,
222  typename ExpressionType,
223  typename LSFESpaceType,
224  typename LSVectorType >
225 IntegrateMatrixElementLSAdapted<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType, LSFESpaceType, LSVectorType>::
226 IntegrateMatrixElementLSAdapted (const std::shared_ptr<MeshType>& mesh,
227  const QRAdapter_Type& QRAdapter,
228  const std::shared_ptr<TestSpaceType>& testSpace,
229  const std::shared_ptr<SolutionSpaceType>& solutionSpace,
230  const ExpressionType& expression)
231  : M_mesh (mesh),
232  M_QRAdapter (QRAdapter),
233  M_testSpace (testSpace),
234  M_solutionSpace (solutionSpace),
235  M_evaluation (expression),
236 
237  M_globalCFE_unadapted (new ETCurrentFE<3, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), QRAdapter.standardQR() ) ),
238  M_globalCFE_adapted (new ETCurrentFE<3, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), QRAdapter.standardQR() ) ),
239 
240  M_testCFE_unadapted (new ETCurrentFE<3, TestSpaceType::S_fieldDim> (testSpace->refFE(), testSpace->geoMap(), QRAdapter.standardQR() ) ),
241  M_testCFE_adapted (new ETCurrentFE<3, TestSpaceType::S_fieldDim> (testSpace->refFE(), testSpace->geoMap(), QRAdapter.standardQR() ) ),
242 
243  M_solutionCFE_unadapted (new ETCurrentFE<3, SolutionSpaceType::S_fieldDim> (solutionSpace->refFE(), testSpace->geoMap(), QRAdapter.standardQR() ) ),
244  M_solutionCFE_adapted (new ETCurrentFE<3, SolutionSpaceType::S_fieldDim> (solutionSpace->refFE(), testSpace->geoMap(), QRAdapter.standardQR() ) ),
245 
246  M_elementalMatrix (TestSpaceType::S_fieldDim * testSpace->refFE().nbDof(),
247  SolutionSpaceType::S_fieldDim * solutionSpace->refFE().nbDof() )
248 {
249  M_evaluation.setQuadrature ( M_QRAdapter.standardQR() );
250  M_evaluation.setGlobalCFE ( M_globalCFE_unadapted );
251  M_evaluation.setTestCFE ( M_testCFE_unadapted );
252  M_evaluation.setSolutionCFE ( M_solutionCFE_unadapted );
253 }
254 
255 template < typename MeshType,
256  typename TestSpaceType,
257  typename SolutionSpaceType,
258  typename ExpressionType,
259  typename LSFESpaceType,
260  typename LSVectorType >
261 IntegrateMatrixElementLSAdapted<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType, LSFESpaceType, LSVectorType>::
262 IntegrateMatrixElementLSAdapted (const IntegrateMatrixElementLSAdapted < MeshType, TestSpaceType, SolutionSpaceType,
263  ExpressionType, LSFESpaceType, LSVectorType > & integrator)
264  : M_mesh (integrator.M_mesh),
265  M_QRAdapter (integrator.M_QRAdapter),
266  M_testSpace (integrator.M_testSpace),
267  M_solutionSpace (integrator.M_solutionSpace),
268  M_evaluation (integrator.M_evaluation),
269 
270  M_globalCFE_unadapted (new ETCurrentFE<3, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), M_QRAdapter.standardQR() ) ),
271  M_globalCFE_adapted (new ETCurrentFE<3, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), M_QRAdapter.standardQR() ) ),
272 
273  M_testCFE_unadapted (new ETCurrentFE<3, TestSpaceType::S_fieldDim> (M_testSpace->refFE(), M_testSpace->geoMap(), M_QRAdapter.standardQR() ) ),
274  M_testCFE_adapted (new ETCurrentFE<3, TestSpaceType::S_fieldDim> (M_testSpace->refFE(), M_testSpace->geoMap(), M_QRAdapter.standardQR() ) ),
275 
276  M_solutionCFE_unadapted (new ETCurrentFE<3, SolutionSpaceType::S_fieldDim> (M_solutionSpace->refFE(), M_solutionSpace->geoMap(), M_QRAdapter.standardQR() ) ),
277  M_solutionCFE_adapted (new ETCurrentFE<3, SolutionSpaceType::S_fieldDim> (M_solutionSpace->refFE(), M_solutionSpace->geoMap(), M_QRAdapter.standardQR() ) ),
278 
279 
280  M_elementalMatrix (integrator.M_elementalMatrix)
281 {
282  M_evaluation.setQuadrature (M_QRAdapter.standardQR() );
283  M_evaluation.setGlobalCFE ( M_globalCFE_unadapted );
284  M_evaluation.setTestCFE ( M_testCFE_unadapted );
285  M_evaluation.setSolutionCFE ( M_solutionCFE_unadapted );
286 }
287 
288 template < typename MeshType,
289  typename TestSpaceType,
290  typename SolutionSpaceType,
291  typename ExpressionType,
292  typename LSFESpaceType,
293  typename LSVectorType >
294 IntegrateMatrixElementLSAdapted<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType, LSFESpaceType, LSVectorType>::
295 ~IntegrateMatrixElementLSAdapted()
296 {
297  delete M_globalCFE_unadapted;
298  delete M_globalCFE_adapted;
299  delete M_testCFE_unadapted;
300  delete M_testCFE_adapted;
301  delete M_solutionCFE_unadapted;
302  delete M_solutionCFE_adapted;
303 }
304 
305 // ===================================================
306 // Methods
307 // ===================================================
308 
309 template < typename MeshType,
310  typename TestSpaceType,
311  typename SolutionSpaceType,
312  typename ExpressionType,
313  typename LSFESpaceType,
314  typename LSVectorType >
315 void
316 IntegrateMatrixElementLSAdapted<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType, LSFESpaceType, LSVectorType>::
317 check (std::ostream& out)
318 {
319  out << " Checking the integration : " << std::endl;
320  M_evaluation.display (out);
321  out << std::endl;
322  out << " Elemental matrix : " << std::endl;
323  M_elementalMatrix.showMe (out);
324  out << std::endl;
325 }
326 
327 template < typename MeshType,
328  typename TestSpaceType,
329  typename SolutionSpaceType,
330  typename ExpressionType,
331  typename LSFESpaceType,
332  typename LSVectorType >
333 template <typename MatrixType>
334 void
335 IntegrateMatrixElementLSAdapted<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType, LSFESpaceType, LSVectorType>::
336 addTo (MatrixType& mat)
337 {
338  const UInt nbElements (M_mesh->numElements() );
339  const UInt nbQuadPt_unadapted (M_QRAdapter.standardQR().nbQuadPt() );
340  const UInt nbTestDof (M_testSpace->refFE().nbDof() );
341  const UInt nbSolutionDof (M_solutionSpace->refFE().nbDof() );
342 
343  bool isPreviousAdapted (true);
344 
345  for (UInt iElement (0); iElement < nbElements; ++iElement)
346  {
347  // Zeros out the matrix
348  M_elementalMatrix.zero();
349 
350  // Update the adapter
351  M_QRAdapter.update (iElement);
352 
353 
354  if ( M_QRAdapter.isAdaptedElement() )
355  {
356  // Reset the QR in the CurrentFEs
357  M_globalCFE_adapted->setQuadratureRule ( M_QRAdapter.adaptedQR() );
358  M_testCFE_adapted->setQuadratureRule ( M_QRAdapter.adaptedQR() );
359  M_solutionCFE_adapted->setQuadratureRule ( M_QRAdapter.adaptedQR() );
360 
361  // Reset the Evaluation
362  M_evaluation.setGlobalCFE ( M_globalCFE_adapted );
363  M_evaluation.setTestCFE ( M_testCFE_adapted );
364  M_evaluation.setSolutionCFE ( M_solutionCFE_adapted );
365 
366  M_evaluation.setQuadrature ( M_QRAdapter.adaptedQR() );
367 
368  // Update the currentFEs
369  M_globalCFE_adapted->update (M_mesh->element (iElement), evaluation_Type::S_globalUpdateFlag | ET_UPDATE_WDET);
370  M_testCFE_adapted->update (M_mesh->element (iElement), evaluation_Type::S_testUpdateFlag);
371  M_solutionCFE_adapted->update (M_mesh->element (iElement), evaluation_Type::S_solutionUpdateFlag);
372 
373  // Update the evaluation
374  M_evaluation.update (iElement);
375 
376  // Loop on the blocks
377 
378  for (UInt iblock (0); iblock < TestSpaceType::S_fieldDim; ++iblock)
379  {
380  for (UInt jblock (0); jblock < SolutionSpaceType::S_fieldDim; ++jblock)
381  {
382 
383  // Set the row global indices in the local matrix
384  for (UInt i (0); i < nbTestDof; ++i)
385  {
386  M_elementalMatrix.setRowIndex
387  (i + iblock * nbTestDof,
388  M_testSpace->dof().localToGlobalMap (iElement, i) + iblock * M_testSpace->dof().numTotalDof() );
389  }
390 
391  // Set the column global indices in the local matrix
392  for (UInt j (0); j < nbSolutionDof; ++j)
393  {
394  M_elementalMatrix.setColumnIndex
395  (j + jblock * nbSolutionDof,
396  M_solutionSpace->dof().localToGlobalMap (iElement, j) + jblock * M_solutionSpace->dof().numTotalDof() );
397  }
398 
399  for (UInt iQuadPt (0); iQuadPt < M_QRAdapter.adaptedQR().nbQuadPt(); ++iQuadPt)
400  {
401  for (UInt i (0); i < nbTestDof; ++i)
402  {
403  for (UInt j (0); j < nbSolutionDof; ++j)
404  {
405  M_elementalMatrix.element (i + iblock * nbTestDof, j + jblock * nbSolutionDof) +=
406  M_evaluation.value_qij (iQuadPt, i + iblock * nbTestDof, j + jblock * nbSolutionDof)
407  * M_globalCFE_adapted->wDet (iQuadPt);
408 
409  }
410  }
411  }
412  }
413  }
414 
415  isPreviousAdapted = true;
416  }
417  else
418  {
419  if (isPreviousAdapted)
420  {
421  // Reset the Evaluation
422  M_evaluation.setGlobalCFE ( M_globalCFE_unadapted );
423  M_evaluation.setTestCFE ( M_testCFE_unadapted );
424  M_evaluation.setSolutionCFE ( M_solutionCFE_unadapted );
425 
426  M_evaluation.setQuadrature ( M_QRAdapter.standardQR() );
427 
428  }
429 
430  // Update the currentFEs
431  M_globalCFE_unadapted->update (M_mesh->element (iElement), evaluation_Type::S_globalUpdateFlag | ET_UPDATE_WDET);
432  M_testCFE_unadapted->update (M_mesh->element (iElement), evaluation_Type::S_testUpdateFlag);
433  M_solutionCFE_unadapted->update (M_mesh->element (iElement), evaluation_Type::S_solutionUpdateFlag);
434 
435  // Update the evaluation
436  M_evaluation.update (iElement);
437 
438  // Loop on the blocks
439 
440  for (UInt iblock (0); iblock < TestSpaceType::S_fieldDim; ++iblock)
441  {
442  for (UInt jblock (0); jblock < SolutionSpaceType::S_fieldDim; ++jblock)
443  {
444 
445  // Set the row global indices in the local matrix
446  for (UInt i (0); i < nbTestDof; ++i)
447  {
448  M_elementalMatrix.setRowIndex
449  (i + iblock * nbTestDof,
450  M_testSpace->dof().localToGlobalMap (iElement, i) + iblock * M_testSpace->dof().numTotalDof() );
451  }
452 
453  // Set the column global indices in the local matrix
454  for (UInt j (0); j < nbSolutionDof; ++j)
455  {
456  M_elementalMatrix.setColumnIndex
457  (j + jblock * nbSolutionDof,
458  M_solutionSpace->dof().localToGlobalMap (iElement, j) + jblock * M_solutionSpace->dof().numTotalDof() );
459  }
460 
461  for (UInt iQuadPt (0); iQuadPt < nbQuadPt_unadapted; ++iQuadPt)
462  {
463  for (UInt i (0); i < nbTestDof; ++i)
464  {
465  for (UInt j (0); j < nbSolutionDof; ++j)
466  {
467  M_elementalMatrix.element (i + iblock * nbTestDof, j + jblock * nbSolutionDof) +=
468  M_evaluation.value_qij (iQuadPt, i + iblock * nbTestDof, j + jblock * nbSolutionDof)
469  * M_globalCFE_unadapted->wDet (iQuadPt);
470 
471  }
472  }
473  }
474  }
475  }
476 
477 
478  isPreviousAdapted = false;
479  }
480 
481 
482  M_elementalMatrix.pushToGlobal (mat);
483  }
484 }
485 
486 
487 } // Namespace ExpressionAssembly
488 
489 } // Namespace LifeV
490 
491 #endif // INTEGRATE_MATRIX_ELEMENT_LS_ADAPTED_HPP
ASSERT
#define ASSERT(X, A)
Definition: LifeAssert.hpp:90