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 "XFEMFuncs.h"
11 :
12 : #include "MooseError.h"
13 : #include "Conversion.h"
14 :
15 : using namespace libMesh;
16 :
17 : namespace Xfem
18 : {
19 :
20 : void
21 0 : dunavant_rule2(const Real * wts,
22 : const Real * a,
23 : const Real * b,
24 : const unsigned int * permutation_ids,
25 : unsigned int n_wts,
26 : std::vector<Point> & points,
27 : std::vector<Real> & weights)
28 : {
29 : // see libmesh/src/quadrature/quadrature_gauss.C
30 : // Figure out how many total points by summing up the entries
31 : // in the permutation_ids array, and resize the _points and _weights
32 : // vectors appropriately.
33 : unsigned int total_pts = 0;
34 0 : for (unsigned int p = 0; p < n_wts; ++p)
35 0 : total_pts += permutation_ids[p];
36 :
37 : // Resize point and weight vectors appropriately.
38 0 : points.resize(total_pts);
39 0 : weights.resize(total_pts);
40 :
41 : // Always insert into the points & weights vector relative to the offset
42 : unsigned int offset = 0;
43 0 : for (unsigned int p = 0; p < n_wts; ++p)
44 : {
45 0 : switch (permutation_ids[p])
46 : {
47 : case 1:
48 : {
49 : // The point has only a single permutation (the centroid!)
50 : // So we don't even need to look in the a or b arrays.
51 0 : points[offset + 0] = Point(1.0L / 3.0L, 1.0L / 3.0L);
52 0 : weights[offset + 0] = wts[p];
53 :
54 0 : offset += 1;
55 0 : break;
56 : }
57 :
58 0 : case 3:
59 : {
60 : // For this type of rule, don't need to look in the b array.
61 0 : points[offset + 0] = Point(a[p], a[p]); // (a,a)
62 0 : points[offset + 1] = Point(a[p], 1.L - 2.L * a[p]); // (a,1-2a)
63 0 : points[offset + 2] = Point(1.L - 2.L * a[p], a[p]); // (1-2a,a)
64 :
65 0 : for (unsigned int j = 0; j < 3; ++j)
66 0 : weights[offset + j] = wts[p];
67 :
68 0 : offset += 3;
69 0 : break;
70 : }
71 :
72 0 : case 6:
73 : {
74 : // This type of point uses all 3 arrays...
75 0 : points[offset + 0] = Point(a[p], b[p]);
76 0 : points[offset + 1] = Point(b[p], a[p]);
77 0 : points[offset + 2] = Point(a[p], 1.L - a[p] - b[p]);
78 0 : points[offset + 3] = Point(1.L - a[p] - b[p], a[p]);
79 0 : points[offset + 4] = Point(b[p], 1.L - a[p] - b[p]);
80 0 : points[offset + 5] = Point(1.L - a[p] - b[p], b[p]);
81 :
82 0 : for (unsigned int j = 0; j < 6; ++j)
83 0 : weights[offset + j] = wts[p];
84 :
85 0 : offset += 6;
86 0 : break;
87 : }
88 :
89 0 : default:
90 0 : mooseError("Unknown permutation id: ", permutation_ids[p], "!");
91 : }
92 : }
93 0 : }
94 :
95 : void
96 585056 : stdQuadr2D(unsigned int nen, unsigned int iord, std::vector<std::vector<Real>> & sg2)
97 : {
98 : // Purpose: get Guass integration points for 2D quad and tri elems
99 : // N.B. only works for n_qp <= 6
100 :
101 585056 : Real lr4[4] = {-1.0, 1.0, -1.0, 1.0}; // libmesh order
102 585056 : Real lz4[4] = {-1.0, -1.0, 1.0, 1.0};
103 585056 : Real lr9[9] = {-1.0, 0.0, 1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0}; // libmesh order
104 585056 : Real lz9[9] = {-1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0};
105 585056 : Real lw9[9] = {25.0, 40.0, 25.0, 40.0, 64.0, 40.0, 25.0, 40.0, 25.0};
106 :
107 585056 : if (nen == 4) // 2d quad element
108 : {
109 204 : if (iord == 1) // 1-point Gauss
110 : {
111 0 : sg2.resize(1);
112 0 : sg2[0].resize(3);
113 0 : sg2[0][0] = 0.0;
114 0 : sg2[0][1] = 0.0;
115 0 : sg2[0][2] = 4.0;
116 : }
117 204 : else if (iord == 2) // 2x2-point Gauss
118 : {
119 204 : sg2.resize(4);
120 1020 : for (unsigned int i = 0; i < 4; ++i)
121 816 : sg2[i].resize(3);
122 1020 : for (unsigned int i = 0; i < 4; ++i)
123 : {
124 816 : sg2[i][0] = (1 / sqrt(3)) * lr4[i];
125 816 : sg2[i][1] = (1 / sqrt(3)) * lz4[i];
126 816 : sg2[i][2] = 1.0;
127 : }
128 : }
129 0 : else if (iord == 3) // 3x3-point Gauss
130 : {
131 0 : sg2.resize(9);
132 0 : for (unsigned int i = 0; i < 9; ++i)
133 0 : sg2[i].resize(3);
134 0 : for (unsigned int i = 0; i < 9; ++i)
135 : {
136 0 : sg2[i][0] = sqrt(0.6) * lr9[i];
137 0 : sg2[i][1] = sqrt(0.6) * lz9[i];
138 0 : sg2[i][2] = (1.0 / 81.0) * lw9[i];
139 : }
140 : }
141 : else
142 0 : mooseError("Invalid quadrature order = " + Moose::stringify(iord) + " for quad elements");
143 : }
144 584852 : else if (nen == 3) // triangle
145 : {
146 584852 : if (iord == 1) // one-point Gauss
147 : {
148 584420 : sg2.resize(1);
149 584420 : sg2[0].resize(4);
150 584420 : sg2[0][0] = 1.0 / 3.0;
151 584420 : sg2[0][1] = 1.0 / 3.0;
152 584420 : sg2[0][2] = 1.0 / 3.0;
153 584420 : sg2[0][3] = 0.5;
154 : }
155 432 : else if (iord == 2) // three-point Gauss
156 : {
157 432 : sg2.resize(3);
158 1728 : for (unsigned int i = 0; i < 3; ++i)
159 1296 : sg2[i].resize(4);
160 432 : sg2[0][0] = 2.0 / 3.0;
161 432 : sg2[0][1] = 1.0 / 6.0;
162 432 : sg2[0][2] = 1.0 / 6.0;
163 432 : sg2[0][3] = 1.0 / 6.0;
164 432 : sg2[1][0] = 1.0 / 6.0;
165 432 : sg2[1][1] = 2.0 / 3.0;
166 432 : sg2[1][2] = 1.0 / 6.0;
167 432 : sg2[1][3] = 1.0 / 6.0;
168 432 : sg2[2][0] = 1.0 / 6.0;
169 432 : sg2[2][1] = 1.0 / 6.0;
170 432 : sg2[2][2] = 2.0 / 3.0;
171 432 : sg2[2][3] = 1.0 / 6.0;
172 : }
173 0 : else if (iord == 3) // four-point Gauss
174 : {
175 0 : sg2.resize(4);
176 0 : for (unsigned int i = 0; i < 4; ++i)
177 0 : sg2[i].resize(4);
178 0 : sg2[0][0] = 1.5505102572168219018027159252941e-01;
179 0 : sg2[0][1] = 1.7855872826361642311703513337422e-01;
180 0 : sg2[0][2] = 1.0 - sg2[0][0] - sg2[0][1];
181 0 : sg2[0][3] = 1.5902069087198858469718450103758e-01;
182 :
183 0 : sg2[1][0] = 6.4494897427831780981972840747059e-01;
184 0 : sg2[1][1] = 7.5031110222608118177475598324603e-02;
185 0 : sg2[1][2] = 1.0 - sg2[1][0] - sg2[1][1];
186 0 : sg2[1][3] = 9.0979309128011415302815498962418e-02;
187 :
188 0 : sg2[2][0] = 1.5505102572168219018027159252941e-01;
189 0 : sg2[2][1] = 6.6639024601470138670269327409637e-01;
190 0 : sg2[2][2] = 1.0 - sg2[2][0] - sg2[2][1];
191 0 : sg2[2][3] = 1.5902069087198858469718450103758e-01;
192 :
193 0 : sg2[3][0] = 6.4494897427831780981972840747059e-01;
194 0 : sg2[3][1] = 2.8001991549907407200279599420481e-01;
195 0 : sg2[3][2] = 1.0 - sg2[3][0] - sg2[3][1];
196 0 : sg2[3][3] = 9.0979309128011415302815498962418e-02;
197 : }
198 0 : else if (iord == 4) // six-point Guass
199 : {
200 : const unsigned int n_wts = 2;
201 0 : const Real wts[n_wts] = {1.1169079483900573284750350421656140e-01L,
202 : 5.4975871827660933819163162450105264e-02L};
203 :
204 0 : const Real a[n_wts] = {4.4594849091596488631832925388305199e-01L,
205 : 9.1576213509770743459571463402201508e-02L};
206 :
207 0 : const Real b[n_wts] = {0., 0.}; // not used
208 0 : const unsigned int permutation_ids[n_wts] = {3, 3};
209 :
210 : std::vector<Point> points;
211 : std::vector<Real> weights;
212 0 : dunavant_rule2(wts, a, b, permutation_ids, n_wts, points, weights); // 6 total points
213 :
214 0 : sg2.resize(6);
215 0 : for (unsigned int i = 0; i < 6; ++i)
216 0 : sg2[i].resize(4);
217 0 : for (unsigned int i = 0; i < 6; ++i)
218 : {
219 0 : sg2[i][0] = points[i](0);
220 0 : sg2[i][1] = points[i](1);
221 0 : sg2[i][2] = 1.0 - points[i](0) - points[i](1);
222 0 : sg2[i][3] = weights[i];
223 : }
224 0 : }
225 : else
226 0 : mooseError("Invalid quadrature order = " + Moose::stringify(iord) + " for triangle elements");
227 : }
228 : else
229 0 : mooseError("Invalid 2D element type");
230 585056 : }
231 :
232 : void
233 0 : wissmannPoints(unsigned int nqp, std::vector<std::vector<Real>> & wss)
234 : {
235 0 : if (nqp == 6)
236 : {
237 0 : wss.resize(6);
238 0 : for (unsigned int i = 0; i < 6; ++i)
239 0 : wss[i].resize(3);
240 0 : wss[0][0] = 0.0;
241 0 : wss[0][1] = 0.0;
242 0 : wss[0][2] = 1.1428571428571428;
243 :
244 0 : wss[1][0] = 0.0;
245 0 : wss[1][1] = 9.6609178307929590e-01;
246 0 : wss[1][2] = 4.3956043956043956e-01;
247 :
248 0 : wss[2][0] = 8.5191465330460049e-01;
249 0 : wss[2][1] = 4.5560372783619284e-01;
250 0 : wss[2][2] = 5.6607220700753210e-01;
251 :
252 0 : wss[3][0] = -wss[2][0];
253 0 : wss[3][1] = wss[2][1];
254 0 : wss[3][2] = wss[2][2];
255 :
256 0 : wss[4][0] = 6.3091278897675402e-01;
257 0 : wss[4][1] = -7.3162995157313452e-01;
258 0 : wss[4][2] = 6.4271900178367668e-01;
259 :
260 0 : wss[5][0] = -wss[4][0];
261 0 : wss[5][1] = wss[4][1];
262 0 : wss[5][2] = wss[4][2];
263 : }
264 : else
265 0 : mooseError("Unknown Wissmann quadrature type");
266 0 : }
267 :
268 : void
269 587960 : shapeFunc2D(unsigned int nen,
270 : std::vector<Real> & ss,
271 : std::vector<Point> & xl,
272 : std::vector<std::vector<Real>> & shp,
273 : Real & xsj,
274 : bool natl_flg)
275 : {
276 : // Get shape functions and derivatives
277 587960 : Real s[4] = {-0.5, 0.5, 0.5, -0.5};
278 587960 : Real t[4] = {-0.5, -0.5, 0.5, 0.5};
279 :
280 587960 : if (nen == 4) // quad element
281 : {
282 2244 : Real xs[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
283 : Real sx[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
284 11220 : for (unsigned int i = 0; i < 4; ++i)
285 : {
286 8976 : shp[i][2] = (0.5 + s[i] * ss[0]) * (0.5 + t[i] * ss[1]);
287 8976 : shp[i][0] = s[i] * (0.5 + t[i] * ss[1]);
288 8976 : shp[i][1] = t[i] * (0.5 + s[i] * ss[0]);
289 : }
290 6732 : for (unsigned int i = 0; i < 2; ++i) // x, y
291 : {
292 13464 : for (unsigned int j = 0; j < 2; ++j) // xi, eta
293 : {
294 8976 : xs[i][j] = 0.0;
295 44880 : for (unsigned int k = 0; k < nen; ++k)
296 35904 : xs[i][j] += xl[k](i) * shp[k][j];
297 : }
298 : }
299 2244 : xsj = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; // det(j)
300 2244 : if (natl_flg == false) // get global derivatives
301 : {
302 0 : Real temp = 1.0 / xsj;
303 0 : sx[0][0] = xs[1][1] * temp; // inv(j)
304 0 : sx[1][1] = xs[0][0] * temp;
305 0 : sx[0][1] = -xs[0][1] * temp;
306 0 : sx[1][0] = -xs[1][0] * temp;
307 0 : for (unsigned int i = 0; i < nen; ++i)
308 : {
309 0 : temp = shp[i][0] * sx[0][0] + shp[i][1] * sx[1][0];
310 0 : shp[i][1] = shp[i][0] * sx[0][1] + shp[i][1] * sx[1][1];
311 0 : shp[i][0] = temp;
312 : }
313 : }
314 : }
315 585716 : else if (nen == 3) // triangle element
316 : {
317 : // x1*(y2 - y3) + x2*(y3 - y1) + x3*(y1 - y2)
318 : Point x13 = xl[2] - xl[0];
319 : Point x23 = xl[2] - xl[1];
320 585716 : Point cross_prod = x13.cross(x23);
321 585716 : xsj = cross_prod.norm();
322 : Real xsjr = 1.0;
323 585716 : if (xsj != 0.0)
324 585716 : xsjr = 1.0 / xsj;
325 : // xsj *= 0.5; // we do not have this 0.5 here because in stdQuad2D the sum of all weights in
326 : // tri is 0.5
327 585716 : shp[0][2] = ss[0];
328 585716 : shp[1][2] = ss[1];
329 585716 : shp[2][2] = ss[2];
330 585716 : if (natl_flg == false) // need global drivatives
331 : {
332 0 : shp[0][0] = (xl[1](1) - xl[2](1)) * xsjr;
333 0 : shp[0][1] = (xl[2](0) - xl[1](0)) * xsjr;
334 0 : shp[1][0] = (xl[2](1) - xl[0](1)) * xsjr;
335 0 : shp[1][1] = (xl[0](0) - xl[2](0)) * xsjr;
336 0 : shp[2][0] = (xl[0](1) - xl[1](1)) * xsjr;
337 0 : shp[2][1] = (xl[1](0) - xl[0](0)) * xsjr;
338 : }
339 : else
340 : {
341 585716 : shp[0][0] = 1.0;
342 585716 : shp[0][1] = 0.0;
343 585716 : shp[1][0] = 0.0;
344 585716 : shp[1][1] = 1.0;
345 585716 : shp[2][0] = -1.0;
346 585716 : shp[2][1] = -1.0;
347 : }
348 : }
349 : else
350 0 : mooseError("ShapeFunc2D only works for linear quads and tris!");
351 587960 : }
352 :
353 : double
354 10440064 : r8vec_norm(int n, double a[])
355 : {
356 : // John Burkardt geometry.cpp
357 : double v = 0.0;
358 41760256 : for (int i = 0; i < n; ++i)
359 31320192 : v = v + a[i] * a[i];
360 10440064 : v = std::sqrt(v);
361 10440064 : return v;
362 : }
363 :
364 : void
365 0 : r8vec_copy(int n, double a1[], double a2[])
366 : {
367 : // John Burkardt geometry.cpp
368 0 : for (int i = 0; i < n; ++i)
369 0 : a2[i] = a1[i];
370 0 : return;
371 : }
372 :
373 : bool
374 5220032 : r8vec_eq(int n, double a1[], double a2[])
375 : {
376 : // John Burkardt geometry.cpp
377 10349776 : for (int i = 0; i < n; ++i)
378 10349776 : if (a1[i] != a2[i])
379 : return false;
380 : return true;
381 : }
382 :
383 : double
384 5220032 : r8vec_dot_product(int n, double a1[], double a2[])
385 : {
386 : // John Burkardt geometry.cpp
387 : double value = 0.0;
388 20880128 : for (int i = 0; i < n; ++i)
389 15660096 : value += a1[i] * a2[i];
390 5220032 : return value;
391 : }
392 :
393 : bool
394 5220032 : line_exp_is_degenerate_nd(int dim_num, double p1[], double p2[])
395 : {
396 : // John Burkardt geometry.cpp
397 : bool value;
398 5220032 : value = r8vec_eq(dim_num, p1, p2);
399 5220032 : return value;
400 : }
401 :
402 : int
403 5220032 : plane_normal_line_exp_int_3d(
404 : double pp[3], double normal[3], double p1[3], double p2[3], double pint[3])
405 : {
406 : // John Burkardt geometry.cpp
407 : // Parameters:
408 : //
409 : // Input, double PP[3], a point on the plane.
410 : //
411 : // Input, double NORMAL[3], a normal vector to the plane.
412 : //
413 : // Input, double P1[3], P2[3], two distinct points on the line.
414 : //
415 : // Output, double PINT[3], the coordinates of a
416 : // common point of the plane and line, when IVAL is 1 or 2.
417 : //
418 : // Output, integer PLANE_NORMAL_LINE_EXP_INT_3D, the kind of intersection;
419 : // 0, the line and plane seem to be parallel and separate;
420 : // 1, the line and plane intersect at a single point;
421 : // 2, the line and plane seem to be parallel and joined.
422 : #define DIM_NUM 3
423 :
424 : double direction[DIM_NUM];
425 : int ival;
426 : double temp;
427 : double temp2;
428 : //
429 : // Make sure the line is not degenerate.
430 5220032 : if (line_exp_is_degenerate_nd(DIM_NUM, p1, p2))
431 0 : mooseError("PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error! The line is degenerate.");
432 : //
433 : // Make sure the plane normal vector is a unit vector.
434 5220032 : temp = r8vec_norm(DIM_NUM, normal);
435 5220032 : if (temp == 0.0)
436 0 : mooseError("PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error! The normal vector of the plane is "
437 : "degenerate.");
438 :
439 20880128 : for (unsigned int i = 0; i < DIM_NUM; ++i)
440 15660096 : normal[i] = normal[i] / temp;
441 : //
442 : // Determine the unit direction vector of the line.
443 20880128 : for (unsigned int i = 0; i < DIM_NUM; ++i)
444 15660096 : direction[i] = p2[i] - p1[i];
445 5220032 : temp = r8vec_norm(DIM_NUM, direction);
446 :
447 20880128 : for (unsigned int i = 0; i < DIM_NUM; ++i)
448 15660096 : direction[i] = direction[i] / temp;
449 : //
450 : // If the normal and direction vectors are orthogonal, then
451 : // we have a special case to deal with.
452 5220032 : if (r8vec_dot_product(DIM_NUM, normal, direction) == 0.0)
453 : {
454 : temp = 0.0;
455 3517536 : for (unsigned int i = 0; i < DIM_NUM; ++i)
456 2638152 : temp = temp + normal[i] * (p1[i] - pp[i]);
457 :
458 879384 : if (temp == 0.0)
459 : {
460 : ival = 2;
461 0 : r8vec_copy(DIM_NUM, p1, pint);
462 : }
463 : else
464 : {
465 : ival = 0;
466 3517536 : for (unsigned int i = 0; i < DIM_NUM; ++i)
467 2638152 : pint[i] = 1.0e20; // dummy huge value
468 : }
469 879384 : return ival;
470 : }
471 : //
472 : // Determine the distance along the direction vector to the intersection point.
473 : temp = 0.0;
474 17362592 : for (unsigned int i = 0; i < DIM_NUM; ++i)
475 13021944 : temp = temp + normal[i] * (pp[i] - p1[i]);
476 : temp2 = 0.0;
477 17362592 : for (unsigned int i = 0; i < DIM_NUM; ++i)
478 13021944 : temp2 = temp2 + normal[i] * direction[i];
479 :
480 : ival = 1;
481 17362592 : for (unsigned int i = 0; i < DIM_NUM; ++i)
482 13021944 : pint[i] = p1[i] + temp * direction[i] / temp2;
483 :
484 : return ival;
485 : #undef DIM_NUM
486 : }
487 :
488 : double
489 54720 : polyhedron_volume_3d(
490 : double coord[], int order_max, int face_num, int node[], int /*node_num*/, int order[])
491 : //
492 : // Purpose:
493 : //
494 : // POLYHEDRON_VOLUME_3D computes the volume of a polyhedron in 3D.
495 : //
496 : // Licensing:
497 : //
498 : // This code is distributed under the GNU LGPL license.
499 : //
500 : // Modified:
501 : //
502 : // 21 August 2003
503 : //
504 : // Author:
505 : //
506 : // John Burkardt
507 : //
508 : // Parameters:
509 : //
510 : // Input, double COORD[NODE_NUM*3], the 3D coordinates of the vertices.
511 : // The vertices may be listed in any order.
512 : //
513 : // Input, int ORDER_MAX, the maximum number of vertices that make
514 : // up a face of the polyhedron.
515 : //
516 : // Input, int FACE_NUM, the number of faces of the polyhedron.
517 : //
518 : // Input, int NODE[FACE_NUM*ORDER_MAX]. Face I is defined by
519 : // the vertices NODE(I,1) through NODE(I,ORDER(I)). These vertices
520 : // are listed in neighboring order.
521 : //
522 : // Input, int NODE_NUM, the number of points stored in COORD.
523 : //
524 : // Input, int ORDER[FACE_NUM], the number of vertices making up
525 : // each face.
526 : //
527 : // Output, double POLYHEDRON_VOLUME_3D, the volume of the polyhedron.
528 : //
529 : {
530 : #define DIM_NUM 3
531 :
532 : int face;
533 : int n1;
534 : int n2;
535 : int n3;
536 : double term;
537 : int v;
538 : double volume;
539 : double x1;
540 : double x2;
541 : double x3;
542 : double y1;
543 : double y2;
544 : double y3;
545 : double z1;
546 : double z2;
547 : double z3;
548 : //
549 : volume = 0.0;
550 : //
551 : // Triangulate each face.
552 : //
553 330581 : for (face = 0; face < face_num; face++)
554 : {
555 275861 : n3 = node[order[face] - 1 + face * order_max];
556 275861 : x3 = coord[0 + n3 * 3];
557 275861 : y3 = coord[1 + n3 * 3];
558 275861 : z3 = coord[2 + n3 * 3];
559 :
560 722649 : for (v = 0; v < order[face] - 2; v++)
561 : {
562 446788 : n1 = node[v + face * order_max];
563 446788 : x1 = coord[0 + n1 * 3];
564 446788 : y1 = coord[1 + n1 * 3];
565 446788 : z1 = coord[2 + n1 * 3];
566 :
567 446788 : n2 = node[v + 1 + face * order_max];
568 446788 : x2 = coord[0 + n2 * 3];
569 446788 : y2 = coord[1 + n2 * 3];
570 446788 : z2 = coord[2 + n2 * 3];
571 :
572 446788 : term =
573 446788 : x1 * y2 * z3 - x1 * y3 * z2 + x2 * y3 * z1 - x2 * y1 * z3 + x3 * y1 * z2 - x3 * y2 * z1;
574 :
575 446788 : volume = volume + term;
576 : }
577 : }
578 :
579 54720 : volume = volume / 6.0;
580 :
581 54720 : return volume;
582 : #undef DIM_NUM
583 : }
584 :
585 : void
586 54720 : i4vec_zero(int n, int a[])
587 : //
588 : // Purpose:
589 : //
590 : // I4VEC_ZERO zeroes an I4VEC.
591 : //
592 : // Licensing:
593 : //
594 : // This code is distributed under the GNU LGPL license.
595 : //
596 : // Modified:
597 : //
598 : // 01 August 2005
599 : //
600 : // Author:
601 : //
602 : // John Burkardt
603 : //
604 : // Parameters:
605 : //
606 : // Input, int N, the number of entries in the vector.
607 : //
608 : // Output, int A[N], a vector of zeroes.
609 : //
610 : {
611 : int i;
612 :
613 1108536 : for (i = 0; i < n; i++)
614 : {
615 1053816 : a[i] = 0;
616 : }
617 54720 : return;
618 : }
619 :
620 : void
621 26667374 : normalizePoint(Point & p)
622 : {
623 26667374 : Real len = p.norm();
624 26667374 : if (len > tol)
625 26518614 : p = (1.0 / len) * p;
626 : else
627 : p.zero();
628 26667374 : }
629 :
630 : void
631 11194 : normalizePoint(EFAPoint & p)
632 : {
633 11194 : Real len = p.norm();
634 11194 : if (len > tol)
635 11194 : p *= (1.0 / len);
636 11194 : }
637 :
638 : double
639 0 : r8_acos(double c)
640 : //
641 : // Purpose:
642 : //
643 : // R8_ACOS computes the arc cosine function, with argument truncation.
644 : //
645 : // Discussion:
646 : //
647 : // If you call your system ACOS routine with an input argument that is
648 : // outside the range [-1.0, 1.0 ], you may get an unpleasant surprise.
649 : // This routine truncates arguments outside the range.
650 : //
651 : // Licensing:
652 : //
653 : // This code is distributed under the GNU LGPL license.
654 : //
655 : // Modified:
656 : //
657 : // 13 June 2002
658 : //
659 : // Author:
660 : //
661 : // John Burkardt
662 : //
663 : // Parameters:
664 : //
665 : // Input, double C, the argument, the cosine of an angle.
666 : //
667 : // Output, double R8_ACOS, an angle whose cosine is C.
668 : //
669 : {
670 : #define PI 3.141592653589793
671 :
672 : double value;
673 :
674 0 : if (c <= -1.0)
675 : {
676 : value = PI;
677 : }
678 0 : else if (1.0 <= c)
679 : {
680 : value = 0.0;
681 : }
682 : else
683 : {
684 0 : value = acos(c);
685 : }
686 0 : return value;
687 : #undef PI
688 : }
689 :
690 : double
691 0 : angle_rad_3d(double p1[3], double p2[3], double p3[3])
692 : //
693 : // Purpose:
694 : //
695 : // ANGLE_RAD_3D returns the angle between two vectors in 3D.
696 : //
697 : // Discussion:
698 : //
699 : // The routine always computes the SMALLER of the two angles between
700 : // two vectors. Thus, if the vectors make an (exterior) angle of 200
701 : // degrees, the (interior) angle of 160 is reported.
702 : //
703 : // X dot Y = Norm(X) * Norm(Y) * Cos ( Angle(X,Y) )
704 : //
705 : // Licensing:
706 : //
707 : // This code is distributed under the GNU LGPL license.
708 : //
709 : // Modified:
710 : //
711 : // 20 June 2005
712 : //
713 : // Author:
714 : //
715 : // John Burkardt
716 : //
717 : // Parameters:
718 : //
719 : // Input, double P1[3], P2[3], P3[3], points defining an angle.
720 : // The rays are P1 - P2 and P3 - P2.
721 : //
722 : // Output, double ANGLE_RAD_3D, the angle between the two vectors, in radians.
723 : // This value will always be between 0 and PI. If either vector has
724 : // zero length, then the angle is returned as zero.
725 : //
726 : {
727 : #define DIM_NUM 3
728 :
729 : double dot;
730 : int i;
731 : double v1norm;
732 : double v2norm;
733 : double value;
734 :
735 : v1norm = 0.0;
736 0 : for (i = 0; i < DIM_NUM; i++)
737 : {
738 0 : v1norm = v1norm + pow(p1[i] - p2[i], 2);
739 : }
740 0 : v1norm = sqrt(v1norm);
741 :
742 0 : if (v1norm == 0.0)
743 : {
744 : value = 0.0;
745 : return value;
746 : }
747 :
748 : v2norm = 0.0;
749 0 : for (i = 0; i < DIM_NUM; i++)
750 : {
751 0 : v2norm = v2norm + pow(p3[i] - p2[i], 2);
752 : }
753 0 : v2norm = sqrt(v2norm);
754 :
755 0 : if (v2norm == 0.0)
756 : {
757 : value = 0.0;
758 : return value;
759 : }
760 :
761 : dot = 0.0;
762 0 : for (i = 0; i < DIM_NUM; i++)
763 : {
764 0 : dot = dot + (p1[i] - p2[i]) * (p3[i] - p2[i]);
765 : }
766 :
767 0 : value = r8_acos(dot / (v1norm * v2norm));
768 :
769 0 : return value;
770 : #undef DIM_NUM
771 : }
772 :
773 : bool
774 6158340 : intersectSegmentWithCutLine(const Point & segment_point1,
775 : const Point & segment_point2,
776 : const std::pair<Point, Point> & cutting_line_points,
777 : const Real & cutting_line_fraction,
778 : Real & segment_intersection_fraction)
779 : {
780 : // Use the algorithm described here to determine whether a line segment is intersected
781 : // by a cutting line, and to compute the fraction along that line where the intersection
782 : // occurs:
783 : // http://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
784 :
785 : bool cut_segment = false;
786 : Point seg_dir = segment_point2 - segment_point1;
787 : Point cut_dir = cutting_line_points.second - cutting_line_points.first;
788 : Point cut_start_to_seg_start = segment_point1 - cutting_line_points.first;
789 :
790 6158340 : Real cut_dir_cross_seg_dir = crossProduct2D(cut_dir, seg_dir);
791 :
792 6158340 : if (std::abs(cut_dir_cross_seg_dir) > Xfem::tol)
793 : {
794 : // Fraction of the distance along the cutting segment where it intersects the edge segment
795 5172768 : Real cut_int_frac = crossProduct2D(cut_start_to_seg_start, seg_dir) / cut_dir_cross_seg_dir;
796 :
797 5172768 : if (cut_int_frac >= 0.0 && cut_int_frac <= cutting_line_fraction)
798 : { // Cutting segment intersects the line of the edge segment, but the intersection point may be
799 : // outside the segment
800 318939 : Real int_frac = crossProduct2D(cut_start_to_seg_start, cut_dir) / cut_dir_cross_seg_dir;
801 318939 : if (int_frac >= 0.0 && int_frac <= 1.0)
802 : {
803 : cut_segment = true;
804 57736 : segment_intersection_fraction = int_frac;
805 : }
806 : }
807 : }
808 6158340 : return cut_segment;
809 : }
810 :
811 : Real
812 11650047 : crossProduct2D(const Point & point_a, const Point & point_b)
813 : {
814 11650047 : return (point_a(0) * point_b(1) - point_b(0) * point_a(1));
815 : }
816 :
817 : Real
818 66649022 : pointSegmentDistance(const Point & x0, const Point & x1, const Point & x2, Point & xp)
819 : {
820 : Point dx = x2 - x1;
821 : Real m2 = dx * dx;
822 66649022 : if (m2 == 0)
823 0 : mooseError("In XFEMFuncs::pointSegmentDistance(), x0 and x1 should be two different points.");
824 : // find parameter coordinate of closest point on segment
825 66649022 : Real s12 = (x2 - x0) * dx / m2;
826 66649022 : if (s12 < 0)
827 : s12 = 0;
828 45089819 : else if (s12 > 1)
829 : s12 = 1;
830 : // and find the distance
831 66649022 : xp = s12 * x1 + (1 - s12) * x2;
832 66649022 : return std::sqrt((x0 - xp) * (x0 - xp));
833 : }
834 :
835 : Real
836 33611052 : pointTriangleDistance(const Point & x0,
837 : const Point & x1,
838 : const Point & x2,
839 : const Point & x3,
840 : Point & xp,
841 : unsigned int & region)
842 : {
843 : Point x13 = x1 - x3, x23 = x2 - x3, x03 = x0 - x3;
844 : Real m13 = x13 * x13, m23 = x23 * x23, d = x13 * x23;
845 33611052 : Real invdet = 1.0 / std::max(m13 * m23 - d * d, 1e-30);
846 : Real a = x13 * x03, b = x23 * x03;
847 :
848 33611052 : Real w23 = invdet * (m23 * a - d * b);
849 33611052 : Real w31 = invdet * (m13 * b - d * a);
850 33611052 : Real w12 = 1 - w23 - w31;
851 33611052 : if (w23 >= 0 && w31 >= 0 && w12 >= 0)
852 : { // if we're inside the triangle
853 286541 : region = 0;
854 286541 : xp = w23 * x1 + w31 * x2 + w12 * x3;
855 286541 : return std::sqrt((x0 - xp) * (x0 - xp));
856 : }
857 : else
858 : {
859 16827898 : if (w23 > 0) // this rules out edge 2-3 for us
860 : {
861 : Point xp1, xp2;
862 16827898 : Real distance_12 = pointSegmentDistance(x0, x1, x2, xp1);
863 16827898 : Real distance_13 = pointSegmentDistance(x0, x1, x3, xp2);
864 16827898 : Real distance_1 = std::sqrt((x0 - x1) * (x0 - x1));
865 16827898 : if (std::min(distance_12, distance_13) < distance_1)
866 : {
867 5879152 : if (distance_12 < distance_13)
868 : {
869 1439890 : region = 4;
870 1439890 : xp = xp1;
871 1439890 : return distance_12;
872 : }
873 : else
874 : {
875 4439262 : region = 6;
876 4439262 : xp = xp2;
877 4439262 : return distance_13;
878 : }
879 : }
880 : else
881 : {
882 10948746 : region = 1;
883 10948746 : xp = x1;
884 10948746 : return distance_1;
885 : }
886 : }
887 16496613 : else if (w31 > 0) // this rules out edge 1-3
888 : {
889 : Point xp1, xp2;
890 12180003 : Real distance_12 = pointSegmentDistance(x0, x1, x2, xp1);
891 12180003 : Real distance_23 = pointSegmentDistance(x0, x2, x3, xp2);
892 12180003 : Real distance_2 = std::sqrt((x0 - x2) * (x0 - x2));
893 12180003 : if (std::min(distance_12, distance_23) < distance_2)
894 : {
895 4124766 : if (distance_12 < distance_23)
896 : {
897 0 : region = 4;
898 0 : xp = xp1;
899 0 : return distance_12;
900 : }
901 : else
902 : {
903 4124766 : region = 5;
904 4124766 : xp = xp2;
905 4124766 : return distance_23;
906 : }
907 : }
908 : else
909 : {
910 8055237 : region = 2;
911 8055237 : xp = x2;
912 8055237 : return distance_2;
913 : }
914 : }
915 : else // w12 must be >0, ruling out edge 1-2
916 : {
917 : Point xp1, xp2;
918 4316610 : Real distance_23 = pointSegmentDistance(x0, x2, x3, xp1);
919 4316610 : Real distance_31 = pointSegmentDistance(x0, x3, x1, xp2);
920 4316610 : Real distance_3 = std::sqrt((x0 - x3) * (x0 - x3));
921 4316610 : if (std::min(distance_23, distance_31) < distance_3)
922 : {
923 0 : if (distance_23 < distance_31)
924 : {
925 0 : region = 5;
926 0 : xp = xp1;
927 0 : return distance_23;
928 : }
929 : else
930 : {
931 0 : region = 6;
932 0 : xp = xp2;
933 0 : return distance_31;
934 : }
935 : }
936 : else
937 : {
938 4316610 : region = 3;
939 4316610 : xp = x3;
940 4316610 : return distance_3;
941 : }
942 : }
943 : }
944 : mooseError("Cannot find closest location in XFEMFuncs::pointTriangleDistance().");
945 : }
946 :
947 : bool
948 4210272 : intersectWithEdge(const Point & p1,
949 : const Point & p2,
950 : const std::vector<Point> & vertices,
951 : Point & pint)
952 : {
953 : bool has_intersection = false;
954 :
955 4210272 : if (vertices.size() != 3)
956 0 : mooseError("The number of vertices of cutting element must be 3.");
957 :
958 4210272 : Plane elem_plane(vertices[0], vertices[1], vertices[2]);
959 4210272 : Point point = vertices[0];
960 4210272 : Point normal = elem_plane.unit_normal(point);
961 :
962 4210272 : std::array<Real, 3> plane_point = {{point(0), point(1), point(2)}};
963 4210272 : std::array<Real, 3> planenormal = {{normal(0), normal(1), normal(2)}};
964 4210272 : std::array<Real, 3> edge_point1 = {{p1(0), p1(1), p1(2)}};
965 4210272 : std::array<Real, 3> edge_point2 = {{p2(0), p2(1), p2(2)}};
966 4210272 : std::array<Real, 3> cut_point = {{0.0, 0.0, 0.0}};
967 :
968 4210272 : if (Xfem::plane_normal_line_exp_int_3d(
969 : &plane_point[0], &planenormal[0], &edge_point1[0], &edge_point2[0], &cut_point[0]) == 1)
970 : {
971 3865840 : Point temp_p(cut_point[0], cut_point[1], cut_point[2]);
972 3865840 : if (isInsideCutPlane(vertices, temp_p) && isInsideEdge(p1, p2, temp_p))
973 : {
974 724 : pint = temp_p;
975 : has_intersection = true;
976 : }
977 : }
978 :
979 4210272 : return has_intersection;
980 4210272 : }
981 :
982 : bool
983 8828 : isInsideEdge(const Point & p1, const Point & p2, const Point & p)
984 : {
985 : Real dotp1 = (p1 - p) * (p2 - p1);
986 : Real dotp2 = (p2 - p) * (p2 - p1);
987 8828 : return (dotp1 * dotp2 <= 0.0);
988 : }
989 :
990 : Real
991 672 : getRelativePosition(const Point & p1, const Point & p2, const Point & p)
992 : {
993 672 : Real full_len = (p2 - p1).norm();
994 672 : Real len_p1_p = (p - p1).norm();
995 672 : return len_p1_p / full_len;
996 : }
997 :
998 : bool
999 3865840 : isInsideCutPlane(const std::vector<Point> & vertices, const Point & p)
1000 : {
1001 3865840 : unsigned int n_node = vertices.size();
1002 :
1003 3865840 : if (n_node != 3)
1004 0 : mooseError("The number of vertices of cutting element must be 3.");
1005 :
1006 3865840 : Plane elem_plane(vertices[0], vertices[1], vertices[2]);
1007 3865840 : Point normal = elem_plane.unit_normal(vertices[0]);
1008 :
1009 : bool inside = false;
1010 : unsigned int counter = 0;
1011 :
1012 15463360 : for (unsigned int i = 0; i < n_node; ++i)
1013 : {
1014 11597520 : unsigned int iplus1 = (i < n_node - 1 ? i + 1 : 0);
1015 11597520 : Point middle2p = p - 0.5 * (vertices[i] + vertices[iplus1]);
1016 : const Point side_tang = vertices[iplus1] - vertices[i];
1017 11597520 : Point side_norm = side_tang.cross(normal);
1018 :
1019 11597520 : normalizePoint(middle2p);
1020 11597520 : normalizePoint(side_norm);
1021 :
1022 11597520 : if (middle2p * side_norm <= 0)
1023 6380708 : counter += 1;
1024 : }
1025 :
1026 3865840 : if (counter == n_node)
1027 : inside = true;
1028 3865840 : return inside;
1029 3865840 : }
1030 :
1031 : } // namespace XFEM
|
