MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
mesh_operators.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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{
86 threshold = 0.0;
88 marked_elements.SetSize(0);
90
91 const long long num_elements = mesh.GetGlobalNE();
92 if (num_elements >= max_elements) { return STOP; }
93
94 const int NE = mesh.GetNE();
95 const Vector &local_err = estimator.GetLocalErrors();
96 MFEM_ASSERT(local_err.Size() == NE, "invalid size of local_err");
97
98 const real_t total_err = GetNorm(local_err, mesh);
99 if (total_err <= total_err_goal) { return STOP; }
100
101 if (total_norm_p < infinity())
102 {
103 threshold = std::max((real_t) (total_err * total_fraction *
104 std::pow(num_elements, -1.0/total_norm_p)),
106 }
107 else
108 {
109 threshold = std::max(total_err * total_fraction, local_err_goal);
110 }
111
112 for (int el = 0; el < NE; el++)
113 {
114 if (local_err(el) > threshold)
115 {
116 marked_elements.Append(Refinement(el));
117 }
118 }
119
120 if (aniso_estimator)
121 {
122 const Array<int> &aniso_flags = aniso_estimator->GetAnisotropicFlags();
123 if (aniso_flags.Size() > 0)
124 {
125 for (int i = 0; i < marked_elements.Size(); i++)
126 {
127 Refinement &ref = marked_elements[i];
128 ref.ref_type = aniso_flags[ref.index];
129 }
130 }
131 }
132
134 if (num_marked_elements == 0LL) { return STOP; }
135
137 return CONTINUE + REFINED;
138}
139
141{
143 current_sequence = -1;
145 // marked_elements.SetSize(0); // not necessary
146}
147
148
150{
151 if (mesh.Conforming()) { return NONE; }
152
153 const Vector &local_err = estimator.GetLocalErrors();
154 bool derefs = mesh.DerefineByError(local_err, threshold, nc_limit, op);
155
156 return derefs ? CONTINUE + DEREFINED : NONE;
157}
158
159
161{
162 int max_it = 1;
163 return PreprocessMesh(mesh, max_it);
164}
165
167{
168 int rank = 0;
169 MFEM_VERIFY(max_it > 0, "max_it must be strictly positive")
170
171 int dim = mesh.Dimension();
172 L2_FECollection l2fec(order, dim);
173 FiniteElementSpace* l2fes = NULL;
174
175 bool par = false;
176 GridFunction *gf = NULL;
177
178#ifdef MFEM_USE_MPI
179 ParMesh* pmesh = dynamic_cast<ParMesh*>(&mesh);
180 if (pmesh && pmesh->Nonconforming())
181 {
182 par = true;
183 l2fes = new ParFiniteElementSpace(pmesh, &l2fec);
184 gf = new ParGridFunction(static_cast<ParFiniteElementSpace*>(l2fes));
185 }
186#endif
187 if (!par)
188 {
189 l2fes = new FiniteElementSpace(&mesh, &l2fec);
190 gf = new GridFunction(l2fes);
191 }
192
193 // If custom integration rule has not been set,
194 // then use the default integration rule
195 if (!irs)
196 {
197 int order_quad = 2*order + 3;
198 for (int i=0; i < Geometry::NumGeom; ++i)
199 {
200 ir_default[i] = &(IntRules.Get(i, order_quad));
201 }
202 irs = ir_default;
203 }
204
205 for (int i = 0; i < max_it; i++)
206 {
207 // Compute number of elements and L2-norm of f.
208 int NE = mesh.GetNE();
209 int globalNE = 0;
210 real_t norm_of_coeff = 0.0;
211 if (par)
212 {
213#ifdef MFEM_USE_MPI
214 globalNE = pmesh->GetGlobalNE();
215 norm_of_coeff = ComputeGlobalLpNorm(2.0,*coeff,*pmesh,irs);
216#endif
217 }
218 else
219 {
220 globalNE = NE;
221 norm_of_coeff = ComputeLpNorm(2.0,*coeff,mesh,irs);
222 }
223
224 // Compute average L2-norm of f
225 real_t av_norm_of_coeff = norm_of_coeff / sqrt(globalNE);
226
227 // Compute element-wise L2-norms of (I - Π) f
228 Vector element_norms_of_fine_scale(NE);
229 gf->SetSpace(l2fes);
231 gf->ComputeElementL2Errors(*coeff,element_norms_of_fine_scale,irs);
232
233 // Define osc_K(f) := || h ⋅ (I - Π) f ||_K and select elements
234 // for refinement based on threshold. Also record relative osc(f).
235 global_osc = 0.0;
239 element_oscs = 0.0;
240 for (int j = 0; j < NE; j++)
241 {
242 real_t h = mesh.GetElementSize(j);
243 real_t element_osc = h * element_norms_of_fine_scale(j);
244 if ( element_osc > threshold * av_norm_of_coeff )
245 {
247 }
248 element_oscs(j) = element_osc/(norm_of_coeff + 1e-10);
249 global_osc += element_osc*element_osc;
250 }
251#ifdef MFEM_USE_MPI
252 if (par)
253 {
254 MPI_Comm comm = pmesh->GetComm();
255 MPI_Allreduce(MPI_IN_PLACE, &global_osc, 1, MPITypeMap<real_t>::mpi_type,
256 MPI_SUM, comm);
257 MPI_Comm_rank(comm, &rank);
258 }
259#endif
260 global_osc = sqrt(global_osc)/(norm_of_coeff + 1e-10);
261
262 // Exit if the global threshold or maximum number of elements is reached.
263 if (global_osc < threshold || globalNE > max_elements)
264 {
265 if (global_osc > threshold && globalNE > max_elements && rank == 0 &&
267 {
268 MFEM_WARNING("Reached maximum number of elements "
269 "before resolving data to tolerance.");
270 }
271 delete l2fes;
272 delete gf;
273 return STOP;
274 }
275
276 // Refine elements.
278 l2fes->Update(false);
279 gf->Update();
280
281 }
282 delete l2fes;
283 delete gf;
284 return CONTINUE + REFINED;
285
286}
287
289{
291 global_osc = 0.0;
292 coeff = NULL;
293 irs = NULL;
294}
295
296
298{
299#ifdef MFEM_USE_MPI
300 ParMesh *pmesh = dynamic_cast<ParMesh*>(&mesh);
301 if (pmesh && pmesh->Nonconforming())
302 {
303 pmesh->Rebalance();
304 return CONTINUE + REBALANCED;
305 }
306#endif
307 return NONE;
308}
309
310
311} // namespace mfem
The AnisotropicErrorEstimator class is the base class for all error estimators that compute one non-n...
Definition: estimators.hpp:65
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:697
int Size() const
Return the logical size of the array.
Definition: array.hpp:144
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition: array.hpp:769
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh once.
virtual void Reset()
Reset.
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
Base class for all element based error estimators.
Definition: estimators.hpp:42
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:220
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3439
static const int NumGeom
Definition: geom.hpp:42
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
Definition: gridfunc.hpp:629
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:164
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2346
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition: gridfunc.cpp:195
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:1006
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:330
Array< MeshOperator * > sequence
MeshOperators sequence, owned by us.
virtual void Reset()
Reset all MeshOperators in the sequence.
virtual int ApplyImpl(Mesh &mesh)
Apply the MeshOperatorSequence.
virtual ~MeshOperatorSequence()
Delete all operators from 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:56
virtual long long ReduceInt(int value) const
Utility function: sum integers from all processors (Allreduce).
Definition: mesh.hpp:2439
bool Conforming() const
Definition: mesh.hpp:2228
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:10558
bool Nonconforming() const
Definition: mesh.hpp:2229
long GetSequence() const
Definition: mesh.hpp:2242
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
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:10304
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:1255
Abstract parallel finite element space.
Definition: pfespace.hpp:29
Class for parallel grid function.
Definition: pgridfunc.hpp:33
Class for parallel meshes.
Definition: pmesh.hpp:34
MPI_Comm GetComm() const
Definition: pmesh.hpp:402
void Rebalance()
Definition: pmesh.cpp:4009
virtual int ApplyImpl(Mesh &mesh)
Rebalance a parallel mesh (only non-conforming parallel meshes are supported).
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh.
ErrorEstimator & estimator
Array< Refinement > marked_elements
ThresholdRefiner(ErrorEstimator &est)
Construct a ThresholdRefiner using the given ErrorEstimator.
virtual int ApplyImpl(Mesh &mesh)
Apply the operator to the mesh.
AnisotropicErrorEstimator * aniso_estimator
virtual void Reset()
Reset the associated estimator.
real_t GetNorm(const Vector &local_err, Mesh &mesh) const
ErrorEstimator & estimator
Vector data type.
Definition: vector.hpp:80
void Destroy()
Destroy a vector.
Definition: vector.hpp:615
int Size() const
Returns the size of the vector.
Definition: vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:538
real_t Normlp(real_t p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:872
int dim
Definition: ex24.cpp:53
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:45
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:450
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:486
Helper struct to convert a C++ type to an MPI type.
char ref_type
refinement XYZ bit mask (7 = full isotropic)
Definition: ncmesh.hpp:40
int index
Mesh element number.
Definition: ncmesh.hpp:39