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)
 
 ~XFEMCutElem2D ()
 
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
 
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 23 of file XFEMCutElem2D.h.

Constructor & Destructor Documentation

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

Definition at line 22 of file XFEMCutElem2D.C.

25  : XFEMCutElem(elem, n_qpoints), _efa_elem2d(CEMelem, true)
26 {
28 }
virtual void computePhysicalVolumeFraction()
Definition: XFEMCutElem2D.C:60
XFEMCutElem(Elem *elem, unsigned int n_qpoints)
Definition: XFEMCutElem.C:20
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
XFEMCutElem2D::~XFEMCutElem2D ( )

Definition at line 30 of file XFEMCutElem2D.C.

30 {}

Member Function Documentation

void XFEMCutElem2D::computeMomentFittingWeights ( )
virtual

Implements XFEMCutElem.

Definition at line 77 of file XFEMCutElem2D.C.

78 {
79  // Purpose: calculate new weights via moment-fitting method
80  std::vector<Point> elem_nodes(_n_nodes, Point(0.0, 0.0, 0.0));
81  std::vector<std::vector<Real>> wsg;
82 
83  for (unsigned int i = 0; i < _n_nodes; ++i)
84  elem_nodes[i] = (*_nodes[i]);
85 
86  if (_efa_elem2d.isPartial() && _n_qpoints <= 6) // ONLY work for <= 6 q_points
87  {
88  std::vector<std::vector<Real>> tsg;
89  getPhysicalQuadraturePoints(tsg); // get tsg - QPs within partial element
91  _n_nodes, _n_qpoints, elem_nodes, tsg, wsg); // get wsg - QPs from moment-fitting
92  _new_weights.resize(wsg.size(), 1.0);
93  for (unsigned int i = 0; i < wsg.size(); ++i)
94  _new_weights[i] = wsg[i][2]; // weight multiplier
95  }
96  else
98 }
virtual bool isPartial() const
Definition: EFAElement2D.C:206
unsigned int _n_qpoints
Definition: XFEMCutElem.h:36
Real _physical_volfrac
Definition: XFEMCutElem.h:41
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
Definition: XFEMCutElem.h:43
void getPhysicalQuadraturePoints(std::vector< std::vector< Real >> &tsg)
unsigned int _n_nodes
Definition: XFEMCutElem.h:35
std::vector< Node * > _nodes
Definition: XFEMCutElem.h:37
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
void XFEMCutElem2D::computePhysicalVolumeFraction ( )
virtual

Implements XFEMCutElem.

Definition at line 60 of file XFEMCutElem2D.C.

Referenced by XFEMCutElem2D().

61 {
62  Real frag_vol = 0.0;
63 
64  // Calculate area of entire element and fragment using the formula:
65  // A = 1/2 sum_{i=0}^{n-1} (x_i y_{i+1} - x_{i+1} y{i})
66 
67  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
68  {
69  Point edge_p1 = getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(0));
70  Point edge_p2 = getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(1));
71  frag_vol += 0.5 * (edge_p1(0) - edge_p2(0)) * (edge_p1(1) + edge_p2(1));
72  }
73  _physical_volfrac = frag_vol / _elem_volume;
74 }
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
Real _physical_volfrac
Definition: XFEMCutElem.h:41
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:33
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
Real _elem_volume
Definition: XFEMCutElem.h:40
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
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 XFEMCutElem2D::getCrackTipOriginAndDirection ( unsigned  tip_id,
Point &  origin,
Point &  direction 
) const
virtual

Implements XFEMCutElem.

Definition at line 152 of file XFEMCutElem2D.C.

155 {
156  // TODO: two cut plane case is not working
157  std::vector<EFANode *> cut_line_nodes;
158  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
159  {
161  {
162  std::vector<EFANode *> node_line(2, NULL);
163  node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
164  node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
165  if (node_line[1]->id() == tip_id)
166  {
167  cut_line_nodes.push_back(node_line[0]);
168  cut_line_nodes.push_back(node_line[1]);
169  }
170  else if (node_line[0]->id() == tip_id)
171  {
172  node_line[1] = node_line[0];
173  node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
174  cut_line_nodes.push_back(node_line[0]);
175  cut_line_nodes.push_back(node_line[1]);
176  }
177  }
178  }
179  if (cut_line_nodes.size() == 0)
180  mooseError("no cut line found in this element");
181 
182  Point cut_line_p1 = getNodeCoordinates(cut_line_nodes[0]);
183  Point cut_line_p2 = getNodeCoordinates(cut_line_nodes[1]);
184  Point cut_line = cut_line_p2 - cut_line_p1;
185  Real len = std::sqrt(cut_line.norm_sq());
186  cut_line /= len;
187  origin = cut_line_p2;
188  direction = Point(cut_line(0), cut_line(1), 0.0);
189 }
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:33
bool isEdgeInterior(unsigned int edge_id) const
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
Point XFEMCutElem2D::getCutPlaneNormal ( unsigned int  plane_id,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 123 of file XFEMCutElem2D.C.

Referenced by getIntersectionInfo().

124 {
125  Point normal(0.0, 0.0, 0.0);
126  std::vector<std::vector<EFANode *>> cut_line_nodes;
127  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
128  {
130  {
131  std::vector<EFANode *> node_line(2, NULL);
132  node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
133  node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
134  cut_line_nodes.push_back(node_line);
135  }
136  }
137  if (cut_line_nodes.size() == 0)
138  mooseError("no cut line found in this element");
139  if (plane_id < cut_line_nodes.size()) // valid plane_id
140  {
141  Point cut_line_p1 = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
142  Point cut_line_p2 = getNodeCoordinates(cut_line_nodes[plane_id][1], displaced_mesh);
143  Point cut_line = cut_line_p2 - cut_line_p1;
144  Real len = std::sqrt(cut_line.norm_sq());
145  cut_line /= len;
146  normal = Point(cut_line(1), -cut_line(0), 0.0);
147  }
148  return normal;
149 }
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:33
bool isEdgeInterior(unsigned int edge_id) const
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
Point XFEMCutElem2D::getCutPlaneOrigin ( unsigned int  plane_id,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 101 of file XFEMCutElem2D.C.

102 {
103  Point orig(0.0, 0.0, 0.0);
104  std::vector<std::vector<EFANode *>> cut_line_nodes;
105  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
106  {
108  {
109  std::vector<EFANode *> node_line(2, NULL);
110  node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
111  node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
112  cut_line_nodes.push_back(node_line);
113  }
114  }
115  if (cut_line_nodes.size() == 0)
116  mooseError("no cut line found in this element");
117  if (plane_id < cut_line_nodes.size()) // valid plane_id
118  orig = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
119  return orig;
120 }
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:33
bool isEdgeInterior(unsigned int edge_id) const
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
const EFAElement * XFEMCutElem2D::getEFAElement ( ) const
virtual

Implements XFEMCutElem.

Definition at line 208 of file XFEMCutElem2D.C.

209 {
210  return &_efa_elem2d;
211 }
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
void XFEMCutElem2D::getFragmentFaces ( std::vector< std::vector< Point >> &  frag_faces,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 192 of file XFEMCutElem2D.C.

194 {
195  frag_faces.clear();
196  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
197  {
198  std::vector<Point> edge_points(2, Point(0.0, 0.0, 0.0));
199  edge_points[0] =
200  getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(0), displaced_mesh);
201  edge_points[1] =
202  getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(1), displaced_mesh);
203  frag_faces.push_back(edge_points);
204  }
205 }
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:33
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
void XFEMCutElem2D::getIntersectionInfo ( unsigned int  plane_id,
Point &  normal,
std::vector< Point > &  intersectionPoints,
MeshBase *  displaced_mesh = NULL 
) const
virtual

Implements XFEMCutElem.

Definition at line 413 of file XFEMCutElem2D.C.

417 {
418  intersectionPoints.resize(2); // 2D: the intersection line is straight and can be represented by
419  // two ending points.(may have issues with double cuts case!)
420 
421  std::vector<std::vector<EFANode *>> cut_line_nodes;
422  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
423  {
425  {
426  std::vector<EFANode *> node_line(2, NULL);
427  node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
428  node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
429  cut_line_nodes.push_back(node_line);
430  }
431  }
432  if (cut_line_nodes.size() == 0)
433  mooseError("No cut line found in this element");
434 
435  if (plane_id < cut_line_nodes.size()) // valid plane_id
436  {
437  intersectionPoints[0] = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
438  intersectionPoints[1] = getNodeCoordinates(cut_line_nodes[plane_id][1], displaced_mesh);
439  }
440 
441  normal = getCutPlaneNormal(plane_id, displaced_mesh);
442 }
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:33
bool isEdgeInterior(unsigned int edge_id) const
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
Real XFEMCutElem::getMomentFittingWeight ( unsigned int  i_qp) const
inherited
Point XFEMCutElem2D::getNodeCoordinates ( EFANode node,
MeshBase *  displaced_mesh = NULL 
) const
privatevirtual

Implements XFEMCutElem.

Definition at line 33 of file XFEMCutElem2D.C.

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

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_elem2d.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), 0.0);
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 }
virtual void getMasterInfo(EFANode *node, std::vector< EFANode * > &master_nodes, std::vector< double > &master_weights) const
Definition: EFAElement2D.C:298
std::vector< Node * > _nodes
Definition: XFEMCutElem.h:37
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
void XFEMCutElem2D::getPhysicalQuadraturePoints ( std::vector< std::vector< Real >> &  tsg)
private

Definition at line 224 of file XFEMCutElem2D.C.

Referenced by computeMomentFittingWeights().

225 {
226  // Get the coords for parial element nodes
228  unsigned int nnd_pe = frag->numEdges();
229  std::vector<Point> frag_points(nnd_pe, Point(0.0, 0.0, 0.0)); // nodal coord of partial elem
230  Real jac = 0.0;
231 
232  for (unsigned int j = 0; j < nnd_pe; ++j)
233  frag_points[j] = getNodeCoordinates(frag->getEdge(j)->getNode(0));
234 
235  // Get centroid coords for partial elements
236  Point xcrd(0.0, 0.0, 0.0);
237  for (unsigned int j = 0; j < nnd_pe; ++j)
238  xcrd += frag_points[j];
239  xcrd *= (1.0 / nnd_pe);
240 
241  // Get tsg - the physical coords of Gaussian Q-points for partial elements
242  if ((nnd_pe == 3) || (nnd_pe == 4)) // partial element is a triangle or quad
243  {
244  std::vector<std::vector<Real>> sg2;
245  std::vector<std::vector<Real>> shape(nnd_pe, std::vector<Real>(3, 0.0));
246  Xfem::stdQuadr2D(nnd_pe, 2, sg2);
247  for (unsigned int l = 0; l < sg2.size(); ++l)
248  {
249  Xfem::shapeFunc2D(nnd_pe, sg2[l], frag_points, shape, jac, true); // Get shape
250  std::vector<Real> tsg_line(3, 0.0);
251  for (unsigned int k = 0; k < nnd_pe; ++k)
252  {
253  tsg_line[0] += shape[k][2] * frag_points[k](0);
254  tsg_line[1] += shape[k][2] * frag_points[k](1);
255  }
256  if (nnd_pe == 3) // tri partial elem
257  tsg_line[2] = sg2[l][3] * jac; // total weights
258  else // quad partial elem
259  tsg_line[2] = sg2[l][2] * jac; // total weights
260  tsg.push_back(tsg_line);
261  }
262  }
263  else if (nnd_pe >= 5) // partial element is a polygon
264  {
265  for (unsigned int j = 0; j < nnd_pe; ++j) // loop all sub-tris
266  {
267  std::vector<std::vector<Real>> shape(3, std::vector<Real>(3, 0.0));
268  std::vector<Point> subtri_points(3, Point(0.0, 0.0, 0.0)); // sub-tri nodal coords
269 
270  int jplus1 = j < nnd_pe - 1 ? j + 1 : 0;
271  subtri_points[0] = xcrd;
272  subtri_points[1] = frag_points[j];
273  subtri_points[2] = frag_points[jplus1];
274 
275  std::vector<std::vector<Real>> sg2;
276  Xfem::stdQuadr2D(3, 2, sg2); // get sg2
277  for (unsigned int l = 0; l < sg2.size(); ++l) // loop all int pts on a sub-tri
278  {
279  Xfem::shapeFunc2D(3, sg2[l], subtri_points, shape, jac, true); // Get shape
280  std::vector<Real> tsg_line(3, 0.0);
281  for (unsigned int k = 0; k < 3; ++k) // loop sub-tri nodes
282  {
283  tsg_line[0] += shape[k][2] * subtri_points[k](0);
284  tsg_line[1] += shape[k][2] * subtri_points[k](1);
285  }
286  tsg_line[2] = sg2[l][3] * jac;
287  tsg.push_back(tsg_line);
288  }
289  }
290  }
291  else
292  mooseError("Invalid partial element!");
293 }
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
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:274
virtual Point getNodeCoordinates(EFANode *node, MeshBase *displaced_mesh=NULL) const
Definition: XFEMCutElem2D.C:33
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
void stdQuadr2D(unsigned int nen, unsigned int iord, std::vector< std::vector< Real >> &sg2)
Definition: XFEMFuncs.C:101
EFAEdge * getEdge(unsigned int edge_id) const
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
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 XFEMCutElem2D::numCutPlanes ( ) const
virtual

Implements XFEMCutElem.

Definition at line 214 of file XFEMCutElem2D.C.

215 {
216  unsigned int counter = 0;
217  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
219  counter += 1;
220  return counter;
221 }
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
bool isEdgeInterior(unsigned int edge_id) const
static unsigned int counter
EFAElement2D _efa_elem2d
Definition: XFEMCutElem2D.h:30
void XFEMCutElem::setQuadraturePointsAndWeights ( const std::vector< Point > &  qp_points,
const std::vector< Real > &  qp_weights 
)
inherited
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 296 of file XFEMCutElem2D.C.

Referenced by computeMomentFittingWeights().

301 {
302  // Get physical coords for the new six-point rule
303  std::vector<std::vector<Real>> shape(nen, std::vector<Real>(3, 0.0));
304  std::vector<std::vector<Real>> wss;
305 
306  if (nen == 4)
307  {
308  wss.resize(_qp_points.size());
309  for (unsigned int i = 0; i < _qp_points.size(); i++)
310  {
311  wss[i].resize(3);
312  wss[i][0] = _qp_points[i](0);
313  wss[i][1] = _qp_points[i](1);
314  wss[i][2] = _qp_weights[i];
315  }
316  }
317  else if (nen == 3)
318  {
319  wss.resize(_qp_points.size());
320  for (unsigned int i = 0; i < _qp_points.size(); i++)
321  {
322  wss[i].resize(4);
323  wss[i][0] = _qp_points[i](0);
324  wss[i][1] = _qp_points[i](1);
325  wss[i][2] = 1.0 - _qp_points[i](0) - _qp_points[i](1);
326  wss[i][3] = _qp_weights[i];
327  }
328  }
329  else
330  mooseError("Invalid element");
331 
332  wsg.resize(wss.size());
333  for (unsigned int i = 0; i < wsg.size(); ++i)
334  wsg[i].resize(3, 0.0);
335  Real jac = 0.0;
336  std::vector<Real> old_weights(wss.size(), 0.0);
337 
338  for (unsigned int l = 0; l < wsg.size(); ++l)
339  {
340  Xfem::shapeFunc2D(nen, wss[l], elem_nodes, shape, jac, true); // Get shape
341  if (nen == 4) // 2D quad elem
342  old_weights[l] = wss[l][2] * jac; // weights for total element
343  else if (nen == 3) // 2D triangle elem
344  old_weights[l] = wss[l][3] * jac;
345  else
346  mooseError("Invalid element!");
347  for (unsigned int k = 0; k < nen; ++k) // physical coords of Q-pts
348  {
349  wsg[l][0] += shape[k][2] * elem_nodes[k](0);
350  wsg[l][1] += shape[k][2] * elem_nodes[k](1);
351  }
352  }
353 
354  // Compute weights via moment fitting
355  Real * A;
356  A = new Real[wsg.size() * wsg.size()];
357  unsigned ind = 0;
358  for (unsigned int i = 0; i < wsg.size(); ++i)
359  {
360  A[ind] = 1.0; // const
361  if (nqp > 1)
362  A[1 + ind] = wsg[i][0]; // x
363  if (nqp > 2)
364  A[2 + ind] = wsg[i][1]; // y
365  if (nqp > 3)
366  A[3 + ind] = wsg[i][0] * wsg[i][1]; // x*y
367  if (nqp > 4)
368  A[4 + ind] = wsg[i][0] * wsg[i][0]; // x^2
369  if (nqp > 5)
370  A[5 + ind] = wsg[i][1] * wsg[i][1]; // y^2
371  if (nqp > 6)
372  mooseError("Q-points of more than 6 are not allowed now!");
373  ind = ind + nqp;
374  }
375 
376  Real * b;
377  b = new Real[wsg.size()];
378  for (unsigned int i = 0; i < wsg.size(); ++i)
379  b[i] = 0.0;
380  for (unsigned int i = 0; i < tsg.size(); ++i)
381  {
382  b[0] += tsg[i][2];
383  if (nqp > 1)
384  b[1] += tsg[i][2] * tsg[i][0];
385  if (nqp > 2)
386  b[2] += tsg[i][2] * tsg[i][1];
387  if (nqp > 3)
388  b[3] += tsg[i][2] * tsg[i][0] * tsg[i][1];
389  if (nqp > 4)
390  b[4] += tsg[i][2] * tsg[i][0] * tsg[i][0];
391  if (nqp > 5)
392  b[5] += tsg[i][2] * tsg[i][1] * tsg[i][1];
393  if (nqp > 6)
394  mooseError("Q-points of more than 6 are not allowed now!");
395  }
396 
397  int nrhs = 1;
398  int info;
399  int n = wsg.size();
400  std::vector<int> ipiv(n);
401 
402  LAPACKgesv_(&n, &nrhs, A, &n, &ipiv[0], b, &n, &info);
403 
404  for (unsigned int i = 0; i < wsg.size(); ++i)
405  wsg[i][2] = b[i] / old_weights[i]; // get the multiplier
406 
407  // delete arrays
408  delete[] A;
409  delete[] b;
410 }
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:274
std::vector< Point > _qp_points
Definition: XFEMCutElem.h:38
std::vector< Real > _qp_weights
Definition: XFEMCutElem.h:39

Member Data Documentation

EFAElement2D XFEMCutElem2D::_efa_elem2d
private
Real XFEMCutElem::_elem_volume
protectedinherited
bool XFEMCutElem::_have_weights
protectedinherited
unsigned int XFEMCutElem::_n_nodes
protectedinherited

Definition at line 35 of file XFEMCutElem.h.

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

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

Definition at line 38 of file XFEMCutElem.h.

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

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

Definition at line 39 of file XFEMCutElem.h.

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


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