LifeV
IntegrateVectorFaceID.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 INTEGRATE_VECTOR_FACE_ID_HPP
37 #define INTEGRATE_VECTOR_FACE_ID_HPP
38 
39 #include <lifev/core/LifeV.hpp>
40 
41 #include <lifev/eta/fem/QuadratureBoundary.hpp>
42 #include <lifev/eta/fem/ETCurrentFE.hpp>
43 #include <lifev/eta/fem/ETCurrentBDFE.hpp>
44 #include <lifev/eta/fem/MeshGeometricMap.hpp>
45 
46 #include <lifev/eta/expression/ExpressionToEvaluation.hpp>
47 
48 #include <lifev/eta/array/ETVectorElemental.hpp>
49 
50 #include <boost/shared_ptr.hpp>
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 assemble a vector
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 vector 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 template < typename MeshType, typename TestSpaceType, typename ExpressionType>
70 class IntegrateVectorFaceID
71 {
72 public:
73 
74  //! @name Public Types
75  //@{
76 
77  //! Type of the Evaluation
78  typedef typename ExpressionToEvaluation < ExpressionType,
79  TestSpaceType::field_dim,
80  0,
81  3 >::evaluation_Type evaluation_Type;
82 
83  //@}
84 
85 
86  //! @name Constructors, destructor
87  //@{
88 
89  //! Full data constructor
90  IntegrateVectorFaceID (const std::shared_ptr<MeshType>& mesh,
91  const UInt boundaryID,
92  const QuadratureBoundary& quadratureBD,
93  const std::shared_ptr<TestSpaceType>& testSpace,
94  const ExpressionType& expression);
95 
96  //! Copy constructor
97  IntegrateVectorFaceID ( const IntegrateVectorFaceID < MeshType, TestSpaceType, ExpressionType>& integrator);
98 
99  //! Destructor
100  ~IntegrateVectorFaceID();
101 
102  //@}
103 
104 
105  //! @name Operator
106  //@{
107 
108  //! Operator wrapping the addTo method
109  template <typename VectorType>
110  inline void operator>> (VectorType& vector)
111  {
112  addTo (vector);
113  }
114 
115  //! Operator wrapping the addTo method (for shared_ptr)
116  template <typename VectorType>
117  inline void operator>> (std::shared_ptr<VectorType> vector)
118  {
119  addTo (*vector);
120  }
121 
122 
123  //@}
124 
125  //! @name Methods
126  //@{
127 
128  //! Ouput method
129  void check (std::ostream& out = std::cout);
130 
131  //! Method that performs the assembly
132  /*!
133  The loop over the elements is located right
134  in this method. Everything for the assembly is then
135  performed: update the values, update the local vector,
136  sum over the quadrature nodes, assemble in the global
137  vector.
138  */
139  template <typename VectorType>
140  void addTo (VectorType& vec);
141 
142  //@}
143 
144 private:
145 
146  //! @name Private Methods
147  //@{
148 
149  // No default constructor
150  IntegrateVectorFaceID();
151 
152  //@}
153 
154  // Pointer on the mesh
155  std::shared_ptr<MeshType> M_mesh;
156 
157  // Identifier for the boundary
158  UInt M_boundaryId;
159 
160  // Quadrature to be used
161  QuadratureBoundary M_quadratureBoundary;
162 
163  // Shared pointer on the Space
164  std::shared_ptr<TestSpaceType> M_testSpace;
165 
166  // Tree to compute the values for the assembly
167  evaluation_Type M_evaluation;
168 
169  std::vector<ETCurrentBDFE<3>*> M_globalCFE;
170  std::vector<ETCurrentFE<3, TestSpaceType::field_dim>*> M_testCFE;
171 
172  ETVectorElemental M_elementalVector;
173 };
174 
175 
176 // ===================================================
177 // IMPLEMENTATION
178 // ===================================================
179 
180 // ===================================================
181 // Constructors & Destructor
182 // ===================================================
183 
184 template < typename MeshType, typename TestSpaceType, typename ExpressionType>
185 IntegrateVectorFaceID < MeshType, TestSpaceType, ExpressionType>::
186 IntegrateVectorFaceID (const std::shared_ptr<MeshType>& mesh,
187  const UInt boundaryID,
188  const QuadratureBoundary& quadratureBD,
189  const std::shared_ptr<TestSpaceType>& testSpace,
190  const ExpressionType& expression)
191  : M_mesh (mesh),
192  M_boundaryId (boundaryID),
193  M_quadratureBoundary (quadratureBD),
194  M_testSpace (testSpace),
195  M_evaluation (expression),
196 
197  M_globalCFE (4),
198  M_testCFE (4),
199 
200  M_elementalVector (TestSpaceType::field_dim * testSpace->refFE().nbDof() )
201 {
202  for (UInt i (0); i < 4; ++i)
203  {
204  M_globalCFE[i] = new ETCurrentBDFE<3> (geometricMapFromMesh<MeshType>()
205  , M_quadratureBoundary.qr (i) );
206  M_testCFE[i] = new ETCurrentFE<3, TestSpaceType::field_dim> (testSpace->refFE()
207  , testSpace->geoMap()
208  , M_quadratureBoundary.qr (i) );
209  }
210 
211  // Set the tangent on the different faces
212  std::vector< VectorSmall<3> > t0 (2, VectorSmall<3> (0.0, 0.0, 0.0) );
213  t0[0][0] = 1;
214  t0[0][1] = 0;
215  t0[0][2] = 0;
216  t0[1][0] = 0;
217  t0[1][1] = 1;
218  t0[1][2] = 0;
219  std::vector< VectorSmall<3> > t1 (2, VectorSmall<3> (0.0, 0.0, 0.0) );
220  t1[0][0] = 0;
221  t1[0][1] = 0;
222  t1[0][2] = 1;
223  t1[1][0] = 1;
224  t1[1][1] = 0;
225  t1[1][2] = 0;
226  std::vector< VectorSmall<3> > t2 (2, VectorSmall<3> (0.0, 0.0, 0.0) );
227  //t2[0][0]=-1/std::sqrt(6); t2[0][1]=-1/std::sqrt(6); t2[0][2]=2/std::sqrt(6);
228  //t2[1][0]=-1/std::sqrt(2); t2[1][1]=1/std::sqrt(2); t2[1][2]=0;
229  t2[0][0] = -1;
230  t2[0][1] = 0;
231  t2[0][2] = 1;
232  t2[1][0] = -1;
233  t2[1][1] = 1;
234  t2[1][2] = 0;
235 
236  std::vector< VectorSmall<3> > t3 (2, VectorSmall<3> (0.0, 0.0, 0.0) );
237  t3[0][0] = 0;
238  t3[0][1] = 1;
239  t3[0][2] = 0;
240  t3[1][0] = 0;
241  t3[1][1] = 0;
242  t3[1][2] = 1;
243 
244  M_globalCFE[0]->setRefTangents (t0);
245  M_globalCFE[1]->setRefTangents (t1);
246  M_globalCFE[2]->setRefTangents (t2);
247  M_globalCFE[3]->setRefTangents (t3);
248 
249 
250  M_evaluation.setQuadrature (M_quadratureBoundary.qr (0) );
251  M_evaluation.setGlobalCFE (M_globalCFE[0]);
252  M_evaluation.setTestCFE (M_testCFE[0]);
253 }
254 
255 
256 template < typename MeshType, typename TestSpaceType, typename ExpressionType>
257 IntegrateVectorFaceID < MeshType, TestSpaceType, ExpressionType>::
258 IntegrateVectorFaceID ( const IntegrateVectorFaceID < MeshType, TestSpaceType, ExpressionType>& integrator)
259  : M_mesh (integrator.M_mesh),
260  M_boundaryId (integrator.M_boundaryId),
261  M_quadratureBoundary (integrator.M_quadratureBoundary),
262  M_testSpace (integrator.M_testSpace),
263  M_evaluation (integrator.M_evaluation),
264 
265  M_globalCFE (4),
266  M_testCFE (4),
267 
268  M_elementalVector (integrator.M_elementalVector)
269 {
270  for (UInt i (0); i < 4; ++i)
271  {
272  M_globalCFE[i] = new ETCurrentBDFE<3> (geometricMapFromMesh<MeshType>()
273  , M_quadratureBoundary.qr (i) );
274  M_testCFE[i] = new ETCurrentFE<3, TestSpaceType::field_dim> (M_testSpace->refFE()
275  , M_testSpace->geoMap()
276  , M_quadratureBoundary.qr (i) );
277  }
278 
279  // Set the tangent on the different faces
280  std::vector< VectorSmall<3> > t0 (2, VectorSmall<3> (0.0, 0.0, 0.0) );
281  t0[0][0] = 1;
282  t0[0][1] = 0;
283  t0[0][2] = 0;
284  t0[1][0] = 0;
285  t0[1][1] = 1;
286  t0[1][2] = 0;
287  std::vector< VectorSmall<3> > t1 (2, VectorSmall<3> (0.0, 0.0, 0.0) );
288  t1[0][0] = 0;
289  t1[0][1] = 0;
290  t1[0][2] = 1;
291  t1[1][0] = 1;
292  t1[1][1] = 0;
293  t1[1][2] = 0;
294  std::vector< VectorSmall<3> > t2 (2, VectorSmall<3> (0.0, 0.0, 0.0) );
295  //t2[0][0]=-1/std::sqrt(6); t2[0][1]=-1/std::sqrt(6); t2[0][2]=2/std::sqrt(6);
296  //t2[1][0]=-1/std::sqrt(2); t2[1][1]=1/std::sqrt(2); t2[1][2]=0;
297  t2[0][0] = -1;
298  t2[0][1] = 0;
299  t2[0][2] = 1;
300  t2[1][0] = -1;
301  t2[1][1] = 1;
302  t2[1][2] = 0;
303 
304  std::vector< VectorSmall<3> > t3 (2, VectorSmall<3> (0.0, 0.0, 0.0) );
305  t3[0][0] = 0;
306  t3[0][1] = 1;
307  t3[0][2] = 0;
308  t3[1][0] = 0;
309  t3[1][1] = 0;
310  t3[1][2] = 1;
311 
312  M_globalCFE[0]->setRefTangents (t0);
313  M_globalCFE[1]->setRefTangents (t1);
314  M_globalCFE[2]->setRefTangents (t2);
315  M_globalCFE[3]->setRefTangents (t3);
316 
317 
318  M_evaluation.setQuadrature (M_quadratureBoundary.qr (0) );
319  M_evaluation.setGlobalCFE (M_globalCFE[0]);
320  M_evaluation.setTestCFE (M_testCFE[0]);
321 }
322 
323 
324 template < typename MeshType, typename TestSpaceType, typename ExpressionType>
325 IntegrateVectorFaceID < MeshType, TestSpaceType, ExpressionType>::
326 ~IntegrateVectorFaceID()
327 {
328  for (UInt i (0); i < 4; ++i)
329  {
330  delete M_globalCFE[i];
331  delete M_testCFE[i];
332  }
333 }
334 
335 // ===================================================
336 // Methods
337 // ===================================================
338 
339 template < typename MeshType, typename TestSpaceType, typename ExpressionType>
340 void
341 IntegrateVectorFaceID < MeshType, TestSpaceType, ExpressionType>::
342 check (std::ostream& out)
343 {
344  out << " Checking the integration : " << std::endl;
345  M_evaluation.display (out);
346  out << std::endl;
347  out << " Elemental vector : " << std::endl;
348  M_elementalVector.showMe (out);
349  out << std::endl;
350 }
351 
352 template < typename MeshType, typename TestSpaceType, typename ExpressionType>
353 template <typename VectorType>
354 void
355 IntegrateVectorFaceID < MeshType, TestSpaceType, ExpressionType>::
356 addTo (VectorType& vec)
357 {
358  UInt nbBoundaryFaces (M_mesh->numBFaces() );
359  UInt nbTestDof (M_testSpace->refFE().nbDof() );
360 
361  /*CurrentBoundaryFE origFE(M_testSpace->refFE().boundaryFE(),
362  M_testSpace->geoMap().boundaryMap(),
363  quadRuleTria3pt);*/
364 
365  for (UInt iFace (0); iFace < nbBoundaryFaces; ++iFace)
366  {
367  // Check the identifier
368  if ( M_mesh->face (iFace).markerID() != M_boundaryId )
369  {
370  continue;
371  }
372 
373  //origFE.updateMeasNormalQuadPt(M_mesh->face(iFace));
374 
375  // Zeros out the elemental vector
376  M_elementalVector.zero();
377 
378  // Get the number of the face in the adjacent element
379  UInt faceIDinAdjacentElement (M_mesh->face (iFace).firstAdjacentElementPosition() );
380 
381  // Get the ID of the adjacent element
382  UInt adjacentElementID (M_mesh->face (iFace).firstAdjacentElementIdentity() );
383 
384  // Update the currentFEs
385  M_globalCFE[faceIDinAdjacentElement]
386  ->update (M_mesh->element (adjacentElementID) );
387  M_testCFE[faceIDinAdjacentElement]
388  ->update (M_mesh->element (adjacentElementID), evaluation_Type::S_testUpdateFlag);
389 
390  // Update the evaluation
391  M_evaluation.setQuadrature (M_quadratureBoundary.qr (faceIDinAdjacentElement) );
392  M_evaluation.setGlobalCFE (M_globalCFE[faceIDinAdjacentElement]);
393  M_evaluation.setTestCFE (M_testCFE[faceIDinAdjacentElement]);
394 
395  M_evaluation.update (adjacentElementID);
396 
397  // Loop on the blocks
398  for (UInt iblock (0); iblock < TestSpaceType::field_dim; ++iblock)
399  {
400  // Set the row global indices in the local vector
401  for (UInt i (0); i < nbTestDof; ++i)
402  {
403  M_elementalVector.setRowIndex
404  (i + iblock * nbTestDof,
405  M_testSpace->dof().localToGlobalMap (adjacentElementID, i) + iblock * M_testSpace->dof().numTotalDof() );
406  }
407 
408  // Make the assembly
409  for (UInt iQuadPt (0); iQuadPt < M_quadratureBoundary.qr (faceIDinAdjacentElement).nbQuadPt(); ++iQuadPt)
410  {
411  for (UInt i (0); i < nbTestDof; ++i)
412  {
413  M_elementalVector.element (i + iblock * nbTestDof) +=
414  M_evaluation.value_qi (iQuadPt, i + iblock * nbTestDof)
415  * M_globalCFE[faceIDinAdjacentElement]->M_wMeas[iQuadPt];
416 
417  }
418  }
419  }
420 
421  M_elementalVector.pushToGlobal (vec);
422  }
423 }
424 
425 
426 } // Namespace ExpressionAssembly
427 
428 } // Namespace LifeV
429 
430 #endif