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