MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
hpref.cpp
Go to the documentation of this file.
1// Serial hp-refinement example
2//
3// Compile with: make hpref
4//
5// Sample runs: hpref -dim 2 -n 1000
6// hpref -dim 3 -n 500
7//
8// Description: This example demonstrates h- and p-refinement in a serial
9// finite element discretization of the Poisson problem (cf. ex1)
10// -Delta u = 1 with homogeneous Dirichlet boundary conditions.
11// Refinements are performed iteratively, each iteration having h-
12// or p-refinements. For simplicity, we randomly choose the
13// elements and the type of refinement, for each iteration. In
14// practice, these choices may be made in a problem-dependent way,
15// but this example serves only to illustrate the capabilities of
16// hp-refinement.
17//
18// We recommend viewing Example 1 before viewing this example.
19
20#include "mfem.hpp"
21#include <fstream>
22#include <iostream>
23
24using namespace std;
25using namespace mfem;
26
28
29// Deterministic function for "random" integers.
30int DetRand(int & seed)
31{
32 seed++;
33 return int(std::abs(1.0e5 * sin(seed * 1.1234 * M_PI)));
34}
35
36void f_exact(const Vector &x, Vector &f);
37
38int main(int argc, char *argv[])
39{
40 // 1. Parse command-line options.
41 int order = 1;
42 const char *device_config = "cpu";
43 bool visualization = true;
44 int numIter = 0;
45 int dim = 2;
46 bool deterministic = true;
47 bool projectSolution = false;
48
49 OptionsParser args(argc, argv);
50 args.AddOption(&order, "-o", "--order",
51 "Finite element order (polynomial degree) or -1 for"
52 " isoparametric space.");
53 args.AddOption(&device_config, "-d", "--device",
54 "Device configuration string, see Device::Configure().");
55 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
56 "--no-visualization",
57 "Enable or disable GLVis visualization.");
58 args.AddOption(&numIter, "-n", "--num-iter", "Number of hp-ref iterations");
59 args.AddOption(&dim, "-dim", "--dim", "Mesh dimension (2 or 3)");
60 args.AddOption(&deterministic, "-det", "--deterministic", "-not-det",
61 "--not-deterministic",
62 "Whether to use deterministic random refinements");
63 args.AddOption(&projectSolution, "-proj", "--project-solution", "-no-proj",
64 "--no-project",
65 "Whether to project a coefficient to solution");
66 args.Parse();
67 if (!args.Good())
68 {
69 args.PrintUsage(cout);
70 return 1;
71 }
72 args.PrintOptions(cout);
73
74 // 2. Enable hardware devices such as GPUs, and programming models such as
75 // CUDA, OCCA, RAJA and OpenMP based on command line options.
76 Device device(device_config);
77 device.Print();
78
79 // 3. Construct a uniform coarse mesh on all processors.
80 Mesh mesh;
81 if (dim == 3)
82 {
84 }
85 else
86 {
88 }
89
90 mesh.EnsureNCMesh();
91
92 // 4. Define a finite element space on the mesh. Here we use continuous
93 // Lagrange finite elements of the specified order. If order < 1, we
94 // instead use an isoparametric/isogeometric space.
96 bool delete_fec;
97 if (order > 0)
98 {
99 fec = new H1_FECollection(order, dim);
100 delete_fec = true;
101 }
102 else if (mesh.GetNodes())
103 {
104 fec = mesh.GetNodes()->OwnFEC();
105 delete_fec = false;
106 cout << "Using isoparametric FEs: " << fec->Name() << endl;
107 }
108 else
109 {
110 fec = new H1_FECollection(order = 1, dim);
111 delete_fec = true;
112 }
113
114 const int fespaceDim = projectSolution ? dim : 1;
115 FiniteElementSpace fespace(&mesh, fec, fespaceDim);
116
117 // 5. Iteratively perform h- and p-refinements.
118 int numH = 0;
119 int numP = 0;
120 int seed = 0;
121
122 const std::vector<char> hp_char = {'h', 'p'};
123
124 for (int iter=0; iter<numIter; ++iter)
125 {
126 const int r1 = deterministic ? DetRand(seed) : rand();
127 const int r2 = deterministic ? DetRand(seed) : rand();
128 const int elem = r1 % mesh.GetNE();
129 const int hp = r2 % 2;
130
131 cout << "hp-refinement iteration " << iter << ": "
132 << hp_char[hp] << "-refinement" << endl;
133
134 if (hp == 1)
135 {
136 // p-ref
138 refs.Append(pRefinement(elem, 1)); // Increase the element order by 1
139 fespace.PRefineAndUpdate(refs);
140 numP++;
141 }
142 else
143 {
144 // h-ref
146 refs.Append(Refinement(elem));
147 mesh.GeneralRefinement(refs);
148 fespace.Update(false);
149 numH++;
150 }
151 }
152
153 const int size = fespace.GetTrueVSize();
154 cout << "Number of finite element unknowns: " << size << endl;
155
156 const int maxP = fespace.GetMaxElementOrder();
157 cout << "Total number of h-refinements: " << numH
158 << "\nTotal number of p-refinements: " << numP
159 << "\nMaximum order " << maxP << "\n";
160
161 GridFunction x(&fespace);
162 Vector X;
163
164 if (projectSolution)
165 {
167 x.ProjectCoefficient(vec_coef);
168
169 X.SetSize(fespace.GetTrueVSize());
170
171 fespace.GetHpRestrictionMatrix()->Mult(x, X);
172 fespace.GetProlongationMatrix()->Mult(X, x);
173
174 // Compute and print the L^2 norm of the error.
175 const real_t error = x.ComputeL2Error(vec_coef);
176 cout << "\n|| E_h - E ||_{L^2} = " << error << '\n' << endl;
177 }
178 else
179 {
180 // 6. Determine the list of essential boundary dofs. In this example, the
181 // boundary conditions are defined by marking all the boundary attributes
182 // from the mesh as essential (Dirichlet) and converting them to a list of
183 // true dofs.
184 Array<int> ess_tdof_list;
185 if (mesh.bdr_attributes.Size())
186 {
187 Array<int> ess_bdr(mesh.bdr_attributes.Max());
188 ess_bdr = 1;
189 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
190 }
191
192 // 7. Set up the linear form b(.) which corresponds to the right-hand side of
193 // the FEM linear system, which in this case is (1,phi_i) where phi_i are
194 // the basis functions in fespace.
195 LinearForm b(&fespace);
196 ConstantCoefficient one(1.0);
197 b.AddDomainIntegrator(new DomainLFIntegrator(one));
198 b.Assemble();
199
200 // 8. Define the solution vector x as a finite element grid function
201 // corresponding to fespace. Initialize x with initial guess of zero,
202 // which satisfies the boundary conditions.
203 x = 0.0;
204
205 // 9. Set up the bilinear form a(.,.) on the finite element space
206 // corresponding to the Laplacian operator -Delta, by adding the diffusion
207 // domain integrator.
208 BilinearForm a(&fespace);
209 a.AddDomainIntegrator(new DiffusionIntegrator(one));
210
211 // 10. Assemble the bilinear form and the corresponding linear system,
212 // applying any necessary transformations such as: assembly, eliminating
213 // boundary conditions, applying conforming constraints for non-conforming
214 // AMR, static condensation, etc.
215 a.Assemble();
216
217 OperatorPtr A;
218 Vector B;
219 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
220
221 // 11. Solve the linear system A X = B.
222 {
223#ifndef MFEM_USE_SUITESPARSE
224 // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
225 GSSmoother M((SparseMatrix&)(*A));
226 PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
227#else
228 // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
229 UMFPackSolver umf_solver;
230 umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
231 umf_solver.SetOperator(*A);
232 umf_solver.Mult(B, X);
233#endif
234 }
235
236 // 12. Recover the grid function corresponding to X.
237 a.RecoverFEMSolution(X, b, x);
238 }
239
240 if (fespaceDim == 1)
241 {
242 const real_t h1error = CheckH1Continuity(x);
243 cout << "H1 continuity error " << h1error << endl;
244 MFEM_VERIFY(h1error < 1.0e-12, "H1 continuity is not satisfied");
245 }
246
247 L2_FECollection fecL2(0, dim);
248 FiniteElementSpace l2fespace(&mesh, &fecL2);
249 GridFunction xo(&l2fespace);
250 xo = 0.0;
251
252 for (int e=0; e<mesh.GetNE(); ++e)
253 {
254 const int p_elem = fespace.GetElementOrder(e);
255 Array<int> dofs;
256 l2fespace.GetElementDofs(e, dofs);
257 MFEM_VERIFY(dofs.Size() == 1, "");
258 xo[dofs[0]] = p_elem;
259 }
260
261 // 13. Save the refined mesh and the solution. This output can be viewed later
262 // using GLVis: "glvis -m refined.mesh -g sol.gf".
263 std::unique_ptr<GridFunction> vis_x = x.ProlongateToMaxOrder();
264 ofstream mesh_ofs("refined.mesh");
265 mesh_ofs.precision(8);
266 mesh.Print(mesh_ofs);
267 ofstream sol_ofs("sol.gf");
268 sol_ofs.precision(8);
269 vis_x->Save(sol_ofs);
270 ofstream order_ofs("order.gf");
271 order_ofs.precision(8);
272 xo.Save(order_ofs);
273
274 // 14. Send the solution by socket to a GLVis server.
275 if (visualization)
276 {
277 char vishost[] = "localhost";
278 int visport = 19916;
279 socketstream sol_sock(vishost, visport);
280 sol_sock.precision(8);
281 sol_sock << "solution\n" << mesh << *vis_x << flush;
282 }
283
284 // 15. Free the used memory.
285 if (delete_fec)
286 {
287 delete fec;
288 }
289
290 return 0;
291}
292
294{
295 const FiniteElementSpace *fes = x.FESpace();
296 Mesh *mesh = fes->GetMesh();
297
298 const int dim = mesh->Dimension();
299
300 // Following the example of KellyErrorEstimator::ComputeEstimates(),
301 // we loop over interior faces and compute their error contributions.
302 real_t errorMax = 0.0;
303 for (int f = 0; f < mesh->GetNumFaces(); f++)
304 {
305 if (mesh->FaceIsInterior(f))
306 {
307 int Inf1, Inf2, NCFace;
308 mesh->GetFaceInfos(f, &Inf1, &Inf2, &NCFace);
309
310 auto FT = mesh->GetFaceElementTransformations(f);
311
312 const int faceOrder = dim == 3 ? fes->GetFaceOrder(f) :
313 fes->GetEdgeOrder(f);
314 auto &int_rule = IntRules.Get(FT->FaceGeom, 2 * faceOrder);
315 const auto nip = int_rule.GetNPoints();
316
317 // Convention
318 // * Conforming face: Face side with smaller element id handles
319 // the integration
320 // * Non-conforming face: The slave handles the integration.
321 // See FaceInfo documentation for details.
322 bool isNCSlave = FT->Elem2No >= 0 && NCFace >= 0;
323 bool isConforming = FT->Elem2No >= 0 && NCFace == -1;
324 if ((FT->Elem1No < FT->Elem2No && isConforming) || isNCSlave)
325 {
326 for (int i = 0; i < nip; i++)
327 {
328 const auto &fip = int_rule.IntPoint(i);
330
331 FT->Loc1.Transform(fip, ip);
332 const real_t v1 = x.GetValue(FT->Elem1No, ip);
333
334 FT->Loc2.Transform(fip, ip);
335 const real_t v2 = x.GetValue(FT->Elem2No, ip);
336
337 const real_t err_i = std::abs(v1 - v2);
338 errorMax = std::max(errorMax, err_i);
339 }
340 }
341 }
342 }
343
344 return errorMax;
345}
346
347void f_exact(const Vector &x, Vector &f)
348{
349 constexpr real_t freq = 1.0;
350 constexpr real_t kappa = freq * M_PI;
351
352 if (x.Size() == 3)
353 {
354 f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
355 f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
356 f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
357 }
358 else
359 {
360 f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
361 f(1) = (1. + kappa * kappa) * sin(kappa * x(0));
362 if (x.Size() == 3) { f(2) = 0.0; }
363 }
364}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
int Size() const
Return the logical size of the array.
Definition array.hpp:147
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:830
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
A coefficient that is constant across space and time.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:297
Class for domain integration .
Definition lininteg.hpp:106
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
virtual const char * Name() const
Definition fe_coll.hpp:79
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:851
DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
Definition fespace.cpp:3513
int GetEdgeOrder(int edge, int variant=0) const
Definition fespace.cpp:3356
int GetFaceOrder(int face, int variant=0) const
Returns the polynomial degree of the i'th face finite element.
Definition fespace.cpp:3372
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:658
virtual const Operator * GetProlongationMatrix() const
Definition fespace.hpp:727
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition fespace.cpp:4157
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
Definition fespace.cpp:221
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
virtual void PRefineAndUpdate(const Array< pRefinement > &refs, bool want_transfer=true)
Definition fespace.cpp:4276
virtual const SparseMatrix * GetHpRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.hpp:754
virtual int GetMaxElementOrder() const
Return the maximum polynomial order over all elements.
Definition fespace.hpp:705
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual real_t GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition gridfunc.cpp:446
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Returns ||exsol - u_h||_L2 for scalar or vector H1 or L2 elements.
FiniteElementSpace * FESpace()
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
std::unique_ptr< GridFunction > ProlongateToMaxOrder() const
Return a GridFunction with the values of this, prolongated to the maximum order of all elements in th...
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:275
Class for integration point with weight.
Definition intrules.hpp:35
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:346
Vector with associated FE space and LinearFormIntegrators.
Mesh data type.
Definition mesh.hpp:64
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition mesh.cpp:1004
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:290
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
Definition mesh.cpp:1463
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition mesh.cpp:6531
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition mesh.cpp:10883
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2433
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
Definition mesh.hpp:1462
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, real_t sx=1.0, real_t sy=1.0, real_t sz=1.0, bool sfc_ordering=true)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz hexahedra if type =...
Definition mesh.cpp:4481
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:9294
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition mesh.cpp:10951
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, real_t sx=1.0, real_t sy=1.0, bool sfc_ordering=true)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny quadrilaterals if type = QUADRILATER...
Definition mesh.cpp:4471
Pointer to an Operator of a specified type.
Definition handle.hpp:34
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
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.
Data type sparse matrix.
Definition sparsemat.hpp:51
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication.
Direct sparse solver using UMFPACK.
Definition solvers.hpp:1121
real_t Control[UMFPACK_CONTROL]
Definition solvers.hpp:1131
void SetOperator(const Operator &op) override
Factorize the given Operator op which must be a SparseMatrix.
Definition solvers.cpp:3218
void Mult(const Vector &b, Vector &x) const override
Direct solution of the linear system using UMFPACK.
Definition solvers.cpp:3313
A general vector function coefficient.
Vector data type.
Definition vector.hpp:82
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
real_t kappa
Definition ex24.cpp:54
int dim
Definition ex24.cpp:53
real_t freq
Definition ex24.cpp:54
int main()
int DetRand(int &seed)
Definition hpref.cpp:30
real_t CheckH1Continuity(GridFunction &x)
Definition hpref.cpp:293
void f_exact(const Vector &x, Vector &f)
Definition hpref.cpp:347
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:949
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:492
const char vishost[]