Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : #ifndef LIBMESH_VARIATIONAL_SMOOTHER_SYSTEM_H
19 : #define LIBMESH_VARIATIONAL_SMOOTHER_SYSTEM_H
20 :
21 : // libMesh includes
22 : #include "libmesh/enum_fe_family.h"
23 : #include "libmesh/enum_quadrature_type.h"
24 : #include "libmesh/fem_function_base.h"
25 : #include "libmesh/fem_system.h"
26 : #include "libmesh/libmesh_common.h"
27 :
28 : // C++ includes
29 : #include <map>
30 : #include <memory>
31 :
32 : namespace libMesh
33 : {
34 :
35 : /**
36 : * Struct to hold smoother-relevant information about the mesh quality.
37 : */
38 216 : struct MeshQualityInfo
39 : {
40 : // dof_id_types will hold the id of the Elem where the metric occurs
41 : // IMPORTANT: the Real should be the first entry of the pair so that taking
42 : // the min/max accross processors will compare the numeric values instead of
43 : // the element ids.
44 :
45 : std::pair<Real, dof_id_type> max_elem_distortion{std::numeric_limits<Real>::lowest(),
46 : DofObject::invalid_id};
47 : std::pair<Real, dof_id_type> min_elem_distortion{std::numeric_limits<Real>::max(),
48 : DofObject::invalid_id};
49 : Real total_distortion = 0.;
50 :
51 : std::pair<Real, dof_id_type> max_elem_dilation{std::numeric_limits<Real>::lowest(),
52 : DofObject::invalid_id};
53 : std::pair<Real, dof_id_type> min_elem_dilation{std::numeric_limits<Real>::max(),
54 : DofObject::invalid_id};
55 : Real total_dilation = 0.;
56 :
57 : std::pair<Real, dof_id_type> max_elem_combined{std::numeric_limits<Real>::lowest(),
58 : DofObject::invalid_id};
59 : std::pair<Real, dof_id_type> min_elem_combined{std::numeric_limits<Real>::max(),
60 : DofObject::invalid_id};
61 : Real total_combined = 0.;
62 :
63 : std::pair<Real, dof_id_type> max_elem_det_S{std::numeric_limits<Real>::lowest(),
64 : DofObject::invalid_id};
65 : std::pair<Real, dof_id_type> min_elem_det_S{std::numeric_limits<Real>::max(),
66 : DofObject::invalid_id};
67 : Real total_det_S = 0.;
68 : Real max_qp_det_S = std::numeric_limits<Real>::lowest();
69 : Real min_qp_det_S = std::numeric_limits<Real>::max();
70 :
71 : bool mesh_is_tangled = false;
72 : bool initialized = false;
73 : };
74 :
75 : // FEMSystem, TimeSolver and NewtonSolver will handle most tasks,
76 : // but we must specify element residuals
77 360 : class VariationalSmootherSystem : public libMesh::FEMSystem
78 : {
79 : /**
80 : * This is an FEMSystem to solve the optimization probelem posed by the
81 : * VariationalMeshSmoother class.
82 : *
83 : * The residual is coded as the gradient of the distortion-dilation metric, and
84 : * the jacobian as analytically coded as the Hessian of the metric.
85 : *
86 : * The nodes of the system mesh are updated during the solve.
87 : */
88 : public:
89 2556 : VariationalSmootherSystem(libMesh::EquationSystems & es,
90 : const std::string & name,
91 : const unsigned int number)
92 2556 : : libMesh::FEMSystem(es, name, number),
93 2412 : _verbosity(0),
94 2412 : _epsilon_squared(TOLERANCE),
95 2412 : _epsilon_squared_assembly(0.),
96 2412 : _ref_vol(0.),
97 2412 : _dilation_weight(0.5),
98 2412 : _untangling_solve(false),
99 2628 : _quadrature_type(QGAUSS)
100 2556 : {}
101 :
102 : // Default destructor
103 : ~VariationalSmootherSystem() override;
104 :
105 : /**
106 : * Assembly method to update the mesh based on the smoother solve.
107 : */
108 : virtual void assembly (bool get_residual,
109 : bool get_jacobian,
110 : bool apply_heterogeneous_constraints = false,
111 : bool apply_no_constraints = false) override;
112 :
113 72 : Real & get_dilation_weight() { return _dilation_weight; }
114 :
115 : /**
116 : * Set the quadrature rule type used to integrate the distortion-dilation
117 : * metric over each element. The default (\p QGAUSS) samples the element
118 : * interior only, which can miss degeneracies localized at element corners
119 : * (e.g. an element collapsing toward one of its nodes). A vertex-sampling
120 : * rule such as \p QTRAP, \p QSIMPSON, \p QNODAL, or \p QGAUSS_LOBATTO
121 : * evaluates the metric at the element nodes, so a folding corner drives the
122 : * metric up (and is flagged as tangled) as it should be. The number of
123 : * points is set from the default order for the mesh's FE type plus the
124 : * System's \p extra_quadrature_order. Not every quadrature type is defined
125 : * for every element type (e.g. \p QGAUSS_LOBATTO is only available for
126 : * tensor-product elements), so choose a type compatible with the mesh.
127 : */
128 2556 : void set_quadrature_type(QuadratureType qt) { _quadrature_type = qt; }
129 :
130 : /**
131 : * Get the quadrature rule type used to integrate the distortion-dilation
132 : * metric. See set_quadrature_type().
133 : */
134 : QuadratureType get_quadrature_type() const { return _quadrature_type; }
135 :
136 : /**
137 : * Solves the system to smooth the mesh. If the mesh is initially tangled,
138 : * a solve is first performed to untangle the mesh, followed by a solve to
139 : * smooth the mesh.
140 : */
141 : virtual void solve() override;
142 :
143 : /**
144 : * Get the jacobians (and determinants) of the target-to-reference element mapping.
145 : * @param target_elem Target element.
146 : * @param femcontext Context used to build mapping.
147 : * @param jacobian Vector in which to store the jacobians for each quadrature point.
148 : * @param jacobian_dets Vector in which to store the determinant of the jacobians
149 : * for each quadrature point.
150 : */
151 : static void get_target_to_reference_jacobian(const Elem * const target_elem,
152 : const FEMContext & femcontext,
153 : std::vector<RealTensor> & jacobians,
154 : std::vector<Real> & jacobian_dets);
155 :
156 : /**
157 : * Getter for the _mesh_info attribute. If this attribute has not yet been
158 : * initialized, compute_mesh_quality_info is called to initialize it.
159 : */
160 : const MeshQualityInfo & get_mesh_info();
161 :
162 : /*
163 : * Computes information about the mesh quality and sets the _mesh_info attribute.
164 : */
165 : void compute_mesh_quality_info();
166 :
167 : /*
168 : * Sets the verbosity of the object.
169 : */
170 2556 : void set_verbosity(const unsigned int verbosity) { _verbosity = verbosity; }
171 :
172 : protected:
173 :
174 : // System initialization
175 : virtual void init_data () override;
176 :
177 : // Context initialization
178 : virtual void init_context (libMesh::DiffContext & context) override;
179 :
180 : // Element residual and jacobian calculations
181 : // Time dependent parts
182 : virtual bool element_time_derivative (bool request_jacobian,
183 : libMesh::DiffContext & context) override;
184 :
185 : /* Computes the element reference volume used in the dilation metric
186 : * The reference value is set to the averaged value of all elements' average
187 : * |J|. Also computes any applicable target element inverse Jacobians. Target
188 : * elements are relavant when the reference element does not minimize the
189 : * distortion metric.
190 : */
191 : void prepare_for_smoothing();
192 :
193 : /**
194 : * Verbosity setting.
195 : * The verbosity levels and the corresponding information output are as
196 : * follows:
197 : *
198 : * verbosity = 0 : No information.
199 : *
200 : * 0 < verbosity : Prints:
201 : * - Initial mesh quality information
202 : * - The untangled mesh quality (if applicable)
203 : * - The smoothed mesh quality
204 : *
205 : * 10 < verbosity: Prints:
206 : * - The reference volume used for the dilation metric
207 : *
208 : * 50 < verbosity: Prints:
209 : * - Mesh quality information after each assembly
210 : *
211 : * 90 < verbosity: Prints:
212 : * - Quality information about each element after each assembly
213 : *
214 : */
215 : unsigned int _verbosity;
216 :
217 : /**
218 : * The small nonzero constant to prevent zero denominators (degenerate meshes only)
219 : */
220 : const Real _epsilon_squared;
221 :
222 : /**
223 : * Epsilon squared value determined at runtime during each assembly. The value
224 : * depends on whether the mesh is tangled.
225 : */
226 : Real _epsilon_squared_assembly;
227 :
228 : /**
229 : * The reference volume for each element
230 : */
231 : Real _ref_vol;
232 :
233 : /**
234 : * The relative weight to give the dilation metric. The distortion metric is given weight 1 - _dilation_weight.
235 : */
236 : Real _dilation_weight;
237 :
238 : /* Map to hold target qp-dependent element target-to-reference mapping
239 : * Jacobians, if any
240 : */
241 : std::map<ElemType, std::vector<RealTensor>> _target_jacobians;
242 :
243 : /*
244 : * Map to hold the determinants of _target_jacobians.
245 : */
246 : std::map<ElemType, std::vector<Real>> _target_jacobian_dets;
247 :
248 : /**
249 : * Information about the mesh quality.
250 : */
251 : MeshQualityInfo _mesh_info;
252 :
253 : /**
254 : * Flag to indicate if the current solve is to untangle or smooth
255 : */
256 : bool _untangling_solve;
257 :
258 : /**
259 : * The quadrature rule type used to integrate the distortion-dilation metric.
260 : * See set_quadrature_type().
261 : */
262 : QuadratureType _quadrature_type;
263 : };
264 :
265 : } // namespace libMesh
266 :
267 : #endif // LIBMESH_VARIATIONAL_SMOOTHER_SYSTEM_H
|
