www.mooseframework.org
Public Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
XFEMCutElem3D Class Reference

#include <XFEMCutElem3D.h>

Inheritance diagram for XFEMCutElem3D:
[legend]

Public Member Functions

 XFEMCutElem3D (Elem *elem, const EFAElement3D *const CEMelem, unsigned int n_qpoints)
 
 ~XFEMCutElem3D ()
 
virtual void computePhysicalVolumeFraction ()
 
virtual void computeMomentFittingWeights ()
 
virtual Point getCutPlaneOrigin (unsigned int plane_id, MeshBase *displaced_mesh=NULL) const
 
virtual Point getCutPlaneNormal (unsigned int plane_id, MeshBase *displaced_mesh=NULL) const
 
virtual void getCrackTipOriginAndDirection (unsigned tip_id, Point &origin, Point &direction) const
 
virtual void getFragmentFaces (std::vector< std::vector< Point >> &frag_faces, MeshBase *displaced_mesh=NULL) const
 
virtual const EFAElementgetEFAElement () const
 
virtual unsigned int numCutPlanes () const
 
virtual void getIntersectionInfo (unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, MeshBase *displaced_mesh=NULL) const
 
void setQuadraturePointsAndWeights (const std::vector< Point > &qp_points, const std::vector< Real > &qp_weights)
 
Real getPhysicalVolumeFraction () const
 
Real getMomentFittingWeight (unsigned int i_qp) const
 
void getWeightMultipliers (MooseArray< Real > &weights, QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
 
void computeXFEMWeights (QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
 
bool isPointPhysical (const Point &p) const
 

Protected Attributes

unsigned int _n_nodes
 
unsigned int _n_qpoints
 
std::vector< Node * > _nodes
 
std::vector< Point > _qp_points
 
std::vector< Real > _qp_weights
 
Real _elem_volume
 
Real _physical_volfrac
 
bool _have_weights
 
std::vector< Real > _new_weights
 

Private Member Functions

virtual Point getNodeCoordinates (EFANode *node, MeshBase *displaced_mesh=NULL) const
 

Private Attributes

EFAElement3D _efa_elem3d
 

Detailed Description

Definition at line 23 of file XFEMCutElem3D.h.

Constructor & Destructor Documentation

XFEMCutElem3D::XFEMCutElem3D ( Elem *  elem,
const EFAElement3D *const  CEMelem,
unsigned int  n_qpoints 
)

Definition at line 21 of file XFEMCutElem3D.C.

24  : XFEMCutElem(elem, n_qpoints), _efa_elem3d(CEMelem, true)
25 {
28 }
EFAElement3D _efa_elem3d
Definition: XFEMCutElem3D.h:30
virtual void computeMomentFittingWeights()
virtual void computePhysicalVolumeFraction()
Definition: XFEMCutElem3D.C:60
XFEMCutElem(Elem *elem, unsigned int n_qpoints)
Definition: XFEMCutElem.C:20
XFEMCutElem3D::~XFEMCutElem3D ( )

Definition at line 30 of file XFEMCutElem3D.C.

30 {}

Member Function Documentation

void XFEMCutElem3D::computeMomentFittingWeights ( )
virtual

Implements XFEMCutElem.

Definition at line 105 of file XFEMCutElem3D.C.

Referenced by XFEMCutElem3D().

106 {
107  // TODO: 3D moment fitting method nod coded yet - use volume fraction for now
109 }
unsigned int _n_qpoints
Definition: XFEMCutElem.h:36
Real _physical_volfrac
Definition: XFEMCutElem.h:41
std::vector< Real > _new_weights
Definition: XFEMCutElem.h:43
void XFEMCutElem3D::computePhysicalVolumeFraction ( )
virtual

Implements XFEMCutElem.

Definition at line 60 of file XFEMCutElem3D.C.

Referenced by XFEMCutElem3D().

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 }
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
EFAElement3D _efa_elem3d
Definition: XFEMCutElem3D.h:30
Real _physical_volfrac
Definition: XFEMCutElem.h:41
EFAFragment3D * getFragment(unsigned int frag_id) const
void i4vec_zero(int n, int a[])
Definition: XFEMFuncs.C:592
Real _elem_volume
Definition: XFEMCutElem.h:40
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem3D.C:33
void XFEMCutElem::computeXFEMWeights ( QBase *  qrule,
Xfem::XFEM_QRULE  xfem_qrule,
const MooseArray< Point > &  q_points 
)
inherited

Definition at line 51 of file XFEMCutElem.C.

Referenced by XFEMCutElem::getWeightMultipliers().

54 {
55  _new_weights.clear();
56 
57  _new_weights.resize(qrule->n_points(), 1.0);
58 
59  switch (xfem_qrule)
60  {
61  case Xfem::VOLFRAC:
62  {
64  Real volfrac = getPhysicalVolumeFraction();
65  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
66  _new_weights[qp] = volfrac;
67  break;
68  }
70  {
71  // These are the coordinates in parametric coordinates
72  _qp_points = qrule->get_points();
73  _qp_weights = qrule->get_weights();
74 
76 
77  // Blend weights from moment fitting and volume fraction to avoid negative weights
78  Real alpha = 1.0;
79  if (*std::min_element(_new_weights.begin(), _new_weights.end()) < 0.0)
80  {
81  // One or more of the weights computed by moment fitting is negative.
82  // Blend moment and volume fraction weights to keep them nonnegative.
83  // Find the largest value of alpha that will keep all of the weights nonnegative.
84  for (unsigned int i = 0; i < _n_qpoints; ++i)
85  {
86  const Real denominator = _physical_volfrac - _new_weights[i];
87  if (denominator > 0.0) // Negative values would give a negative value for alpha, which
88  // must be between 0 and 1.
89  {
90  const Real alpha_i = _physical_volfrac / denominator;
91  if (alpha_i < alpha)
92  alpha = alpha_i;
93  }
94  }
95  for (unsigned int i = 0; i < _n_qpoints; ++i)
96  {
97  _new_weights[i] = alpha * _new_weights[i] + (1.0 - alpha) * _physical_volfrac;
98  if (_new_weights[i] < 0.0) // We can end up with small (roundoff) negative weights
99  _new_weights[i] = 0.0;
100  }
101  }
102  break;
103  }
104  case Xfem::DIRECT: // remove q-points outside the partial element's physical domain
105  {
106  bool nonzero = false;
107  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
108  {
109  // q_points contains quadrature point locations in physical coordinates
110  if (isPointPhysical(q_points[qp]))
111  {
112  nonzero = true;
113  _new_weights[qp] = 1.0;
114  }
115  else
116  _new_weights[qp] = 0.0;
117  }
118  if (!nonzero)
119  {
120  // Set the weights to a small value (1e-3) to avoid having DOFs
121  // with zeros on the diagonal, which occurs for nodes that are
122  // connected to no physical material.
123  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
124  _new_weights[qp] = 1e-3;
125  }
126  break;
127  }
128  default:
129  mooseError("Undefined option for XFEM_QRULE");
130  }
131  _have_weights = true;
132 }
unsigned int _n_qpoints
Definition: XFEMCutElem.h:36
std::vector< Point > _qp_points
Definition: XFEMCutElem.h:38
Real _physical_volfrac
Definition: XFEMCutElem.h:41
Real getPhysicalVolumeFraction() const
Definition: XFEMCutElem.C:31
bool isPointPhysical(const Point &p) const
Definition: XFEMCutElem.C:135
std::vector< Real > _qp_weights
Definition: XFEMCutElem.h:39
std::vector< Real > _new_weights
Definition: XFEMCutElem.h:43
bool _have_weights
Definition: XFEMCutElem.h:42
virtual void computePhysicalVolumeFraction()=0
virtual void computeMomentFittingWeights()=0
void XFEMCutElem3D::getCrackTipOriginAndDirection ( unsigned  tip_id,
Point &  origin,
Point &  direction 
) const
virtual

Implements XFEMCutElem.

Definition at line 185 of file XFEMCutElem3D.C.

188 {
189  // TODO: not implemented for 3D
190  mooseError("getCrackTipOriginAndDirection not yet implemented for XFEMCutElem3D");
191 }
Point XFEMCutElem3D::getCutPlaneNormal ( unsigned int  plane_id,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 143 of file XFEMCutElem3D.C.

Referenced by getIntersectionInfo().

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 }
unsigned int numNodes() const
Definition: EFAFace.C:85
EFAElement3D _efa_elem3d
Definition: XFEMCutElem3D.h:30
bool isFaceInterior(unsigned int face_id) const
EFAFragment3D * getFragment(unsigned int frag_id) const
EFAFace * getFace(unsigned int face_id) const
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
Point XFEMCutElem3D::getCutPlaneOrigin ( unsigned int  plane_id,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 112 of file XFEMCutElem3D.C.

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 }
unsigned int numNodes() const
Definition: EFAFace.C:85
EFAElement3D _efa_elem3d
Definition: XFEMCutElem3D.h:30
bool isFaceInterior(unsigned int face_id) const
EFAFragment3D * getFragment(unsigned int frag_id) const
EFAFace * getFace(unsigned int face_id) const
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
const EFAElement * XFEMCutElem3D::getEFAElement ( ) const
virtual

Implements XFEMCutElem.

Definition at line 202 of file XFEMCutElem3D.C.

203 {
204  return &_efa_elem3d;
205 }
EFAElement3D _efa_elem3d
Definition: XFEMCutElem3D.h:30
void XFEMCutElem3D::getFragmentFaces ( std::vector< std::vector< Point >> &  frag_faces,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 194 of file XFEMCutElem3D.C.

196 {
197  // TODO: not implemented for 3D
198  mooseError("getFragmentFaces not yet implemented for XFEMCutElem3D");
199 }
void XFEMCutElem3D::getIntersectionInfo ( unsigned int  plane_id,
Point &  normal,
std::vector< Point > &  intersectionPoints,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 218 of file XFEMCutElem3D.C.

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
EFAElement3D _efa_elem3d
Definition: XFEMCutElem3D.h:30
bool isFaceInterior(unsigned int face_id) const
EFAFragment3D * getFragment(unsigned int frag_id) const
EFAFace * getFace(unsigned int face_id) const
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const
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
Real XFEMCutElem::getMomentFittingWeight ( unsigned int  i_qp) const
inherited
Point XFEMCutElem3D::getNodeCoordinates ( EFANode node,
MeshBase *  displaced_mesh = NULL 
) const
privatevirtual

Implements XFEMCutElem.

Definition at line 33 of file XFEMCutElem3D.C.

Referenced by computePhysicalVolumeFraction(), getCutPlaneNormal(), getCutPlaneOrigin(), and getIntersectionInfo().

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 }
EFAElement3D _efa_elem3d
Definition: XFEMCutElem3D.h:30
virtual void getMasterInfo(EFANode *node, std::vector< EFANode * > &master_nodes, std::vector< double > &master_weights) const
Definition: EFAElement3D.C:367
std::vector< Node * > _nodes
Definition: XFEMCutElem.h:37
Real XFEMCutElem::getPhysicalVolumeFraction ( ) const
inherited

Definition at line 31 of file XFEMCutElem.C.

Referenced by XFEMCutElem::computeXFEMWeights(), and XFEM::getPhysicalVolumeFraction().

32 {
33  return _physical_volfrac;
34 }
Real _physical_volfrac
Definition: XFEMCutElem.h:41
void XFEMCutElem::getWeightMultipliers ( MooseArray< Real > &  weights,
QBase *  qrule,
Xfem::XFEM_QRULE  xfem_qrule,
const MooseArray< Point > &  q_points 
)
inherited

Definition at line 37 of file XFEMCutElem.C.

Referenced by XFEM::getXFEMWeights().

41 {
42  if (!_have_weights)
43  computeXFEMWeights(qrule, xfem_qrule, q_points);
44 
45  weights.resize(_new_weights.size());
46  for (unsigned int qp = 0; qp < _new_weights.size(); ++qp)
47  weights[qp] = _new_weights[qp];
48 }
std::vector< Real > _new_weights
Definition: XFEMCutElem.h:43
bool _have_weights
Definition: XFEMCutElem.h:42
void computeXFEMWeights(QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
Definition: XFEMCutElem.C:51
bool XFEMCutElem::isPointPhysical ( const Point &  p) const
inherited

Definition at line 135 of file XFEMCutElem.C.

Referenced by XFEMCutElem::computeXFEMWeights().

136 {
137  // determine whether the point is inside the physical domain of a partial element
138  bool physical_flag = true;
139  unsigned int n_cut_planes = numCutPlanes();
140  for (unsigned int plane_id = 0; plane_id < n_cut_planes; ++plane_id)
141  {
142  Point origin = getCutPlaneOrigin(plane_id);
143  Point normal = getCutPlaneNormal(plane_id);
144  Point origin2qp = p - origin;
145  if (origin2qp * normal > 0.0)
146  {
147  physical_flag = false; // Point outside pysical domain
148  break;
149  }
150  }
151  return physical_flag;
152 }
virtual Point getCutPlaneOrigin(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const =0
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const =0
virtual unsigned int numCutPlanes() const =0
unsigned int XFEMCutElem3D::numCutPlanes ( ) const
virtual

Implements XFEMCutElem.

Definition at line 208 of file XFEMCutElem3D.C.

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 }
EFAElement3D _efa_elem3d
Definition: XFEMCutElem3D.h:30
bool isFaceInterior(unsigned int face_id) const
EFAFragment3D * getFragment(unsigned int frag_id) const
unsigned int numFaces() const
static unsigned int counter
void XFEMCutElem::setQuadraturePointsAndWeights ( const std::vector< Point > &  qp_points,
const std::vector< Real > &  qp_weights 
)
inherited

Member Data Documentation

EFAElement3D XFEMCutElem3D::_efa_elem3d
private
Real XFEMCutElem::_elem_volume
protectedinherited
bool XFEMCutElem::_have_weights
protectedinherited
unsigned int XFEMCutElem::_n_nodes
protectedinherited
unsigned int XFEMCutElem::_n_qpoints
protectedinherited
std::vector<Real> XFEMCutElem::_new_weights
protectedinherited
std::vector<Node *> XFEMCutElem::_nodes
protectedinherited
Real XFEMCutElem::_physical_volfrac
protectedinherited
std::vector<Point> XFEMCutElem::_qp_points
protectedinherited
std::vector<Real> XFEMCutElem::_qp_weights
protectedinherited

The documentation for this class was generated from the following files: