libMesh
variational_smoother_system.C
Go to the documentation of this file.
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 
34 // C++ includes
35 #include <functional> // std::reference_wrapper
36 
37 namespace libMesh
38 {
39 
41 
42 void VariationalSmootherSystem::assembly (bool get_residual,
43  bool get_jacobian,
44  bool apply_heterogeneous_constraints,
45  bool apply_no_constraints)
46 {
47  // Update the mesh based on the current variable values
48  auto & mesh = this->get_mesh();
49  this->solution->close();
50 
51  for (auto * node : mesh.local_node_ptr_range())
52  {
53  for (const auto d : make_range(mesh.mesh_dimension()))
54  {
55  const auto dof_id = node->dof_number(this->number(), d, 0);
56  // Update mesh
57  (*node)(d) = libmesh_real((*current_local_solution)(dof_id));
58  }
59  }
60 
61  SyncNodalPositions sync_object(mesh);
62  Parallel::sync_dofobject_data_by_id (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);
63 
64  FEMSystem::assembly(get_residual, get_jacobian, apply_heterogeneous_constraints, apply_no_constraints);
65 }
66 
68 {
69  auto & mesh = this->get_mesh();
70  const auto & elem_orders = mesh.elem_default_orders();
71  libmesh_error_msg_if(elem_orders.size() != 1,
72  "The variational smoother cannot be used for mixed-order meshes!");
73  const auto fe_order = *elem_orders.begin();
74  // Add a variable for each dimension of the mesh
75  // "r0" for x, "r1" for y, "r2" for z
76  for (const auto & d : make_range(mesh.mesh_dimension()))
77  this->add_variable ("r" + std::to_string(d), fe_order);
78 
79  // Do the parent's initialization after variables are defined
81 
82  // Set the current_local_solution to the current mesh for the initial guess
83  this->solution->close();
84 
85  for (auto * node : mesh.local_node_ptr_range())
86  {
87  for (const auto d : make_range(mesh.mesh_dimension()))
88  {
89  const auto dof_id = node->dof_number(this->number(), d, 0);
90  // Update solution
91  (*solution).set(dof_id, (*node)(d));
92  }
93  }
94 
95  this->prepare_for_smoothing();
96 }
97 
99 {
100  std::unique_ptr<DiffContext> con = this->build_context();
101  FEMContext & femcontext = cast_ref<FEMContext &>(*con);
102  this->init_context(femcontext);
103 
104  const auto & mesh = this->get_mesh();
105  const auto dim = mesh.mesh_dimension();
106 
107  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  const auto & fe_map = femcontext.get_element_fe(0)->get_fe_map();
113  const auto & JxW = fe_map.get_JxW();
114 
115  std::map<ElemType, std::vector<Real>> target_elem_inverse_jacobian_dets;
116 
117  for (const auto * elem : mesh.active_local_element_ptr_range())
118  {
119  femcontext.pre_fe_reinit(*this, elem);
120  femcontext.elem_fe_reinit();
121 
122  // Add target element info, if applicable
123  if (_target_inverse_jacobians.find(elem->type()) == _target_inverse_jacobians.end())
124  {
125  // Create FEMap to compute target_element mapping information
126  FEMap fe_map_target;
127 
128  // pre-request mapping derivatives
129  const auto & dxyzdxi = fe_map_target.get_dxyzdxi();
130  const auto & dxyzdeta = fe_map_target.get_dxyzdeta();
131  // const auto & dxyzdzeta = fe_map_target.get_dxyzdzeta();
132 
133  const auto & qrule_points = femcontext.get_element_qrule().get_points();
134  const auto & qrule_weights = femcontext.get_element_qrule().get_weights();
135  const auto nq_points = femcontext.get_element_qrule().n_points();
136 
137  // If the target element is the reference element, Jacobian matrix is
138  // identity, det of inverse is 1. This will only be overwritten if a
139  // different target elemen is explicitly specified.
140  target_elem_inverse_jacobian_dets[elem->type()] =
141  std::vector<Real>(nq_points, 1.0);
142 
143  const auto type_str = Utility::enum_to_string(elem->type());
144 
145  // Elems deriving from Tri
146  if (type_str.compare(0, 3, "TRI") == 0)
147  {
148 
149  // The target element will be an equilateral triangle with area equal to
150  // the area of the reference element.
151 
152  // First, let's define the nodal locations of the vertices
153  const Real sqrt_3 = std::sqrt(3.);
154 
155  std::vector<Point> equilateral_points{
156  Point(0., 0.),
157  Point(1., 0.),
158  Point(0.5, 0.5 * sqrt_3)
159  };
160 
161  // Target element
162  const auto target_elem = Elem::build(elem->type());
163 
164  // Area of the reference element
165  const auto ref_area = target_elem->reference_elem()->volume();
166 
167  switch (elem->type())
168  {
169  case TRI3:
170  {
171  // Nothing to do here, vertices already defined in equilateral_points
172  break;
173  }
174 
175  case TRI6:
176  {
177  // Define the midpoint nodes of the equilateral triangle
178  equilateral_points.emplace_back(0.50, 0. );
179  equilateral_points.emplace_back(0.75, 0.25 * sqrt_3);
180  equilateral_points.emplace_back(0.25, 0.25 * sqrt_3);
181 
182  break;
183  }
184 
185  default:
186  libmesh_error_msg("Unsupported triangular element: " << type_str);
187  break;
188  }
189 
190  // Equilateral triangle side length preserving area of the reference element
191  const auto side_length = std::sqrt(4. / sqrt_3 * ref_area);
192 
193  std::vector<std::unique_ptr<Node>> owned_nodes;
194 
195  // Set nodes of target element to form an equilateral triangle
196  for (const dof_id_type node_id : index_range(equilateral_points))
197  {
198  // Scale the nodal positions to conserve area, store the pointer to keep it alive
199  owned_nodes.emplace_back(
200  Node::build(equilateral_points[node_id] * side_length, node_id));
201  target_elem->set_node(node_id, owned_nodes.back().get());
202  }
203 
204  // build map
205  fe_map_target.init_reference_to_physical_map(dim, qrule_points, target_elem.get());
206  fe_map_target.compute_map(dim, qrule_weights, target_elem.get(), /*d2phi=*/false);
207 
208  // Yes, triangle-to-triangle mappings have constant Jacobians, but we
209  // will keep things general for now
210  _target_inverse_jacobians[target_elem->type()] =
211  std::vector<RealTensor>(nq_points);
212  for (const auto qp : make_range(nq_points))
213  {
214  const RealTensor H_inv =
215  RealTensor(dxyzdxi[qp](0), dxyzdeta[qp](0), 0, dxyzdxi[qp](1),
216  dxyzdeta[qp](1), 0, 0, 0, 1)
217  .inverse();
218 
219  _target_inverse_jacobians[target_elem->type()][qp] = H_inv;
220  target_elem_inverse_jacobian_dets[target_elem->type()][qp] =
221  H_inv.det();
222  }
223  }// if Tri
224  }// if find == end()
225 
226  // Reference volume computation
227  Real elem_integrated_det_J = 0.;
228  for (const auto qp : index_range(JxW))
229  elem_integrated_det_J +=
230  JxW[qp] * target_elem_inverse_jacobian_dets[elem->type()][qp];
231  const auto ref_elem_vol = elem->reference_elem()->volume();
232  elem_averaged_det_J_sum += elem_integrated_det_J / ref_elem_vol;
233 
234  } // for elem
235 
236  // Get contributions from elements on other processors
237  mesh.comm().sum(elem_averaged_det_J_sum);
238 
239  _ref_vol = elem_averaged_det_J_sum / mesh.n_active_elem();
240 }
241 
243 {
244  FEMContext & c = cast_ref<FEMContext &>(context);
245 
246  FEBase * my_fe = nullptr;
247 
248  // Now make sure we have requested all the data
249  // we need to build the system.
250 
251  // We might have a multi-dimensional mesh
252  const std::set<unsigned char> & elem_dims =
253  c.elem_dimensions();
254 
255  for (const auto & dim : elem_dims)
256  {
257  c.get_element_fe( 0, my_fe, dim );
258  my_fe->get_nothing();
259 
260  auto & fe_map = my_fe->get_fe_map();
261  fe_map.get_dxyzdxi();
262  fe_map.get_dxyzdeta();
263  fe_map.get_dxyzdzeta();
264  fe_map.get_JxW();
265 
266  c.get_side_fe( 0, my_fe, dim );
267  my_fe->get_nothing();
268  }
269 
270  FEMSystem::init_context(context);
271 }
272 
273 
275  DiffContext & context)
276 {
277  FEMContext & c = cast_ref<FEMContext &>(context);
278 
279  const Elem & elem = c.get_elem();
280 
281  // First we get some references to cell-specific data that
282  // will be used to assemble the linear system.
283 
284  const auto & fe_map = c.get_element_fe(0)->get_fe_map();
285 
286  const auto & dxyzdxi = fe_map.get_dxyzdxi();
287  const auto & dxyzdeta = fe_map.get_dxyzdeta();
288  const auto & dxyzdzeta = fe_map.get_dxyzdzeta();
289 
290  const auto & dphidxi_map = fe_map.get_dphidxi_map();
291  const auto & dphideta_map = fe_map.get_dphideta_map();
292  const auto & dphidzeta_map = fe_map.get_dphidzeta_map();
293 
294  unsigned int dim = c.get_dim();
295 
296  unsigned int x_var = 0, y_var = 1, z_var = 2;
297  // In the case of lower dimensions, we will not access z (1D/2D) or y (1D)
298  if (dim < 3)
299  {
300  z_var = 0;
301  if (dim < 2)
302  y_var = 0;
303  }
304 
305  // The subvectors and submatrices we need to fill:
306  // system residual
307  std::reference_wrapper<DenseSubVector<Number>> F[3] =
308  {
309  c.get_elem_residual(x_var),
310  c.get_elem_residual(y_var),
311  c.get_elem_residual(z_var)
312  };
313  // system jacobian
314  std::reference_wrapper<DenseSubMatrix<Number>> K[3][3] =
315  {
316  {c.get_elem_jacobian(x_var, x_var), c.get_elem_jacobian(x_var, y_var), c.get_elem_jacobian(x_var, z_var)},
317  {c.get_elem_jacobian(y_var, x_var), c.get_elem_jacobian(y_var, y_var), c.get_elem_jacobian(y_var, z_var)},
318  {c.get_elem_jacobian(z_var, x_var), c.get_elem_jacobian(z_var, y_var), c.get_elem_jacobian(z_var, z_var)}
319  };
320 
321  // Quadrature info
322  const auto quad_weights = c.get_element_qrule().get_weights();
323 
324  const auto distortion_weight = 1. - _dilation_weight;
325 
326  // Integrate the distortion-dilation metric over the reference element
327  for (const auto qp : index_range(quad_weights))
328  {
329  // Compute quantities needed to evaluate the distortion-dilation metric
330  // and its gradient and Hessian.
331  // The metric will be minimized when it's gradient with respect to the node
332  // locations (R) is zero. For Newton's method, minimizing the gradient
333  // requires computation of the gradient's Jacobian with respect to R.
334  // The Jacobian of the distortion-dilation metric's gradientis the Hessian
335  // of the metric.
336  //
337  // Note that the term "Jacobian" has two meanings in this mesh smoothing
338  // application. The first meaning refers to the Jacobian w.r.t R of the
339  // gradient w.r.t. R, (i.e., the Hessian of the metric). This is the K
340  // variable defined above. This is also the Jacobian the 'request_jacobian'
341  // variable refers to. The second mesning refers to the Jacobian of the
342  // physical-to-reference element mapping. This is the Jacobian used to
343  // compute the distorion-dilation metric.
344  //
345  // Grab the physical-to-reference mapping Jacobian matrix (i.e., "S") at this qp
346  RealTensor S;
347  // RealTensors are always 3x3, so we will fill any dimensions above dim
348  // with 1s on the diagonal. This indicates a 1 to 1 relationship between
349  // the physical and reference elements in these extra dimensions.
350  switch (dim)
351  {
352  case 1:
353  S = RealTensor(
354  dxyzdxi[qp](0), 0, 0,
355  0, 1, 0,
356  0, 0, 1
357  );
358  break;
359 
360  case 2:
361  S = RealTensor(
362  dxyzdxi[qp](0), dxyzdeta[qp](0), 0,
363  dxyzdxi[qp](1), dxyzdeta[qp](1), 0,
364  0, 0, 1
365  );
366  break;
367 
368  case 3:
369  S = RealTensor(
370  dxyzdxi[qp],
371  dxyzdeta[qp],
372  dxyzdzeta[qp]
373  ).transpose(); // Note the transposition!
374  break;
375 
376  default:
377  libmesh_error_msg("Unsupported dimension.");
378  }
379 
380  // Apply any applicable target element transformations
381  if (_target_inverse_jacobians.find(elem.type()) !=
383  S = S * _target_inverse_jacobians[elem.type()][qp];
384 
385  // Compute quantities needed for the smoothing algorithm
386 
387  // determinant
388  const Real det = S.det();
389  const Real det_sq = det * det;
390  const Real det_cube = det_sq * det;
391 
392  const Real ref_vol_sq = _ref_vol * _ref_vol;
393 
394  // trace of S^T * S
395  // DO NOT USE RealTensor.tr for the trace, it will NOT be correct for
396  // 1D and 2D meshes because of our hack of putting 1s in the diagonal of
397  // S for the extra dimensions
398  const RealTensor ST_S = S.transpose() * S;
399  Real tr = 0.0;
400  for (const auto i : make_range(dim))
401  tr += ST_S(i, i);
402  const Real tr_div_dim = tr / dim;
403 
404  // Precompute pow(tr_div_dim, 0.5 * dim - x) for x = 0, 1, 2
405  const Real half_dim = 0.5 * dim;
406  const std::vector<Real> trace_powers{
407  std::pow(tr_div_dim, half_dim),
408  std::pow(tr_div_dim, half_dim - 1.),
409  std::pow(tr_div_dim, half_dim - 2.),
410  };
411 
412  // inverse of S
413  const RealTensor S_inv = S.inverse();
414  // inverse transpose of S
415  const RealTensor S_inv_T = S_inv.transpose();
416 
417  // Identity matrix
418  const RealTensor I(1, 0, 0,
419  0, 1, 0,
420  0, 0, 1);
421 
422  // The chi function allows us to handle degenerate elements
423  // When the element is not degenerate (i.e., det(S) > 0), we can set
424  // _epsilon_squared to zero to get chi(det(S)) = det(S)
425  const Real epsilon_sq = det > 0. ? 0. : _epsilon_squared;
426  const Real chi = 0.5 * (det + std::sqrt(epsilon_sq + det_sq));
427  const Real chi_sq = chi * chi;
428  const Real sqrt_term = std::sqrt(epsilon_sq + det_sq);
429  // dchi(x) / dx
430  const Real chi_prime = 0.5 * (1. + det / sqrt_term);
431  const Real chi_prime_sq = chi_prime * chi_prime;
432  // d2chi(x) / dx2
433  const Real chi_2prime = 0.5 * (1. / sqrt_term - det_sq / Utility::pow<3>(sqrt_term));
434 
435  // Distortion metric (beta)
436  //const Real beta = trace_powers[0] / chi;
437  const RealTensor dbeta_dS = (trace_powers[1] / chi) * S - (trace_powers[0] / chi_sq * chi_prime * det) * S_inv_T;
438 
439  // Dilation metric (mu)
440  //const Real mu = 0.5 * (_ref_vol + det_sq / _ref_vol) / chi;
441  // We represent d mu / dS as alpha(S) * S^-T, where alpha is a scalar function
442  const Real alpha = (-chi_prime * det_cube + 2. * det_sq * chi - ref_vol_sq * det * chi_prime) / (2. * _ref_vol * chi_sq);
443  const RealTensor dmu_dS = alpha * S_inv_T;
444 
445  // Combined metric (E)
446  //const Real E = distortion_weight * beta + _dilation_weight * mu;
447  const RealTensor dE_dS = distortion_weight * dbeta_dS + _dilation_weight * dmu_dS;
448 
449  // This vector is useful in computing dS/dR below
450  std::vector<std::vector<std::vector<Real>>> dphi_maps = {dphidxi_map, dphideta_map, dphidzeta_map};
451 
452  // Compute residual (i.e., the gradient of the combined metric w.r.t node locations)
453  // Recall that when the gradient (residual) is zero, the combined metric is minimized
454  for (const auto l : elem.node_index_range())
455  {
456  for (const auto var_id : make_range(dim))
457  {
458  // Build dS/dR, the derivative of the physical-to-reference mapping Jacobin w.r.t. the mesh node locations
459  RealTensor dS_dR = RealTensor(0);
460  for (const auto jj : make_range(dim))
461  dS_dR(var_id, jj) = dphi_maps[jj][l][qp];
462 
463  // Residual contribution. The contraction of dE/dS and dS/dR gives us
464  // the gradient we are looking for, dE/dR
465  F[var_id](l) += quad_weights[qp] * dE_dS.contract(dS_dR);
466  }// for var_id
467  }// for l
468 
469 
470  if (request_jacobian)
471  {
472  // Compute jacobian of the smoothing system (i.e., the Hessian of the
473  // combined metric w.r.t. the mesh node locations)
474 
475  // Precompute coefficients to be applied to each component of tensor
476  // products in the loops below. At first glance, these coefficients look
477  // like gibberish, but everything in the Hessian has been verified by
478  // taking finite differences of the gradient. We should probably write
479  // down the derivations of the gradient and Hessian somewhere...
480 
481  // Recall that above, dbeta_dS takes the form:
482  // d(beta)/dS = c1(S) * S - c2(S) * S_inv_T,
483  // where c1 and c2 are scalar-valued functions.
484  const std::vector<Real> d2beta_dS2_coefs_times_distortion_weight = {
485  //Part 1: scaler coefficients of d(c1 * S) / dS
486  //
487  // multiplies I[i,a] x I[j,b]
488  (trace_powers[1] / chi) * distortion_weight,
489  // multiplies S[a,b] x S[i,j]
490  (((dim - 2.) / dim) * trace_powers[2] / chi) * distortion_weight,
491  // multiplies S_inv[b,a] * S[i,j]
492  (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
493  //
494  //Part 2: scaler coefficients of d(-c2 * S_inv_T) / dS
495  //
496  // multiplies S[a,b] x S_inv[j,i]
497  (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
498  // multiplies S_inv[b,a] x S_inv[j,i]
499  (trace_powers[0] * (det / chi_sq)
500  * ((2. * chi_prime_sq / chi - chi_2prime) * det - chi_prime)) * distortion_weight,
501  // multiplies S_inv[b,i] x S_inv[j,a]
502  ((trace_powers[0] / chi_sq) * chi_prime * det) * distortion_weight,
503  };
504 
505  // d alpha / dS has the form c(S) * S^-T, where c is the scalar coefficient defined below
506  const Real dalpha_dS_coef_times_dilation_weight = ((det / (2. * _ref_vol * chi_sq))
507  * (-4. * _ref_vol * alpha * chi * chi_prime - chi_2prime * det_cube
508  - chi_prime * det_sq + 4 * chi * det - ref_vol_sq * (chi_prime + det * chi_2prime))) * _dilation_weight;
509 
510  // This is also useful to precompute
511  const Real alpha_times_dilation_weight = alpha * _dilation_weight;
512 
513  /*
514 
515  To increase the efficiency of the Jacobian computation, we take
516  advantage of l-p symmetry, ij-ab symmetry, and the sparsity pattern of
517  dS_dR. We also factor all possible multipliers out of inner loops for
518  efficiency. The result is code that is more difficult to read. For
519  clarity, consult the pseudo-code below in this comment.
520 
521  for (const auto l: elem.node_index_range()) // Contribution to Hessian
522  from node l
523  {
524  for (const auto var_id1 : make_range(dim)) // Contribution from each
525  x/y/z component of node l
526  {
527  // Build dS/dR_l, the derivative of the physical-to-reference
528  mapping Jacobin w.r.t.
529  // the l-th node
530  RealTensor dS_dR_l = RealTensor(0);
531  for (const auto ii : make_range(dim))
532  dS_dR_l(var_id1, ii) = dphi_maps[ii][l][qp];
533 
534  for (const auto p: elem.node_index_range()) // Contribution to
535  Hessian from node p
536  {
537  for (const auto var_id2 : make_range(dim)) // Contribution from
538  each x/y/z component of node p
539  {
540  // Build dS/dR_l, the derivative of the physical-to-reference
541  mapping Jacobin w.r.t.
542  // the p-th node
543  RealTensor dS_dR_p = RealTensor(0);
544  for (const auto jj : make_range(dim))
545  dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
546 
547  Real d2beta_dR2 = 0.;
548  Real d2mu_dR2 = 0.;
549  // Perform tensor contraction
550  for (const auto i : make_range(dim))
551  {
552  for (const auto j : make_range(dim))
553  {
554  for (const auto a : make_range(dim))
555  {
556  for (const auto b : make_range(dim))
557  {
558  // Nasty tensor products to be multiplied by
559  d2beta_dS2_coefs to get d2(beta) / dS2 const std::vector<Real>
560  d2beta_dS2_tensor_contributions =
561  {
562  I(i,a) * I(j,b),
563  S(a,b) * S(i,j),
564  S_inv(b,a) * S(i,j),
565  S(a,b) * S_inv(j,i),
566  S_inv(b,a) * S_inv(j,i),
567  S_inv(j,a) * S_inv(b,i),
568  };
569 
570  // Combine precomputed coefficients with tensor products
571  to get d2(beta) / dS2 Real d2beta_dS2 = 0.; for (const auto comp_id :
572  index_range(d2beta_dS2_coefs))
573  {
574  const Real contribution = d2beta_dS2_coefs[comp_id] *
575  d2beta_dS2_tensor_contributions[comp_id]; d2beta_dS2 += contribution;
576  }
577 
578  // Incorporate tensor product portion to get d2(mu) /
579  dS2 const Real d2mu_dS2 = dalpha_dS_coef * S_inv(b,a) * S_inv(j,i) -
580  alpha * S_inv(b,i) * S_inv(j,a);
581 
582  // Chain rule to change d/dS to d/dR
583  d2beta_dR2 += d2beta_dS2 * dS_dR_l(a, b) * dS_dR_p(i,
584  j); d2mu_dR2 += d2mu_dS2 * dS_dR_l(a, b) * dS_dR_p(i, j);
585 
586  }// for b
587  }// for a
588  }// for j
589  }// for i, end tensor contraction
590 
591  // Jacobian contribution
592  K[var_id1][var_id2](l, p) += quad_weights[qp] * ((1. -
593  _dilation_weight) * d2beta_dR2 + _dilation_weight * d2mu_dR2);
594 
595  }// for var_id2
596  }// for p
597  }// for var_id1
598  }// for l
599 
600  End pseudo-code, begin efficient code
601 
602  */
603 
604  for (const auto l: elem.node_index_range()) // Contribution to Hessian from node l
605  {
606  for (const auto var_id1 : make_range(dim)) // Contribution from each x/y/z component of node l
607  {
608  // Build dS/dR_l, the derivative of the physical-to-reference mapping Jacobin w.r.t.
609  // the l-th node
610  RealTensor dS_dR_l = RealTensor(0);
611  for (const auto ii : make_range(dim))
612  dS_dR_l(var_id1, ii) = dphi_maps[ii][l][qp];
613 
614  // Jacobian is symmetric, only need to loop over lower triangular portion
615  for (const auto p: make_range(l + 1)) // Contribution to Hessian from node p
616  {
617  for (const auto var_id2 : make_range(dim)) // Contribution from each x/y/z component of node p
618  {
619  // Build dS/dR_l, the derivative of the physical-to-reference mapping Jacobin w.r.t.
620  // the p-th node
621  RealTensor dS_dR_p = RealTensor(0);
622  for (const auto jj : make_range(dim))
623  dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
624 
625  Real d2E_dR2 = 0.;
626  // Perform tensor contraction
627  for (const auto i : make_range(dim))
628  {
629 
630  for (const auto j : make_range(dim))
631  {
632 
633  const auto S_ij = S(i,j);
634  const auto S_inv_ji = S_inv(j,i);
635 
636  // Apply the stuff only depending on i and j before entering a, b loops
637  // beta
638  const std::vector<Real> d2beta_dS2_coefs_ij_applied{
639  d2beta_dS2_coefs_times_distortion_weight[1] * S_ij,
640  d2beta_dS2_coefs_times_distortion_weight[2] * S_ij,
641  d2beta_dS2_coefs_times_distortion_weight[3] * S_inv_ji,
642  d2beta_dS2_coefs_times_distortion_weight[4] * S_inv_ji,
643  };
644  // mu
645  const auto dalpha_dS_coef_ij_applied = dalpha_dS_coef_times_dilation_weight * S_inv_ji;
646 
647  Real d2E_dSdR_l = 0.;
648  Real d2E_dSdR_p = 0.;
649 
650  for (const auto a : make_range(i + 1))
651  {
652 
653  // If this condition is met, both the ijab and abij
654  // contributions to the Jacobian are zero due to the
655  // spasity patterns of dS_dR_l and dS_dR_p and this
656  // iteration may be skipped
657  if (!(a == var_id1 && i == var_id2) &&
658  !(a == var_id2 && i == var_id1))
659  continue;
660 
661  const auto S_inv_ja = S_inv(j,a);
662 
663  const Real d2beta_dS2_coef_ia_applied = d2beta_dS2_coefs_times_distortion_weight[0] * I(i,a);
664  const Real d2beta_dS2_coef_ja_applied = d2beta_dS2_coefs_times_distortion_weight[5] * S_inv_ja;
665  const Real alpha_ja_applied = alpha_times_dilation_weight * S_inv_ja;
666 
667  const auto b_limit = (a == i) ? j + 1 : dim;
668  for (const auto b : make_range(b_limit))
669  {
670 
671  // Combine precomputed coefficients with tensor products
672  // to get d2(beta) / dS2
673  Real d2beta_dS2_times_distortion_weight = (
674  d2beta_dS2_coef_ia_applied * I(j,b) +
675  d2beta_dS2_coefs_ij_applied[0] * S(a,b) +
676  d2beta_dS2_coefs_ij_applied[1] * S_inv(b,a) +
677  d2beta_dS2_coefs_ij_applied[2] * S(a,b) +
678  d2beta_dS2_coefs_ij_applied[3] * S_inv(b,a) +
679  d2beta_dS2_coef_ja_applied * S_inv(b,i)
680  );
681 
682  // Incorporate tensor product portion to get d2(mu) /
683  // dS2
684  const Real d2mu_dS2_times_dilation_weight = dalpha_dS_coef_ij_applied * S_inv(b,a) - alpha_ja_applied * S_inv(b,i);
685 
686 
687  // Chain rule to change d/dS to d/dR
688  const auto d2E_dS2 =
689  d2beta_dS2_times_distortion_weight +
690  d2mu_dS2_times_dilation_weight;
691 
692  // if !(a == var_id1 (next line) && i == var_id2
693  // (outside 'a' loop)), dS_dR_l(p) multiplier is zero
694  d2E_dSdR_l += d2E_dS2 * dS_dR_l(a, b);
695 
696  if (!(i == a && j == b))
697  // if !(a == var_id2 (next line) && i == var_id1
698  // (outside 'a' loop)), dS_dR_p(l) multiplier is zero
699  d2E_dSdR_p += d2E_dS2 * dS_dR_p(a, b);
700 
701  } // for b
702  } // for a
703  d2E_dR2 +=
704  d2E_dSdR_l * dS_dR_p(i, j) + d2E_dSdR_p * dS_dR_l(i, j);
705  }// for j
706  }// for i, end tensor contraction
707 
708  // Jacobian contribution
709  const Real jacobian_contribution = quad_weights[qp] * d2E_dR2;
710  K[var_id1][var_id2](l, p) += jacobian_contribution;
711  // Jacobian is symmetric, add contribution to p,l entry
712  // Don't get the diagonal twice!
713  if (p < l)
714  // Note the transposition of var_id1 and var_id2 as these are also jacobian indices
715  K[var_id2][var_id1](p, l) += jacobian_contribution;
716 
717  }// for var_id2
718  }// for p
719  }// for var_id1
720  }// for l
721  }
722  } // end of the quadrature point qp-loop
723 
724  return request_jacobian;
725 }
726 
727 } // namespace libMesh
T libmesh_real(T a)
virtual std::unique_ptr< DiffContext > build_context() override
Builds a FEMContext object with enough information to do evaluations on each element.
Definition: fem_system.C:1346
virtual void init_data() override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:274
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:55
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
Definition: fem_context.h:317
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
Definition: fem_context.C:1683
virtual void compute_map(const unsigned int dim, const std::vector< Real > &qw, const Elem *elem, bool calculate_d2phi)
Compute the jacobian and some other additional data fields.
Definition: fe_map.C:1439
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:908
unsigned int dim
virtual void init_context(DiffContext &) override
Definition: fem_system.C:1370
TypeTensor< T > transpose() const
Definition: type_tensor.h:1050
std::map< ElemType, std::vector< RealTensor > > _target_inverse_jacobians
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) override
Assembly method to update the mesh based on the smoother solve.
MeshBase & mesh
const std::vector< Real > & get_weights() const
Definition: quadrature.h:168
RealTensorValue RealTensor
The libMesh namespace provides an interface to certain functionality in the library.
const MeshBase & get_mesh() const
Definition: system.h:2358
Real _dilation_weight
The relative weight to give the dilation metric. The distortion metric is given weight 1 - _dilation_...
unsigned char get_dim() const
Accessor for cached mesh dimension.
Definition: fem_context.h:936
virtual void elem_fe_reinit(const std::vector< Point > *const pts=nullptr)
Reinitializes interior FE objects on the current geometric element.
Definition: fem_context.C:1477
unsigned int number() const
Definition: system.h:2350
T pow(const T &x)
Definition: utility.h:328
const std::vector< RealGradient > & get_dxyzdxi() const
Definition: fe_map.h:239
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
virtual void init_data() override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: fem_system.C:873
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:62
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1357
unsigned int n_points() const
Definition: quadrature.h:131
virtual bool element_time_derivative(bool request_jacobian, libMesh::DiffContext &context) override
Adds the time derivative contribution on elem to elem_residual.
void sync_dofobject_data_by_id(const Communicator &comm, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
Request data about a range of ghost dofobjects uniquely identified by their id.
CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
Multiply 2 tensors together to return a scalar, i.e.
Definition: type_tensor.h:1268
std::string enum_to_string(const T e)
void init_reference_to_physical_map(const std::vector< Point > &qp, const Elem *elem)
Definition: fe_map.C:109
Real _ref_vol
The reference volume for each element.
static std::unique_ptr< Node > build(const Node &n)
Definition: node.h:315
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:242
const Real _epsilon_squared
The small nonzero constant to prevent zero denominators (degenerate elements only) ...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< Point > & get_points() const
Definition: quadrature.h:156
TypeTensor< T > inverse() const
Definition: type_tensor.h:1080
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1605
virtual void init_context(libMesh::DiffContext &context) override
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:277
const std::vector< RealGradient > & get_dxyzdeta() const
Definition: fe_map.h:247
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) override
Prepares matrix or rhs for matrix assembly.
Definition: fem_system.C:880
Class contained in FE that encapsulates mapping (i.e.
Definition: fe_map.h:47
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
This class forms the foundation from which generic finite elements may be derived.
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:802
const std::set< unsigned char > & elem_dimensions() const
Definition: fem_context.h:951
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
uint8_t dof_id_type
Definition: id_types.h:67