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