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 412452 : 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 412452 : const auto msm_elem_order = msm_elem->default_order();
33 412452 : const auto msm_elem_type = msm_elem->type();
34 :
35 : // Get normal to linearized element, could store and query but computation is easy
36 412452 : const Point e1 = msm_elem->point(0) - msm_elem->point(1);
37 412452 : const Point e2 = msm_elem->point(2) - msm_elem->point(1);
38 412452 : const Point normal = e2.cross(e1).unit();
39 :
40 : // Get sub-elem (for second order meshes, otherwise trivial)
41 412452 : const auto sub_elem = msm_elem->get_extra_integer(sub_elem_index);
42 412452 : const ElemType primal_type = primal_elem->type();
43 :
44 641448 : auto get_sub_elem_node_indices = [primal_type, sub_elem]() -> std::vector<unsigned int>
45 : {
46 412452 : switch (primal_type)
47 : {
48 4032 : case TRI3:
49 4032 : return {{0, 1, 2}};
50 179424 : case QUAD4:
51 179424 : return {{0, 1, 2, 3}};
52 13056 : case TRI6:
53 : case TRI7:
54 13056 : switch (sub_elem)
55 : {
56 3264 : case 0:
57 3264 : return {{0, 3, 5}};
58 3264 : case 1:
59 3264 : return {{3, 4, 5}};
60 3264 : case 2:
61 3264 : return {{3, 1, 4}};
62 3264 : case 3:
63 3264 : return {{5, 4, 2}};
64 0 : default:
65 0 : mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
66 : }
67 174756 : case QUAD8:
68 174756 : switch (sub_elem)
69 : {
70 27318 : case 0:
71 27318 : return {{0, 4, 7}};
72 30426 : case 1:
73 30426 : return {{4, 1, 5}};
74 27348 : case 2:
75 27348 : return {{5, 2, 6}};
76 28740 : case 3:
77 28740 : return {{7, 6, 3}};
78 60924 : case 4:
79 60924 : return {{4, 5, 6, 7}};
80 0 : default:
81 0 : mooseError("get_sub_elem_nodes: Invalid sub_elem: ", sub_elem);
82 : }
83 41184 : case QUAD9:
84 41184 : switch (sub_elem)
85 : {
86 10032 : case 0:
87 10032 : return {{0, 4, 8, 7}};
88 10596 : case 1:
89 10596 : return {{4, 1, 5, 8}};
90 9960 : case 2:
91 9960 : return {{8, 5, 2, 6}};
92 10596 : case 3:
93 10596 : 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 412452 : };
103 :
104 : // Transforms quadrature point from first order sub-elements (in case of second-order)
105 : // to primal element
106 4005048 : auto transform_qp = [primal_type, sub_elem](const Real nu, const Real xi)
107 : {
108 2369436 : switch (primal_type)
109 : {
110 16128 : case TRI3:
111 16128 : return Point(nu, xi, 0);
112 717696 : case QUAD4:
113 717696 : return Point(nu, xi, 0);
114 124032 : case TRI6:
115 : case TRI7:
116 124032 : switch (sub_elem)
117 : {
118 31008 : case 0:
119 31008 : return Point(0.5 * nu, 0.5 * xi, 0);
120 31008 : case 1:
121 31008 : return Point(0.5 * (1 - xi), 0.5 * (nu + xi), 0);
122 31008 : case 2:
123 31008 : return Point(0.5 * (1 + nu), 0.5 * xi, 0);
124 31008 : case 3:
125 31008 : 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 1223292 : case QUAD8:
130 1223292 : switch (sub_elem)
131 : {
132 191226 : case 0:
133 191226 : return Point(nu - 1, xi - 1, 0);
134 212982 : case 1:
135 212982 : return Point(nu + xi, xi - 1, 0);
136 191436 : case 2:
137 191436 : return Point(1 - xi, nu + xi, 0);
138 201180 : case 3:
139 201180 : return Point(nu - 1, nu + xi, 0);
140 426468 : case 4:
141 426468 : 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 288288 : case QUAD9:
146 288288 : switch (sub_elem)
147 : {
148 70224 : case 0:
149 70224 : return Point(0.5 * (nu - 1), 0.5 * (xi - 1), 0);
150 74172 : case 1:
151 74172 : return Point(0.5 * (nu + 1), 0.5 * (xi - 1), 0);
152 69720 : case 2:
153 69720 : return Point(0.5 * (nu + 1), 0.5 * (xi + 1), 0);
154 74172 : case 3:
155 74172 : 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 412452 : };
165 :
166 36029232 : auto sub_element_type = [primal_type, sub_elem]()
167 : {
168 26025072 : switch (primal_type)
169 : {
170 800172 : case TRI3:
171 : case TRI6:
172 : case TRI7:
173 800172 : return TRI3;
174 15220740 : case QUAD4:
175 : case QUAD9:
176 15220740 : return QUAD4;
177 10004160 : case QUAD8:
178 10004160 : switch (sub_elem)
179 : {
180 4771872 : case 0:
181 : case 1:
182 : case 2:
183 : case 3:
184 4771872 : return TRI3;
185 5232288 : case 4:
186 5232288 : 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 412452 : };
196 :
197 : // Get sub-elem node indices
198 412452 : auto sub_elem_node_indices = get_sub_elem_node_indices();
199 :
200 : // Loop through quadrature points on msm_elem
201 2781888 : for (auto qp : make_range(qrule_msm.n_points()))
202 : {
203 : // Get physical point on msm_elem to project
204 2369436 : Point x0;
205 9477744 : for (auto n : make_range(msm_elem->n_nodes()))
206 7108308 : x0 += Moose::fe_lagrange_2D_shape(msm_elem_type,
207 : msm_elem_order,
208 : n,
209 7108308 : static_cast<const TypeVector<Real> &>(qrule_msm.qp(qp))) *
210 14216616 : msm_elem->point(n);
211 :
212 : // Use msm_elem quadrature point as initial guess
213 : // (will be correct for aligned meshes)
214 2369436 : Dual2 xi1{};
215 2369436 : xi1.value() = qrule_msm.qp(qp)(0);
216 2369436 : xi1.derivatives()[0] = 1.0;
217 2369436 : Dual2 xi2{};
218 2369436 : xi2.value() = qrule_msm.qp(qp)(1);
219 2369436 : xi2.derivatives()[1] = 1.0;
220 2369436 : VectorValue<Dual2> xi(xi1, xi2, 0);
221 2369436 : 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 6970605 : VectorValue<Dual2> x1;
228 32995677 : for (auto n : make_range(sub_elem_node_indices.size()))
229 52050144 : x1 += Moose::fe_lagrange_2D_shape(sub_element_type(), FIRST, n, xi) *
230 52050144 : primal_elem->point(sub_elem_node_indices[n]);
231 6970605 : auto u = x1 - x0;
232 :
233 13941210 : VectorValue<Dual2> F(u(1) * normal(2) - u(2) * normal(1),
234 13941210 : u(2) * normal(0) - u(0) * normal(2),
235 20911815 : u(0) * normal(1) - u(1) * normal(0));
236 :
237 6970605 : 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 6970605 : if (!u.is_zero() && u.norm().value() > 1.0)
244 330408 : projection_tolerance *= u.norm().value();
245 :
246 6970605 : if (MetaPhysicL::raw_value(F).norm() < projection_tolerance)
247 2369436 : break;
248 :
249 4601169 : RealEigenMatrix J(3, 2);
250 9202338 : J << F(0).derivatives()[0], F(0).derivatives()[1], F(1).derivatives()[0],
251 4601169 : F(1).derivatives()[1], F(2).derivatives()[0], F(2).derivatives()[1];
252 4601169 : RealEigenVector f(3);
253 4601169 : f << F(0).value(), F(1).value(), F(2).value();
254 4601169 : const RealEigenVector dxi = -J.colPivHouseholderQr().solve(f);
255 :
256 4601169 : xi(0) += dxi(0);
257 4601169 : xi(1) += dxi(1);
258 16310646 : } while (++current_iterate < max_iterates);
259 :
260 2369436 : if (current_iterate < max_iterates)
261 : {
262 : // Transfer quadrature point from sub-element to element and store
263 2369436 : 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 2369436 : auto & qp_back = q_pts.back();
270 2369436 : if (primal_elem->type() == TRI3 || primal_elem->type() == TRI6 || primal_elem->type() == TRI7)
271 : {
272 280320 : if (qp_back(0) < -TOLERANCE || qp_back(1) < -TOLERANCE ||
273 140160 : qp_back(0) + qp_back(1) > (1 + TOLERANCE))
274 0 : mooseException("Quadrature point: ", qp_back, " out of bounds, truncating.");
275 : }
276 2517564 : else if (primal_elem->type() == QUAD4 || primal_elem->type() == QUAD8 ||
277 288288 : primal_elem->type() == QUAD9)
278 : {
279 4458552 : if (qp_back(0) < (-1 - TOLERANCE) || qp_back(0) > (1 + TOLERANCE) ||
280 4458552 : 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 2369436 : }
294 412452 : }
295 : }
296 : }
|
