LifeV
IntegrateVectorElementLSAdapted.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 IntegrateVectorElementLSAdapted class.
31 
32  @date 06/2011
33  @author Samuel Quinodoz <samuel.quinodoz@epfl.ch>
34  */
35 
36 #ifndef INTEGRATE_VECTOR_ELEMENT_LS_ADAPTED_HPP
37 #define INTEGRATE_VECTOR_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/ETVectorElemental.hpp>
48 
49 #include <boost/shared_ptr.hpp>
50 
51 
52 
53 namespace LifeV
54 {
55 
56 namespace ExpressionAssembly
57 {
58 
59 //! The class to actually perform the loop over the elements to assemble a vector using special QR.
60 /*!
61  @author Samuel Quinodoz <samuel.quinodoz@epfl.ch>
62 
63  This class is used to store the data required for the assembly of a vector and
64  perform that assembly with a loop over the elements, and then, for each elements,
65  using the Evaluation corresponding to the Expression (this convertion is done
66  within a typedef).
67 
68  The quadrature used is given by the structure adapting the quadrature to the
69  situation.
70  */
71 template < typename MeshType,
72  typename TestSpaceType,
73  typename ExpressionType,
74  typename LSFESpaceType,
75  typename LSVectorType >
76 class IntegrateVectorElementLSAdapted
77 {
78 public:
79 
80  //! @name Public Types
81  //@{
82 
83  //! Type of the Evaluation
84  typedef typename ExpressionToEvaluation < ExpressionType,
85  TestSpaceType::S_fieldDim,
86  0,
87  3 >::evaluation_Type evaluation_Type;
88 
89  typedef LevelSetQRAdapter<LSFESpaceType, LSVectorType> QRAdapter_Type;
90 
91  //@}
92 
93 
94  //! @name Constructors, destructor
95  //@{
96 
97  //! Full data constructor
98  IntegrateVectorElementLSAdapted (const std::shared_ptr<MeshType>& mesh,
99  const QRAdapter_Type& QRAdapter,
100  const std::shared_ptr<TestSpaceType>& testSpace,
101  const ExpressionType& expression);
102 
103  //! Copy constructor
104  IntegrateVectorElementLSAdapted ( const IntegrateVectorElementLSAdapted
105  < MeshType, TestSpaceType, ExpressionType, LSFESpaceType, LSVectorType>& integrator);
106 
107  //! Destructor
108  ~IntegrateVectorElementLSAdapted();
109 
110  //@}
111 
112 
113  //! @name Operator
114  //@{
115 
116  //! Operator wrapping the addTo method
117  template <typename VectorType>
118  inline void operator>> (VectorType& vector)
119  {
120  addTo (vector);
121  }
122 
123  template <typename VectorType>
124  inline void operator>> (std::shared_ptr<VectorType> vector)
125  {
126  addTo (*vector);
127  }
128 
129  //@}
130 
131  //! @name Methods
132  //@{
133 
134  //! Ouput method
135  void check (std::ostream& out = std::cout);
136 
137  //! Method that performs the assembly
138  /*!
139  The loop over the elements is located right
140  in this method. Everything for the assembly is then
141  performed: update the values, update the local vector,
142  sum over the quadrature nodes, assemble in the global
143  vector.
144  */
145  template <typename VectorType>
146  void addTo (VectorType& vec);
147 
148  template <typename VectorType>
149  void addTo (std::shared_ptr<VectorType> vec)
150  {
151  addTo (*vec);
152  }
153 
154  //@}
155 
156 private:
157 
158  //! @name Private Methods
159  //@{
160 
161  // No default constructor
162  IntegrateVectorElementLSAdapted();
163 
164  //@}
165 
166  // Pointer on the mesh
167  std::shared_ptr<MeshType> M_mesh;
168 
169  // Quadrature adapter
170  QRAdapter_Type M_QRAdapter;
171 
172  // Shared pointer on the Space
173  std::shared_ptr<TestSpaceType> M_testSpace;
174 
175  // Tree to compute the values for the assembly
176  evaluation_Type M_evaluation;
177 
178  // CurrentFEs are duplicated to consider the
179  // adapted and unadapted quadrature rules.
180  ETCurrentFE<3, 1>* M_globalCFE_unadapted;
181  ETCurrentFE<3, 1>* M_globalCFE_adapted;
182 
183  ETCurrentFE<3, TestSpaceType::S_fieldDim>* M_testCFE_unadapted;
184  ETCurrentFE<3, TestSpaceType::S_fieldDim>* M_testCFE_adapted;
185 
186  //ETVectorElemental<1> M_elementalVector;
187  ETVectorElemental M_elementalVector;
188 };
189 
190 
191 // ===================================================
192 // IMPLEMENTATION
193 // ===================================================
194 
195 // ===================================================
196 // Constructors & Destructor
197 // ===================================================
198 
199 template < typename MeshType,
200  typename TestSpaceType,
201  typename ExpressionType,
202  typename LSFESpaceType,
203  typename LSVectorType >
204 IntegrateVectorElementLSAdapted < MeshType, TestSpaceType, ExpressionType, LSFESpaceType, LSVectorType>::
205 IntegrateVectorElementLSAdapted (const std::shared_ptr<MeshType>& mesh,
206  const QRAdapter_Type& QRAdapter,
207  const std::shared_ptr<TestSpaceType>& testSpace,
208  const ExpressionType& expression)
209  : M_mesh (mesh),
210  M_QRAdapter (QRAdapter),
211  M_testSpace (testSpace),
212  M_evaluation (expression),
213 
214  M_globalCFE_unadapted (new ETCurrentFE<3, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), M_QRAdapter.standardQR() ) ),
215  M_globalCFE_adapted (new ETCurrentFE<3, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), M_QRAdapter.standardQR() ) ),
216 
217  M_testCFE_unadapted (new ETCurrentFE<3, TestSpaceType::S_fieldDim> (testSpace->refFE(), testSpace->geoMap(), M_QRAdapter.standardQR() ) ),
218  M_testCFE_adapted (new ETCurrentFE<3, TestSpaceType::S_fieldDim> (testSpace->refFE(), testSpace->geoMap(), M_QRAdapter.standardQR() ) ),
219 
220  M_elementalVector (TestSpaceType::S_fieldDim * testSpace->refFE().nbDof() )
221 {
222  M_evaluation.setQuadrature (M_QRAdapter.standardQR() );
223  M_evaluation.setGlobalCFE (M_globalCFE_unadapted);
224  M_evaluation.setTestCFE (M_testCFE_unadapted);
225 }
226 
227 
228 template < typename MeshType,
229  typename TestSpaceType,
230  typename ExpressionType,
231  typename LSFESpaceType,
232  typename LSVectorType >
233 IntegrateVectorElementLSAdapted < MeshType, TestSpaceType, ExpressionType, LSFESpaceType, LSVectorType>::
234 IntegrateVectorElementLSAdapted ( const IntegrateVectorElementLSAdapted
235  < MeshType, TestSpaceType, ExpressionType, LSFESpaceType, LSVectorType>& integrator)
236  : M_mesh (integrator.M_mesh),
237  M_QRAdapter (integrator.M_QRAdapter),
238  M_testSpace (integrator.M_testSpace),
239  M_evaluation (integrator.M_evaluation),
240 
241  M_globalCFE_unadapted (new ETCurrentFE<3, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), M_QRAdapter.standardQR() ) ),
242  M_globalCFE_adapted (new ETCurrentFE<3, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), M_QRAdapter.standardQR() ) ),
243 
244  M_testCFE_unadapted (new ETCurrentFE<3, TestSpaceType::S_fieldDim> (M_testSpace->refFE(), M_testSpace->geoMap(), M_QRAdapter.standardQR() ) ),
245  M_testCFE_adapted (new ETCurrentFE<3, TestSpaceType::S_fieldDim> (M_testSpace->refFE(), M_testSpace->geoMap(), M_QRAdapter.standardQR() ) ),
246 
247  M_elementalVector (integrator.M_elementalVector)
248 {
249  M_evaluation.setQuadrature (M_QRAdapter.standardQR() );
250  M_evaluation.setGlobalCFE (M_globalCFE_unadapted);
251  M_evaluation.setTestCFE (M_testCFE_unadapted);
252 }
253 
254 
255 template < typename MeshType,
256  typename TestSpaceType,
257  typename ExpressionType,
258  typename LSFESpaceType,
259  typename LSVectorType >
260 IntegrateVectorElementLSAdapted < MeshType, TestSpaceType, ExpressionType, LSFESpaceType, LSVectorType>::
261 ~IntegrateVectorElementLSAdapted()
262 {
263  delete M_globalCFE_unadapted;
264  delete M_globalCFE_adapted;
265  delete M_testCFE_unadapted;
266  delete M_testCFE_adapted;
267 }
268 
269 // ===================================================
270 // Methods
271 // ===================================================
272 
273 template < typename MeshType,
274  typename TestSpaceType,
275  typename ExpressionType,
276  typename LSFESpaceType,
277  typename LSVectorType >
278 void
279 IntegrateVectorElementLSAdapted < MeshType, TestSpaceType, ExpressionType, LSFESpaceType, LSVectorType>::
280 check (std::ostream& out)
281 {
282  out << " Checking the integration : " << std::endl;
283  M_evaluation.display (out);
284  out << std::endl;
285  out << " Elemental vector : " << std::endl;
286  M_elementalVector.showMe (out);
287  out << std::endl;
288 }
289 
290 template < typename MeshType,
291  typename TestSpaceType,
292  typename ExpressionType,
293  typename LSFESpaceType,
294  typename LSVectorType >
295 template <typename VectorType>
296 void
297 IntegrateVectorElementLSAdapted < MeshType, TestSpaceType, ExpressionType, LSFESpaceType, LSVectorType>::
298 addTo (VectorType& vec)
299 {
300  const UInt nbElements (M_mesh->numElements() );
301  const UInt nbQuadPt_unadapted (M_QRAdapter.standardQR().nbQuadPt() );
302  const UInt nbTestDof (M_testSpace->refFE().nbDof() );
303 
304  // Even if this is not the case, this
305  // prevents many potentially buggy behaviours
306  // and does not harm.
307  bool isPreviousAdapted (true);
308 
309  for (UInt iElement (0); iElement < nbElements; ++iElement)
310  {
311  // Zeros out the elemental vector
312  M_elementalVector.zero();
313 
314  // Update the quadrature rule
315  M_QRAdapter.update (iElement);
316 
317 
318  if (M_QRAdapter.isAdaptedElement() )
319  {
320  // Change the QR
321  M_evaluation.setQuadrature ( M_QRAdapter.adaptedQR() );
322  M_globalCFE_adapted->setQuadratureRule ( M_QRAdapter.adaptedQR() );
323  M_testCFE_adapted->setQuadratureRule ( M_QRAdapter.adaptedQR() );
324 
325  // Reset the CFEs in the evaluation
326  M_evaluation.setGlobalCFE ( M_globalCFE_adapted );
327  M_evaluation.setTestCFE ( M_testCFE_adapted );
328 
329  // Update the currentFEs
330  M_globalCFE_adapted->update (M_mesh->element (iElement), evaluation_Type::S_globalUpdateFlag | ET_UPDATE_WDET);
331  M_testCFE_adapted->update (M_mesh->element (iElement), evaluation_Type::S_testUpdateFlag);
332 
333  // Update the evaluation
334  M_evaluation.update (iElement);
335 
336  // Loop on the blocks
337 
338  for (UInt iblock (0); iblock < TestSpaceType::S_fieldDim; ++iblock)
339  {
340  // Set the row global indices in the local vector
341  for (UInt i (0); i < nbTestDof; ++i)
342  {
343  M_elementalVector.setRowIndex
344  (i + iblock * nbTestDof,
345  M_testSpace->dof().localToGlobalMap (iElement, i) + iblock * M_testSpace->dof().numTotalDof() );
346  }
347 
348  // Make the assembly
349  for (UInt iQuadPt (0); iQuadPt < M_QRAdapter.adaptedQR().nbQuadPt(); ++iQuadPt)
350  {
351  for (UInt i (0); i < nbTestDof; ++i)
352  {
353  M_elementalVector.element (i + iblock * nbTestDof) +=
354  M_evaluation.value_qi (iQuadPt, i + iblock * nbTestDof)
355  * M_globalCFE_adapted->wDet (iQuadPt);
356  }
357  }
358  }
359 
360  M_elementalVector.pushToGlobal (vec);
361 
362  isPreviousAdapted = true;
363 
364  }
365  else
366  {
367  if (isPreviousAdapted)
368  {
369  M_evaluation.setQuadrature ( M_QRAdapter.standardQR() );
370  M_evaluation.setGlobalCFE ( M_globalCFE_unadapted );
371  M_evaluation.setTestCFE ( M_testCFE_unadapted );
372  }
373 
374  // Update the currentFEs
375  M_globalCFE_unadapted->update (M_mesh->element (iElement), evaluation_Type::S_globalUpdateFlag | ET_UPDATE_WDET);
376  M_testCFE_unadapted->update (M_mesh->element (iElement), evaluation_Type::S_testUpdateFlag);
377 
378  // Update the evaluation
379  M_evaluation.update (iElement);
380 
381  // Loop on the blocks
382 
383  for (UInt iblock (0); iblock < TestSpaceType::S_fieldDim; ++iblock)
384  {
385  // Set the row global indices in the local vector
386  for (UInt i (0); i < nbTestDof; ++i)
387  {
388  M_elementalVector.setRowIndex
389  (i + iblock * nbTestDof,
390  M_testSpace->dof().localToGlobalMap (iElement, i) + iblock * M_testSpace->dof().numTotalDof() );
391  }
392 
393  // Make the assembly
394  for (UInt iQuadPt (0); iQuadPt < nbQuadPt_unadapted; ++iQuadPt)
395  {
396  for (UInt i (0); i < nbTestDof; ++i)
397  {
398  M_elementalVector.element (i + iblock * nbTestDof) +=
399  M_evaluation.value_qi (iQuadPt, i + iblock * nbTestDof)
400  * M_globalCFE_unadapted->wDet (iQuadPt);
401  }
402  }
403  }
404 
405  M_elementalVector.pushToGlobal (vec);
406 
407  isPreviousAdapted = false;
408 
409  }
410  }
411 }
412 
413 
414 } // Namespace ExpressionAssembly
415 
416 } // Namespace LifeV
417 
418 #endif // INTEGRATE_VECTOR_ELEMENT_LS_ADAPTED_HPP