Line data Source code
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 821524 : 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 821524 : const auto msm_elem_order = msm_elem->default_order();
33 821524 : const auto msm_elem_type = msm_elem->type();
34 :
35 : // Get normal to linearized element, could store and query but computation is easy
36 821524 : const Point e1 = msm_elem->point(0) - msm_elem->point(1);
37 821524 : const Point e2 = msm_elem->point(2) - msm_elem->point(1);
38 821524 : const Point normal = e2.cross(e1).unit();
39 :
40 : // Get sub-elem (for second order meshes, otherwise trivial)
41 821524 : const auto sub_elem = msm_elem->get_extra_integer(sub_elem_index);
42 821524 : const ElemType primal_type = primal_elem->type();
43 :
44 821524 : auto get_sub_elem_node_indices = [primal_type, sub_elem]() -> std::vector<unsigned int>
45 : {
46 821524 : switch (primal_type)
47 : {
48 3584 : case TRI3:
49 10752 : return {{0, 1, 2}};
50 295536 : case QUAD4:
51 886608 : return {{0, 1, 2, 3}};
52 17152 : case TRI6:
53 : case TRI7:
54 17152 : switch (sub_elem)
55 : {
56 4288 : case 0:
57 12864 : return {{0, 3, 5}};
58 4288 : case 1:
59 12864 : return {{3, 4, 5}};
60 4288 : case 2:
61 12864 : return {{3, 1, 4}};
62 4288 : case 3:
63 12864 : return {{5, 4, 2}};
64 0 : default:
65 0 : mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
66 : }
67 430408 : case QUAD8:
68 430408 : switch (sub_elem)
69 : {
70 68308 : case 0:
71 204924 : return {{0, 4, 7}};
72 78428 : case 1:
73 235284 : return {{4, 1, 5}};
74 68483 : case 2:
75 205449 : return {{5, 2, 6}};
76 66724 : case 3:
77 200172 : return {{7, 6, 3}};
78 148465 : case 4:
79 445395 : return {{4, 5, 6, 7}};
80 0 : default:
81 0 : mooseError("get_sub_elem_nodes: Invalid sub_elem: ", sub_elem);
82 : }
83 74844 : case QUAD9:
84 74844 : switch (sub_elem)
85 : {
86 18063 : case 0:
87 54189 : return {{0, 4, 8, 7}};
88 19305 : case 1:
89 57915 : return {{4, 1, 5, 8}};
90 18063 : case 2:
91 54189 : return {{8, 5, 2, 6}};
92 19413 : case 3:
93 58239 : return {{7, 8, 6, 3}};
94 0 : default:
95 0 : mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
96 : }
97 0 : default:
98 0 : mooseError("get_sub_elem_indices: Face element type: ",
99 0 : libMesh::Utility::enum_to_string<ElemType>(primal_type),
100 : " invalid for 3D mortar");
101 : }
102 821524 : };
103 :
104 : // Transforms quadrature point from first order sub-elements (in case of second-order)
105 : // to primal element
106 4896188 : auto transform_qp = [primal_type, sub_elem](const Real nu, const Real xi)
107 : {
108 4896188 : switch (primal_type)
109 : {
110 14336 : case TRI3:
111 14336 : return Point(nu, xi, 0);
112 1182144 : case QUAD4:
113 1182144 : return Point(nu, xi, 0);
114 162944 : case TRI6:
115 : case TRI7:
116 162944 : switch (sub_elem)
117 : {
118 40736 : case 0:
119 40736 : return Point(0.5 * nu, 0.5 * xi, 0);
120 40736 : case 1:
121 40736 : return Point(0.5 * (1 - xi), 0.5 * (nu + xi), 0);
122 40736 : case 2:
123 40736 : return Point(0.5 * (1 + nu), 0.5 * xi, 0);
124 40736 : case 3:
125 40736 : return Point(0.5 * nu, 0.5 * (1 + xi), 0);
126 0 : default:
127 0 : mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
128 : }
129 3012856 : case QUAD8:
130 3012856 : switch (sub_elem)
131 : {
132 478156 : case 0:
133 478156 : return Point(nu - 1, xi - 1, 0);
134 548996 : case 1:
135 548996 : return Point(nu + xi, xi - 1, 0);
136 479381 : case 2:
137 479381 : return Point(1 - xi, nu + xi, 0);
138 467068 : case 3:
139 467068 : return Point(nu - 1, nu + xi, 0);
140 1039255 : case 4:
141 1039255 : return Point(0.5 * (nu - xi), 0.5 * (nu + xi), 0);
142 0 : default:
143 0 : mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
144 : }
145 523908 : case QUAD9:
146 523908 : switch (sub_elem)
147 : {
148 126441 : case 0:
149 126441 : return Point(0.5 * (nu - 1), 0.5 * (xi - 1), 0);
150 135135 : case 1:
151 135135 : return Point(0.5 * (nu + 1), 0.5 * (xi - 1), 0);
152 126441 : case 2:
153 126441 : return Point(0.5 * (nu + 1), 0.5 * (xi + 1), 0);
154 135891 : case 3:
155 135891 : return Point(0.5 * (nu - 1), 0.5 * (xi + 1), 0);
156 0 : default:
157 0 : mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
158 : }
159 0 : default:
160 0 : mooseError("transform_qp: Face element type: ",
161 0 : libMesh::Utility::enum_to_string<ElemType>(primal_type),
162 : " invalid for 3D mortar");
163 : }
164 821524 : };
165 :
166 51735877 : auto sub_element_type = [primal_type, sub_elem]()
167 : {
168 51735877 : switch (primal_type)
169 : {
170 1010784 : case TRI3:
171 : case TRI6:
172 : case TRI7:
173 1010784 : return TRI3;
174 26332296 : case QUAD4:
175 : case QUAD9:
176 26332296 : return QUAD4;
177 24392797 : case QUAD8:
178 24392797 : switch (sub_elem)
179 : {
180 11819925 : case 0:
181 : case 1:
182 : case 2:
183 : case 3:
184 11819925 : return TRI3;
185 12572872 : case 4:
186 12572872 : return QUAD4;
187 0 : default:
188 0 : mooseError("sub_element_type: Invalid sub_elem: ", sub_elem);
189 : }
190 0 : default:
191 0 : mooseError("sub_element_type: Face element type: ",
192 0 : libMesh::Utility::enum_to_string<ElemType>(primal_type),
193 : " invalid for 3D mortar");
194 : }
195 821524 : };
196 :
197 : // Get sub-elem node indices
198 821524 : auto sub_elem_node_indices = get_sub_elem_node_indices();
199 :
200 : // Loop through quadrature points on msm_elem
201 5717712 : for (auto qp : make_range(qrule_msm.n_points()))
202 : {
203 : // Get physical point on msm_elem to project
204 4896188 : Point x0;
205 19584752 : for (auto n : make_range(msm_elem->n_nodes()))
206 14688564 : x0 += Moose::fe_lagrange_2D_shape(msm_elem_type,
207 : msm_elem_order,
208 : n,
209 14688564 : static_cast<const TypeVector<Real> &>(qrule_msm.qp(qp))) *
210 29377128 : msm_elem->point(n);
211 :
212 : // Use msm_elem quadrature point as initial guess
213 : // (will be correct for aligned meshes)
214 4896188 : Dual2 xi1{};
215 4896188 : xi1.value() = qrule_msm.qp(qp)(0);
216 4896188 : xi1.derivatives()[0] = 1.0;
217 4896188 : Dual2 xi2{};
218 4896188 : xi2.value() = qrule_msm.qp(qp)(1);
219 4896188 : xi2.derivatives()[1] = 1.0;
220 4896188 : VectorValue<Dual2> xi(xi1, xi2, 0);
221 4896188 : 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 14003195 : VectorValue<Dual2> x1;
228 65739072 : for (auto n : make_range(sub_elem_node_indices.size()))
229 103471754 : x1 += Moose::fe_lagrange_2D_shape(sub_element_type(), FIRST, n, xi) *
230 103471754 : primal_elem->point(sub_elem_node_indices[n]);
231 14003195 : auto u = x1 - x0;
232 :
233 28006390 : VectorValue<Dual2> F(u(1) * normal(2) - u(2) * normal(1),
234 28006390 : u(2) * normal(0) - u(0) * normal(2),
235 42009585 : u(0) * normal(1) - u(1) * normal(0));
236 :
237 14003195 : 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 14003195 : if (!u.is_zero() && u.norm().value() > 1.0)
244 320594 : projection_tolerance *= u.norm().value();
245 :
246 14003195 : if (MetaPhysicL::raw_value(F).norm() < projection_tolerance)
247 4896188 : break;
248 :
249 9107007 : RealEigenMatrix J(3, 2);
250 18214014 : J << F(0).derivatives()[0], F(0).derivatives()[1], F(1).derivatives()[0],
251 9107007 : F(1).derivatives()[1], F(2).derivatives()[0], F(2).derivatives()[1];
252 9107007 : RealEigenVector f(3);
253 9107007 : f << F(0).value(), F(1).value(), F(2).value();
254 9107007 : const RealEigenVector dxi = -J.colPivHouseholderQr().solve(f);
255 :
256 9107007 : xi(0) += dxi(0);
257 9107007 : xi(1) += dxi(1);
258 32902578 : } while (++current_iterate < max_iterates);
259 :
260 4896188 : if (current_iterate < max_iterates)
261 : {
262 : // Transfer quadrature point from sub-element to element and store
263 4896188 : 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 4896188 : auto & qp_back = q_pts.back();
270 4896188 : if (primal_elem->type() == TRI3 || primal_elem->type() == TRI6 || primal_elem->type() == TRI7)
271 : {
272 354560 : if (qp_back(0) < -TOLERANCE || qp_back(1) < -TOLERANCE ||
273 177280 : qp_back(0) + qp_back(1) > (1 + TOLERANCE))
274 0 : mooseException("Quadrature point: ", qp_back, " out of bounds, truncating.");
275 : }
276 5242816 : else if (primal_elem->type() == QUAD4 || primal_elem->type() == QUAD8 ||
277 523908 : primal_elem->type() == QUAD9)
278 : {
279 9437816 : if (qp_back(0) < (-1 - TOLERANCE) || qp_back(0) > (1 + TOLERANCE) ||
280 9437816 : qp_back(1) < (-1 - TOLERANCE) || qp_back(1) > (1 + TOLERANCE))
281 0 : mooseException("Quadrature point: ", qp_back, " out of bounds, truncating");
282 : }
283 : }
284 : else
285 : {
286 0 : mooseError("Newton iteration for mortar quadrature mapping msm element: ",
287 0 : msm_elem->id(),
288 : " to elem: ",
289 0 : primal_elem->id(),
290 : " didn't converge. MSM element volume: ",
291 0 : msm_elem->volume());
292 : }
293 4896188 : }
294 821524 : }
295 : }
296 : }
|
