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 "libmesh/fe_interface.h"
15 : #include "metaphysicl/dualnumberarray.h"
16 : #include "Eigen/Dense"
17 :
18 : using MetaPhysicL::NumberArray;
19 :
20 : typedef DualNumber<Real, NumberArray<2, Real>> Dual2;
21 :
22 : namespace Moose
23 : {
24 : namespace Mortar
25 : {
26 : std::vector<unsigned int>
27 1049890 : getMortarSubElementNodeIndices(const Elem & parent_elem, const unsigned int sub_elem)
28 : {
29 1049890 : if (sub_elem >= parent_elem.n_sub_elem())
30 0 : mooseError("Invalid 3D mortar sub-element index ",
31 : sub_elem,
32 : " for parent element ",
33 0 : parent_elem.id(),
34 : " of type ",
35 0 : libMesh::Utility::enum_to_string<ElemType>(parent_elem.type()),
36 : ", which has ",
37 0 : parent_elem.n_sub_elem(),
38 : " sub-elements.");
39 :
40 1049890 : switch (parent_elem.type())
41 : {
42 9131 : case TRI3:
43 27393 : return {0, 1, 2};
44 437731 : case QUAD4:
45 1313193 : return {0, 1, 2, 3};
46 34192 : case TRI6:
47 : case TRI7:
48 34192 : switch (sub_elem)
49 : {
50 8548 : case 0:
51 25644 : return {0, 3, 5};
52 8548 : case 1:
53 25644 : return {3, 4, 5};
54 8548 : case 2:
55 25644 : return {3, 1, 4};
56 8548 : case 3:
57 25644 : return {5, 4, 2};
58 0 : default:
59 0 : mooseError("Invalid 3D mortar triangular sub-element index ", sub_elem, ".");
60 : }
61 485104 : case QUAD8:
62 485104 : switch (sub_elem)
63 : {
64 79559 : case 0:
65 238677 : return {0, 4, 7};
66 89031 : case 1:
67 267093 : return {4, 1, 5};
68 79559 : case 2:
69 238677 : return {5, 2, 6};
70 77921 : case 3:
71 233763 : return {7, 6, 3};
72 159034 : case 4:
73 477102 : return {4, 5, 6, 7};
74 0 : default:
75 0 : mooseError("Invalid 3D mortar QUAD8 sub-element index ", sub_elem, ".");
76 : }
77 83732 : case QUAD9:
78 83732 : switch (sub_elem)
79 : {
80 20285 : case 0:
81 60855 : return {0, 4, 8, 7};
82 21527 : case 1:
83 64581 : return {4, 1, 5, 8};
84 20285 : case 2:
85 60855 : return {8, 5, 2, 6};
86 21635 : case 3:
87 64905 : return {7, 8, 6, 3};
88 0 : default:
89 0 : mooseError("Invalid 3D mortar QUAD9 sub-element index ", sub_elem, ".");
90 : }
91 0 : default:
92 0 : mooseError("Parent face element ",
93 0 : parent_elem.id(),
94 : " has unsupported type ",
95 0 : libMesh::Utility::enum_to_string<ElemType>(parent_elem.type()),
96 : " for 3D mortar sub-element topology.");
97 : }
98 : }
99 :
100 : void
101 44034 : mapQPoints3dFromReference(const Elem & mortar_segment_elem,
102 : const MortarSegmentReferencePoints & reference_points,
103 : const QBase & qrule_msm,
104 : std::vector<Point> & secondary_q_pts,
105 : std::vector<Point> & primary_q_pts)
106 : {
107 : mooseAssert(mortar_segment_elem.type() == TRI3,
108 : "Reference interpolation expects triangular mortar segments.");
109 44034 : const FEType fe_type(FIRST, LAGRANGE);
110 :
111 352272 : for (const auto qp : make_range(qrule_msm.n_points()))
112 : {
113 308238 : Point secondary_qp;
114 308238 : Point primary_qp;
115 :
116 1232952 : for (const auto n : index_range(reference_points.secondary_reference_points))
117 : {
118 : const auto phi =
119 924714 : FEInterface::shape(fe_type, &mortar_segment_elem, n, qrule_msm.qp(qp), false);
120 924714 : secondary_qp += phi * reference_points.secondary_reference_points[n];
121 924714 : primary_qp += phi * reference_points.primary_reference_points[n];
122 : }
123 :
124 308238 : secondary_q_pts.push_back(secondary_qp);
125 308238 : primary_q_pts.push_back(primary_qp);
126 : }
127 44034 : }
128 :
129 : void
130 911564 : projectQPoints3d(const Elem * const msm_elem,
131 : const Elem * const primal_elem,
132 : const unsigned int sub_elem_index,
133 : const QBase & qrule_msm,
134 : std::vector<Point> & q_pts)
135 : {
136 911564 : const auto msm_elem_order = msm_elem->default_order();
137 911564 : const auto msm_elem_type = msm_elem->type();
138 :
139 : // Get normal to linearized element, could store and query but computation is easy
140 911564 : const Point e1 = msm_elem->point(0) - msm_elem->point(1);
141 911564 : const Point e2 = msm_elem->point(2) - msm_elem->point(1);
142 911564 : const Point normal = e2.cross(e1).unit();
143 :
144 : // Get sub-elem (for second order meshes, otherwise trivial)
145 911564 : const auto sub_elem = msm_elem->get_extra_integer(sub_elem_index);
146 911564 : const ElemType primal_type = primal_elem->type();
147 :
148 : // Transforms quadrature point from first order sub-elements (in case of second-order)
149 : // to primal element
150 5259276 : auto transform_qp = [primal_type, sub_elem](const Real nu, const Real xi)
151 : {
152 5259276 : switch (primal_type)
153 : {
154 14336 : case TRI3:
155 14336 : return Point(nu, xi, 0);
156 1538400 : case QUAD4:
157 1538400 : return Point(nu, xi, 0);
158 162944 : case TRI6:
159 : case TRI7:
160 162944 : switch (sub_elem)
161 : {
162 40736 : case 0:
163 40736 : return Point(0.5 * nu, 0.5 * xi, 0);
164 40736 : case 1:
165 40736 : return Point(0.5 * (1 - xi), 0.5 * (nu + xi), 0);
166 40736 : case 2:
167 40736 : return Point(0.5 * (1 + nu), 0.5 * xi, 0);
168 40736 : case 3:
169 40736 : return Point(0.5 * nu, 0.5 * (1 + xi), 0);
170 0 : default:
171 0 : mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
172 : }
173 3019688 : case QUAD8:
174 3019688 : switch (sub_elem)
175 : {
176 481705 : case 0:
177 481705 : return Point(nu - 1, xi - 1, 0);
178 548009 : case 1:
179 548009 : return Point(nu + xi, xi - 1, 0);
180 481705 : case 2:
181 481705 : return Point(1 - xi, nu + xi, 0);
182 470239 : case 3:
183 470239 : return Point(nu - 1, nu + xi, 0);
184 1038030 : case 4:
185 1038030 : return Point(0.5 * (nu - xi), 0.5 * (nu + xi), 0);
186 0 : default:
187 0 : mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
188 : }
189 523908 : case QUAD9:
190 523908 : switch (sub_elem)
191 : {
192 126441 : case 0:
193 126441 : return Point(0.5 * (nu - 1), 0.5 * (xi - 1), 0);
194 135135 : case 1:
195 135135 : return Point(0.5 * (nu + 1), 0.5 * (xi - 1), 0);
196 126441 : case 2:
197 126441 : return Point(0.5 * (nu + 1), 0.5 * (xi + 1), 0);
198 135891 : case 3:
199 135891 : return Point(0.5 * (nu - 1), 0.5 * (xi + 1), 0);
200 0 : default:
201 0 : mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
202 : }
203 0 : default:
204 0 : mooseError("transform_qp: Face element type: ",
205 0 : libMesh::Utility::enum_to_string<ElemType>(primal_type),
206 : " invalid for 3D mortar");
207 : }
208 911564 : };
209 :
210 57476744 : auto sub_element_type = [primal_type, sub_elem]()
211 : {
212 57476744 : switch (primal_type)
213 : {
214 1010784 : case TRI3:
215 : case TRI6:
216 : case TRI7:
217 1010784 : return TRI3;
218 32040844 : case QUAD4:
219 : case QUAD9:
220 32040844 : return QUAD4;
221 24425116 : case QUAD8:
222 24425116 : switch (sub_elem)
223 : {
224 11867520 : case 0:
225 : case 1:
226 : case 2:
227 : case 3:
228 11867520 : return TRI3;
229 12557596 : case 4:
230 12557596 : return QUAD4;
231 0 : default:
232 0 : mooseError("sub_element_type: Invalid sub_elem: ", sub_elem);
233 : }
234 0 : default:
235 0 : mooseError("sub_element_type: Face element type: ",
236 0 : libMesh::Utility::enum_to_string<ElemType>(primal_type),
237 : " invalid for 3D mortar");
238 : }
239 911564 : };
240 :
241 : // Get sub-elem node indices
242 911564 : const auto sub_elem_node_indices = getMortarSubElementNodeIndices(*primal_elem, sub_elem);
243 :
244 : // Loop through quadrature points on msm_elem
245 6170840 : for (auto qp : make_range(qrule_msm.n_points()))
246 : {
247 : // Get physical point on msm_elem to project
248 5259276 : Point x0;
249 21037104 : for (auto n : make_range(msm_elem->n_nodes()))
250 15777828 : x0 += Moose::fe_lagrange_2D_shape(msm_elem_type,
251 : msm_elem_order,
252 : n,
253 15777828 : static_cast<const TypeVector<Real> &>(qrule_msm.qp(qp))) *
254 31555656 : msm_elem->point(n);
255 :
256 : // Use msm_elem quadrature point as initial guess
257 : // (will be correct for aligned meshes)
258 5259276 : Dual2 xi1{};
259 5259276 : xi1.value() = qrule_msm.qp(qp)(0);
260 5259276 : xi1.derivatives()[0] = 1.0;
261 5259276 : Dual2 xi2{};
262 5259276 : xi2.value() = qrule_msm.qp(qp)(1);
263 5259276 : xi2.derivatives()[1] = 1.0;
264 5259276 : VectorValue<Dual2> xi(xi1, xi2, 0);
265 5259276 : unsigned int current_iterate = 0, max_iterates = 10;
266 :
267 : // Project qp from mortar segments to first order sub-elements (elements in case of first order
268 : // geometry)
269 : do
270 : {
271 15442378 : VectorValue<Dual2> x1;
272 72919122 : for (auto n : make_range(sub_elem_node_indices.size()))
273 114953488 : x1 += Moose::fe_lagrange_2D_shape(sub_element_type(), FIRST, n, xi) *
274 114953488 : primal_elem->point(sub_elem_node_indices[n]);
275 15442378 : auto u = x1 - x0;
276 :
277 30884756 : VectorValue<Dual2> F(u(1) * normal(2) - u(2) * normal(1),
278 30884756 : u(2) * normal(0) - u(0) * normal(2),
279 46327134 : u(0) * normal(1) - u(1) * normal(0));
280 :
281 15442378 : Real projection_tolerance(1e-10);
282 :
283 : // Normalize tolerance with quantities involved in the projection.
284 : // Absolute projection tolerance is loosened for displacements larger than those on the order
285 : // of one. Tightening the tolerance for displacements of smaller orders causes this tolerance
286 : // to not be reached in a number of tests.
287 15442378 : if (!u.is_zero() && u.norm().value() > 1.0)
288 320594 : projection_tolerance *= u.norm().value();
289 :
290 15442378 : if (MetaPhysicL::raw_value(F).norm() < projection_tolerance)
291 5259276 : break;
292 :
293 10183102 : RealEigenMatrix J(3, 2);
294 20366204 : J << F(0).derivatives()[0], F(0).derivatives()[1], F(1).derivatives()[0],
295 10183102 : F(1).derivatives()[1], F(2).derivatives()[0], F(2).derivatives()[1];
296 10183102 : RealEigenVector f(3);
297 10183102 : f << F(0).value(), F(1).value(), F(2).value();
298 10183102 : const RealEigenVector dxi = -J.colPivHouseholderQr().solve(f);
299 :
300 10183102 : xi(0) += dxi(0);
301 10183102 : xi(1) += dxi(1);
302 36144032 : } while (++current_iterate < max_iterates);
303 :
304 5259276 : if (current_iterate < max_iterates)
305 : {
306 : // Transfer quadrature point from sub-element to element and store
307 5259276 : q_pts.push_back(transform_qp(xi(0).value(), xi(1).value()));
308 :
309 : // The following checks if quadrature point falls in correct domain.
310 : // On small mortar segment elements with very distorted elements this can fail, instead of
311 : // erroring simply truncate quadrature point, these points typically have very small
312 : // contributions to integrals
313 5259276 : auto & qp_back = q_pts.back();
314 5259276 : if (primal_elem->type() == TRI3 || primal_elem->type() == TRI6 || primal_elem->type() == TRI7)
315 : {
316 354560 : if (qp_back(0) < -TOLERANCE || qp_back(1) < -TOLERANCE ||
317 177280 : qp_back(0) + qp_back(1) > (1 + TOLERANCE))
318 0 : mooseException("Quadrature point: ", qp_back, " out of bounds, truncating.");
319 : }
320 5605904 : else if (primal_elem->type() == QUAD4 || primal_elem->type() == QUAD8 ||
321 523908 : primal_elem->type() == QUAD9)
322 : {
323 10163992 : if (qp_back(0) < (-1 - TOLERANCE) || qp_back(0) > (1 + TOLERANCE) ||
324 10163992 : qp_back(1) < (-1 - TOLERANCE) || qp_back(1) > (1 + TOLERANCE))
325 0 : mooseException("Quadrature point: ", qp_back, " out of bounds, truncating");
326 : }
327 : }
328 : else
329 : {
330 0 : mooseError("Newton iteration for mortar quadrature mapping msm element: ",
331 0 : msm_elem->id(),
332 : " to elem: ",
333 0 : primal_elem->id(),
334 : " didn't converge. MSM element volume: ",
335 0 : msm_elem->volume());
336 : }
337 5259276 : }
338 911564 : }
339 : }
340 : }
|
