MFEM  v4.6.0
Finite element discretization library
lsf_integral.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // ---------------------------------------------------------------------
13 // Miniapp: Integral of implicit domains defined by a level set function
14 // ---------------------------------------------------------------------
15 // The miniapp demonstrates an interface to the Algoim library for computing
16 // volumetric and surface integrals over domains and surfaces defined implicitly
17 // by a level set function. The miniapp requires MFEM to be built with
18 // Blitz and Algoim libraries (see INSTALL).
19 //
20 // Compile with: make lsf_integral
21 //
22 // Sample runs:
23 //
24 // Evaluates surface and volumetric integral for a circle with radius 1
25 // lsf_integral -m ../../data/star-q3.mesh
26 //
27 // Evaluates surface and volumetric integral for a level set defined
28 // by y=0.5-(0.1*sin(3*pi*x+pi/2))
29 // lsf_integral -ls 2 -m ../../data/inline-quad.mesh -rs 2 -o 3 -ao 3
30 
31 #include "mfem.hpp"
32 #include "integ_algoim.hpp"
33 
34 using namespace mfem;
35 using namespace std;
36 
37 //Level set function for sphere in 3D and circle in 2D
38 double sphere_ls(const Vector &x)
39 {
40  double r2= x*x;
41  return -sqrt(r2)+1.0;//the radius is 1.0
42 }
43 
44 //Level set function for a sinusoidal wave.
45 //Resulting zero isocontour is at y=0.5-(0.1*sin(3*pi*x+pi/2))
46 double sinusoidal_ls(const Vector &x)
47 {
48  double a1 = 20., a2 = 2., a3 = 3.;
49  return tanh(a1*(x(1)-0.5) + a2*sin(a3*(x(0)-0.5)*M_PI));
50 }
51 
52 int main(int argc, char *argv[])
53 {
54  // Parse command-line options
55  const char *mesh_file = "../../data/star-q3.mesh";
56  int ser_ref_levels = 1;
57  int order = 2;
58  int iorder = 2; //MFEM integration integration points
59  int aorder = 2; //Algoim integration integration points
60  bool visualization = true;
61  int print_level = 0;
62  int ls_type = 1;
63 
64  OptionsParser args(argc, argv);
65  args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
66  args.AddOption(&ser_ref_levels,
67  "-rs",
68  "--refine-serial",
69  "Number of times to refine the mesh uniformly in serial.");
70  args.AddOption(&order,
71  "-o",
72  "--order",
73  "Order (degree) of the finite elements.");
74  args.AddOption(&iorder,
75  "-io",
76  "--iorder",
77  "MFEM Integration order.");
78  args.AddOption(&aorder,
79  "-ao",
80  "--aorder",
81  "Algoim Integration order.");
82  args.AddOption(&ls_type,
83  "-ls",
84  "--ls-type",
85  "Level set type: 1: circle, 2 sinusoidal wave");
86  args.AddOption(&visualization,
87  "-vis",
88  "--visualization",
89  "-no-vis",
90  "--no-visualization",
91  "Enable or disable GLVis visualization.");
92  args.AddOption((&print_level), "-prt", "--print-level", "Print level.");
93  args.Parse();
94  if (!args.Good())
95  {
96  args.PrintUsage(std::cout);
97  return 1;
98  }
99  args.PrintOptions(std::cout);
100 
101  // Read the (serial) mesh from the given mesh file.
102  Mesh *mesh = new Mesh(mesh_file, 1, 1);
103  int dim = mesh->Dimension();
104  // Refine the mesh in serial to increase the resolution. In this example
105  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
106  // a command-line parameter.
107  for (int lev = 0; lev < ser_ref_levels; lev++)
108  {
109  mesh->UniformRefinement();
110  }
111 
112  // Define the finite element space for the level-set function.
113  H1_FECollection fec(order, dim);
114  FiniteElementSpace fespace(mesh, &fec, 1, Ordering::byVDIM);
115  int glob_size = fespace.GetTrueVSize();
116  std::cout << "Number of finite element unknowns: " << glob_size << std::endl;
117 
118  // Define the level set grid function
119  GridFunction x(&fespace);
120 
121  // Define the level set coefficient
122  Coefficient *ls_coeff = nullptr;
123  if (ls_type == 1)
124  {
125  ls_coeff=new FunctionCoefficient(sphere_ls);
126  }
127  else if (ls_type == 2)
128  {
129  ls_coeff=new FunctionCoefficient(sinusoidal_ls);
130  }
131  else
132  {
133  MFEM_ABORT("Level set coefficient not defined");
134  }
135 
136  // Project the coefficient onto the LS grid function
137  x.ProjectCoefficient(*ls_coeff);
138 
139  if (visualization)
140  {
141  osockstream sock(19916, "localhost");
142  sock << "solution\n";
143  mesh->Print(sock);
144  x.Save(sock);
145  sock.send();
146  sock << "window_title 'Level set'\n"
147  << "window_geometry "
148  << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
149  << "keys jRmmclA" << endl;
150  }
151 
152  // Exact volume and area
153  double exact_volume = -10, exact_area = -10;
154  if (ls_type == 1)
155  {
156  if (strncmp(mesh_file,"../../data/star-q3.mesh",100) == 0)
157  {
158  exact_volume = M_PI;
159  exact_area = M_PI*2;
160  }
161  else if (strncmp(mesh_file, "../../data/inline-quad.mesh",100) == 0)
162  {
163  exact_volume = M_PI/4;
164  exact_area = M_PI/2;
165  }
166  }
167  else if (ls_type == 2)
168  {
169  if (strncmp(mesh_file, "../../data/inline-quad.mesh",100) == 0)
170  {
171  exact_volume = 0.5;
172  exact_area = 1.194452300992437;
173  }
174  }
175  (void)(&exact_area); // suppress a warning
176 
178  const IntegrationRule* ir=nullptr;
179 
180  // Integration with algoim
181  double vol=0.0;
182 
183 #ifdef MFEM_USE_ALGOIM
184  double area=0.0;
185  DenseMatrix bmat; //gradients of the shape functions in isoparametric space
186  DenseMatrix pmat; //gradients of the shape functions in physical space
187  Vector inormal; //normal to the level set in isoparametric space
188  Vector tnormal; //normal to the level set in physical space
189  Vector lsfun; //level set function restricted to an element
190  DofTransformation *doftrans;
191  Array<int> vdofs;
192  for (int i=0; i<fespace.GetNE(); i++)
193  {
194  const FiniteElement* el=fespace.GetFE(i);
195 
196  //get the element transformation
197  trans = fespace.GetElementTransformation(i);
198 
199  //extract the element vector from the level-set
200  doftrans = fespace.GetElementVDofs(i,vdofs);
201  x.GetSubVector(vdofs, lsfun);
202 
203  //construct Algoim integration object
204  AlgoimIntegrationRule air(aorder, *el, *trans, lsfun);
205 
206  //compute the volume contribution from the element
207  ir = air.GetVolumeIntegrationRule();
208  for (int j = 0; j < ir->GetNPoints(); j++)
209  {
210  const IntegrationPoint &ip = ir->IntPoint(j);
211  trans->SetIntPoint(&ip);
212  vol += ip.weight * trans->Weight();
213  }
214 
215  //compute the perimeter/area contribution from the element
216  bmat.SetSize(el->GetDof(),el->GetDim());
217  pmat.SetSize(el->GetDof(),el->GetDim());
218  inormal.SetSize(el->GetDim());
219  tnormal.SetSize(el->GetDim());
220 
221  ir = air.GetSurfaceIntegrationRule();
222  for (int j = 0; j < ir->GetNPoints(); j++)
223  {
224  const IntegrationPoint &ip = ir->IntPoint(j);
225  trans->SetIntPoint(&ip);
226 
227  el->CalcDShape(ip,bmat);
228  Mult(bmat, trans->AdjugateJacobian(), pmat);
229  //compute the normal to the LS in isoparametric space
230  bmat.MultTranspose(lsfun,inormal);
231  //compute the normal to the LS in physical space
232  pmat.MultTranspose(lsfun,tnormal);
233  area += ip.weight * tnormal.Norml2() / inormal.Norml2();
234  }
235  }
236 
237  if (exact_volume > 0)
238  {
239  std::cout<<"Algoim Volume="<<vol<<" Error="<<vol-exact_volume<<std::endl;
240  std::cout<<"Algoim Area="<<area<<" Error="<<area-exact_area<<std::endl;
241  }
242  else
243  {
244  std::cout<<"Algoim Volume="<<vol<<std::endl;
245  std::cout<<"Algoim Area="<<area<<std::endl;
246  }
247 #endif
248 
249  //Perform standard MFEM integration
250  vol=0.0;
251  for (int i=0; i<fespace.GetNE(); i++)
252  {
253  const FiniteElement* el=fespace.GetFE(i);
254  //get the element transformation
255  trans = fespace.GetElementTransformation(i);
256 
257  //compute the volume contribution from the element
258  ir=&IntRules.Get(el->GetGeomType(), iorder);
259  for (int j = 0; j < ir->GetNPoints(); j++)
260  {
261  const IntegrationPoint &ip = ir->IntPoint(j);
262  trans->SetIntPoint(&ip);
263  double vlsf=x.GetValue(*trans,ip);
264  if (vlsf>=0.0)
265  {
266  vol += ip.weight * trans->Weight();
267  }
268  }
269  }
270 
271  if (exact_volume > 0.0)
272  {
273  std::cout<<"MFEM Volume="<<vol<<" Error="<<vol-exact_volume<<std::endl;
274  }
275  else
276  {
277  std::cout<<"MFEM Volume="<<vol<<std::endl;
278  }
279 
280  ParaViewDataCollection dacol("ParaViewLSF", mesh);
281  dacol.SetLevelsOfDetail(order);
282  dacol.RegisterField("LSF",&x);
283  dacol.SetTime(1.0);
284  dacol.SetCycle(1);
285  dacol.Save();
286 
287  delete ls_coeff;
288  delete mesh;
289 
290  return 0;
291 }
Abstract class for all finite elements.
Definition: fe_base.hpp:233
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
const IntegrationRule * GetSurfaceIntegrationRule()
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
Helper class for ParaView visualization data.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
const IntegrationRule * GetVolumeIntegrationRule()
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:579
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: gridfunc.cpp:449
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2841
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
double sphere_ls(const Vector &x)
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
int main(int argc, char *argv[])
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:214
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
void SetTime(double t)
Set physical time (for time-dependent simulations)
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:736
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
Class for integration point with weight.
Definition: intrules.hpp:31
int dim
Definition: ex24.cpp:53
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2415
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:835
void SetLevelsOfDetail(int levels_of_detail_)
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
virtual void Save() override
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3696
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
double sinusoidal_ls(const Vector &x)
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:769