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 : #include "libmesh/variational_smoother_system.h"
19 :
20 : #include "libmesh/elem.h"
21 : #include "libmesh/face_tri3.h"
22 : #include "libmesh/face_tri6.h"
23 : #include "libmesh/fe_base.h"
24 : #include "libmesh/fe_interface.h"
25 : #include "libmesh/fem_context.h"
26 : #include "libmesh/mesh.h"
27 : #include "libmesh/node_range.h"
28 : #include "libmesh/numeric_vector.h"
29 : #include "libmesh/parallel_ghost_sync.h"
30 : #include "libmesh/quadrature.h"
31 : #include "libmesh/string_to_enum.h"
32 : #include "libmesh/utility.h"
33 : #include "libmesh/enum_to_string.h"
34 : #include <libmesh/reference_elem.h>
35 :
36 : // C++ includes
37 : #include <functional> // std::reference_wrapper
38 :
39 : namespace libMesh
40 : {
41 :
42 : /*
43 : * Gets the dof_id_type value corresponding to the minimum of the Real value.
44 : */
45 140008 : void communicate_pair_min(std::pair<Real, dof_id_type> & pair, const Parallel::Communicator & comm)
46 : {
47 : // Get rank where minimum occurs
48 : unsigned int rank;
49 140008 : comm.minloc(pair.first, rank);
50 140008 : comm.broadcast(pair.second, rank);
51 140008 : }
52 :
53 : /*
54 : * Gets the dof_id_type value corresponding to the maximum of the Real value.
55 : */
56 140008 : void communicate_pair_max(std::pair<Real, dof_id_type> & pair, const Parallel::Communicator & comm)
57 : {
58 : // Get rank where minimum occurs
59 : unsigned int rank;
60 140008 : comm.maxloc(pair.first, rank);
61 140008 : comm.broadcast(pair.second, rank);
62 140008 : }
63 :
64 : /**
65 : * Function to prevent dividing by zero for degenerate elements
66 : */
67 6802624 : Real chi_epsilon(const Real & x, const Real epsilon_squared)
68 : {
69 7369296 : return 0.5 * (x + std::sqrt(epsilon_squared + Utility::pow<2>(x)));
70 : }
71 :
72 : /**
73 : * Given an fe_map, element dimension, and quadrature point index, returns the
74 : * Jacobian of the physical-to-reference mapping.
75 : */
76 6829137 : RealTensor get_jacobian_at_qp(const FEMap & fe_map,
77 : const unsigned int & dim,
78 : const unsigned int & qp)
79 : {
80 6829137 : libmesh_error_msg_if(dim > 3, "Unsupported dimension.");
81 :
82 : // RealTensors are always 3x3, so we will fill any dimensions above dim
83 : // with 1s on the diagonal. This indicates a 1 to 1 relationship between
84 : // the physical and reference elements in these extra dimensions.
85 :
86 6829137 : const auto & dxyzdxi = dim >= 1 ? fe_map.get_dxyzdxi()[qp] : RealGradient(1, 0, 0);
87 6829137 : const auto & dxyzdeta = dim >= 2 ? fe_map.get_dxyzdeta()[qp] : RealGradient(0, 1, 0);
88 6829137 : const auto & dxyzdzeta = dim >= 3 ? fe_map.get_dxyzdzeta()[qp] : RealGradient(0, 0, 1);
89 :
90 6829137 : return RealTensor(dxyzdxi, dxyzdeta, dxyzdzeta).transpose(); // Note the transposition!
91 : }
92 :
93 : /**
94 : * Compute the trace of a dim-dimensional matrix.
95 : */
96 6802624 : Real trace(const RealTensor & A, const unsigned int & dim)
97 : {
98 566672 : Real tr = 0.0;
99 26823120 : for (const auto i : make_range(dim))
100 20020496 : tr += A(i, i);
101 :
102 6802624 : return tr;
103 : }
104 :
105 7236 : VariationalSmootherSystem::~VariationalSmootherSystem () = default;
106 :
107 32446 : void VariationalSmootherSystem::assembly (bool get_residual,
108 : bool get_jacobian,
109 : bool apply_heterogeneous_constraints,
110 : bool apply_no_constraints)
111 : {
112 : // Update the mesh based on the current variable values
113 1816 : auto & mesh = this->get_mesh();
114 32446 : this->solution->close();
115 :
116 32446 : const auto dim = mesh.mesh_dimension();
117 :
118 : Threads::parallel_for_each
119 32446 : (mesh.local_node_stored_range(),
120 1838412 : [this, dim](Node * node)
121 : {
122 1838412 : for (const auto d : make_range(dim))
123 : {
124 1362664 : const auto dof_id = node->dof_number(this->number(), d, 0);
125 : // Update mesh
126 1362664 : (*node)(d) = libmesh_real((*current_local_solution)(dof_id));
127 : }
128 475748 : });
129 :
130 32446 : SyncNodalPositions sync_object(mesh);
131 63984 : Parallel::sync_dofobject_data_by_id (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);
132 :
133 : // Compute and update mesh quality information
134 32446 : compute_mesh_quality_info();
135 32446 : const bool is_tangled = _mesh_info.mesh_is_tangled;
136 :
137 : // Update _epsilon_squared_assembly based on whether we are untangling or
138 : // smoothing
139 32446 : if (_untangling_solve)
140 : {
141 172 : const Real & min_S = _mesh_info.min_qp_det_S;
142 : // This component is based on what Larisa did in the original code.
143 6318 : const Real variable_component = 100. * Utility::pow<2>(_ref_vol * min_S);
144 6318 : _epsilon_squared_assembly = _epsilon_squared + variable_component;
145 : }
146 :
147 : else
148 26128 : _epsilon_squared_assembly = 0.;
149 :
150 32446 : FEMSystem::assembly(get_residual, get_jacobian, apply_heterogeneous_constraints, apply_no_constraints);
151 :
152 32446 : if (_untangling_solve && !is_tangled)
153 : {
154 : // Mesh is untangled, artificially reduce residual by factor 0.9 to tell
155 : // the solver we are on the right track. The mesh may become re-tangled in
156 : // subsequent nonlinear and line search iterations, so we will not use
157 : // this reduction factor in the case of re-tangulation. This approach
158 : // should drive the solver to eventually favor untangled solutions.
159 4472 : rhs->close();
160 4472 : (*rhs) *= 0.9;
161 : }
162 32446 : }
163 :
164 2556 : void VariationalSmootherSystem::solve()
165 : {
166 2556 : const auto & mesh_info = get_mesh_info();
167 2556 : if (_verbosity > 10)
168 4 : libMesh::out << "Initial " << mesh_info << std::endl;
169 2556 : if (mesh_info.mesh_is_tangled)
170 : {
171 : // Untangling solve
172 142 : _untangling_solve = true;
173 :
174 : // Untangling seems to work better using only the distortion metric.
175 : // For a mixed metric, I've seen it *technically* untangle a mesh to a
176 : // still-suboptimal mesh (i.e., a local minima), but not been able to
177 : // smooth it well because the smoothest solution is on the other side of
178 : // a tangulation barrier.
179 142 : const auto dilation_weight = _dilation_weight;
180 142 : _dilation_weight = 0.;
181 :
182 142 : if (_verbosity > 10)
183 0 : libMesh::out << "Untangling the mesh" << std::endl;
184 142 : FEMSystem::solve();
185 :
186 : // Reset the dilation weight
187 142 : _dilation_weight = dilation_weight;
188 :
189 142 : if (_verbosity > 10)
190 0 : libMesh::out << "Untangled " << mesh_info << std::endl;
191 : }
192 :
193 : // Smoothing solve
194 2556 : _untangling_solve = false;
195 2556 : if (_verbosity > 10)
196 4 : libMesh::out << "Smoothing the mesh" << std::endl;
197 2556 : FEMSystem::solve();
198 2556 : libmesh_error_msg_if(mesh_info.mesh_is_tangled, "The smoothing solve tangled the mesh!");
199 2556 : if (_verbosity > 10)
200 4 : libMesh::out << "Smoothed " << mesh_info << std::endl;
201 2556 : }
202 :
203 2556 : void VariationalSmootherSystem::init_data ()
204 : {
205 144 : auto & mesh = this->get_mesh();
206 72 : const auto & elem_orders = mesh.elem_default_orders();
207 2556 : libmesh_error_msg_if(elem_orders.size() != 1,
208 : "The variational smoother cannot be used for mixed-order meshes!");
209 2556 : const auto fe_order = *elem_orders.begin();
210 : // Add a variable for each dimension of the mesh
211 : // "r0" for x, "r1" for y, "r2" for z
212 8946 : for (const auto & d : make_range(mesh.mesh_dimension()))
213 12600 : this->add_variable ("r" + std::to_string(d), fe_order);
214 :
215 : // Do the parent's initialization after variables are defined
216 2556 : FEMSystem::init_data();
217 :
218 : // Set the current_local_solution to the current mesh for the initial guess
219 2556 : this->solution->close();
220 :
221 77772 : for (auto * node : mesh.local_node_ptr_range())
222 : {
223 154332 : for (const auto d : make_range(mesh.mesh_dimension()))
224 : {
225 114660 : const auto dof_id = node->dof_number(this->number(), d, 0);
226 : // Update solution
227 114660 : (*solution).set(dof_id, (*node)(d));
228 : }
229 2412 : }
230 :
231 2556 : this->prepare_for_smoothing();
232 2556 : }
233 :
234 2556 : void VariationalSmootherSystem::prepare_for_smoothing()
235 : {
236 : // If this method has already been called to set _ref_vol, just return.
237 2628 : if (std::abs(_ref_vol) > TOLERANCE * TOLERANCE)
238 0 : return;
239 :
240 2628 : std::unique_ptr<DiffContext> con = this->build_context();
241 72 : FEMContext & femcontext = cast_ref<FEMContext &>(*con);
242 2556 : this->init_context(femcontext);
243 :
244 144 : const auto & mesh = this->get_mesh();
245 :
246 2556 : Real elem_averaged_det_S_sum = 0.;
247 :
248 : // Make pre-requests before reinit() for efficiency in
249 : // --enable-deprecated builds, and to avoid errors in
250 : // --disable-deprecated builds.
251 72 : const auto & fe_map = femcontext.get_element_fe(0)->get_fe_map();
252 72 : const auto & JxW = fe_map.get_JxW();
253 :
254 40968 : for (const auto * elem : mesh.active_local_element_ptr_range())
255 : {
256 35928 : femcontext.pre_fe_reinit(*this, elem);
257 35928 : femcontext.elem_fe_reinit();
258 :
259 : // Add target element info, if applicable
260 38922 : if (_target_jacobians.find(elem->type()) == _target_jacobians.end())
261 : {
262 2094 : const auto [target_elem, target_nodes] = ReferenceElem::ideal_target(elem->type());
263 2094 : get_target_to_reference_jacobian(target_elem.get(),
264 : femcontext,
265 2094 : _target_jacobians[elem->type()],
266 4188 : _target_jacobian_dets[elem->type()]);
267 : }// if find == end()
268 :
269 : // Reference volume computation
270 2994 : Real elem_integrated_det_S = 0.;
271 337464 : for (const auto qp : index_range(JxW))
272 326664 : elem_integrated_det_S += JxW[qp] / _target_jacobian_dets[elem->type()][qp];
273 35928 : const auto ref_elem_vol = elem->reference_elem()->volume();
274 35928 : elem_averaged_det_S_sum += elem_integrated_det_S / ref_elem_vol;
275 :
276 2412 : } // for elem
277 :
278 : // Get contributions from elements on other processors
279 2556 : mesh.comm().sum(elem_averaged_det_S_sum);
280 :
281 2556 : _ref_vol = elem_averaged_det_S_sum / mesh.n_active_elem();
282 2556 : if (_verbosity > 10)
283 8 : libMesh::out << "Reference volume: " << _ref_vol << std::endl;
284 2412 : }
285 :
286 70004 : void VariationalSmootherSystem::init_context(DiffContext & context)
287 : {
288 1960 : FEMContext & c = cast_ref<FEMContext &>(context);
289 :
290 1960 : FEBase * my_fe = nullptr;
291 :
292 : // Now make sure we have requested all the data
293 : // we need to build the system.
294 :
295 : // We might have a multi-dimensional mesh
296 : const std::set<unsigned char> & elem_dims =
297 1960 : c.elem_dimensions();
298 :
299 140008 : for (const auto & dim : elem_dims)
300 : {
301 70004 : c.get_element_fe( 0, my_fe, dim );
302 1960 : my_fe->get_nothing();
303 :
304 1960 : auto & fe_map = my_fe->get_fe_map();
305 1960 : fe_map.get_dxyzdxi();
306 1960 : fe_map.get_dxyzdeta();
307 1960 : fe_map.get_dxyzdzeta();
308 1960 : fe_map.get_JxW();
309 :
310 : // Mesh may be tangled, allow negative Jacobians
311 1960 : fe_map.set_jacobian_tolerance(std::numeric_limits<Real>::lowest());
312 :
313 70004 : c.get_side_fe( 0, my_fe, dim );
314 1960 : my_fe->get_nothing();
315 : }
316 :
317 70004 : FEMSystem::init_context(context);
318 :
319 : // Override the default (interior-only) Gauss rule if a different quadrature
320 : // type has been requested. A vertex-sampling rule (e.g. QTRAP, QSIMPSON,
321 : // QNODAL, QGAUSS_LOBATTO) lets the metric "see" degeneracies localized at
322 : // element corners that interior Gauss points miss. Do this after
323 : // FEMSystem::init_context() so it is not overwritten. The System's
324 : // extra_quadrature_order is respected for the point count.
325 70004 : if (_quadrature_type != QGAUSS)
326 1562 : c.use_quadrature_rules(_quadrature_type, this->extra_quadrature_order);
327 70004 : }
328 :
329 :
330 383672 : bool VariationalSmootherSystem::element_time_derivative (bool request_jacobian,
331 : DiffContext & context)
332 : {
333 31946 : FEMContext & c = cast_ref<FEMContext &>(context);
334 :
335 63892 : const Elem & elem = c.get_elem();
336 :
337 383672 : unsigned int dim = c.get_dim();
338 :
339 31946 : unsigned int x_var = 0, y_var = 1, z_var = 2;
340 : // In the case of lower dimensions, we will not access z (1D/2D) or y (1D)
341 383672 : if (dim < 3)
342 : {
343 3672 : z_var = 0;
344 44384 : if (dim < 2)
345 100 : y_var = 0;
346 : }
347 :
348 : // The subvectors and submatrices we need to fill:
349 : // system residual
350 : std::reference_wrapper<DenseSubVector<Number>> F[3] =
351 : {
352 : c.get_elem_residual(x_var),
353 : c.get_elem_residual(y_var),
354 : c.get_elem_residual(z_var)
355 31946 : };
356 : // system jacobian
357 : std::reference_wrapper<DenseSubMatrix<Number>> K[3][3] =
358 : {
359 : {c.get_elem_jacobian(x_var, x_var), c.get_elem_jacobian(x_var, y_var), c.get_elem_jacobian(x_var, z_var)},
360 : {c.get_elem_jacobian(y_var, x_var), c.get_elem_jacobian(y_var, y_var), c.get_elem_jacobian(y_var, z_var)},
361 : {c.get_elem_jacobian(z_var, x_var), c.get_elem_jacobian(z_var, y_var), c.get_elem_jacobian(z_var, z_var)}
362 31946 : };
363 :
364 : // Quadrature info
365 31946 : const auto & quad_weights = c.get_element_qrule().get_weights();
366 :
367 383672 : const auto distortion_weight = 1. - _dilation_weight;
368 :
369 : // Get some references to cell-specific data that
370 : // will be used to assemble the linear system.
371 :
372 31946 : const auto & fe_map = c.get_element_fe(0)->get_fe_map();
373 :
374 31946 : const auto & dphidxi_map = fe_map.get_dphidxi_map();
375 31946 : const auto & dphideta_map = fe_map.get_dphideta_map();
376 31946 : const auto & dphidzeta_map = fe_map.get_dphidzeta_map();
377 :
378 383672 : const auto & target_jacobian_dets = _target_jacobian_dets[elem.type()];
379 383672 : const auto & target_jacobians = _target_jacobians[elem.type()];
380 :
381 : // Integrate the distortion-dilation metric over the reference element
382 3634216 : for (const auto qp : index_range(quad_weights))
383 : {
384 : // Compute quantities needed to evaluate the distortion-dilation metric
385 : // and its gradient and Hessian.
386 : // The metric will be minimized when it's gradient with respect to the node
387 : // locations (R) is zero. For Newton's method, minimizing the gradient
388 : // requires computation of the gradient's Jacobian with respect to R.
389 : // The Jacobian of the distortion-dilation metric's gradientis the Hessian
390 : // of the metric.
391 :
392 : // Transform quad weight from reference element to quad weight for target element
393 3521316 : const auto quad_weight = quad_weights[qp] / target_jacobian_dets[qp];
394 :
395 : // Note that the term "Jacobian" has two meanings in this mesh smoothing
396 : // application. The first meaning refers to the Jacobian w.r.t R of the
397 : // gradient w.r.t. R, (i.e., the Hessian of the metric). This is the K
398 : // variable defined above. This is also the Jacobian the 'request_jacobian'
399 : // variable refers to. The second meaning refers to the Jacobian of the
400 : // physical-to-target element mapping. This is the Jacobian used to
401 : // compute the distorion-dilation metric.
402 : //
403 : // Grab the physical-to-reference mapping Jacobian matrix (i.e., "S") at this qp
404 3521316 : RealTensor S = get_jacobian_at_qp(fe_map, dim, qp);
405 :
406 : // Apply target element transformation to get the physical-to-target jacobian
407 3250544 : S *= _target_jacobians[elem.type()][qp];
408 :
409 : // Compute quantities needed for the smoothing algorithm
410 :
411 : // determinant
412 3250544 : const Real det = S.det();
413 3250544 : const Real det_sq = det * det;
414 3250544 : const Real det_cube = det_sq * det;
415 :
416 3250544 : const Real ref_vol_sq = _ref_vol * _ref_vol;
417 :
418 : // trace of S^T * S
419 : // DO NOT USE RealTensor.tr for the trace, it will NOT be correct for
420 : // 1D and 2D meshes because of our hack of putting 1s in the diagonal of
421 : // S for the extra dimensions (see get_jacobian_at_qp)
422 3250544 : const auto tr = trace(S.transpose() * S, dim);
423 3250544 : const Real tr_div_dim = tr / dim;
424 :
425 : // Precompute pow(tr_div_dim, 0.5 * dim - x) for x = 0, 1, 2
426 3250544 : const Real half_dim = 0.5 * dim;
427 : const std::vector<Real> trace_powers{
428 3250544 : std::pow(tr_div_dim, half_dim),
429 3250544 : std::pow(tr_div_dim, half_dim - 1.),
430 3250544 : std::pow(tr_div_dim, half_dim - 2.),
431 3521316 : };
432 :
433 : // inverse of S
434 3521316 : const RealTensor S_inv = S.inverse();
435 : // inverse transpose of S
436 541544 : const RealTensor S_inv_T = S_inv.transpose();
437 :
438 : // Identity matrix
439 5418000 : const RealTensor I(1, 0, 0,
440 5418000 : 0, 1, 0,
441 3250544 : 0, 0, 1);
442 :
443 : // The chi function allows us to handle degenerate elements
444 3250544 : const auto chi = chi_epsilon(det, _epsilon_squared_assembly);
445 3250544 : const Real chi_sq = chi * chi;
446 3250544 : const Real sqrt_term = std::sqrt(_epsilon_squared_assembly + det_sq);
447 : // dchi(x) / dx
448 3250544 : const Real chi_prime = 0.5 * (1. + det / sqrt_term);
449 3250544 : const Real chi_prime_sq = chi_prime * chi_prime;
450 : // d2chi(x) / dx2
451 3250544 : const Real chi_2prime = 0.5 * (1. / sqrt_term - det_sq / Utility::pow<3>(sqrt_term));
452 :
453 : // Distortion metric (beta)
454 : //const Real beta = trace_powers[0] / chi;
455 4333632 : const RealTensor dbeta_dS = (trace_powers[1] / chi) * S - (trace_powers[0] / chi_sq * chi_prime * det) * S_inv_T;
456 :
457 : // Dilation metric (mu)
458 : //const Real mu = 0.5 * (_ref_vol + det_sq / _ref_vol) / chi;
459 : // We represent d mu / dS as alpha(S) * S^-T, where alpha is a scalar function
460 3250544 : const Real alpha = (-chi_prime * det_cube + 2. * det_sq * chi - ref_vol_sq * det * chi_prime) / (2. * _ref_vol * chi_sq);
461 3250544 : const RealTensor dmu_dS = alpha * S_inv_T;
462 :
463 : // Combined metric (E)
464 : //const Real E = distortion_weight * beta + _dilation_weight * mu;
465 4062860 : const RealTensor dE_dS = distortion_weight * dbeta_dS + _dilation_weight * dmu_dS;
466 :
467 : // This vector is useful in computing dS/dR below
468 14085264 : std::vector<std::vector<std::vector<Real>>> dphi_maps = {dphidxi_map, dphideta_map, dphidzeta_map};
469 :
470 : // Compute residual (i.e., the gradient of the combined metric w.r.t node locations)
471 : // Recall that when the gradient (residual) is zero, the combined metric is minimized
472 32701888 : for (const auto l : elem.node_index_range())
473 : {
474 117065664 : for (const auto var_id : make_range(dim))
475 : {
476 : // Build dS/dR, the derivative of the physical-to-target mapping Jacobin w.r.t. the mesh
477 : // node locations
478 80313980 : RealTensor dS_dR = RealTensor(0);
479 348993696 : for (const auto jj : make_range(dim))
480 348499008 : dS_dR(var_id, jj) = dphi_maps[jj][l][qp];
481 94914660 : dS_dR *= target_jacobians[qp];
482 :
483 : // Residual contribution. The contraction of dE/dS and dS/dR gives us
484 : // the gradient we are looking for, dE/dR
485 94914660 : F[var_id](l) += quad_weight * dE_dS.contract(dS_dR);
486 : }// for var_id
487 : }// for l
488 :
489 :
490 3250544 : if (request_jacobian)
491 : {
492 : // Compute jacobian of the smoothing system (i.e., the Hessian of the
493 : // combined metric w.r.t. the mesh node locations)
494 :
495 : // Precompute coefficients to be applied to each component of tensor
496 : // products in the loops below. At first glance, these coefficients look
497 : // like gibberish, but everything in the Hessian has been verified by
498 : // taking finite differences of the gradient. We should probably write
499 : // down the derivations of the gradient and Hessian somewhere...
500 :
501 : // Recall that above, dbeta_dS takes the form:
502 : // d(beta)/dS = c1(S) * S - c2(S) * S_inv_T,
503 : // where c1 and c2 are scalar-valued functions.
504 : const std::vector<Real> d2beta_dS2_coefs_times_distortion_weight = {
505 : //Part 1: scaler coefficients of d(c1 * S) / dS
506 : //
507 : // multiplies I[i,a] x I[j,b]
508 1599752 : (trace_powers[1] / chi) * distortion_weight,
509 : // multiplies S[a,b] x S[i,j]
510 1733054 : (((dim - 2.) / dim) * trace_powers[2] / chi) * distortion_weight,
511 : // multiplies S_inv[b,a] * S[i,j]
512 1733054 : (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
513 : //
514 : //Part 2: scaler coefficients of d(-c2 * S_inv_T) / dS
515 : //
516 : // multiplies S[a,b] x S_inv[j,i]
517 1733054 : (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
518 : // multiplies S_inv[b,a] x S_inv[j,i]
519 1733054 : (trace_powers[0] * (det / chi_sq)
520 1599752 : * ((2. * chi_prime_sq / chi - chi_2prime) * det - chi_prime)) * distortion_weight,
521 : // multiplies S_inv[b,i] x S_inv[j,a]
522 1733054 : ((trace_powers[0] / chi_sq) * chi_prime * det) * distortion_weight,
523 1733054 : };
524 :
525 : // d alpha / dS has the form c(S) * S^-T, where c is the scalar coefficient defined below
526 1866356 : const Real dalpha_dS_coef_times_dilation_weight = ((det / (2. * _ref_vol * chi_sq))
527 1866356 : * (-4. * _ref_vol * alpha * chi * chi_prime - chi_2prime * det_cube
528 1866356 : - chi_prime * det_sq + 4 * chi * det - ref_vol_sq * (chi_prime + det * chi_2prime))) * _dilation_weight;
529 :
530 : // This is also useful to precompute
531 1599752 : const Real alpha_times_dilation_weight = alpha * _dilation_weight;
532 :
533 : /*
534 :
535 : To increase the efficiency of the Jacobian computation, we take
536 : advantage of l-p symmetry, ij-ab symmetry, and the sparsity pattern of
537 : dS_dR. We also factor all possible multipliers out of inner loops for
538 : efficiency. The result is code that is more difficult to read. For
539 : clarity, consult the pseudo-code below in this comment.
540 :
541 : for (const auto l: elem.node_index_range()) // Contribution to Hessian
542 : from node l
543 : {
544 : for (const auto var_id1 : make_range(dim)) // Contribution from each
545 : x/y/z component of node l
546 : {
547 : // Build dS/dR_l, the derivative of the physical-to-target
548 : mapping Jacobin w.r.t.
549 : // the l-th node
550 : RealTensor dS_dR_l = RealTensor(0);
551 : for (const auto ii : make_range(dim))
552 : dS_dR_l(var_id1, ii) = dphi_maps[ii][l][qp];
553 : dS_dR_l *= target_jacobians[qp];
554 :
555 : for (const auto p: elem.node_index_range()) // Contribution to
556 : Hessian from node p
557 : {
558 : for (const auto var_id2 : make_range(dim)) // Contribution from
559 : each x/y/z component of node p
560 : {
561 : // Build dS/dR_l, the derivative of the physical-to-target
562 : mapping Jacobin w.r.t.
563 : // the p-th node
564 : RealTensor dS_dR_p = RealTensor(0);
565 : for (const auto jj : make_range(dim))
566 : dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
567 : dS_dR_p *= target_jacobians[qp];
568 :
569 : Real d2beta_dR2 = 0.;
570 : Real d2mu_dR2 = 0.;
571 : // Perform tensor contraction
572 : for (const auto i : make_range(dim))
573 : {
574 : for (const auto j : make_range(dim))
575 : {
576 : for (const auto a : make_range(dim))
577 : {
578 : for (const auto b : make_range(dim))
579 : {
580 : // Nasty tensor products to be multiplied by
581 : d2beta_dS2_coefs to get d2(beta) / dS2 const std::vector<Real>
582 : d2beta_dS2_tensor_contributions =
583 : {
584 : I(i,a) * I(j,b),
585 : S(a,b) * S(i,j),
586 : S_inv(b,a) * S(i,j),
587 : S(a,b) * S_inv(j,i),
588 : S_inv(b,a) * S_inv(j,i),
589 : S_inv(j,a) * S_inv(b,i),
590 : };
591 :
592 : // Combine precomputed coefficients with tensor products
593 : to get d2(beta) / dS2 Real d2beta_dS2 = 0.; for (const auto comp_id :
594 : index_range(d2beta_dS2_coefs))
595 : {
596 : const Real contribution = d2beta_dS2_coefs[comp_id] *
597 : d2beta_dS2_tensor_contributions[comp_id]; d2beta_dS2 += contribution;
598 : }
599 :
600 : // Incorporate tensor product portion to get d2(mu) /
601 : dS2 const Real d2mu_dS2 = dalpha_dS_coef * S_inv(b,a) * S_inv(j,i) -
602 : alpha * S_inv(b,i) * S_inv(j,a);
603 :
604 : // Chain rule to change d/dS to d/dR
605 : d2beta_dR2 += d2beta_dS2 * dS_dR_l(a, b) * dS_dR_p(i,
606 : j); d2mu_dR2 += d2mu_dS2 * dS_dR_l(a, b) * dS_dR_p(i, j);
607 :
608 : }// for b
609 : }// for a
610 : }// for j
611 : }// for i, end tensor contraction
612 :
613 : // Jacobian contribution
614 : K[var_id1][var_id2](l, p) += quad_weight * ((1. -
615 : _dilation_weight) * d2beta_dR2 + _dilation_weight * d2mu_dR2);
616 :
617 : }// for var_id2
618 : }// for p
619 : }// for var_id1
620 : }// for l
621 :
622 : End pseudo-code, begin efficient code
623 :
624 : */
625 :
626 15992080 : for (const auto l: elem.node_index_range()) // Contribution to Hessian from node l
627 : {
628 57252864 : for (const auto var_id1 : make_range(dim)) // Contribution from each x/y/z component of node l
629 : {
630 : // Build dS/dR_l, the derivative of the physical-to-target mapping Jacobin w.r.t.
631 : // the l-th node
632 42860536 : RealTensor dS_dR_l = RealTensor(0);
633 170817168 : for (const auto ii : make_range(dim))
634 170608160 : dS_dR_l(var_id1, ii) = dphi_maps[ii][l][qp];
635 46432162 : dS_dR_l *= target_jacobians[qp];
636 :
637 : // Jacobian is symmetric, only need to loop over lower triangular portion
638 343146728 : for (const auto p: make_range(l + 1)) // Contribution to Hessian from node p
639 : {
640 1199363232 : for (const auto var_id2 : make_range(dim)) // Contribution from each x/y/z component of node p
641 : {
642 : // Build dS/dR_l, the derivative of the physical-to-target mapping Jacobin w.r.t.
643 : // the p-th node
644 899077040 : RealTensor dS_dR_p = RealTensor(0);
645 3592759776 : for (const auto jj : make_range(dim))
646 3591573568 : dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
647 973999700 : dS_dR_p *= target_jacobians[qp];
648 :
649 74922660 : Real d2E_dR2 = 0.;
650 : // Perform tensor contraction
651 3592759776 : for (const auto i : make_range(dim))
652 : {
653 :
654 10767648864 : for (const auto j : make_range(dim))
655 : {
656 :
657 8073966128 : const auto S_ij = S(i,j);
658 8073966128 : const auto S_inv_ji = S_inv(j,i);
659 :
660 : // Apply the stuff only depending on i and j before entering a, b loops
661 : // beta
662 : const std::vector<Real> d2beta_dS2_coefs_ij_applied{
663 8073966128 : d2beta_dS2_coefs_times_distortion_weight[1] * S_ij,
664 8746794932 : d2beta_dS2_coefs_times_distortion_weight[2] * S_ij,
665 8746794932 : d2beta_dS2_coefs_times_distortion_weight[3] * S_inv_ji,
666 8746794932 : d2beta_dS2_coefs_times_distortion_weight[4] * S_inv_ji,
667 8073966128 : };
668 : // mu
669 8073966128 : const auto dalpha_dS_coef_ij_applied = dalpha_dS_coef_times_dilation_weight * S_inv_ji;
670 :
671 672828804 : Real d2E_dSdR_l = 0.;
672 672828804 : Real d2E_dSdR_p = 0.;
673 :
674 24214823648 : for (const auto a : make_range(i + 1))
675 : {
676 :
677 : // If this condition is met, both the ijab and abij
678 : // contributions to the Jacobian are zero due to the
679 : // spasity patterns of dS_dR_l and dS_dR_p and this
680 : // iteration may be skipped
681 16140857520 : if (!(a == var_id1 && i == var_id2) &&
682 13522834512 : !(a == var_id2 && i == var_id1))
683 12326578592 : continue;
684 :
685 2693682736 : const auto S_inv_ja = S_inv(j,a);
686 :
687 2693682736 : const Real d2beta_dS2_coef_ia_applied = d2beta_dS2_coefs_times_distortion_weight[0] * I(i,a);
688 2693682736 : const Real d2beta_dS2_coef_ja_applied = d2beta_dS2_coefs_times_distortion_weight[5] * S_inv_ja;
689 2693682736 : const Real alpha_ja_applied = alpha_times_dilation_weight * S_inv_ja;
690 :
691 2693682736 : const auto b_limit = (a == i) ? j + 1 : dim;
692 9870346016 : for (const auto b : make_range(b_limit))
693 : {
694 :
695 : // Combine precomputed coefficients with tensor products
696 : // to get d2(beta) / dS2
697 : Real d2beta_dS2_times_distortion_weight = (
698 7774717060 : d2beta_dS2_coef_ia_applied * I(j,b) +
699 7176663280 : d2beta_dS2_coefs_ij_applied[0] * S(a,b) +
700 7176663280 : d2beta_dS2_coefs_ij_applied[1] * S_inv(b,a) +
701 7774717060 : d2beta_dS2_coefs_ij_applied[2] * S(a,b) +
702 7176663280 : d2beta_dS2_coefs_ij_applied[3] * S_inv(b,a) +
703 7176663280 : d2beta_dS2_coef_ja_applied * S_inv(b,i)
704 7176663280 : );
705 :
706 : // Incorporate tensor product portion to get d2(mu) /
707 : // dS2
708 7176663280 : const Real d2mu_dS2_times_dilation_weight = dalpha_dS_coef_ij_applied * S_inv(b,a) - alpha_ja_applied * S_inv(b,i);
709 :
710 :
711 : // Chain rule to change d/dS to d/dR
712 7176663280 : const auto d2E_dS2 =
713 : d2beta_dS2_times_distortion_weight +
714 : d2mu_dS2_times_dilation_weight;
715 :
716 : // if !(a == var_id1 (next line) && i == var_id2
717 : // (outside 'a' loop)), dS_dR_l(p) multiplier is zero
718 7176663280 : d2E_dSdR_l += d2E_dS2 * dS_dR_l(a, b);
719 :
720 7176663280 : if (!(i == a && j == b))
721 : // if !(a == var_id2 (next line) && i == var_id1
722 : // (outside 'a' loop)), dS_dR_p(l) multiplier is zero
723 6277586240 : d2E_dSdR_p += d2E_dS2 * dS_dR_p(a, b);
724 :
725 : } // for b
726 : } // for a
727 8073966128 : d2E_dR2 +=
728 8073966128 : d2E_dSdR_l * dS_dR_p(i, j) + d2E_dSdR_p * dS_dR_l(i, j);
729 : }// for j
730 : }// for i, end tensor contraction
731 :
732 : // Jacobian contribution
733 899077040 : const Real jacobian_contribution = quad_weight * d2E_dR2;
734 899077040 : K[var_id1][var_id2](l, p) += jacobian_contribution;
735 : // Jacobian is symmetric, add contribution to p,l entry
736 : // Don't get the diagonal twice!
737 899077040 : if (p < l)
738 : // Note the transposition of var_id1 and var_id2 as these are also jacobian indices
739 835380186 : K[var_id2][var_id1](p, l) += jacobian_contribution;
740 :
741 : }// for var_id2
742 : }// for p
743 : }// for var_id1
744 : }// for l
745 : }
746 2709000 : } // end of the quadrature point qp-loop
747 :
748 383672 : return request_jacobian;
749 2709000 : }
750 :
751 2698 : const MeshQualityInfo & VariationalSmootherSystem::get_mesh_info()
752 : {
753 2698 : if (!_mesh_info.initialized)
754 2556 : compute_mesh_quality_info();
755 :
756 2698 : return _mesh_info;
757 : }
758 :
759 35002 : void VariationalSmootherSystem::compute_mesh_quality_info()
760 : {
761 : // If the reference volume has not yet been computed, compute it.
762 35982 : if (std::abs(_ref_vol) < TOLERANCE * TOLERANCE)
763 0 : prepare_for_smoothing();
764 :
765 35982 : std::unique_ptr<DiffContext> con = this->build_context();
766 980 : FEMContext & femcontext = cast_ref<FEMContext &>(*con);
767 35002 : this->init_context(femcontext);
768 :
769 1960 : const auto & mesh = this->get_mesh();
770 35002 : const auto dim = mesh.mesh_dimension();
771 35002 : const Real half_dim = 0.5 * dim;
772 35002 : const auto distortion_weight = 1. - _dilation_weight;
773 :
774 : // Make pre-requests before reinit() for efficiency in
775 : // --enable-deprecated builds, and to avoid errors in
776 : // --disable-deprecated builds.
777 980 : const auto & fe_map = femcontext.get_element_fe(0)->get_fe_map();
778 980 : const auto & quad_weights = femcontext.get_element_qrule().get_weights();
779 980 : const auto & JxW = fe_map.get_JxW();
780 980 : fe_map.get_dxyzdxi();
781 980 : fe_map.get_dxyzdeta();
782 980 : fe_map.get_dxyzdzeta();
783 :
784 35002 : MeshQualityInfo info;
785 :
786 838344 : for (const auto * elem : mesh.active_local_element_ptr_range())
787 : {
788 419600 : femcontext.pre_fe_reinit(*this, elem);
789 419600 : femcontext.elem_fe_reinit();
790 :
791 : // Element-integrated quantities
792 34940 : Real det_S_int = 0.;
793 34940 : Real beta_int = 0.;
794 34940 : Real mu_int = 0.;
795 34940 : Real combined_int = 0.;
796 :
797 419600 : const auto & target_jacobian_dets = _target_jacobian_dets[elem->type()];
798 :
799 3971680 : for (const auto qp : index_range(JxW))
800 : {
801 3847980 : det_S_int += JxW[qp] / _target_jacobian_dets[elem->type()][qp];
802 3847980 : const auto quad_weight = quad_weights[qp] / target_jacobian_dets[qp];
803 :
804 : // Grab the physical-to-reference mapping Jacobian matrix (i.e., "S") at this qp
805 3552080 : RealTensor S = get_jacobian_at_qp(fe_map, dim, qp);
806 :
807 : // Apply target element transformation to get the physical-to-target jacobian
808 3552080 : S *= _target_jacobians[elem->type()][qp];
809 :
810 : // Determinant of S
811 3552080 : const auto det = S.det();
812 3552080 : const auto det_sq = det * det;
813 :
814 3552080 : if (det > info.max_qp_det_S)
815 110122 : info.max_qp_det_S = det;
816 3441958 : else if (det < info.min_qp_det_S)
817 130512 : info.min_qp_det_S = det;
818 :
819 3552080 : if (det < TOLERANCE * TOLERANCE)
820 1332 : info.mesh_is_tangled = true;
821 :
822 : // trace of S^T * S
823 3552080 : const auto tr = trace(S.transpose() * S, dim);
824 :
825 : // The chi function allows us to handle degenerate elements
826 3552080 : const auto chi = chi_epsilon(det, _epsilon_squared_assembly);
827 :
828 : // distortion
829 3552080 : const Real beta = std::pow(tr / dim, half_dim) / chi;
830 3552080 : beta_int += beta * quad_weight;
831 :
832 : // dilation
833 3552080 : const Real mu = 0.5 * (_ref_vol + det_sq / _ref_vol) / chi;
834 3552080 : mu_int += mu * quad_weight;
835 :
836 : // combined
837 3552080 : const Real E = distortion_weight * beta + _dilation_weight * mu;
838 3552080 : combined_int += E * quad_weight;
839 : }
840 :
841 419600 : info.total_det_S += det_S_int;
842 419600 : if (det_S_int > info.max_elem_det_S.first)
843 5545 : info.max_elem_det_S = std::make_pair(det_S_int, elem->id());
844 419600 : if (det_S_int < info.min_elem_det_S.first)
845 6692 : info.min_elem_det_S = std::make_pair(det_S_int, elem->id());
846 :
847 419600 : info.total_distortion += beta_int;
848 419600 : if (beta_int > info.max_elem_distortion.first)
849 5745 : info.max_elem_distortion = std::make_pair(beta_int, elem->id());
850 419600 : if (beta_int < info.min_elem_distortion.first)
851 6147 : info.min_elem_distortion = std::make_pair(beta_int, elem->id());
852 :
853 419600 : info.total_dilation += mu_int;
854 419600 : if (mu_int > info.max_elem_dilation.first)
855 5181 : info.max_elem_dilation = std::make_pair(mu_int, elem->id());
856 419600 : if (mu_int < info.min_elem_dilation.first)
857 5986 : info.min_elem_dilation = std::make_pair(mu_int, elem->id());
858 :
859 419600 : info.total_combined += combined_int;
860 419600 : if (combined_int > info.max_elem_combined.first)
861 5627 : info.max_elem_combined = std::make_pair(combined_int, elem->id());
862 419600 : if (combined_int < info.min_elem_combined.first)
863 6130 : info.min_elem_combined = std::make_pair(combined_int, elem->id());
864 :
865 419600 : if (_verbosity > 90)
866 : {
867 51216 : libMesh::out << "Elem " << elem->id() << " quality:" << std::endl
868 17072 : << " distortion-dilation metric: " << combined_int << std::endl
869 17072 : << " distortion metric: " << beta_int << std::endl
870 17072 : << " dilation metric: " << mu_int << std::endl
871 17072 : << " det(S): " << det_S_int << std::endl;
872 : }
873 :
874 33042 : } // for elem
875 :
876 : // Get contributions from elements on other processors
877 35002 : communicate_pair_max(info.max_elem_det_S, mesh.comm());
878 35002 : communicate_pair_min(info.min_elem_det_S, mesh.comm());
879 35002 : mesh.comm().max(info.max_qp_det_S);
880 35002 : mesh.comm().min(info.min_qp_det_S);
881 35002 : mesh.comm().sum(info.total_det_S);
882 :
883 35002 : communicate_pair_max(info.max_elem_distortion, mesh.comm());
884 35002 : communicate_pair_min(info.min_elem_distortion, mesh.comm());
885 35002 : mesh.comm().sum(info.total_distortion);
886 :
887 35002 : communicate_pair_max(info.max_elem_dilation, mesh.comm());
888 35002 : communicate_pair_min(info.min_elem_dilation, mesh.comm());
889 35002 : mesh.comm().sum(info.total_dilation);
890 :
891 35002 : communicate_pair_max(info.max_elem_combined, mesh.comm());
892 35002 : communicate_pair_min(info.min_elem_combined, mesh.comm());
893 35002 : mesh.comm().sum(info.total_combined);
894 :
895 35002 : mesh.comm().max(info.mesh_is_tangled);
896 :
897 35002 : info.initialized = true;
898 :
899 35002 : _mesh_info = info;
900 :
901 35002 : if (_verbosity > 50)
902 44 : libMesh::out << info;
903 35002 : }
904 :
905 2094 : void VariationalSmootherSystem::get_target_to_reference_jacobian(
906 : const Elem * const target_elem,
907 : const FEMContext & femcontext,
908 : std::vector<RealTensor> & jacobians,
909 : std::vector<Real> & jacobian_dets)
910 : {
911 :
912 2094 : const auto dim = target_elem->dim();
913 :
914 78 : const auto & qrule_points = femcontext.get_element_qrule().get_points();
915 78 : const auto & qrule_weights = femcontext.get_element_qrule().get_weights();
916 78 : const auto nq_points = femcontext.get_element_qrule().n_points();
917 :
918 : // If the target element is the reference element, Jacobian matrix is
919 : // identity, det of inverse is 1. These will only be overwritten if a
920 : // different target element is explicitly specified.
921 2172 : jacobians = std::vector<RealTensor>(nq_points, RealTensor(
922 : 1., 0., 0.,
923 : 0., 1., 0.,
924 78 : 0., 0., 1.));
925 2094 : jacobian_dets = std::vector<Real>(nq_points, 1.0);
926 :
927 : // Don't use "if (*target_elem == *(target_elem->reference_elem()))" here, it
928 : // only compares global node ids, not the node locations themselves.
929 78 : bool target_equals_reference = true;
930 2094 : const auto * ref_elem = target_elem->reference_elem();
931 23365 : for (const auto local_id : make_range(target_elem->n_nodes()))
932 21271 : target_equals_reference &= target_elem->node_ref(local_id) == ref_elem->node_ref(local_id);
933 2094 : if (target_equals_reference)
934 802 : return;
935 :
936 : // Create FEMap to compute target_element mapping information
937 1380 : FEMap fe_map_target;
938 :
939 : // pre-request mapping derivatives
940 44 : fe_map_target.get_dxyzdxi();
941 44 : fe_map_target.get_dxyzdeta();
942 44 : fe_map_target.get_dxyzdzeta();
943 :
944 : // build map
945 1292 : fe_map_target.init_reference_to_physical_map(dim, qrule_points, target_elem);
946 1292 : fe_map_target.compute_map(dim, qrule_weights, target_elem, /*d2phi=*/false);
947 :
948 27805 : for (const auto qp : make_range(nq_points))
949 : {
950 : // We use Larisa's H notation to denote the reference-to-target jacobian
951 26513 : RealTensor H = get_jacobian_at_qp(fe_map_target, dim, qp);
952 :
953 : // The target-to-reference jacobian is the inverse of the
954 : // reference-to-target jacobian
955 26513 : jacobians[qp] = H.inverse();
956 27567 : jacobian_dets[qp] = jacobians[qp].det();
957 : }
958 1204 : }
959 :
960 1846 : std::ostream &operator<<(std::ostream &os, const MeshQualityInfo & info)
961 : {
962 1794 : os << "Mesh quality info:" << std::endl
963 1794 : << " Mesh distortion-dilation metric: "
964 1846 : << info.total_combined << std::endl
965 1794 : << " Mesh distortion metric: "
966 1846 : << info.total_distortion << std::endl
967 1794 : << " Mesh dilation metric: "
968 1846 : << info.total_dilation << std::endl
969 1794 : << " Max distortion-dilation is in elem "
970 1846 : << info.max_elem_combined.second << ": "
971 1846 : << info.max_elem_combined.first << std::endl
972 1794 : << " Max distortion is in elem "
973 1846 : << info.max_elem_distortion.second << ": "
974 1846 : << info.max_elem_distortion.first << std::endl
975 1794 : << " Max dilation is in elem "
976 1846 : << info.max_elem_dilation.second << ": "
977 1846 : << info.max_elem_dilation.first << std::endl
978 1794 : << " Max det(S) is in elem "
979 1846 : << info.max_elem_det_S.second << ": "
980 1846 : << info.max_elem_det_S.first << std::endl
981 1794 : << " Min distortion-dilation is in elem "
982 1846 : << info.min_elem_combined.second << ": "
983 1846 : << info.min_elem_combined.first << std::endl
984 1794 : << " Min distortion is in elem "
985 1846 : << info.min_elem_distortion.second << ": "
986 1846 : << info.min_elem_distortion.first << std::endl
987 1794 : << " Min dilation is in elem "
988 1846 : << info.min_elem_dilation.second << ": "
989 1846 : << info.min_elem_dilation.first << std::endl
990 1794 : << " Min det(S) is in elem "
991 1846 : << info.min_elem_det_S.second << ": "
992 1846 : << info.min_elem_det_S.first << std::endl
993 1846 : << " Max qp det(S): " << info.max_qp_det_S << std::endl
994 1846 : << " Min qp det(S): " << info.min_qp_det_S << std::endl
995 1846 : << " Mesh-integrated det(S): " << info.total_det_S << std::endl
996 1846 : << " Tangled: " << info.mesh_is_tangled << std::endl;
997 :
998 1846 : return os;
999 : }
1000 :
1001 : } // namespace libMesh
|
