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 378080 : 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 378080 : const auto msm_elem_order = msm_elem->default_order();
33 378080 : const auto msm_elem_type = msm_elem->type();
34 :
35 : // Get normal to linearized element, could store and query but computation is easy
36 378080 : const Point e1 = msm_elem->point(0) - msm_elem->point(1);
37 378080 : const Point e2 = msm_elem->point(2) - msm_elem->point(1);
38 378080 : const Point normal = e2.cross(e1).unit();
39 :
40 : // Get sub-elem (for second order meshes, otherwise trivial)
41 378080 : const auto sub_elem = msm_elem->get_extra_integer(sub_elem_index);
42 378080 : const ElemType primal_type = primal_elem->type();
43 :
44 592560 : auto get_sub_elem_node_indices = [primal_type, sub_elem]() -> std::vector<unsigned int>
45 : {
46 378080 : switch (primal_type)
47 : {
48 3584 : case TRI3:
49 3584 : return {{0, 1, 2}};
50 160016 : case QUAD4:
51 160016 : return {{0, 1, 2, 3}};
52 11776 : case TRI6:
53 : case TRI7:
54 11776 : switch (sub_elem)
55 : {
56 2944 : case 0:
57 2944 : return {{0, 3, 5}};
58 2944 : case 1:
59 2944 : return {{3, 4, 5}};
60 2944 : case 2:
61 2944 : return {{3, 1, 4}};
62 2944 : case 3:
63 2944 : return {{5, 4, 2}};
64 0 : default:
65 0 : mooseError("get_sub_elem_indices: Invalid sub_elem: ", sub_elem);
66 : }
67 161520 : case QUAD8:
68 161520 : switch (sub_elem)
69 : {
70 25272 : case 0:
71 25272 : return {{0, 4, 7}};
72 28200 : case 1:
73 28200 : return {{4, 1, 5}};
74 25302 : case 2:
75 25302 : return {{5, 2, 6}};
76 26472 : case 3:
77 26472 : return {{7, 6, 3}};
78 56274 : case 4:
79 56274 : 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 378080 : };
103 :
104 : // Transforms quadrature point from first order sub-elements (in case of second-order)
105 : // to primal element
106 3716000 : auto transform_qp = [primal_type, sub_elem](const Real nu, const Real xi)
107 : {
108 2185200 : switch (primal_type)
109 : {
110 14336 : case TRI3:
111 14336 : return Point(nu, xi, 0);
112 640064 : case QUAD4:
113 640064 : return Point(nu, xi, 0);
114 111872 : case TRI6:
115 : case TRI7:
116 111872 : switch (sub_elem)
117 : {
118 27968 : case 0:
119 27968 : return Point(0.5 * nu, 0.5 * xi, 0);
120 27968 : case 1:
121 27968 : return Point(0.5 * (1 - xi), 0.5 * (nu + xi), 0);
122 27968 : case 2:
123 27968 : return Point(0.5 * (1 + nu), 0.5 * xi, 0);
124 27968 : case 3:
125 27968 : 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 1130640 : case QUAD8:
130 1130640 : switch (sub_elem)
131 : {
132 176904 : case 0:
133 176904 : return Point(nu - 1, xi - 1, 0);
134 197400 : case 1:
135 197400 : return Point(nu + xi, xi - 1, 0);
136 177114 : case 2:
137 177114 : return Point(1 - xi, nu + xi, 0);
138 185304 : case 3:
139 185304 : return Point(nu - 1, nu + xi, 0);
140 393918 : case 4:
141 393918 : 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 378080 : };
165 :
166 33164692 : auto sub_element_type = [primal_type, sub_elem]()
167 : {
168 23923786 : switch (primal_type)
169 : {
170 720480 : case TRI3:
171 : case TRI6:
172 : case TRI7:
173 720480 : return TRI3;
174 13962400 : case QUAD4:
175 : case QUAD9:
176 13962400 : return QUAD4;
177 9240906 : case QUAD8:
178 9240906 : switch (sub_elem)
179 : {
180 4411962 : case 0:
181 : case 1:
182 : case 2:
183 : case 3:
184 4411962 : return TRI3;
185 4828944 : case 4:
186 4828944 : 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 378080 : };
196 :
197 : // Get sub-elem node indices
198 378080 : auto sub_elem_node_indices = get_sub_elem_node_indices();
199 :
200 : // Loop through quadrature points on msm_elem
201 2563280 : for (auto qp : make_range(qrule_msm.n_points()))
202 : {
203 : // Get physical point on msm_elem to project
204 2185200 : Point x0;
205 8740800 : for (auto n : make_range(msm_elem->n_nodes()))
206 6555600 : x0 += Moose::fe_lagrange_2D_shape(msm_elem_type,
207 : msm_elem_order,
208 : n,
209 6555600 : static_cast<const TypeVector<Real> &>(qrule_msm.qp(qp))) *
210 13111200 : msm_elem->point(n);
211 :
212 : // Use msm_elem quadrature point as initial guess
213 : // (will be correct for aligned meshes)
214 2185200 : Dual2 xi1{};
215 2185200 : xi1.value() = qrule_msm.qp(qp)(0);
216 2185200 : xi1.derivatives()[0] = 1.0;
217 2185200 : Dual2 xi2{};
218 2185200 : xi2.value() = qrule_msm.qp(qp)(1);
219 2185200 : xi2.derivatives()[1] = 1.0;
220 2185200 : VectorValue<Dual2> xi(xi1, xi2, 0);
221 2185200 : 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 6408650 : VectorValue<Dual2> x1;
228 30332436 : for (auto n : make_range(sub_elem_node_indices.size()))
229 47847572 : x1 += Moose::fe_lagrange_2D_shape(sub_element_type(), FIRST, n, xi) *
230 47847572 : primal_elem->point(sub_elem_node_indices[n]);
231 6408650 : auto u = x1 - x0;
232 :
233 12817300 : VectorValue<Dual2> F(u(1) * normal(2) - u(2) * normal(1),
234 12817300 : u(2) * normal(0) - u(0) * normal(2),
235 19225950 : u(0) * normal(1) - u(1) * normal(0));
236 :
237 6408650 : 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 6408650 : if (!u.is_zero() && u.norm().value() > 1.0)
244 308498 : projection_tolerance *= u.norm().value();
245 :
246 6408650 : if (MetaPhysicL::raw_value(F).norm() < projection_tolerance)
247 2185200 : break;
248 :
249 4223450 : RealEigenMatrix J(3, 2);
250 8446900 : J << F(0).derivatives()[0], F(0).derivatives()[1], F(1).derivatives()[0],
251 4223450 : F(1).derivatives()[1], F(2).derivatives()[0], F(2).derivatives()[1];
252 4223450 : RealEigenVector f(3);
253 4223450 : f << F(0).value(), F(1).value(), F(2).value();
254 4223450 : const RealEigenVector dxi = -J.colPivHouseholderQr().solve(f);
255 :
256 4223450 : xi(0) += dxi(0);
257 4223450 : xi(1) += dxi(1);
258 15002500 : } while (++current_iterate < max_iterates);
259 :
260 2185200 : if (current_iterate < max_iterates)
261 : {
262 : // Transfer quadrature point from sub-element to element and store
263 2185200 : 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 2185200 : auto & qp_back = q_pts.back();
270 2185200 : if (primal_elem->type() == TRI3 || primal_elem->type() == TRI6 || primal_elem->type() == TRI7)
271 : {
272 252416 : if (qp_back(0) < -TOLERANCE || qp_back(1) < -TOLERANCE ||
273 126208 : qp_back(0) + qp_back(1) > (1 + TOLERANCE))
274 0 : mooseException("Quadrature point: ", qp_back, " out of bounds, truncating.");
275 : }
276 2347280 : else if (primal_elem->type() == QUAD4 || primal_elem->type() == QUAD8 ||
277 288288 : primal_elem->type() == QUAD9)
278 : {
279 4117984 : if (qp_back(0) < (-1 - TOLERANCE) || qp_back(0) > (1 + TOLERANCE) ||
280 4117984 : 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 2185200 : }
294 378080 : }
295 : }
296 : }
|
