56 return 1. - (pow(X(0), 2.) + pow(X(1), 2.));
58 return 1. - (pow(X(0) / 1.5, 2.) + pow(X(1) / .75, 2.));
60 return 1. - (pow(X(0), 2.) + pow(X(1), 2.) + pow(X(2), 2.));
62 return 1. - (pow(X(0) / 1.5, 2.) + pow(X(1) / .75, 2.) + pow(X(2) / .5, 2.));
76 return 3. * pow(X(0), 2.) - pow(X(1), 2.);
80 return 4. - 3. * pow(X(0), 2.) + 2. * pow(X(1), 2.) - pow(X(2), 2.);
98 return 7.26633616541076;
100 return 40. / 3. * M_PI;
102 return 9.90182151329315;
118 return 9. / 8. * M_PI;
122 return 3. / 4. * M_PI;
128#ifdef MFEM_USE_LAPACK
156 SIntegrationRule(
int Order,
Coefficient& LvlSet,
int lsOrder,
Mesh* mesh)
164 MFIRs.GetSurfaceIntegrationRule(Tr, ir);
175 MFIRs.GetSurfaceWeights(Tr, ir, w);
176 SurfaceWeights.
SetCol(0, w);
198 for (
int elem = 1; elem < mesh->
GetNE(); elem++)
201 MFIRs.GetSurfaceIntegrationRule(Tr, ir);
202 MFIRs.GetSurfaceWeights(Tr, ir, w);
203 SurfaceWeights.
SetCol(elem, w);
224 void SetElementinclSurfaceWeight(
int Element)
231 cout << intp.
x <<
" " <<
Element << endl;
259 ~SIntegrationRule() {}
287 CIntegrationRule(
int Order,
Coefficient& LvlSet,
int lsOrder,
Mesh* mesh)
295 MFIRs.GetVolumeIntegrationRule(Tr, ir);
324 for (
int elem = 1; elem < mesh->
GetNE(); elem++)
327 MFIRs.GetVolumeIntegrationRule(Tr, ir);
337 Weights(2 * ip, elem) = ir.
IntPoint(ip).
x;
363 ~CIntegrationRule() {}
378 SIntegrationRule* SIntRule;
398 SIntegrationRule* ir)
421 SIntRule->SetElementinclSurfaceWeight(Tr.
ElementNo);
423 for (
int ip = 0; ip < SIntRule->GetNPoints(); ip++)
427 el.
CalcShape(SIntRule->IntPoint(ip), shape);
428 add(elvect, SIntRule->IntPoint(ip).weight * val, shape, elvect);
446 CIntegrationRule* CIntRule;
466 CIntegrationRule* ir)
491 for (
int ip = 0; ip < CIntRule->GetNPoints(); ip++)
495 * Q.
Eval(Tr, CIntRule->IntPoint(ip));
497 add(elvect, CIntRule->IntPoint(ip).weight * val, shape, elvect);
503int main(
int argc,
char *argv[])
505#ifndef MFEM_USE_LAPACK
506 cout <<
"MFEM must be built with LAPACK for this example." << endl;
507 return MFEM_SKIP_RETURN_VALUE;
512 const char *inttype =
"surface2d";
513 bool visualization =
true;
517 args.
AddOption(&order,
"-o",
"--order",
"Order of quadrature rule");
518 args.
AddOption(&ref_levels,
"-r",
"--refine",
"Number of meh refinements");
519 args.
AddOption(&inttype,
"-i",
"--integrationtype",
520 "IntegrationType to demonstrate");
521 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
522 "--no-visualization",
523 "Enable or disable GLVis visualization.");
526 if (strcmp(inttype,
"volumetric1d") == 0
527 || strcmp(inttype,
"Volumetric1D") == 0)
531 else if (strcmp(inttype,
"surface2d") == 0
532 || strcmp(inttype,
"Surface2D") == 0)
536 else if (strcmp(inttype,
"volumetric2d") == 0
537 || strcmp(inttype,
"Volumetric2D") == 0)
541 else if (strcmp(inttype,
"surface3d") == 0
542 || strcmp(inttype,
"Surface3d") == 0)
546 else if (strcmp(inttype,
"volumetric3d") == 0
547 || strcmp(inttype,
"Volumetric3d") == 0)
556 mesh =
new Mesh(
"../data/inline-segment.mesh");
561 mesh =
new Mesh(2, 4, 1, 0, 2);
572 mesh =
new Mesh(3, 8, 1, 0, 3);
581 mesh->
AddHex(0,1,2,3,4,5,6,7);
585 for (
int lev = 0; lev < ref_levels; lev++)
601 SIntegrationRule* sir =
new SIntegrationRule(order, levelset, 2, mesh);
602 CIntegrationRule* cir = NULL;
607 cir =
new CIntegrationRule(order, levelset, 2, mesh);
627 int nbasis = 2 * (order + 1) + (
int)(order * (order + 1) / 2);
634 cout <<
"============================================" << endl;
635 cout <<
"Mesh size dx: ";
638 cout << 3.2 / pow(2., (
real_t)ref_levels) << endl;
642 cout << .25 / pow(2., (
real_t)ref_levels) << endl;
647 cout <<
"Number of div free basis functions: " << nbasis << endl;
648 cout <<
"Number of quadrature points: " << ir.
GetNPoints() << endl;
650 cout << scientific << setprecision(2);
651 cout <<
"============================================" << endl;
652 cout <<
"Computed value of surface integral: " << surface.
Sum() << endl;
653 cout <<
"True value of surface integral: " <<
Surface() << endl;
654 cout <<
"Absolute Error (Surface): ";
655 cout << abs(surface.
Sum() -
Surface()) << endl;
656 cout <<
"Relative Error (Surface): ";
662 cout <<
"--------------------------------------------" << endl;
663 cout <<
"Computed value of volume integral: " << volume.
Sum() << endl;
664 cout <<
"True value of volume integral: " <<
Volume() << endl;
665 cout <<
"Absolute Error (Volume): ";
666 cout << abs(volume.
Sum() -
Volume()) << endl;
667 cout <<
"Relative Error (Volume): ";
670 cout <<
"============================================" << endl;
683 sol_sock.precision(8);
684 sol_sock <<
"solution\n" << *mesh << lgf << flush;
685 sol_sock <<
"keys pppppppppppppppppppppppppppcmmlRj\n";
686 sol_sock <<
"levellines " << 0. <<
" " << 0. <<
" " << 1 <<
"\n" << flush;
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
Data type dense matrix using column-major storage.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void SetCol(int c, const real_t *col)
Abstract data type element.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Abstract class for all finite elements.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
A general function coefficient.
Class for grid function - Vector with associated FE space.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Arbitrary order H1-conforming (continuous) finite elements.
Class for integration point with weight.
Class for an integration rule - an Array of IntegrationPoint.
int GetNPoints() const
Returns the number of the points in the integration rule.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Container class for integration rules.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
int GetNE() const
Returns number of elements.
int Dimension() const
Dimension of the reference space used within the elements.
void FinalizeHexMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a hexahedral Mesh.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
int AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8, int attr=1)
Adds a hexahedron to the mesh given by 8 vertices v1 through v8.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Class for subdomain IntegrationRules by means of moment-fitting.
void ParseCheck(std::ostream &out=mfem::out)
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 'var' to receive the value. Enable/disable tags are used to set the bool...
real_t Sum() const
Return the sum of the vector entries.
void SetSize(int s)
Resize the vector to size s.
real_t Surface()
Analytic surface integral.
real_t integrand(const Vector &X)
Function that should be integrated.
real_t Volume()
Analytic volume integral over subdomain with positive level-set.
real_t lvlset(const Vector &X)
Level-set function defining the implicit interface.
IntegrationType
Integration rule the example should demonstrate.
real_t u(const Vector &xvec)
void add(const Vector &v1, const Vector &v2, Vector &v)