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