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