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 1608 : VariationalSmootherSystem::~VariationalSmootherSystem () = default;
42 :
43 16898 : void VariationalSmootherSystem::assembly (bool get_residual,
44 : bool get_jacobian,
45 : bool apply_heterogeneous_constraints,
46 : bool apply_no_constraints)
47 : {
48 : // Update the mesh based on the current variable values
49 952 : auto & mesh = this->get_mesh();
50 16898 : this->solution->close();
51 :
52 1349756 : for (auto * node : mesh.local_node_ptr_range())
53 : {
54 2750160 : for (const auto d : make_range(mesh.mesh_dimension()))
55 : {
56 2032104 : const auto dof_id = node->dof_number(this->number(), d, 0);
57 : // Update mesh
58 2032104 : (*node)(d) = libmesh_real((*current_local_solution)(dof_id));
59 : }
60 15946 : }
61 :
62 16898 : SyncNodalPositions sync_object(mesh);
63 33320 : Parallel::sync_dofobject_data_by_id (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);
64 :
65 16898 : FEMSystem::assembly(get_residual, get_jacobian, apply_heterogeneous_constraints, apply_no_constraints);
66 16898 : }
67 :
68 852 : void VariationalSmootherSystem::init_data ()
69 : {
70 48 : auto & mesh = this->get_mesh();
71 24 : const auto & elem_orders = mesh.elem_default_orders();
72 852 : libmesh_error_msg_if(elem_orders.size() != 1,
73 : "The variational smoother cannot be used for mixed-order meshes!");
74 852 : const auto fe_order = *elem_orders.begin();
75 : // Add a variable for each dimension of the mesh
76 : // "r0" for x, "r1" for y, "r2" for z
77 2627 : for (const auto & d : make_range(mesh.mesh_dimension()))
78 3500 : this->add_variable ("r" + std::to_string(d), fe_order);
79 :
80 : // Do the parent's initialization after variables are defined
81 852 : FEMSystem::init_data();
82 :
83 : // Set the current_local_solution to the current mesh for the initial guess
84 852 : this->solution->close();
85 :
86 88074 : for (auto * node : mesh.local_node_ptr_range())
87 : {
88 184644 : for (const auto d : make_range(mesh.mesh_dimension()))
89 : {
90 137520 : const auto dof_id = node->dof_number(this->number(), d, 0);
91 : // Update solution
92 137520 : (*solution).set(dof_id, (*node)(d));
93 : }
94 804 : }
95 :
96 852 : this->prepare_for_smoothing();
97 852 : }
98 :
99 852 : void VariationalSmootherSystem::prepare_for_smoothing()
100 : {
101 876 : std::unique_ptr<DiffContext> con = this->build_context();
102 24 : FEMContext & femcontext = cast_ref<FEMContext &>(*con);
103 852 : this->init_context(femcontext);
104 :
105 48 : const auto & mesh = this->get_mesh();
106 :
107 852 : Real elem_averaged_det_J_sum = 0.;
108 :
109 : // Make pre-requests before reinit() for efficiency in
110 : // --enable-deprecated builds, and to avoid errors in
111 : // --disable-deprecated builds.
112 24 : const auto & fe_map = femcontext.get_element_fe(0)->get_fe_map();
113 24 : const auto & JxW = fe_map.get_JxW();
114 :
115 24 : std::map<ElemType, std::vector<Real>> target_jacobians_dets;
116 :
117 10260 : for (const auto * elem : mesh.active_local_element_ptr_range())
118 : {
119 8580 : femcontext.pre_fe_reinit(*this, elem);
120 8580 : femcontext.elem_fe_reinit();
121 :
122 : // Add target element info, if applicable
123 9295 : if (_target_jacobians.find(elem->type()) == _target_jacobians.end())
124 : {
125 758 : const auto [target_elem, target_nodes] = get_target_elem(elem->type());
126 758 : get_target_to_reference_jacobian(target_elem.get(),
127 : femcontext,
128 758 : _target_jacobians[elem->type()],
129 1516 : target_jacobians_dets[elem->type()]);
130 : }// if find == end()
131 :
132 : // Reference volume computation
133 715 : Real elem_integrated_det_J = 0.;
134 153960 : for (const auto qp : index_range(JxW))
135 157495 : elem_integrated_det_J += JxW[qp] * target_jacobians_dets[elem->type()][qp];
136 8580 : const auto ref_elem_vol = elem->reference_elem()->volume();
137 8580 : elem_averaged_det_J_sum += elem_integrated_det_J / ref_elem_vol;
138 :
139 804 : } // for elem
140 :
141 : // Get contributions from elements on other processors
142 852 : mesh.comm().sum(elem_averaged_det_J_sum);
143 :
144 852 : _ref_vol = elem_averaged_det_J_sum / mesh.n_active_elem();
145 852 : }
146 :
147 19178 : void VariationalSmootherSystem::init_context(DiffContext & context)
148 : {
149 976 : FEMContext & c = cast_ref<FEMContext &>(context);
150 :
151 976 : FEBase * my_fe = nullptr;
152 :
153 : // Now make sure we have requested all the data
154 : // we need to build the system.
155 :
156 : // We might have a multi-dimensional mesh
157 : const std::set<unsigned char> & elem_dims =
158 976 : c.elem_dimensions();
159 :
160 38356 : for (const auto & dim : elem_dims)
161 : {
162 19178 : c.get_element_fe( 0, my_fe, dim );
163 976 : my_fe->get_nothing();
164 :
165 976 : auto & fe_map = my_fe->get_fe_map();
166 976 : fe_map.get_dxyzdxi();
167 976 : fe_map.get_dxyzdeta();
168 976 : fe_map.get_dxyzdzeta();
169 976 : fe_map.get_JxW();
170 :
171 19178 : c.get_side_fe( 0, my_fe, dim );
172 976 : my_fe->get_nothing();
173 : }
174 :
175 19178 : FEMSystem::init_context(context);
176 19178 : }
177 :
178 :
179 164160 : bool VariationalSmootherSystem::element_time_derivative (bool request_jacobian,
180 : DiffContext & context)
181 : {
182 13680 : FEMContext & c = cast_ref<FEMContext &>(context);
183 :
184 27360 : const Elem & elem = c.get_elem();
185 :
186 : // First we get some references to cell-specific data that
187 : // will be used to assemble the linear system.
188 :
189 13680 : const auto & fe_map = c.get_element_fe(0)->get_fe_map();
190 :
191 13680 : const auto & dxyzdxi = fe_map.get_dxyzdxi();
192 13680 : const auto & dxyzdeta = fe_map.get_dxyzdeta();
193 13680 : const auto & dxyzdzeta = fe_map.get_dxyzdzeta();
194 :
195 13680 : const auto & dphidxi_map = fe_map.get_dphidxi_map();
196 13680 : const auto & dphideta_map = fe_map.get_dphideta_map();
197 13680 : const auto & dphidzeta_map = fe_map.get_dphidzeta_map();
198 :
199 164160 : unsigned int dim = c.get_dim();
200 :
201 13680 : unsigned int x_var = 0, y_var = 1, z_var = 2;
202 : // In the case of lower dimensions, we will not access z (1D/2D) or y (1D)
203 164160 : if (dim < 3)
204 : {
205 7180 : z_var = 0;
206 86160 : if (dim < 2)
207 180 : y_var = 0;
208 : }
209 :
210 : // The subvectors and submatrices we need to fill:
211 : // system residual
212 : std::reference_wrapper<DenseSubVector<Number>> F[3] =
213 : {
214 : c.get_elem_residual(x_var),
215 : c.get_elem_residual(y_var),
216 : c.get_elem_residual(z_var)
217 13680 : };
218 : // system jacobian
219 : std::reference_wrapper<DenseSubMatrix<Number>> K[3][3] =
220 : {
221 : {c.get_elem_jacobian(x_var, x_var), c.get_elem_jacobian(x_var, y_var), c.get_elem_jacobian(x_var, z_var)},
222 : {c.get_elem_jacobian(y_var, x_var), c.get_elem_jacobian(y_var, y_var), c.get_elem_jacobian(y_var, z_var)},
223 : {c.get_elem_jacobian(z_var, x_var), c.get_elem_jacobian(z_var, y_var), c.get_elem_jacobian(z_var, z_var)}
224 13680 : };
225 :
226 : // Quadrature info
227 164160 : const auto quad_weights = c.get_element_qrule().get_weights();
228 :
229 164160 : const auto distortion_weight = 1. - _dilation_weight;
230 :
231 : // Integrate the distortion-dilation metric over the reference element
232 2413440 : for (const auto qp : index_range(quad_weights))
233 : {
234 : // Compute quantities needed to evaluate the distortion-dilation metric
235 : // and its gradient and Hessian.
236 : // The metric will be minimized when it's gradient with respect to the node
237 : // locations (R) is zero. For Newton's method, minimizing the gradient
238 : // requires computation of the gradient's Jacobian with respect to R.
239 : // The Jacobian of the distortion-dilation metric's gradientis the Hessian
240 : // of the metric.
241 : //
242 : // Note that the term "Jacobian" has two meanings in this mesh smoothing
243 : // application. The first meaning refers to the Jacobian w.r.t R of the
244 : // gradient w.r.t. R, (i.e., the Hessian of the metric). This is the K
245 : // variable defined above. This is also the Jacobian the 'request_jacobian'
246 : // variable refers to. The second mesning refers to the Jacobian of the
247 : // physical-to-reference element mapping. This is the Jacobian used to
248 : // compute the distorion-dilation metric.
249 : //
250 : // Grab the physical-to-reference mapping Jacobian matrix (i.e., "S") at this qp
251 374880 : RealTensor S;
252 : // RealTensors are always 3x3, so we will fill any dimensions above dim
253 : // with 1s on the diagonal. This indicates a 1 to 1 relationship between
254 : // the physical and reference elements in these extra dimensions.
255 2249280 : switch (dim)
256 : {
257 5880 : case 1:
258 5880 : S = RealTensor(
259 980 : dxyzdxi[qp](0), 0, 0,
260 : 0, 1, 0,
261 : 0, 0, 1
262 : );
263 5880 : break;
264 :
265 422400 : case 2:
266 528000 : S = RealTensor(
267 70400 : dxyzdxi[qp](0), dxyzdeta[qp](0), 0,
268 140800 : dxyzdxi[qp](1), dxyzdeta[qp](1), 0,
269 : 0, 0, 1
270 : );
271 422400 : break;
272 :
273 303500 : case 3:
274 1972750 : S = RealTensor(
275 303500 : dxyzdxi[qp],
276 303500 : dxyzdeta[qp],
277 303500 : dxyzdzeta[qp]
278 303500 : ).transpose(); // Note the transposition!
279 1821000 : break;
280 :
281 0 : default:
282 0 : libmesh_error_msg("Unsupported dimension.");
283 : }
284 :
285 : // Apply any applicable target element transformations
286 2436720 : if (_target_jacobians.find(elem.type()) != _target_jacobians.end())
287 2249280 : S = S * _target_jacobians[elem.type()][qp];
288 :
289 : // Compute quantities needed for the smoothing algorithm
290 :
291 : // determinant
292 2249280 : const Real det = S.det();
293 2249280 : const Real det_sq = det * det;
294 2249280 : const Real det_cube = det_sq * det;
295 :
296 2249280 : const Real ref_vol_sq = _ref_vol * _ref_vol;
297 :
298 : // trace of S^T * S
299 : // DO NOT USE RealTensor.tr for the trace, it will NOT be correct for
300 : // 1D and 2D meshes because of our hack of putting 1s in the diagonal of
301 : // S for the extra dimensions
302 2624160 : const RealTensor ST_S = S.transpose() * S;
303 187440 : Real tr = 0.0;
304 8562960 : for (const auto i : make_range(dim))
305 6313680 : tr += ST_S(i, i);
306 2249280 : const Real tr_div_dim = tr / dim;
307 :
308 : // Precompute pow(tr_div_dim, 0.5 * dim - x) for x = 0, 1, 2
309 2249280 : const Real half_dim = 0.5 * dim;
310 : const std::vector<Real> trace_powers{
311 2249280 : std::pow(tr_div_dim, half_dim),
312 2249280 : std::pow(tr_div_dim, half_dim - 1.),
313 2249280 : std::pow(tr_div_dim, half_dim - 2.),
314 2436720 : };
315 :
316 : // inverse of S
317 2436720 : const RealTensor S_inv = S.inverse();
318 : // inverse transpose of S
319 374880 : const RealTensor S_inv_T = S_inv.transpose();
320 :
321 : // Identity matrix
322 3748800 : const RealTensor I(1, 0, 0,
323 3748800 : 0, 1, 0,
324 2249280 : 0, 0, 1);
325 :
326 : // The chi function allows us to handle degenerate elements
327 : // When the element is not degenerate (i.e., det(S) > 0), we can set
328 : // _epsilon_squared to zero to get chi(det(S)) = det(S)
329 2249280 : const Real epsilon_sq = det > 0. ? 0. : _epsilon_squared;
330 2249280 : const Real chi = 0.5 * (det + std::sqrt(epsilon_sq + det_sq));
331 2249280 : const Real chi_sq = chi * chi;
332 2249280 : const Real sqrt_term = std::sqrt(epsilon_sq + det_sq);
333 : // dchi(x) / dx
334 2249280 : const Real chi_prime = 0.5 * (1. + det / sqrt_term);
335 2249280 : const Real chi_prime_sq = chi_prime * chi_prime;
336 : // d2chi(x) / dx2
337 2249280 : const Real chi_2prime = 0.5 * (1. / sqrt_term - det_sq / Utility::pow<3>(sqrt_term));
338 :
339 : // Distortion metric (beta)
340 : //const Real beta = trace_powers[0] / chi;
341 2999040 : const RealTensor dbeta_dS = (trace_powers[1] / chi) * S - (trace_powers[0] / chi_sq * chi_prime * det) * S_inv_T;
342 :
343 : // Dilation metric (mu)
344 : //const Real mu = 0.5 * (_ref_vol + det_sq / _ref_vol) / chi;
345 : // We represent d mu / dS as alpha(S) * S^-T, where alpha is a scalar function
346 2249280 : const Real alpha = (-chi_prime * det_cube + 2. * det_sq * chi - ref_vol_sq * det * chi_prime) / (2. * _ref_vol * chi_sq);
347 2249280 : const RealTensor dmu_dS = alpha * S_inv_T;
348 :
349 : // Combined metric (E)
350 : //const Real E = distortion_weight * beta + _dilation_weight * mu;
351 2811600 : const RealTensor dE_dS = distortion_weight * dbeta_dS + _dilation_weight * dmu_dS;
352 :
353 : // This vector is useful in computing dS/dR below
354 11621280 : std::vector<std::vector<std::vector<Real>>> dphi_maps = {dphidxi_map, dphideta_map, dphidzeta_map};
355 :
356 : // Compute residual (i.e., the gradient of the combined metric w.r.t node locations)
357 : // Recall that when the gradient (residual) is zero, the combined metric is minimized
358 47079720 : for (const auto l : elem.node_index_range())
359 : {
360 177392880 : for (const auto var_id : make_range(dim))
361 : {
362 : // Build dS/dR, the derivative of the physical-to-reference mapping Jacobin w.r.t. the mesh node locations
363 121515570 : RealTensor dS_dR = RealTensor(0);
364 526424880 : for (const auto jj : make_range(dim))
365 525149920 : dS_dR(var_id, jj) = dphi_maps[jj][l][qp];
366 :
367 : // Residual contribution. The contraction of dE/dS and dS/dR gives us
368 : // the gradient we are looking for, dE/dR
369 276171750 : F[var_id](l) += quad_weights[qp] * dE_dS.contract(dS_dR);
370 : }// for var_id
371 : }// for l
372 :
373 :
374 2249280 : if (request_jacobian)
375 : {
376 : // Compute jacobian of the smoothing system (i.e., the Hessian of the
377 : // combined metric w.r.t. the mesh node locations)
378 :
379 : // Precompute coefficients to be applied to each component of tensor
380 : // products in the loops below. At first glance, these coefficients look
381 : // like gibberish, but everything in the Hessian has been verified by
382 : // taking finite differences of the gradient. We should probably write
383 : // down the derivations of the gradient and Hessian somewhere...
384 :
385 : // Recall that above, dbeta_dS takes the form:
386 : // d(beta)/dS = c1(S) * S - c2(S) * S_inv_T,
387 : // where c1 and c2 are scalar-valued functions.
388 : const std::vector<Real> d2beta_dS2_coefs_times_distortion_weight = {
389 : //Part 1: scaler coefficients of d(c1 * S) / dS
390 : //
391 : // multiplies I[i,a] x I[j,b]
392 1124640 : (trace_powers[1] / chi) * distortion_weight,
393 : // multiplies S[a,b] x S[i,j]
394 1218360 : (((dim - 2.) / dim) * trace_powers[2] / chi) * distortion_weight,
395 : // multiplies S_inv[b,a] * S[i,j]
396 1218360 : (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
397 : //
398 : //Part 2: scaler coefficients of d(-c2 * S_inv_T) / dS
399 : //
400 : // multiplies S[a,b] x S_inv[j,i]
401 1218360 : (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
402 : // multiplies S_inv[b,a] x S_inv[j,i]
403 1218360 : (trace_powers[0] * (det / chi_sq)
404 1124640 : * ((2. * chi_prime_sq / chi - chi_2prime) * det - chi_prime)) * distortion_weight,
405 : // multiplies S_inv[b,i] x S_inv[j,a]
406 1218360 : ((trace_powers[0] / chi_sq) * chi_prime * det) * distortion_weight,
407 1218360 : };
408 :
409 : // d alpha / dS has the form c(S) * S^-T, where c is the scalar coefficient defined below
410 1312080 : const Real dalpha_dS_coef_times_dilation_weight = ((det / (2. * _ref_vol * chi_sq))
411 1312080 : * (-4. * _ref_vol * alpha * chi * chi_prime - chi_2prime * det_cube
412 1312080 : - chi_prime * det_sq + 4 * chi * det - ref_vol_sq * (chi_prime + det * chi_2prime))) * _dilation_weight;
413 :
414 : // This is also useful to precompute
415 1124640 : const Real alpha_times_dilation_weight = alpha * _dilation_weight;
416 :
417 : /*
418 :
419 : To increase the efficiency of the Jacobian computation, we take
420 : advantage of l-p symmetry, ij-ab symmetry, and the sparsity pattern of
421 : dS_dR. We also factor all possible multipliers out of inner loops for
422 : efficiency. The result is code that is more difficult to read. For
423 : clarity, consult the pseudo-code below in this comment.
424 :
425 : for (const auto l: elem.node_index_range()) // Contribution to Hessian
426 : from node l
427 : {
428 : for (const auto var_id1 : make_range(dim)) // Contribution from each
429 : x/y/z component of node l
430 : {
431 : // Build dS/dR_l, the derivative of the physical-to-reference
432 : mapping Jacobin w.r.t.
433 : // the l-th node
434 : RealTensor dS_dR_l = RealTensor(0);
435 : for (const auto ii : make_range(dim))
436 : dS_dR_l(var_id1, ii) = dphi_maps[ii][l][qp];
437 :
438 : for (const auto p: elem.node_index_range()) // Contribution to
439 : Hessian from node p
440 : {
441 : for (const auto var_id2 : make_range(dim)) // Contribution from
442 : each x/y/z component of node p
443 : {
444 : // Build dS/dR_l, the derivative of the physical-to-reference
445 : mapping Jacobin w.r.t.
446 : // the p-th node
447 : RealTensor dS_dR_p = RealTensor(0);
448 : for (const auto jj : make_range(dim))
449 : dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
450 :
451 : Real d2beta_dR2 = 0.;
452 : Real d2mu_dR2 = 0.;
453 : // Perform tensor contraction
454 : for (const auto i : make_range(dim))
455 : {
456 : for (const auto j : make_range(dim))
457 : {
458 : for (const auto a : make_range(dim))
459 : {
460 : for (const auto b : make_range(dim))
461 : {
462 : // Nasty tensor products to be multiplied by
463 : d2beta_dS2_coefs to get d2(beta) / dS2 const std::vector<Real>
464 : d2beta_dS2_tensor_contributions =
465 : {
466 : I(i,a) * I(j,b),
467 : S(a,b) * S(i,j),
468 : S_inv(b,a) * S(i,j),
469 : S(a,b) * S_inv(j,i),
470 : S_inv(b,a) * S_inv(j,i),
471 : S_inv(j,a) * S_inv(b,i),
472 : };
473 :
474 : // Combine precomputed coefficients with tensor products
475 : to get d2(beta) / dS2 Real d2beta_dS2 = 0.; for (const auto comp_id :
476 : index_range(d2beta_dS2_coefs))
477 : {
478 : const Real contribution = d2beta_dS2_coefs[comp_id] *
479 : d2beta_dS2_tensor_contributions[comp_id]; d2beta_dS2 += contribution;
480 : }
481 :
482 : // Incorporate tensor product portion to get d2(mu) /
483 : dS2 const Real d2mu_dS2 = dalpha_dS_coef * S_inv(b,a) * S_inv(j,i) -
484 : alpha * S_inv(b,i) * S_inv(j,a);
485 :
486 : // Chain rule to change d/dS to d/dR
487 : d2beta_dR2 += d2beta_dS2 * dS_dR_l(a, b) * dS_dR_p(i,
488 : j); d2mu_dR2 += d2mu_dS2 * dS_dR_l(a, b) * dS_dR_p(i, j);
489 :
490 : }// for b
491 : }// for a
492 : }// for j
493 : }// for i, end tensor contraction
494 :
495 : // Jacobian contribution
496 : K[var_id1][var_id2](l, p) += quad_weights[qp] * ((1. -
497 : _dilation_weight) * d2beta_dR2 + _dilation_weight * d2mu_dR2);
498 :
499 : }// for var_id2
500 : }// for p
501 : }// for var_id1
502 : }// for l
503 :
504 : End pseudo-code, begin efficient code
505 :
506 : */
507 :
508 23539860 : for (const auto l: elem.node_index_range()) // Contribution to Hessian from node l
509 : {
510 88696440 : for (const auto var_id1 : make_range(dim)) // Contribution from each x/y/z component of node l
511 : {
512 : // Build dS/dR_l, the derivative of the physical-to-reference mapping Jacobin w.r.t.
513 : // the l-th node
514 66281220 : RealTensor dS_dR_l = RealTensor(0);
515 263212440 : for (const auto ii : make_range(dim))
516 262574960 : dS_dR_l(var_id1, ii) = dphi_maps[ii][l][qp];
517 :
518 : // Jacobian is symmetric, only need to loop over lower triangular portion
519 900014460 : for (const auto p: make_range(l + 1)) // Contribution to Hessian from node p
520 : {
521 3329246880 : for (const auto var_id2 : make_range(dim)) // Contribution from each x/y/z component of node p
522 : {
523 : // Build dS/dR_l, the derivative of the physical-to-reference mapping Jacobin w.r.t.
524 : // the p-th node
525 2495513640 : RealTensor dS_dR_p = RealTensor(0);
526 9970714080 : for (const auto jj : make_range(dim))
527 9966933920 : dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
528 :
529 207959470 : Real d2E_dR2 = 0.;
530 : // Perform tensor contraction
531 9970714080 : for (const auto i : make_range(dim))
532 : {
533 :
534 29878152480 : for (const auto j : make_range(dim))
535 : {
536 :
537 22402952040 : const auto S_ij = S(i,j);
538 22402952040 : const auto S_inv_ji = S_inv(j,i);
539 :
540 : // Apply the stuff only depending on i and j before entering a, b loops
541 : // beta
542 : const std::vector<Real> d2beta_dS2_coefs_ij_applied{
543 22402952040 : d2beta_dS2_coefs_times_distortion_weight[1] * S_ij,
544 24269864710 : d2beta_dS2_coefs_times_distortion_weight[2] * S_ij,
545 24269864710 : d2beta_dS2_coefs_times_distortion_weight[3] * S_inv_ji,
546 24269864710 : d2beta_dS2_coefs_times_distortion_weight[4] * S_inv_ji,
547 22402952040 : };
548 : // mu
549 22402952040 : const auto dalpha_dS_coef_ij_applied = dalpha_dS_coef_times_dilation_weight * S_inv_ji;
550 :
551 1866912670 : Real d2E_dSdR_l = 0.;
552 1866912670 : Real d2E_dSdR_p = 0.;
553 :
554 67186222680 : for (const auto a : make_range(i + 1))
555 : {
556 :
557 : // If this condition is met, both the ijab and abij
558 : // contributions to the Jacobian are zero due to the
559 : // spasity patterns of dS_dR_l and dS_dR_p and this
560 : // iteration may be skipped
561 44783270640 : if (!(a == var_id1 && i == var_id2) &&
562 37518384350 : !(a == var_id2 && i == var_id1))
563 34199064350 : continue;
564 :
565 7475200440 : const auto S_inv_ja = S_inv(j,a);
566 :
567 7475200440 : const Real d2beta_dS2_coef_ia_applied = d2beta_dS2_coefs_times_distortion_weight[0] * I(i,a);
568 7475200440 : const Real d2beta_dS2_coef_ja_applied = d2beta_dS2_coefs_times_distortion_weight[5] * S_inv_ja;
569 7475200440 : const Real alpha_ja_applied = alpha_times_dilation_weight * S_inv_ja;
570 :
571 7475200440 : const auto b_limit = (a == i) ? j + 1 : dim;
572 27388309080 : for (const auto b : make_range(b_limit))
573 : {
574 :
575 : // Combine precomputed coefficients with tensor products
576 : // to get d2(beta) / dS2
577 : Real d2beta_dS2_times_distortion_weight = (
578 21572534360 : d2beta_dS2_coef_ia_applied * I(j,b) +
579 19913108640 : d2beta_dS2_coefs_ij_applied[0] * S(a,b) +
580 19913108640 : d2beta_dS2_coefs_ij_applied[1] * S_inv(b,a) +
581 21572534360 : d2beta_dS2_coefs_ij_applied[2] * S(a,b) +
582 19913108640 : d2beta_dS2_coefs_ij_applied[3] * S_inv(b,a) +
583 19913108640 : d2beta_dS2_coef_ja_applied * S_inv(b,i)
584 19913108640 : );
585 :
586 : // Incorporate tensor product portion to get d2(mu) /
587 : // dS2
588 19913108640 : const Real d2mu_dS2_times_dilation_weight = dalpha_dS_coef_ij_applied * S_inv(b,a) - alpha_ja_applied * S_inv(b,i);
589 :
590 :
591 : // Chain rule to change d/dS to d/dR
592 19913108640 : const auto d2E_dS2 =
593 : d2beta_dS2_times_distortion_weight +
594 : d2mu_dS2_times_dilation_weight;
595 :
596 : // if !(a == var_id1 (next line) && i == var_id2
597 : // (outside 'a' loop)), dS_dR_l(p) multiplier is zero
598 19913108640 : d2E_dSdR_l += d2E_dS2 * dS_dR_l(a, b);
599 :
600 19913108640 : if (!(i == a && j == b))
601 : // if !(a == var_id2 (next line) && i == var_id1
602 : // (outside 'a' loop)), dS_dR_p(l) multiplier is zero
603 17417595000 : d2E_dSdR_p += d2E_dS2 * dS_dR_p(a, b);
604 :
605 : } // for b
606 : } // for a
607 22402952040 : d2E_dR2 +=
608 22402952040 : d2E_dSdR_l * dS_dR_p(i, j) + d2E_dSdR_p * dS_dR_l(i, j);
609 : }// for j
610 : }// for i, end tensor contraction
611 :
612 : // Jacobian contribution
613 2495513640 : const Real jacobian_contribution = quad_weights[qp] * d2E_dR2;
614 2495513640 : K[var_id1][var_id2](l, p) += jacobian_contribution;
615 : // Jacobian is symmetric, add contribution to p,l entry
616 : // Don't get the diagonal twice!
617 2495513640 : if (p < l)
618 : // Note the transposition of var_id1 and var_id2 as these are also jacobian indices
619 2490130955 : K[var_id2][var_id1](p, l) += jacobian_contribution;
620 :
621 : }// for var_id2
622 : }// for p
623 : }// for var_id1
624 : }// for l
625 : }
626 : } // end of the quadrature point qp-loop
627 :
628 177840 : return request_jacobian;
629 1874400 : }
630 :
631 : std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
632 758 : VariationalSmootherSystem::get_target_elem(const ElemType & type)
633 : {
634 : // Build target element
635 782 : auto target_elem = Elem::build(type);
636 :
637 : // Volume of reference element
638 758 : const auto ref_vol = target_elem->reference_elem()->volume();
639 :
640 : // Update the nodes of the target element, depending on type
641 24 : const Real sqrt_3 = std::sqrt(Real(3));
642 72 : std::vector<std::unique_ptr<Node>> owned_nodes;
643 :
644 782 : const auto type_str = Utility::enum_to_string(type);
645 :
646 : // Elems deriving from Tri
647 758 : if (type_str.compare(0, 3, "TRI") == 0)
648 : {
649 :
650 : // The target element will be an equilateral triangle with area equal to
651 : // the area of the reference element.
652 :
653 : // Equilateral triangle side length preserving area of the reference element
654 213 : const auto side_length = std::sqrt(4. / sqrt_3 * ref_vol);
655 :
656 : // Define the nodal locations of the vertices
657 6 : const auto & s = side_length;
658 : // x y node_id
659 219 : owned_nodes.emplace_back(Node::build(Point(0., 0.), 0));
660 213 : owned_nodes.emplace_back(Node::build(Point(s, 0.), 1));
661 219 : owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s), 2));
662 :
663 213 : switch (type)
664 : {
665 4 : case TRI3: {
666 : // Nothing to do here, vertices already added above
667 4 : break;
668 : }
669 :
670 71 : case TRI6: {
671 : // Define the midpoint nodes of the equilateral triangle
672 : // x y node_id
673 71 : owned_nodes.emplace_back(Node::build(Point(0.50 * s, 0.00), 3));
674 73 : owned_nodes.emplace_back(Node::build(Point(0.75 * s, 0.25 * sqrt_3 * s), 4));
675 73 : owned_nodes.emplace_back(Node::build(Point(0.25 * s, 0.25 * sqrt_3 * s), 5));
676 :
677 71 : break;
678 : }
679 :
680 0 : default:
681 0 : libmesh_error_msg("Unsupported triangular element: " << type_str);
682 : break;
683 : }
684 : } // if Tri
685 :
686 :
687 : // Set the target_elem equal to the reference elem
688 : else
689 7265 : for (const auto & node : target_elem->reference_elem()->node_ref_range())
690 7094 : owned_nodes.emplace_back(Node::build(node, node.id()));
691 :
692 : // Set nodes of target element
693 8312 : for (const auto & node_ptr : owned_nodes)
694 7774 : target_elem->set_node(node_ptr->id(), node_ptr.get());
695 :
696 24 : libmesh_assert(relative_fuzzy_equals(target_elem->volume(), ref_vol));
697 :
698 782 : return std::make_pair(std::move(target_elem), std::move(owned_nodes));
699 710 : }
700 :
701 758 : void VariationalSmootherSystem::get_target_to_reference_jacobian(
702 : const Elem * const target_elem,
703 : const FEMContext & femcontext,
704 : std::vector<RealTensor> & jacobians,
705 : std::vector<Real> & jacobian_dets)
706 : {
707 :
708 758 : const auto dim = target_elem->dim();
709 :
710 24 : const auto & qrule_points = femcontext.get_element_qrule().get_points();
711 24 : const auto & qrule_weights = femcontext.get_element_qrule().get_weights();
712 24 : const auto nq_points = femcontext.get_element_qrule().n_points();
713 :
714 : // If the target element is the reference element, Jacobian matrix is
715 : // identity, det of inverse is 1. These will only be overwritten if a
716 : // different target element is explicitly specified.
717 782 : jacobians = std::vector<RealTensor>(nq_points, RealTensor(
718 : 1., 0., 0.,
719 : 0., 1., 0.,
720 24 : 0., 0., 1.));
721 758 : jacobian_dets = std::vector<Real>(nq_points, 1.0);
722 :
723 : // Don't use "if (*target_elem == *(target_elem->reference_elem()))" here, it
724 : // only compares global node ids, not the node locations themselves.
725 24 : bool target_equals_reference = true;
726 758 : const auto * ref_elem = target_elem->reference_elem();
727 8312 : for (const auto local_id : make_range(target_elem->n_nodes()))
728 7554 : target_equals_reference &= target_elem->node_ref(local_id) == ref_elem->node_ref(local_id);
729 758 : if (target_equals_reference)
730 545 : return;
731 :
732 : // Create FEMap to compute target_element mapping information
733 225 : FEMap fe_map_target;
734 :
735 : // pre-request mapping derivatives
736 6 : const auto & dxyzdxi = fe_map_target.get_dxyzdxi();
737 6 : const auto & dxyzdeta = fe_map_target.get_dxyzdeta();
738 6 : const auto & dxyzdzeta = fe_map_target.get_dxyzdzeta();
739 :
740 : // build map
741 213 : fe_map_target.init_reference_to_physical_map(dim, qrule_points, target_elem);
742 213 : fe_map_target.compute_map(dim, qrule_weights, target_elem, /*d2phi=*/false);
743 :
744 1278 : for (const auto qp : make_range(nq_points))
745 : {
746 : // We use Larisa's H notation to denote the reference-to-target jacobian
747 30 : RealTensor H;
748 1065 : switch (dim)
749 : {
750 0 : case 1: {
751 0 : H = RealTensor(
752 0 : dxyzdxi[qp](0), 0., 0.,
753 : 0., 1., 0.,
754 : 0., 0., 1.);
755 :
756 0 : break;
757 : }
758 1065 : case 2: {
759 1155 : H = RealTensor(
760 60 : dxyzdxi[qp](0), dxyzdeta[qp](0), 0.,
761 1125 : dxyzdxi[qp](1), dxyzdeta[qp](1), 0.,
762 : 0., 0., 1.);
763 :
764 1065 : break;
765 : }
766 0 : case 3: {
767 0 : H = RealTensor(dxyzdxi[qp], dxyzdeta[qp], dxyzdzeta[qp]).transpose();
768 :
769 0 : break;
770 : }
771 :
772 0 : default:
773 0 : libmesh_error_msg("Unsupported mesh dimension " << dim);
774 : }
775 :
776 : // The target-to-reference jacobian is the inverse of the
777 : // reference-to-target jacobian
778 1065 : jacobians[qp] = H.inverse();
779 1095 : jacobian_dets[qp] = jacobians[qp].det();
780 : }
781 201 : }
782 :
783 : } // namespace libMesh
|
