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 target element for a given element type.
145 : * @param type Element type
146 : * @return a std::pair containing the target element for type and the
147 : * corresponding nodes that must be kept in scope while the target element is
148 : * used.
149 : */
150 : static std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
151 : get_target_elem(const ElemType & type);
152 :
153 : /**
154 : * Get the jacobians (and determinants) of the target-to-reference element mapping.
155 : * @param target_elem Target element.
156 : * @param femcontext Context used to build mapping.
157 : * @param jacobian Vector in which to store the jacobians for each quadrature point.
158 : * @param jacobian_dets Vector in which to store the determinant of the jacobians
159 : * for each quadrature point.
160 : */
161 : static void get_target_to_reference_jacobian(const Elem * const target_elem,
162 : const FEMContext & femcontext,
163 : std::vector<RealTensor> & jacobians,
164 : std::vector<Real> & jacobian_dets);
165 :
166 : /**
167 : * Getter for the _mesh_info attribute. If this attribute has not yet been
168 : * initialized, compute_mesh_quality_info is called to initialize it.
169 : */
170 : const MeshQualityInfo & get_mesh_info();
171 :
172 : /*
173 : * Computes information about the mesh quality and sets the _mesh_info attribute.
174 : */
175 : void compute_mesh_quality_info();
176 :
177 : /*
178 : * Sets the verbosity of the object.
179 : */
180 2556 : void set_verbosity(const unsigned int verbosity) { _verbosity = verbosity; }
181 :
182 : protected:
183 :
184 : // System initialization
185 : virtual void init_data () override;
186 :
187 : // Context initialization
188 : virtual void init_context (libMesh::DiffContext & context) override;
189 :
190 : // Element residual and jacobian calculations
191 : // Time dependent parts
192 : virtual bool element_time_derivative (bool request_jacobian,
193 : libMesh::DiffContext & context) override;
194 :
195 : /* Computes the element reference volume used in the dilation metric
196 : * The reference value is set to the averaged value of all elements' average
197 : * |J|. Also computes any applicable target element inverse Jacobians. Target
198 : * elements are relavant when the reference element does not minimize the
199 : * distortion metric.
200 : */
201 : void prepare_for_smoothing();
202 :
203 : /**
204 : * Verbosity setting.
205 : * The verbosity levels and the corresponding information output are as
206 : * follows:
207 : *
208 : * verbosity = 0 : No information.
209 : *
210 : * 0 < verbosity : Prints:
211 : * - Initial mesh quality information
212 : * - The untangled mesh quality (if applicable)
213 : * - The smoothed mesh quality
214 : *
215 : * 10 < verbosity: Prints:
216 : * - The reference volume used for the dilation metric
217 : *
218 : * 50 < verbosity: Prints:
219 : * - Mesh quality information after each assembly
220 : *
221 : * 90 < verbosity: Prints:
222 : * - Quality information about each element after each assembly
223 : *
224 : */
225 : unsigned int _verbosity;
226 :
227 : /**
228 : * The small nonzero constant to prevent zero denominators (degenerate meshes only)
229 : */
230 : const Real _epsilon_squared;
231 :
232 : /**
233 : * Epsilon squared value determined at runtime during each assembly. The value
234 : * depends on whether the mesh is tangled.
235 : */
236 : Real _epsilon_squared_assembly;
237 :
238 : /**
239 : * The reference volume for each element
240 : */
241 : Real _ref_vol;
242 :
243 : /**
244 : * The relative weight to give the dilation metric. The distortion metric is given weight 1 - _dilation_weight.
245 : */
246 : Real _dilation_weight;
247 :
248 : /* Map to hold target qp-dependent element target-to-reference mapping
249 : * Jacobians, if any
250 : */
251 : std::map<ElemType, std::vector<RealTensor>> _target_jacobians;
252 :
253 : /*
254 : * Map to hold the determinants of _target_jacobians.
255 : */
256 : std::map<ElemType, std::vector<Real>> _target_jacobian_dets;
257 :
258 : /**
259 : * Information about the mesh quality.
260 : */
261 : MeshQualityInfo _mesh_info;
262 :
263 : /**
264 : * Flag to indicate if the current solve is to untangle or smooth
265 : */
266 : bool _untangling_solve;
267 :
268 : /**
269 : * The quadrature rule type used to integrate the distortion-dilation metric.
270 : * See set_quadrature_type().
271 : */
272 : QuadratureType _quadrature_type;
273 : };
274 :
275 : } // namespace libMesh
276 :
277 : #endif // LIBMESH_VARIATIONAL_SMOOTHER_SYSTEM_H
|
