Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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/fem_function_base.h"
24 : #include "libmesh/fem_system.h"
25 : #include "libmesh/libmesh_common.h"
26 :
27 : // C++ includes
28 : #include <map>
29 : #include <memory>
30 :
31 : namespace libMesh
32 : {
33 :
34 : /**
35 : * Struct to hold smoother-relevant information about the mesh quality.
36 : */
37 204 : struct MeshQualityInfo
38 : {
39 : // dof_id_types will hold the id of the Elem where the metric occurs
40 : // IMPORTANT: the Real should be the first entry of the pair so that taking
41 : // the min/max accross processors will compare the numeric values instead of
42 : // the element ids.
43 :
44 : std::pair<Real, dof_id_type> max_elem_distortion{std::numeric_limits<Real>::lowest(),
45 : DofObject::invalid_id};
46 : std::pair<Real, dof_id_type> min_elem_distortion{std::numeric_limits<Real>::max(),
47 : DofObject::invalid_id};
48 : Real total_distortion = 0.;
49 :
50 : std::pair<Real, dof_id_type> max_elem_dilation{std::numeric_limits<Real>::lowest(),
51 : DofObject::invalid_id};
52 : std::pair<Real, dof_id_type> min_elem_dilation{std::numeric_limits<Real>::max(),
53 : DofObject::invalid_id};
54 : Real total_dilation = 0.;
55 :
56 : std::pair<Real, dof_id_type> max_elem_combined{std::numeric_limits<Real>::lowest(),
57 : DofObject::invalid_id};
58 : std::pair<Real, dof_id_type> min_elem_combined{std::numeric_limits<Real>::max(),
59 : DofObject::invalid_id};
60 : Real total_combined = 0.;
61 :
62 : std::pair<Real, dof_id_type> max_elem_det_S{std::numeric_limits<Real>::lowest(),
63 : DofObject::invalid_id};
64 : std::pair<Real, dof_id_type> min_elem_det_S{std::numeric_limits<Real>::max(),
65 : DofObject::invalid_id};
66 : Real total_det_S = 0.;
67 : Real max_qp_det_S = std::numeric_limits<Real>::lowest();
68 : Real min_qp_det_S = std::numeric_limits<Real>::max();
69 :
70 : bool mesh_is_tangled = false;
71 : bool initialized = false;
72 : };
73 :
74 : // FEMSystem, TimeSolver and NewtonSolver will handle most tasks,
75 : // but we must specify element residuals
76 340 : class VariationalSmootherSystem : public libMesh::FEMSystem
77 : {
78 : /**
79 : * This is an FEMSystem to solve the optimization probelem posed by the
80 : * VariationalMeshSmoother class.
81 : *
82 : * The residual is coded as the gradient of the distortion-dilation metric, and
83 : * the jacobian as analytically coded as the Hessian of the metric.
84 : *
85 : * The nodes of the system mesh are updated during the solve.
86 : */
87 : public:
88 2414 : VariationalSmootherSystem(libMesh::EquationSystems & es,
89 : const std::string & name,
90 : const unsigned int number)
91 2414 : : libMesh::FEMSystem(es, name, number),
92 2278 : _verbosity(0),
93 2278 : _epsilon_squared(TOLERANCE),
94 2278 : _epsilon_squared_assembly(0.),
95 2278 : _ref_vol(0.),
96 2278 : _dilation_weight(0.5),
97 2482 : _untangling_solve(false)
98 2414 : {}
99 :
100 : // Default destructor
101 : ~VariationalSmootherSystem() override;
102 :
103 : /**
104 : * Assembly method to update the mesh based on the smoother solve.
105 : */
106 : virtual void assembly (bool get_residual,
107 : bool get_jacobian,
108 : bool apply_heterogeneous_constraints = false,
109 : bool apply_no_constraints = false) override;
110 :
111 68 : Real & get_dilation_weight() { return _dilation_weight; }
112 :
113 : /**
114 : * Solves the system to smooth the mesh. If the mesh is initially tangled,
115 : * a solve is first performed to untangle the mesh, followed by a solve to
116 : * smooth the mesh.
117 : */
118 : virtual void solve() override;
119 :
120 : /**
121 : * Get the target element for a given element type.
122 : * @param type Element type
123 : * @return a std::pair containing the target element for type and the
124 : * corresponding nodes that must be kept in scope while the target element is
125 : * used.
126 : */
127 : static std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
128 : get_target_elem(const ElemType & type);
129 :
130 : /**
131 : * Get the jacobians (and determinants) of the target-to-reference element mapping.
132 : * @param target_elem Target element.
133 : * @param femcontext Context used to build mapping.
134 : * @param jacobian Vector in which to store the jacobians for each quadrature point.
135 : * @param jacobian_dets Vector in which to store the determinant of the jacobians
136 : * for each quadrature point.
137 : */
138 : static void get_target_to_reference_jacobian(const Elem * const target_elem,
139 : const FEMContext & femcontext,
140 : std::vector<RealTensor> & jacobians,
141 : std::vector<Real> & jacobian_dets);
142 :
143 : /**
144 : * Getter for the _mesh_info attribute. If this attribute has not yet been
145 : * initialized, compute_mesh_quality_info is called to initialize it.
146 : */
147 : const MeshQualityInfo & get_mesh_info();
148 :
149 : /*
150 : * Computes information about the mesh quality and sets the _mesh_info attribute.
151 : */
152 : void compute_mesh_quality_info();
153 :
154 : /*
155 : * Sets the verbosity of the object.
156 : */
157 2414 : void set_verbosity(const unsigned int verbosity) { _verbosity = verbosity; }
158 :
159 : protected:
160 :
161 : // System initialization
162 : virtual void init_data () override;
163 :
164 : // Context initialization
165 : virtual void init_context (libMesh::DiffContext & context) override;
166 :
167 : // Element residual and jacobian calculations
168 : // Time dependent parts
169 : virtual bool element_time_derivative (bool request_jacobian,
170 : libMesh::DiffContext & context) override;
171 :
172 : /* Computes the element reference volume used in the dilation metric
173 : * The reference value is set to the averaged value of all elements' average
174 : * |J|. Also computes any applicable target element inverse Jacobians. Target
175 : * elements are relavant when the reference element does not minimize the
176 : * distortion metric.
177 : */
178 : void prepare_for_smoothing();
179 :
180 : /**
181 : * Verbosity setting.
182 : * The verbosity levels and the corresponding information output are as
183 : * follows:
184 : *
185 : * verbosity = 0 : No information.
186 : *
187 : * 0 < verbosity : Prints:
188 : * - Initial mesh quality information
189 : * - The untangled mesh quality (if applicable)
190 : * - The smoothed mesh quality
191 : *
192 : * 10 < verbosity: Prints:
193 : * - The reference volume used for the dilation metric
194 : *
195 : * 50 < verbosity: Prints:
196 : * - Mesh quality information after each assembly
197 : *
198 : * 90 < verbosity: Prints:
199 : * - Quality information about each element after each assembly
200 : *
201 : */
202 : unsigned int _verbosity;
203 :
204 : /**
205 : * The small nonzero constant to prevent zero denominators (degenerate meshes only)
206 : */
207 : const Real _epsilon_squared;
208 :
209 : /**
210 : * Epsilon squared value determined at runtime during each assembly. The value
211 : * depends on whether the mesh is tangled.
212 : */
213 : Real _epsilon_squared_assembly;
214 :
215 : /**
216 : * The reference volume for each element
217 : */
218 : Real _ref_vol;
219 :
220 : /**
221 : * The relative weight to give the dilation metric. The distortion metric is given weight 1 - _dilation_weight.
222 : */
223 : Real _dilation_weight;
224 :
225 : /* Map to hold target qp-dependent element target-to-reference mapping
226 : * Jacobians, if any
227 : */
228 : std::map<ElemType, std::vector<RealTensor>> _target_jacobians;
229 :
230 : /*
231 : * Map to hold the determinants of _target_jacobians.
232 : */
233 : std::map<ElemType, std::vector<Real>> _target_jacobian_dets;
234 :
235 : /**
236 : * Information about the mesh quality.
237 : */
238 : MeshQualityInfo _mesh_info;
239 :
240 : /**
241 : * Flag to indicate if the current solve is to untangle or smooth
242 : */
243 : bool _untangling_solve;
244 : };
245 :
246 : } // namespace libMesh
247 :
248 : #endif // LIBMESH_VARIATIONAL_SMOOTHER_SYSTEM_H
|
