MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
ex30p.cpp
Go to the documentation of this file.
1// MFEM Example 30 - Parallel Version
2//
3// Compile with: make ex30p
4//
5// Sample runs: mpirun -np 4 ex30p -m ../data/square-disc.mesh -o 1
6// mpirun -np 4 ex30p -m ../data/square-disc.mesh -o 2
7// mpirun -np 4 ex30p -m ../data/square-disc.mesh -o 2 -me 1e+4
8// mpirun -np 4 ex30p -m ../data/square-disc-nurbs.mesh -o 2
9// mpirun -np 4 ex30p -m ../data/star.mesh -o 2 -eo 4
10// mpirun -np 4 ex30p -m ../data/fichera.mesh -o 2 -me 1e+5 -e 5e-2
11// mpirun -np 4 ex30p -m ../data/disc-nurbs.mesh -o 2
12// mpirun -np 4 ex30p -m ../data/ball-nurbs.mesh -o 2 -eo 3 -e 5e-2 -me 1e+5
13// mpirun -np 4 ex30p -m ../data/star-surf.mesh -o 2
14// mpirun -np 4 ex30p -m ../data/square-disc-surf.mesh -o 2
15// mpirun -np 4 ex30p -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 // 0. Initialize MPI and HYPRE.
90 Mpi::Init(argc, argv);
91 int num_procs = Mpi::WorldSize();
92 int myid = Mpi::WorldRank();
94
95 // 1. Parse command-line options.
96 const char *mesh_file = "../data/star.mesh";
97 int order = 1;
98 int nc_limit = 1;
99 int max_elems = 1e5;
100 real_t double_max_elems = real_t(max_elems);
101 bool visualization = true;
102 bool nc_simplices = true;
103 real_t osc_threshold = 1e-3;
104 int enriched_order = 5;
105
106 OptionsParser args(argc, argv);
107 args.AddOption(&mesh_file, "-m", "--mesh",
108 "Mesh file to use.");
109 args.AddOption(&order, "-o", "--order",
110 "Finite element order (polynomial degree).");
111 args.AddOption(&nc_limit, "-l", "--nc-limit",
112 "Maximum level of hanging nodes.");
113 args.AddOption(&double_max_elems, "-me", "--max-elems",
114 "Stop after reaching this many elements.");
115 args.AddOption(&osc_threshold, "-e", "--error",
116 "relative data oscillation threshold.");
117 args.AddOption(&enriched_order, "-eo", "--enriched_order",
118 "Enriched quadrature order.");
119 args.AddOption(&nc_simplices, "-ns", "--nonconforming-simplices",
120 "-cs", "--conforming-simplices",
121 "For simplicial meshes, enable/disable nonconforming"
122 " refinement");
123 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
124 "--no-visualization",
125 "Enable or disable GLVis visualization.");
126 args.Parse();
127 if (!args.Good())
128 {
129 if (myid == 0)
130 {
131 args.PrintUsage(cout);
132 }
133 return 1;
134 }
135 if (myid == 0)
136 {
137 args.PrintOptions(cout);
138 }
139
140 max_elems = int(double_max_elems);
141 Mesh mesh(mesh_file, 1, 1);
142
143 // 2. Since a NURBS mesh can currently only be refined uniformly, we need to
144 // convert it to a piecewise-polynomial curved mesh. First we refine the
145 // NURBS mesh a bit and then project the curvature to quadratic Nodes.
146 if (mesh.NURBSext)
147 {
148 for (int i = 0; i < 2; i++)
149 {
150 mesh.UniformRefinement();
151 }
152 mesh.SetCurvature(2);
153 }
154
155 // 3. Make sure the mesh is in the non-conforming mode to enable local
156 // refinement of quadrilaterals/hexahedra. Simplices can be refined either
157 // in conforming or in non-conforming mode. The conforming mode however
158 // does not support dynamic partitioning.
159 mesh.EnsureNCMesh(nc_simplices);
160
161 // 4. Define a parallel mesh by partitioning the serial mesh. Once the
162 // parallel mesh is defined, the serial mesh can be deleted.
163 ParMesh pmesh(MPI_COMM_WORLD, mesh);
164 mesh.Clear();
165
166 // 5. Define functions and refiner.
170 CoefficientRefiner coeffrefiner(affine_coeff,order);
171
172 // 6. Connect to GLVis.
173 char vishost[] = "localhost";
174 int visport = 19916;
175 socketstream sol_sock;
176 if (visualization)
177 {
178 sol_sock.open(vishost, visport);
179 }
180
181 // 7. Define custom integration rule (optional).
183 int order_quad = 2*order + enriched_order;
184 for (int i = 0; i < Geometry::NumGeom; ++i)
185 {
186 irs[i] = &(IntRules.Get(i, order_quad));
187 }
188
189 // 8. Apply custom refiner settings.
190 coeffrefiner.SetIntRule(irs);
191 coeffrefiner.SetMaxElements(max_elems);
192 coeffrefiner.SetThreshold(osc_threshold);
193 coeffrefiner.SetNCLimit(nc_limit);
194 coeffrefiner.PrintWarnings();
195
196 // 9. Preprocess mesh to control osc (piecewise-affine function). This is
197 // mostly just a verification check. The oscillation should be zero if the
198 // function is mesh-conforming and order > 0.
199 coeffrefiner.PreprocessMesh(pmesh);
200
201 int globalNE = pmesh.GetGlobalNE();
202 real_t osc = coeffrefiner.GetOsc();
203 if (myid == 0)
204 {
205 mfem::out << "\n";
206 mfem::out << "Function 0 (affine) \n";
207 mfem::out << "Number of Elements " << globalNE << "\n";
208 mfem::out << "Osc error " << osc << "\n";
209 }
210
211 // 10. Preprocess mesh to control osc (jump function).
212 coeffrefiner.ResetCoefficient(jump_coeff);
213 coeffrefiner.PreprocessMesh(pmesh);
214
215 globalNE = pmesh.GetGlobalNE();
216 osc = coeffrefiner.GetOsc();
217 if (myid == 0)
218 {
219 mfem::out << "\n";
220 mfem::out << "Function 1 (discontinuous) \n";
221 mfem::out << "Number of Elements " << globalNE << "\n";
222 mfem::out << "Osc error " << osc << "\n";
223 }
224
225 // 11. Preprocess mesh to control osc (singular function).
226 coeffrefiner.ResetCoefficient(singular_coeff);
227 coeffrefiner.PreprocessMesh(pmesh);
228
229 globalNE = pmesh.GetGlobalNE();
230 osc = coeffrefiner.GetOsc();
231 if (myid == 0)
232 {
233 mfem::out << "\n";
234 mfem::out << "Function 2 (singular) \n";
235 mfem::out << "Number of Elements " << globalNE << "\n";
236 mfem::out << "Osc error " << osc << "\n";
237 }
238
239 if (visualization)
240 {
241 sol_sock.precision(8);
242 sol_sock << "parallel " << num_procs << " " << myid << "\n";
243 sol_sock << "mesh\n" << pmesh << flush;
244 }
245
246 return 0;
247}
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
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
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
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition mesh.hpp:1255
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition mesh.cpp:10626
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
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
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.
Class for parallel meshes.
Definition pmesh.hpp:34
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 ex30p.cpp:45
real_t singular_function(const Vector &p)
Definition ex30p.cpp:73
real_t jump_function(const Vector &p)
Definition ex30p.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)