19 #include "libmesh/mesh.h" 20 #include "libmesh/elem.h" 21 #include "libmesh/node.h" 25 unsigned int n_qpoints,
27 :
XFEMCutElem(elem, n_qpoints, n_sides), _efa_elem3d(CEMelem, true)
38 Point node_coor(0.0, 0.0, 0.0);
39 std::vector<EFANode *> master_nodes;
40 std::vector<Point> master_points;
41 std::vector<double> master_weights;
44 for (
unsigned int i = 0; i < master_nodes.size(); ++i)
48 Node * node =
_nodes[master_nodes[i]->id()];
50 node = displaced_mesh->node_ptr(node->id());
51 Point node_p((*node)(0), (*node)(1), (*node)(2));
52 master_points.push_back(node_p);
57 for (
unsigned int i = 0; i < master_nodes.size(); ++i)
58 node_coor += master_weights[i] * master_points[i];
68 std::vector<std::vector<unsigned int>> frag_face_indices;
69 std::vector<EFANode *> frag_nodes;
71 int face_num = frag_face_indices.size();
72 int node_num = frag_nodes.size();
75 int * order =
new int[face_num];
76 for (
int i = 0; i < face_num; ++i)
78 if (frag_face_indices[i].size() > (
unsigned int)order_max)
79 order_max = frag_face_indices[i].size();
80 order[i] = frag_face_indices[i].size();
83 double * coord =
new double[3 * node_num];
84 for (
unsigned int i = 0; i < frag_nodes.size(); ++i)
87 coord[3 * i + 0] = p(0);
88 coord[3 * i + 1] = p(1);
89 coord[3 * i + 2] = p(2);
92 int * node =
new int[face_num * order_max];
94 for (
int i = 0; i < face_num; ++i)
95 for (
unsigned int j = 0;
j < frag_face_indices[i].size(); ++
j)
96 node[order_max * i +
j] = frag_face_indices[i][
j];
110 Real frag_surf = 0.0;
112 std::vector<std::vector<unsigned int>> frag_face_indices;
113 std::vector<EFANode *> frag_nodes;
115 int face_num = frag_face_indices.size();
118 bool contains_all =
true;
121 for (
int i = 0; i < face_num; ++i)
124 for (
unsigned int j = 0;
j < frag_face_indices[i].size(); ++
j)
126 EFANode * efa_node = frag_nodes[frag_face_indices[i][
j]];
129 contains_all =
false;
136 for (
unsigned int j = 0;
j < frag_face_indices[i].size(); ++
j)
138 unsigned int m = ((
j + 1) == frag_face_indices[i].size()) ? 0 :
j + 1;
142 frag_surf += 0.5 * (edge_p1(0) - edge_p2(0)) * (edge_p1(1) + edge_p2(1));
161 Point orig(0.0, 0.0, 0.0);
162 std::vector<std::vector<EFANode *>> cut_plane_nodes;
168 std::vector<EFANode *> node_line;
169 for (
unsigned int j = 0;
j < face->
numNodes(); ++
j)
170 node_line.push_back(face->
getNode(
j));
171 cut_plane_nodes.push_back(node_line);
174 if (cut_plane_nodes.size() == 0)
175 mooseError(
"no cut plane found in this element");
176 if (plane_id < cut_plane_nodes.size())
178 std::vector<Point> cut_plane_points;
179 for (
unsigned int i = 0; i < cut_plane_nodes[plane_id].size(); ++i)
180 cut_plane_points.push_back(
getNodeCoordinates(cut_plane_nodes[plane_id][i], displaced_mesh));
182 for (
unsigned int i = 0; i < cut_plane_points.size(); ++i)
183 orig += cut_plane_points[i];
184 orig /= cut_plane_points.size();
192 Point normal(0.0, 0.0, 0.0);
193 std::vector<std::vector<EFANode *>> cut_plane_nodes;
199 std::vector<EFANode *> node_line;
200 for (
unsigned int j = 0;
j < face->
numNodes(); ++
j)
201 node_line.push_back(face->
getNode(
j));
202 cut_plane_nodes.push_back(node_line);
205 if (cut_plane_nodes.size() == 0)
206 mooseError(
"no cut plane found in this element");
207 if (plane_id < cut_plane_nodes.size())
209 std::vector<Point> cut_plane_points;
210 for (
unsigned int i = 0; i < cut_plane_nodes[plane_id].size(); ++i)
211 cut_plane_points.push_back(
getNodeCoordinates(cut_plane_nodes[plane_id][i], displaced_mesh));
213 Point
center(0.0, 0.0, 0.0);
214 for (
unsigned int i = 0; i < cut_plane_points.size(); ++i)
215 center += cut_plane_points[i];
216 center /= cut_plane_points.size();
218 for (
unsigned int i = 0; i < cut_plane_points.size(); ++i)
220 unsigned int iplus1 = i < cut_plane_points.size() - 1 ? i + 1 : 0;
221 Point ray1 = cut_plane_points[i] -
center;
222 Point ray2 = cut_plane_points[iplus1] -
center;
223 normal += ray1.cross(ray2);
225 normal /= cut_plane_points.size();
237 mooseError(
"getCrackTipOriginAndDirection not yet implemented for XFEMCutElem3D");
245 mooseError(
"getFragmentFaces not yet implemented for XFEMCutElem3D");
257 unsigned int counter = 0;
267 std::vector<Point> & intersectionPoints,
268 MeshBase * displaced_mesh)
const 270 intersectionPoints.clear();
271 std::vector<std::vector<EFANode *>> cut_plane_nodes;
277 std::vector<EFANode *> node_line;
278 for (
unsigned int j = 0;
j < face->
numNodes(); ++
j)
279 node_line.push_back(face->
getNode(
j));
280 cut_plane_nodes.push_back(node_line);
283 if (cut_plane_nodes.size() == 0)
284 mooseError(
"No cut plane found in this element");
286 if (plane_id < cut_plane_nodes.size())
288 intersectionPoints.resize(cut_plane_nodes[plane_id].size());
289 for (
unsigned int i = 0; i < cut_plane_nodes[plane_id].size(); ++i)
290 intersectionPoints[i] =
getNodeCoordinates(cut_plane_nodes[plane_id][i], displaced_mesh);
EFAFragment3D * getFragment(unsigned int frag_id) const
double polyhedron_volume_3d(double coord[], int order_max, int face_num, int node[], int node_num, int order[])
EFAFace * getFace(unsigned int face_id) const
std::vector< Real > _elem_side_area
bool containsNode(const EFANode *node) const
virtual const EFAElement * getEFAElement() const
void mooseError(Args &&... args)
virtual void getIntersectionInfo(unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, MeshBase *displaced_mesh=nullptr) const
virtual void computeMomentFittingWeights()
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=nullptr) const
virtual void computePhysicalFaceAreaFraction(unsigned int side)
Computes the surface area fraction of the element side.
EFANode * getNode(unsigned int node_id) const
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
virtual void getCrackTipOriginAndDirection(unsigned tip_id, Point &origin, Point &direction) const
bool isFaceInterior(unsigned int face_id) const
void getNodeInfo(std::vector< std::vector< unsigned int >> &face_node_indices, std::vector< EFANode *> &nodes) const
void i4vec_zero(int n, int a[])
EFAFace * getFace(unsigned int face_id) const
unsigned int numNodes() const
virtual Point getCutPlaneOrigin(unsigned int plane_id, MeshBase *displaced_mesh=nullptr) const
virtual void getMasterInfo(EFANode *node, std::vector< EFANode *> &master_nodes, std::vector< double > &master_weights) const
std::vector< Real > _new_weights
quadrature weights from volume fraction and moment fitting
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=nullptr) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void computePhysicalVolumeFraction()
Computes the volume fraction of the element fragment.
std::vector< Node * > _nodes
virtual unsigned int numCutPlanes() const
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
void normalizePoint(Point &p)
virtual void getFragmentFaces(std::vector< std::vector< Point >> &frag_faces, MeshBase *displaced_mesh=nullptr) const
std::vector< Real > _physical_areafrac
static const std::string center
unsigned int numFaces() const
XFEMCutElem3D(Elem *elem, const EFAElement3D *const CEMelem, unsigned int n_qpoints, unsigned int n_sides)
Constructor initializes XFEMCutElem3D object.