MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
mesh_operators.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#include "mesh_operators.hpp"
13#include "pmesh.hpp"
14
15namespace mfem
16{
17
19{
20 // delete in reverse order
21 for (int i = sequence.Size()-1; i >= 0; i--)
22 {
23 delete sequence[i];
24 }
25}
26
28{
29 if (sequence.Size() == 0) { return NONE; }
30next_step:
31 step = (step + 1) % sequence.Size();
32 bool last = (step == sequence.Size() - 1);
33 int mod = sequence[step]->ApplyImpl(mesh);
34 switch (mod & MASK_ACTION)
35 {
36 case NONE: if (last) { return NONE; } goto next_step;
37 case CONTINUE: return last ? mod : (REPEAT | (mod & MASK_INFO));
38 case STOP: return STOP;
39 case REPEAT: --step; return mod;
40 }
41 return NONE;
42}
43
45{
46 for (int i = 0; i < sequence.Size(); i++)
47 {
48 sequence[i]->Reset();
49 }
50 step = 0;
51}
52
53
55 : estimator(est)
56{
59 total_err_goal = 0.0;
60 total_fraction = 0.5;
61 local_err_goal = 0.0;
62 max_elements = std::numeric_limits<long long>::max();
63
64 threshold = 0.0;
67
68 non_conforming = -1;
69 nc_limit = 0;
70}
71
72real_t ThresholdRefiner::GetNorm(const Vector &local_err, Mesh &mesh) const
73{
74#ifdef MFEM_USE_MPI
75 ParMesh *pmesh = dynamic_cast<ParMesh*>(&mesh);
76 if (pmesh)
77 {
78 return ParNormlp(local_err, total_norm_p, pmesh->GetComm());
79 }
80#endif
81 return local_err.Normlp(total_norm_p);
82}
83
85 Array<Refinement> & refinements)
86{
87 threshold = 0.0;
89 refinements.SetSize(0);
91
92 const long long num_elements = mesh.GetGlobalNE();
93 if (num_elements >= max_elements) { return STOP; }
94
95 const int NE = mesh.GetNE();
96 const Vector &local_err = estimator.GetLocalErrors();
97 MFEM_ASSERT(local_err.Size() == NE, "invalid size of local_err");
98
99 const real_t total_err = GetNorm(local_err, mesh);
100 if (total_err <= total_err_goal) { return STOP; }
101
102 if (total_norm_p < infinity())
103 {
104 threshold = std::max((real_t) (total_err * total_fraction *
105 std::pow(num_elements, -1.0/total_norm_p)),
107 }
108 else
109 {
110 threshold = std::max(total_err * total_fraction, local_err_goal);
111 }
112
113 for (int el = 0; el < NE; el++)
114 {
115 if (local_err(el) > threshold)
116 {
117 refinements.Append(Refinement(el));
118 }
119 }
120
121 if (aniso_estimator)
122 {
123 const Array<int> &aniso_flags = aniso_estimator->GetAnisotropicFlags();
124 if (aniso_flags.Size() > 0)
125 {
126 for (int i = 0; i < refinements.Size(); i++)
127 {
128 Refinement &ref = refinements[i];
129 ref.SetType(aniso_flags[ref.index]);
130 }
131 }
132 }
133
134 return NONE;
135}
136
138{
139 const int action = MarkWithoutRefining(mesh, marked_elements);
140 if (action == STOP) { return STOP; }
141
143 if (num_marked_elements == 0LL) { return STOP; }
144
146 return CONTINUE + REFINED;
147}
148
150{
152 current_sequence = -1;
154 // marked_elements.SetSize(0); // not necessary
155}
156
157
159{
160 if (mesh.Conforming()) { return NONE; }
161
162 const Vector &local_err = estimator.GetLocalErrors();
163 bool derefs = mesh.DerefineByError(local_err, threshold, nc_limit, op);
164
165 return derefs ? CONTINUE + DEREFINED : NONE;
166}
167
168
170{
171 int max_it = 1;
172 return PreprocessMesh(mesh, max_it);
173}
174
176{
177 int rank = 0;
178 MFEM_VERIFY(max_it > 0, "max_it must be strictly positive")
179
180 int dim = mesh.Dimension();
181 L2_FECollection l2fec(order, dim);
182 FiniteElementSpace* l2fes = NULL;
183
184 bool par = false;
185 GridFunction *gf = NULL;
186
187#ifdef MFEM_USE_MPI
188 ParMesh* pmesh = dynamic_cast<ParMesh*>(&mesh);
189 if (pmesh && pmesh->Nonconforming())
190 {
191 par = true;
192 l2fes = new ParFiniteElementSpace(pmesh, &l2fec);
193 gf = new ParGridFunction(static_cast<ParFiniteElementSpace*>(l2fes));
194 }
195#endif
196 if (!par)
197 {
198 l2fes = new FiniteElementSpace(&mesh, &l2fec);
199 gf = new GridFunction(l2fes);
200 }
201
202 // If custom integration rule has not been set,
203 // then use the default integration rule
204 if (!irs)
205 {
206 int order_quad = 2*order + 3;
207 for (int i=0; i < Geometry::NumGeom; ++i)
208 {
209 ir_default[i] = &(IntRules.Get(i, order_quad));
210 }
211 irs = ir_default;
212 }
213
214 for (int i = 0; i < max_it; i++)
215 {
216 // Compute number of elements and L2-norm of f.
217 int NE = mesh.GetNE();
218 int globalNE = 0;
219 real_t norm_of_coeff = 0.0;
220 if (par)
221 {
222#ifdef MFEM_USE_MPI
223 globalNE = pmesh->GetGlobalNE();
224 norm_of_coeff = ComputeGlobalLpNorm(2.0,*coeff,*pmesh,irs);
225#endif
226 }
227 else
228 {
229 globalNE = NE;
230 norm_of_coeff = ComputeLpNorm(2.0,*coeff,mesh,irs);
231 }
232
233 // Compute average L2-norm of f
234 real_t av_norm_of_coeff = norm_of_coeff / sqrt(globalNE);
235
236 // Compute element-wise L2-norms of (I - Π) f
237 Vector element_norms_of_fine_scale(NE);
238 gf->SetSpace(l2fes);
240 gf->ComputeElementL2Errors(*coeff,element_norms_of_fine_scale,irs);
241
242 // Define osc_K(f) := || h ⋅ (I - Π) f ||_K and select elements
243 // for refinement based on threshold. Also record relative osc(f).
244 global_osc = 0.0;
248 element_oscs = 0.0;
249 for (int j = 0; j < NE; j++)
250 {
251 real_t h = mesh.GetElementSize(j);
252 real_t element_osc = h * element_norms_of_fine_scale(j);
253 if ( element_osc > threshold * av_norm_of_coeff )
254 {
256 }
257 element_oscs(j) = element_osc/(norm_of_coeff + 1e-10);
258 global_osc += element_osc*element_osc;
259 }
260#ifdef MFEM_USE_MPI
261 if (par)
262 {
263 MPI_Comm comm = pmesh->GetComm();
264 MPI_Allreduce(MPI_IN_PLACE, &global_osc, 1, MPITypeMap<real_t>::mpi_type,
265 MPI_SUM, comm);
266 MPI_Comm_rank(comm, &rank);
267 }
268#endif
269 global_osc = sqrt(global_osc)/(norm_of_coeff + 1e-10);
270
271 // Exit if the global threshold or maximum number of elements is reached.
273 {
274 if (global_osc > threshold && globalNE > max_elements && rank == 0 &&
276 {
277 MFEM_WARNING("Reached maximum number of elements "
278 "before resolving data to tolerance.");
279 }
280 delete l2fes;
281 delete gf;
282 return STOP;
283 }
284
285 // Refine elements.
287 l2fes->Update(false);
288 gf->Update();
289
290 }
291 delete l2fes;
292 delete gf;
293 return CONTINUE + REFINED;
294
295}
296
298{
300 global_osc = 0.0;
301 coeff = NULL;
302 irs = NULL;
303}
304
305
307{
308#ifdef MFEM_USE_MPI
309 ParMesh *pmesh = dynamic_cast<ParMesh*>(&mesh);
310 if (pmesh && pmesh->Nonconforming())
311 {
312 pmesh->Rebalance();
313 return CONTINUE + REBALANCED;
314 }
315#endif
316 return NONE;
317}
318
319
320} // namespace mfem
The AnisotropicErrorEstimator class is the base class for all error estimators that compute one non-n...
virtual const Array< int > & GetAnisotropicFlags()=0
Get an Array<int> with anisotropic flags for all mesh elements.
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
const IntegrationRule * ir_default[Geometry::NumGeom]
virtual int PreprocessMesh(Mesh &mesh, int max_it)
Apply the operator to the mesh max_it times or until tolerance achieved.
const IntegrationRule ** irs
void Reset() override
Reset.
int ApplyImpl(Mesh &mesh) override
Apply the operator to the mesh once.
Base class for all element based error estimators.
virtual void Reset()=0
Force recomputation of the estimates on the next call to GetLocalErrors.
virtual const Vector & GetLocalErrors()=0
Get a Vector with all element errors.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
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
static const int NumGeom
Definition geom.hpp:46
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void ComputeElementL2Errors(Coefficient &exsol, Vector &error, const IntegrationRule *irs[]=NULL) const
Returns ||u_ex - u_h||_L2 elementwise for H1 or L2 elements.
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition gridfunc.cpp:167
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition gridfunc.cpp:225
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
Array< MeshOperator * > sequence
MeshOperators sequence, owned by us.
int ApplyImpl(Mesh &mesh) override
Apply the MeshOperatorSequence.
virtual ~MeshOperatorSequence()
Delete all operators from the sequence.
void Reset() override
Reset all MeshOperators in the sequence.
@ MASK_ACTION
bit mask for all "action" bits
@ STOP
a stopping criterion was satisfied
@ MASK_INFO
bit mask for all "info" bits
@ REFINED
the mesh was refined
@ REBALANCED
the mesh was rebalanced
@ DEREFINED
the mesh was de-refined
Mesh data type.
Definition mesh.hpp:64
virtual long long ReduceInt(int value) const
Utility function: sum integers from all processors (Allreduce).
Definition mesh.hpp:2589
bool Conforming() const
Return a bool indicating whether this mesh is conforming.
Definition mesh.hpp:2365
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition mesh.cpp:10883
bool Nonconforming() const
Return a bool indicating whether this mesh is nonconforming.
Definition mesh.hpp:2367
long GetSequence() const
Definition mesh.hpp:2387
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
real_t GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
Definition mesh.cpp:107
bool DerefineByError(Array< real_t > &elem_error, real_t threshold, int nc_limit=0, int op=1)
Definition mesh.cpp:10629
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition mesh.hpp:1311
Abstract parallel finite element space.
Definition pfespace.hpp:29
Class for parallel grid function.
Definition pgridfunc.hpp:50
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
void Rebalance()
Definition pmesh.cpp:4008
int ApplyImpl(Mesh &mesh) override
Rebalance a parallel mesh (only non-conforming parallel meshes are supported).
int ApplyImpl(Mesh &mesh) override
Apply the operator to the mesh.
Array< Refinement > marked_elements
ThresholdRefiner(ErrorEstimator &est)
Construct a ThresholdRefiner using the given ErrorEstimator.
AnisotropicErrorEstimator * aniso_estimator
void Reset() override
Reset the associated estimator.
int MarkWithoutRefining(Mesh &mesh, Array< Refinement > &refinements)
Set the array refinements of elements to refine, without refining.
real_t GetNorm(const Vector &local_err, Mesh &mesh) const
ErrorEstimator & estimator
int ApplyImpl(Mesh &mesh) override
Apply the operator to the mesh.
Vector data type.
Definition vector.hpp:82
void Destroy()
Destroy a vector.
Definition vector.hpp:635
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 Normlp(real_t p) const
Returns the l_p norm of the vector.
Definition vector.cpp:998
int dim
Definition ex24.cpp:53
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition vector.hpp:47
real_t ParNormlp(const Vector &vec, real_t p, MPI_Comm comm)
Compute the l_p norm of the Vector which is split without overlap across the given communicator.
Definition hypre.cpp:479
real_t ComputeGlobalLpNorm(real_t p, Coefficient &coeff, ParMesh &pmesh, const IntegrationRule *irs[])
Compute the global Lp norm of a function f. .
real_t ComputeLpNorm(real_t p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
Compute the Lp norm of a function f. .
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:492
Helper struct to convert a C++ type to an MPI type.
int index
Mesh element number.
Definition ncmesh.hpp:42
void SetType(char type, real_t scale=0.5)
Set the type and scale, assuming the element is already set.
Definition ncmesh.cpp:596