MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
ref321.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11//
12// -----------------------------------------------------------------
13// 3:1 Refinement Miniapp: Perform 3:1 anisotropic mesh refinements
14// -----------------------------------------------------------------
15//
16// This miniapp performs random 3:1 refinements of a quadrilateral or hexahedral
17// mesh. A diffusion equation is solved in an H1 finite element space defined on
18// the refined mesh, and its continuity is verified.
19//
20// Compile with: make ref321
21//
22// Sample runs: ref321 -mm -dim 2 -o 2 -r 100
23// ref321 -mm -dim 3 -o 2 -r 100
24// ref321 -m ../../data/star.mesh -o 2 -r 100
25
26#include "mfem.hpp"
27#include <fstream>
28#include <iostream>
29
30using namespace std;
31using namespace mfem;
32
34
35// Find the two children of parent element `elem` after its refinement in one
36// direction.
37void FindChildren(const Mesh & mesh, int elem, Array<int> & children)
38{
40 MFEM_ASSERT(mesh.GetNE() == cf.embeddings.Size(), "");
41
42 // Note that row `elem` of the table constructed by cf.MakeCoarseToFineTable
43 // is an alternative to this global loop, but constructing the table is also
44 // a global operation with global storage.
45 for (int i = 0; i < mesh.GetNE(); i++)
46 {
47 const int p = cf.embeddings[i].parent;
48 if (p == elem)
49 {
50 children.Append(i);
51 }
52 }
53}
54
55// Refine 3:1 via 2 refinements with scalings 2/3 and 1/2.
56void Refine31(Mesh & mesh, int elem, char type)
57{
58 Array<Refinement> refs; // Refinement is defined in ncmesh.hpp
59 refs.Append(Refinement(elem, type, 2.0 / 3.0));
60 mesh.GeneralRefinement(refs);
61
62 // Find the elements with parent `elem`
63 Array<int> children;
64 FindChildren(mesh, elem, children);
65 MFEM_ASSERT(children.Size() == 2, "");
66
67 const int elem1 = children[0];
68
69 refs.SetSize(0);
70 refs.Append(Refinement(elem1, type)); // Default scaling of 0.5
71 mesh.GeneralRefinement(refs);
72}
73
74// Deterministic, somewhat random integer generator
75int MyRand(int & s)
76{
77 s++;
78 const double a = 1000 * sin(s * 1.1234 * M_PI);
79 return int(std::abs(a));
80}
81
82// Randomly select elements for 3:1 refinements in random directions.
83void TestAnisoRefRandom(int iter, int dim, Mesh & mesh)
84{
85 int seed = 0;
86 for (int i = 0; i < iter; i++)
87 {
88 const int elem = MyRand(seed) % mesh.GetNE();
89 const int t = MyRand(seed) % dim;
90 auto type = t == 0 ? Refinement::X :
91 (t == 1 ? Refinement::Y : Refinement::Z);
92 Refine31(mesh, elem, type);
93 }
94
95 mesh.EnsureNodes();
96 mesh.SetScaledNCMesh();
97}
98
99int main(int argc, char *argv[])
100{
101 // 1. Parse command-line options.
102 const char *mesh_file = "../../data/star.mesh";
103 int order = 1;
104 bool visualization = true;
105 bool makeMesh = false;
106 int num_refs = 1;
107 int tdim = 2; // Mesh dimension
108
109 OptionsParser args(argc, argv);
110 args.AddOption(&mesh_file, "-m", "--mesh",
111 "Mesh file to use.");
112 args.AddOption(&order, "-o", "--order",
113 "Finite element order (polynomial degree) or -1 for"
114 " isoparametric space.");
115 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
116 "--no-visualization",
117 "Enable or disable GLVis visualization.");
118 args.AddOption(&makeMesh, "-mm", "--make-mesh", "-no-mm",
119 "--no-make-mesh", "Create Cartesian mesh");
120 args.AddOption(&tdim, "-dim", "--dimension", "Dimension for Cartesian mesh");
121 args.AddOption(&num_refs, "-r", "--refs", "Number of 3:1 refinements");
122 args.Parse();
123 if (!args.Good())
124 {
125 args.PrintUsage(cout);
126 return 1;
127 }
128 args.PrintOptions(cout);
129
130 // 2. Create or read the mesh from the given mesh file.
131 Mesh mesh;
132 if (makeMesh)
133 {
134 mesh = tdim == 3 ? Mesh::MakeCartesian3D(2, 2, 2, Element::HEXAHEDRON) :
136 }
137 else
138 {
139 mesh = Mesh::LoadFromFile(mesh_file, 1, 1);
140 }
141
142 const int dim = mesh.Dimension();
143
144 // 3. Randomly perform 3:1 refinements in the mesh.
145 TestAnisoRefRandom(num_refs, tdim, mesh);
146
147 // 4. Define a finite element space on the mesh. Here we use continuous
148 // Lagrange finite elements of the specified order.
149 H1_FECollection fec(order, dim);
150 FiniteElementSpace fespace(&mesh, &fec);
151 cout << "Number of finite element unknowns: "
152 << fespace.GetTrueVSize() << endl;
153
154 // 5. Define the solution vector x as a finite element grid function
155 // corresponding to fespace. Solve the Poisson problem, as in ex1.
156 GridFunction x(&fespace);
157
158 {
159 x = 0.0;
160 LinearForm b(&fespace);
161 ConstantCoefficient one(1.0);
162 b.AddDomainIntegrator(new DomainLFIntegrator(one));
163 b.Assemble();
164
165 BilinearForm a(&fespace);
166 a.AddDomainIntegrator(new DiffusionIntegrator());
167 a.Assemble();
168
169 OperatorPtr A;
170 Vector B, X;
171 Array<int> ess_tdof_list;
172 if (mesh.bdr_attributes.Size())
173 {
174 Array<int> ess_bdr(mesh.bdr_attributes.Max());
175 ess_bdr = 1;
176 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
177 }
178
179 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
180
181 GSSmoother M((SparseMatrix&)(*A));
182 PCG(*A, M, B, X, 1, 2000, 1e-12, 0.0);
183 a.RecoverFEMSolution(X, b, x);
184 }
185
186 // 6. Verify the continuity of the projected function in H1.
187 const real_t h1err = CheckH1Continuity(x);
188 cout << "Error of H1 continuity: " << h1err << endl;
189 MFEM_VERIFY(h1err < 1.0e-7, "");
190
191 // 7. Save the refined mesh and the solution. This output can be viewed later
192 // using GLVis: "glvis -m ref321.mesh -g sol.gf".
193 ofstream mesh_ofs("ref321.mesh");
194 mesh_ofs.precision(8);
195 mesh.Print(mesh_ofs);
196 ofstream sol_ofs("sol.gf");
197 sol_ofs.precision(8);
198 x.Save(sol_ofs);
199
200 // 8. Send the solution by socket to a GLVis server.
201 if (visualization)
202 {
203 char vishost[] = "localhost";
204 int visport = 19916;
205 socketstream sol_sock(vishost, visport);
206 sol_sock.precision(8);
207 sol_sock << "solution\n" << mesh << x << flush;
208 }
209
210 return 0;
211}
212
214{
215 const FiniteElementSpace *fes = x.FESpace();
216 Mesh *mesh = fes->GetMesh();
217
218 const int dim = mesh->Dimension();
219
220 // Following the example of KellyErrorEstimator::ComputeEstimates(), we loop
221 // over interior faces and then shared faces.
222
223 // Compute error contribution from local interior faces
224 real_t errorMax = 0.0;
225 for (int f = 0; f < mesh->GetNumFaces(); f++)
226 {
227 if (mesh->FaceIsInterior(f))
228 {
229 int Inf1, Inf2, NCFace;
230 mesh->GetFaceInfos(f, &Inf1, &Inf2, &NCFace);
231
232 auto FT = mesh->GetFaceElementTransformations(f);
233
234 const int faceOrder = dim == 3 ? fes->GetFaceOrder(f) :
235 fes->GetEdgeOrder(f);
236 auto &int_rule = IntRules.Get(FT->FaceGeom, 2 * faceOrder);
237 const auto nip = int_rule.GetNPoints();
238
239 // Convention:
240 // * Conforming face: Face side with smaller element id handles the
241 // integration
242 // * Non-conforming face: The slave handles the integration.
243 // See FaceInfo documentation for details.
244 bool isNCSlave = FT->Elem2No >= 0 && NCFace >= 0;
245 bool isConforming = FT->Elem2No >= 0 && NCFace == -1;
246 if ((FT->Elem1No < FT->Elem2No && isConforming) || isNCSlave)
247 {
248 for (int i = 0; i < nip; i++)
249 {
250 const auto &fip = int_rule.IntPoint(i);
252
253 FT->Loc1.Transform(fip, ip);
254 const real_t v1 = x.GetValue(FT->Elem1No, ip);
255
256 FT->Loc2.Transform(fip, ip);
257 const real_t v2 = x.GetValue(FT->Elem2No, ip);
258
259 const real_t err_i = std::abs(v1 - v2);
260 errorMax = std::max(errorMax, err_i);
261 }
262 }
263 }
264 }
265
266 return errorMax;
267}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:758
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.
Class for domain integration .
Definition lininteg.hpp:106
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
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
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
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.
FiniteElementSpace * FESpace()
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.
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
void EnsureNodes()
Make sure that the mesh has valid nodes, i.e. its geometry is described by a vector finite element gr...
Definition mesh.cpp:6432
void SetScaledNCMesh()
Definition mesh.hpp:2371
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
static Mesh LoadFromFile(const std::string &filename, int generate_edges=0, int refine=1, bool fix_orientation=true)
Definition mesh.cpp:4453
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition mesh.hpp:299
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
const CoarseFineTransformations & GetRefinementTransforms() const
Definition ncmesh.cpp:5153
Pointer to an Operator of a specified type.
Definition handle.hpp:34
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
Vector data type.
Definition vector.hpp:82
int dim
Definition ex24.cpp:53
int main()
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[]
real_t p(const Vector &x, real_t t)
void FindChildren(const Mesh &mesh, int elem, Array< int > &children)
Definition ref321.cpp:37
real_t CheckH1Continuity(GridFunction &x)
Definition ref321.cpp:213
void Refine31(Mesh &mesh, int elem, char type)
Definition ref321.cpp:56
int MyRand(int &s)
Definition ref321.cpp:75
void TestAnisoRefRandom(int iter, int dim, Mesh &mesh)
Definition ref321.cpp:83
Defines the coarse-fine transformations of all fine elements.
Definition ncmesh.hpp:90
Array< Embedding > embeddings
Fine element positions in their parents.
Definition ncmesh.hpp:92