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

Classes

struct  CutEdge
 Data structure defining a cut on an element edge. More...
 
struct  CutElemInfo
 Information about a cut element. 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...
 

Typedefs

typedef std::array< std::unordered_map< unsigned int, std::string >, 2 > CachedMaterialProperties
 Convenient typedef for local storage of stateful material properties. 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 intersectSegmentWithCutLine (const Point &segment_point1, const Point &segment_point2, const std::pair< Point, Point > &cutting_line_points, const Real &cutting_line_fraction, Real &segment_intersection_fraction)
 Determine whether a line segment is intersected by a cutting line, and compute the fraction along that line where the intersection occurs. More...
 
Real crossProduct2D (const Point &point_a, const Point &point_b)
 Compute the cross product of two vectors, provided as Point objects, which have nonzero components only in the x,y plane. More...
 
Real pointSegmentDistance (const Point &x0, const Point &x1, const Point &x2, Point &xp)
 Calculate the signed distance from a point to a line segment. More...
 
Real pointTriangleDistance (const Point &x0, const Point &x1, const Point &x2, const Point &x3, Point &xp, unsigned int &region)
 Calculate the signed distance from a point to a triangle. More...
 
bool intersectWithEdge (const Point &p1, const Point &p2, const std::vector< Point > &vertices, Point &pint)
 check if a line intersects with an element defined by vertices calculate the distance from a point to triangle. More...
 
bool isInsideEdge (const Point &p1, const Point &p2, const Point &p)
 check if point is inside the straight edge p1-p2 More...
 
Real getRelativePosition (const Point &p1, const Point &p2, const Point &p)
 Get the relative position of p from p1 respect to the total length of the line segment. More...
 
bool isInsideCutPlane (const std::vector< Point > &vertices, const Point &p)
 Check if point p is inside a plane. More...
 
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
 

Typedef Documentation

◆ CachedMaterialProperties

typedef std::array<std::unordered_map<unsigned int, std::string>, 2> Xfem::CachedMaterialProperties

Convenient typedef for local storage of stateful material properties.

The first (component 0) entry in the CachedMaterialProperties is a map for old material properties. The second (component 1) entry is a map for older material properties. The second entry will be empty if the material storage doesn't have older material properties.

Definition at line 50 of file XFEM.h.

Enumeration Type Documentation

◆ XFEM_CUTPLANE_QUANTITY

Enumerator
ORIGIN_X 
ORIGIN_Y 
ORIGIN_Z 
NORMAL_X 
NORMAL_Y 
NORMAL_Z 

Definition at line 27 of file XFEM.h.

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

◆ XFEM_QRULE

Enumerator
VOLFRAC 
MOMENT_FITTING 
DIRECT 

Definition at line 37 of file XFEM.h.

38 {
39  VOLFRAC,
41  DIRECT
42 };

Function Documentation

◆ angle_rad_3d()

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

Definition at line 691 of file XFEMFuncs.C.

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  for (i = 0; i < DIM_NUM; i++)
737  {
738  v1norm = v1norm + pow(p1[i] - p2[i], 2);
739  }
740  v1norm = sqrt(v1norm);
741 
742  if (v1norm == 0.0)
743  {
744  value = 0.0;
745  return value;
746  }
747 
748  v2norm = 0.0;
749  for (i = 0; i < DIM_NUM; i++)
750  {
751  v2norm = v2norm + pow(p3[i] - p2[i], 2);
752  }
753  v2norm = sqrt(v2norm);
754 
755  if (v2norm == 0.0)
756  {
757  value = 0.0;
758  return value;
759  }
760 
761  dot = 0.0;
762  for (i = 0; i < DIM_NUM; i++)
763  {
764  dot = dot + (p1[i] - p2[i]) * (p3[i] - p2[i]);
765  }
766 
767  value = r8_acos(dot / (v1norm * v2norm));
768 
769  return value;
770 #undef DIM_NUM
771 }
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
double r8_acos(double c)
Definition: XFEMFuncs.C:639

◆ crossProduct2D()

Real Xfem::crossProduct2D ( const Point point_a,
const Point point_b 
)

Compute the cross product of two vectors, provided as Point objects, which have nonzero components only in the x,y plane.

Because the vectors both lie in the x,y plane, the only nonzero component of the cross product vector is in the z direction, and that is returned as a scalar value by this function.

Parameters
point_aFirst vector in cross product
point_bSecond vector in cross product
Returns
z component of cross product vector

Definition at line 812 of file XFEMFuncs.C.

Referenced by intersectSegmentWithCutLine().

813 {
814  return (point_a(0) * point_b(1) - point_b(0) * point_a(1));
815 }

◆ 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 }
void mooseError(Args &&... args)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ getRelativePosition()

Real Xfem::getRelativePosition ( const Point p1,
const Point p2,
const Point p 
)

Get the relative position of p from p1 respect to the total length of the line segment.

Parameters
p1,p2End points of the line segment
pPoint coordinate
Returns
the relative position of p from p1

Definition at line 991 of file XFEMFuncs.C.

Referenced by InterfaceMeshCut3DUserObject::cutElementByGeometry(), and CrackMeshCut3DUserObject::cutElementByGeometry().

992 {
993  Real full_len = (p2 - p1).norm();
994  Real len_p1_p = (p - p1).norm();
995  return len_p1_p / full_len;
996 }
auto norm(const T &a) -> decltype(std::abs(a))
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ 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 }

◆ intersectSegmentWithCutLine()

bool Xfem::intersectSegmentWithCutLine ( const Point segment_point1,
const Point segment_point2,
const std::pair< Point, Point > &  cutting_line_points,
const Real cutting_line_fraction,
Real segment_intersection_fraction 
)

Determine whether a line segment is intersected by a cutting line, and compute the fraction along that line where the intersection occurs.

Parameters
segment_point1Point at one end of the line segment
segment_point1Point at other end of the line segment
cutting_line_pointsPair of points that define the cutting line
cutting_line_fractionFractional distance from the start to end point of the cutting line over which the cutting line is currently active
segment_intersection_fractionFrictional distance along the cut segment from segment_point1 where the intersection occurs
Returns
true if the segment is intersected, false if it is not

Definition at line 774 of file XFEMFuncs.C.

Referenced by XFEMCrackGrowthIncrement2DCut::cutElementByCrackGrowthIncrement(), GeometricCut2DUserObject::cutElementByGeometry(), InterfaceMeshCut2DUserObject::cutElementByGeometry(), MeshCut2DUserObjectBase::cutElementByGeometry(), GeometricCut2DUserObject::cutFragmentByGeometry(), and MeshCut2DUserObjectBase::cutFragmentByGeometry().

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  Real cut_dir_cross_seg_dir = crossProduct2D(cut_dir, seg_dir);
791 
792  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  Real cut_int_frac = crossProduct2D(cut_start_to_seg_start, seg_dir) / cut_dir_cross_seg_dir;
796 
797  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  Real int_frac = crossProduct2D(cut_start_to_seg_start, cut_dir) / cut_dir_cross_seg_dir;
801  if (int_frac >= 0.0 && int_frac <= 1.0)
802  {
803  cut_segment = true;
804  segment_intersection_fraction = int_frac;
805  }
806  }
807  }
808  return cut_segment;
809 }
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Real crossProduct2D(const Point &point_a, const Point &point_b)
Compute the cross product of two vectors, provided as Point objects, which have nonzero components on...
Definition: XFEMFuncs.C:812
static const double tol
Definition: XFEMFuncs.h:21
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ intersectWithEdge()

bool Xfem::intersectWithEdge ( const Point p1,
const Point p2,
const std::vector< Point > &  vertices,
Point pint 
)

check if a line intersects with an element defined by vertices calculate the distance from a point to triangle.

Parameters
p1,p2End points of the line segment
verticesVertices of two-node element.
pintIntersection point
Returns
true if a line intersects with an element

Definition at line 948 of file XFEMFuncs.C.

Referenced by InterfaceMeshCut3DUserObject::cutElementByGeometry(), and CrackMeshCut3DUserObject::cutElementByGeometry().

952 {
953  bool has_intersection = false;
954 
955  if (vertices.size() != 3)
956  mooseError("The number of vertices of cutting element must be 3.");
957 
958  Plane elem_plane(vertices[0], vertices[1], vertices[2]);
959  Point point = vertices[0];
960  Point normal = elem_plane.unit_normal(point);
961 
962  std::array<Real, 3> plane_point = {{point(0), point(1), point(2)}};
963  std::array<Real, 3> planenormal = {{normal(0), normal(1), normal(2)}};
964  std::array<Real, 3> edge_point1 = {{p1(0), p1(1), p1(2)}};
965  std::array<Real, 3> edge_point2 = {{p2(0), p2(1), p2(2)}};
966  std::array<Real, 3> cut_point = {{0.0, 0.0, 0.0}};
967 
969  &plane_point[0], &planenormal[0], &edge_point1[0], &edge_point2[0], &cut_point[0]) == 1)
970  {
971  Point temp_p(cut_point[0], cut_point[1], cut_point[2]);
972  if (isInsideCutPlane(vertices, temp_p) && isInsideEdge(p1, p2, temp_p))
973  {
974  pint = temp_p;
975  has_intersection = true;
976  }
977  }
978 
979  return has_intersection;
980 }
void mooseError(Args &&... args)
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:403
bool isInsideCutPlane(const std::vector< Point > &vertices, const Point &p)
Check if point p is inside a plane.
Definition: XFEMFuncs.C:999
bool isInsideEdge(const Point &p1, const Point &p2, const Point &p)
check if point is inside the straight edge p1-p2
Definition: XFEMFuncs.C:983

◆ isInsideCutPlane()

bool Xfem::isInsideCutPlane ( const std::vector< Point > &  vertices,
const Point p 
)

Check if point p is inside a plane.

Parameters
verticesVertices of the plane
pPoint coordinate
Returns
true if point p is inside a plane

Definition at line 999 of file XFEMFuncs.C.

Referenced by intersectWithEdge().

1000 {
1001  unsigned int n_node = vertices.size();
1002 
1003  if (n_node != 3)
1004  mooseError("The number of vertices of cutting element must be 3.");
1005 
1006  Plane elem_plane(vertices[0], vertices[1], vertices[2]);
1007  Point normal = elem_plane.unit_normal(vertices[0]);
1008 
1009  bool inside = false;
1010  unsigned int counter = 0;
1011 
1012  for (unsigned int i = 0; i < n_node; ++i)
1013  {
1014  unsigned int iplus1 = (i < n_node - 1 ? i + 1 : 0);
1015  Point middle2p = p - 0.5 * (vertices[i] + vertices[iplus1]);
1016  const Point side_tang = vertices[iplus1] - vertices[i];
1017  Point side_norm = side_tang.cross(normal);
1018 
1019  normalizePoint(middle2p);
1020  normalizePoint(side_norm);
1021 
1022  if (middle2p * side_norm <= 0)
1023  counter += 1;
1024  }
1025 
1026  if (counter == n_node)
1027  inside = true;
1028  return inside;
1029 }
void mooseError(Args &&... args)
TypeVector< typename CompareTypes< Real, T2 >::supertype > cross(const TypeVector< T2 > &v) const
void normalizePoint(Point &p)
Definition: XFEMFuncs.C:621

◆ isInsideEdge()

bool Xfem::isInsideEdge ( const Point p1,
const Point p2,
const Point p 
)

check if point is inside the straight edge p1-p2

Parameters
p1,p2End points of the line segment
pPoint coordinate
Returns
true if a point is inside the edge p1-p2

Definition at line 983 of file XFEMFuncs.C.

Referenced by intersectWithEdge().

984 {
985  Real dotp1 = (p1 - p) * (p2 - p1);
986  Real dotp2 = (p2 - p) * (p2 - p1);
987  return (dotp1 * dotp2 <= 0.0);
988 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ 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 }
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
bool r8vec_eq(int n, double a1[], double a2[])
Definition: XFEMFuncs.C:374

◆ normalizePoint() [1/2]

void Xfem::normalizePoint ( Point p)

Definition at line 621 of file XFEMFuncs.C.

Referenced by CircleCutUserObject::CircleCutUserObject(), EFAElement2D::createChild(), EFAElement3D::createChild(), EllipseCutUserObject::EllipseCutUserObject(), XFEMCutElem3D::getCutPlaneNormal(), RectangleCutUserObject::isInsideCutPlane(), isInsideCutPlane(), CrackMeshCut3DUserObject::isInsideCutPlane(), and RectangleCutUserObject::RectangleCutUserObject().

622 {
623  Real len = p.norm();
624  if (len > tol)
625  p = (1.0 / len) * p;
626  else
627  p.zero();
628 }
auto norm() const -> decltype(std::norm(Real()))
const double tol
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ 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
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
double norm()
Definition: EFAPoint.C:84

◆ 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 CrackMeshCut3DUserObject::findIntersection(), GeometricCut3DUserObject::intersectWithEdge(), intersectWithEdge(), and CrackMeshCut3DUserObject::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 }
void mooseError(Args &&... args)
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

◆ pointSegmentDistance()

Real Xfem::pointSegmentDistance ( const Point x0,
const Point x1,
const Point x2,
Point xp 
)

Calculate the signed distance from a point to a line segment.

Positive values are on the side of the line segment's normal (using standard conventions).

Parameters
x1,x2Coordinates of line segment end points
x0Coordinate of the point
xpClosest point coordinate on the line segment
Returns
Distance from a point x0 to a line segment defined by x1-x2

Definition at line 818 of file XFEMFuncs.C.

Referenced by pointTriangleDistance().

819 {
820  Point dx = x2 - x1;
821  Real m2 = dx * dx;
822  if (m2 == 0)
823  mooseError("In XFEMFuncs::pointSegmentDistance(), x0 and x1 should be two different points.");
824  // find parameter coordinate of closest point on segment
825  Real s12 = (x2 - x0) * dx / m2;
826  if (s12 < 0)
827  s12 = 0;
828  else if (s12 > 1)
829  s12 = 1;
830  // and find the distance
831  xp = s12 * x1 + (1 - s12) * x2;
832  return std::sqrt((x0 - xp) * (x0 - xp));
833 }
void mooseError(Args &&... args)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ pointTriangleDistance()

Real Xfem::pointTriangleDistance ( const Point x0,
const Point x1,
const Point x2,
const Point x3,
Point xp,
unsigned int region 
)

Calculate the signed distance from a point to a triangle.

Positive values are on the side of the triangle's normal (using standard conventions).

Parameters
x1,x2,x3Coordinates of triangle vertices
x0Coordinate of the point
xpClosest point coordinate on the triangle
regionThe seven regions where the closest point could be located
Returns
distance from a point x0 to a triangle defined by x1-x2-x3

Definition at line 836 of file XFEMFuncs.C.

Referenced by InterfaceMeshCut3DUserObject::calculateSignedDistance().

842 {
843  Point x13 = x1 - x3, x23 = x2 - x3, x03 = x0 - x3;
844  Real m13 = x13 * x13, m23 = x23 * x23, d = x13 * x23;
845  Real invdet = 1.0 / std::max(m13 * m23 - d * d, 1e-30);
846  Real a = x13 * x03, b = x23 * x03;
847 
848  Real w23 = invdet * (m23 * a - d * b);
849  Real w31 = invdet * (m13 * b - d * a);
850  Real w12 = 1 - w23 - w31;
851  if (w23 >= 0 && w31 >= 0 && w12 >= 0)
852  { // if we're inside the triangle
853  region = 0;
854  xp = w23 * x1 + w31 * x2 + w12 * x3;
855  return std::sqrt((x0 - xp) * (x0 - xp));
856  }
857  else
858  {
859  if (w23 > 0) // this rules out edge 2-3 for us
860  {
861  Point xp1, xp2;
862  Real distance_12 = pointSegmentDistance(x0, x1, x2, xp1);
863  Real distance_13 = pointSegmentDistance(x0, x1, x3, xp2);
864  Real distance_1 = std::sqrt((x0 - x1) * (x0 - x1));
865  if (std::min(distance_12, distance_13) < distance_1)
866  {
867  if (distance_12 < distance_13)
868  {
869  region = 4;
870  xp = xp1;
871  return distance_12;
872  }
873  else
874  {
875  region = 6;
876  xp = xp2;
877  return distance_13;
878  }
879  }
880  else
881  {
882  region = 1;
883  xp = x1;
884  return distance_1;
885  }
886  }
887  else if (w31 > 0) // this rules out edge 1-3
888  {
889  Point xp1, xp2;
890  Real distance_12 = pointSegmentDistance(x0, x1, x2, xp1);
891  Real distance_23 = pointSegmentDistance(x0, x2, x3, xp2);
892  Real distance_2 = std::sqrt((x0 - x2) * (x0 - x2));
893  if (std::min(distance_12, distance_23) < distance_2)
894  {
895  if (distance_12 < distance_23)
896  {
897  region = 4;
898  xp = xp1;
899  return distance_12;
900  }
901  else
902  {
903  region = 5;
904  xp = xp2;
905  return distance_23;
906  }
907  }
908  else
909  {
910  region = 2;
911  xp = x2;
912  return distance_2;
913  }
914  }
915  else // w12 must be >0, ruling out edge 1-2
916  {
917  Point xp1, xp2;
918  Real distance_23 = pointSegmentDistance(x0, x2, x3, xp1);
919  Real distance_31 = pointSegmentDistance(x0, x3, x1, xp2);
920  Real distance_3 = std::sqrt((x0 - x3) * (x0 - x3));
921  if (std::min(distance_23, distance_31) < distance_3)
922  {
923  if (distance_23 < distance_31)
924  {
925  region = 5;
926  xp = xp1;
927  return distance_23;
928  }
929  else
930  {
931  region = 6;
932  xp = xp2;
933  return distance_31;
934  }
935  }
936  else
937  {
938  region = 3;
939  xp = x3;
940  return distance_3;
941  }
942  }
943  }
944  mooseError("Cannot find closest location in XFEMFuncs::pointTriangleDistance().");
945 }
Real pointSegmentDistance(const Point &x0, const Point &x1, const Point &x2, Point &xp)
Calculate the signed distance from a point to a line segment.
Definition: XFEMFuncs.C:818
void mooseError(Args &&... args)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ 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 }
static const std::string v
Definition: NS.h:82

◆ r8_acos()

double Xfem::r8_acos ( double  c)

Definition at line 639 of file XFEMFuncs.C.

Referenced by angle_rad_3d().

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

◆ 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 }
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)

◆ 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 }
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
static const std::string v
Definition: NS.h:82

◆ 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 }
auto norm() const -> decltype(std::norm(Real()))
void mooseError(Args &&... args)
TypeVector< typename CompareTypes< Real, T2 >::supertype > cross(const TypeVector< T2 > &v) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
static const std::string k
Definition: NS.h:124

◆ 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 mooseError(Args &&... args)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
std::string stringify(const T &t)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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 }
void mooseError(Args &&... args)

Variable Documentation

◆ tol

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