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);
931 else if (type_str.compare(0, 7,
"PYRAMID") == 0)
944 const Real sqrt_2 = std::sqrt(
Real(2));
946 const auto side_length =
std::pow(3. * sqrt_2 * ref_vol, 1. / 3.);
948 const auto target_height = side_length / sqrt_2;
950 const auto & s = side_length;
951 const auto & h = target_height;
962 const auto & on = owned_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));
980 owned_nodes.emplace_back(
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));
995 libmesh_error_msg(
"Unsupported pyramid element: " << type_str);
1001 for (
const auto & node : target_elem->reference_elem()->node_ref_range())
1002 owned_nodes.emplace_back(
Node::build(node, node.id()));
1005 for (
const auto & node_ptr : owned_nodes)
1006 target_elem->set_node(node_ptr->id(), node_ptr.get());
1010 return std::make_pair(std::move(target_elem), std::move(owned_nodes));
1014 const Elem *
const target_elem,
1016 std::vector<RealTensor> & jacobians,
1017 std::vector<Real> & jacobian_dets)
1020 const auto dim = target_elem->
dim();
1029 jacobians = std::vector<RealTensor>(nq_points,
RealTensor(
1033 jacobian_dets = std::vector<Real>(nq_points, 1.0);
1037 bool target_equals_reference =
true;
1040 target_equals_reference &= target_elem->
node_ref(local_id) == ref_elem->node_ref(local_id);
1041 if (target_equals_reference)
1045 FEMap fe_map_target;
1054 fe_map_target.
compute_map(
dim, qrule_weights, target_elem,
false);
1064 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.