24 const unsigned int * permutation_ids,
26 std::vector<Point> & points,
27 std::vector<Real> & weights)
33 unsigned int total_pts = 0;
34 for (
unsigned int p = 0; p < n_wts; ++p)
35 total_pts += permutation_ids[p];
38 points.resize(total_pts);
39 weights.resize(total_pts);
42 unsigned int offset = 0;
43 for (
unsigned int p = 0; p < n_wts; ++p)
45 switch (permutation_ids[p])
51 points[offset + 0] =
Point(1.0L / 3.0L, 1.0L / 3.0L);
52 weights[offset + 0] = wts[p];
61 points[offset + 0] =
Point(
a[p],
a[p]);
62 points[offset + 1] =
Point(
a[p], 1.L - 2.L *
a[p]);
63 points[offset + 2] =
Point(1.L - 2.L *
a[p],
a[p]);
65 for (
unsigned int j = 0;
j < 3; ++
j)
66 weights[offset +
j] = wts[p];
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]);
82 for (
unsigned int j = 0;
j < 6; ++
j)
83 weights[offset +
j] = wts[p];
90 mooseError(
"Unknown permutation id: ", permutation_ids[p],
"!");
96 stdQuadr2D(
unsigned int nen,
unsigned int iord, std::vector<std::vector<Real>> & sg2)
101 Real lr4[4] = {-1.0, 1.0, -1.0, 1.0};
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};
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};
120 for (
unsigned int i = 0; i < 4; ++i)
122 for (
unsigned int i = 0; i < 4; ++i)
124 sg2[i][0] = (1 /
sqrt(3)) * lr4[i];
125 sg2[i][1] = (1 /
sqrt(3)) * lz4[i];
132 for (
unsigned int i = 0; i < 9; ++i)
134 for (
unsigned int i = 0; i < 9; ++i)
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];
150 sg2[0][0] = 1.0 / 3.0;
151 sg2[0][1] = 1.0 / 3.0;
152 sg2[0][2] = 1.0 / 3.0;
158 for (
unsigned int i = 0; i < 3; ++i)
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;
176 for (
unsigned int i = 0; i < 4; ++i)
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;
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;
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;
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;
200 const unsigned int n_wts = 2;
201 const Real wts[n_wts] = {1.1169079483900573284750350421656140e-01L,
202 5.4975871827660933819163162450105264e-02L};
204 const Real a[n_wts] = {4.4594849091596488631832925388305199e-01L,
205 9.1576213509770743459571463402201508e-02L};
207 const Real b[n_wts] = {0., 0.};
208 const unsigned int permutation_ids[n_wts] = {3, 3};
210 std::vector<Point> points;
211 std::vector<Real> weights;
215 for (
unsigned int i = 0; i < 6; ++i)
217 for (
unsigned int i = 0; i < 6; ++i)
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];
238 for (
unsigned int i = 0; i < 6; ++i)
242 wss[0][2] = 1.1428571428571428;
245 wss[1][1] = 9.6609178307929590e-01;
246 wss[1][2] = 4.3956043956043956e-01;
248 wss[2][0] = 8.5191465330460049e-01;
249 wss[2][1] = 4.5560372783619284e-01;
250 wss[2][2] = 5.6607220700753210e-01;
252 wss[3][0] = -wss[2][0];
253 wss[3][1] = wss[2][1];
254 wss[3][2] = wss[2][2];
256 wss[4][0] = 6.3091278897675402e-01;
257 wss[4][1] = -7.3162995157313452e-01;
258 wss[4][2] = 6.4271900178367668e-01;
260 wss[5][0] = -wss[4][0];
261 wss[5][1] = wss[4][1];
262 wss[5][2] = wss[4][2];
265 mooseError(
"Unknown Wissmann quadrature type");
270 std::vector<Real> & ss,
271 std::vector<Point> & xl,
272 std::vector<std::vector<Real>> & shp,
277 Real s[4] = {-0.5, 0.5, 0.5, -0.5};
278 Real t[4] = {-0.5, -0.5, 0.5, 0.5};
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)
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]);
290 for (
unsigned int i = 0; i < 2; ++i)
292 for (
unsigned int j = 0;
j < 2; ++
j)
295 for (
unsigned int k = 0;
k < nen; ++
k)
296 xs[i][
j] += xl[
k](i) * shp[
k][
j];
299 xsj = xs[0][0] * xs[1][1] - xs[0][1] * xs[1][0];
300 if (natl_flg ==
false)
302 Real temp = 1.0 / xsj;
303 sx[0][0] = xs[1][1] * temp;
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)
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];
318 Point x13 = xl[2] - xl[0];
319 Point x23 = xl[2] - xl[1];
321 xsj = cross_prod.
norm();
330 if (natl_flg ==
false)
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;
350 mooseError(
"ShapeFunc2D only works for linear quads and tris!");
358 for (
int i = 0; i < n; ++i)
368 for (
int i = 0; i < n; ++i)
377 for (
int i = 0; i < n; ++i)
388 for (
int i = 0; i < n; ++i)
389 value += a1[i] * a2[i];
404 double pp[3],
double normal[3],
double p1[3],
double p2[3],
double pint[3])
424 double direction[DIM_NUM];
431 mooseError(
"PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error! The line is degenerate.");
436 mooseError(
"PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error! The normal vector of the plane is " 439 for (
unsigned int i = 0; i < DIM_NUM; ++i)
440 normal[i] = normal[i] / temp;
443 for (
unsigned int i = 0; i < DIM_NUM; ++i)
444 direction[i] = p2[i] - p1[i];
447 for (
unsigned int i = 0; i < DIM_NUM; ++i)
448 direction[i] = direction[i] / temp;
455 for (
unsigned int i = 0; i < DIM_NUM; ++i)
456 temp = temp + normal[i] * (p1[i] - pp[i]);
466 for (
unsigned int i = 0; i < DIM_NUM; ++i)
474 for (
unsigned int i = 0; i < DIM_NUM; ++i)
475 temp = temp + normal[i] * (pp[i] - p1[i]);
477 for (
unsigned int i = 0; i < DIM_NUM; ++i)
478 temp2 = temp2 + normal[i] * direction[i];
481 for (
unsigned int i = 0; i < DIM_NUM; ++i)
482 pint[i] = p1[i] + temp * direction[i] / temp2;
490 double coord[],
int order_max,
int face_num,
int node[],
int ,
int order[])
553 for (face = 0; face < face_num; face++)
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];
560 for (
v = 0;
v < order[face] - 2;
v++)
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];
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];
573 x1 * y2 * z3 - x1 * y3 * z2 + x2 * y3 * z1 - x2 * y1 * z3 + x3 * y1 * z2 - x3 * y2 * z1;
575 volume = volume + term;
579 volume = volume / 6.0;
613 for (i = 0; i < n; i++)
670 #define PI 3.141592653589793 736 for (i = 0; i < DIM_NUM; i++)
738 v1norm = v1norm +
pow(p1[i] - p2[i], 2);
740 v1norm =
sqrt(v1norm);
749 for (i = 0; i < DIM_NUM; i++)
751 v2norm = v2norm +
pow(p3[i] - p2[i], 2);
753 v2norm =
sqrt(v2norm);
762 for (i = 0; i < DIM_NUM; i++)
764 dot = dot + (p1[i] - p2[i]) * (p3[i] - p2[i]);
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)
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;
795 Real cut_int_frac =
crossProduct2D(cut_start_to_seg_start, seg_dir) / cut_dir_cross_seg_dir;
797 if (cut_int_frac >= 0.0 && cut_int_frac <= cutting_line_fraction)
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)
804 segment_intersection_fraction = int_frac;
814 return (point_a(0) * point_b(1) - point_b(0) * point_a(1));
823 mooseError(
"In XFEMFuncs::pointSegmentDistance(), x0 and x1 should be two different points.");
825 Real s12 = (x2 - x0) * dx / m2;
831 xp = s12 * x1 + (1 - s12) * x2;
841 unsigned int & region)
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;
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)
854 xp = w23 * x1 + w31 * x2 + w12 * x3;
865 if (std::min(distance_12, distance_13) < distance_1)
867 if (distance_12 < distance_13)
893 if (std::min(distance_12, distance_23) < distance_2)
895 if (distance_12 < distance_23)
921 if (std::min(distance_23, distance_31) < distance_3)
923 if (distance_23 < distance_31)
944 mooseError(
"Cannot find closest location in XFEMFuncs::pointTriangleDistance().");
950 const std::vector<Point> & vertices,
953 bool has_intersection =
false;
955 if (vertices.size() != 3)
956 mooseError(
"The number of vertices of cutting element must be 3.");
958 Plane elem_plane(vertices[0], vertices[1], vertices[2]);
959 Point point = vertices[0];
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}};
969 &plane_point[0], &planenormal[0], &edge_point1[0], &edge_point2[0], &cut_point[0]) == 1)
971 Point temp_p(cut_point[0], cut_point[1], cut_point[2]);
975 has_intersection =
true;
979 return has_intersection;
985 Real dotp1 = (p1 - p) * (p2 - p1);
986 Real dotp2 = (p2 - p) * (p2 - p1);
987 return (dotp1 * dotp2 <= 0.0);
995 return len_p1_p / full_len;
1001 unsigned int n_node = vertices.size();
1004 mooseError(
"The number of vertices of cutting element must be 3.");
1006 Plane elem_plane(vertices[0], vertices[1], vertices[2]);
1009 bool inside =
false;
1010 unsigned int counter = 0;
1012 for (
unsigned int i = 0; i < n_node; ++i)
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];
1022 if (middle2p * side_norm <= 0)
1026 if (counter == n_node)
double polyhedron_volume_3d(double coord[], int order_max, int face_num, int node[], int node_num, int order[])
void shapeFunc2D(unsigned int nen, std::vector< Real > &ss, std::vector< Point > &xl, std::vector< std::vector< Real >> &shp, Real &xsj, bool natl_flg)
auto norm() const -> decltype(std::norm(Real()))
Real pointSegmentDistance(const Point &x0, const Point &x1, const Point &x2, Point &xp)
Calculate the signed distance from a point to a line segment.
void mooseError(Args &&... args)
void wissmannPoints(unsigned int nqp, std::vector< std::vector< Real >> &wss)
bool line_exp_is_degenerate_nd(int dim_num, double p1[], double p2[])
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
The following methods are specializations for using the Parallel::packed_range_* routines for a vecto...
double r8vec_norm(int n, double a[])
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
int plane_normal_line_exp_int_3d(double pp[3], double normal[3], double p1[3], double p2[3], double pint[3])
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.
void i4vec_zero(int n, int a[])
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
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...
bool isInsideCutPlane(const std::vector< Point > &vertices, const Point &p)
Check if point p is inside a plane.
TypeVector< typename CompareTypes< Real, T2 >::supertype > cross(const TypeVector< T2 > &v) const
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...
std::string stringify(const T &t)
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 tha...
bool r8vec_eq(int n, double a1[], double a2[])
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::string v
void stdQuadr2D(unsigned int nen, unsigned int iord, std::vector< std::vector< Real >> &sg2)
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)
double r8vec_dot_product(int n, double a1[], double a2[])
bool isInsideEdge(const Point &p1, const Point &p2, const Point &p)
check if point is inside the straight edge p1-p2
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
void normalizePoint(Point &p)
double angle_rad_3d(double p1[3], double p2[3], double p3[3])
virtual Point unit_normal(const Point &p) const override
Real pointTriangleDistance(const Point &x0, const Point &x1, const Point &x2, const Point &x3, Point &xp, unsigned int ®ion)
Calculate the signed distance from a point to a triangle.
static const std::string k
void r8vec_copy(int n, double a1[], double a2[])