LifeV
LevelSetBDQRAdapter.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  @author Samuel Quinodoz <samuel.quinodoz@epfl.ch>
32  @date 15 Dec 2011
33 
34  */
35 
36 #ifndef LEVELSETBDQRADAPTER_H
37 #define LEVELSETBDQRADAPTER_H 1
38 
39 #include <lifev/core/LifeV.hpp>
40 
41 #include <lifev/eta/fem/QuadratureBoundary.hpp>
42 #include <lifev/eta/fem/ETCurrentFE.hpp>
43 
44 #include <boost/shared_ptr.hpp>
45 
46 namespace LifeV
47 {
48 
49 
50 // Cannot inherit from the volume version because it would
51 // yield the volume weight, not the surface weight.
52 
53 template< typename FESpaceType, typename VectorType >
54 class LevelSetBDQRAdapter
55 {
56 public:
57 
58  //! @name Public Types
59  //@{
60 
61  typedef std::shared_ptr<FESpaceType> FESpaceType_Ptr;
62 
63  //@}
64 
65 
66  //! @name Constructor & Destructor
67  //@{
68 
69  LevelSetBDQRAdapter ( FESpaceType_Ptr fespace, const VectorType& vect, const QuadratureBoundary& qrbd)
70  : M_lsFESpace (fespace),
71  M_lsValue (vect, Repeated),
72  M_qrBd (qrbd),
73  M_currentFE (fespace->refFE(), getGeometricMap (*fespace->mesh() ), qrbd.qr (0) ),
74  M_isAdaptedElement (false),
75  M_adaptedQrBd (new QuadratureBoundary (qrbd) )
76  {};
77 
78  LevelSetBDQRAdapter ( const LevelSetBDQRAdapter& lsbdqra)
79  : M_lsFESpace ( lsbdqra.M_lsFESpace),
80  M_lsValue ( lsbdqra.M_lsValue),
81  M_qrBd ( lsbdqra.M_qrBd),
82  M_currentFE ( lsbdqra.M_currentFE),
83  M_isAdaptedElement ( lsbdqra.M_isAdaptedElement),
84  M_adaptedQrBd (new QuadratureBoundary (lsbdqra.M_qrBd) )
85  {};
86 
87  ~LevelSetBDQRAdapter() {}
88 
89  //@}
90 
91 
92  //! @name Methods
93  //@{
94 
95  void update (UInt elementID, UInt localFaceID);
96 
97  //@}
98 
99 
100  //! @name Set Methods
101  //@{
102 
103  //@}
104 
105 
106  //! @name Get Methods
107  //@{
108 
109  QuadratureRule adaptedBdQR ( const UInt i) const
110  {
111  if (this->M_isAdaptedElement)
112  {
113  return M_adaptedQrBd->qr (i);
114  }
115  return M_qrBd.qr (i);
116  }
117 
118  //@}
119 
120 private:
121 
122  LevelSetBDQRAdapter();
123 
124 
125  FESpaceType_Ptr M_lsFESpace;
126 
127  VectorType M_lsValue;
128 
129  QuadratureBoundary M_qrBd;
130 
131  ETCurrentFE<FESpaceType::space_dim, 1> M_currentFE;
132 
133  bool M_isAdaptedElement;
134 
135  std::shared_ptr<QuadratureBoundary> M_adaptedQrBd;
136 
137 };
138 
139 template<typename FESpaceType, typename VectorType>
140 LevelSetBDQRAdapter<FESpaceType, VectorType>
141 adapt (std::shared_ptr<FESpaceType> fespace, const VectorType& vector, const QuadratureBoundary& qrbd)
142 {
143  return LevelSetBDQRAdapter<FESpaceType, VectorType> (fespace, vector, qrbd);
144 }
145 
146 
147 template<typename FESpaceType, typename VectorType>
148 void
149 LevelSetBDQRAdapter<FESpaceType, VectorType>::
150 update (UInt elementID, UInt localFaceID)
151 {
152  // They are ordered so that the "mapping" in
153  // QuadratureBoundary is right.
154 
155  std::vector<UInt> myLID (3);
156 
157  if (localFaceID == 0)
158  {
159  myLID[0] = 0;
160  myLID[1] = 1;
161  myLID[2] = 2;
162  }
163  else if (localFaceID == 1)
164  {
165  myLID[0] = 0;
166  myLID[1] = 1;
167  myLID[2] = 3;
168  }
169  else if (localFaceID == 2)
170  {
171  myLID[0] = 3;
172  myLID[1] = 1;
173  myLID[2] = 2;
174  }
175  else
176  {
177  myLID[0] = 0;
178  myLID[1] = 3;
179  myLID[2] = 2;
180  }
181 
182  // Compute the global IDs
183  std::vector<UInt> myGID (3);
184 
185  for (UInt i (0); i < 3; ++i)
186  {
187  myGID[i] = M_lsFESpace->dof().localToGlobalMap (elementID, myLID[i]);
188  }
189 
190  // Check the level set values
191  std::vector<UInt> localPositive;
192  std::vector<UInt> localNegative;
193  std::vector<Real> localPositiveValue;
194  std::vector<Real> localNegativeValue;
195 
196  for (UInt i (0); i < 3; ++i)
197  {
198  if (M_lsValue[myGID[i]] >= 0)
199  {
200  localPositive.push_back (i);
201  localPositiveValue.push_back (M_lsValue[myGID[i]]);
202  }
203  else
204  {
205  localNegative.push_back (i);
206  localNegativeValue.push_back (M_lsValue[myGID[i]]);
207  }
208  }
209 
210  // Check if the element is actually cut by the interface.
211  if (localPositive.empty() || localNegative.empty() )
212  {
213  M_isAdaptedElement = false;
214  *M_adaptedQrBd = M_qrBd;
215  }
216  else
217  {
218 
219  M_isAdaptedElement = true;
220 
221  std::vector<std::vector<Real> >refCoor (3, std::vector<Real> (2) );
222  refCoor[0][0] = 0.0;
223  refCoor[0][1] = 0.0;
224  refCoor[1][0] = 1.0;
225  refCoor[1][1] = 0.0;
226  refCoor[2][0] = 0.0;
227  refCoor[2][1] = 1.0;
228 
229  QuadratureRule qrRef ("none", TRIANGLE, 3, 0, 0);
230 
231  // Disinguish the two cases
232  if (localPositive.size() == 1)
233  {
234  // Do nothing
235  }
236  else
237  {
238  // Exchange positive and negative values
239  localPositive.swap (localNegative);
240  localPositiveValue.swap (localNegativeValue);
241  }
242 
243 
244  // Two negative case
245  Real lambda1 = -localPositiveValue[0] / (localNegativeValue[0] - localPositiveValue[0]);
246 
247  std::vector < Real > I1 (2);
248  I1[0] = refCoor[localPositive[0]][0] + lambda1 * (refCoor[localNegative[0]][0] - refCoor[localPositive[0]][0]);
249  I1[1] = refCoor[localPositive[0]][1] + lambda1 * (refCoor[localNegative[0]][1] - refCoor[localPositive[0]][1]);
250 
251  Real lambda2 = -localPositiveValue[0] / (localNegativeValue[1] - localPositiveValue[0]);
252 
253  std::vector < Real > I2 (2);
254  I2[0] = refCoor[localPositive[0]][0] + lambda2 * (refCoor[localNegative[1]][0] - refCoor[localPositive[0]][0]);
255  I2[1] = refCoor[localPositive[0]][1] + lambda2 * (refCoor[localNegative[1]][1] - refCoor[localPositive[0]][1]);
256 
257  // First triangle
258  {
259  std::vector<Real> P0 (2);
260  P0[0] = refCoor[localPositive[0]][0];
261  P0[1] = refCoor[localPositive[0]][1];
262 
263  std::vector<Real> v1 (2);
264  v1[0] = I1[0] - refCoor[localPositive[0]][0];
265  v1[1] = I1[1] - refCoor[localPositive[0]][1];
266 
267  std::vector<Real> v2 (2);
268  v2[0] = I2[0] - refCoor[localPositive[0]][0];
269  v2[1] = I2[1] - refCoor[localPositive[0]][1];
270 
271  Real wRel (std::abs (v1[0]*v2[1] - v1[1]*v2[0]) );
272 
273  for (UInt iq (0); iq < M_qrBd.qr (0).nbQuadPt(); ++iq)
274  {
275  Real x (P0[0] + M_qrBd.qr (0).quadPointCoor (iq, 0) *v1[0] + M_qrBd.qr (0).quadPointCoor (iq, 1) *v2[0]);
276  Real y (P0[1] + M_qrBd.qr (0).quadPointCoor (iq, 0) *v1[1] + M_qrBd.qr (0).quadPointCoor (iq, 1) *v2[1]);
277 
278  qrRef.addPoint (QuadraturePoint (x, y, M_qrBd.qr (0).weight (iq) *wRel) );
279  }
280  }
281 
282  // Second triangle
283  {
284  std::vector<Real> P0 (2);
285  P0[0] = refCoor[localNegative[0]][0];
286  P0[1] = refCoor[localNegative[0]][1];
287 
288  std::vector<Real> v1 (2);
289  v1[0] = I1[0] - refCoor[localNegative[0]][0];
290  v1[1] = I1[1] - refCoor[localNegative[0]][1];
291 
292  std::vector<Real> v2 (2);
293  v2[0] = I2[0] - refCoor[localNegative[0]][0];
294  v2[1] = I2[1] - refCoor[localNegative[0]][1];
295 
296  Real wRel (std::abs (v1[0]*v2[1] - v1[1]*v2[0]) );
297 
298  for (UInt iq (0); iq < M_qrBd.qr (0).nbQuadPt(); ++iq)
299  {
300  Real x (P0[0] + M_qrBd.qr (0).quadPointCoor (iq, 0) *v1[0] + M_qrBd.qr (0).quadPointCoor (iq, 1) *v2[0]);
301  Real y (P0[1] + M_qrBd.qr (0).quadPointCoor (iq, 0) *v1[1] + M_qrBd.qr (0).quadPointCoor (iq, 1) *v2[1]);
302 
303  qrRef.addPoint (QuadraturePoint (x, y, M_qrBd.qr (0).weight (iq) *wRel) );
304  }
305  }
306 
307  // Third triangle
308  {
309  std::vector<Real> P0 (2);
310  P0[0] = refCoor[localNegative[0]][0];
311  P0[1] = refCoor[localNegative[0]][1];
312 
313  std::vector<Real> v1 (2);
314  v1[0] = I2[0] - refCoor[localNegative[0]][0];
315  v1[1] = I2[1] - refCoor[localNegative[0]][1];
316 
317  std::vector<Real> v2 (2);
318  v2[0] = refCoor[localNegative[1]][0] - refCoor[localNegative[0]][0];
319  v2[1] = refCoor[localNegative[1]][1] - refCoor[localNegative[0]][1];
320 
321  Real wRel (std::abs (v1[0]*v2[1] - v1[1]*v2[0]) );
322 
323  for (UInt iq (0); iq < M_qrBd.qr (0).nbQuadPt(); ++iq)
324  {
325  Real x (P0[0] + M_qrBd.qr (0).quadPointCoor (iq, 0) *v1[0] + M_qrBd.qr (0).quadPointCoor (iq, 1) *v2[0]);
326  Real y (P0[1] + M_qrBd.qr (0).quadPointCoor (iq, 0) *v1[1] + M_qrBd.qr (0).quadPointCoor (iq, 1) *v2[1]);
327 
328  qrRef.addPoint (QuadraturePoint (x, y, M_qrBd.qr (0).weight (iq) *wRel) );
329  }
330  }
331 
332  // Call the make tretra
333  M_adaptedQrBd.reset (new QuadratureBoundary (buildTetraBDQR (qrRef) ) );
334  }
335 }
336 
337 
338 } // Namespace LifeV
339 
340 #endif /* LEVELSETBDQRADAPTER_H */