13 #include "libmesh/enum_to_string.h" 14 #include "metaphysicl/dualnumberarray.h" 15 #include "Eigen/Dense" 17 using MetaPhysicL::NumberArray;
19 typedef DualNumber<Real, NumberArray<2, Real>>
Dual2;
27 const Elem *
const primal_elem,
28 const unsigned int sub_elem_index,
29 const QBase & qrule_msm,
30 std::vector<Point> & q_pts)
32 const auto msm_elem_order = msm_elem->default_order();
33 const auto msm_elem_type = msm_elem->type();
36 const Point e1 = msm_elem->point(0) - msm_elem->point(1);
37 const Point e2 = msm_elem->point(2) - msm_elem->point(1);
38 const Point normal = e2.cross(e1).unit();
41 const auto sub_elem = msm_elem->get_extra_integer(sub_elem_index);
42 const ElemType primal_type = primal_elem->type();
44 auto get_sub_elem_node_indices = [primal_type, sub_elem]() -> std::vector<unsigned int>
51 return {{0, 1, 2, 3}};
65 mooseError(
"get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
79 return {{4, 5, 6, 7}};
81 mooseError(
"get_sub_elem_nodes: Invalid sub_elem: ", sub_elem);
87 return {{0, 4, 8, 7}};
89 return {{4, 1, 5, 8}};
91 return {{8, 5, 2, 6}};
93 return {{7, 8, 6, 3}};
95 mooseError(
"get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
98 mooseError(
"get_sub_elem_indices: Face element type: ",
99 libMesh::Utility::enum_to_string<ElemType>(primal_type),
100 " invalid for 3D mortar");
106 auto transform_qp = [primal_type, sub_elem](
const Real nu,
const Real xi)
111 return Point(nu, xi, 0);
113 return Point(nu, xi, 0);
119 return Point(0.5 * nu, 0.5 * xi, 0);
121 return Point(0.5 * (1 - xi), 0.5 * (nu + xi), 0);
123 return Point(0.5 * (1 + nu), 0.5 * xi, 0);
125 return Point(0.5 * nu, 0.5 * (1 + xi), 0);
127 mooseError(
"get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
133 return Point(nu - 1, xi - 1, 0);
135 return Point(nu + xi, xi - 1, 0);
137 return Point(1 - xi, nu + xi, 0);
139 return Point(nu - 1, nu + xi, 0);
141 return Point(0.5 * (nu - xi), 0.5 * (nu + xi), 0);
143 mooseError(
"get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
149 return Point(0.5 * (nu - 1), 0.5 * (xi - 1), 0);
151 return Point(0.5 * (nu + 1), 0.5 * (xi - 1), 0);
153 return Point(0.5 * (nu + 1), 0.5 * (xi + 1), 0);
155 return Point(0.5 * (nu - 1), 0.5 * (xi + 1), 0);
157 mooseError(
"get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
160 mooseError(
"transform_qp: Face element type: ",
161 libMesh::Utility::enum_to_string<ElemType>(primal_type),
162 " invalid for 3D mortar");
166 auto sub_element_type = [primal_type, sub_elem]()
188 mooseError(
"sub_element_type: Invalid sub_elem: ", sub_elem);
191 mooseError(
"sub_element_type: Face element type: ",
192 libMesh::Utility::enum_to_string<ElemType>(primal_type),
193 " invalid for 3D mortar");
198 auto sub_elem_node_indices = get_sub_elem_node_indices();
201 for (
auto qp :
make_range(qrule_msm.n_points()))
205 for (
auto n :
make_range(msm_elem->n_nodes()))
215 xi1.value() = qrule_msm.qp(qp)(0);
216 xi1.derivatives()[0] = 1.0;
218 xi2.value() = qrule_msm.qp(qp)(1);
219 xi2.derivatives()[1] = 1.0;
220 VectorValue<Dual2> xi(xi1, xi2, 0);
221 unsigned int current_iterate = 0, max_iterates = 10;
227 VectorValue<Dual2> x1;
228 for (
auto n :
make_range(sub_elem_node_indices.size()))
230 primal_elem->point(sub_elem_node_indices[n]);
233 VectorValue<Dual2> F(u(1) * normal(2) - u(2) * normal(1),
234 u(2) * normal(0) - u(0) * normal(2),
235 u(0) * normal(1) - u(1) * normal(0));
237 Real projection_tolerance(1e-10);
243 if (!u.is_zero() && u.norm().value() > 1.0)
244 projection_tolerance *= u.norm().value();
250 J << F(0).derivatives()[0], F(0).derivatives()[1], F(1).derivatives()[0],
251 F(1).derivatives()[1], F(2).derivatives()[0], F(2).derivatives()[1];
253 f << F(0).value(), F(1).value(), F(2).value();
258 }
while (++current_iterate < max_iterates);
260 if (current_iterate < max_iterates)
263 q_pts.push_back(transform_qp(xi(0).
value(), xi(1).
value()));
269 auto & qp_back = q_pts.back();
270 if (primal_elem->type() ==
TRI3 || primal_elem->type() ==
TRI6 || primal_elem->type() ==
TRI7)
272 if (qp_back(0) < -TOLERANCE || qp_back(1) < -TOLERANCE ||
273 qp_back(0) + qp_back(1) > (1 + TOLERANCE))
274 mooseException(
"Quadrature point: ", qp_back,
" out of bounds, truncating.");
276 else if (primal_elem->type() ==
QUAD4 || primal_elem->type() ==
QUAD8 ||
277 primal_elem->type() ==
QUAD9)
279 if (qp_back(0) < (-1 - TOLERANCE) || qp_back(0) > (1 + TOLERANCE) ||
280 qp_back(1) < (-1 - TOLERANCE) || qp_back(1) > (1 + TOLERANCE))
281 mooseException(
"Quadrature point: ", qp_back,
" out of bounds, truncating");
286 mooseError(
"Newton iteration for mortar quadrature mapping msm element: ",
290 " didn't converge. MSM element volume: ",
T fe_lagrange_2D_shape(const libMesh::ElemType type, const Order order, const unsigned int i, const VectorType< T > &p)
DualNumber< Real, NumberArray< 2, Real > > Dual2
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
Eigen::Matrix< Real, Eigen::Dynamic, Eigen::Dynamic > RealEigenMatrix
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
Eigen::Matrix< Real, Eigen::Dynamic, 1 > RealEigenVector
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
void projectQPoints3d(const Elem *msm_elem, const Elem *primal_elem, unsigned int sub_elem_index, const QBase &qrule_msm, std::vector< Point > &q_pts)
3D projection operator for mapping qpoints on mortar segments to secondary or primary elements ...