XFEMFuncs.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7
8 /****************************************************************/
9 /* This file also contains modified functions from libMesh and */
10 /* the geometry library of John Burkardt */
11 /* http://people.sc.fsu.edu/~jburkardt/ */
13 /****************************************************************/
14
15 #include "XFEMFuncs.h"
16
17 #include "MooseError.h"
18 #include "Conversion.h"
19
20 using namespace libMesh;
21
22 namespace Xfem
23 {
24
25 void
26 dunavant_rule2(const Real * wts,
27  const Real * a,
28  const Real * b,
29  const unsigned int * permutation_ids,
30  unsigned int n_wts,
31  std::vector<Point> & points,
32  std::vector<Real> & weights)
33 {
35  // Figure out how many total points by summing up the entries
36  // in the permutation_ids array, and resize the _points and _weights
37  // vectors appropriately.
38  unsigned int total_pts = 0;
39  for (unsigned int p = 0; p < n_wts; ++p)
40  total_pts += permutation_ids[p];
41
42  // Resize point and weight vectors appropriately.
43  points.resize(total_pts);
44  weights.resize(total_pts);
45
46  // Always insert into the points & weights vector relative to the offset
47  unsigned int offset = 0;
48  for (unsigned int p = 0; p < n_wts; ++p)
49  {
50  switch (permutation_ids[p])
51  {
52  case 1:
53  {
54  // The point has only a single permutation (the centroid!)
55  // So we don't even need to look in the a or b arrays.
56  points[offset + 0] = Point(1.0L / 3.0L, 1.0L / 3.0L);
57  weights[offset + 0] = wts[p];
58
59  offset += 1;
60  break;
61  }
62
63  case 3:
64  {
65  // For this type of rule, don't need to look in the b array.
66  points[offset + 0] = Point(a[p], a[p]); // (a,a)
67  points[offset + 1] = Point(a[p], 1.L - 2.L * a[p]); // (a,1-2a)
68  points[offset + 2] = Point(1.L - 2.L * a[p], a[p]); // (1-2a,a)
69
70  for (unsigned int j = 0; j < 3; ++j)
71  weights[offset + j] = wts[p];
72
73  offset += 3;
74  break;
75  }
76
77  case 6:
78  {
79  // This type of point uses all 3 arrays...
80  points[offset + 0] = Point(a[p], b[p]);
81  points[offset + 1] = Point(b[p], a[p]);
82  points[offset + 2] = Point(a[p], 1.L - a[p] - b[p]);
83  points[offset + 3] = Point(1.L - a[p] - b[p], a[p]);
84  points[offset + 4] = Point(b[p], 1.L - a[p] - b[p]);
85  points[offset + 5] = Point(1.L - a[p] - b[p], b[p]);
86
87  for (unsigned int j = 0; j < 6; ++j)
88  weights[offset + j] = wts[p];
89
90  offset += 6;
91  break;
92  }
93
94  default:
95  mooseError("Unknown permutation id: ", permutation_ids[p], "!");
96  }
97  }
98 }
99
100 void
101 stdQuadr2D(unsigned int nen, unsigned int iord, std::vector<std::vector<Real>> & sg2)
102 {
103  // Purpose: get Guass integration points for 2D quad and tri elems
104  // N.B. only works for n_qp <= 6
105
106  Real lr4[4] = {-1.0, 1.0, -1.0, 1.0}; // libmesh order
107  Real lz4[4] = {-1.0, -1.0, 1.0, 1.0};
108  Real lr9[9] = {-1.0, 0.0, 1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0}; // libmesh order
109  Real lz9[9] = {-1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0};
110  Real lw9[9] = {25.0, 40.0, 25.0, 40.0, 64.0, 40.0, 25.0, 40.0, 25.0};
111
112  if (nen == 4) // 2d quad element
113  {
114  if (iord == 1) // 1-point Gauss
115  {
116  sg2.resize(1);
117  sg2[0].resize(3);
118  sg2[0][0] = 0.0;
119  sg2[0][1] = 0.0;
120  sg2[0][2] = 4.0;
121  }
122  else if (iord == 2) // 2x2-point Gauss
123  {
124  sg2.resize(4);
125  for (unsigned int i = 0; i < 4; ++i)
126  sg2[i].resize(3);
127  for (unsigned int i = 0; i < 4; ++i)
128  {
129  sg2[i][0] = (1 / sqrt(3)) * lr4[i];
130  sg2[i][1] = (1 / sqrt(3)) * lz4[i];
131  sg2[i][2] = 1.0;
132  }
133  }
134  else if (iord == 3) // 3x3-point Gauss
135  {
136  sg2.resize(9);
137  for (unsigned int i = 0; i < 9; ++i)
138  sg2[i].resize(3);
139  for (unsigned int i = 0; i < 9; ++i)
140  {
141  sg2[i][0] = sqrt(0.6) * lr9[i];
142  sg2[i][1] = sqrt(0.6) * lz9[i];
143  sg2[i][2] = (1.0 / 81.0) * lw9[i];
144  }
145  }
146  else
147  mooseError("Invalid quadrature order = " + Moose::stringify(iord) + " for quad elements");
148  }
149  else if (nen == 3) // triangle
150  {
151  if (iord == 1) // one-point Gauss
152  {
153  sg2.resize(1);
154  sg2[0].resize(4);
155  sg2[0][0] = 1.0 / 3.0;
156  sg2[0][1] = 1.0 / 3.0;
157  sg2[0][2] = 1.0 / 3.0;
158  sg2[0][3] = 0.5;
159  }
160  else if (iord == 2) // three-point Gauss
161  {
162  sg2.resize(3);
163  for (unsigned int i = 0; i < 3; ++i)
164  sg2[i].resize(4);
165  sg2[0][0] = 2.0 / 3.0;
166  sg2[0][1] = 1.0 / 6.0;
167  sg2[0][2] = 1.0 / 6.0;
168  sg2[0][3] = 1.0 / 6.0;
169  sg2[1][0] = 1.0 / 6.0;
170  sg2[1][1] = 2.0 / 3.0;
171  sg2[1][2] = 1.0 / 6.0;
172  sg2[1][3] = 1.0 / 6.0;
173  sg2[2][0] = 1.0 / 6.0;
174  sg2[2][1] = 1.0 / 6.0;
175  sg2[2][2] = 2.0 / 3.0;
176  sg2[2][3] = 1.0 / 6.0;
177  }
178  else if (iord == 3) // four-point Gauss
179  {
180  sg2.resize(4);
181  for (unsigned int i = 0; i < 4; ++i)
182  sg2[i].resize(4);
183  sg2[0][0] = 1.5505102572168219018027159252941e-01;
184  sg2[0][1] = 1.7855872826361642311703513337422e-01;
185  sg2[0][2] = 1.0 - sg2[0][0] - sg2[0][1];
186  sg2[0][3] = 1.5902069087198858469718450103758e-01;
187
188  sg2[1][0] = 6.4494897427831780981972840747059e-01;
189  sg2[1][1] = 7.5031110222608118177475598324603e-02;
190  sg2[1][2] = 1.0 - sg2[1][0] - sg2[1][1];
191  sg2[1][3] = 9.0979309128011415302815498962418e-02;
192
193  sg2[2][0] = 1.5505102572168219018027159252941e-01;
194  sg2[2][1] = 6.6639024601470138670269327409637e-01;
195  sg2[2][2] = 1.0 - sg2[2][0] - sg2[2][1];
196  sg2[2][3] = 1.5902069087198858469718450103758e-01;
197
198  sg2[3][0] = 6.4494897427831780981972840747059e-01;
199  sg2[3][1] = 2.8001991549907407200279599420481e-01;
200  sg2[3][2] = 1.0 - sg2[3][0] - sg2[3][1];
201  sg2[3][3] = 9.0979309128011415302815498962418e-02;
202  }
203  else if (iord == 4) // six-point Guass
204  {
205  const unsigned int n_wts = 2;
206  const Real wts[n_wts] = {1.1169079483900573284750350421656140e-01L,
207  5.4975871827660933819163162450105264e-02L};
208
209  const Real a[n_wts] = {4.4594849091596488631832925388305199e-01L,
210  9.1576213509770743459571463402201508e-02L};
211
212  const Real b[n_wts] = {0., 0.}; // not used
213  const unsigned int permutation_ids[n_wts] = {3, 3};
214
215  std::vector<Point> points;
216  std::vector<Real> weights;
217  dunavant_rule2(wts, a, b, permutation_ids, n_wts, points, weights); // 6 total points
218
219  sg2.resize(6);
220  for (unsigned int i = 0; i < 6; ++i)
221  sg2[i].resize(4);
222  for (unsigned int i = 0; i < 6; ++i)
223  {
224  sg2[i][0] = points[i](0);
225  sg2[i][1] = points[i](1);
226  sg2[i][2] = 1.0 - points[i](0) - points[i](1);
227  sg2[i][3] = weights[i];
228  }
229  }
230  else
231  mooseError("Invalid quadrature order = " + Moose::stringify(iord) + " for triangle elements");
232  }
233  else
234  mooseError("Invalid 2D element type");
235 }
236
237 void
238 wissmannPoints(unsigned int nqp, std::vector<std::vector<Real>> & wss)
239 {
240  if (nqp == 6)
241  {
242  wss.resize(6);
243  for (unsigned int i = 0; i < 6; ++i)
244  wss[i].resize(3);
245  wss[0][0] = 0.0;
246  wss[0][1] = 0.0;
247  wss[0][2] = 1.1428571428571428;
248
249  wss[1][0] = 0.0;
250  wss[1][1] = 9.6609178307929590e-01;
251  wss[1][2] = 4.3956043956043956e-01;
252
253  wss[2][0] = 8.5191465330460049e-01;
254  wss[2][1] = 4.5560372783619284e-01;
255  wss[2][2] = 5.6607220700753210e-01;
256
257  wss[3][0] = -wss[2][0];
258  wss[3][1] = wss[2][1];
259  wss[3][2] = wss[2][2];
260
261  wss[4][0] = 6.3091278897675402e-01;
262  wss[4][1] = -7.3162995157313452e-01;
263  wss[4][2] = 6.4271900178367668e-01;
264
265  wss[5][0] = -wss[4][0];
266  wss[5][1] = wss[4][1];
267  wss[5][2] = wss[4][2];
268  }
269  else
271 }
272
273 void
274 shapeFunc2D(unsigned int nen,
275  std::vector<Real> & ss,
276  std::vector<Point> & xl,
277  std::vector<std::vector<Real>> & shp,
278  Real & xsj,
279  bool natl_flg)
280 {
281  // Get shape functions and derivatives
282  Real s[4] = {-0.5, 0.5, 0.5, -0.5};
283  Real t[4] = {-0.5, -0.5, 0.5, 0.5};
284
285  if (nen == 4) // quad element
286  {
287  Real xs[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
288  Real sx[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
289  for (unsigned int i = 0; i < 4; ++i)
290  {
291  shp[i][2] = (0.5 + s[i] * ss[0]) * (0.5 + t[i] * ss[1]);
292  shp[i][0] = s[i] * (0.5 + t[i] * ss[1]);
293  shp[i][1] = t[i] * (0.5 + s[i] * ss[0]);
294  }
295  for (unsigned int i = 0; i < 2; ++i) // x, y
296  {
297  for (unsigned int j = 0; j < 2; ++j) // xi, eta
298  {
299  xs[i][j] = 0.0;
300  for (unsigned int k = 0; k < nen; ++k)
301  xs[i][j] += xl[k](i) * shp[k][j];
302  }
303  }
304  xsj = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; // det(j)
305  if (natl_flg == false) // get global derivatives
306  {
307  Real temp = 1.0 / xsj;
308  sx[0][0] = xs[1][1] * temp; // inv(j)
309  sx[1][1] = xs[0][0] * temp;
310  sx[0][1] = -xs[0][1] * temp;
311  sx[1][0] = -xs[1][0] * temp;
312  for (unsigned int i = 0; i < nen; ++i)
313  {
314  temp = shp[i][0] * sx[0][0] + shp[i][1] * sx[1][0];
315  shp[i][1] = shp[i][0] * sx[0][1] + shp[i][1] * sx[1][1];
316  shp[i][0] = temp;
317  }
318  }
319  }
320  else if (nen == 3) // triangle element
321  {
322  // x1*(y2 - y3) + x2*(y3 - y1) + x3*(y1 - y2)
323  Point x13 = xl[2] - xl[0];
324  Point x23 = xl[2] - xl[1];
325  Point cross_prod = x13.cross(x23);
326  xsj = cross_prod.norm();
327  Real xsjr = 1.0;
328  if (xsj != 0.0)
329  xsjr = 1.0 / xsj;
330  // xsj *= 0.5; // we do not have this 0.5 here because in stdQuad2D the sum of all weights in
331  // tri is 0.5
332  shp[0][2] = ss[0];
333  shp[1][2] = ss[1];
334  shp[2][2] = ss[2];
335  if (natl_flg == false) // need global drivatives
336  {
337  shp[0][0] = (xl[1](1) - xl[2](1)) * xsjr;
338  shp[0][1] = (xl[2](0) - xl[1](0)) * xsjr;
339  shp[1][0] = (xl[2](1) - xl[0](1)) * xsjr;
340  shp[1][1] = (xl[0](0) - xl[2](0)) * xsjr;
341  shp[2][0] = (xl[0](1) - xl[1](1)) * xsjr;
342  shp[2][1] = (xl[1](0) - xl[0](0)) * xsjr;
343  }
344  else
345  {
346  shp[0][0] = 1.0;
347  shp[0][1] = 0.0;
348  shp[1][0] = 0.0;
349  shp[1][1] = 1.0;
350  shp[2][0] = -1.0;
351  shp[2][1] = -1.0;
352  }
353  }
354  else
355  mooseError("ShapeFunc2D only works for linear quads and tris!");
356 }
357
358 double
359 r8vec_norm(int n, double a[])
360 {
361  // John Burkardt geometry.cpp
362  double v = 0.0;
363  for (int i = 0; i < n; ++i)
364  v = v + a[i] * a[i];
365  v = std::sqrt(v);
366  return v;
367 }
368
369 void
370 r8vec_copy(int n, double a1[], double a2[])
371 {
372  // John Burkardt geometry.cpp
373  for (int i = 0; i < n; ++i)
374  a2[i] = a1[i];
375  return;
376 }
377
378 bool
379 r8vec_eq(int n, double a1[], double a2[])
380 {
381  // John Burkardt geometry.cpp
382  for (int i = 0; i < n; ++i)
383  if (a1[i] != a2[i])
384  return false;
385  return true;
386 }
387
388 double
389 r8vec_dot_product(int n, double a1[], double a2[])
390 {
391  // John Burkardt geometry.cpp
392  double value = 0.0;
393  for (int i = 0; i < n; ++i)
394  value += a1[i] * a2[i];
395  return value;
396 }
397
398 bool
399 line_exp_is_degenerate_nd(int dim_num, double p1[], double p2[])
400 {
401  // John Burkardt geometry.cpp
402  bool value;
403  value = r8vec_eq(dim_num, p1, p2);
404  return value;
405 }
406
407 int
409  double pp[3], double normal[3], double p1[3], double p2[3], double pint[3])
410 {
411 // John Burkardt geometry.cpp
412 // Parameters:
413 //
414 // Input, double PP[3], a point on the plane.
415 //
416 // Input, double NORMAL[3], a normal vector to the plane.
417 //
418 // Input, double P1[3], P2[3], two distinct points on the line.
419 //
420 // Output, double PINT[3], the coordinates of a
421 // common point of the plane and line, when IVAL is 1 or 2.
422 //
423 // Output, integer PLANE_NORMAL_LINE_EXP_INT_3D, the kind of intersection;
424 // 0, the line and plane seem to be parallel and separate;
425 // 1, the line and plane intersect at a single point;
426 // 2, the line and plane seem to be parallel and joined.
427 #define DIM_NUM 3
428
429  double direction[DIM_NUM];
430  int ival;
431  double temp;
432  double temp2;
433  //
434  // Make sure the line is not degenerate.
435  if (line_exp_is_degenerate_nd(DIM_NUM, p1, p2))
436  mooseError("PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error! The line is degenerate.");
437  //
438  // Make sure the plane normal vector is a unit vector.
439  temp = r8vec_norm(DIM_NUM, normal);
440  if (temp == 0.0)
441  mooseError("PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error! The normal vector of the plane is "
442  "degenerate.");
443
444  for (unsigned int i = 0; i < DIM_NUM; ++i)
445  normal[i] = normal[i] / temp;
446  //
447  // Determine the unit direction vector of the line.
448  for (unsigned int i = 0; i < DIM_NUM; ++i)
449  direction[i] = p2[i] - p1[i];
450  temp = r8vec_norm(DIM_NUM, direction);
451
452  for (unsigned int i = 0; i < DIM_NUM; ++i)
453  direction[i] = direction[i] / temp;
454  //
455  // If the normal and direction vectors are orthogonal, then
456  // we have a special case to deal with.
457  if (r8vec_dot_product(DIM_NUM, normal, direction) == 0.0)
458  {
459  temp = 0.0;
460  for (unsigned int i = 0; i < DIM_NUM; ++i)
461  temp = temp + normal[i] * (p1[i] - pp[i]);
462
463  if (temp == 0.0)
464  {
465  ival = 2;
466  r8vec_copy(DIM_NUM, p1, pint);
467  }
468  else
469  {
470  ival = 0;
471  for (unsigned int i = 0; i < DIM_NUM; ++i)
472  pint[i] = 1.0e20; // dummy huge value
473  }
474  return ival;
475  }
476  //
477  // Determine the distance along the direction vector to the intersection point.
478  temp = 0.0;
479  for (unsigned int i = 0; i < DIM_NUM; ++i)
480  temp = temp + normal[i] * (pp[i] - p1[i]);
481  temp2 = 0.0;
482  for (unsigned int i = 0; i < DIM_NUM; ++i)
483  temp2 = temp2 + normal[i] * direction[i];
484
485  ival = 1;
486  for (unsigned int i = 0; i < DIM_NUM; ++i)
487  pint[i] = p1[i] + temp * direction[i] / temp2;
488
489  return ival;
490 #undef DIM_NUM
491 }
492
493 double
495  double coord[], int order_max, int face_num, int node[], int /*node_num*/, int order[])
496 //****************************************************************************80
497 //
498 // Purpose:
499 //
500 // POLYHEDRON_VOLUME_3D computes the volume of a polyhedron in 3D.
501 //
502 // Licensing:
503 //
505 //
506 // Modified:
507 //
508 // 21 August 2003
509 //
510 // Author:
511 //
512 // John Burkardt
513 //
514 // Parameters:
515 //
516 // Input, double COORD[NODE_NUM*3], the 3D coordinates of the vertices.
517 // The vertices may be listed in any order.
518 //
519 // Input, int ORDER_MAX, the maximum number of vertices that make
520 // up a face of the polyhedron.
521 //
522 // Input, int FACE_NUM, the number of faces of the polyhedron.
523 //
524 // Input, int NODE[FACE_NUM*ORDER_MAX]. Face I is defined by
525 // the vertices NODE(I,1) through NODE(I,ORDER(I)). These vertices
526 // are listed in neighboring order.
527 //
528 // Input, int NODE_NUM, the number of points stored in COORD.
529 //
530 // Input, int ORDER[FACE_NUM], the number of vertices making up
531 // each face.
532 //
533 // Output, double POLYHEDRON_VOLUME_3D, the volume of the polyhedron.
534 //
535 {
536 #define DIM_NUM 3
537
538  int face;
539  int n1;
540  int n2;
541  int n3;
542  double term;
543  int v;
544  double volume;
545  double x1;
546  double x2;
547  double x3;
548  double y1;
549  double y2;
550  double y3;
551  double z1;
552  double z2;
553  double z3;
554  //
555  volume = 0.0;
556  //
557  // Triangulate each face.
558  //
559  for (face = 0; face < face_num; face++)
560  {
561  n3 = node[order[face] - 1 + face * order_max];
562  x3 = coord[0 + n3 * 3];
563  y3 = coord[1 + n3 * 3];
564  z3 = coord[2 + n3 * 3];
565
566  for (v = 0; v < order[face] - 2; v++)
567  {
568  n1 = node[v + face * order_max];
569  x1 = coord[0 + n1 * 3];
570  y1 = coord[1 + n1 * 3];
571  z1 = coord[2 + n1 * 3];
572
573  n2 = node[v + 1 + face * order_max];
574  x2 = coord[0 + n2 * 3];
575  y2 = coord[1 + n2 * 3];
576  z2 = coord[2 + n2 * 3];
577
578  term =
579  x1 * y2 * z3 - x1 * y3 * z2 + x2 * y3 * z1 - x2 * y1 * z3 + x3 * y1 * z2 - x3 * y2 * z1;
580
581  volume = volume + term;
582  }
583  }
584
585  volume = volume / 6.0;
586
587  return volume;
588 #undef DIM_NUM
589 }
590
591 void
592 i4vec_zero(int n, int a[])
593 //****************************************************************************80
594 //
595 // Purpose:
596 //
597 // I4VEC_ZERO zeroes an I4VEC.
598 //
599 // Licensing:
600 //
602 //
603 // Modified:
604 //
605 // 01 August 2005
606 //
607 // Author:
608 //
609 // John Burkardt
610 //
611 // Parameters:
612 //
613 // Input, int N, the number of entries in the vector.
614 //
615 // Output, int A[N], a vector of zeroes.
616 //
617 {
618  int i;
619
620  for (i = 0; i < n; i++)
621  {
622  a[i] = 0;
623  }
624  return;
625 }
626
627 void
628 normalizePoint(Point & p)
629 {
630  Real len = p.norm();
631  if (len > tol)
632  p = (1.0 / len) * p;
633  else
634  p.zero();
635 }
636
637 void
639 {
640  Real len = p.norm();
641  if (len > tol)
642  p *= (1.0 / len);
643 }
644
645 } // namespace XFEM
double polyhedron_volume_3d(double coord[], int order_max, int face_num, int node[], int node_num, int order[])
Definition: XFEMFuncs.C:494
void shapeFunc2D(unsigned int nen, std::vector< Real > &ss, std::vector< Point > &xl, std::vector< std::vector< Real >> &shp, Real &xsj, bool natl_flg)
Definition: XFEMFuncs.C:274
void wissmannPoints(unsigned int nqp, std::vector< std::vector< Real >> &wss)
Definition: XFEMFuncs.C:238
bool line_exp_is_degenerate_nd(int dim_num, double p1[], double p2[])
Definition: XFEMFuncs.C:399
double r8vec_norm(int n, double a[])
Definition: XFEMFuncs.C:359
int plane_normal_line_exp_int_3d(double pp[3], double normal[3], double p1[3], double p2[3], double pint[3])
Definition: XFEMFuncs.C:408
void i4vec_zero(int n, int a[])
Definition: XFEMFuncs.C:592
Definition: XFEM.h:24
static const double tol
Definition: XFEMFuncs.h:26
bool r8vec_eq(int n, double a1[], double a2[])
Definition: XFEMFuncs.C:379
void stdQuadr2D(unsigned int nen, unsigned int iord, std::vector< std::vector< Real >> &sg2)
Definition: XFEMFuncs.C:101
double norm()
Definition: EFAPoint.C:80
void dunavant_rule2(const Real *wts, const Real *a, const Real *b, const unsigned int *permutation_ids, unsigned int n_wts, std::vector< Point > &points, std::vector< Real > &weights)
Definition: XFEMFuncs.C:26
double r8vec_dot_product(int n, double a1[], double a2[])
Definition: XFEMFuncs.C:389
void normalizePoint(Point &p)
Definition: XFEMFuncs.C:628
void r8vec_copy(int n, double a1[], double a2[])
Definition: XFEMFuncs.C:370