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 114 : 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 190 : 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 1349 : VariationalSmootherSystem(libMesh::EquationSystems & es,
89 : const std::string & name,
90 : const unsigned int number)
91 1349 : : libMesh::FEMSystem(es, name, number),
92 1273 : _epsilon_squared(TOLERANCE),
93 1273 : _epsilon_squared_assembly(0.),
94 1273 : _ref_vol(0.),
95 1273 : _dilation_weight(0.5),
96 1387 : _untangling_solve(false)
97 1349 : {}
98 :
99 : // Default destructor
100 : ~VariationalSmootherSystem() override;
101 :
102 : /**
103 : * Assembly method to update the mesh based on the smoother solve.
104 : */
105 : virtual void assembly (bool get_residual,
106 : bool get_jacobian,
107 : bool apply_heterogeneous_constraints = false,
108 : bool apply_no_constraints = false) override;
109 :
110 38 : Real & get_dilation_weight() { return _dilation_weight; }
111 :
112 : /**
113 : * Solves the system to smooth the mesh. If the mesh is initially tangled,
114 : * a solve is first performed to untangle the mesh, followed by a solve to
115 : * smooth the mesh.
116 : */
117 : virtual void solve() override;
118 :
119 : /**
120 : * Get the target element for a given element type.
121 : * @param type Element type
122 : * @return a std::pair containing the target element for type and the
123 : * corresponding nodes that must be kept in scope while the target element is
124 : * used.
125 : */
126 : static std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
127 : get_target_elem(const ElemType & type);
128 :
129 : /**
130 : * Get the jacobians (and determinants) of the target-to-reference element mapping.
131 : * @param target_elem Target element.
132 : * @param femcontext Context used to build mapping.
133 : * @param jacobian Vector in which to store the jacobians for each quadrature point.
134 : * @param jacobian_dets Vector in which to store the determinant of the jacobians
135 : * for each quadrature point.
136 : */
137 : static void get_target_to_reference_jacobian(const Elem * const target_elem,
138 : const FEMContext & femcontext,
139 : std::vector<RealTensor> & jacobians,
140 : std::vector<Real> & jacobian_dets);
141 :
142 : /**
143 : * Getter for the _mesh_info attribute. If this attribute has not yet been
144 : * initialized, compute_mesh_quality_info is called to initialize it.
145 : */
146 : const MeshQualityInfo & get_mesh_info();
147 :
148 : /*
149 : * Computes information about the mesh quality and sets the _mesh_info attribute.
150 : */
151 : void compute_mesh_quality_info();
152 :
153 : protected:
154 :
155 : // System initialization
156 : virtual void init_data () override;
157 :
158 : // Context initialization
159 : virtual void init_context (libMesh::DiffContext & context) override;
160 :
161 : // Element residual and jacobian calculations
162 : // Time dependent parts
163 : virtual bool element_time_derivative (bool request_jacobian,
164 : libMesh::DiffContext & context) override;
165 :
166 : /* Computes the element reference volume used in the dilation metric
167 : * The reference value is set to the averaged value of all elements' average
168 : * |J|. Also computes any applicable target element inverse Jacobians. Target
169 : * elements are relavant when the reference element does not minimize the
170 : * distortion metric.
171 : */
172 : void prepare_for_smoothing();
173 :
174 : /**
175 : * The small nonzero constant to prevent zero denominators (degenerate meshes only)
176 : */
177 : const Real _epsilon_squared;
178 :
179 : /**
180 : * Epsilon squared value determined at runtime during each assembly. The value
181 : * depends on whether the mesh is tangled.
182 : */
183 : Real _epsilon_squared_assembly;
184 :
185 : /**
186 : * The reference volume for each element
187 : */
188 : Real _ref_vol;
189 :
190 : /**
191 : * The relative weight to give the dilation metric. The distortion metric is given weight 1 - _dilation_weight.
192 : */
193 : Real _dilation_weight;
194 :
195 : /* Map to hold target qp-dependent element target-to-reference mapping
196 : * Jacobians, if any
197 : */
198 : std::map<ElemType, std::vector<RealTensor>> _target_jacobians;
199 :
200 : /*
201 : * Map to hold the determinants of _target_jacobians.
202 : */
203 : std::map<ElemType, std::vector<Real>> _target_jacobian_dets;
204 :
205 : /**
206 : * Information about the mesh quality.
207 : */
208 : MeshQualityInfo _mesh_info;
209 :
210 : /**
211 : * Flag to indicate if the current solve is to untangle or smooth
212 : */
213 : bool _untangling_solve;
214 : };
215 :
216 : } // namespace libMesh
217 :
218 : #endif // LIBMESH_VARIATIONAL_SMOOTHER_SYSTEM_H
|
