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