MFEM  v4.6.0
Finite element discretization library
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 
41 using namespace std;
42 using namespace mfem;
43 
44 // Piecewise-affine function which is sometimes mesh-conforming
45 double affine_function(const Vector &p)
46 {
47  double 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
59 double jump_function(const Vector &p)
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].
73 double singular_function(const Vector &p)
74 {
75  double x = p(0), y = p(1);
76  double alpha = 1000.0;
77  double xc = 0.75, yc = 0.5;
78  double r0 = 0.7;
79  double r = sqrt(pow(x - xc,2.0) + pow(y - yc,2.0));
80  double num = - ( alpha - pow(alpha,3) * (pow(r,2) - pow(r0,2)) );
81  double 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 = max(denom,1e-8);
84  return num / denom;
85 }
86 
87 int 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();
93  Hypre::Init();
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  double double_max_elems = double(max_elems);
101  bool visualization = true;
102  bool nc_simplices = true;
103  double 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.
167  FunctionCoefficient affine_coeff(affine_function);
169  FunctionCoefficient singular_coeff(singular_function);
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).
182  const IntegrationRule *irs[Geometry::NumGeom];
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  double 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 }
int visport
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
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 main(int argc, char *argv[])
Definition: ex30p.cpp:87
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void ResetCoefficient(Coefficient &coeff_)
Reset the function f.
double jump_function(const Vector &p)
Definition: ex30p.cpp:59
double affine_function(const Vector &p)
Definition: ex30p.cpp:45
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.
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
void SetMaxElements(long long max_elements_)
Set the maximum number of elements stopping criterion: stop when the input mesh has num_elements >= m...
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:1115
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9888
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
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:5635
void SetNCLimit(int nc_limit_)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited). The default value is...
double singular_function(const Vector &p)
Definition: ex30p.cpp:73
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
Refinement operator to control data oscillation.
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
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
void SetThreshold(double threshold_)
Set the refinement threshold. The default value is 1.0e-2.
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:678
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
Class for parallel meshes.
Definition: pmesh.hpp:32