www.mooseframework.org
Classes | Enumerations | Functions | Variables
Xfem Namespace Reference

Classes

struct  CutEdge
 Data structure defining a cut on an element edge. More...
 
struct  CutFace
 Data structure defining a cut through a face. More...
 
struct  CutNode
 Data structure defining a cut through a node. More...
 
struct  GeomMarkedElemInfo2D
 Data structure describing geometrically described cut through 2D element. More...
 
struct  GeomMarkedElemInfo3D
 Data structure describing geometrically described cut through 3D element. More...
 

Enumerations

enum  XFEM_CUTPLANE_QUANTITY {
  ORIGIN_X, ORIGIN_Y, ORIGIN_Z, NORMAL_X,
  NORMAL_Y, NORMAL_Z
}
 
enum  XFEM_QRULE { VOLFRAC, MOMENT_FITTING, DIRECT }
 

Functions

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)
 
void stdQuadr2D (unsigned int nen, unsigned int iord, std::vector< std::vector< Real >> &sg2)
 
void wissmannPoints (unsigned int nqp, std::vector< std::vector< Real >> &wss)
 
void shapeFunc2D (unsigned int nen, std::vector< Real > &ss, std::vector< Point > &xl, std::vector< std::vector< Real >> &shp, Real &xsj, bool natl_flg)
 
double r8vec_norm (int n, double a[])
 
void r8vec_copy (int n, double a1[], double a2[])
 
bool r8vec_eq (int n, double a1[], double a2[])
 
double r8vec_dot_product (int n, double a1[], double a2[])
 
bool line_exp_is_degenerate_nd (int dim_num, double p1[], double p2[])
 
int plane_normal_line_exp_int_3d (double pp[3], double normal[3], double p1[3], double p2[3], double pint[3])
 
double polyhedron_volume_3d (double coord[], int order_max, int face_num, int node[], int node_num, int order[])
 
void i4vec_zero (int n, int a[])
 
void normalizePoint (Point &p)
 
void normalizePoint (EFAPoint &p)
 
double r8_acos (double c)
 
double angle_rad_3d (double p1[3], double p2[3], double p3[3])
 
bool operator< (const CutEdge &lhs, const CutEdge &rhs)
 Operator < for two CutEdge Objects Needed to allow the use of std::set<CutEdge> More...
 

Variables

static const double tol = 1.0e-10
 

Enumeration Type Documentation

◆ XFEM_CUTPLANE_QUANTITY

Enumerator
ORIGIN_X 
ORIGIN_Y 
ORIGIN_Z 
NORMAL_X 
NORMAL_Y 
NORMAL_Z 

Definition at line 28 of file XFEM.h.

29 {
30  ORIGIN_X,
31  ORIGIN_Y,
32  ORIGIN_Z,
33  NORMAL_X,
34  NORMAL_Y,
35  NORMAL_Z
36 };

◆ XFEM_QRULE

Enumerator
VOLFRAC 
MOMENT_FITTING 
DIRECT 

Definition at line 38 of file XFEM.h.

39 {
40  VOLFRAC,
42  DIRECT
43 };

Function Documentation

◆ angle_rad_3d()

double Xfem::angle_rad_3d ( double  p1[3],
double  p2[3],
double  p3[3] 
)

Definition at line 692 of file XFEMFuncs.C.

Referenced by MeshCut3DUserObject::findBoundaryNodes().

728 {
729 #define DIM_NUM 3
730 
731  double dot;
732  int i;
733  double v1norm;
734  double v2norm;
735  double value;
736 
737  v1norm = 0.0;
738  for (i = 0; i < DIM_NUM; i++)
739  {
740  v1norm = v1norm + pow(p1[i] - p2[i], 2);
741  }
742  v1norm = sqrt(v1norm);
743 
744  if (v1norm == 0.0)
745  {
746  value = 0.0;
747  return value;
748  }
749 
750  v2norm = 0.0;
751  for (i = 0; i < DIM_NUM; i++)
752  {
753  v2norm = v2norm + pow(p3[i] - p2[i], 2);
754  }
755  v2norm = sqrt(v2norm);
756 
757  if (v2norm == 0.0)
758  {
759  value = 0.0;
760  return value;
761  }
762 
763  dot = 0.0;
764  for (i = 0; i < DIM_NUM; i++)
765  {
766  dot = dot + (p1[i] - p2[i]) * (p3[i] - p2[i]);
767  }
768 
769  value = r8_acos(dot / (v1norm * v2norm));
770 
771  return value;
772 #undef DIM_NUM
773 }
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
double r8_acos(double c)
Definition: XFEMFuncs.C:639

◆ dunavant_rule2()

void Xfem::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 at line 21 of file XFEMFuncs.C.

Referenced by stdQuadr2D().

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  for (unsigned int p = 0; p < n_wts; ++p)
35  total_pts += permutation_ids[p];
36 
37  // Resize point and weight vectors appropriately.
38  points.resize(total_pts);
39  weights.resize(total_pts);
40 
41  // Always insert into the points & weights vector relative to the offset
42  unsigned int offset = 0;
43  for (unsigned int p = 0; p < n_wts; ++p)
44  {
45  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  points[offset + 0] = Point(1.0L / 3.0L, 1.0L / 3.0L);
52  weights[offset + 0] = wts[p];
53 
54  offset += 1;
55  break;
56  }
57 
58  case 3:
59  {
60  // For this type of rule, don't need to look in the b array.
61  points[offset + 0] = Point(a[p], a[p]); // (a,a)
62  points[offset + 1] = Point(a[p], 1.L - 2.L * a[p]); // (a,1-2a)
63  points[offset + 2] = Point(1.L - 2.L * a[p], a[p]); // (1-2a,a)
64 
65  for (unsigned int j = 0; j < 3; ++j)
66  weights[offset + j] = wts[p];
67 
68  offset += 3;
69  break;
70  }
71 
72  case 6:
73  {
74  // This type of point uses all 3 arrays...
75  points[offset + 0] = Point(a[p], b[p]);
76  points[offset + 1] = Point(b[p], a[p]);
77  points[offset + 2] = Point(a[p], 1.L - a[p] - b[p]);
78  points[offset + 3] = Point(1.L - a[p] - b[p], a[p]);
79  points[offset + 4] = Point(b[p], 1.L - a[p] - b[p]);
80  points[offset + 5] = Point(1.L - a[p] - b[p], b[p]);
81 
82  for (unsigned int j = 0; j < 6; ++j)
83  weights[offset + j] = wts[p];
84 
85  offset += 6;
86  break;
87  }
88 
89  default:
90  mooseError("Unknown permutation id: ", permutation_ids[p], "!");
91  }
92  }
93 }

◆ i4vec_zero()

void Xfem::i4vec_zero ( int  n,
int  a[] 
)

Definition at line 586 of file XFEMFuncs.C.

Referenced by XFEMCutElem3D::computePhysicalVolumeFraction().

610 {
611  int i;
612 
613  for (i = 0; i < n; i++)
614  {
615  a[i] = 0;
616  }
617  return;
618 }

◆ line_exp_is_degenerate_nd()

bool Xfem::line_exp_is_degenerate_nd ( int  dim_num,
double  p1[],
double  p2[] 
)

Definition at line 394 of file XFEMFuncs.C.

Referenced by plane_normal_line_exp_int_3d().

395 {
396  // John Burkardt geometry.cpp
397  bool value;
398  value = r8vec_eq(dim_num, p1, p2);
399  return value;
400 }
bool r8vec_eq(int n, double a1[], double a2[])
Definition: XFEMFuncs.C:374

◆ normalizePoint() [1/2]

void Xfem::normalizePoint ( Point &  p)

◆ normalizePoint() [2/2]

void Xfem::normalizePoint ( EFAPoint p)

Definition at line 631 of file XFEMFuncs.C.

632 {
633  Real len = p.norm();
634  if (len > tol)
635  p *= (1.0 / len);
636 }
const double tol
Definition: Setup.h:19
double norm()
Definition: EFAPoint.C:82

◆ operator<()

bool Xfem::operator< ( const CutEdge lhs,
const CutEdge rhs 
)
inline

Operator < for two CutEdge Objects Needed to allow the use of std::set<CutEdge>

Parameters
lhsCutEdge object on the left side of the comparison
rhsCutEdge object on the right side of the comparison
Returns
bool true if lhs < rhs

Definition at line 47 of file GeometricCutUserObject.h.

48 {
49  if (lhs._id1 != rhs._id1)
50  return lhs._id1 < rhs._id1;
51  else
52  return lhs._id2 < rhs._id2;
53 }

◆ plane_normal_line_exp_int_3d()

int Xfem::plane_normal_line_exp_int_3d ( double  pp[3],
double  normal[3],
double  p1[3],
double  p2[3],
double  pint[3] 
)

Definition at line 403 of file XFEMFuncs.C.

Referenced by MeshCut3DUserObject::findIntersection(), GeometricCut3DUserObject::intersectWithEdge(), and MeshCut3DUserObject::intersectWithEdge().

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  if (line_exp_is_degenerate_nd(DIM_NUM, p1, p2))
431  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  temp = r8vec_norm(DIM_NUM, normal);
435  if (temp == 0.0)
436  mooseError("PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error! The normal vector of the plane is "
437  "degenerate.");
438 
439  for (unsigned int i = 0; i < DIM_NUM; ++i)
440  normal[i] = normal[i] / temp;
441  //
442  // Determine the unit direction vector of the line.
443  for (unsigned int i = 0; i < DIM_NUM; ++i)
444  direction[i] = p2[i] - p1[i];
445  temp = r8vec_norm(DIM_NUM, direction);
446 
447  for (unsigned int i = 0; i < DIM_NUM; ++i)
448  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  if (r8vec_dot_product(DIM_NUM, normal, direction) == 0.0)
453  {
454  temp = 0.0;
455  for (unsigned int i = 0; i < DIM_NUM; ++i)
456  temp = temp + normal[i] * (p1[i] - pp[i]);
457 
458  if (temp == 0.0)
459  {
460  ival = 2;
461  r8vec_copy(DIM_NUM, p1, pint);
462  }
463  else
464  {
465  ival = 0;
466  for (unsigned int i = 0; i < DIM_NUM; ++i)
467  pint[i] = 1.0e20; // dummy huge value
468  }
469  return ival;
470  }
471  //
472  // Determine the distance along the direction vector to the intersection point.
473  temp = 0.0;
474  for (unsigned int i = 0; i < DIM_NUM; ++i)
475  temp = temp + normal[i] * (pp[i] - p1[i]);
476  temp2 = 0.0;
477  for (unsigned int i = 0; i < DIM_NUM; ++i)
478  temp2 = temp2 + normal[i] * direction[i];
479 
480  ival = 1;
481  for (unsigned int i = 0; i < DIM_NUM; ++i)
482  pint[i] = p1[i] + temp * direction[i] / temp2;
483 
484  return ival;
485 #undef DIM_NUM
486 }
bool line_exp_is_degenerate_nd(int dim_num, double p1[], double p2[])
Definition: XFEMFuncs.C:394
double r8vec_norm(int n, double a[])
Definition: XFEMFuncs.C:354
double r8vec_dot_product(int n, double a1[], double a2[])
Definition: XFEMFuncs.C:384
void r8vec_copy(int n, double a1[], double a2[])
Definition: XFEMFuncs.C:365

◆ polyhedron_volume_3d()

double Xfem::polyhedron_volume_3d ( double  coord[],
int  order_max,
int  face_num,
int  node[],
int  node_num,
int  order[] 
)

Definition at line 489 of file XFEMFuncs.C.

Referenced by XFEMCutElem3D::computePhysicalVolumeFraction().

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  for (face = 0; face < face_num; face++)
554  {
555  n3 = node[order[face] - 1 + face * order_max];
556  x3 = coord[0 + n3 * 3];
557  y3 = coord[1 + n3 * 3];
558  z3 = coord[2 + n3 * 3];
559 
560  for (v = 0; v < order[face] - 2; v++)
561  {
562  n1 = node[v + face * order_max];
563  x1 = coord[0 + n1 * 3];
564  y1 = coord[1 + n1 * 3];
565  z1 = coord[2 + n1 * 3];
566 
567  n2 = node[v + 1 + face * order_max];
568  x2 = coord[0 + n2 * 3];
569  y2 = coord[1 + n2 * 3];
570  z2 = coord[2 + n2 * 3];
571 
572  term =
573  x1 * y2 * z3 - x1 * y3 * z2 + x2 * y3 * z1 - x2 * y1 * z3 + x3 * y1 * z2 - x3 * y2 * z1;
574 
575  volume = volume + term;
576  }
577  }
578 
579  volume = volume / 6.0;
580 
581  return volume;
582 #undef DIM_NUM
583 }

◆ r8_acos()

double Xfem::r8_acos ( double  c)

Definition at line 639 of file XFEMFuncs.C.

Referenced by angle_rad_3d().

670 {
671 #define PI 3.141592653589793
672 
673  double value;
674 
675  if (c <= -1.0)
676  {
677  value = PI;
678  }
679  else if (1.0 <= c)
680  {
681  value = 0.0;
682  }
683  else
684  {
685  value = acos(c);
686  }
687  return value;
688 #undef PI
689 }

◆ r8vec_copy()

void Xfem::r8vec_copy ( int  n,
double  a1[],
double  a2[] 
)

Definition at line 365 of file XFEMFuncs.C.

Referenced by plane_normal_line_exp_int_3d().

366 {
367  // John Burkardt geometry.cpp
368  for (int i = 0; i < n; ++i)
369  a2[i] = a1[i];
370  return;
371 }

◆ r8vec_dot_product()

double Xfem::r8vec_dot_product ( int  n,
double  a1[],
double  a2[] 
)

Definition at line 384 of file XFEMFuncs.C.

Referenced by plane_normal_line_exp_int_3d().

385 {
386  // John Burkardt geometry.cpp
387  double value = 0.0;
388  for (int i = 0; i < n; ++i)
389  value += a1[i] * a2[i];
390  return value;
391 }

◆ r8vec_eq()

bool Xfem::r8vec_eq ( int  n,
double  a1[],
double  a2[] 
)

Definition at line 374 of file XFEMFuncs.C.

Referenced by line_exp_is_degenerate_nd().

375 {
376  // John Burkardt geometry.cpp
377  for (int i = 0; i < n; ++i)
378  if (a1[i] != a2[i])
379  return false;
380  return true;
381 }

◆ r8vec_norm()

double Xfem::r8vec_norm ( int  n,
double  a[] 
)

Definition at line 354 of file XFEMFuncs.C.

Referenced by plane_normal_line_exp_int_3d().

355 {
356  // John Burkardt geometry.cpp
357  double v = 0.0;
358  for (int i = 0; i < n; ++i)
359  v = v + a[i] * a[i];
360  v = std::sqrt(v);
361  return v;
362 }

◆ shapeFunc2D()

void Xfem::shapeFunc2D ( unsigned int  nen,
std::vector< Real > &  ss,
std::vector< Point > &  xl,
std::vector< std::vector< Real >> &  shp,
Real &  xsj,
bool  natl_flg 
)

Definition at line 269 of file XFEMFuncs.C.

Referenced by XFEMCutElem2D::getPhysicalQuadraturePoints(), XFEM::getXFEMqRuleOnSurface(), and XFEMCutElem2D::solveMomentFitting().

275 {
276  // Get shape functions and derivatives
277  Real s[4] = {-0.5, 0.5, 0.5, -0.5};
278  Real t[4] = {-0.5, -0.5, 0.5, 0.5};
279 
280  if (nen == 4) // quad element
281  {
282  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  for (unsigned int i = 0; i < 4; ++i)
285  {
286  shp[i][2] = (0.5 + s[i] * ss[0]) * (0.5 + t[i] * ss[1]);
287  shp[i][0] = s[i] * (0.5 + t[i] * ss[1]);
288  shp[i][1] = t[i] * (0.5 + s[i] * ss[0]);
289  }
290  for (unsigned int i = 0; i < 2; ++i) // x, y
291  {
292  for (unsigned int j = 0; j < 2; ++j) // xi, eta
293  {
294  xs[i][j] = 0.0;
295  for (unsigned int k = 0; k < nen; ++k)
296  xs[i][j] += xl[k](i) * shp[k][j];
297  }
298  }
299  xsj = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0]; // det(j)
300  if (natl_flg == false) // get global derivatives
301  {
302  Real temp = 1.0 / xsj;
303  sx[0][0] = xs[1][1] * temp; // inv(j)
304  sx[1][1] = xs[0][0] * temp;
305  sx[0][1] = -xs[0][1] * temp;
306  sx[1][0] = -xs[1][0] * temp;
307  for (unsigned int i = 0; i < nen; ++i)
308  {
309  temp = shp[i][0] * sx[0][0] + shp[i][1] * sx[1][0];
310  shp[i][1] = shp[i][0] * sx[0][1] + shp[i][1] * sx[1][1];
311  shp[i][0] = temp;
312  }
313  }
314  }
315  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  Point cross_prod = x13.cross(x23);
321  xsj = cross_prod.norm();
322  Real xsjr = 1.0;
323  if (xsj != 0.0)
324  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  shp[0][2] = ss[0];
328  shp[1][2] = ss[1];
329  shp[2][2] = ss[2];
330  if (natl_flg == false) // need global drivatives
331  {
332  shp[0][0] = (xl[1](1) - xl[2](1)) * xsjr;
333  shp[0][1] = (xl[2](0) - xl[1](0)) * xsjr;
334  shp[1][0] = (xl[2](1) - xl[0](1)) * xsjr;
335  shp[1][1] = (xl[0](0) - xl[2](0)) * xsjr;
336  shp[2][0] = (xl[0](1) - xl[1](1)) * xsjr;
337  shp[2][1] = (xl[1](0) - xl[0](0)) * xsjr;
338  }
339  else
340  {
341  shp[0][0] = 1.0;
342  shp[0][1] = 0.0;
343  shp[1][0] = 0.0;
344  shp[1][1] = 1.0;
345  shp[2][0] = -1.0;
346  shp[2][1] = -1.0;
347  }
348  }
349  else
350  mooseError("ShapeFunc2D only works for linear quads and tris!");
351 }

◆ stdQuadr2D()

void Xfem::stdQuadr2D ( unsigned int  nen,
unsigned int  iord,
std::vector< std::vector< Real >> &  sg2 
)

Definition at line 96 of file XFEMFuncs.C.

Referenced by XFEMCutElem2D::getPhysicalQuadraturePoints(), and XFEM::getXFEMqRuleOnSurface().

97 {
98  // Purpose: get Guass integration points for 2D quad and tri elems
99  // N.B. only works for n_qp <= 6
100 
101  Real lr4[4] = {-1.0, 1.0, -1.0, 1.0}; // libmesh order
102  Real lz4[4] = {-1.0, -1.0, 1.0, 1.0};
103  Real lr9[9] = {-1.0, 0.0, 1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0}; // libmesh order
104  Real lz9[9] = {-1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0};
105  Real lw9[9] = {25.0, 40.0, 25.0, 40.0, 64.0, 40.0, 25.0, 40.0, 25.0};
106 
107  if (nen == 4) // 2d quad element
108  {
109  if (iord == 1) // 1-point Gauss
110  {
111  sg2.resize(1);
112  sg2[0].resize(3);
113  sg2[0][0] = 0.0;
114  sg2[0][1] = 0.0;
115  sg2[0][2] = 4.0;
116  }
117  else if (iord == 2) // 2x2-point Gauss
118  {
119  sg2.resize(4);
120  for (unsigned int i = 0; i < 4; ++i)
121  sg2[i].resize(3);
122  for (unsigned int i = 0; i < 4; ++i)
123  {
124  sg2[i][0] = (1 / sqrt(3)) * lr4[i];
125  sg2[i][1] = (1 / sqrt(3)) * lz4[i];
126  sg2[i][2] = 1.0;
127  }
128  }
129  else if (iord == 3) // 3x3-point Gauss
130  {
131  sg2.resize(9);
132  for (unsigned int i = 0; i < 9; ++i)
133  sg2[i].resize(3);
134  for (unsigned int i = 0; i < 9; ++i)
135  {
136  sg2[i][0] = sqrt(0.6) * lr9[i];
137  sg2[i][1] = sqrt(0.6) * lz9[i];
138  sg2[i][2] = (1.0 / 81.0) * lw9[i];
139  }
140  }
141  else
142  mooseError("Invalid quadrature order = " + Moose::stringify(iord) + " for quad elements");
143  }
144  else if (nen == 3) // triangle
145  {
146  if (iord == 1) // one-point Gauss
147  {
148  sg2.resize(1);
149  sg2[0].resize(4);
150  sg2[0][0] = 1.0 / 3.0;
151  sg2[0][1] = 1.0 / 3.0;
152  sg2[0][2] = 1.0 / 3.0;
153  sg2[0][3] = 0.5;
154  }
155  else if (iord == 2) // three-point Gauss
156  {
157  sg2.resize(3);
158  for (unsigned int i = 0; i < 3; ++i)
159  sg2[i].resize(4);
160  sg2[0][0] = 2.0 / 3.0;
161  sg2[0][1] = 1.0 / 6.0;
162  sg2[0][2] = 1.0 / 6.0;
163  sg2[0][3] = 1.0 / 6.0;
164  sg2[1][0] = 1.0 / 6.0;
165  sg2[1][1] = 2.0 / 3.0;
166  sg2[1][2] = 1.0 / 6.0;
167  sg2[1][3] = 1.0 / 6.0;
168  sg2[2][0] = 1.0 / 6.0;
169  sg2[2][1] = 1.0 / 6.0;
170  sg2[2][2] = 2.0 / 3.0;
171  sg2[2][3] = 1.0 / 6.0;
172  }
173  else if (iord == 3) // four-point Gauss
174  {
175  sg2.resize(4);
176  for (unsigned int i = 0; i < 4; ++i)
177  sg2[i].resize(4);
178  sg2[0][0] = 1.5505102572168219018027159252941e-01;
179  sg2[0][1] = 1.7855872826361642311703513337422e-01;
180  sg2[0][2] = 1.0 - sg2[0][0] - sg2[0][1];
181  sg2[0][3] = 1.5902069087198858469718450103758e-01;
182 
183  sg2[1][0] = 6.4494897427831780981972840747059e-01;
184  sg2[1][1] = 7.5031110222608118177475598324603e-02;
185  sg2[1][2] = 1.0 - sg2[1][0] - sg2[1][1];
186  sg2[1][3] = 9.0979309128011415302815498962418e-02;
187 
188  sg2[2][0] = 1.5505102572168219018027159252941e-01;
189  sg2[2][1] = 6.6639024601470138670269327409637e-01;
190  sg2[2][2] = 1.0 - sg2[2][0] - sg2[2][1];
191  sg2[2][3] = 1.5902069087198858469718450103758e-01;
192 
193  sg2[3][0] = 6.4494897427831780981972840747059e-01;
194  sg2[3][1] = 2.8001991549907407200279599420481e-01;
195  sg2[3][2] = 1.0 - sg2[3][0] - sg2[3][1];
196  sg2[3][3] = 9.0979309128011415302815498962418e-02;
197  }
198  else if (iord == 4) // six-point Guass
199  {
200  const unsigned int n_wts = 2;
201  const Real wts[n_wts] = {1.1169079483900573284750350421656140e-01L,
202  5.4975871827660933819163162450105264e-02L};
203 
204  const Real a[n_wts] = {4.4594849091596488631832925388305199e-01L,
205  9.1576213509770743459571463402201508e-02L};
206 
207  const Real b[n_wts] = {0., 0.}; // not used
208  const unsigned int permutation_ids[n_wts] = {3, 3};
209 
210  std::vector<Point> points;
211  std::vector<Real> weights;
212  dunavant_rule2(wts, a, b, permutation_ids, n_wts, points, weights); // 6 total points
213 
214  sg2.resize(6);
215  for (unsigned int i = 0; i < 6; ++i)
216  sg2[i].resize(4);
217  for (unsigned int i = 0; i < 6; ++i)
218  {
219  sg2[i][0] = points[i](0);
220  sg2[i][1] = points[i](1);
221  sg2[i][2] = 1.0 - points[i](0) - points[i](1);
222  sg2[i][3] = weights[i];
223  }
224  }
225  else
226  mooseError("Invalid quadrature order = " + Moose::stringify(iord) + " for triangle elements");
227  }
228  else
229  mooseError("Invalid 2D element type");
230 }
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:21

◆ wissmannPoints()

void Xfem::wissmannPoints ( unsigned int  nqp,
std::vector< std::vector< Real >> &  wss 
)

Definition at line 233 of file XFEMFuncs.C.

234 {
235  if (nqp == 6)
236  {
237  wss.resize(6);
238  for (unsigned int i = 0; i < 6; ++i)
239  wss[i].resize(3);
240  wss[0][0] = 0.0;
241  wss[0][1] = 0.0;
242  wss[0][2] = 1.1428571428571428;
243 
244  wss[1][0] = 0.0;
245  wss[1][1] = 9.6609178307929590e-01;
246  wss[1][2] = 4.3956043956043956e-01;
247 
248  wss[2][0] = 8.5191465330460049e-01;
249  wss[2][1] = 4.5560372783619284e-01;
250  wss[2][2] = 5.6607220700753210e-01;
251 
252  wss[3][0] = -wss[2][0];
253  wss[3][1] = wss[2][1];
254  wss[3][2] = wss[2][2];
255 
256  wss[4][0] = 6.3091278897675402e-01;
257  wss[4][1] = -7.3162995157313452e-01;
258  wss[4][2] = 6.4271900178367668e-01;
259 
260  wss[5][0] = -wss[4][0];
261  wss[5][1] = wss[4][1];
262  wss[5][2] = wss[4][2];
263  }
264  else
265  mooseError("Unknown Wissmann quadrature type");
266 }

Variable Documentation

◆ tol

const double Xfem::tol = 1.0e-10
static