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