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

#include <XFEMCutElem2D.h>

Inheritance diagram for XFEMCutElem2D:
[legend]

Public Member Functions

 XFEMCutElem2D (Elem *elem, const EFAElement2D *const CEMelem, unsigned int n_qpoints, unsigned int n_sides)
 Constructor initializes XFEMCutElem2D object. More...
 
 ~XFEMCutElem2D ()
 
virtual void computePhysicalVolumeFraction ()
 Computes the volume fraction of the element fragment. More...
 
virtual void computePhysicalFaceAreaFraction (unsigned int side)
 Computes the surface area fraction of the element side. More...
 
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
 Returns the volume fraction of the element fragment. More...
 
Real getPhysicalFaceAreaFraction (unsigned int side) const
 Returns the surface area fraction of the element side. More...
 
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 getFaceWeightMultipliers (MooseArray< Real > &face_weights, QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points, unsigned int side)
 
void computeXFEMWeights (QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
 Computes integration weights for the cut element. More...
 
void computeXFEMFaceWeights (QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points, unsigned int side)
 Computes face integration weights for the cut element side. More...
 
bool isPointPhysical (const Point &p) const
 

Protected Attributes

unsigned int _n_nodes
 
unsigned int _n_qpoints
 
unsigned int _n_sides
 
std::vector< Node * > _nodes
 
std::vector< Point > _qp_points
 
std::vector< Real > _qp_weights
 
Real _elem_volume
 
std::vector< Real > _elem_side_area
 
Real _physical_volfrac
 
std::vector< Real > _physical_areafrac
 
bool _have_weights
 
std::vector< bool > _have_face_weights
 
std::vector< Real > _new_weights
 quadrature weights from volume fraction and moment fitting More...
 
std::vector< std::vector< Real > > _new_face_weights
 face quadrature weights from surface area fraction More...
 

Private Member Functions

virtual Point getNodeCoordinates (EFANode *node, MeshBase *displaced_mesh=NULL) const
 
void getPhysicalQuadraturePoints (std::vector< std::vector< Real >> &tsg)
 
void solveMomentFitting (unsigned int nen, unsigned int nqp, std::vector< Point > &elem_nodes, std::vector< std::vector< Real >> &tsg, std::vector< std::vector< Real >> &wsg)
 

Private Attributes

EFAElement2D _efa_elem2d
 

Detailed Description

Definition at line 25 of file XFEMCutElem2D.h.

Constructor & Destructor Documentation

◆ XFEMCutElem2D()

XFEMCutElem2D::XFEMCutElem2D ( Elem *  elem,
const EFAElement2D *const  CEMelem,
unsigned int  n_qpoints,
unsigned int  n_sides 
)

Constructor initializes XFEMCutElem2D object.

Parameters
elemThe element on which XFEMCutElem2D is built
CEMelemThe EFAFragment2D object that belongs to XFEMCutElem2D
n_qpointsThe number of quadrature points
n_sidesThe number of sides which the element has

Definition at line 24 of file XFEMCutElem2D.C.

28  : XFEMCutElem(elem, n_qpoints, n_sides), _efa_elem2d(CEMelem, true)
29 {
31 }
virtual void computePhysicalVolumeFraction()
Computes the volume fraction of the element fragment.
Definition: XFEMCutElem2D.C:63
XFEMCutElem(Elem *elem, unsigned int n_qpoints, unsigned int n_sides)
Constructor initializes XFEMCutElem object.
Definition: XFEMCutElem.C:22
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ ~XFEMCutElem2D()

XFEMCutElem2D::~XFEMCutElem2D ( )

Definition at line 33 of file XFEMCutElem2D.C.

33 {}

Member Function Documentation

◆ computeMomentFittingWeights()

void XFEMCutElem2D::computeMomentFittingWeights ( )
virtual

Implements XFEMCutElem.

Definition at line 105 of file XFEMCutElem2D.C.

106 {
107  // Purpose: calculate new weights via moment-fitting method
108  std::vector<Point> elem_nodes(_n_nodes, Point(0.0, 0.0, 0.0));
109  std::vector<std::vector<Real>> wsg;
110 
111  for (unsigned int i = 0; i < _n_nodes; ++i)
112  elem_nodes[i] = (*_nodes[i]);
113 
114  if (_efa_elem2d.isPartial() && _n_qpoints <= 6) // ONLY work for <= 6 q_points
115  {
116  std::vector<std::vector<Real>> tsg;
117  getPhysicalQuadraturePoints(tsg); // get tsg - QPs within partial element
119  _n_nodes, _n_qpoints, elem_nodes, tsg, wsg); // get wsg - QPs from moment-fitting
120  _new_weights.resize(wsg.size(), 1.0);
121  for (unsigned int i = 0; i < wsg.size(); ++i)
122  _new_weights[i] = wsg[i][2]; // weight multiplier
123  }
124  else
126 }
unsigned int _n_qpoints
Definition: XFEMCutElem.h:44
Real _physical_volfrac
Definition: XFEMCutElem.h:51
void solveMomentFitting(unsigned int nen, unsigned int nqp, std::vector< Point > &elem_nodes, std::vector< std::vector< Real >> &tsg, std::vector< std::vector< Real >> &wsg)
std::vector< Real > _new_weights
quadrature weights from volume fraction and moment fitting
Definition: XFEMCutElem.h:56
void getPhysicalQuadraturePoints(std::vector< std::vector< Real >> &tsg)
unsigned int _n_nodes
Definition: XFEMCutElem.h:43
std::vector< Node * > _nodes
Definition: XFEMCutElem.h:46
virtual bool isPartial() const
Definition: EFAElement2D.C:208
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ computePhysicalFaceAreaFraction()

void XFEMCutElem2D::computePhysicalFaceAreaFraction ( unsigned int  side)
virtual

Computes the surface area fraction of the element side.

Parameters
sideThe side of the element

find a fragment edge which is covered by element side

Implements XFEMCutElem.

Definition at line 80 of file XFEMCutElem2D.C.

81 {
82  Real frag_surf = 0.0;
83 
84  EFAEdge * edge = _efa_elem2d.getEdge(side);
85 
86  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
87  {
88  EFANode * node_1 = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
89  EFANode * node_2 = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
90 
92  if (edge->containsNode(node_1) && edge->containsNode(node_2))
93  {
94  Point edge_p1 = getNodeCoordinates(node_1);
95  Point edge_p2 = getNodeCoordinates(node_2);
96  frag_surf = (edge_p1 - edge_p2).norm();
97  _physical_areafrac[side] = frag_surf / _elem_side_area[side];
98  return;
99  }
100  }
101  _physical_areafrac[side] = 1.0;
102 }
std::vector< Real > _elem_side_area
Definition: XFEMCutElem.h:50
unsigned int numEdges() const
EFAEdge * getEdge(unsigned int edge_id) const
EFAFragment2D * getFragment(unsigned int frag_id) const
bool containsNode(const EFANode *node) const
Definition: EFAEdge.C:377
std::vector< Real > _physical_areafrac
Definition: XFEMCutElem.h:52
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:36
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:179
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ computePhysicalVolumeFraction()

void XFEMCutElem2D::computePhysicalVolumeFraction ( )
virtual

Computes the volume fraction of the element fragment.

Implements XFEMCutElem.

Definition at line 63 of file XFEMCutElem2D.C.

Referenced by XFEMCutElem2D().

64 {
65  Real frag_vol = 0.0;
66 
67  // Calculate area of entire element and fragment using the formula:
68  // A = 1/2 sum_{i=0}^{n-1} (x_i y_{i+1} - x_{i+1} y{i})
69 
70  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
71  {
72  Point edge_p1 = getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(0));
73  Point edge_p2 = getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(1));
74  frag_vol += 0.5 * (edge_p1(0) - edge_p2(0)) * (edge_p1(1) + edge_p2(1));
75  }
76  _physical_volfrac = frag_vol / _elem_volume;
77 }
unsigned int numEdges() const
Real _physical_volfrac
Definition: XFEMCutElem.h:51
EFAFragment2D * getFragment(unsigned int frag_id) const
Real _elem_volume
Definition: XFEMCutElem.h:49
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:36
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:179
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ computeXFEMFaceWeights()

void XFEMCutElem::computeXFEMFaceWeights ( QBase *  qrule,
Xfem::XFEM_QRULE  xfem_qrule,
const MooseArray< Point > &  q_points,
unsigned int  side 
)
inherited

Computes face integration weights for the cut element side.

Parameters
qruleThe standard MOOSE face quadrature rule
xfem_qruleThe integration scheme for the cut element (We use surface area fraction only)
q_pointsThe quadrature points for the element side
sideThe side of the element

Definition at line 90 of file XFEMCutElem.C.

Referenced by XFEMCutElem::getFaceWeightMultipliers().

94 {
95  _new_face_weights[side].clear();
96  _new_face_weights[side].resize(qrule->n_points(), 1.0);
97 
99  Real surffrac = getPhysicalFaceAreaFraction(side);
100  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
101  _new_face_weights[side][qp] = surffrac;
102 
103  _have_face_weights[side] = true;
104 }
std::vector< std::vector< Real > > _new_face_weights
face quadrature weights from surface area fraction
Definition: XFEMCutElem.h:58
std::vector< bool > _have_face_weights
Definition: XFEMCutElem.h:54
virtual void computePhysicalFaceAreaFraction(unsigned int side)=0
Computes the surface area fraction of the element side.
Real getPhysicalFaceAreaFraction(unsigned int side) const
Returns the surface area fraction of the element side.
Definition: XFEMCutElem.C:55

◆ computeXFEMWeights()

void XFEMCutElem::computeXFEMWeights ( QBase *  qrule,
Xfem::XFEM_QRULE  xfem_qrule,
const MooseArray< Point > &  q_points 
)
inherited

Computes integration weights for the cut element.

Parameters
qruleThe standard MOOSE quadrature rule
xfem_qruleThe integration scheme for the cut element
q_pointsThe quadrature points for the element

Definition at line 107 of file XFEMCutElem.C.

Referenced by XFEMCutElem::getWeightMultipliers().

110 {
111  _new_weights.clear();
112 
113  _new_weights.resize(qrule->n_points(), 1.0);
114 
115  switch (xfem_qrule)
116  {
117  case Xfem::VOLFRAC:
118  {
120  Real volfrac = getPhysicalVolumeFraction();
121  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
122  _new_weights[qp] = volfrac;
123  break;
124  }
126  {
127  // These are the coordinates in parametric coordinates
128  _qp_points = qrule->get_points();
129  _qp_weights = qrule->get_weights();
130 
132 
133  // Blend weights from moment fitting and volume fraction to avoid negative weights
134  Real alpha = 1.0;
135  if (*std::min_element(_new_weights.begin(), _new_weights.end()) < 0.0)
136  {
137  // One or more of the weights computed by moment fitting is negative.
138  // Blend moment and volume fraction weights to keep them nonnegative.
139  // Find the largest value of alpha that will keep all of the weights nonnegative.
140  for (unsigned int i = 0; i < _n_qpoints; ++i)
141  {
142  const Real denominator = _physical_volfrac - _new_weights[i];
143  if (denominator > 0.0) // Negative values would give a negative value for alpha, which
144  // must be between 0 and 1.
145  {
146  const Real alpha_i = _physical_volfrac / denominator;
147  if (alpha_i < alpha)
148  alpha = alpha_i;
149  }
150  }
151  for (unsigned int i = 0; i < _n_qpoints; ++i)
152  {
153  _new_weights[i] = alpha * _new_weights[i] + (1.0 - alpha) * _physical_volfrac;
154  if (_new_weights[i] < 0.0) // We can end up with small (roundoff) negative weights
155  _new_weights[i] = 0.0;
156  }
157  }
158  break;
159  }
160  case Xfem::DIRECT: // remove q-points outside the partial element's physical domain
161  {
162  bool nonzero = false;
163  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
164  {
165  // q_points contains quadrature point locations in physical coordinates
166  if (isPointPhysical(q_points[qp]))
167  {
168  nonzero = true;
169  _new_weights[qp] = 1.0;
170  }
171  else
172  _new_weights[qp] = 0.0;
173  }
174  if (!nonzero)
175  {
176  // Set the weights to a small value (1e-3) to avoid having DOFs
177  // with zeros on the diagonal, which occurs for nodes that are
178  // connected to no physical material.
179  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
180  _new_weights[qp] = 1e-3;
181  }
182  break;
183  }
184  default:
185  mooseError("Undefined option for XFEM_QRULE");
186  }
187  _have_weights = true;
188 }
unsigned int _n_qpoints
Definition: XFEMCutElem.h:44
std::vector< Point > _qp_points
Definition: XFEMCutElem.h:47
Real getPhysicalVolumeFraction() const
Returns the volume fraction of the element fragment.
Definition: XFEMCutElem.C:49
Real _physical_volfrac
Definition: XFEMCutElem.h:51
std::vector< Real > _qp_weights
Definition: XFEMCutElem.h:48
std::vector< Real > _new_weights
quadrature weights from volume fraction and moment fitting
Definition: XFEMCutElem.h:56
bool _have_weights
Definition: XFEMCutElem.h:53
virtual void computePhysicalVolumeFraction()=0
Computes the volume fraction of the element fragment.
virtual void computeMomentFittingWeights()=0
bool isPointPhysical(const Point &p) const
Definition: XFEMCutElem.C:191

◆ getCrackTipOriginAndDirection()

void XFEMCutElem2D::getCrackTipOriginAndDirection ( unsigned  tip_id,
Point &  origin,
Point &  direction 
) const
virtual

Implements XFEMCutElem.

Definition at line 180 of file XFEMCutElem2D.C.

183 {
184  // TODO: two cut plane case is not working
185  std::vector<EFANode *> cut_line_nodes;
186  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
187  {
189  {
190  std::vector<EFANode *> node_line(2, NULL);
191  node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
192  node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
193  if (node_line[1]->id() == tip_id)
194  {
195  cut_line_nodes.push_back(node_line[0]);
196  cut_line_nodes.push_back(node_line[1]);
197  }
198  else if (node_line[0]->id() == tip_id)
199  {
200  node_line[1] = node_line[0];
201  node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
202  cut_line_nodes.push_back(node_line[0]);
203  cut_line_nodes.push_back(node_line[1]);
204  }
205  }
206  }
207  if (cut_line_nodes.size() == 0)
208  mooseError("no cut line found in this element");
209 
210  Point cut_line_p1 = getNodeCoordinates(cut_line_nodes[0]);
211  Point cut_line_p2 = getNodeCoordinates(cut_line_nodes[1]);
212  Point cut_line = cut_line_p2 - cut_line_p1;
213  Real len = std::sqrt(cut_line.norm_sq());
214  cut_line /= len;
215  origin = cut_line_p2;
216  direction = Point(cut_line(0), cut_line(1), 0.0);
217 }
unsigned int numEdges() const
bool isEdgeInterior(unsigned int edge_id) const
EFAFragment2D * getFragment(unsigned int frag_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:36
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:179
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ getCutPlaneNormal()

Point XFEMCutElem2D::getCutPlaneNormal ( unsigned int  plane_id,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 151 of file XFEMCutElem2D.C.

Referenced by getIntersectionInfo().

152 {
153  Point normal(0.0, 0.0, 0.0);
154  std::vector<std::vector<EFANode *>> cut_line_nodes;
155  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
156  {
158  {
159  std::vector<EFANode *> node_line(2, NULL);
160  node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
161  node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
162  cut_line_nodes.push_back(node_line);
163  }
164  }
165  if (cut_line_nodes.size() == 0)
166  mooseError("no cut line found in this element");
167  if (plane_id < cut_line_nodes.size()) // valid plane_id
168  {
169  Point cut_line_p1 = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
170  Point cut_line_p2 = getNodeCoordinates(cut_line_nodes[plane_id][1], displaced_mesh);
171  Point cut_line = cut_line_p2 - cut_line_p1;
172  Real len = std::sqrt(cut_line.norm_sq());
173  cut_line /= len;
174  normal = Point(cut_line(1), -cut_line(0), 0.0);
175  }
176  return normal;
177 }
unsigned int numEdges() const
bool isEdgeInterior(unsigned int edge_id) const
EFAFragment2D * getFragment(unsigned int frag_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:36
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:179
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ getCutPlaneOrigin()

Point XFEMCutElem2D::getCutPlaneOrigin ( unsigned int  plane_id,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 129 of file XFEMCutElem2D.C.

130 {
131  Point orig(0.0, 0.0, 0.0);
132  std::vector<std::vector<EFANode *>> cut_line_nodes;
133  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
134  {
136  {
137  std::vector<EFANode *> node_line(2, NULL);
138  node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
139  node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
140  cut_line_nodes.push_back(node_line);
141  }
142  }
143  if (cut_line_nodes.size() == 0)
144  mooseError("no cut line found in this element");
145  if (plane_id < cut_line_nodes.size()) // valid plane_id
146  orig = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
147  return orig;
148 }
unsigned int numEdges() const
bool isEdgeInterior(unsigned int edge_id) const
EFAFragment2D * getFragment(unsigned int frag_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:36
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:179
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ getEFAElement()

const EFAElement * XFEMCutElem2D::getEFAElement ( ) const
virtual

Implements XFEMCutElem.

Definition at line 236 of file XFEMCutElem2D.C.

237 {
238  return &_efa_elem2d;
239 }
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ getFaceWeightMultipliers()

void XFEMCutElem::getFaceWeightMultipliers ( MooseArray< Real > &  face_weights,
QBase *  qrule,
Xfem::XFEM_QRULE  xfem_qrule,
const MooseArray< Point > &  q_points,
unsigned int  side 
)
inherited

Definition at line 75 of file XFEMCutElem.C.

Referenced by XFEM::getXFEMFaceWeights().

80 {
81  if (!_have_face_weights[side])
82  computeXFEMFaceWeights(qrule, xfem_qrule, q_points, side);
83 
84  face_weights.resize(_new_face_weights[side].size());
85  for (unsigned int qp = 0; qp < _new_face_weights[side].size(); ++qp)
86  face_weights[qp] = _new_face_weights[side][qp];
87 }
std::vector< std::vector< Real > > _new_face_weights
face quadrature weights from surface area fraction
Definition: XFEMCutElem.h:58
std::vector< bool > _have_face_weights
Definition: XFEMCutElem.h:54
void computeXFEMFaceWeights(QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points, unsigned int side)
Computes face integration weights for the cut element side.
Definition: XFEMCutElem.C:90

◆ getFragmentFaces()

void XFEMCutElem2D::getFragmentFaces ( std::vector< std::vector< Point >> &  frag_faces,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 220 of file XFEMCutElem2D.C.

222 {
223  frag_faces.clear();
224  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
225  {
226  std::vector<Point> edge_points(2, Point(0.0, 0.0, 0.0));
227  edge_points[0] =
228  getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(0), displaced_mesh);
229  edge_points[1] =
230  getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(1), displaced_mesh);
231  frag_faces.push_back(edge_points);
232  }
233 }
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:36
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:179
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ getIntersectionInfo()

void XFEMCutElem2D::getIntersectionInfo ( unsigned int  plane_id,
Point &  normal,
std::vector< Point > &  intersectionPoints,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 441 of file XFEMCutElem2D.C.

445 {
446  intersectionPoints.resize(2); // 2D: the intersection line is straight and can be represented by
447  // two ending points.(may have issues with double cuts case!)
448 
449  std::vector<std::vector<EFANode *>> cut_line_nodes;
450  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
451  {
453  {
454  std::vector<EFANode *> node_line(2, NULL);
455  node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
456  node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
457  cut_line_nodes.push_back(node_line);
458  }
459  }
460  if (cut_line_nodes.size() == 0)
461  mooseError("No cut line found in this element");
462 
463  if (plane_id < cut_line_nodes.size()) // valid plane_id
464  {
465  intersectionPoints[0] = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
466  intersectionPoints[1] = getNodeCoordinates(cut_line_nodes[plane_id][1], displaced_mesh);
467  }
468 
469  normal = getCutPlaneNormal(plane_id, displaced_mesh);
470 }
unsigned int numEdges() const
bool isEdgeInterior(unsigned int edge_id) const
EFAFragment2D * getFragment(unsigned int frag_id) const
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:36
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:179
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ getMomentFittingWeight()

Real XFEMCutElem::getMomentFittingWeight ( unsigned int  i_qp) const
inherited

◆ getNodeCoordinates()

Point XFEMCutElem2D::getNodeCoordinates ( EFANode node,
MeshBase *  displaced_mesh = NULL 
) const
privatevirtual

Implements XFEMCutElem.

Definition at line 36 of file XFEMCutElem2D.C.

Referenced by computePhysicalFaceAreaFraction(), computePhysicalVolumeFraction(), getCrackTipOriginAndDirection(), getCutPlaneNormal(), getCutPlaneOrigin(), getFragmentFaces(), getIntersectionInfo(), and getPhysicalQuadraturePoints().

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_elem2d.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), 0.0);
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 }
std::vector< Node * > _nodes
Definition: XFEMCutElem.h:46
virtual void getMasterInfo(EFANode *node, std::vector< EFANode *> &master_nodes, std::vector< double > &master_weights) const
Definition: EFAElement2D.C:300
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ getPhysicalFaceAreaFraction()

Real XFEMCutElem::getPhysicalFaceAreaFraction ( unsigned int  side) const
inherited

Returns the surface area fraction of the element side.

Parameters
sideThe side of the element

Definition at line 55 of file XFEMCutElem.C.

Referenced by XFEMCutElem::computeXFEMFaceWeights().

56 {
57  return _physical_areafrac[side];
58 }
std::vector< Real > _physical_areafrac
Definition: XFEMCutElem.h:52

◆ getPhysicalQuadraturePoints()

void XFEMCutElem2D::getPhysicalQuadraturePoints ( std::vector< std::vector< Real >> &  tsg)
private

Definition at line 252 of file XFEMCutElem2D.C.

Referenced by computeMomentFittingWeights().

253 {
254  // Get the coords for parial element nodes
256  unsigned int nnd_pe = frag->numEdges();
257  std::vector<Point> frag_points(nnd_pe, Point(0.0, 0.0, 0.0)); // nodal coord of partial elem
258  Real jac = 0.0;
259 
260  for (unsigned int j = 0; j < nnd_pe; ++j)
261  frag_points[j] = getNodeCoordinates(frag->getEdge(j)->getNode(0));
262 
263  // Get centroid coords for partial elements
264  Point xcrd(0.0, 0.0, 0.0);
265  for (unsigned int j = 0; j < nnd_pe; ++j)
266  xcrd += frag_points[j];
267  xcrd *= (1.0 / nnd_pe);
268 
269  // Get tsg - the physical coords of Gaussian Q-points for partial elements
270  if ((nnd_pe == 3) || (nnd_pe == 4)) // partial element is a triangle or quad
271  {
272  std::vector<std::vector<Real>> sg2;
273  std::vector<std::vector<Real>> shape(nnd_pe, std::vector<Real>(3, 0.0));
274  Xfem::stdQuadr2D(nnd_pe, 2, sg2);
275  for (unsigned int l = 0; l < sg2.size(); ++l)
276  {
277  Xfem::shapeFunc2D(nnd_pe, sg2[l], frag_points, shape, jac, true); // Get shape
278  std::vector<Real> tsg_line(3, 0.0);
279  for (unsigned int k = 0; k < nnd_pe; ++k)
280  {
281  tsg_line[0] += shape[k][2] * frag_points[k](0);
282  tsg_line[1] += shape[k][2] * frag_points[k](1);
283  }
284  if (nnd_pe == 3) // tri partial elem
285  tsg_line[2] = sg2[l][3] * jac; // total weights
286  else // quad partial elem
287  tsg_line[2] = sg2[l][2] * jac; // total weights
288  tsg.push_back(tsg_line);
289  }
290  }
291  else if (nnd_pe >= 5) // partial element is a polygon
292  {
293  for (unsigned int j = 0; j < nnd_pe; ++j) // loop all sub-tris
294  {
295  std::vector<std::vector<Real>> shape(3, std::vector<Real>(3, 0.0));
296  std::vector<Point> subtri_points(3, Point(0.0, 0.0, 0.0)); // sub-tri nodal coords
297 
298  int jplus1 = j < nnd_pe - 1 ? j + 1 : 0;
299  subtri_points[0] = xcrd;
300  subtri_points[1] = frag_points[j];
301  subtri_points[2] = frag_points[jplus1];
302 
303  std::vector<std::vector<Real>> sg2;
304  Xfem::stdQuadr2D(3, 2, sg2); // get sg2
305  for (unsigned int l = 0; l < sg2.size(); ++l) // loop all int pts on a sub-tri
306  {
307  Xfem::shapeFunc2D(3, sg2[l], subtri_points, shape, jac, true); // Get shape
308  std::vector<Real> tsg_line(3, 0.0);
309  for (unsigned int k = 0; k < 3; ++k) // loop sub-tri nodes
310  {
311  tsg_line[0] += shape[k][2] * subtri_points[k](0);
312  tsg_line[1] += shape[k][2] * subtri_points[k](1);
313  }
314  tsg_line[2] = sg2[l][3] * jac;
315  tsg.push_back(tsg_line);
316  }
317  }
318  }
319  else
320  mooseError("Invalid partial element!");
321 }
void shapeFunc2D(unsigned int nen, std::vector< Real > &ss, std::vector< Point > &xl, std::vector< std::vector< Real >> &shp, Real &xsj, bool natl_flg)
Definition: XFEMFuncs.C:269
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
void stdQuadr2D(unsigned int nen, unsigned int iord, std::vector< std::vector< Real >> &sg2)
Definition: XFEMFuncs.C:96
EFAEdge * getEdge(unsigned int edge_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:36
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:179
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ getPhysicalVolumeFraction()

Real XFEMCutElem::getPhysicalVolumeFraction ( ) const
inherited

Returns the volume fraction of the element fragment.

Definition at line 49 of file XFEMCutElem.C.

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

50 {
51  return _physical_volfrac;
52 }
Real _physical_volfrac
Definition: XFEMCutElem.h:51

◆ getWeightMultipliers()

void XFEMCutElem::getWeightMultipliers ( MooseArray< Real > &  weights,
QBase *  qrule,
Xfem::XFEM_QRULE  xfem_qrule,
const MooseArray< Point > &  q_points 
)
inherited

Definition at line 61 of file XFEMCutElem.C.

Referenced by XFEM::getXFEMWeights().

65 {
66  if (!_have_weights)
67  computeXFEMWeights(qrule, xfem_qrule, q_points);
68 
69  weights.resize(_new_weights.size());
70  for (unsigned int qp = 0; qp < _new_weights.size(); ++qp)
71  weights[qp] = _new_weights[qp];
72 }
std::vector< Real > _new_weights
quadrature weights from volume fraction and moment fitting
Definition: XFEMCutElem.h:56
bool _have_weights
Definition: XFEMCutElem.h:53
void computeXFEMWeights(QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
Computes integration weights for the cut element.
Definition: XFEMCutElem.C:107

◆ isPointPhysical()

bool XFEMCutElem::isPointPhysical ( const Point &  p) const
inherited

Definition at line 191 of file XFEMCutElem.C.

Referenced by XFEMCutElem::computeXFEMWeights(), XFEM::healMesh(), and XFEM::isPointInsidePhysicalDomain().

192 {
193  // determine whether the point is inside the physical domain of a partial element
194  bool physical_flag = true;
195  unsigned int n_cut_planes = numCutPlanes();
196  for (unsigned int plane_id = 0; plane_id < n_cut_planes; ++plane_id)
197  {
198  Point origin = getCutPlaneOrigin(plane_id);
199  Point normal = getCutPlaneNormal(plane_id);
200  Point origin2qp = p - origin;
201  if (origin2qp * normal > 0.0)
202  {
203  physical_flag = false; // Point outside pysical domain
204  break;
205  }
206  }
207  return physical_flag;
208 }
virtual Point getCutPlaneOrigin(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const =0
virtual unsigned int numCutPlanes() const =0
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const =0

◆ numCutPlanes()

unsigned int XFEMCutElem2D::numCutPlanes ( ) const
virtual

Implements XFEMCutElem.

Definition at line 242 of file XFEMCutElem2D.C.

243 {
244  unsigned int counter = 0;
245  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
247  counter += 1;
248  return counter;
249 }
unsigned int numEdges() const
bool isEdgeInterior(unsigned int edge_id) const
EFAFragment2D * getFragment(unsigned int frag_id) const
static unsigned int counter
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:42

◆ setQuadraturePointsAndWeights()

void XFEMCutElem::setQuadraturePointsAndWeights ( const std::vector< Point > &  qp_points,
const std::vector< Real > &  qp_weights 
)
inherited

◆ solveMomentFitting()

void XFEMCutElem2D::solveMomentFitting ( unsigned int  nen,
unsigned int  nqp,
std::vector< Point > &  elem_nodes,
std::vector< std::vector< Real >> &  tsg,
std::vector< std::vector< Real >> &  wsg 
)
private

Definition at line 324 of file XFEMCutElem2D.C.

Referenced by computeMomentFittingWeights().

329 {
330  // Get physical coords for the new six-point rule
331  std::vector<std::vector<Real>> shape(nen, std::vector<Real>(3, 0.0));
332  std::vector<std::vector<Real>> wss;
333 
334  if (nen == 4)
335  {
336  wss.resize(_qp_points.size());
337  for (unsigned int i = 0; i < _qp_points.size(); i++)
338  {
339  wss[i].resize(3);
340  wss[i][0] = _qp_points[i](0);
341  wss[i][1] = _qp_points[i](1);
342  wss[i][2] = _qp_weights[i];
343  }
344  }
345  else if (nen == 3)
346  {
347  wss.resize(_qp_points.size());
348  for (unsigned int i = 0; i < _qp_points.size(); i++)
349  {
350  wss[i].resize(4);
351  wss[i][0] = _qp_points[i](0);
352  wss[i][1] = _qp_points[i](1);
353  wss[i][2] = 1.0 - _qp_points[i](0) - _qp_points[i](1);
354  wss[i][3] = _qp_weights[i];
355  }
356  }
357  else
358  mooseError("Invalid element");
359 
360  wsg.resize(wss.size());
361  for (unsigned int i = 0; i < wsg.size(); ++i)
362  wsg[i].resize(3, 0.0);
363  Real jac = 0.0;
364  std::vector<Real> old_weights(wss.size(), 0.0);
365 
366  for (unsigned int l = 0; l < wsg.size(); ++l)
367  {
368  Xfem::shapeFunc2D(nen, wss[l], elem_nodes, shape, jac, true); // Get shape
369  if (nen == 4) // 2D quad elem
370  old_weights[l] = wss[l][2] * jac; // weights for total element
371  else if (nen == 3) // 2D triangle elem
372  old_weights[l] = wss[l][3] * jac;
373  else
374  mooseError("Invalid element!");
375  for (unsigned int k = 0; k < nen; ++k) // physical coords of Q-pts
376  {
377  wsg[l][0] += shape[k][2] * elem_nodes[k](0);
378  wsg[l][1] += shape[k][2] * elem_nodes[k](1);
379  }
380  }
381 
382  // Compute weights via moment fitting
383  Real * A;
384  A = new Real[wsg.size() * wsg.size()];
385  unsigned ind = 0;
386  for (unsigned int i = 0; i < wsg.size(); ++i)
387  {
388  A[ind] = 1.0; // const
389  if (nqp > 1)
390  A[1 + ind] = wsg[i][0]; // x
391  if (nqp > 2)
392  A[2 + ind] = wsg[i][1]; // y
393  if (nqp > 3)
394  A[3 + ind] = wsg[i][0] * wsg[i][1]; // x*y
395  if (nqp > 4)
396  A[4 + ind] = wsg[i][0] * wsg[i][0]; // x^2
397  if (nqp > 5)
398  A[5 + ind] = wsg[i][1] * wsg[i][1]; // y^2
399  if (nqp > 6)
400  mooseError("Q-points of more than 6 are not allowed now!");
401  ind = ind + nqp;
402  }
403 
404  Real * b;
405  b = new Real[wsg.size()];
406  for (unsigned int i = 0; i < wsg.size(); ++i)
407  b[i] = 0.0;
408  for (unsigned int i = 0; i < tsg.size(); ++i)
409  {
410  b[0] += tsg[i][2];
411  if (nqp > 1)
412  b[1] += tsg[i][2] * tsg[i][0];
413  if (nqp > 2)
414  b[2] += tsg[i][2] * tsg[i][1];
415  if (nqp > 3)
416  b[3] += tsg[i][2] * tsg[i][0] * tsg[i][1];
417  if (nqp > 4)
418  b[4] += tsg[i][2] * tsg[i][0] * tsg[i][0];
419  if (nqp > 5)
420  b[5] += tsg[i][2] * tsg[i][1] * tsg[i][1];
421  if (nqp > 6)
422  mooseError("Q-points of more than 6 are not allowed now!");
423  }
424 
425  int nrhs = 1;
426  int info;
427  int n = wsg.size();
428  std::vector<int> ipiv(n);
429 
430  LAPACKgesv_(&n, &nrhs, A, &n, &ipiv[0], b, &n, &info);
431 
432  for (unsigned int i = 0; i < wsg.size(); ++i)
433  wsg[i][2] = b[i] / old_weights[i]; // get the multiplier
434 
435  // delete arrays
436  delete[] A;
437  delete[] b;
438 }
void shapeFunc2D(unsigned int nen, std::vector< Real > &ss, std::vector< Point > &xl, std::vector< std::vector< Real >> &shp, Real &xsj, bool natl_flg)
Definition: XFEMFuncs.C:269
std::vector< Point > _qp_points
Definition: XFEMCutElem.h:47
std::vector< Real > _qp_weights
Definition: XFEMCutElem.h:48

Member Data Documentation

◆ _efa_elem2d

EFAElement2D XFEMCutElem2D::_efa_elem2d
private

◆ _elem_side_area

std::vector<Real> XFEMCutElem::_elem_side_area
protectedinherited

◆ _elem_volume

Real XFEMCutElem::_elem_volume
protectedinherited

◆ _have_face_weights

std::vector<bool> XFEMCutElem::_have_face_weights
protectedinherited

◆ _have_weights

bool XFEMCutElem::_have_weights
protectedinherited

◆ _n_nodes

unsigned int XFEMCutElem::_n_nodes
protectedinherited

Definition at line 43 of file XFEMCutElem.h.

Referenced by computeMomentFittingWeights(), and XFEMCutElem::XFEMCutElem().

◆ _n_qpoints

unsigned int XFEMCutElem::_n_qpoints
protectedinherited

◆ _n_sides

unsigned int XFEMCutElem::_n_sides
protectedinherited

Definition at line 45 of file XFEMCutElem.h.

Referenced by XFEMCutElem::XFEMCutElem().

◆ _new_face_weights

std::vector<std::vector<Real> > XFEMCutElem::_new_face_weights
protectedinherited

face quadrature weights from surface area fraction

Definition at line 58 of file XFEMCutElem.h.

Referenced by XFEMCutElem::computeXFEMFaceWeights(), and XFEMCutElem::getFaceWeightMultipliers().

◆ _new_weights

std::vector<Real> XFEMCutElem::_new_weights
protectedinherited

quadrature weights from volume fraction and moment fitting

Definition at line 56 of file XFEMCutElem.h.

Referenced by computeMomentFittingWeights(), XFEMCutElem3D::computeMomentFittingWeights(), XFEMCutElem::computeXFEMWeights(), and XFEMCutElem::getWeightMultipliers().

◆ _nodes

std::vector<Node *> XFEMCutElem::_nodes
protectedinherited

◆ _physical_areafrac

std::vector<Real> XFEMCutElem::_physical_areafrac
protectedinherited

◆ _physical_volfrac

Real XFEMCutElem::_physical_volfrac
protectedinherited

◆ _qp_points

std::vector<Point> XFEMCutElem::_qp_points
protectedinherited

Definition at line 47 of file XFEMCutElem.h.

Referenced by XFEMCutElem::computeXFEMWeights(), and solveMomentFitting().

◆ _qp_weights

std::vector<Real> XFEMCutElem::_qp_weights
protectedinherited

Definition at line 48 of file XFEMCutElem.h.

Referenced by XFEMCutElem::computeXFEMWeights(), and solveMomentFitting().


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