MFEM  v4.5.2
Finite element discretization library
integ_algoim.cpp
Go to the documentation of this file.
1 #include "integ_algoim.hpp"
2 
3 
4 #ifdef MFEM_USE_ALGOIM
5 
6 
7 namespace mfem
8 {
9 
12  const Vector &lsfun)
13 {
14  int_order=o;
15  vir=nullptr;
16  sir=nullptr;
17 
18  if (el.GetGeomType()==Geometry::Type::SQUARE)
19  {
21  }
22  else if (el.GetGeomType()==Geometry::Type::CUBE)
23  {
24  pe=new H1Pos_HexahedronElement(el.GetOrder());
25  }
26  else
27  {
28  MFEM_ABORT("Currently MFEM + Algoim supports only quads and hexes.");
29  }
30 
31  //change the basis of the level-set function
32  //from Lagrangian to Bernstein (positive)
33  lsvec.SetSize(pe->GetDof());
34  DenseMatrix T(pe->GetDof());
35  pe->Project(el,trans,T);
36  T.Mult(lsfun,lsvec);
37 }
38 
39 
41 {
42  if (vir!=nullptr) {return vir;}
43 
44  const int dim=pe->GetDim();
45  int np1d=int_order/2+1;
46  if (dim==2)
47  {
48  LevelSet2D ls(pe,lsvec);
49  auto q = Algoim::quadGen<2>(ls,Algoim::BoundingBox<double,2>(0.0,1.0),
50  -1, -1, np1d);
51 
52  vir=new IntegrationRule(q.nodes.size());
53  vir->SetOrder(int_order);
54  for (size_t i=0; i<q.nodes.size(); i++)
55  {
56  IntegrationPoint& ip=vir->IntPoint(i);
57  ip.Set2w(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].w);
58  }
59  }
60  else
61  {
62  LevelSet3D ls(pe,lsvec);
63  auto q = Algoim::quadGen<3>(ls,Algoim::BoundingBox<double,3>(0.0,1.0),
64  -1, -1, np1d);
65 
66  vir=new IntegrationRule(q.nodes.size());
67  vir->SetOrder(int_order);
68  for (size_t i=0; i<q.nodes.size(); i++)
69  {
70  IntegrationPoint& ip=vir->IntPoint(i);
71  ip.Set(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].x(2),q.nodes[i].w);
72  }
73  }
74 
75  return vir;
76 }
77 
79 {
80  if (sir!=nullptr) {return sir;}
81 
82  int np1d=int_order/2+1;
83  const int dim=pe->GetDim();
84  if (dim==2)
85  {
86  LevelSet2D ls(pe,lsvec);
87  auto q = Algoim::quadGen<2>(ls,Algoim::BoundingBox<double,2>(0.0,1.0),
88  2, -1, np1d);
89 
90  sir=new IntegrationRule(q.nodes.size());
91  sir->SetOrder(int_order);
92  for (size_t i=0; i<q.nodes.size(); i++)
93  {
94  IntegrationPoint& ip=sir->IntPoint(i);
95  ip.Set2w(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].w);
96  }
97  }
98  else
99  {
100  LevelSet3D ls(pe,lsvec);
101  auto q = Algoim::quadGen<3>(ls,Algoim::BoundingBox<double,3>(0.0,1.0),
102  3, -1, np1d);
103 
104  sir=new IntegrationRule(q.nodes.size());
105  sir->SetOrder(int_order);
106  for (size_t i=0; i<q.nodes.size(); i++)
107  {
108  IntegrationPoint& ip=sir->IntPoint(i);
109  ip.Set(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].x(2),q.nodes[i].w);
110  }
111  }
112 
113  return sir;
114 }
115 
116 }
117 
118 #endif
119 
120 
Abstract class for all finite elements.
Definition: fe_base.hpp:232
void Set(const double *p, const int dim)
Definition: intrules.hpp:37
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()
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
const IntegrationRule * GetVolumeIntegrationRule()
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:319
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:310
Arbitrary order H1 elements in 3D utilizing the Bernstein basis on a cube.
Definition: fe_pos.hpp:159
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:322
void SetOrder(const int order)
Sets the order of the integration rule. This is only for keeping order information, it does not alter any data in the IntegrationRule.
Definition: intrules.hpp:244
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition: fe_pos.cpp:24
Class for integration point with weight.
Definition: intrules.hpp:25
AlgoimIntegrationRule(int o, const FiniteElement &el, ElementTransformation &trans, const Vector &lsfun)
int dim
Definition: ex24.cpp:53
Arbitrary order H1 elements in 2D utilizing the Bernstein basis on a square.
Definition: fe_pos.hpp:140
Vector data type.
Definition: vector.hpp:60
void Set2w(const double x1, const double x2, const double w)
Definition: intrules.hpp:75
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