www.mooseframework.org
XFEMCutElem3D.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "XFEMCutElem3D.h"
11 
12 #include "EFANode.h"
13 #include "EFAFace.h"
14 #include "EFAFragment3D.h"
15 #include "EFAFuncs.h"
16 #include "XFEMFuncs.h"
17 #include "MooseError.h"
18 
19 #include "libmesh/mesh.h"
20 #include "libmesh/elem.h"
21 #include "libmesh/node.h"
22 
24  const EFAElement3D * const CEMelem,
25  unsigned int n_qpoints,
26  unsigned int n_sides)
27  : XFEMCutElem(elem, n_qpoints, n_sides), _efa_elem3d(CEMelem, true)
28 {
31 }
32 
34 
35 Point
36 XFEMCutElem3D::getNodeCoordinates(EFANode * CEMnode, MeshBase * displaced_mesh) const
37 {
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;
42 
43  _efa_elem3d.getMasterInfo(CEMnode, master_nodes, master_weights);
44  for (unsigned int i = 0; i < master_nodes.size(); ++i)
45  {
46  if (master_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
47  {
48  Node * node = _nodes[master_nodes[i]->id()];
49  if (displaced_mesh)
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);
53  }
54  else
55  mooseError("master nodes must be local");
56  }
57  for (unsigned int i = 0; i < master_nodes.size(); ++i)
58  node_coor += master_weights[i] * master_points[i];
59  return node_coor;
60 }
61 
62 void
64 {
65  Real frag_vol = 0.0;
66 
67  // collect fragment info needed by polyhedron_volume_3d()
68  std::vector<std::vector<unsigned int>> frag_face_indices;
69  std::vector<EFANode *> frag_nodes;
70  _efa_elem3d.getFragment(0)->getNodeInfo(frag_face_indices, frag_nodes);
71  int face_num = frag_face_indices.size();
72  int node_num = frag_nodes.size();
73 
74  int order_max = 0;
75  int * order = new int[face_num];
76  for (int i = 0; i < face_num; ++i)
77  {
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();
81  }
82 
83  double * coord = new double[3 * node_num];
84  for (unsigned int i = 0; i < frag_nodes.size(); ++i)
85  {
86  Point p = getNodeCoordinates(frag_nodes[i]);
87  coord[3 * i + 0] = p(0);
88  coord[3 * i + 1] = p(1);
89  coord[3 * i + 2] = p(2);
90  }
91 
92  int * node = new int[face_num * order_max];
93  Xfem::i4vec_zero(face_num * order_max, node);
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];
97 
98  // compute fragment volume and volume fraction
99  frag_vol = Xfem::polyhedron_volume_3d(coord, order_max, face_num, node, node_num, order);
100  _physical_volfrac = frag_vol / _elem_volume;
101 
102  delete[] order;
103  delete[] coord;
104  delete[] node;
105 }
106 
107 void
109 {
110  Real frag_surf = 0.0;
111 
112  std::vector<std::vector<unsigned int>> frag_face_indices;
113  std::vector<EFANode *> frag_nodes;
114  _efa_elem3d.getFragment(0)->getNodeInfo(frag_face_indices, frag_nodes);
115  int face_num = frag_face_indices.size();
116 
117  EFAFace * efa_face = _efa_elem3d.getFace(side);
118  bool contains_all = true;
119 
121  for (int i = 0; i < face_num; ++i)
122  {
123  contains_all = true;
124  for (unsigned int j = 0; j < frag_face_indices[i].size(); ++j)
125  {
126  EFANode * efa_node = frag_nodes[frag_face_indices[i][j]];
127  if (!efa_face->containsNode(efa_node))
128  {
129  contains_all = false;
130  break;
131  }
132  }
133 
134  if (contains_all)
135  {
136  for (unsigned int j = 0; j < frag_face_indices[i].size(); ++j)
137  {
138  unsigned int m = ((j + 1) == frag_face_indices[i].size()) ? 0 : j + 1;
139  Point edge_p1 = getNodeCoordinates(frag_nodes[frag_face_indices[i][j]]);
140  Point edge_p2 = getNodeCoordinates(frag_nodes[frag_face_indices[i][m]]);
141 
142  frag_surf += 0.5 * (edge_p1(0) - edge_p2(0)) * (edge_p1(1) + edge_p2(1));
143  }
144  _physical_areafrac[side] = std::abs(frag_surf) / _elem_side_area[side];
145  return;
146  }
147  }
148  _physical_areafrac[side] = 1.0;
149 }
150 
151 void
153 {
154  // TODO: 3D moment fitting method nod coded yet - use volume fraction for now
156 }
157 
158 Point
159 XFEMCutElem3D::getCutPlaneOrigin(unsigned int plane_id, MeshBase * displaced_mesh) const
160 {
161  Point orig(0.0, 0.0, 0.0);
162  std::vector<std::vector<EFANode *>> cut_plane_nodes;
163  for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
164  {
166  {
167  EFAFace * face = _efa_elem3d.getFragment(0)->getFace(i);
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);
172  }
173  }
174  if (cut_plane_nodes.size() == 0)
175  mooseError("no cut plane found in this element");
176  if (plane_id < cut_plane_nodes.size()) // valid plane_id
177  {
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));
181 
182  for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
183  orig += cut_plane_points[i];
184  orig /= cut_plane_points.size();
185  }
186  return orig;
187 }
188 
189 Point
190 XFEMCutElem3D::getCutPlaneNormal(unsigned int plane_id, MeshBase * displaced_mesh) const
191 {
192  Point normal(0.0, 0.0, 0.0);
193  std::vector<std::vector<EFANode *>> cut_plane_nodes;
194  for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
195  {
197  {
198  EFAFace * face = _efa_elem3d.getFragment(0)->getFace(i);
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);
203  }
204  }
205  if (cut_plane_nodes.size() == 0)
206  mooseError("no cut plane found in this element");
207  if (plane_id < cut_plane_nodes.size()) // valid plane_id
208  {
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));
212 
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();
217 
218  for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
219  {
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);
224  }
225  normal /= cut_plane_points.size();
226  }
227  Xfem::normalizePoint(normal);
228  return normal;
229 }
230 
231 void
233  Point & /*origin*/,
234  Point & /*direction*/) const
235 {
236  // TODO: not implemented for 3D
237  mooseError("getCrackTipOriginAndDirection not yet implemented for XFEMCutElem3D");
238 }
239 
240 void
241 XFEMCutElem3D::getFragmentFaces(std::vector<std::vector<Point>> & /*frag_faces*/,
242  MeshBase * /*displaced_mesh*/) const
243 {
244  // TODO: not implemented for 3D
245  mooseError("getFragmentFaces not yet implemented for XFEMCutElem3D");
246 }
247 
248 const EFAElement *
250 {
251  return &_efa_elem3d;
252 }
253 
254 unsigned int
256 {
257  unsigned int counter = 0;
258  for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
260  counter += 1;
261  return counter;
262 }
263 
264 void
266  Point & normal,
267  std::vector<Point> & intersectionPoints,
268  MeshBase * displaced_mesh) const
269 {
270  intersectionPoints.clear();
271  std::vector<std::vector<EFANode *>> cut_plane_nodes;
272  for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
273  {
275  {
276  EFAFace * face = _efa_elem3d.getFragment(0)->getFace(i);
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);
281  }
282  }
283  if (cut_plane_nodes.size() == 0)
284  mooseError("No cut plane found in this element");
285 
286  if (plane_id < cut_plane_nodes.size()) // valid plane_id
287  {
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);
291  }
292 
293  normal = getCutPlaneNormal(plane_id, displaced_mesh);
294 }
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[])
Definition: XFEMFuncs.C:489
EFAFace * getFace(unsigned int face_id) const
std::vector< Real > _elem_side_area
Definition: XFEMCutElem.h:49
EFAElement3D _efa_elem3d
Definition: XFEMCutElem3D.h:41
bool containsNode(const EFANode *node) const
Definition: EFAFace.C:362
virtual const EFAElement * getEFAElement() const
unsigned int _n_qpoints
Definition: XFEMCutElem.h:43
void mooseError(Args &&... args)
virtual void getIntersectionInfo(unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, MeshBase *displaced_mesh=nullptr) const
Real _physical_volfrac
Definition: XFEMCutElem.h:50
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
Definition: EFAFace.C:99
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[])
Definition: XFEMFuncs.C:586
EFAFace * getFace(unsigned int face_id) const
unsigned int numNodes() const
Definition: EFAFace.C:87
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
Definition: EFAElement3D.C:373
std::vector< Real > _new_weights
quadrature weights from volume fraction and moment fitting
Definition: XFEMCutElem.h:55
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=nullptr) const
Definition: XFEMCutElem3D.C:36
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void computePhysicalVolumeFraction()
Computes the volume fraction of the element fragment.
Definition: XFEMCutElem3D.C:63
std::vector< Node * > _nodes
Definition: XFEMCutElem.h:45
Real _elem_volume
Definition: XFEMCutElem.h:48
virtual unsigned int numCutPlanes() const
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
void normalizePoint(Point &p)
Definition: XFEMFuncs.C:621
virtual void getFragmentFaces(std::vector< std::vector< Point >> &frag_faces, MeshBase *displaced_mesh=nullptr) const
std::vector< Real > _physical_areafrac
Definition: XFEMCutElem.h:51
static const std::string center
Definition: NS.h:28
unsigned int numFaces() const
XFEMCutElem3D(Elem *elem, const EFAElement3D *const CEMelem, unsigned int n_qpoints, unsigned int n_sides)
Constructor initializes XFEMCutElem3D object.
Definition: XFEMCutElem3D.C:23