LifeV
LevelSetQRAdapter.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 LifeV.
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 License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  @file
29  @brief Adapter for the quadrature on a given level set.
30 
31  This class is intended to be used in conjonction with the ETA framework
32  (this module must be available to use this class). Based on the values
33  of a level function, it computes (via a piecewise linear reconstruction
34  of the interface) a quadrature rule which is adapted to the 0 level set.
35 
36  @author Samuel Quinodoz <samuel.quinodoz@epfl.ch>
37  @date 15 Dec 2011
38 
39  */
40 
41 
42 #ifndef LEVELSETQRADAPTER_H
43 #define LEVELSETQRADAPTER_H 1
44 
45 #include <lifev/core/LifeV.hpp>
46 
47 #ifdef LIFEV_HAS_ETA
48 
49 #include <lifev/eta/fem/ETCurrentFE.hpp>
50 #include <lifev/core/fem/GeometricMap.hpp>
51 #include <lifev/core/fem/QuadratureRule.hpp>
52 
53 #include <lifev/eta/fem/QRAdapterBase.hpp>
54 
55 #include <boost/shared_ptr.hpp>
56 
57 namespace LifeV
58 {
59 
60 
61 /*!
62  The class LevelSetQRAdapter
63 
64  P1 reconstruction of the interface
65  and quadrature splitting
66 
67  3D only for the moment and P1
68 
69  */
70 template< typename FESpaceType, typename VectorType >
71 class LevelSetQRAdapter : public QRAdapterBase< LevelSetQRAdapter<FESpaceType, VectorType> >
72 {
73 public:
74 
75  //! @name Public Types
76  //@{
77 
78  typedef std::shared_ptr<FESpaceType> FESpaceType_Ptr;
79 
80  typedef QRAdapterBase< LevelSetQRAdapter<FESpaceType, VectorType> > base_Type;
81 
82  //@}
83 
84 
85  //! @name Constructor & Destructor
86  //@{
87 
88  //! Constructor with both FESpace and level set values
89  LevelSetQRAdapter (FESpaceType_Ptr fespace, const VectorType& vect, const QuadratureRule& qr);
90 
91  //! Copy constructor
92  LevelSetQRAdapter (const LevelSetQRAdapter<FESpaceType, VectorType>& lsqra);
93 
94  //! Simple destructor
95  ~LevelSetQRAdapter() {}
96 
97  //@}
98 
99 
100  //! @name Operators
101  //@{
102 
103  //@}
104 
105 
106  //! @name Methods
107  //@{
108 
109  //! Update this structure with the current element
110  /*!
111  This checks if the element is crossed by the interface
112  and if so, computes the adapted quadrature rule.
113  */
114  void update (UInt elementID);
115 
116 
117  //@}
118 
119 
120  //! @name Set Methods
121  //@{
122 
123  //@}
124 
125 
126  //! @name Get Methods
127  //@{
128 
129  //! Is the current element crossed by the interface
130  bool isAdaptedElement() const
131  {
132  return M_isAdaptedElement;
133  }
134 
135  //! Getter for the non-adapted quadrature
136  const QuadratureRule& standardQR() const
137  {
138  return M_qr;
139  }
140 
141  //! Getter for the adapted quadrature
142  const QuadratureRule& adaptedQR() const
143  {
144  return *M_adaptedQR;
145  }
146 
147 
148  //@}
149 
150 private:
151 
152  //! @name Private Methods
153  //@{
154 
155  /*! No default constructor as there is no default
156  constructor in the ETCurrentFE
157  */
158  LevelSetQRAdapter();
159 
160  //@}
161 
162  // Finite element space for the level set values
163  FESpaceType_Ptr M_lsFESpace;
164 
165  // Value of the level set
166  VectorType M_lsValue;
167 
168  // Quadrature rule
169  QuadratureRule M_qr;
170 
171  // CurrentFE for the level set
172 
173  ETCurrentFE<FESpaceType::space_dim, 1> M_currentFE;
174 
175  // Boolean indicating if the element is crossed by the interface
176  bool M_isAdaptedElement;
177 
178  // Adapted Quadrature Rule
179  std::shared_ptr<QuadratureRule> M_adaptedQR;
180 
181 };
182 
183 template<typename FESpaceType, typename VectorType>
184 LevelSetQRAdapter<FESpaceType, VectorType>
185 adapt (std::shared_ptr<FESpaceType> fespace, const VectorType& vector, const QuadratureRule& qr)
186 {
187  return LevelSetQRAdapter<FESpaceType, VectorType> (fespace, vector, qr);
188 }
189 
190 
191 template< typename FESpaceType, typename VectorType >
192 LevelSetQRAdapter<FESpaceType, VectorType>::
193 LevelSetQRAdapter (FESpaceType_Ptr fespace, const VectorType& vect, const QuadratureRule& qr)
194  :
195  base_Type(),
196  M_lsFESpace (fespace),
197  M_lsValue (vect, Repeated),
198  M_qr (qr),
199  M_currentFE (fespace->refFE(), getGeometricMap (*fespace->mesh() ), qr),
200  M_isAdaptedElement (false),
201  M_adaptedQR (new QuadratureRule (qr) )
202 {
203 
204  ASSERT ( FESpaceType::field_dim == 1, "Quadrature adaptation only for scalar fields!");
205  ASSERT ( FESpaceType::space_dim == 3, "Quadrature adaptation only 3D cases for the moment!");
206 }
207 
208 template< typename FESpaceType, typename VectorType >
209 LevelSetQRAdapter<FESpaceType, VectorType>::
210 LevelSetQRAdapter (const LevelSetQRAdapter<FESpaceType, VectorType>& lsqra)
211  :
212  base_Type(),
213  M_lsFESpace ( lsqra.M_lsFESpace ),
214  M_lsValue ( lsqra.M_lsValue ),
215  M_qr ( lsqra.M_qr ),
216  M_currentFE ( lsqra.M_currentFE ),
217  M_isAdaptedElement ( lsqra.M_isAdaptedElement ),
218  M_adaptedQR ( new QuadratureRule (*lsqra.M_adaptedQR) )
219 {}
220 
221 
222 template< typename FESpaceType, typename VectorType >
223 void
224 LevelSetQRAdapter<FESpaceType, VectorType>::
225 update (UInt elementID)
226 {
227  //! Check that the level set values are Repeated
228  ASSERT ( M_lsValue.mapType() == Repeated, "Internal error: level values are Unique!");
229  ASSERT ( M_lsFESpace != 0, "Internal error: empty pointer to the Space");
230 
231  // First, check the values of the level set and detect sign changes
232 
233  std::vector<UInt> localPositive;
234  std::vector<UInt> localNegative;
235  std::vector<Real> localLSValue (4, 0.0);
236 
237  for (UInt iDof (0); iDof < 4; ++iDof)
238  {
239  UInt myGlobalID ( M_lsFESpace->dof().localToGlobalMap (elementID, iDof) );
240  localLSValue[iDof] = M_lsValue[myGlobalID];
241 
242  if (localLSValue[iDof] >= 0)
243  {
244  localPositive.push_back (iDof);
245  }
246  else
247  {
248  localNegative.push_back (iDof);
249  }
250  }
251 
252  // Check if adaptation is needed
253  if ( localPositive.empty() || localNegative.empty() )
254  {
255  // No adaptation
256  M_isAdaptedElement = false;
257  *M_adaptedQR = M_qr;
258  }
259  else
260  {
261  // Adaptation
262  M_isAdaptedElement = true;
263 
264  // Define the reference tetrahedra
265  std::vector< VectorSmall<3> > refPts (4);
266  refPts[1][0] = 1;
267  refPts[2][1] = 1;
268  refPts[3][2] = 1;
269 
270  // Define the CurrentFE that will be used to map the
271  // quadrature from its reference shape to the part
272  // of the tetrahedra
273  ETCurrentFE<3, 1> QRMapper (feTetraP1, getGeometricMap (*M_lsFESpace->mesh() ), M_qr);
274 
275  // Reset the adapted quadrature
276  M_adaptedQR.reset (new QuadratureRule);
277  M_adaptedQR->setDimensionShape (3, TETRA);
278 
279  // Constant
280  UInt QuadratureSize (M_qr.nbQuadPt() );
281 
282  // If there's only one negative value
283  // the tetrahedra has to be cut into 4
284  // pieces
285  if (localPositive.size() == 1)
286  {
287  // Localize the intersections
288  std::vector< VectorSmall<3> > intersections (3);
289 
290  for (UInt i (0); i < 3; ++i)
291  {
292  UInt pos (localPositive[0]);
293  UInt neg (localNegative[i]);
294  Real mu ( localLSValue[pos] / (localLSValue[pos] - localLSValue[neg]) );
295  VectorSmall<3> I ( refPts[pos] + mu * ( refPts[neg] - refPts[pos] ) );
296  intersections[i] = I;
297  }
298 
299 
300  // Build the first tetra
301  std::vector< VectorSmall<3> > tetra1 (4);
302  tetra1[0] = refPts[localPositive[0]];
303  tetra1[1] = intersections[0];
304  tetra1[2] = intersections[1];
305  tetra1[3] = intersections[2];
306 
307  // Map the QR
308  QRMapper.update ( tetra1, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
309 
310  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
311  {
312  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
313  QRMapper.quadNode (iQuadPt, 1),
314  QRMapper.quadNode (iQuadPt, 2),
315  std::abs (QRMapper.wDet (iQuadPt) ) ) );
316  }
317 
318  // Build the second tetra
319  std::vector< VectorSmall<3> > tetra2 (4);
320  tetra2[0] = intersections[0];
321  tetra2[1] = intersections[1];
322  tetra2[2] = intersections[2];
323  tetra2[3] = refPts[localNegative[0]];
324 
325  // Map the QR
326  QRMapper.update ( tetra2, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
327 
328  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
329  {
330  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
331  QRMapper.quadNode (iQuadPt, 1),
332  QRMapper.quadNode (iQuadPt, 2),
333  std::abs (QRMapper.wDet (iQuadPt) ) ) );
334  }
335 
336  // Build the third tetra
337  std::vector< VectorSmall<3> > tetra3 (4);
338  tetra3[0] = refPts[localNegative[0]];
339  tetra3[1] = refPts[localNegative[1]];
340  tetra3[2] = refPts[localNegative[2]];
341  tetra3[3] = intersections[2];
342 
343  // Map the QR
344  QRMapper.update ( tetra3, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
345 
346  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
347  {
348  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
349  QRMapper.quadNode (iQuadPt, 1),
350  QRMapper.quadNode (iQuadPt, 2),
351  std::abs (QRMapper.wDet (iQuadPt) ) ) );
352  }
353 
354 
355  // Build the fourth tetra
356  std::vector< VectorSmall<3> > tetra4 (4);
357  tetra4[0] = refPts[localNegative[0]];
358  tetra4[1] = refPts[localNegative[1]];
359  tetra4[2] = intersections[1];
360  tetra4[3] = intersections[2];
361 
362  // Map the QR
363  QRMapper.update ( tetra4, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
364 
365  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
366  {
367  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
368  QRMapper.quadNode (iQuadPt, 1),
369  QRMapper.quadNode (iQuadPt, 2),
370  std::abs (QRMapper.wDet (iQuadPt) ) ) );
371  }
372 
373  }
374  else if ( localPositive.size() == 2)
375  {
376 
377  // Localize the intersections
378  std::vector< VectorSmall<3> > intersections (4);
379 
380  for (UInt i (0); i < 2; ++i)
381  {
382  for (UInt j (0); j < 2; ++j)
383  {
384  UInt pos (localPositive[i]);
385  UInt neg (localNegative[j]);
386  Real mu ( localLSValue[pos] / (localLSValue[pos] - localLSValue[neg]) );
387  VectorSmall<3> I ( refPts[pos] + mu * ( refPts[neg] - refPts[pos] ) );
388  intersections[i + 2 * j] = I;
389  }
390  }
391 
392 
393  // Build the first tetra
394  std::vector< VectorSmall<3> > tetra1 (4);
395  tetra1[0] = refPts[localPositive[0]];
396  tetra1[1] = intersections[0];
397  tetra1[2] = intersections[1];
398  tetra1[3] = intersections[2];
399 
400  // Map the QR
401  QRMapper.update ( tetra1, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
402 
403  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
404  {
405  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
406  QRMapper.quadNode (iQuadPt, 1),
407  QRMapper.quadNode (iQuadPt, 2),
408  std::abs (QRMapper.wDet (iQuadPt) ) ) );
409  }
410 
411  // Build the second tetra
412  std::vector< VectorSmall<3> > tetra2 (4);
413  tetra2[0] = intersections[1];
414  tetra2[1] = intersections[2];
415  tetra2[2] = intersections[3];
416  tetra2[3] = refPts[localPositive[0]];
417 
418  // Map the QR
419  QRMapper.update ( tetra2, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
420 
421  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
422  {
423  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
424  QRMapper.quadNode (iQuadPt, 1),
425  QRMapper.quadNode (iQuadPt, 2),
426  std::abs (QRMapper.wDet (iQuadPt) ) ) );
427  }
428 
429  // Build the third tetra
430  std::vector< VectorSmall<3> > tetra3 (4);
431  tetra3[0] = refPts[localPositive[0]];
432  tetra3[1] = refPts[localPositive[1]];
433  tetra3[2] = intersections[1];
434  tetra3[3] = intersections[3];
435 
436  // Map the QR
437  QRMapper.update ( tetra3, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
438 
439  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
440  {
441  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
442  QRMapper.quadNode (iQuadPt, 1),
443  QRMapper.quadNode (iQuadPt, 2),
444  std::abs (QRMapper.wDet (iQuadPt) ) ) );
445  }
446 
447  // Build the fourth tetra
448  std::vector< VectorSmall<3> > tetra4 (4);
449  tetra4[0] = refPts[localNegative[0]];
450  tetra4[1] = intersections[0];
451  tetra4[2] = intersections[1];
452  tetra4[3] = intersections[2];
453 
454  // Map the QR
455  QRMapper.update ( tetra4, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
456 
457  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
458  {
459  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
460  QRMapper.quadNode (iQuadPt, 1),
461  QRMapper.quadNode (iQuadPt, 2),
462  std::abs (QRMapper.wDet (iQuadPt) ) ) );
463  }
464 
465  // Build the fifth tetra
466  std::vector< VectorSmall<3> > tetra5 (4);
467  tetra5[0] = refPts[localNegative[0]];
468  tetra5[1] = intersections[1];
469  tetra5[2] = intersections[2];
470  tetra5[3] = intersections[3];
471 
472  // Map the QR
473  QRMapper.update ( tetra5, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
474 
475  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
476  {
477  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
478  QRMapper.quadNode (iQuadPt, 1),
479  QRMapper.quadNode (iQuadPt, 2),
480  std::abs (QRMapper.wDet (iQuadPt) ) ) );
481  }
482 
483  // Build the sixth tetra
484  std::vector< VectorSmall<3> > tetra6 (4);
485  tetra6[0] = refPts[localNegative[0]];
486  tetra6[1] = refPts[localNegative[1]];
487  tetra6[2] = intersections[2];
488  tetra6[3] = intersections[3];
489 
490  // Map the QR
491  QRMapper.update ( tetra6, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
492 
493  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
494  {
495  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
496  QRMapper.quadNode (iQuadPt, 1),
497  QRMapper.quadNode (iQuadPt, 2),
498  std::abs (QRMapper.wDet (iQuadPt) ) ) );
499  }
500 
501 
502 
503  }
504  else
505  {
506  ASSERT ( localPositive.size() == 3, "Internal inconsistency");
507 
508  // Localize the intersections
509  std::vector< VectorSmall<3> > intersections (3);
510 
511  for (UInt i (0); i < 3; ++i)
512  {
513  UInt pos (localPositive[i]);
514  UInt neg (localNegative[0]);
515  Real mu ( localLSValue[pos] / (localLSValue[pos] - localLSValue[neg]) );
516  VectorSmall<3> I ( refPts[pos] + mu * ( refPts[neg] - refPts[pos] ) );
517  intersections[i] = I;
518  }
519 
520 
521  // Build the first tetra
522  std::vector< VectorSmall<3> > tetra1 (4);
523  tetra1[0] = refPts[localNegative[0]];
524  tetra1[1] = intersections[0];
525  tetra1[2] = intersections[1];
526  tetra1[3] = intersections[2];
527 
528  // Map the QR
529  QRMapper.update ( tetra1, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
530 
531  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
532  {
533  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
534  QRMapper.quadNode (iQuadPt, 1),
535  QRMapper.quadNode (iQuadPt, 2),
536  std::abs (QRMapper.wDet (iQuadPt) ) ) );
537  }
538 
539 
540  // Build the second tetra
541  std::vector< VectorSmall<3> > tetra2 (4);
542  tetra2[0] = intersections[0];
543  tetra2[1] = intersections[1];
544  tetra2[2] = intersections[2];
545  tetra2[3] = refPts[localPositive[0]];
546 
547  // Map the QR
548  QRMapper.update ( tetra2, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
549 
550  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
551  {
552  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
553  QRMapper.quadNode (iQuadPt, 1),
554  QRMapper.quadNode (iQuadPt, 2),
555  std::abs (QRMapper.wDet (iQuadPt) ) ) );
556  }
557 
558  // Build the third tetra
559  std::vector< VectorSmall<3> > tetra3 (4);
560  tetra3[0] = refPts[localPositive[0]];
561  tetra3[1] = refPts[localPositive[1]];
562  tetra3[2] = refPts[localPositive[2]];
563  tetra3[3] = intersections[2];
564 
565  // Map the QR
566  QRMapper.update ( tetra3, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
567 
568  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
569  {
570  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
571  QRMapper.quadNode (iQuadPt, 1),
572  QRMapper.quadNode (iQuadPt, 2),
573  std::abs (QRMapper.wDet (iQuadPt) ) ) );
574  }
575 
576  // Build the fourth tetra
577  std::vector< VectorSmall<3> > tetra4 (4);
578  tetra4[0] = refPts[localPositive[0]];
579  tetra4[1] = refPts[localPositive[1]];
580  tetra4[2] = intersections[1];
581  tetra4[3] = intersections[2];
582 
583  // Map the QR
584  QRMapper.update ( tetra4, ET_UPDATE_QUAD_NODE | ET_UPDATE_WDET);
585 
586  for (UInt iQuadPt (0); iQuadPt < QuadratureSize; ++iQuadPt)
587  {
588  M_adaptedQR->addPoint (QuadraturePoint ( QRMapper.quadNode (iQuadPt, 0),
589  QRMapper.quadNode (iQuadPt, 1),
590  QRMapper.quadNode (iQuadPt, 2),
591  std::abs (QRMapper.wDet (iQuadPt) ) ) );
592  }
593  }
594  }
595 }
596 
597 
598 } // Namespace LifeV
599 
600 #endif /* LEVELSETQRADAPTER_H */
601 
602 #endif /* LIFEV_HAS_ETA */