www.mooseframework.org
XFEMCutElem3D.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "XFEMCutElem3D.h"
9 
10 #include "EFANode.h"
11 #include "EFAFace.h"
12 #include "EFAFragment3D.h"
13 #include "EFAFuncs.h"
14 #include "XFEMFuncs.h"
15 #include "MooseError.h"
16 
17 #include "libmesh/mesh.h"
18 #include "libmesh/elem.h"
19 #include "libmesh/node.h"
20 
22  const EFAElement3D * const CEMelem,
23  unsigned int n_qpoints)
24  : XFEMCutElem(elem, n_qpoints), _efa_elem3d(CEMelem, true)
25 {
28 }
29 
31 
32 Point
33 XFEMCutElem3D::getNodeCoordinates(EFANode * CEMnode, MeshBase * displaced_mesh) const
34 {
35  Point node_coor(0.0, 0.0, 0.0);
36  std::vector<EFANode *> master_nodes;
37  std::vector<Point> master_points;
38  std::vector<double> master_weights;
39 
40  _efa_elem3d.getMasterInfo(CEMnode, master_nodes, master_weights);
41  for (unsigned int i = 0; i < master_nodes.size(); ++i)
42  {
43  if (master_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
44  {
45  Node * node = _nodes[master_nodes[i]->id()];
46  if (displaced_mesh)
47  node = displaced_mesh->node_ptr(node->id());
48  Point node_p((*node)(0), (*node)(1), (*node)(2));
49  master_points.push_back(node_p);
50  }
51  else
52  mooseError("master nodes must be local");
53  }
54  for (unsigned int i = 0; i < master_nodes.size(); ++i)
55  node_coor += master_weights[i] * master_points[i];
56  return node_coor;
57 }
58 
59 void
61 {
62  Real frag_vol = 0.0;
63 
64  // collect fragment info needed by polyhedron_volume_3d()
65  std::vector<std::vector<unsigned int>> frag_face_indices;
66  std::vector<EFANode *> frag_nodes;
67  _efa_elem3d.getFragment(0)->getNodeInfo(frag_face_indices, frag_nodes);
68  int face_num = frag_face_indices.size();
69  int node_num = frag_nodes.size();
70 
71  int order_max = 0;
72  int * order = new int[face_num];
73  for (int i = 0; i < face_num; ++i)
74  {
75  if (frag_face_indices[i].size() > (unsigned int)order_max)
76  order_max = frag_face_indices[i].size();
77  order[i] = frag_face_indices[i].size();
78  }
79 
80  double * coord = new double[3 * node_num];
81  for (unsigned int i = 0; i < frag_nodes.size(); ++i)
82  {
83  Point p = getNodeCoordinates(frag_nodes[i]);
84  coord[3 * i + 0] = p(0);
85  coord[3 * i + 1] = p(1);
86  coord[3 * i + 2] = p(2);
87  }
88 
89  int * node = new int[face_num * order_max];
90  Xfem::i4vec_zero(face_num * order_max, node);
91  for (int i = 0; i < face_num; ++i)
92  for (unsigned int j = 0; j < frag_face_indices[i].size(); ++j)
93  node[order_max * i + j] = frag_face_indices[i][j];
94 
95  // compute fragment volume and volume fraction
96  frag_vol = Xfem::polyhedron_volume_3d(coord, order_max, face_num, node, node_num, order);
97  _physical_volfrac = frag_vol / _elem_volume;
98 
99  delete[] order;
100  delete[] coord;
101  delete[] node;
102 }
103 
104 void
106 {
107  // TODO: 3D moment fitting method nod coded yet - use volume fraction for now
109 }
110 
111 Point
112 XFEMCutElem3D::getCutPlaneOrigin(unsigned int plane_id, MeshBase * displaced_mesh) const
113 {
114  Point orig(0.0, 0.0, 0.0);
115  std::vector<std::vector<EFANode *>> cut_plane_nodes;
116  for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
117  {
119  {
120  EFAFace * face = _efa_elem3d.getFragment(0)->getFace(i);
121  std::vector<EFANode *> node_line;
122  for (unsigned int j = 0; j < face->numNodes(); ++j)
123  node_line.push_back(face->getNode(j));
124  cut_plane_nodes.push_back(node_line);
125  }
126  }
127  if (cut_plane_nodes.size() == 0)
128  mooseError("no cut plane found in this element");
129  if (plane_id < cut_plane_nodes.size()) // valid plane_id
130  {
131  std::vector<Point> cut_plane_points;
132  for (unsigned int i = 0; i < cut_plane_nodes[plane_id].size(); ++i)
133  cut_plane_points.push_back(getNodeCoordinates(cut_plane_nodes[plane_id][i], displaced_mesh));
134 
135  for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
136  orig += cut_plane_points[i];
137  orig /= cut_plane_points.size();
138  }
139  return orig;
140 }
141 
142 Point
143 XFEMCutElem3D::getCutPlaneNormal(unsigned int plane_id, MeshBase * displaced_mesh) const
144 {
145  Point normal(0.0, 0.0, 0.0);
146  std::vector<std::vector<EFANode *>> cut_plane_nodes;
147  for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
148  {
150  {
151  EFAFace * face = _efa_elem3d.getFragment(0)->getFace(i);
152  std::vector<EFANode *> node_line;
153  for (unsigned int j = 0; j < face->numNodes(); ++j)
154  node_line.push_back(face->getNode(j));
155  cut_plane_nodes.push_back(node_line);
156  }
157  }
158  if (cut_plane_nodes.size() == 0)
159  mooseError("no cut plane found in this element");
160  if (plane_id < cut_plane_nodes.size()) // valid plane_id
161  {
162  std::vector<Point> cut_plane_points;
163  for (unsigned int i = 0; i < cut_plane_nodes[plane_id].size(); ++i)
164  cut_plane_points.push_back(getNodeCoordinates(cut_plane_nodes[plane_id][i], displaced_mesh));
165 
166  Point center(0.0, 0.0, 0.0);
167  for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
168  center += cut_plane_points[i];
169  center /= cut_plane_points.size();
170 
171  for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
172  {
173  unsigned int iplus1 = i < cut_plane_points.size() - 1 ? i + 1 : 0;
174  Point ray1 = cut_plane_points[i] - center;
175  Point ray2 = cut_plane_points[iplus1] - center;
176  normal += ray1.cross(ray2);
177  }
178  normal /= cut_plane_points.size();
179  }
180  Xfem::normalizePoint(normal);
181  return normal;
182 }
183 
184 void
186  Point & /*origin*/,
187  Point & /*direction*/) const
188 {
189  // TODO: not implemented for 3D
190  mooseError("getCrackTipOriginAndDirection not yet implemented for XFEMCutElem3D");
191 }
192 
193 void
194 XFEMCutElem3D::getFragmentFaces(std::vector<std::vector<Point>> & /*frag_faces*/,
195  MeshBase * /*displaced_mesh*/) const
196 {
197  // TODO: not implemented for 3D
198  mooseError("getFragmentFaces not yet implemented for XFEMCutElem3D");
199 }
200 
201 const EFAElement *
203 {
204  return &_efa_elem3d;
205 }
206 
207 unsigned int
209 {
210  unsigned int counter = 0;
211  for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
213  counter += 1;
214  return counter;
215 }
216 
217 void
219  Point & normal,
220  std::vector<Point> & intersectionPoints,
221  MeshBase * displaced_mesh) const
222 {
223  intersectionPoints.clear();
224  std::vector<std::vector<EFANode *>> cut_plane_nodes;
225  for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
226  {
228  {
229  EFAFace * face = _efa_elem3d.getFragment(0)->getFace(i);
230  std::vector<EFANode *> node_line;
231  for (unsigned int j = 0; j < face->numNodes(); ++j)
232  node_line.push_back(face->getNode(j));
233  cut_plane_nodes.push_back(node_line);
234  }
235  }
236  if (cut_plane_nodes.size() == 0)
237  mooseError("No cut plane found in this element");
238 
239  if (plane_id < cut_plane_nodes.size()) // valid plane_id
240  {
241  intersectionPoints.resize(cut_plane_nodes[plane_id].size());
242  for (unsigned int i = 0; i < cut_plane_nodes[plane_id].size(); ++i)
243  intersectionPoints[i] = getNodeCoordinates(cut_plane_nodes[plane_id][i], displaced_mesh);
244  }
245 
246  normal = getCutPlaneNormal(plane_id, displaced_mesh);
247 }
unsigned int numNodes() const
Definition: EFAFace.C:85
void getNodeInfo(std::vector< std::vector< unsigned int >> &face_node_indices, std::vector< EFANode * > &nodes) const
double polyhedron_volume_3d(double coord[], int order_max, int face_num, int node[], int node_num, int order[])
Definition: XFEMFuncs.C:494
virtual unsigned int numCutPlanes() const
EFAElement3D _efa_elem3d
Definition: XFEMCutElem3D.h:30
bool isFaceInterior(unsigned int face_id) const
unsigned int _n_qpoints
Definition: XFEMCutElem.h:36
virtual void getIntersectionInfo(unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, MeshBase *displaced_mesh=NULL) const
Real _physical_volfrac
Definition: XFEMCutElem.h:41
virtual void computeMomentFittingWeights()
EFAFragment3D * getFragment(unsigned int frag_id) const
void i4vec_zero(int n, int a[])
Definition: XFEMFuncs.C:592
virtual Point getCutPlaneOrigin(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const
virtual void getCrackTipOriginAndDirection(unsigned tip_id, Point &origin, Point &direction) const
EFAFace * getFace(unsigned int face_id) const
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const
std::vector< Real > _new_weights
Definition: XFEMCutElem.h:43
virtual void getMasterInfo(EFANode *node, std::vector< EFANode * > &master_nodes, std::vector< double > &master_weights) const
Definition: EFAElement3D.C:367
virtual void getFragmentFaces(std::vector< std::vector< Point >> &frag_faces, MeshBase *displaced_mesh=NULL) const
XFEMCutElem3D(Elem *elem, const EFAElement3D *const CEMelem, unsigned int n_qpoints)
Definition: XFEMCutElem3D.C:21
virtual void computePhysicalVolumeFraction()
Definition: XFEMCutElem3D.C:60
std::vector< Node * > _nodes
Definition: XFEMCutElem.h:37
virtual const EFAElement * getEFAElement() const
Real _elem_volume
Definition: XFEMCutElem.h:40
void normalizePoint(Point &p)
Definition: XFEMFuncs.C:628
unsigned int numFaces() const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem3D.C:33
EFANode * getNode(unsigned int node_id) const
Definition: EFAFace.C:97
static unsigned int counter