18 #include "libmesh/variational_smoother_system.h" 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> 48 comm.
minloc(pair.first, rank);
59 comm.
maxloc(pair.first, rank);
68 return 0.5 * (x + std::sqrt(epsilon_squared + Utility::pow<2>(x)));
76 const unsigned int &
dim,
77 const unsigned int & qp)
98 dxyzdxi(0), dxyzdeta(0), 0,
99 dxyzdxi(1), dxyzdeta(1), 0,
113 libmesh_error_msg(
"Unsupported dimension.");
133 bool apply_heterogeneous_constraints,
134 bool apply_no_constraints)
140 for (
auto * node :
mesh.local_node_ptr_range())
144 const auto dof_id = node->dof_number(this->
number(), d, 0);
163 const Real variable_component = 100. * Utility::pow<2>(
_ref_vol * min_S);
170 FEMSystem::assembly(get_residual, get_jacobian, apply_heterogeneous_constraints, apply_no_constraints);
187 if (mesh_info.mesh_is_tangled)
209 libmesh_error_msg_if(
get_mesh_info().mesh_is_tangled,
"The smoothing solve tangled the 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();
230 for (
auto * node :
mesh.local_node_ptr_range())
234 const auto dof_id = node->dof_number(this->
number(), d, 0);
236 (*solution).set(dof_id, (*node)(d));
250 FEMContext & femcontext = cast_ref<FEMContext &>(*con);
255 Real elem_averaged_det_S_sum = 0.;
261 const auto & JxW = fe_map.get_JxW();
263 for (
const auto * elem :
mesh.active_local_element_ptr_range())
271 const auto [target_elem, target_nodes] =
get_target_elem(elem->type());
279 Real elem_integrated_det_S = 0.;
282 const auto ref_elem_vol = elem->reference_elem()->volume();
283 elem_averaged_det_S_sum += elem_integrated_det_S / ref_elem_vol;
288 mesh.comm().sum(elem_averaged_det_S_sum);
290 _ref_vol = elem_averaged_det_S_sum /
mesh.n_active_elem();
295 FEMContext & c = cast_ref<FEMContext &>(context);
303 const std::set<unsigned char> & elem_dims =
306 for (
const auto &
dim : elem_dims)
309 my_fe->get_nothing();
311 auto & fe_map = my_fe->get_fe_map();
312 fe_map.get_dxyzdxi();
313 fe_map.get_dxyzdeta();
314 fe_map.get_dxyzdzeta();
318 fe_map.set_jacobian_tolerance(std::numeric_limits<Real>::lowest());
321 my_fe->get_nothing();
331 FEMContext & c = cast_ref<FEMContext &>(context);
337 unsigned int x_var = 0, y_var = 1, z_var = 2;
348 std::reference_wrapper<DenseSubVector<Number>> F[3] =
355 std::reference_wrapper<DenseSubMatrix<Number>> K[3][3] =
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();
405 const Real det_sq = det * det;
406 const Real det_cube = det_sq * det;
415 const Real tr_div_dim = tr /
dim;
418 const Real half_dim = 0.5 *
dim;
419 const std::vector<Real> trace_powers{
421 std::pow(tr_div_dim, half_dim - 1.),
422 std::pow(tr_div_dim, half_dim - 2.),
437 const Real chi_sq = chi * chi;
440 const Real chi_prime = 0.5 * (1. + det / sqrt_term);
441 const Real chi_prime_sq = chi_prime * chi_prime;
443 const Real chi_2prime = 0.5 * (1. / sqrt_term - det_sq / Utility::pow<3>(sqrt_term));
447 const RealTensor dbeta_dS = (trace_powers[1] / chi) * S - (trace_powers[0] / chi_sq * chi_prime * det) * S_inv_T;
452 const Real alpha = (-chi_prime * det_cube + 2. * det_sq * chi - ref_vol_sq * det * chi_prime) / (2. *
_ref_vol * chi_sq);
460 std::vector<std::vector<std::vector<Real>>> dphi_maps = {dphidxi_map, dphideta_map, dphidzeta_map};
464 for (
const auto l : elem.node_index_range())
471 dS_dR(var_id, jj) = dphi_maps[jj][l][qp];
475 F[var_id](l) += quad_weights[qp] * dE_dS.
contract(dS_dR);
480 if (request_jacobian)
494 const std::vector<Real> d2beta_dS2_coefs_times_distortion_weight = {
498 (trace_powers[1] / chi) * distortion_weight,
500 (((
dim - 2.) /
dim) * trace_powers[2] / chi) * distortion_weight,
502 (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
507 (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
509 (trace_powers[0] * (det / chi_sq)
510 * ((2. * chi_prime_sq / chi - chi_2prime) * det - chi_prime)) * distortion_weight,
512 ((trace_powers[0] / chi_sq) * chi_prime * det) * distortion_weight,
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;
614 for (
const auto l: elem.node_index_range())
622 dS_dR_l(var_id1, ii) = dphi_maps[ii][l][qp];
633 dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
643 const auto S_ij = S(i,j);
644 const auto S_inv_ji = S_inv(j,i);
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,
655 const auto dalpha_dS_coef_ij_applied = dalpha_dS_coef_times_dilation_weight * S_inv_ji;
657 Real d2E_dSdR_l = 0.;
658 Real d2E_dSdR_p = 0.;
667 if (!(a == var_id1 && i == var_id2) &&
668 !(a == var_id2 && i == var_id1))
671 const auto S_inv_ja = S_inv(j,a);
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;
677 const auto b_limit = (a == i) ? j + 1 :
dim;
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)
694 const Real d2mu_dS2_times_dilation_weight = dalpha_dS_coef_ij_applied * S_inv(b,a) - alpha_ja_applied * S_inv(b,i);
699 d2beta_dS2_times_distortion_weight +
700 d2mu_dS2_times_dilation_weight;
704 d2E_dSdR_l += d2E_dS2 * dS_dR_l(a, b);
706 if (!(i == a && j == b))
709 d2E_dSdR_p += d2E_dS2 * dS_dR_p(a, b);
714 d2E_dSdR_l * dS_dR_p(i, j) + d2E_dSdR_p * dS_dR_l(i, j);
719 const Real jacobian_contribution = quad_weights[qp] * d2E_dR2;
720 K[var_id1][var_id2](l, p) += jacobian_contribution;
725 K[var_id2][var_id1](p, l) += jacobian_contribution;
734 return request_jacobian;
752 FEMContext & femcontext = cast_ref<FEMContext &>(*con);
756 const auto dim =
mesh.mesh_dimension();
757 const Real half_dim = 0.5 *
dim;
764 const auto & JxW = fe_map.get_JxW();
765 fe_map.get_dxyzdxi();
766 fe_map.get_dxyzdeta();
767 fe_map.get_dxyzdzeta();
771 for (
const auto * elem :
mesh.active_local_element_ptr_range())
780 Real combined_int = 0.;
785 det_S_int += det_SxW;
794 const auto det = S.
det();
795 const auto det_sq = det * det;
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;
803 info.mesh_is_tangled =
true;
813 beta_int += beta * det_SxW;
817 mu_int += mu * det_SxW;
821 combined_int += E * det_SxW;
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());
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());
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());
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());
859 mesh.comm().sum(
info.total_distortion);
863 mesh.comm().sum(
info.total_dilation);
867 mesh.comm().sum(
info.total_combined);
869 mesh.comm().max(
info.mesh_is_tangled);
874 std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
881 const auto ref_vol = target_elem->reference_elem()->volume();
884 const Real sqrt_3 = std::sqrt(
Real(3));
885 std::vector<std::unique_ptr<Node>> owned_nodes;
890 if (type_str.compare(0, 3,
"TRI") == 0)
897 const auto side_length = std::sqrt(4. / sqrt_3 * ref_vol);
900 const auto & s = side_length;
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));
924 libmesh_error_msg(
"Unsupported triangular element: " << type_str);
930 else if (type_str.compare(0, 5,
"PRISM") == 0)
948 const auto side_length =
std::pow(16. * ref_vol / 3., 1. / 3.);
950 const auto target_height = 0.25 * side_length * sqrt_3;
952 const auto & s = side_length;
953 const auto & h = target_height;
957 owned_nodes.emplace_back(
Node::build(
Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
960 owned_nodes.emplace_back(
Node::build(
Point(0.5 * s, 0.5 * sqrt_3 * s, h), 5));
965 const auto & on = owned_nodes;
979 owned_nodes.emplace_back(
Node::build(
Point((*on[0] + *on[1] + *on[3] + *on[4]) / 4.), 15));
980 owned_nodes.emplace_back(
Node::build(
Point((*on[1] + *on[2] + *on[4] + *on[5]) / 4.), 16));
981 owned_nodes.emplace_back(
Node::build(
Point((*on[0] + *on[2] + *on[3] + *on[5]) / 4.), 17));
986 owned_nodes.emplace_back(
Node::build(
Point((*on[0] + *on[1] + *on[2]) / 3.), 18));
987 owned_nodes.emplace_back(
Node::build(
Point((*on[3] + *on[4] + *on[5]) / 3.), 19));
991 owned_nodes.emplace_back(
Node::build(
Point((*on[9] + *on[10] + *on[11]) / 3.), 20));
998 libmesh_error_msg(
"Unsupported prism element: " << type_str);
1003 else if (type_str.compare(0, 7,
"PYRAMID") == 0)
1016 const Real sqrt_2 = std::sqrt(
Real(2));
1018 const auto side_length =
std::pow(3. * sqrt_2 * ref_vol, 1. / 3.);
1020 const auto target_height = side_length / sqrt_2;
1022 const auto & s = side_length;
1023 const auto & h = target_height;
1034 const auto & on = owned_nodes;
1038 owned_nodes.emplace_back(
Node::build((*on[0] + *on[1]) / 2., 5));
1039 owned_nodes.emplace_back(
Node::build((*on[1] + *on[2]) / 2., 6));
1040 owned_nodes.emplace_back(
Node::build((*on[2] + *on[3]) / 2., 7));
1041 owned_nodes.emplace_back(
Node::build((*on[3] + *on[0]) / 2., 8));
1052 owned_nodes.emplace_back(
1058 owned_nodes.emplace_back(
Node::build(
Point((*on[0] + *on[1] + *on[4]) / 3.), 14));
1059 owned_nodes.emplace_back(
Node::build(
Point((*on[1] + *on[2] + *on[4]) / 3.), 15));
1060 owned_nodes.emplace_back(
Node::build(
Point((*on[2] + *on[3] + *on[4]) / 3.), 16));
1061 owned_nodes.emplace_back(
Node::build(
Point((*on[3] + *on[0] + *on[4]) / 3.), 17));
1067 libmesh_error_msg(
"Unsupported pyramid element: " << type_str);
1073 for (
const auto & node : target_elem->reference_elem()->node_ref_range())
1074 owned_nodes.emplace_back(
Node::build(node, node.id()));
1077 for (
const auto & node_ptr : owned_nodes)
1078 target_elem->set_node(node_ptr->id(), node_ptr.get());
1082 return std::make_pair(std::move(target_elem), std::move(owned_nodes));
1086 const Elem *
const target_elem,
1088 std::vector<RealTensor> & jacobians,
1089 std::vector<Real> & jacobian_dets)
1092 const auto dim = target_elem->
dim();
1101 jacobians = std::vector<RealTensor>(nq_points,
RealTensor(
1105 jacobian_dets = std::vector<Real>(nq_points, 1.0);
1109 bool target_equals_reference =
true;
1112 target_equals_reference &= target_elem->
node_ref(local_id) == ref_elem->node_ref(local_id);
1113 if (target_equals_reference)
1117 FEMap fe_map_target;
1126 fe_map_target.
compute_map(
dim, qrule_weights, target_elem,
false);
1136 jacobian_dets[qp] = jacobians[qp].det();
virtual std::unique_ptr< DiffContext > build_context() override
Builds a FEMContext object with enough information to do evaluations on each element.
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.
This class provides all data required for a physics package (e.g.
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...
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
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.
void minloc(T &r, unsigned int &min_id) const
static constexpr Real TOLERANCE
const Elem & get_elem() const
Accessor for current Elem object.
virtual void init_context(DiffContext &) override
TypeTensor< T > transpose() const
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.
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.
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
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
The libMesh namespace provides an interface to certain functionality in the library.
const MeshBase & get_mesh() const
Real _dilation_weight
The relative weight to give the dilation metric.
unsigned char get_dim() const
Accessor for cached mesh dimension.
virtual void elem_fe_reinit(const std::vector< Point > *const pts=nullptr)
Reinitializes interior FE objects on the current geometric element.
Struct to hold smoother-relevant information about the mesh quality.
unsigned int number() const
const Node & node_ref(const unsigned int i) const
virtual unsigned int n_nodes() const =0
const std::vector< RealGradient > & get_dxyzdxi() const
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
const MeshQualityInfo & get_mesh_info()
Getter for the _mesh_info attribute.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
virtual void init_data() override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
This class provides all data required for a physics package (e.g.
bool _untangling_solve
Flag to indicate if the current solve is to untangle or smooth.
const Elem * reference_elem() const
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.
unsigned int n_points() const
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'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.
std::string enum_to_string(const T e)
void init_reference_to_physical_map(const std::vector< Point > &qp, const Elem *elem)
Real _ref_vol
The reference volume for each element.
static std::unique_ptr< Node > build(const Node &n)
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
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
virtual unsigned short dim() const =0
TypeTensor< T > inverse() const
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...
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.
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...
virtual void solve() override
Invokes the solver associated with the system.
const std::vector< RealGradient > & get_dxyzdeta() const
void compute_mesh_quality_info()
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.
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.
std::map< ElemType, std::vector< RealTensor > > _target_jacobians
Class contained in FE that encapsulates mapping (i.e.
A Point defines a location in LIBMESH_DIM dimensional Real space.
~VariationalSmootherSystem() override
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
void prepare_for_smoothing()
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.
const std::set< unsigned char > & elem_dimensions() const
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.