LifeV
ComputeFineScaleVelocity.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 computation of fine scale velocity inside an element.
31 
32  @date 06/2011
33  @author Davide Forti <davide.forti@epfl.ch>
34  */
35 
36 #ifndef COMPUTE_FINESCALEVELOCITY_ELEMENT_HPP
37 #define COMPUTE_FINESCALEVELOCITY_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 #include <lifev/eta/expression/EvaluationPhiI.hpp>
48 
49 #include <lifev/eta/array/ETVectorElemental.hpp>
50 
51 
52 
53 
54 namespace LifeV
55 {
56 
57 namespace ExpressionAssembly
58 {
59 
60 //! The class to actually perform the loop over the elements to compute the fine scale velocity at the center of each element
61 /*!
62  @author Davide Forti <davide.forti@epfl.ch>
63 
64  */
65 template < typename MeshType, typename TestSpaceType, typename ExpressionType, typename QRAdapterType>
66 class ComputeFineScaleVelocity
67 {
68 public:
69 
70  //! @name Public Types
71  //@{
72 
73  //! Type of the Evaluation
74  typedef typename ExpressionToEvaluation < ExpressionType,
75  TestSpaceType::field_dim,
76  0,
77  MeshType::S_geoDimensions >::evaluation_Type evaluation_Type;
78 
79  //@}
80 
81 
82  //! @name Constructors, destructor
83  //@{
84 
85  //! Full data constructor
86  ComputeFineScaleVelocity (const std::shared_ptr<MeshType>& mesh,
87  const QRAdapterType& qrAdapter,
88  const std::shared_ptr<TestSpaceType>& testSpace,
89  const ExpressionType& expression,
90  const UInt offset = 0);
91 
92  //! Copy constructor
93  ComputeFineScaleVelocity ( const IntegrateVectorElement < MeshType, TestSpaceType, ExpressionType, QRAdapterType>& integrator);
94 
95  //! Destructor
96  ~ComputeFineScaleVelocity();
97 
98  //@}
99 
100 
101  //! @name Operator
102  //@{
103 
104  //! Operator wrapping the addTo method
105  template <typename VectorType>
106  inline void operator>> (VectorType& vector)
107  {
108  compute (vector);
109  }
110 
111  //@}
112 
113  //! @name Methods
114  //@{
115 
116  //! Ouput method
117  void check (std::ostream& out = std::cout);
118 
119  //! Method that performs the assembly
120  /*!
121  The loop over the elements is located right
122  in this method. Everything for the assembly is then
123  performed: update the values, update the local vector,
124  sum over the quadrature nodes, assemble in the global
125  vector.
126  */
127  template <typename VectorType>
128  void compute (VectorType& vec);
129 
130  //@}
131 
132 private:
133 
134  //! @name Private Methods
135  //@{
136 
137  // No default constructor
138  ComputeFineScaleVelocity();
139 
140  //@}
141 
142  // Pointer on the mesh
143  std::shared_ptr<MeshType> M_mesh;
144 
145  // Quadrature to be used
146  QRAdapterType M_qrAdapter;
147 
148  // Shared pointer on the Space
149  std::shared_ptr<TestSpaceType> M_testSpace;
150 
151  // Tree to compute the values for the assembly
152  evaluation_Type M_evaluation;
153 
154  ETCurrentFE<MeshType::S_geoDimensions, 1>* M_globalCFE_std;
155  ETCurrentFE<MeshType::S_geoDimensions, 1>* M_globalCFE_adapted;
156 
157  ETCurrentFE<TestSpaceType::space_dim, TestSpaceType::field_dim>* M_testCFE_std;
158  ETCurrentFE<TestSpaceType::space_dim, TestSpaceType::field_dim>* M_testCFE_adapted;
159 
160  ETVectorElemental M_elementalVector;
161 
162  // Offset
163  UInt M_offset;
164 };
165 
166 
167 // ===================================================
168 // IMPLEMENTATION
169 // ===================================================
170 
171 // ===================================================
172 // Constructors & Destructor
173 // ===================================================
174 
175 template < typename MeshType, typename TestSpaceType, typename ExpressionType, typename QRAdapterType>
176 ComputeFineScaleVelocity < MeshType, TestSpaceType, ExpressionType, QRAdapterType>::
177 ComputeFineScaleVelocity (const std::shared_ptr<MeshType>& mesh,
178  const QRAdapterType& qrAdapter,
179  const std::shared_ptr<TestSpaceType>& testSpace,
180  const ExpressionType& expression,
181  const UInt offset)
182  : M_mesh (mesh),
183  M_qrAdapter (qrAdapter),
184  M_testSpace (testSpace),
185  M_evaluation (expression),
186 
187  M_testCFE_std (new ETCurrentFE<TestSpaceType::space_dim, TestSpaceType::field_dim> (testSpace->refFE(), testSpace->geoMap(), qrAdapter.standardQR() ) ),
188  M_testCFE_adapted (new ETCurrentFE<TestSpaceType::space_dim, TestSpaceType::field_dim> (testSpace->refFE(), testSpace->geoMap(), qrAdapter.standardQR() ) ),
189 
190  M_elementalVector (TestSpaceType::field_dim * testSpace->refFE().nbDof() ),
191 
192  M_offset (offset)
193 {
194  switch (MeshType::geoShape_Type::BasRefSha::S_shape)
195  {
196  case LINE:
197  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feSegP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
198  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feSegP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
199  break;
200  case TRIANGLE:
201  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTriaP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
202  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTriaP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
203  break;
204  case QUAD:
205  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feQuadQ0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
206  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feQuadQ0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
207  break;
208  case TETRA:
209  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
210  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
211  break;
212  case HEXA:
213  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feHexaQ0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
214  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feHexaQ0, geometricMapFromMesh<MeshType>(), qrAdapter.standardQR() );
215  break;
216  default:
217  ERROR_MSG ("Unrecognized element shape");
218  }
219  M_evaluation.setQuadrature ( qrAdapter.standardQR() );
220 
221  M_evaluation.setGlobalCFE (M_globalCFE_std);
222  M_evaluation.setTestCFE (M_testCFE_std);
223 }
224 
225 
226 template < typename MeshType, typename TestSpaceType, typename ExpressionType, typename QRAdapterType>
227 ComputeFineScaleVelocity < MeshType, TestSpaceType, ExpressionType, QRAdapterType>::
228 ComputeFineScaleVelocity ( const IntegrateVectorElement < MeshType, TestSpaceType, ExpressionType, QRAdapterType>& integrator)
229  : M_mesh (integrator.M_mesh),
230  M_qrAdapter (integrator.M_qrAdapter),
231  M_testSpace (integrator.M_testSpace),
232  M_evaluation (integrator.M_evaluation),
233 
234  M_testCFE_std (new ETCurrentFE<TestSpaceType::space_dim, TestSpaceType::field_dim> (M_testSpace->refFE(), M_testSpace->geoMap(), integrator.M_qrAdapter.standardQR() ) ),
235  M_testCFE_adapted (new ETCurrentFE<TestSpaceType::space_dim, TestSpaceType::field_dim> (M_testSpace->refFE(), M_testSpace->geoMap(), integrator.M_qrAdapter.standardQR() ) ),
236 
237  M_elementalVector (integrator.M_elementalVector),
238  M_offset (integrator.M_offset)
239 {
240  switch (MeshType::geoShape_Type::BasRefSha::S_shape)
241  {
242  case LINE:
243  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feSegP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
244  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feSegP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
245  break;
246  case TRIANGLE:
247  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTriaP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
248  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTriaP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
249  break;
250  case QUAD:
251  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feQuadQ0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
252  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feQuadQ0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
253  break;
254  case TETRA:
255  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
256  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feTetraP0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
257  break;
258  case HEXA:
259  M_globalCFE_std = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feHexaQ0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
260  M_globalCFE_adapted = new ETCurrentFE<MeshType::S_geoDimensions, 1> (feHexaQ0, geometricMapFromMesh<MeshType>(), integrator.M_qrAdapter.standardQR() );
261  break;
262  default:
263  ERROR_MSG ("Unrecognized element shape");
264  }
265  M_evaluation.setQuadrature (integrator.M_qrAdapter.standardQR() );
266  M_evaluation.setGlobalCFE (M_globalCFE_std);
267  M_evaluation.setTestCFE (M_testCFE_std);
268 }
269 
270 
271 template < typename MeshType, typename TestSpaceType, typename ExpressionType, typename QRAdapterType>
272 ComputeFineScaleVelocity < MeshType, TestSpaceType, ExpressionType, QRAdapterType>::
273 ~ComputeFineScaleVelocity()
274 {
275  delete M_globalCFE_std;
276  delete M_globalCFE_adapted;
277  delete M_testCFE_std;
278  delete M_testCFE_adapted;
279 }
280 
281 // ===================================================
282 // Methods
283 // ===================================================
284 
285 template < typename MeshType, typename TestSpaceType, typename ExpressionType, typename QRAdapterType>
286 void
287 ComputeFineScaleVelocity < MeshType, TestSpaceType, ExpressionType, QRAdapterType>::
288 check (std::ostream& out)
289 {
290  out << " Checking the integration : " << std::endl;
291  M_evaluation.display (out);
292  out << std::endl;
293  out << " Elemental vector : " << std::endl;
294  M_elementalVector.showMe (out);
295  out << std::endl;
296 }
297 
298 template < typename MeshType, typename TestSpaceType, typename ExpressionType, typename QRAdapterType>
299 template <typename VectorType>
300 void
301 ComputeFineScaleVelocity < MeshType, TestSpaceType, ExpressionType, QRAdapterType>::
302 compute (VectorType& vec)
303 {
304  UInt nbElements (M_mesh->numElements() );
305  UInt nbQuadPt_std (M_qrAdapter.standardQR().nbQuadPt() );
306  UInt nbTestDof (M_testSpace->refFE().nbDof() );
307 
308  // Defaulted to true for security
309  bool isPreviousAdapted (true);
310 
311  VectorSmall<3> elementalVelocities;
312 
313  for (UInt iElement (0); iElement < nbElements; ++iElement)
314  {
315 
316  // Update the quadrature rule adapter
317  M_qrAdapter.update (iElement);
318 
319  // Check if the last one was adapted
320  if (isPreviousAdapted)
321  {
322  M_evaluation.setQuadrature ( M_qrAdapter.standardQR() );
323  M_evaluation.setGlobalCFE ( M_globalCFE_std );
324  M_evaluation.setTestCFE ( M_testCFE_std );
325 
326  isPreviousAdapted = false;
327  }
328 
329  // Update the currentFEs
330  M_globalCFE_std->update (M_mesh->element (iElement), evaluation_Type::S_globalUpdateFlag | ET_UPDATE_WDET);
331  M_testCFE_std->update (M_mesh->element (iElement), evaluation_Type::S_testUpdateFlag);
332 
333  // Update the evaluation
334  M_evaluation.update (iElement);
335 
336  elementalVelocities = M_evaluation.value_qi (0, 0);
337 
338  for (UInt i (0); i < 3; ++i)
339  {
340  vec[M_mesh->element(iElement).id() + i * M_mesh->numGlobalElements() ] = elementalVelocities(i);
341  elementalVelocities(i) = 0.0;
342  }
343  }
344 
345 }
346 
347 } // Namespace ExpressionAssembly
348 
349 } // Namespace LifeV
350 
351 #endif
ERROR_MSG
#define ERROR_MSG(A)
Definition: LifeAssert.hpp:69