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 578628 : 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 578628 : comm.minloc(pair.first, rank);
49 578628 : comm.broadcast(pair.second, rank);
50 578628 : }
51 :
52 : /*
53 : * Gets the dof_id_type value corresponding to the maximum of the Real value.
54 : */
55 578628 : 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 578628 : comm.maxloc(pair.first, rank);
60 578628 : comm.broadcast(pair.second, rank);
61 578628 : }
62 :
63 : /**
64 : * Function to prevent dividing by zero for degenerate elements
65 : */
66 19561504 : Real chi_epsilon(const Real & x, const Real epsilon_squared)
67 : {
68 21196776 : 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 19579428 : RealTensor get_jacobian_at_qp(const FEMap & fe_map,
76 : const unsigned int & dim,
77 : const unsigned int & qp)
78 : {
79 19579428 : const auto & dxyzdxi = fe_map.get_dxyzdxi()[qp];
80 3261388 : const auto & dxyzdeta = fe_map.get_dxyzdeta()[qp];
81 3261388 : 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 19579428 : switch (dim)
87 : {
88 6528 : case 1:
89 : return RealTensor(
90 : dxyzdxi(0), 0, 0,
91 : 0, 1, 0,
92 : 0, 0, 1
93 544 : );
94 : break;
95 :
96 428905 : case 2:
97 : return RealTensor(
98 : dxyzdxi(0), dxyzdeta(0), 0,
99 : dxyzdxi(1), dxyzdeta(1), 0,
100 : 0, 0, 1
101 35502 : );
102 : break;
103 :
104 19143995 : case 3:
105 3178544 : return RealTensor(
106 : dxyzdxi,
107 : dxyzdeta,
108 : dxyzdzeta
109 1589272 : ).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 19561504 : Real trace(const RealTensor & A, const unsigned int & dim)
121 : {
122 1624520 : Real tr = 0.0;
123 77804992 : for (const auto i : make_range(dim))
124 58243488 : tr += A(i, i);
125 :
126 19561504 : return tr;
127 : }
128 :
129 5025 : VariationalSmootherSystem::~VariationalSmootherSystem () = default;
130 :
131 140965 : 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 7980 : auto & mesh = this->get_mesh();
138 140965 : this->solution->close();
139 :
140 3706430 : for (auto * node : mesh.local_node_ptr_range())
141 : {
142 7420030 : for (const auto d : make_range(mesh.mesh_dimension()))
143 : {
144 5549370 : const auto dof_id = node->dof_number(this->number(), d, 0);
145 : // Update mesh
146 5549370 : (*node)(d) = libmesh_real((*current_local_solution)(dof_id));
147 : }
148 132985 : }
149 :
150 140965 : SyncNodalPositions sync_object(mesh);
151 277924 : 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 140965 : compute_mesh_quality_info();
155 140965 : 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 140965 : 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 134653 : _epsilon_squared_assembly = 0.;
169 :
170 140965 : FEMSystem::assembly(get_residual, get_jacobian, apply_heterogeneous_constraints, apply_no_constraints);
171 :
172 140965 : 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 140965 : }
183 :
184 1775 : void VariationalSmootherSystem::solve()
185 : {
186 1775 : const auto & mesh_info = get_mesh_info();
187 1775 : if (mesh_info.mesh_is_tangled)
188 : {
189 : // Untangling solve
190 142 : _untangling_solve = true;
191 :
192 : // Untangling seems to work better using only the distortion metric.
193 : // For a mixed metric, I've seen it *technically* untangle a mesh to a
194 : // still-suboptimal mesh (i.e., a local minima), but not been able to
195 : // smooth it well because the smoothest solution is on the other side of
196 : // a tangulation barrier.
197 142 : const auto dilation_weight = _dilation_weight;
198 142 : _dilation_weight = 0.;
199 :
200 142 : FEMSystem::solve();
201 :
202 : // Reset the dilation weight
203 142 : _dilation_weight = dilation_weight;
204 : }
205 :
206 : // Smoothing solve
207 1775 : _untangling_solve = false;
208 1775 : FEMSystem::solve();
209 1775 : libmesh_error_msg_if(get_mesh_info().mesh_is_tangled, "The smoothing solve tangled the mesh!");
210 1775 : }
211 :
212 1775 : void VariationalSmootherSystem::init_data ()
213 : {
214 100 : auto & mesh = this->get_mesh();
215 50 : const auto & elem_orders = mesh.elem_default_orders();
216 1775 : libmesh_error_msg_if(elem_orders.size() != 1,
217 : "The variational smoother cannot be used for mixed-order meshes!");
218 1775 : const auto fe_order = *elem_orders.begin();
219 : // Add a variable for each dimension of the mesh
220 : // "r0" for x, "r1" for y, "r2" for z
221 6248 : for (const auto & d : make_range(mesh.mesh_dimension()))
222 8820 : this->add_variable ("r" + std::to_string(d), fe_order);
223 :
224 : // Do the parent's initialization after variables are defined
225 1775 : FEMSystem::init_data();
226 :
227 : // Set the current_local_solution to the current mesh for the initial guess
228 1775 : this->solution->close();
229 :
230 43342 : for (auto * node : mesh.local_node_ptr_range())
231 : {
232 84768 : for (const auto d : make_range(mesh.mesh_dimension()))
233 : {
234 63036 : const auto dof_id = node->dof_number(this->number(), d, 0);
235 : // Update solution
236 63036 : (*solution).set(dof_id, (*node)(d));
237 : }
238 1675 : }
239 :
240 1775 : this->prepare_for_smoothing();
241 1775 : }
242 :
243 1775 : void VariationalSmootherSystem::prepare_for_smoothing()
244 : {
245 : // If this method has already been called to set _ref_vol, just return.
246 1825 : if (std::abs(_ref_vol) > TOLERANCE * TOLERANCE)
247 0 : return;
248 :
249 1775 : std::unique_ptr<DiffContext> con = this->build_context();
250 50 : FEMContext & femcontext = cast_ref<FEMContext &>(*con);
251 1775 : this->init_context(femcontext);
252 :
253 100 : const auto & mesh = this->get_mesh();
254 :
255 1775 : Real elem_averaged_det_S_sum = 0.;
256 :
257 : // Make pre-requests before reinit() for efficiency in
258 : // --enable-deprecated builds, and to avoid errors in
259 : // --disable-deprecated builds.
260 50 : const auto & fe_map = femcontext.get_element_fe(0)->get_fe_map();
261 50 : const auto & JxW = fe_map.get_JxW();
262 :
263 14084 : for (const auto * elem : mesh.active_local_element_ptr_range())
264 : {
265 10584 : femcontext.pre_fe_reinit(*this, elem);
266 10584 : femcontext.elem_fe_reinit();
267 :
268 : // Add target element info, if applicable
269 11466 : if (_target_jacobians.find(elem->type()) == _target_jacobians.end())
270 : {
271 1314 : const auto [target_elem, target_nodes] = get_target_elem(elem->type());
272 1314 : get_target_to_reference_jacobian(target_elem.get(),
273 : femcontext,
274 1314 : _target_jacobians[elem->type()],
275 2628 : _target_jacobian_dets[elem->type()]);
276 : }// if find == end()
277 :
278 : // Reference volume computation
279 882 : Real elem_integrated_det_S = 0.;
280 134376 : for (const auto qp : index_range(JxW))
281 134108 : elem_integrated_det_S += JxW[qp] * _target_jacobian_dets[elem->type()][qp];
282 10584 : const auto ref_elem_vol = elem->reference_elem()->volume();
283 10584 : elem_averaged_det_S_sum += elem_integrated_det_S / ref_elem_vol;
284 :
285 1675 : } // for elem
286 :
287 : // Get contributions from elements on other processors
288 1775 : mesh.comm().sum(elem_averaged_det_S_sum);
289 :
290 1775 : _ref_vol = elem_averaged_det_S_sum / mesh.n_active_elem();
291 1675 : }
292 :
293 295309 : void VariationalSmootherSystem::init_context(DiffContext & context)
294 : {
295 10796 : FEMContext & c = cast_ref<FEMContext &>(context);
296 :
297 10796 : FEBase * my_fe = nullptr;
298 :
299 : // Now make sure we have requested all the data
300 : // we need to build the system.
301 :
302 : // We might have a multi-dimensional mesh
303 : const std::set<unsigned char> & elem_dims =
304 10796 : c.elem_dimensions();
305 :
306 590618 : for (const auto & dim : elem_dims)
307 : {
308 295309 : c.get_element_fe( 0, my_fe, dim );
309 10796 : my_fe->get_nothing();
310 :
311 10796 : auto & fe_map = my_fe->get_fe_map();
312 10796 : fe_map.get_dxyzdxi();
313 10796 : fe_map.get_dxyzdeta();
314 10796 : fe_map.get_dxyzdzeta();
315 10796 : fe_map.get_JxW();
316 :
317 : // Mesh may be tangled, allow negative Jacobians
318 10796 : fe_map.set_jacobian_tolerance(std::numeric_limits<Real>::lowest());
319 :
320 295309 : c.get_side_fe( 0, my_fe, dim );
321 10796 : my_fe->get_nothing();
322 : }
323 :
324 295309 : FEMSystem::init_context(context);
325 295309 : }
326 :
327 :
328 817286 : bool VariationalSmootherSystem::element_time_derivative (bool request_jacobian,
329 : DiffContext & context)
330 : {
331 67364 : FEMContext & c = cast_ref<FEMContext &>(context);
332 :
333 135746 : const Elem & elem = c.get_elem();
334 :
335 817286 : unsigned int dim = c.get_dim();
336 :
337 67364 : unsigned int x_var = 0, y_var = 1, z_var = 2;
338 : // In the case of lower dimensions, we will not access z (1D/2D) or y (1D)
339 817286 : if (dim < 3)
340 : {
341 4148 : z_var = 0;
342 50064 : if (dim < 2)
343 100 : y_var = 0;
344 : }
345 :
346 : // The subvectors and submatrices we need to fill:
347 : // system residual
348 : std::reference_wrapper<DenseSubVector<Number>> F[3] =
349 : {
350 : c.get_elem_residual(x_var),
351 : c.get_elem_residual(y_var),
352 : c.get_elem_residual(z_var)
353 67364 : };
354 : // system jacobian
355 : std::reference_wrapper<DenseSubMatrix<Number>> K[3][3] =
356 : {
357 : {c.get_elem_jacobian(x_var, x_var), c.get_elem_jacobian(x_var, y_var), c.get_elem_jacobian(x_var, z_var)},
358 : {c.get_elem_jacobian(y_var, x_var), c.get_elem_jacobian(y_var, y_var), c.get_elem_jacobian(y_var, z_var)},
359 : {c.get_elem_jacobian(z_var, x_var), c.get_elem_jacobian(z_var, y_var), c.get_elem_jacobian(z_var, z_var)}
360 67364 : };
361 :
362 : // Quadrature info
363 817286 : const auto quad_weights = c.get_element_qrule().get_weights();
364 :
365 817286 : const auto distortion_weight = 1. - _dilation_weight;
366 :
367 : // Get some references to cell-specific data that
368 : // will be used to assemble the linear system.
369 :
370 67364 : const auto & fe_map = c.get_element_fe(0)->get_fe_map();
371 :
372 67364 : const auto & dphidxi_map = fe_map.get_dphidxi_map();
373 67364 : const auto & dphideta_map = fe_map.get_dphideta_map();
374 67364 : const auto & dphidzeta_map = fe_map.get_dphidzeta_map();
375 :
376 : // Integrate the distortion-dilation metric over the reference element
377 10472566 : 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 : // Note that the term "Jacobian" has two meanings in this mesh smoothing
388 : // application. The first meaning refers to the Jacobian w.r.t R of the
389 : // gradient w.r.t. R, (i.e., the Hessian of the metric). This is the K
390 : // variable defined above. This is also the Jacobian the 'request_jacobian'
391 : // variable refers to. The second mesning refers to the Jacobian of the
392 : // physical-to-reference element mapping. This is the Jacobian used to
393 : // compute the distorion-dilation metric.
394 : //
395 : // Grab the physical-to-reference mapping Jacobian matrix (i.e., "S") at this qp
396 10457084 : RealTensor S = get_jacobian_at_qp(fe_map, dim, qp);
397 :
398 : // Apply target element transformation
399 9655280 : S *= _target_jacobians[elem.type()][qp];
400 :
401 : // Compute quantities needed for the smoothing algorithm
402 :
403 : // determinant
404 9655280 : const Real det = S.det();
405 9655280 : const Real det_sq = det * det;
406 9655280 : const Real det_cube = det_sq * det;
407 :
408 9655280 : const Real ref_vol_sq = _ref_vol * _ref_vol;
409 :
410 : // trace of S^T * S
411 : // DO NOT USE RealTensor.tr for the trace, it will NOT be correct for
412 : // 1D and 2D meshes because of our hack of putting 1s in the diagonal of
413 : // S for the extra dimensions (see get_jacobian_at_qp)
414 9655280 : const auto tr = trace(S.transpose() * S, dim);
415 9655280 : const Real tr_div_dim = tr / dim;
416 :
417 : // Precompute pow(tr_div_dim, 0.5 * dim - x) for x = 0, 1, 2
418 9655280 : const Real half_dim = 0.5 * dim;
419 : const std::vector<Real> trace_powers{
420 9655280 : std::pow(tr_div_dim, half_dim),
421 9655280 : std::pow(tr_div_dim, half_dim - 1.),
422 9655280 : std::pow(tr_div_dim, half_dim - 2.),
423 10457084 : };
424 :
425 : // inverse of S
426 10457084 : const RealTensor S_inv = S.inverse();
427 : // inverse transpose of S
428 1603608 : const RealTensor S_inv_T = S_inv.transpose();
429 :
430 : // Identity matrix
431 16092592 : const RealTensor I(1, 0, 0,
432 16092592 : 0, 1, 0,
433 9649904 : 0, 0, 1);
434 :
435 : // The chi function allows us to handle degenerate elements
436 9655280 : const auto chi = chi_epsilon(det, _epsilon_squared_assembly);
437 9655280 : const Real chi_sq = chi * chi;
438 9655280 : const Real sqrt_term = std::sqrt(_epsilon_squared_assembly + det_sq);
439 : // dchi(x) / dx
440 9655280 : const Real chi_prime = 0.5 * (1. + det / sqrt_term);
441 9655280 : const Real chi_prime_sq = chi_prime * chi_prime;
442 : // d2chi(x) / dx2
443 9655280 : const Real chi_2prime = 0.5 * (1. / sqrt_term - det_sq / Utility::pow<3>(sqrt_term));
444 :
445 : // Distortion metric (beta)
446 : //const Real beta = trace_powers[0] / chi;
447 12873248 : const RealTensor dbeta_dS = (trace_powers[1] / chi) * S - (trace_powers[0] / chi_sq * chi_prime * det) * S_inv_T;
448 :
449 : // Dilation metric (mu)
450 : //const Real mu = 0.5 * (_ref_vol + det_sq / _ref_vol) / chi;
451 : // We represent d mu / dS as alpha(S) * S^-T, where alpha is a scalar function
452 9655280 : const Real alpha = (-chi_prime * det_cube + 2. * det_sq * chi - ref_vol_sq * det * chi_prime) / (2. * _ref_vol * chi_sq);
453 9649904 : const RealTensor dmu_dS = alpha * S_inv_T;
454 :
455 : // Combined metric (E)
456 : //const Real E = distortion_weight * beta + _dilation_weight * mu;
457 12066068 : const RealTensor dE_dS = distortion_weight * dbeta_dS + _dilation_weight * dmu_dS;
458 :
459 : // This vector is useful in computing dS/dR below
460 41839088 : std::vector<std::vector<std::vector<Real>>> dphi_maps = {dphidxi_map, dphideta_map, dphidzeta_map};
461 :
462 : // Compute residual (i.e., the gradient of the combined metric w.r.t node locations)
463 : // Recall that when the gradient (residual) is zero, the combined metric is minimized
464 117949280 : for (const auto l : elem.node_index_range())
465 : {
466 432399456 : for (const auto var_id : make_range(dim))
467 : {
468 : // Build dS/dR, the derivative of the physical-to-reference mapping Jacobin w.r.t. the mesh node locations
469 297020924 : RealTensor dS_dR = RealTensor(0);
470 1294884576 : for (const auto jj : make_range(dim))
471 1295284160 : dS_dR(var_id, jj) = dphi_maps[jj][l][qp];
472 :
473 : // Residual contribution. The contraction of dE/dS and dS/dR gives us
474 : // the gradient we are looking for, dE/dR
475 675268804 : F[var_id](l) += quad_weights[qp] * dE_dS.contract(dS_dR);
476 : }// for var_id
477 : }// for l
478 :
479 :
480 9655280 : if (request_jacobian)
481 : {
482 : // Compute jacobian of the smoothing system (i.e., the Hessian of the
483 : // combined metric w.r.t. the mesh node locations)
484 :
485 : // Precompute coefficients to be applied to each component of tensor
486 : // products in the loops below. At first glance, these coefficients look
487 : // like gibberish, but everything in the Hessian has been verified by
488 : // taking finite differences of the gradient. We should probably write
489 : // down the derivations of the gradient and Hessian somewhere...
490 :
491 : // Recall that above, dbeta_dS takes the form:
492 : // d(beta)/dS = c1(S) * S - c2(S) * S_inv_T,
493 : // where c1 and c2 are scalar-valued functions.
494 : const std::vector<Real> d2beta_dS2_coefs_times_distortion_weight = {
495 : //Part 1: scaler coefficients of d(c1 * S) / dS
496 : //
497 : // multiplies I[i,a] x I[j,b]
498 3008296 : (trace_powers[1] / chi) * distortion_weight,
499 : // multiplies S[a,b] x S[i,j]
500 3258334 : (((dim - 2.) / dim) * trace_powers[2] / chi) * distortion_weight,
501 : // multiplies S_inv[b,a] * S[i,j]
502 3258334 : (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
503 : //
504 : //Part 2: scaler coefficients of d(-c2 * S_inv_T) / dS
505 : //
506 : // multiplies S[a,b] x S_inv[j,i]
507 3258334 : (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
508 : // multiplies S_inv[b,a] x S_inv[j,i]
509 3258334 : (trace_powers[0] * (det / chi_sq)
510 3008296 : * ((2. * chi_prime_sq / chi - chi_2prime) * det - chi_prime)) * distortion_weight,
511 : // multiplies S_inv[b,i] x S_inv[j,a]
512 3258334 : ((trace_powers[0] / chi_sq) * chi_prime * det) * distortion_weight,
513 3258334 : };
514 :
515 : // d alpha / dS has the form c(S) * S^-T, where c is the scalar coefficient defined below
516 3509060 : const Real dalpha_dS_coef_times_dilation_weight = ((det / (2. * _ref_vol * chi_sq))
517 3509060 : * (-4. * _ref_vol * alpha * chi * chi_prime - chi_2prime * det_cube
518 3509060 : - chi_prime * det_sq + 4 * chi * det - ref_vol_sq * (chi_prime + det * chi_2prime))) * _dilation_weight;
519 :
520 : // This is also useful to precompute
521 3008296 : const Real alpha_times_dilation_weight = alpha * _dilation_weight;
522 :
523 : /*
524 :
525 : To increase the efficiency of the Jacobian computation, we take
526 : advantage of l-p symmetry, ij-ab symmetry, and the sparsity pattern of
527 : dS_dR. We also factor all possible multipliers out of inner loops for
528 : efficiency. The result is code that is more difficult to read. For
529 : clarity, consult the pseudo-code below in this comment.
530 :
531 : for (const auto l: elem.node_index_range()) // Contribution to Hessian
532 : from node l
533 : {
534 : for (const auto var_id1 : make_range(dim)) // Contribution from each
535 : x/y/z component of node l
536 : {
537 : // Build dS/dR_l, the derivative of the physical-to-reference
538 : mapping Jacobin w.r.t.
539 : // the l-th node
540 : RealTensor dS_dR_l = RealTensor(0);
541 : for (const auto ii : make_range(dim))
542 : dS_dR_l(var_id1, ii) = dphi_maps[ii][l][qp];
543 :
544 : for (const auto p: elem.node_index_range()) // Contribution to
545 : Hessian from node p
546 : {
547 : for (const auto var_id2 : make_range(dim)) // Contribution from
548 : each x/y/z component of node p
549 : {
550 : // Build dS/dR_l, the derivative of the physical-to-reference
551 : mapping Jacobin w.r.t.
552 : // the p-th node
553 : RealTensor dS_dR_p = RealTensor(0);
554 : for (const auto jj : make_range(dim))
555 : dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
556 :
557 : Real d2beta_dR2 = 0.;
558 : Real d2mu_dR2 = 0.;
559 : // Perform tensor contraction
560 : for (const auto i : make_range(dim))
561 : {
562 : for (const auto j : make_range(dim))
563 : {
564 : for (const auto a : make_range(dim))
565 : {
566 : for (const auto b : make_range(dim))
567 : {
568 : // Nasty tensor products to be multiplied by
569 : d2beta_dS2_coefs to get d2(beta) / dS2 const std::vector<Real>
570 : d2beta_dS2_tensor_contributions =
571 : {
572 : I(i,a) * I(j,b),
573 : S(a,b) * S(i,j),
574 : S_inv(b,a) * S(i,j),
575 : S(a,b) * S_inv(j,i),
576 : S_inv(b,a) * S_inv(j,i),
577 : S_inv(j,a) * S_inv(b,i),
578 : };
579 :
580 : // Combine precomputed coefficients with tensor products
581 : to get d2(beta) / dS2 Real d2beta_dS2 = 0.; for (const auto comp_id :
582 : index_range(d2beta_dS2_coefs))
583 : {
584 : const Real contribution = d2beta_dS2_coefs[comp_id] *
585 : d2beta_dS2_tensor_contributions[comp_id]; d2beta_dS2 += contribution;
586 : }
587 :
588 : // Incorporate tensor product portion to get d2(mu) /
589 : dS2 const Real d2mu_dS2 = dalpha_dS_coef * S_inv(b,a) * S_inv(j,i) -
590 : alpha * S_inv(b,i) * S_inv(j,a);
591 :
592 : // Chain rule to change d/dS to d/dR
593 : d2beta_dR2 += d2beta_dS2 * dS_dR_l(a, b) * dS_dR_p(i,
594 : j); d2mu_dR2 += d2mu_dS2 * dS_dR_l(a, b) * dS_dR_p(i, j);
595 :
596 : }// for b
597 : }// for a
598 : }// for j
599 : }// for i, end tensor contraction
600 :
601 : // Jacobian contribution
602 : K[var_id1][var_id2](l, p) += quad_weights[qp] * ((1. -
603 : _dilation_weight) * d2beta_dR2 + _dilation_weight * d2mu_dR2);
604 :
605 : }// for var_id2
606 : }// for p
607 : }// for var_id1
608 : }// for l
609 :
610 : End pseudo-code, begin efficient code
611 :
612 : */
613 :
614 36190016 : for (const auto l: elem.node_index_range()) // Contribution to Hessian from node l
615 : {
616 132356016 : for (const auto var_id1 : make_range(dim)) // Contribution from each x/y/z component of node l
617 : {
618 : // Build dS/dR_l, the derivative of the physical-to-reference mapping Jacobin w.r.t.
619 : // the l-th node
620 99170984 : RealTensor dS_dR_l = RealTensor(0);
621 395963376 : for (const auto ii : make_range(dim))
622 395766816 : dS_dR_l(var_id1, ii) = dphi_maps[ii][l][qp];
623 :
624 : // Jacobian is symmetric, only need to loop over lower triangular portion
625 865478536 : for (const auto p: make_range(l + 1)) // Contribution to Hessian from node p
626 : {
627 3063407136 : for (const auto var_id2 : make_range(dim)) // Contribution from each x/y/z component of node p
628 : {
629 : // Build dS/dR_l, the derivative of the physical-to-reference mapping Jacobin w.r.t.
630 : // the p-th node
631 2297272240 : RealTensor dS_dR_p = RealTensor(0);
632 9184806624 : for (const auto jj : make_range(dim))
633 9185526528 : dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
634 :
635 191754644 : Real d2E_dR2 = 0.;
636 : // Perform tensor contraction
637 9184806624 : for (const auto i : make_range(dim))
638 : {
639 :
640 27543619680 : for (const auto j : make_range(dim))
641 : {
642 :
643 20655915952 : const auto S_ij = S(i,j);
644 20655915952 : const auto S_inv_ji = S_inv(j,i);
645 :
646 : // Apply the stuff only depending on i and j before entering a, b loops
647 : // beta
648 : const std::vector<Real> d2beta_dS2_coefs_ij_applied{
649 20655915952 : d2beta_dS2_coefs_times_distortion_weight[1] * S_ij,
650 22380207972 : d2beta_dS2_coefs_times_distortion_weight[2] * S_ij,
651 22380207972 : d2beta_dS2_coefs_times_distortion_weight[3] * S_inv_ji,
652 22380207972 : d2beta_dS2_coefs_times_distortion_weight[4] * S_inv_ji,
653 20655915952 : };
654 : // mu
655 20655915952 : const auto dalpha_dS_coef_ij_applied = dalpha_dS_coef_times_dilation_weight * S_inv_ji;
656 :
657 1724292020 : Real d2E_dSdR_l = 0.;
658 1724292020 : Real d2E_dSdR_p = 0.;
659 :
660 61960559968 : for (const auto a : make_range(i + 1))
661 : {
662 :
663 : // If this condition is met, both the ijab and abij
664 : // contributions to the Jacobian are zero due to the
665 : // spasity patterns of dS_dR_l and dS_dR_p and this
666 : // iteration may be skipped
667 41304644016 : if (!(a == var_id1 && i == var_id2) &&
668 34605487648 : !(a == var_id2 && i == var_id1))
669 31546458704 : continue;
670 :
671 6887703728 : const auto S_inv_ja = S_inv(j,a);
672 :
673 6887703728 : const Real d2beta_dS2_coef_ia_applied = d2beta_dS2_coefs_times_distortion_weight[0] * I(i,a);
674 6887703728 : const Real d2beta_dS2_coef_ja_applied = d2beta_dS2_coefs_times_distortion_weight[5] * S_inv_ja;
675 6887703728 : const Real alpha_ja_applied = alpha_times_dilation_weight * S_inv_ja;
676 :
677 6887703728 : const auto b_limit = (a == i) ? j + 1 : dim;
678 25248319264 : for (const auto b : make_range(b_limit))
679 : {
680 :
681 : // Combine precomputed coefficients with tensor products
682 : // to get d2(beta) / dS2
683 : Real d2beta_dS2_times_distortion_weight = (
684 19891948260 : d2beta_dS2_coef_ia_applied * I(j,b) +
685 18360615536 : d2beta_dS2_coefs_ij_applied[0] * S(a,b) +
686 18360615536 : d2beta_dS2_coefs_ij_applied[1] * S_inv(b,a) +
687 19891948260 : d2beta_dS2_coefs_ij_applied[2] * S(a,b) +
688 18360615536 : d2beta_dS2_coefs_ij_applied[3] * S_inv(b,a) +
689 18360615536 : d2beta_dS2_coef_ja_applied * S_inv(b,i)
690 18360615536 : );
691 :
692 : // Incorporate tensor product portion to get d2(mu) /
693 : // dS2
694 18360615536 : const Real d2mu_dS2_times_dilation_weight = dalpha_dS_coef_ij_applied * S_inv(b,a) - alpha_ja_applied * S_inv(b,i);
695 :
696 :
697 : // Chain rule to change d/dS to d/dR
698 18360615536 : const auto d2E_dS2 =
699 : d2beta_dS2_times_distortion_weight +
700 : d2mu_dS2_times_dilation_weight;
701 :
702 : // if !(a == var_id1 (next line) && i == var_id2
703 : // (outside 'a' loop)), dS_dR_l(p) multiplier is zero
704 18360615536 : d2E_dSdR_l += d2E_dS2 * dS_dR_l(a, b);
705 :
706 18360615536 : if (!(i == a && j == b))
707 : // if !(a == var_id2 (next line) && i == var_id1
708 : // (outside 'a' loop)), dS_dR_p(l) multiplier is zero
709 16063512640 : d2E_dSdR_p += d2E_dS2 * dS_dR_p(a, b);
710 :
711 : } // for b
712 : } // for a
713 20655915952 : d2E_dR2 +=
714 20655915952 : d2E_dSdR_l * dS_dR_p(i, j) + d2E_dSdR_p * dS_dR_l(i, j);
715 : }// for j
716 : }// for i, end tensor contraction
717 :
718 : // Jacobian contribution
719 2297102896 : const Real jacobian_contribution = quad_weights[qp] * d2E_dR2;
720 2297102896 : K[var_id1][var_id2](l, p) += jacobian_contribution;
721 : // Jacobian is symmetric, add contribution to p,l entry
722 : // Don't get the diagonal twice!
723 2297102896 : if (p < l)
724 : // Note the transposition of var_id1 and var_id2 as these are also jacobian indices
725 2167154682 : K[var_id2][var_id1](p, l) += jacobian_contribution;
726 :
727 : }// for var_id2
728 : }// for p
729 : }// for var_id1
730 : }// for l
731 : }
732 8046296 : } // end of the quadrature point qp-loop
733 :
734 884650 : return request_jacobian;
735 8046296 : }
736 :
737 3692 : const MeshQualityInfo & VariationalSmootherSystem::get_mesh_info()
738 : {
739 3692 : if (!_mesh_info.initialized)
740 3692 : compute_mesh_quality_info();
741 :
742 3692 : return _mesh_info;
743 : }
744 :
745 144657 : void VariationalSmootherSystem::compute_mesh_quality_info()
746 : {
747 : // If the reference volume has not yet been computed, compute it.
748 148735 : if (std::abs(_ref_vol) < TOLERANCE * TOLERANCE)
749 0 : prepare_for_smoothing();
750 :
751 148767 : std::unique_ptr<DiffContext> con = this->build_context();
752 4110 : FEMContext & femcontext = cast_ref<FEMContext &>(*con);
753 144657 : this->init_context(femcontext);
754 :
755 8188 : const auto & mesh = this->get_mesh();
756 144657 : const auto dim = mesh.mesh_dimension();
757 144657 : const Real half_dim = 0.5 * dim;
758 144657 : const auto distortion_weight = 1. - _dilation_weight;
759 :
760 : // Make pre-requests before reinit() for efficiency in
761 : // --enable-deprecated builds, and to avoid errors in
762 : // --disable-deprecated builds.
763 4110 : const auto & fe_map = femcontext.get_element_fe(0)->get_fe_map();
764 4110 : const auto & JxW = fe_map.get_JxW();
765 4110 : fe_map.get_dxyzdxi();
766 4110 : fe_map.get_dxyzdeta();
767 4110 : fe_map.get_dxyzdzeta();
768 :
769 144657 : MeshQualityInfo info;
770 :
771 1823366 : for (const auto * elem : mesh.active_local_element_ptr_range())
772 : {
773 838742 : femcontext.pre_fe_reinit(*this, elem);
774 838742 : femcontext.elem_fe_reinit();
775 :
776 : // Element-integrated quantities
777 69152 : Real det_S_int = 0.;
778 69152 : Real beta_int = 0.;
779 69152 : Real mu_int = 0.;
780 69152 : Real combined_int = 0.;
781 :
782 10744966 : for (const auto qp : index_range(JxW))
783 : {
784 10734316 : const auto det_SxW = JxW[qp] * _target_jacobian_dets[elem->type()][qp];
785 9906224 : det_S_int += det_SxW;
786 :
787 : // Grab the physical-to-reference mapping Jacobian matrix (i.e., "S") at this qp
788 9906224 : RealTensor S = get_jacobian_at_qp(fe_map, dim, qp);
789 :
790 : // Apply target element transformation
791 9906224 : S *= _target_jacobians[elem->type()][qp];
792 :
793 : // Determinant of S
794 9906224 : const auto det = S.det();
795 9906224 : const auto det_sq = det * det;
796 :
797 9906224 : if (det > info.max_qp_det_S)
798 543727 : info.max_qp_det_S = det;
799 9362497 : else if (det < info.min_qp_det_S)
800 567507 : info.min_qp_det_S = det;
801 :
802 9906224 : if (det < TOLERANCE * TOLERANCE)
803 1356 : info.mesh_is_tangled = true;
804 :
805 : // trace of S^T * S
806 9906224 : const auto tr = trace(S.transpose() * S, dim);
807 :
808 : // The chi function allows us to handle degenerate elements
809 9906224 : const auto chi = chi_epsilon(det, _epsilon_squared_assembly);
810 :
811 : // distortion
812 9906224 : const Real beta = std::pow(tr / dim, half_dim) / chi;
813 9906224 : beta_int += beta * det_SxW;
814 :
815 : // dilation
816 9906224 : const Real mu = 0.5 * (_ref_vol + det_sq / _ref_vol) / chi;
817 9906224 : mu_int += mu * det_SxW;
818 :
819 : // combined
820 9906224 : const Real E = distortion_weight * beta + _dilation_weight * mu;
821 9906224 : combined_int += E * det_SxW;
822 : }
823 :
824 838742 : info.total_det_S += det_S_int;
825 838742 : if (det_S_int > info.max_elem_det_S.first)
826 20818 : info.max_elem_det_S = std::make_pair(det_S_int, elem->id());
827 625838 : else if (det_S_int < info.min_elem_det_S.first)
828 16879 : info.min_elem_det_S = std::make_pair(det_S_int, elem->id());
829 :
830 838742 : info.total_distortion += beta_int;
831 838742 : if (beta_int > info.max_elem_distortion.first)
832 20662 : info.max_elem_distortion = std::make_pair(beta_int, elem->id());
833 626793 : else if (beta_int < info.min_elem_distortion.first)
834 17208 : info.min_elem_distortion = std::make_pair(beta_int, elem->id());
835 :
836 838742 : info.total_dilation += mu_int;
837 838742 : if (mu_int > info.max_elem_dilation.first)
838 20793 : info.max_elem_dilation = std::make_pair(mu_int, elem->id());
839 626068 : else if (mu_int < info.min_elem_dilation.first)
840 17063 : info.min_elem_dilation = std::make_pair(mu_int, elem->id());
841 :
842 838742 : info.total_combined += combined_int;
843 838742 : if (combined_int > info.max_elem_combined.first)
844 20698 : info.max_elem_combined = std::make_pair(combined_int, elem->id());
845 626492 : else if (combined_int < info.min_elem_combined.first)
846 17104 : info.min_elem_combined = std::make_pair(combined_int, elem->id());
847 :
848 136469 : } // for elem
849 :
850 : // Get contributions from elements on other processors
851 144657 : communicate_pair_max(info.max_elem_det_S, mesh.comm());
852 144657 : communicate_pair_min(info.min_elem_det_S, mesh.comm());
853 144657 : mesh.comm().max(info.max_qp_det_S);
854 144657 : mesh.comm().min(info.min_qp_det_S);
855 144657 : mesh.comm().sum(info.total_det_S);
856 :
857 144657 : communicate_pair_max(info.max_elem_distortion, mesh.comm());
858 144657 : communicate_pair_min(info.min_elem_distortion, mesh.comm());
859 144657 : mesh.comm().sum(info.total_distortion);
860 :
861 144657 : communicate_pair_max(info.max_elem_dilation, mesh.comm());
862 144657 : communicate_pair_min(info.min_elem_dilation, mesh.comm());
863 144657 : mesh.comm().sum(info.total_dilation);
864 :
865 144657 : communicate_pair_max(info.max_elem_combined, mesh.comm());
866 144657 : communicate_pair_min(info.min_elem_combined, mesh.comm());
867 144657 : mesh.comm().sum(info.total_combined);
868 :
869 144657 : mesh.comm().max(info.mesh_is_tangled);
870 :
871 144657 : _mesh_info = info;
872 144657 : }
873 :
874 : std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
875 1314 : VariationalSmootherSystem::get_target_elem(const ElemType & type)
876 : {
877 : // Build target element
878 1364 : auto target_elem = Elem::build(type);
879 :
880 : // Volume of reference element
881 1314 : const auto ref_vol = target_elem->reference_elem()->volume();
882 :
883 : // Update the nodes of the target element, depending on type
884 50 : const Real sqrt_3 = std::sqrt(Real(3));
885 150 : std::vector<std::unique_ptr<Node>> owned_nodes;
886 :
887 1364 : const auto type_str = Utility::enum_to_string(type);
888 :
889 : // Elems deriving from Tri
890 1314 : if (type_str.compare(0, 3, "TRI") == 0)
891 : {
892 :
893 : // The target element will be an equilateral triangle with area equal to
894 : // the area of the reference element.
895 :
896 : // Equilateral triangle side length preserving area of the reference element
897 193 : const auto side_length = std::sqrt(4. / sqrt_3 * ref_vol);
898 :
899 : // Define the nodal locations of the vertices
900 6 : const auto & s = side_length;
901 : // x y node_id
902 199 : owned_nodes.emplace_back(Node::build(Point(0., 0.), 0));
903 193 : owned_nodes.emplace_back(Node::build(Point(s, 0.), 1));
904 199 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s), 2));
905 :
906 193 : switch (type)
907 : {
908 4 : case TRI3: {
909 : // Nothing to do here, vertices already added above
910 4 : break;
911 : }
912 :
913 55 : case TRI6: {
914 : // Define the midpoint nodes of the equilateral triangle
915 : // x y node_id
916 55 : owned_nodes.emplace_back(Node::build(Point(0.50 * s, 0.00), 3));
917 57 : owned_nodes.emplace_back(Node::build(Point(0.75 * s, 0.25 * sqrt_3 * s), 4));
918 57 : owned_nodes.emplace_back(Node::build(Point(0.25 * s, 0.25 * sqrt_3 * s), 5));
919 :
920 55 : break;
921 : }
922 :
923 0 : default:
924 0 : libmesh_error_msg("Unsupported triangular element: " << type_str);
925 : break;
926 : }
927 : } // if Tri
928 :
929 : // Elems deriving from Prism
930 1121 : else if (type_str.compare(0, 5, "PRISM") == 0)
931 : {
932 :
933 : // The target element will be a prism with an equilateral triangular
934 : // base with volume equal to the volume of the reference element.
935 :
936 : // For an equilateral triangular base with side length s, the
937 : // base area is s^2 * sqrt(3) / 4.
938 : // The prism height that will result in equal face areas is
939 : // s * sqrt(3) / 4. We choose s such that the target element has
940 : // the same volume as the reference element:
941 : // v = (s^2 * sqrt(3) / 4) * (s * sqrt(3) / 4) = 3 * s^3 / 4
942 : // --> s = (16 * v / 3)^(1/3)
943 : // I have no particular motivation for imposing equal face areas,
944 : // so this can be updated if a more `optimal` target prism is
945 : // identified.
946 :
947 : // Side length that preserves the volume of the reference element
948 278 : const auto side_length = std::pow(16. * ref_vol / 3., 1. / 3.);
949 : // Prism height with the property that all faces have equal area
950 278 : const auto target_height = 0.25 * side_length * sqrt_3;
951 :
952 12 : const auto & s = side_length;
953 12 : const auto & h = target_height;
954 : // x y z node_id
955 290 : owned_nodes.emplace_back(Node::build(Point(0., 0., 0.), 0));
956 278 : owned_nodes.emplace_back(Node::build(Point(s, 0., 0.), 1));
957 302 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
958 290 : owned_nodes.emplace_back(Node::build(Point(0., 0., h), 3));
959 290 : owned_nodes.emplace_back(Node::build(Point(s, 0., h), 4));
960 278 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, h), 5));
961 :
962 278 : if (type == PRISM15 || type == PRISM18 || type == PRISM20 || type == PRISM21)
963 : {
964 : // Define the edge midpoint nodes of the prism
965 10 : const auto & on = owned_nodes;
966 217 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 6));
967 217 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 7));
968 217 : owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 8));
969 217 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 9));
970 217 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
971 217 : owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[5]) / 2.), 11));
972 217 : owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
973 217 : owned_nodes.emplace_back(Node::build(Point((*on[4] + *on[5]) / 2.), 13));
974 217 : owned_nodes.emplace_back(Node::build(Point((*on[5] + *on[3]) / 2.), 14));
975 :
976 207 : if (type == PRISM18 || type == PRISM20 || type == PRISM21)
977 : {
978 : // Define the rectangular face midpoint nodes of the prism
979 145 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3] + *on[4]) / 4.), 15));
980 145 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4] + *on[5]) / 4.), 16));
981 145 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3] + *on[5]) / 4.), 17));
982 :
983 137 : if (type == PRISM20 || type == PRISM21)
984 : {
985 : // Define the triangular face midpoint nodes of the prism
986 72 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 18));
987 72 : owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4] + *on[5]) / 3.), 19));
988 :
989 66 : if (type == PRISM21)
990 : // Define the interior point of the prism
991 52 : owned_nodes.emplace_back(Node::build(Point((*on[9] + *on[10] + *on[11]) / 3.), 20));
992 :
993 : }
994 10 : }
995 : }
996 :
997 71 : else if (type != PRISM6)
998 0 : libmesh_error_msg("Unsupported prism element: " << type_str);
999 :
1000 : } // if Prism
1001 :
1002 : // Elems deriving from Pyramid
1003 843 : else if (type_str.compare(0, 7, "PYRAMID") == 0)
1004 : {
1005 :
1006 : // The target element is a pyramid with an square base and
1007 : // equilateral triangular sides with volume equal to the volume of the
1008 : // reference element.
1009 :
1010 : // A pyramid with square base sidelength s and equilateral triangular
1011 : // sides has height h = s / sqrt(2).
1012 : // The volume is v = s^2 h / 3 = s^3 / ( 3 sqrt(2)).
1013 : // Solving for s: s = (3 sqrt(2) v)^(1/3), where v is the volume of the
1014 : // non-optimal reference element.
1015 :
1016 10 : const Real sqrt_2 = std::sqrt(Real(2));
1017 : // Side length that preserves the volume of the reference element
1018 305 : const auto side_length = std::pow(3. * sqrt_2 * ref_vol, 1. / 3.);
1019 : // Pyramid height with the property that all faces are equilateral triangles
1020 305 : const auto target_height = side_length / sqrt_2;
1021 :
1022 10 : const auto & s = side_length;
1023 10 : const auto & h = target_height;
1024 :
1025 : // x y z node_id
1026 315 : owned_nodes.emplace_back(Node::build(Point(0., 0., 0.), 0));
1027 315 : owned_nodes.emplace_back(Node::build(Point(s, 0., 0.), 1));
1028 315 : owned_nodes.emplace_back(Node::build(Point(s, s, 0.), 2));
1029 305 : owned_nodes.emplace_back(Node::build(Point(0., s, 0.), 3));
1030 315 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * s, h), 4));
1031 :
1032 305 : if (type == PYRAMID13 || type == PYRAMID14 || type == PYRAMID18)
1033 : {
1034 8 : const auto & on = owned_nodes;
1035 : // Define the edge midpoint nodes of the pyramid
1036 :
1037 : // Base node to base node midpoint nodes
1038 242 : owned_nodes.emplace_back(Node::build((*on[0] + *on[1]) / 2., 5));
1039 242 : owned_nodes.emplace_back(Node::build((*on[1] + *on[2]) / 2., 6));
1040 242 : owned_nodes.emplace_back(Node::build((*on[2] + *on[3]) / 2., 7));
1041 242 : owned_nodes.emplace_back(Node::build((*on[3] + *on[0]) / 2., 8));
1042 :
1043 : // Base node to apex node midpoint nodes
1044 242 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[4]) / 2.), 9));
1045 242 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
1046 242 : owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[4]) / 2.), 11));
1047 242 : owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
1048 :
1049 234 : if (type == PYRAMID14 || type == PYRAMID18)
1050 : {
1051 : // Define the square face midpoint node of the pyramid
1052 151 : owned_nodes.emplace_back(
1053 169 : Node::build(Point((*on[0] + *on[1] + *on[2] + *on[3]) / 4.), 13));
1054 :
1055 163 : if (type == PYRAMID18)
1056 : {
1057 : // Define the triangular face nodes
1058 96 : owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[4]) / 3.), 14));
1059 96 : owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4]) / 3.), 15));
1060 96 : owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3] + *on[4]) / 3.), 16));
1061 100 : owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[0] + *on[4]) / 3.), 17));
1062 : }
1063 8 : }
1064 : }
1065 :
1066 71 : else if (type != PYRAMID5)
1067 0 : libmesh_error_msg("Unsupported pyramid element: " << type_str);
1068 :
1069 : } // if Pyramid
1070 :
1071 : // Set the target_elem equal to the reference elem
1072 : else
1073 7557 : for (const auto & node : target_elem->reference_elem()->node_ref_range())
1074 7513 : owned_nodes.emplace_back(Node::build(node, node.id()));
1075 :
1076 : // Set nodes of target element
1077 17101 : for (const auto & node_ptr : owned_nodes)
1078 16407 : target_elem->set_node(node_ptr->id(), node_ptr.get());
1079 :
1080 50 : libmesh_assert(relative_fuzzy_equals(target_elem->volume(), ref_vol, TOLERANCE));
1081 :
1082 1364 : return std::make_pair(std::move(target_elem), std::move(owned_nodes));
1083 1214 : }
1084 :
1085 1314 : void VariationalSmootherSystem::get_target_to_reference_jacobian(
1086 : const Elem * const target_elem,
1087 : const FEMContext & femcontext,
1088 : std::vector<RealTensor> & jacobians,
1089 : std::vector<Real> & jacobian_dets)
1090 : {
1091 :
1092 1314 : const auto dim = target_elem->dim();
1093 :
1094 50 : const auto & qrule_points = femcontext.get_element_qrule().get_points();
1095 50 : const auto & qrule_weights = femcontext.get_element_qrule().get_weights();
1096 50 : const auto nq_points = femcontext.get_element_qrule().n_points();
1097 :
1098 : // If the target element is the reference element, Jacobian matrix is
1099 : // identity, det of inverse is 1. These will only be overwritten if a
1100 : // different target element is explicitly specified.
1101 1364 : jacobians = std::vector<RealTensor>(nq_points, RealTensor(
1102 : 1., 0., 0.,
1103 : 0., 1., 0.,
1104 50 : 0., 0., 1.));
1105 1314 : jacobian_dets = std::vector<Real>(nq_points, 1.0);
1106 :
1107 : // Don't use "if (*target_elem == *(target_elem->reference_elem()))" here, it
1108 : // only compares global node ids, not the node locations themselves.
1109 50 : bool target_equals_reference = true;
1110 1314 : const auto * ref_elem = target_elem->reference_elem();
1111 17101 : for (const auto local_id : make_range(target_elem->n_nodes()))
1112 15787 : target_equals_reference &= target_elem->node_ref(local_id) == ref_elem->node_ref(local_id);
1113 1314 : if (target_equals_reference)
1114 538 : return;
1115 :
1116 : // Create FEMap to compute target_element mapping information
1117 832 : FEMap fe_map_target;
1118 :
1119 : // pre-request mapping derivatives
1120 28 : fe_map_target.get_dxyzdxi();
1121 28 : fe_map_target.get_dxyzdeta();
1122 28 : fe_map_target.get_dxyzdzeta();
1123 :
1124 : // build map
1125 776 : fe_map_target.init_reference_to_physical_map(dim, qrule_points, target_elem);
1126 776 : fe_map_target.compute_map(dim, qrule_weights, target_elem, /*d2phi=*/false);
1127 :
1128 18700 : for (const auto qp : make_range(nq_points))
1129 : {
1130 : // We use Larisa's H notation to denote the reference-to-target jacobian
1131 17924 : RealTensor H = get_jacobian_at_qp(fe_map_target, dim, qp);
1132 :
1133 : // The target-to-reference jacobian is the inverse of the
1134 : // reference-to-target jacobian
1135 17924 : jacobians[qp] = H.inverse();
1136 18722 : jacobian_dets[qp] = jacobians[qp].det();
1137 : }
1138 720 : }
1139 :
1140 : } // namespace libMesh
|
