MFEM  v4.5.2
Finite element discretization library
integ_algoim.hpp
Go to the documentation of this file.
1 #ifndef INTEG_ALGOIM_HPP
2 #define INTEG_ALGOIM_HPP
3 
4 #include <mfem.hpp>
5 
6 #ifdef MFEM_USE_ALGOIM
7 #include <algoim_quad.hpp>
8 
9 namespace mfem
10 {
11 
12 
13 //define templated element bases
14 namespace TmplPoly_1D
15 {
16 
17 /// Templated version of CalcBinomTerms
18 template<typename float_type>
19 void CalcBinomTerms(const int p, const float_type x, const float_type y,
20  float_type* u)
21 {
22  if (p == 0)
23  {
24  u[0] = float_type(1.);
25  }
26  else
27  {
28  int i;
29  const int *b = Poly_1D::Binom(p);
30  float_type z = x;
31  for (i = 1; i < p; i++)
32  {
33  u[i] = b[i]*z;
34  z *= x;
35  }
36  u[p] = z;
37  z = y;
38  for (i--; i > 0; i--)
39  {
40  u[i] *= z;
41  z *= y;
42  }
43  u[0] = z;
44  }
45 }
46 
47 /// Templated version of CalcBinomTerms
48 template<typename float_type>
49 void CalcBinomTerms(const int p, const float_type x, const float_type y,
50  float_type* u, float_type* d)
51 {
52  if (p == 0)
53  {
54  u[0] = float_type(1.);
55  d[0] = float_type(0.);
56  }
57  else
58  {
59  int i;
60  const int *b = Poly_1D::Binom(p);
61  const float_type xpy = x + y, ptx = p*x;
62  float_type z = float_type(1.);
63 
64  for (i = 1; i < p; i++)
65  {
66  d[i] = b[i]*z*(i*xpy - ptx);
67  z *= x;
68  u[i] = b[i]*z;
69  }
70  d[p] = p*z;
71  u[p] = z*x;
72  z = float_type(1.);
73  for (i--; i > 0; i--)
74  {
75  d[i] *= z;
76  z *= y;
77  u[i] *= z;
78  }
79  d[0] = -p*z;
80  u[0] = z*y;
81  }
82 
83 }
84 
85 /// Templated evaluation of Bernstein basis
86 template <typename float_type>
87 void CalcBernstein(const int p, const float_type x, float_type *u)
88 {
89  CalcBinomTerms(p, x, 1. - x, u);
90 }
91 
92 
93 /// Templated evaluation of Bernstein basis
94 template <typename float_type>
95 void CalcBernstein(const int p, const float_type x,
96  float_type *u, float_type *d)
97 {
98  CalcBinomTerms(p, x, 1. - x, u, d);
99 }
100 
101 
102 }
103 
104 /// Construct volumetric and surface integration rules for a given element
105 /// using the Algoim library. The volume is define as the positive part of
106 /// a level-set function(LSF) (lsfun argument in the constructor). The surface
107 /// is defined as the zero level-set of the LSF.
109 {
110 public:
111 
112  /// Construct Algoim object using a finite element, its transformation
113  /// and level-set function defined over the element using Lagrangian
114  /// bases. The argument o provides the order of the of the 1D Gaussian
115  /// integration rule used for deriving the vol/surface integration rules.
116  AlgoimIntegrationRule(int o, const FiniteElement &el,
117  ElementTransformation &trans, const Vector &lsfun);
118 
119 
120  /// Destructor of the Algoim object
122  {
123  delete pe;
124  delete vir;
125  delete sir;
126  }
127 
128  /// Returns volumetric integration rule based on the provided
129  /// level-set function.
131 
132  /// Returns surface integration rule based on the provided
133  /// level-set function.
135 
136 
137 private:
138 
139  /// 3D level-set function object required by Algoim.
140  struct LevelSet3D
141  {
142  /// Constructor for 3D level-set function object required by Algoim.
143  LevelSet3D(PositiveTensorFiniteElement* el_, Vector& lsfun_)
144  : el(el_), lsfun(lsfun_) { }
145 
146  /// Returns the value of the LSF for point x.
147  template<typename T>
148  T operator() (const blitz::TinyVector<T,3>& x) const
149  {
150  int el_order=el->GetOrder();
151  T u1[el_order+1];
152  T u2[el_order+1];
153  T u3[el_order+1];
154  TmplPoly_1D::CalcBernstein(el_order, x[0], u1);
155  TmplPoly_1D::CalcBernstein(el_order, x[1], u2);
156  TmplPoly_1D::CalcBernstein(el_order, x[2], u3);
157 
158  const Array<int>& dof_map=el->GetDofMap();
159 
160  T res=T(0.0);
161  for (int oo = 0, kk = 0; kk <= el_order; kk++)
162  for (int jj = 0; jj <= el_order; jj++)
163  for (int ii = 0; ii <= el_order; ii++)
164  {
165  res=res-u1[ii]*u2[jj]*u3[kk]*lsfun(dof_map[oo++]);
166  }
167  return res;
168  }
169 
170  /// Returns the gradients of the LSF for point x.
171  template<typename T>
172  blitz::TinyVector<T,3> grad(const blitz::TinyVector<T,3>& x) const
173  {
174  int el_order=el->GetOrder();
175  T u1[el_order+1];
176  T u2[el_order+1];
177  T u3[el_order+1];
178  T d1[el_order+1];
179  T d2[el_order+1];
180  T d3[el_order+1];
181 
182  TmplPoly_1D::CalcBernstein(el_order,x[0], u1, d1);
183  TmplPoly_1D::CalcBernstein(el_order,x[1], u2, d2);
184  TmplPoly_1D::CalcBernstein(el_order,x[2], u3, d3);
185 
186  blitz::TinyVector<T,3> res(T(0.0),T(0.0),T(0.0));
187 
188  const Array<int>& dof_map=el->GetDofMap();
189 
190  for (int oo = 0, kk = 0; kk <= el_order; kk++)
191  for (int jj = 0; jj <= el_order; jj++)
192  for (int ii = 0; ii <= el_order; ii++)
193  {
194  res[0]=res[0]-d1[ii]*u2[jj]*u3[kk]*lsfun(dof_map[oo]);
195  res[1]=res[1]-u1[ii]*d2[jj]*u3[kk]*lsfun(dof_map[oo]);
196  res[2]=res[2]-u1[ii]*u2[jj]*d3[kk]*lsfun(dof_map[oo]);
197  oo++;
198  }
199 
200  return res;
201  }
202 
203  private:
204  PositiveTensorFiniteElement* el;
205  Vector& lsfun;
206  };
207 
208  /// 2D level-set function object required by Algoim.
209  struct LevelSet2D
210  {
211  /// Constructor for 2D level-set function object required by Algoim.
212  LevelSet2D(PositiveTensorFiniteElement* el_, Vector& lsfun_)
213  :el(el_), lsfun(lsfun_) { }
214 
215  /// Returns the value of the LSF for point x.
216  template<typename T>
217  T operator() (const blitz::TinyVector<T,2>& x) const
218  {
219  int el_order=el->GetOrder();
220  T u1[el_order+1];
221  T u2[el_order+1];
222  TmplPoly_1D::CalcBernstein(el_order, x[0], u1);
223  TmplPoly_1D::CalcBernstein(el_order, x[1], u2);
224 
225  const Array<int>& dof_map=el->GetDofMap();
226 
227  T res=T(0.0);
228 
229  for (int oo = 0, jj = 0; jj <= el_order; jj++)
230  for (int ii = 0; ii <= el_order; ii++)
231  {
232  res=res-u1[ii]*u2[jj]*lsfun(dof_map[oo++]);
233  }
234  return res;
235  }
236 
237  /// Returns the gradients of the LSF for point x.
238  template<typename T>
239  blitz::TinyVector<T,2> grad(const blitz::TinyVector<T,2>& x) const
240  {
241  int el_order=el->GetOrder();
242  T u1[el_order+1];
243  T u2[el_order+1];
244  T d1[el_order+1];
245  T d2[el_order+1];
246 
247  TmplPoly_1D::CalcBernstein(el_order,x[0], u1, d1);
248  TmplPoly_1D::CalcBernstein(el_order,x[1], u2, d2);
249 
250  blitz::TinyVector<T,2> res(T(0.0),T(0.0));
251 
252  const Array<int>& dof_map=el->GetDofMap();
253 
254  for (int oo = 0, jj = 0; jj <= el_order; jj++)
255  for (int ii = 0; ii <= el_order; ii++)
256  {
257  res[0]=res[0]-(d1[ii]*u2[jj])*lsfun(dof_map[oo]);
258  res[1]=res[1]-(u1[ii]*d2[jj])*lsfun(dof_map[oo]);
259  oo++;
260  }
261 
262  return res;
263  }
264 
265 
266  private:
267  PositiveTensorFiniteElement* el;
268  Vector& lsfun;
269  };
270 
271 
272  IntegrationRule* sir; // Surface integration rule. Owned.
273  IntegrationRule* vir; // Volumetric integration rule. Owned.
274  PositiveTensorFiniteElement *pe;
275  Vector lsvec; //level-set in Bernstein bases
276  int int_order; //integration order
277 };
278 
279 }
280 #endif
281 
282 #endif
Abstract class for all finite elements.
Definition: fe_base.hpp:232
void CalcBernstein(const int p, const float_type x, float_type *u)
Templated evaluation of Bernstein basis.
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
const IntegrationRule * GetSurfaceIntegrationRule()
const IntegrationRule * GetVolumeIntegrationRule()
double b
Definition: lissajous.cpp:42
void CalcBinomTerms(const int p, const float_type x, const float_type y, float_type *u)
Templated version of CalcBinomTerms.
static const int * Binom(const int p)
Get a pointer to an array containing the binomial coefficients "p choose k" for k=0,...,p for the given p.
Definition: fe_base.cpp:1930
~AlgoimIntegrationRule()
Destructor of the Algoim object.
AlgoimIntegrationRule(int o, const FiniteElement &el, ElementTransformation &trans, const Vector &lsfun)
Vector data type.
Definition: vector.hpp:60
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:326