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" 44 bool apply_heterogeneous_constraints,
45 bool apply_no_constraints)
51 for (
auto * node :
mesh.local_node_ptr_range())
55 const auto dof_id = node->dof_number(this->
number(), d, 0);
64 FEMSystem::assembly(get_residual, get_jacobian, apply_heterogeneous_constraints, apply_no_constraints);
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();
85 for (
auto * node :
mesh.local_node_ptr_range())
89 const auto dof_id = node->dof_number(this->
number(), d, 0);
91 (*solution).set(dof_id, (*node)(d));
101 FEMContext & femcontext = cast_ref<FEMContext &>(*con);
105 const auto dim =
mesh.mesh_dimension();
107 Real elem_averaged_det_J_sum = 0.;
113 const auto & JxW = fe_map.get_JxW();
115 std::map<ElemType, std::vector<Real>> target_elem_inverse_jacobian_dets;
117 for (
const auto * elem :
mesh.active_local_element_ptr_range())
129 const auto & dxyzdxi = fe_map_target.
get_dxyzdxi();
140 target_elem_inverse_jacobian_dets[elem->type()] =
141 std::vector<Real>(nq_points, 1.0);
146 if (type_str.compare(0, 3,
"TRI") == 0)
153 const Real sqrt_3 = std::sqrt(3.);
155 std::vector<Point> equilateral_points{
158 Point(0.5, 0.5 * sqrt_3)
162 const auto target_elem =
Elem::build(elem->type());
165 const auto ref_area = target_elem->reference_elem()->volume();
167 switch (elem->type())
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);
186 libmesh_error_msg(
"Unsupported triangular element: " << type_str);
191 const auto side_length = std::sqrt(4. / sqrt_3 * ref_area);
193 std::vector<std::unique_ptr<Node>> owned_nodes;
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());
206 fe_map_target.
compute_map(
dim, qrule_weights, target_elem.get(),
false);
211 std::vector<RealTensor>(nq_points);
215 RealTensor(dxyzdxi[qp](0), dxyzdeta[qp](0), 0, dxyzdxi[qp](1),
216 dxyzdeta[qp](1), 0, 0, 0, 1)
220 target_elem_inverse_jacobian_dets[target_elem->type()][qp] =
227 Real elem_integrated_det_J = 0.;
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;
237 mesh.comm().sum(elem_averaged_det_J_sum);
239 _ref_vol = elem_averaged_det_J_sum /
mesh.n_active_elem();
244 FEMContext & c = cast_ref<FEMContext &>(context);
252 const std::set<unsigned char> & elem_dims =
255 for (
const auto &
dim : elem_dims)
258 my_fe->get_nothing();
260 auto & fe_map = my_fe->get_fe_map();
261 fe_map.get_dxyzdxi();
262 fe_map.get_dxyzdeta();
263 fe_map.get_dxyzdzeta();
267 my_fe->get_nothing();
277 FEMContext & c = cast_ref<FEMContext &>(context);
286 const auto & dxyzdxi = fe_map.get_dxyzdxi();
287 const auto & dxyzdeta = fe_map.get_dxyzdeta();
288 const auto & dxyzdzeta = fe_map.get_dxyzdzeta();
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();
296 unsigned int x_var = 0, y_var = 1, z_var = 2;
307 std::reference_wrapper<DenseSubVector<Number>> F[3] =
314 std::reference_wrapper<DenseSubMatrix<Number>> K[3][3] =
354 dxyzdxi[qp](0), 0, 0,
362 dxyzdxi[qp](0), dxyzdeta[qp](0), 0,
363 dxyzdxi[qp](1), dxyzdeta[qp](1), 0,
377 libmesh_error_msg(
"Unsupported dimension.");
389 const Real det_sq = det * det;
390 const Real det_cube = det_sq * det;
402 const Real tr_div_dim = tr /
dim;
405 const Real half_dim = 0.5 *
dim;
406 const std::vector<Real> trace_powers{
408 std::pow(tr_div_dim, half_dim - 1.),
409 std::pow(tr_div_dim, half_dim - 2.),
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);
430 const Real chi_prime = 0.5 * (1. + det / sqrt_term);
431 const Real chi_prime_sq = chi_prime * chi_prime;
433 const Real chi_2prime = 0.5 * (1. / sqrt_term - det_sq / Utility::pow<3>(sqrt_term));
437 const RealTensor dbeta_dS = (trace_powers[1] / chi) * S - (trace_powers[0] / chi_sq * chi_prime * det) * S_inv_T;
442 const Real alpha = (-chi_prime * det_cube + 2. * det_sq * chi - ref_vol_sq * det * chi_prime) / (2. *
_ref_vol * chi_sq);
450 std::vector<std::vector<std::vector<Real>>> dphi_maps = {dphidxi_map, dphideta_map, dphidzeta_map};
454 for (
const auto l : elem.node_index_range())
461 dS_dR(var_id, jj) = dphi_maps[jj][l][qp];
465 F[var_id](l) += quad_weights[qp] * dE_dS.
contract(dS_dR);
470 if (request_jacobian)
484 const std::vector<Real> d2beta_dS2_coefs_times_distortion_weight = {
488 (trace_powers[1] / chi) * distortion_weight,
490 (((
dim - 2.) /
dim) * trace_powers[2] / chi) * distortion_weight,
492 (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
497 (-(trace_powers[1] / chi_sq) * chi_prime * det) * distortion_weight,
499 (trace_powers[0] * (det / chi_sq)
500 * ((2. * chi_prime_sq / chi - chi_2prime) * det - chi_prime)) * distortion_weight,
502 ((trace_powers[0] / chi_sq) * chi_prime * det) * distortion_weight,
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;
604 for (
const auto l: elem.node_index_range())
612 dS_dR_l(var_id1, ii) = dphi_maps[ii][l][qp];
623 dS_dR_p(var_id2, jj) = dphi_maps[jj][p][qp];
633 const auto S_ij = S(i,j);
634 const auto S_inv_ji = S_inv(j,i);
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,
645 const auto dalpha_dS_coef_ij_applied = dalpha_dS_coef_times_dilation_weight * S_inv_ji;
647 Real d2E_dSdR_l = 0.;
648 Real d2E_dSdR_p = 0.;
657 if (!(a == var_id1 && i == var_id2) &&
658 !(a == var_id2 && i == var_id1))
661 const auto S_inv_ja = S_inv(j,a);
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;
667 const auto b_limit = (a == i) ? j + 1 :
dim;
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)
684 const Real d2mu_dS2_times_dilation_weight = dalpha_dS_coef_ij_applied * S_inv(b,a) - alpha_ja_applied * S_inv(b,i);
689 d2beta_dS2_times_distortion_weight +
690 d2mu_dS2_times_dilation_weight;
694 d2E_dSdR_l += d2E_dS2 * dS_dR_l(a, b);
696 if (!(i == a && j == b))
699 d2E_dSdR_p += d2E_dS2 * dS_dR_p(a, b);
704 d2E_dSdR_l * dS_dR_p(i, j) + d2E_dSdR_p * dS_dR_l(i, j);
709 const Real jacobian_contribution = quad_weights[qp] * d2E_dR2;
710 K[var_id1][var_id2](l, p) += jacobian_contribution;
715 K[var_id2][var_id1](p, l) += jacobian_contribution;
724 return request_jacobian;
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...
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.
const Elem & get_elem() const
Accessor for current Elem object.
virtual void init_context(DiffContext &) override
TypeTensor< T > transpose() const
std::map< ElemType, std::vector< RealTensor > > _target_inverse_jacobians
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.
const std::vector< Real > & get_weights() const
RealTensorValue RealTensor
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. The distortion metric is given weight 1 - _dilation_...
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.
unsigned int number() const
const std::vector< RealGradient > & get_dxyzdxi() const
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
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.
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 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.
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 elements only) ...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< Point > & get_points() const
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...
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...
const std::vector< RealGradient > & get_dxyzdeta() const
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.
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
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.