MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
ex38.cpp
Go to the documentation of this file.
1// MFEM Example 38
2//
3// Compile with: make ex38
4//
5// Sample runs:
6// (since all sample runs require LAPACK, the * symbol is used to exclude them
7// from the automatically generated internal MFEM tests).
8// * ex38
9// * ex38 -i volumetric1d
10// * ex38 -i surface2d
11// * ex38 -i surface2d -o 4 -r 5
12// * ex38 -i volumetric2d
13// * ex38 -i volumetric2d -o 4 -r 5
14// * ex38 -i surface3d
15// * ex38 -i surface3d -o 4 -r 5
16// * ex38 -i volumetric3d
17// * ex38 -i volumetric3d -o 4 -r 5
18//
19// Description: This example code demonstrates the use of MFEM to integrate
20// functions over implicit interfaces and subdomains bounded by
21// implicit interfaces.
22//
23// The quadrature rules are constructed by means of moment-fitting.
24// The interface is given by the zero isoline of a level-set
25// function ϕ and the subdomain is given as the domain where ϕ>0
26// holds. The algorithm for construction of the quadrature rules
27// was introduced by Mueller, Kummer and Oberlack [1].
28//
29// This example also showcases how to set up integrators using the
30// integration rules on implicit surfaces and subdomains.
31//
32// [1] Mueller, B., Kummer, F. and Oberlack, M. (2013) Highly accurate surface
33// and volume integration on implicit domains by means of moment-fitting.
34// Int. J. Numer. Meth. Engr. (96) 512-528. DOI:10.1002/nme.4569
35
36#include "mfem.hpp"
37#include <iostream>
38
39using namespace std;
40using namespace mfem;
41
42/// @brief Integration rule the example should demonstrate
47
48/// @brief Level-set function defining the implicit interface
50{
51 switch (itype)
52 {
54 return .55 - X(0);
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.));
63 default:
64 return 1.;
65 }
66}
67
68/// @brief Function that should be integrated
70{
71 switch (itype)
72 {
74 return 1.;
76 return 3. * pow(X(0), 2.) - pow(X(1), 2.);
78 return 1.;
80 return 4. - 3. * pow(X(0), 2.) + 2. * pow(X(1), 2.) - pow(X(2), 2.);
82 return 1.;
83 default:
84 return 0.;
85 }
86}
87
88/// @brief Analytic surface integral
90{
91 switch (itype)
92 {
94 return 1.;
96 return 2. * M_PI;
98 return 7.26633616541076;
100 return 40. / 3. * M_PI;
102 return 9.90182151329315;
103 default:
104 return 0.;
105 }
106}
107
108/// @brief Analytic volume integral over subdomain with positive level-set
110{
111 switch (itype)
112 {
114 return .55;
116 return NAN;
118 return 9. / 8. * M_PI;
120 return NAN;
122 return 3. / 4. * M_PI;
123 default:
124 return 0.;
125 }
126}
127
128#ifdef MFEM_USE_LAPACK
129/**
130 @brief Class for surface IntegrationRule
131
132 This class demonstrates how IntegrationRules computed as CutIntegrationRules
133 can be saved to reduce the impact by computing them from scratch each time.
134 */
135class SIntegrationRule : public IntegrationRule
136{
137protected:
138 /// @brief Space Dimension of the IntegrationRule
139 int dim;
140 /// @brief Column-wise matrix of the quadtrature weights
141 DenseMatrix Weights;
142 /// @brief Column-wise matrix of the transformation weights of the normal
143 DenseMatrix SurfaceWeights;
144
145public:
146 /**
147 @brief Constructor of SIntegrationRule
148
149 The surface integrationRules are computed and saved in the constructor.
150
151 @param [in] Order Order of the IntegrationRule
152 @param [in] LvlSet Level-set defining the implicit interface
153 @param [in] lsOrder Polynomial degree for approx of level-set function
154 @param [in] mesh Pointer to the mesh that is used
155 */
156 SIntegrationRule(int Order, Coefficient& LvlSet, int lsOrder, Mesh* mesh)
157 {
158 dim = mesh->Dimension();
159
161 MomentFittingIntRules MFIRs(Order, LvlSet, lsOrder);
162 mesh->GetElementTransformation(0, &Tr);
164 MFIRs.GetSurfaceIntegrationRule(Tr, ir);
165 if (dim >1)
166 {
167 Weights.SetSize(ir.GetNPoints(), mesh->GetNE());
168 }
169 else
170 {
171 Weights.SetSize(2, mesh->GetNE());
172 }
173 SurfaceWeights.SetSize(ir.GetNPoints(), mesh->GetNE());
174 Vector w;
175 MFIRs.GetSurfaceWeights(Tr, ir, w);
176 SurfaceWeights.SetCol(0, w);
177 SetSize(ir.GetNPoints());
178
179 for (int ip = 0; ip < GetNPoints(); ip++)
180 {
181 IntPoint(ip).index = ip;
182 IntegrationPoint &intp = IntPoint(ip);
183 intp.x = ir.IntPoint(ip).x;
184 intp.y = ir.IntPoint(ip).y;
185 intp.z = ir.IntPoint(ip).z;
186 if (dim > 1)
187 {
188 Weights(ip, 0) = ir.IntPoint(ip).weight;
189 }
190 else
191 {
192 Weights(0, 0) = ir.IntPoint(ip).x;
193 Weights(1, 0) = ir.IntPoint(ip).weight;
194 }
195 }
196
197
198 for (int elem = 1; elem < mesh->GetNE(); elem++)
199 {
200 mesh->GetElementTransformation(elem, &Tr);
201 MFIRs.GetSurfaceIntegrationRule(Tr, ir);
202 MFIRs.GetSurfaceWeights(Tr, ir, w);
203 SurfaceWeights.SetCol(elem, w);
204
205 for (int ip = 0; ip < GetNPoints(); ip++)
206 {
207 if (dim > 1)
208 {
209 Weights(ip, elem) = ir.IntPoint(ip).weight;
210 }
211 else
212 {
213 Weights(0, elem) = ir.IntPoint(ip).x;
214 Weights(1, elem) = ir.IntPoint(ip).weight;
215 }
216 }
217 }
218 }
219
220 /**
221 @brief Set the weights for the given element and multiply them with the
222 transformation of the interface
223 */
224 void SetElementinclSurfaceWeight(int Element)
225 {
226 if (dim == 1)
227 {
228 IntegrationPoint &intp = IntPoint(0);
229 intp.x = Weights(0, Element);
230 intp.weight = Weights(1, Element);
231 cout << intp.x << " " << Element << endl;
232 }
233 else
234 for (int ip = 0; ip < GetNPoints(); ip++)
235 {
236 IntegrationPoint &intp = IntPoint(ip);
237 intp.weight = Weights(ip, Element) * SurfaceWeights(ip, Element);
238 }
239 }
240
241 /// @brief Set the weights for the given element
242 void SetElement(int Element)
243 {
244 if (dim == 1)
245 {
246 IntegrationPoint &intp = IntPoint(0);
247 intp.x = Weights(0, Element);
248 intp.weight = Weights(1, Element);
249 }
250 else
251 for (int ip = 0; ip < GetNPoints(); ip++)
252 {
253 IntegrationPoint &intp = IntPoint(ip);
254 intp.weight = Weights(ip, Element);
255 }
256 }
257
258 /// @brief Destructor of SIntegrationRule
259 ~SIntegrationRule() {}
260};
261
262/**
263 @brief Class for volume IntegrationRule
264
265 This class demonstrates how IntegrationRules computed as CutIntegrationRules
266 can be saved to reduce the impact by computing them from scratch each time.
267 */
268class CIntegrationRule : public IntegrationRule
269{
270protected:
271 /// @brief Space Dimension of the IntegrationRule
272 int dim;
273 /// @brief Column-wise matrix of the quadtrature weights
274 DenseMatrix Weights;
275
276public:
277 /**
278 @brief Constructor of CIntegrationRule
279
280 The volume integrationRules are computed and saved in the constructor.
281
282 @param [in] Order Order of the IntegrationRule
283 @param [in] LvlSet Level-set defining the implicit interface
284 @param [in] lsOrder Polynomial degree for approx of level-set function
285 @param [in] mesh Pointer to the mesh that is used
286 */
287 CIntegrationRule(int Order, Coefficient& LvlSet, int lsOrder, Mesh* mesh)
288 {
289 dim = mesh->Dimension();
290
292 MomentFittingIntRules MFIRs(Order, LvlSet, lsOrder);
293 mesh->GetElementTransformation(0, &Tr);
295 MFIRs.GetVolumeIntegrationRule(Tr, ir);
296 if (dim > 1)
297 {
298 Weights.SetSize(ir.GetNPoints(), mesh->GetNE());
299 }
300 else
301 {
302 Weights.SetSize(2 * ir.GetNPoints(), mesh->GetNE());
303 }
304
305 SetSize(ir.GetNPoints());
306 for (int ip = 0; ip < GetNPoints(); ip++)
307 {
308 IntPoint(ip).index = ip;
309 IntegrationPoint &intp = IntPoint(ip);
310 intp.x = ir.IntPoint(ip).x;
311 intp.y = ir.IntPoint(ip).y;
312 intp.z = ir.IntPoint(ip).z;
313 if (dim > 1)
314 {
315 Weights(ip, 0) = ir.IntPoint(ip).weight;
316 }
317 else
318 {
319 Weights(2 * ip, 0) = ir.IntPoint(ip).x;
320 Weights(2 * ip + 1, 0) = ir.IntPoint(ip).weight;
321 }
322 }
323
324 for (int elem = 1; elem < mesh->GetNE(); elem++)
325 {
326 mesh->GetElementTransformation(elem, &Tr);
327 MFIRs.GetVolumeIntegrationRule(Tr, ir);
328
329 for (int ip = 0; ip < GetNPoints(); ip++)
330 {
331 if (dim > 1)
332 {
333 Weights(ip, elem) = ir.IntPoint(ip).weight;
334 }
335 else
336 {
337 Weights(2 * ip, elem) = ir.IntPoint(ip).x;
338 Weights(2 * ip + 1, elem) = ir.IntPoint(ip).weight;
339 }
340 }
341 }
342 }
343
344 /// @brief Set the weights for the given element
345 void SetElement(int Element)
346 {
347 if (dim == 1)
348 for (int ip = 0; ip < GetNPoints(); ip++)
349 {
350 IntegrationPoint &intp = IntPoint(ip);
351 intp.x = Weights(2 * ip, Element);
352 intp.weight = Weights(2 * ip + 1, Element);
353 }
354 else
355 for (int ip = 0; ip < GetNPoints(); ip++)
356 {
357 IntegrationPoint &intp = IntPoint(ip);
358 intp.weight = Weights(ip, Element);
359 }
360 }
361
362 /// @brief Destructor of CIntegrationRule
363 ~CIntegrationRule() {}
364};
365/**
366 @brief Class for surface linearform integrator
367
368 Integrator to demonstrate the use of the surface integration rule on an
369 implicit surface defined by a level-set.
370 */
371class SurfaceLFIntegrator : public LinearFormIntegrator
372{
373protected:
374 /// @brief vector to evaluate the basis functions
375 Vector shape;
376
377 /// @brief surface integration rule
378 SIntegrationRule* SIntRule;
379
380 /// @brief coefficient representing the level-set defining the interface
381 Coefficient &LevelSet;
382
383 /// @brief coefficient representing the integrand
384 Coefficient &Q;
385
386public:
387 /**
388 @brief Constructor for the surface linear form integrator
389
390 Constructor for the surface linear form integrator to demonstrate the use
391 of the surface integration rule by means of moment-fitting.
392
393 @param [in] q coefficient representing the inegrand
394 @param [in] levelset level-set defining the implicit interfac
395 @param [in] ir surface integrtion rule to be used
396 */
397 SurfaceLFIntegrator(Coefficient &q, Coefficient &levelset,
398 SIntegrationRule* ir)
399 : LinearFormIntegrator(), SIntRule(ir), LevelSet(levelset), Q(q) { }
400
401 /**
402 @brief Assembly of the element vector
403
404 Assemble the element vector of for the right hand side on the element given
405 by the FiniteElement and ElementTransformation.
406
407 @param [in] el finite Element the vector belongs to
408 @param [in] Tr transformation of finite element
409 @param [out] elvect vector containing the
410 */
411 virtual void AssembleRHSElementVect(const FiniteElement &el,
413 Vector &elvect) override
414 {
415 int dof = el.GetDof();
416 shape.SetSize(dof);
417 elvect.SetSize(dof);
418 elvect = 0.;
419
420 // Update the surface integration rule for the current element
421 SIntRule->SetElementinclSurfaceWeight(Tr.ElementNo);
422
423 for (int ip = 0; ip < SIntRule->GetNPoints(); ip++)
424 {
425 Tr.SetIntPoint((&(SIntRule->IntPoint(ip))));
426 real_t val = Tr.Weight() * Q.Eval(Tr, SIntRule->IntPoint(ip));
427 el.CalcShape(SIntRule->IntPoint(ip), shape);
428 add(elvect, SIntRule->IntPoint(ip).weight * val, shape, elvect);
429 }
430 }
431};
432
433/**
434 @brief Class for subdomain linearform integrator
435
436 Integrator to demonstrate the use of the subdomain integration rule within
437 an area defined by an implicit surface defined by a level-set.
438 */
439class SubdomainLFIntegrator : public LinearFormIntegrator
440{
441protected:
442 /// @brief vector to evaluate the basis functions
443 Vector shape;
444
445 /// @brief surface integration rule
446 CIntegrationRule* CIntRule;
447
448 /// @brief coefficient representing the level-set defining the interface
449 Coefficient &LevelSet;
450
451 /// @brief coefficient representing the integrand
452 Coefficient &Q;
453
454public:
455 /**
456 @brief Constructor for the volumetric subdomain linear form integrator
457
458 Constructor for the subdomain linear form integrator to demonstrate the use
459 of the volumetric subdomain integration rule by means of moment-fitting.
460
461 @param [in] q coefficient representing the inegrand
462 @param [in] levelset level-set defining the implicit interfac
463 @param [in] ir subdomain integrtion rule to be used
464 */
465 SubdomainLFIntegrator(Coefficient &q, Coefficient &levelset,
466 CIntegrationRule* ir)
467 : LinearFormIntegrator(), CIntRule(ir), LevelSet(levelset), Q(q) { }
468
469 /**
470 @brief Assembly of the element vector
471
472 Assemble the element vector of for the right hand side on the element given
473 by the FiniteElement and ElementTransformation.
474
475 @param [in] el finite Element the vector belongs to
476 @param [in] Tr transformation of finite element
477 @param [out] elvect vector containing the
478 */
479 virtual void AssembleRHSElementVect(const FiniteElement &el,
481 Vector &elvect) override
482 {
483 int dof = el.GetDof();
484 shape.SetSize(dof);
485 elvect.SetSize(dof);
486 elvect = 0.;
487
488 // Update the subdomain integration rule
489 CIntRule->SetElement(Tr.ElementNo);
490
491 for (int ip = 0; ip < CIntRule->GetNPoints(); ip++)
492 {
493 Tr.SetIntPoint((&(CIntRule->IntPoint(ip))));
494 real_t val = Tr.Weight()
495 * Q.Eval(Tr, CIntRule->IntPoint(ip));
496 el.CalcPhysShape(Tr, shape);
497 add(elvect, CIntRule->IntPoint(ip).weight * val, shape, elvect);
498 }
499 }
500};
501#endif // MFEM_USE_LAPACK
502
503int main(int argc, char *argv[])
504{
505#ifndef MFEM_USE_LAPACK
506 cout << "MFEM must be built with LAPACK for this example." << endl;
507 return MFEM_SKIP_RETURN_VALUE;
508#else
509 // 1. Parse he command-line options.
510 int ref_levels = 3;
511 int order = 2;
512 const char *inttype = "surface2d";
513 bool visualization = true;
515
516 OptionsParser args(argc, argv);
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.");
524 args.ParseCheck();
525
526 if (strcmp(inttype, "volumetric1d") == 0
527 || strcmp(inttype, "Volumetric1D") == 0)
528 {
530 }
531 else if (strcmp(inttype, "surface2d") == 0
532 || strcmp(inttype, "Surface2D") == 0)
533 {
535 }
536 else if (strcmp(inttype, "volumetric2d") == 0
537 || strcmp(inttype, "Volumetric2D") == 0)
538 {
540 }
541 else if (strcmp(inttype, "surface3d") == 0
542 || strcmp(inttype, "Surface3d") == 0)
543 {
545 }
546 else if (strcmp(inttype, "volumetric3d") == 0
547 || strcmp(inttype, "Volumetric3d") == 0)
548 {
550 }
551
552 // 2. Construct and refine the mesh.
553 Mesh *mesh;
555 {
556 mesh = new Mesh("../data/inline-segment.mesh");
557 }
560 {
561 mesh = new Mesh(2, 4, 1, 0, 2);
562 mesh->AddVertex(-1.6,-1.6);
563 mesh->AddVertex(1.6,-1.6);
564 mesh->AddVertex(1.6,1.6);
565 mesh->AddVertex(-1.6,1.6);
566 mesh->AddQuad(0,1,2,3);
567 mesh->FinalizeQuadMesh(1, 0, 1);
568 }
571 {
572 mesh = new Mesh(3, 8, 1, 0, 3);
573 mesh->AddVertex(-1.6,-1.6,-1.6);
574 mesh->AddVertex(1.6,-1.6,-1.6);
575 mesh->AddVertex(1.6,1.6,-1.6);
576 mesh->AddVertex(-1.6,1.6,-1.6);
577 mesh->AddVertex(-1.6,-1.6,1.6);
578 mesh->AddVertex(1.6,-1.6,1.6);
579 mesh->AddVertex(1.6,1.6,1.6);
580 mesh->AddVertex(-1.6,1.6,1.6);
581 mesh->AddHex(0,1,2,3,4,5,6,7);
582 mesh->FinalizeHexMesh(1, 0, 1);
583 }
584
585 for (int lev = 0; lev < ref_levels; lev++)
586 {
587 mesh->UniformRefinement();
588 }
589
590 // 3. Define the necessary finite element space on the mesh.
591 H1_FECollection fe_coll(1, mesh->Dimension());
592 FiniteElementSpace *fespace = new FiniteElementSpace(mesh, &fe_coll);
593
594 // 4. Construction Coefficients for the level set and the integrand.
595 FunctionCoefficient levelset(lvlset);
597
598 // 5. Define the necessary Integration rules on element 0.
600 mesh->GetElementTransformation(0, &Tr);
601 SIntegrationRule* sir = new SIntegrationRule(order, levelset, 2, mesh);
602 CIntegrationRule* cir = NULL;
606 {
607 cir = new CIntegrationRule(order, levelset, 2, mesh);
608 }
609
610 // 6. Define and assemble the linear forms on the finite element space.
611 LinearForm surface(fespace);
612 LinearForm volume(fespace);
613
614 surface.AddDomainIntegrator(new SurfaceLFIntegrator(u, levelset, sir));
615 surface.Assemble();
616
620 {
621 volume.AddDomainIntegrator(new SubdomainLFIntegrator(u, levelset, cir));
622 volume.Assemble();
623 }
624
625 // 7. Print information, computed values and errors to the console.
626 int qorder = 0;
627 int nbasis = 2 * (order + 1) + (int)(order * (order + 1) / 2);
629 IntegrationRule ir = irs.Get(Geometry::SQUARE, qorder);
630 for (; ir.GetNPoints() <= nbasis; qorder++)
631 {
632 ir = irs.Get(Geometry::SQUARE, qorder);
633 }
634 cout << "============================================" << endl;
635 cout << "Mesh size dx: ";
637 {
638 cout << 3.2 / pow(2., (real_t)ref_levels) << endl;
639 }
640 else
641 {
642 cout << .25 / pow(2., (real_t)ref_levels) << endl;
643 }
646 {
647 cout << "Number of div free basis functions: " << nbasis << endl;
648 cout << "Number of quadrature points: " << ir.GetNPoints() << endl;
649 }
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): ";
657 cout << abs(surface.Sum() - Surface()) / Surface() << endl;
661 {
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): ";
668 cout << abs(volume.Sum() - Volume()) / Volume() << endl;
669 }
670 cout << "============================================" << endl;
671
672 // 8. Plot the level-set function on a high order finite element space.
673 if (visualization)
674 {
675 H1_FECollection fe_coll2(5, mesh->Dimension());
676 FiniteElementSpace fespace2(mesh, &fe_coll2);
677 FunctionCoefficient levelset_coeff(levelset);
678 GridFunction lgf(&fespace2);
679 lgf.ProjectCoefficient(levelset_coeff);
680 char vishost[] = "localhost";
681 int visport = 19916;
682 socketstream sol_sock(vishost, visport);
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;
687 }
688
689 delete sir;
690 delete cir;
691 delete fespace;
692 delete mesh;
693 return EXIT_SUCCESS;
694#endif //MFEM_USE_LAPACK
695}
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.
Definition densemat.hpp:24
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
void SetCol(int c, const real_t *col)
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:131
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:93
Abstract data type element.
Definition element.hpp:29
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
Abstract class for all finite elements.
Definition fe_base.hpp:239
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.
Definition fe_base.hpp:329
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 ...
Definition fe_base.cpp:182
A general function coefficient.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
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.
Definition fe_coll.hpp:260
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
Container class for integration rules.
Definition intrules.hpp:416
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
A standard isoparametric element transformation.
Definition eltrans.hpp:363
Abstract base class LinearFormIntegrator.
Definition lininteg.hpp:25
Vector with associated FE space and LinearFormIntegrators.
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Mesh data type.
Definition mesh.hpp:56
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.
Definition mesh.cpp:1743
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition mesh.cpp:1658
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void FinalizeHexMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a hexahedral Mesh.
Definition mesh.cpp:3098
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...
Definition mesh.cpp:357
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition mesh.cpp:2177
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.
Definition mesh.cpp:1806
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
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...
Definition optparser.hpp:82
Vector data type.
Definition vector.hpp:80
real_t Sum() const
Return the sum of the vector entries.
Definition vector.cpp:1283
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
int dim
Definition ex24.cpp:53
real_t Surface()
Analytic surface integral.
Definition ex38.cpp:89
IntegrationType itype
Definition ex38.cpp:46
real_t integrand(const Vector &X)
Function that should be integrated.
Definition ex38.cpp:69
real_t Volume()
Analytic volume integral over subdomain with positive level-set.
Definition ex38.cpp:109
real_t lvlset(const Vector &X)
Level-set function defining the implicit interface.
Definition ex38.cpp:49
IntegrationType
Integration rule the example should demonstrate.
Definition ex38.cpp:43
int main()
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:316
float real_t
Definition config.hpp:43
const char vishost[]