MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
ex30.cpp
Go to the documentation of this file.
1// MFEM Example 30
2//
3// Compile with: make ex30
4//
5// Sample runs: ex30 -m ../data/square-disc.mesh -o 1
6// ex30 -m ../data/square-disc.mesh -o 2
7// ex30 -m ../data/square-disc.mesh -o 2 -me 1e+4
8// ex30 -m ../data/square-disc-nurbs.mesh -o 2
9// ex30 -m ../data/star.mesh -o 2 -eo 4
10// ex30 -m ../data/fichera.mesh -o 2 -me 1e+5 -e 5e-2
11// ex30 -m ../data/disc-nurbs.mesh -o 2
12// ex30 -m ../data/ball-nurbs.mesh -o 2 -eo 3 -e 5e-2 -me 1e+5
13// ex30 -m ../data/star-surf.mesh -o 2
14// ex30 -m ../data/square-disc-surf.mesh -o 2
15// ex30 -m ../data/amr-quad.mesh -l 2
16//
17// Description: This is an example of adaptive mesh refinement preprocessing
18// which lowers the data oscillation [1] to a user-defined
19// relative threshold. There is no PDE being solved.
20//
21// MFEM's capability to work with both conforming and
22// nonconforming meshes is demonstrated in example 6. In some
23// problems, the material data or loading data is not sufficiently
24// resolved on the initial mesh. This missing fine scale data
25// reduces the accuracy of the solution as well as the accuracy of
26// some local error estimators. By preprocessing the mesh before
27// solving the PDE, many issues can be avoided.
28//
29// [1] Morin, P., Nochetto, R. H., & Siebert, K. G. (2000). Data
30// oscillation and convergence of adaptive FEM. SIAM Journal
31// on Numerical Analysis, 38(2), 466-488.
32//
33// [2] Mitchell, W. F. (2013). A collection of 2D elliptic
34// problems for testing adaptive grid refinement algorithms.
35// Applied mathematics and computation, 220, 350-364.
36
37#include "mfem.hpp"
38#include <fstream>
39#include <iostream>
40
41using namespace std;
42using namespace mfem;
43
44// Piecewise-affine function which is sometimes mesh-conforming
46{
47 real_t x = p(0), y = p(1);
48 if (x < 0.0)
49 {
50 return 1.0 + x + y;
51 }
52 else
53 {
54 return 1.0;
55 }
56}
57
58// Piecewise-constant function which is never mesh-conforming
60{
61 if (p.Normlp(2.0) > 0.4 && p.Normlp(2.0) < 0.6)
62 {
63 return 1.0;
64 }
65 else
66 {
67 return 5.0;
68 }
69}
70
71// Singular function derived from the Laplacian of the "steep wavefront" problem
72// in [2].
74{
75 real_t x = p(0), y = p(1);
76 real_t alpha = 1000.0;
77 real_t xc = 0.75, yc = 0.5;
78 real_t r0 = 0.7;
79 real_t r = sqrt(pow(x - xc,2.0) + pow(y - yc,2.0));
80 real_t num = - ( alpha - pow(alpha,3) * (pow(r,2) - pow(r0,2)) );
81 real_t denom = pow(r * ( pow(alpha,2) * pow(r0,2) + pow(alpha,2) * pow(r,2) \
82 - 2 * pow(alpha,2) * r0 * r + 1.0 ),2);
83 denom = std::max(denom, (real_t) 1.0e-8);
84 return num / denom;
85}
86
87int main(int argc, char *argv[])
88{
89 // 1. Parse command-line options.
90 const char *mesh_file = "../data/star.mesh";
91 int order = 1;
92 int nc_limit = 1;
93 int max_elems = 100*1000;
94 real_t double_max_elems = real_t(max_elems);
95 bool visualization = true;
96 real_t osc_threshold = 1e-3;
97 int enriched_order = 5;
98
99 OptionsParser args(argc, argv);
100 args.AddOption(&mesh_file, "-m", "--mesh",
101 "Mesh file to use.");
102 args.AddOption(&order, "-o", "--order",
103 "Finite element order (polynomial degree).");
104 args.AddOption(&nc_limit, "-l", "--nc-limit",
105 "Maximum level of hanging nodes.");
106 args.AddOption(&double_max_elems, "-me", "--max-elems",
107 "Stop after reaching this many elements.");
108 args.AddOption(&osc_threshold, "-e", "--error",
109 "relative data oscillation threshold.");
110 args.AddOption(&enriched_order, "-eo", "--enriched_order",
111 "Enriched quadrature order.");
112 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
113 "--no-visualization",
114 "Enable or disable GLVis visualization.");
115 args.Parse();
116 if (!args.Good())
117 {
118 args.PrintUsage(cout);
119 return 1;
120 }
121 args.PrintOptions(cout);
122
123 max_elems = int(double_max_elems);
124 Mesh mesh(mesh_file, 1, 1);
125
126 // 2. Since a NURBS mesh can currently only be refined uniformly, we need to
127 // convert it to a piecewise-polynomial curved mesh. First we refine the
128 // NURBS mesh a bit and then project the curvature to quadratic Nodes.
129 if (mesh.NURBSext)
130 {
131 for (int i = 0; i < 2; i++)
132 {
133 mesh.UniformRefinement();
134 }
135 mesh.SetCurvature(2);
136 }
137
138 // 3. Define functions and refiner.
142 CoefficientRefiner coeffrefiner(affine_coeff, order);
143
144 // 4. Connect to GLVis.
145 char vishost[] = "localhost";
146 int visport = 19916;
147 socketstream sol_sock;
148 if (visualization)
149 {
150 sol_sock.open(vishost, visport);
151 }
152
153 // 5. Define custom integration rule (optional).
155 int order_quad = 2*order + enriched_order;
156 for (int i = 0; i < Geometry::NumGeom; ++i)
157 {
158 irs[i] = &(IntRules.Get(i, order_quad));
159 }
160
161 // 6. Apply custom refiner settings.
162 coeffrefiner.SetIntRule(irs);
163 coeffrefiner.SetMaxElements(max_elems);
164 coeffrefiner.SetThreshold(osc_threshold);
165 coeffrefiner.SetNCLimit(nc_limit);
166 coeffrefiner.PrintWarnings();
167
168 // 7. Preprocess mesh to control osc (piecewise-affine function). This is
169 // mostly just a verification check. The oscillation should be zero if the
170 // function is mesh-conforming and order > 0.
171 coeffrefiner.PreprocessMesh(mesh);
172
173 mfem::out << "\n";
174 mfem::out << "Function 0 (affine) \n";
175 mfem::out << "Number of Elements " << mesh.GetNE() << "\n";
176 mfem::out << "Osc error " << coeffrefiner.GetOsc() << "\n";
177
178 // 8. Preprocess mesh to control osc (jump function).
179 coeffrefiner.ResetCoefficient(jump_coeff);
180 coeffrefiner.PreprocessMesh(mesh);
181
182 mfem::out << "\n";
183 mfem::out << "Function 1 (discontinuous) \n";
184 mfem::out << "Number of Elements " << mesh.GetNE() << "\n";
185 mfem::out << "Osc error " << coeffrefiner.GetOsc() << "\n";
186
187 // 9. Preprocess mesh to control osc (singular function).
188 coeffrefiner.ResetCoefficient(singular_coeff);
189 coeffrefiner.PreprocessMesh(mesh);
190
191 mfem::out << "\n";
192 mfem::out << "Function 2 (singular) \n";
193 mfem::out << "Number of Elements " << mesh.GetNE() << "\n";
194 mfem::out << "Osc error " << coeffrefiner.GetOsc() << "\n";
195
196 if (visualization)
197 {
198 sol_sock.precision(8);
199 sol_sock << "mesh\n" << mesh << flush;
200 }
201
202 return 0;
203}
Refinement operator to control data oscillation.
void SetNCLimit(int nc_limit_)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited). The default value is...
void SetMaxElements(long long max_elements_)
Set the maximum number of elements stopping criterion: stop when the input mesh has num_elements >= m...
void SetIntRule(const IntegrationRule *irs_[])
virtual int PreprocessMesh(Mesh &mesh, int max_it)
Apply the operator to the mesh max_it times or until tolerance achieved.
void SetThreshold(real_t threshold_)
Set the refinement threshold. The default value is 1.0e-2.
void ResetCoefficient(Coefficient &coeff_)
Reset the function f.
A general function coefficient.
static const int NumGeom
Definition geom.hpp:42
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Mesh data type.
Definition mesh.hpp:56
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition mesh.cpp:6211
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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
bool Good() const
Return true if the command line options were parsed successfully.
Vector data type.
Definition vector.hpp:80
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
const real_t alpha
Definition ex15.cpp:369
real_t affine_function(const Vector &p)
Definition ex30.cpp:45
real_t singular_function(const Vector &p)
Definition ex30.cpp:73
real_t jump_function(const Vector &p)
Definition ex30.cpp:59
int main()
const int visport
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
float real_t
Definition config.hpp:43
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:486
const char vishost[]
real_t p(const Vector &x, real_t t)