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