https://mooseframework.inl.gov
MortarUtils.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "MortarUtils.h"
11 #include "MooseLagrangeHelpers.h"
12 
13 #include "libmesh/enum_to_string.h"
14 #include "metaphysicl/dualnumberarray.h"
15 #include "Eigen/Dense"
16 
17 using MetaPhysicL::NumberArray;
18 
19 typedef DualNumber<Real, NumberArray<2, Real>> Dual2;
20 
21 namespace Moose
22 {
23 namespace Mortar
24 {
25 void
26 projectQPoints3d(const Elem * const msm_elem,
27  const Elem * const primal_elem,
28  const unsigned int sub_elem_index,
29  const QBase & qrule_msm,
30  std::vector<Point> & q_pts)
31 {
32  const auto msm_elem_order = msm_elem->default_order();
33  const auto msm_elem_type = msm_elem->type();
34 
35  // Get normal to linearized element, could store and query but computation is easy
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();
39 
40  // Get sub-elem (for second order meshes, otherwise trivial)
41  const auto sub_elem = msm_elem->get_extra_integer(sub_elem_index);
42  const ElemType primal_type = primal_elem->type();
43 
44  auto get_sub_elem_node_indices = [primal_type, sub_elem]() -> std::vector<unsigned int>
45  {
46  switch (primal_type)
47  {
48  case TRI3:
49  return {{0, 1, 2}};
50  case QUAD4:
51  return {{0, 1, 2, 3}};
52  case TRI6:
53  case TRI7:
54  switch (sub_elem)
55  {
56  case 0:
57  return {{0, 3, 5}};
58  case 1:
59  return {{3, 4, 5}};
60  case 2:
61  return {{3, 1, 4}};
62  case 3:
63  return {{5, 4, 2}};
64  default:
65  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
66  }
67  case QUAD8:
68  switch (sub_elem)
69  {
70  case 0:
71  return {{0, 4, 7}};
72  case 1:
73  return {{4, 1, 5}};
74  case 2:
75  return {{5, 2, 6}};
76  case 3:
77  return {{7, 6, 3}};
78  case 4:
79  return {{4, 5, 6, 7}};
80  default:
81  mooseError("get_sub_elem_nodes: Invalid sub_elem: ", sub_elem);
82  }
83  case QUAD9:
84  switch (sub_elem)
85  {
86  case 0:
87  return {{0, 4, 8, 7}};
88  case 1:
89  return {{4, 1, 5, 8}};
90  case 2:
91  return {{8, 5, 2, 6}};
92  case 3:
93  return {{7, 8, 6, 3}};
94  default:
95  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
96  }
97  default:
98  mooseError("get_sub_elem_indices: Face element type: ",
99  libMesh::Utility::enum_to_string<ElemType>(primal_type),
100  " invalid for 3D mortar");
101  }
102  };
103 
104  // Transforms quadrature point from first order sub-elements (in case of second-order)
105  // to primal element
106  auto transform_qp = [primal_type, sub_elem](const Real nu, const Real xi)
107  {
108  switch (primal_type)
109  {
110  case TRI3:
111  return Point(nu, xi, 0);
112  case QUAD4:
113  return Point(nu, xi, 0);
114  case TRI6:
115  case TRI7:
116  switch (sub_elem)
117  {
118  case 0:
119  return Point(0.5 * nu, 0.5 * xi, 0);
120  case 1:
121  return Point(0.5 * (1 - xi), 0.5 * (nu + xi), 0);
122  case 2:
123  return Point(0.5 * (1 + nu), 0.5 * xi, 0);
124  case 3:
125  return Point(0.5 * nu, 0.5 * (1 + xi), 0);
126  default:
127  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
128  }
129  case QUAD8:
130  switch (sub_elem)
131  {
132  case 0:
133  return Point(nu - 1, xi - 1, 0);
134  case 1:
135  return Point(nu + xi, xi - 1, 0);
136  case 2:
137  return Point(1 - xi, nu + xi, 0);
138  case 3:
139  return Point(nu - 1, nu + xi, 0);
140  case 4:
141  return Point(0.5 * (nu - xi), 0.5 * (nu + xi), 0);
142  default:
143  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
144  }
145  case QUAD9:
146  switch (sub_elem)
147  {
148  case 0:
149  return Point(0.5 * (nu - 1), 0.5 * (xi - 1), 0);
150  case 1:
151  return Point(0.5 * (nu + 1), 0.5 * (xi - 1), 0);
152  case 2:
153  return Point(0.5 * (nu + 1), 0.5 * (xi + 1), 0);
154  case 3:
155  return Point(0.5 * (nu - 1), 0.5 * (xi + 1), 0);
156  default:
157  mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
158  }
159  default:
160  mooseError("transform_qp: Face element type: ",
161  libMesh::Utility::enum_to_string<ElemType>(primal_type),
162  " invalid for 3D mortar");
163  }
164  };
165 
166  auto sub_element_type = [primal_type, sub_elem]()
167  {
168  switch (primal_type)
169  {
170  case TRI3:
171  case TRI6:
172  case TRI7:
173  return TRI3;
174  case QUAD4:
175  case QUAD9:
176  return QUAD4;
177  case QUAD8:
178  switch (sub_elem)
179  {
180  case 0:
181  case 1:
182  case 2:
183  case 3:
184  return TRI3;
185  case 4:
186  return QUAD4;
187  default:
188  mooseError("sub_element_type: Invalid sub_elem: ", sub_elem);
189  }
190  default:
191  mooseError("sub_element_type: Face element type: ",
192  libMesh::Utility::enum_to_string<ElemType>(primal_type),
193  " invalid for 3D mortar");
194  }
195  };
196 
197  // Get sub-elem node indices
198  auto sub_elem_node_indices = get_sub_elem_node_indices();
199 
200  // Loop through quadrature points on msm_elem
201  for (auto qp : make_range(qrule_msm.n_points()))
202  {
203  // Get physical point on msm_elem to project
204  Point x0;
205  for (auto n : make_range(msm_elem->n_nodes()))
206  x0 += Moose::fe_lagrange_2D_shape(msm_elem_type,
207  msm_elem_order,
208  n,
209  static_cast<const TypeVector<Real> &>(qrule_msm.qp(qp))) *
210  msm_elem->point(n);
211 
212  // Use msm_elem quadrature point as initial guess
213  // (will be correct for aligned meshes)
214  Dual2 xi1{};
215  xi1.value() = qrule_msm.qp(qp)(0);
216  xi1.derivatives()[0] = 1.0;
217  Dual2 xi2{};
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;
222 
223  // Project qp from mortar segments to first order sub-elements (elements in case of first order
224  // geometry)
225  do
226  {
227  VectorValue<Dual2> x1;
228  for (auto n : make_range(sub_elem_node_indices.size()))
229  x1 += Moose::fe_lagrange_2D_shape(sub_element_type(), FIRST, n, xi) *
230  primal_elem->point(sub_elem_node_indices[n]);
231  auto u = x1 - x0;
232 
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));
236 
237  Real projection_tolerance(1e-10);
238 
239  // Normalize tolerance with quantities involved in the projection.
240  // Absolute projection tolerance is loosened for displacements larger than those on the order
241  // of one. Tightening the tolerance for displacements of smaller orders causes this tolerance
242  // to not be reached in a number of tests.
243  if (!u.is_zero() && u.norm().value() > 1.0)
244  projection_tolerance *= u.norm().value();
245 
246  if (MetaPhysicL::raw_value(F).norm() < projection_tolerance)
247  break;
248 
249  RealEigenMatrix J(3, 2);
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];
252  RealEigenVector f(3);
253  f << F(0).value(), F(1).value(), F(2).value();
254  const RealEigenVector dxi = -J.colPivHouseholderQr().solve(f);
255 
256  xi(0) += dxi(0);
257  xi(1) += dxi(1);
258  } while (++current_iterate < max_iterates);
259 
260  if (current_iterate < max_iterates)
261  {
262  // Transfer quadrature point from sub-element to element and store
263  q_pts.push_back(transform_qp(xi(0).value(), xi(1).value()));
264 
265  // The following checks if quadrature point falls in correct domain.
266  // On small mortar segment elements with very distorted elements this can fail, instead of
267  // erroring simply truncate quadrature point, these points typically have very small
268  // contributions to integrals
269  auto & qp_back = q_pts.back();
270  if (primal_elem->type() == TRI3 || primal_elem->type() == TRI6 || primal_elem->type() == TRI7)
271  {
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.");
275  }
276  else if (primal_elem->type() == QUAD4 || primal_elem->type() == QUAD8 ||
277  primal_elem->type() == QUAD9)
278  {
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");
282  }
283  }
284  else
285  {
286  mooseError("Newton iteration for mortar quadrature mapping msm element: ",
287  msm_elem->id(),
288  " to elem: ",
289  primal_elem->id(),
290  " didn't converge. MSM element volume: ",
291  msm_elem->volume());
292  }
293  }
294 }
295 }
296 }
T fe_lagrange_2D_shape(const libMesh::ElemType type, const Order order, const unsigned int i, const VectorType< T > &p)
ElemType
QUAD8
DualNumber< Real, NumberArray< 2, Real > > Dual2
Definition: MortarUtils.C:19
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
FIRST
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
TRI3
QUAD4
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
TRI6
Eigen::Matrix< Real, Eigen::Dynamic, Eigen::Dynamic > RealEigenMatrix
Definition: MooseTypes.h:149
auto norm(const T &a) -> decltype(std::abs(a))
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
TRI7
QUAD9
Eigen::Matrix< Real, Eigen::Dynamic, 1 > RealEigenVector
Definition: MooseTypes.h:146
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 ...
Definition: MortarUtils.C:26