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 #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 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  comm.minloc(pair.first, rank);
49  comm.broadcast(pair.second, rank);
50 }
51 
52 /*
53  * Gets the dof_id_type value corresponding to the maximum of the Real value.
54  */
55 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  comm.maxloc(pair.first, rank);
60  comm.broadcast(pair.second, rank);
61 }
62 
66 Real chi_epsilon(const Real & x, const Real epsilon_squared)
67 {
68  return 0.5 * (x + std::sqrt(epsilon_squared + Utility::pow<2>(x)));
69 }
70 
76  const unsigned int & dim,
77  const unsigned int & qp)
78 {
79  const auto & dxyzdxi = fe_map.get_dxyzdxi()[qp];
80  const auto & dxyzdeta = fe_map.get_dxyzdeta()[qp];
81  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  switch (dim)
87  {
88  case 1:
89  return RealTensor(
90  dxyzdxi(0), 0, 0,
91  0, 1, 0,
92  0, 0, 1
93  );
94  break;
95 
96  case 2:
97  return RealTensor(
98  dxyzdxi(0), dxyzdeta(0), 0,
99  dxyzdxi(1), dxyzdeta(1), 0,
100  0, 0, 1
101  );
102  break;
103 
104  case 3:
105  return RealTensor(
106  dxyzdxi,
107  dxyzdeta,
108  dxyzdzeta
109  ).transpose(); // Note the transposition!
110  break;
111 
112  default:
113  libmesh_error_msg("Unsupported dimension.");
114  }
115 }
116 
120 Real trace(const RealTensor & A, const unsigned int & dim)
121 {
122  Real tr = 0.0;
123  for (const auto i : make_range(dim))
124  tr += A(i, i);
125 
126  return tr;
127 }
128 
130 
131 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  auto & mesh = this->get_mesh();
138  this->solution->close();
139 
140  for (auto * node : mesh.local_node_ptr_range())
141  {
142  for (const auto d : make_range(mesh.mesh_dimension()))
143  {
144  const auto dof_id = node->dof_number(this->number(), d, 0);
145  // Update mesh
146  (*node)(d) = libmesh_real((*current_local_solution)(dof_id));
147  }
148  }
149 
150  SyncNodalPositions sync_object(mesh);
151  Parallel::sync_dofobject_data_by_id (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);
152 
153  // Compute and update mesh quality information
155  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  if (_untangling_solve)
160  {
161  const Real & min_S = _mesh_info.min_qp_det_S;
162  // This component is based on what Larisa did in the original code.
163  const Real variable_component = 100. * Utility::pow<2>(_ref_vol * min_S);
164  _epsilon_squared_assembly = _epsilon_squared + variable_component;
165  }
166 
167  else
169 
170  FEMSystem::assembly(get_residual, get_jacobian, apply_heterogeneous_constraints, apply_no_constraints);
171 
172  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  rhs->close();
180  (*rhs) *= 0.9;
181  }
182 }
183 
185 {
186  const auto & mesh_info = get_mesh_info();
187  if (mesh_info.mesh_is_tangled)
188  {
189  // Untangling solve
190  _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  const auto dilation_weight = _dilation_weight;
198  _dilation_weight = 0.;
199 
201 
202  // Reset the dilation weight
203  _dilation_weight = dilation_weight;
204  }
205 
206  // Smoothing solve
207  _untangling_solve = false;
209  libmesh_error_msg_if(get_mesh_info().mesh_is_tangled, "The smoothing solve tangled the mesh!");
210 }
211 
213 {
214  auto & mesh = this->get_mesh();
215  const auto & elem_orders = mesh.elem_default_orders();
216  libmesh_error_msg_if(elem_orders.size() != 1,
217  "The variational smoother cannot be used for mixed-order meshes!");
218  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  for (const auto & d : make_range(mesh.mesh_dimension()))
222  this->add_variable ("r" + std::to_string(d), fe_order);
223 
224  // Do the parent's initialization after variables are defined
226 
227  // Set the current_local_solution to the current mesh for the initial guess
228  this->solution->close();
229 
230  for (auto * node : mesh.local_node_ptr_range())
231  {
232  for (const auto d : make_range(mesh.mesh_dimension()))
233  {
234  const auto dof_id = node->dof_number(this->number(), d, 0);
235  // Update solution
236  (*solution).set(dof_id, (*node)(d));
237  }
238  }
239 
240  this->prepare_for_smoothing();
241 }
242 
244 {
245  // If this method has already been called to set _ref_vol, just return.
246  if (std::abs(_ref_vol) > TOLERANCE * TOLERANCE)
247  return;
248 
249  std::unique_ptr<DiffContext> con = this->build_context();
250  FEMContext & femcontext = cast_ref<FEMContext &>(*con);
251  this->init_context(femcontext);
252 
253  const auto & mesh = this->get_mesh();
254 
255  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  const auto & fe_map = femcontext.get_element_fe(0)->get_fe_map();
261  const auto & JxW = fe_map.get_JxW();
262 
263  for (const auto * elem : mesh.active_local_element_ptr_range())
264  {
265  femcontext.pre_fe_reinit(*this, elem);
266  femcontext.elem_fe_reinit();
267 
268  // Add target element info, if applicable
269  if (_target_jacobians.find(elem->type()) == _target_jacobians.end())
270  {
271  const auto [target_elem, target_nodes] = get_target_elem(elem->type());
272  get_target_to_reference_jacobian(target_elem.get(),
273  femcontext,
274  _target_jacobians[elem->type()],
275  _target_jacobian_dets[elem->type()]);
276  }// if find == end()
277 
278  // Reference volume computation
279  Real elem_integrated_det_S = 0.;
280  for (const auto qp : index_range(JxW))
281  elem_integrated_det_S += JxW[qp] * _target_jacobian_dets[elem->type()][qp];
282  const auto ref_elem_vol = elem->reference_elem()->volume();
283  elem_averaged_det_S_sum += elem_integrated_det_S / ref_elem_vol;
284 
285  } // for elem
286 
287  // Get contributions from elements on other processors
288  mesh.comm().sum(elem_averaged_det_S_sum);
289 
290  _ref_vol = elem_averaged_det_S_sum / mesh.n_active_elem();
291 }
292 
294 {
295  FEMContext & c = cast_ref<FEMContext &>(context);
296 
297  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  c.elem_dimensions();
305 
306  for (const auto & dim : elem_dims)
307  {
308  c.get_element_fe( 0, my_fe, dim );
309  my_fe->get_nothing();
310 
311  auto & fe_map = my_fe->get_fe_map();
312  fe_map.get_dxyzdxi();
313  fe_map.get_dxyzdeta();
314  fe_map.get_dxyzdzeta();
315  fe_map.get_JxW();
316 
317  // Mesh may be tangled, allow negative Jacobians
318  fe_map.set_jacobian_tolerance(std::numeric_limits<Real>::lowest());
319 
320  c.get_side_fe( 0, my_fe, dim );
321  my_fe->get_nothing();
322  }
323 
324  FEMSystem::init_context(context);
325 }
326 
327 
329  DiffContext & context)
330 {
331  FEMContext & c = cast_ref<FEMContext &>(context);
332 
333  const Elem & elem = c.get_elem();
334 
335  unsigned int dim = c.get_dim();
336 
337  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  if (dim < 3)
340  {
341  z_var = 0;
342  if (dim < 2)
343  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  };
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  };
361 
362  // Quadrature info
363  const auto quad_weights = c.get_element_qrule().get_weights();
364 
365  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  const auto & fe_map = c.get_element_fe(0)->get_fe_map();
371 
372  const auto & dphidxi_map = fe_map.get_dphidxi_map();
373  const auto & dphideta_map = fe_map.get_dphideta_map();
374  const auto & dphidzeta_map = fe_map.get_dphidzeta_map();
375 
376  // Integrate the distortion-dilation metric over the reference element
377  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  RealTensor S = get_jacobian_at_qp(fe_map, dim, qp);
397 
398  // Apply target element transformation
399  S *= _target_jacobians[elem.type()][qp];
400 
401  // Compute quantities needed for the smoothing algorithm
402 
403  // determinant
404  const Real det = S.det();
405  const Real det_sq = det * det;
406  const Real det_cube = det_sq * det;
407 
408  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  const auto tr = trace(S.transpose() * S, dim);
415  const Real tr_div_dim = tr / dim;
416 
417  // Precompute pow(tr_div_dim, 0.5 * dim - x) for x = 0, 1, 2
418  const Real half_dim = 0.5 * dim;
419  const std::vector<Real> trace_powers{
420  std::pow(tr_div_dim, half_dim),
421  std::pow(tr_div_dim, half_dim - 1.),
422  std::pow(tr_div_dim, half_dim - 2.),
423  };
424 
425  // inverse of S
426  const RealTensor S_inv = S.inverse();
427  // inverse transpose of S
428  const RealTensor S_inv_T = S_inv.transpose();
429 
430  // Identity matrix
431  const RealTensor I(1, 0, 0,
432  0, 1, 0,
433  0, 0, 1);
434 
435  // The chi function allows us to handle degenerate elements
436  const auto chi = chi_epsilon(det, _epsilon_squared_assembly);
437  const Real chi_sq = chi * chi;
438  const Real sqrt_term = std::sqrt(_epsilon_squared_assembly + det_sq);
439  // dchi(x) / dx
440  const Real chi_prime = 0.5 * (1. + det / sqrt_term);
441  const Real chi_prime_sq = chi_prime * chi_prime;
442  // d2chi(x) / dx2
443  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  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  const Real alpha = (-chi_prime * det_cube + 2. * det_sq * chi - ref_vol_sq * det * chi_prime) / (2. * _ref_vol * chi_sq);
453  const RealTensor dmu_dS = alpha * S_inv_T;
454 
455  // Combined metric (E)
456  //const Real E = distortion_weight * beta + _dilation_weight * mu;
457  const RealTensor dE_dS = distortion_weight * dbeta_dS + _dilation_weight * dmu_dS;
458 
459  // This vector is useful in computing dS/dR below
460  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  for (const auto l : elem.node_index_range())
465  {
466  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  RealTensor dS_dR = RealTensor(0);
470  for (const auto jj : make_range(dim))
471  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  F[var_id](l) += quad_weights[qp] * dE_dS.contract(dS_dR);
476  }// for var_id
477  }// for l
478 
479 
480  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  (trace_powers[1] / chi) * distortion_weight,
499  // multiplies S[a,b] x S[i,j]
500  (((dim - 2.) / dim) * trace_powers[2] / chi) * distortion_weight,
501  // multiplies S_inv[b,a] * S[i,j]
502  (-(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  (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
508  // multiplies S_inv[b,a] x S_inv[j,i]
509  (trace_powers[0] * (det / chi_sq)
510  * ((2. * chi_prime_sq / chi - chi_2prime) * det - chi_prime)) * distortion_weight,
511  // multiplies S_inv[b,i] x S_inv[j,a]
512  ((trace_powers[0] / chi_sq) * chi_prime * det) * distortion_weight,
513  };
514 
515  // d alpha / dS has the form c(S) * S^-T, where c is the scalar coefficient defined below
516  const Real dalpha_dS_coef_times_dilation_weight = ((det / (2. * _ref_vol * chi_sq))
517  * (-4. * _ref_vol * alpha * chi * chi_prime - chi_2prime * det_cube
518  - 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  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  for (const auto l: elem.node_index_range()) // Contribution to Hessian from node l
615  {
616  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  RealTensor dS_dR_l = RealTensor(0);
621  for (const auto ii : make_range(dim))
622  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  for (const auto p: make_range(l + 1)) // Contribution to Hessian from node p
626  {
627  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  RealTensor dS_dR_p = RealTensor(0);
632  for (const auto jj : make_range(dim))
633  dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
634 
635  Real d2E_dR2 = 0.;
636  // Perform tensor contraction
637  for (const auto i : make_range(dim))
638  {
639 
640  for (const auto j : make_range(dim))
641  {
642 
643  const auto S_ij = S(i,j);
644  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  d2beta_dS2_coefs_times_distortion_weight[1] * S_ij,
650  d2beta_dS2_coefs_times_distortion_weight[2] * S_ij,
651  d2beta_dS2_coefs_times_distortion_weight[3] * S_inv_ji,
652  d2beta_dS2_coefs_times_distortion_weight[4] * S_inv_ji,
653  };
654  // mu
655  const auto dalpha_dS_coef_ij_applied = dalpha_dS_coef_times_dilation_weight * S_inv_ji;
656 
657  Real d2E_dSdR_l = 0.;
658  Real d2E_dSdR_p = 0.;
659 
660  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  if (!(a == var_id1 && i == var_id2) &&
668  !(a == var_id2 && i == var_id1))
669  continue;
670 
671  const auto S_inv_ja = S_inv(j,a);
672 
673  const Real d2beta_dS2_coef_ia_applied = d2beta_dS2_coefs_times_distortion_weight[0] * I(i,a);
674  const Real d2beta_dS2_coef_ja_applied = d2beta_dS2_coefs_times_distortion_weight[5] * S_inv_ja;
675  const Real alpha_ja_applied = alpha_times_dilation_weight * S_inv_ja;
676 
677  const auto b_limit = (a == i) ? j + 1 : dim;
678  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  d2beta_dS2_coef_ia_applied * I(j,b) +
685  d2beta_dS2_coefs_ij_applied[0] * S(a,b) +
686  d2beta_dS2_coefs_ij_applied[1] * S_inv(b,a) +
687  d2beta_dS2_coefs_ij_applied[2] * S(a,b) +
688  d2beta_dS2_coefs_ij_applied[3] * S_inv(b,a) +
689  d2beta_dS2_coef_ja_applied * S_inv(b,i)
690  );
691 
692  // Incorporate tensor product portion to get d2(mu) /
693  // dS2
694  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  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  d2E_dSdR_l += d2E_dS2 * dS_dR_l(a, b);
705 
706  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  d2E_dSdR_p += d2E_dS2 * dS_dR_p(a, b);
710 
711  } // for b
712  } // for a
713  d2E_dR2 +=
714  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  const Real jacobian_contribution = quad_weights[qp] * d2E_dR2;
720  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  if (p < l)
724  // Note the transposition of var_id1 and var_id2 as these are also jacobian indices
725  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  } // end of the quadrature point qp-loop
733 
734  return request_jacobian;
735 }
736 
738 {
739  if (!_mesh_info.initialized)
741 
742  return _mesh_info;
743 }
744 
746 {
747  // If the reference volume has not yet been computed, compute it.
748  if (std::abs(_ref_vol) < TOLERANCE * TOLERANCE)
750 
751  std::unique_ptr<DiffContext> con = this->build_context();
752  FEMContext & femcontext = cast_ref<FEMContext &>(*con);
753  this->init_context(femcontext);
754 
755  const auto & mesh = this->get_mesh();
756  const auto dim = mesh.mesh_dimension();
757  const Real half_dim = 0.5 * dim;
758  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  const auto & fe_map = femcontext.get_element_fe(0)->get_fe_map();
764  const auto & JxW = fe_map.get_JxW();
765  fe_map.get_dxyzdxi();
766  fe_map.get_dxyzdeta();
767  fe_map.get_dxyzdzeta();
768 
770 
771  for (const auto * elem : mesh.active_local_element_ptr_range())
772  {
773  femcontext.pre_fe_reinit(*this, elem);
774  femcontext.elem_fe_reinit();
775 
776  // Element-integrated quantities
777  Real det_S_int = 0.;
778  Real beta_int = 0.;
779  Real mu_int = 0.;
780  Real combined_int = 0.;
781 
782  for (const auto qp : index_range(JxW))
783  {
784  const auto det_SxW = JxW[qp] * _target_jacobian_dets[elem->type()][qp];
785  det_S_int += det_SxW;
786 
787  // Grab the physical-to-reference mapping Jacobian matrix (i.e., "S") at this qp
788  RealTensor S = get_jacobian_at_qp(fe_map, dim, qp);
789 
790  // Apply target element transformation
791  S *= _target_jacobians[elem->type()][qp];
792 
793  // Determinant of S
794  const auto det = S.det();
795  const auto det_sq = det * det;
796 
797  if (det > info.max_qp_det_S)
798  info.max_qp_det_S = det;
799  else if (det < info.min_qp_det_S)
800  info.min_qp_det_S = det;
801 
802  if (det < TOLERANCE * TOLERANCE)
803  info.mesh_is_tangled = true;
804 
805  // trace of S^T * S
806  const auto tr = trace(S.transpose() * S, dim);
807 
808  // The chi function allows us to handle degenerate elements
809  const auto chi = chi_epsilon(det, _epsilon_squared_assembly);
810 
811  // distortion
812  const Real beta = std::pow(tr / dim, half_dim) / chi;
813  beta_int += beta * det_SxW;
814 
815  // dilation
816  const Real mu = 0.5 * (_ref_vol + det_sq / _ref_vol) / chi;
817  mu_int += mu * det_SxW;
818 
819  // combined
820  const Real E = distortion_weight * beta + _dilation_weight * mu;
821  combined_int += E * det_SxW;
822  }
823 
824  info.total_det_S += det_S_int;
825  if (det_S_int > info.max_elem_det_S.first)
826  info.max_elem_det_S = std::make_pair(det_S_int, elem->id());
827  else if (det_S_int < info.min_elem_det_S.first)
828  info.min_elem_det_S = std::make_pair(det_S_int, elem->id());
829 
830  info.total_distortion += beta_int;
831  if (beta_int > info.max_elem_distortion.first)
832  info.max_elem_distortion = std::make_pair(beta_int, elem->id());
833  else if (beta_int < info.min_elem_distortion.first)
834  info.min_elem_distortion = std::make_pair(beta_int, elem->id());
835 
836  info.total_dilation += mu_int;
837  if (mu_int > info.max_elem_dilation.first)
838  info.max_elem_dilation = std::make_pair(mu_int, elem->id());
839  else if (mu_int < info.min_elem_dilation.first)
840  info.min_elem_dilation = std::make_pair(mu_int, elem->id());
841 
842  info.total_combined += combined_int;
843  if (combined_int > info.max_elem_combined.first)
844  info.max_elem_combined = std::make_pair(combined_int, elem->id());
845  else if (combined_int < info.min_elem_combined.first)
846  info.min_elem_combined = std::make_pair(combined_int, elem->id());
847 
848  } // for elem
849 
850  // Get contributions from elements on other processors
851  communicate_pair_max(info.max_elem_det_S, mesh.comm());
852  communicate_pair_min(info.min_elem_det_S, mesh.comm());
853  mesh.comm().max(info.max_qp_det_S);
854  mesh.comm().min(info.min_qp_det_S);
855  mesh.comm().sum(info.total_det_S);
856 
857  communicate_pair_max(info.max_elem_distortion, mesh.comm());
858  communicate_pair_min(info.min_elem_distortion, mesh.comm());
859  mesh.comm().sum(info.total_distortion);
860 
861  communicate_pair_max(info.max_elem_dilation, mesh.comm());
862  communicate_pair_min(info.min_elem_dilation, mesh.comm());
863  mesh.comm().sum(info.total_dilation);
864 
865  communicate_pair_max(info.max_elem_combined, mesh.comm());
866  communicate_pair_min(info.min_elem_combined, mesh.comm());
867  mesh.comm().sum(info.total_combined);
868 
869  mesh.comm().max(info.mesh_is_tangled);
870 
871  _mesh_info = info;
872 }
873 
874 std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
876 {
877  // Build target element
878  auto target_elem = Elem::build(type);
879 
880  // Volume of reference element
881  const auto ref_vol = target_elem->reference_elem()->volume();
882 
883  // Update the nodes of the target element, depending on type
884  const Real sqrt_3 = std::sqrt(Real(3));
885  std::vector<std::unique_ptr<Node>> owned_nodes;
886 
887  const auto type_str = Utility::enum_to_string(type);
888 
889  // Elems deriving from Tri
890  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  const auto side_length = std::sqrt(4. / sqrt_3 * ref_vol);
898 
899  // Define the nodal locations of the vertices
900  const auto & s = side_length;
901  // x y node_id
902  owned_nodes.emplace_back(Node::build(Point(0., 0.), 0));
903  owned_nodes.emplace_back(Node::build(Point(s, 0.), 1));
904  owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s), 2));
905 
906  switch (type)
907  {
908  case TRI3: {
909  // Nothing to do here, vertices already added above
910  break;
911  }
912 
913  case TRI6: {
914  // Define the midpoint nodes of the equilateral triangle
915  // x y node_id
916  owned_nodes.emplace_back(Node::build(Point(0.50 * s, 0.00), 3));
917  owned_nodes.emplace_back(Node::build(Point(0.75 * s, 0.25 * sqrt_3 * s), 4));
918  owned_nodes.emplace_back(Node::build(Point(0.25 * s, 0.25 * sqrt_3 * s), 5));
919 
920  break;
921  }
922 
923  default:
924  libmesh_error_msg("Unsupported triangular element: " << type_str);
925  break;
926  }
927  } // if Tri
928 
929 
930  // Elems deriving from Pyramid
931  else if (type_str.compare(0, 7, "PYRAMID") == 0)
932  {
933 
934  // The target element is a pyramid with an square base and
935  // equilateral triangular sides with volume equal to the volume of the
936  // reference element.
937 
938  // A pyramid with square base sidelength s and equilateral triangular
939  // sides has height h = s / sqrt(2).
940  // The volume is v = s^2 h / 3 = s^3 / ( 3 sqrt(2)).
941  // Solving for s: s = (3 sqrt(2) v)^(1/3), where v is the volume of the
942  // non-optimal reference element.
943 
944  const Real sqrt_2 = std::sqrt(Real(2));
945  // Side length that preserves the volume of the reference element
946  const auto side_length = std::pow(3. * sqrt_2 * ref_vol, 1. / 3.);
947  // Pyramid height with the property that all faces are equilateral triangles
948  const auto target_height = side_length / sqrt_2;
949 
950  const auto & s = side_length;
951  const auto & h = target_height;
952 
953  // x y z node_id
954  owned_nodes.emplace_back(Node::build(Point(0., 0., 0.), 0));
955  owned_nodes.emplace_back(Node::build(Point(s, 0., 0.), 1));
956  owned_nodes.emplace_back(Node::build(Point(s, s, 0.), 2));
957  owned_nodes.emplace_back(Node::build(Point(0., s, 0.), 3));
958  owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * s, h), 4));
959 
960  if (type == PYRAMID13 || type == PYRAMID14 || type == PYRAMID18)
961  {
962  const auto & on = owned_nodes;
963  // Define the edge midpoint nodes of the pyramid
964 
965  // Base node to base node midpoint nodes
966  owned_nodes.emplace_back(Node::build((*on[0] + *on[1]) / 2., 5));
967  owned_nodes.emplace_back(Node::build((*on[1] + *on[2]) / 2., 6));
968  owned_nodes.emplace_back(Node::build((*on[2] + *on[3]) / 2., 7));
969  owned_nodes.emplace_back(Node::build((*on[3] + *on[0]) / 2., 8));
970 
971  // Base node to apex node midpoint nodes
972  owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[4]) / 2.), 9));
973  owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
974  owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[4]) / 2.), 11));
975  owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
976 
977  if (type == PYRAMID14 || type == PYRAMID18)
978  {
979  // Define the square face midpoint node of the pyramid
980  owned_nodes.emplace_back(
981  Node::build(Point((*on[0] + *on[1] + *on[2] + *on[3]) / 4.), 13));
982 
983  if (type == PYRAMID18)
984  {
985  // Define the triangular face nodes
986  owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[4]) / 3.), 14));
987  owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4]) / 3.), 15));
988  owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3] + *on[4]) / 3.), 16));
989  owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[0] + *on[4]) / 3.), 17));
990  }
991  }
992  }
993 
994  else if (type != PYRAMID5)
995  libmesh_error_msg("Unsupported pyramid element: " << type_str);
996 
997  } // if Pyramid
998 
999  // Set the target_elem equal to the reference elem
1000  else
1001  for (const auto & node : target_elem->reference_elem()->node_ref_range())
1002  owned_nodes.emplace_back(Node::build(node, node.id()));
1003 
1004  // Set nodes of target element
1005  for (const auto & node_ptr : owned_nodes)
1006  target_elem->set_node(node_ptr->id(), node_ptr.get());
1007 
1008  libmesh_assert(relative_fuzzy_equals(target_elem->volume(), ref_vol, TOLERANCE));
1009 
1010  return std::make_pair(std::move(target_elem), std::move(owned_nodes));
1011 }
1012 
1014  const Elem * const target_elem,
1015  const FEMContext & femcontext,
1016  std::vector<RealTensor> & jacobians,
1017  std::vector<Real> & jacobian_dets)
1018 {
1019 
1020  const auto dim = target_elem->dim();
1021 
1022  const auto & qrule_points = femcontext.get_element_qrule().get_points();
1023  const auto & qrule_weights = femcontext.get_element_qrule().get_weights();
1024  const auto nq_points = femcontext.get_element_qrule().n_points();
1025 
1026  // If the target element is the reference element, Jacobian matrix is
1027  // identity, det of inverse is 1. These will only be overwritten if a
1028  // different target element is explicitly specified.
1029  jacobians = std::vector<RealTensor>(nq_points, RealTensor(
1030  1., 0., 0.,
1031  0., 1., 0.,
1032  0., 0., 1.));
1033  jacobian_dets = std::vector<Real>(nq_points, 1.0);
1034 
1035  // Don't use "if (*target_elem == *(target_elem->reference_elem()))" here, it
1036  // only compares global node ids, not the node locations themselves.
1037  bool target_equals_reference = true;
1038  const auto * ref_elem = target_elem->reference_elem();
1039  for (const auto local_id : make_range(target_elem->n_nodes()))
1040  target_equals_reference &= target_elem->node_ref(local_id) == ref_elem->node_ref(local_id);
1041  if (target_equals_reference)
1042  return;
1043 
1044  // Create FEMap to compute target_element mapping information
1045  FEMap fe_map_target;
1046 
1047  // pre-request mapping derivatives
1048  fe_map_target.get_dxyzdxi();
1049  fe_map_target.get_dxyzdeta();
1050  fe_map_target.get_dxyzdzeta();
1051 
1052  // build map
1053  fe_map_target.init_reference_to_physical_map(dim, qrule_points, target_elem);
1054  fe_map_target.compute_map(dim, qrule_weights, target_elem, /*d2phi=*/false);
1055 
1056  for (const auto qp : make_range(nq_points))
1057  {
1058  // We use Larisa's H notation to denote the reference-to-target jacobian
1059  RealTensor H = get_jacobian_at_qp(fe_map_target, dim, qp);
1060 
1061  // The target-to-reference jacobian is the inverse of the
1062  // reference-to-target jacobian
1063  jacobians[qp] = H.inverse();
1064  jacobian_dets[qp] = jacobians[qp].det();
1065  }
1066 }
1067 
1068 } // 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...
virtual void solve() override
Solves the system to smooth the mesh.
ElemType
Defines an enum for geometric element types.
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
MPI_Info info
void minloc(T &r, unsigned int &min_id) const
static constexpr Real TOLERANCE
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
Real trace(const RealTensor &A, const unsigned int &dim)
Compute the trace of a dim-dimensional matrix.
Real _epsilon_squared_assembly
Epsilon squared value determined at runtime during each assembly.
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
Real chi_epsilon(const Real &x, const Real epsilon_squared)
Function to prevent dividing by zero for degenerate elements.
const std::vector< Real > & get_weights() const
Definition: quadrature.h:168
NumericVector< Number > * rhs
The system matrix.
RealTensorValue RealTensor
static std::pair< std::unique_ptr< Elem >, std::vector< std::unique_ptr< Node > > > get_target_elem(const ElemType &type)
Get the target element for a given element type.
static void get_target_to_reference_jacobian(const Elem *const target_elem, const FEMContext &femcontext, std::vector< RealTensor > &jacobians, std::vector< Real > &jacobian_dets)
Get the jacobians (and determinants) of the target-to-reference element mapping.
const std::vector< RealGradient > & get_dxyzdzeta() const
Definition: fe_map.h:255
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.
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
Struct to hold smoother-relevant information about the mesh quality.
unsigned int number() const
Definition: system.h:2350
T pow(const T &x)
Definition: utility.h:328
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2529
virtual unsigned int n_nodes() const =0
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
const MeshQualityInfo & get_mesh_info()
Getter for the _mesh_info attribute.
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
libmesh_assert(ctx)
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:62
bool _untangling_solve
Flag to indicate if the current solve is to untangle or smooth.
const Elem * reference_elem() const
Definition: elem.C:570
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 maxloc(T &r, unsigned int &max_id) const
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.
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
void communicate_pair_min(std::pair< Real, dof_id_type > &pair, const Parallel::Communicator &comm)
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
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 meshes only)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< ElemType, std::vector< Real > > _target_jacobian_dets
const std::vector< Point > & get_points() const
Definition: quadrature.h:156
virtual unsigned short dim() const =0
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
void communicate_pair_max(std::pair< Real, dof_id_type > &pair, const Parallel::Communicator &comm)
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
virtual void solve() override
Invokes the solver associated with the system.
Definition: fem_system.C:1076
const std::vector< RealGradient > & get_dxyzdeta() const
Definition: fe_map.h:247
bool relative_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
Function to check whether two variables are equal within a relative tolerance.
Definition: fuzzy_equals.h:78
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
std::map< ElemType, std::vector< RealTensor > > _target_jacobians
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
RealTensor get_jacobian_at_qp(const FEMap &fe_map, const unsigned int &dim, const unsigned int &qp)
Given an fe_map, element dimension, and quadrature point index, returns the Jacobian of the physical-...
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
MeshQualityInfo _mesh_info
Information about the mesh quality.