#include <XFEMCutElem2D.h>

Inheritance diagram for XFEMCutElem2D:

Public Member Functions
	XFEMCutElem2D (Elem elem, const EFAElement2D const CEMelem, unsigned int n_qpoints, unsigned int n_sides)
	Constructor initializes XFEMCutElem2D object.

	~XFEMCutElem2D ()

virtual void	computePhysicalVolumeFraction ()
	Computes the volume fraction of the element fragment.

virtual void	computePhysicalFaceAreaFraction (unsigned int side)
	Computes the surface area fraction of the element side.

virtual void	computeMomentFittingWeights ()

virtual Point	getCutPlaneOrigin (unsigned int plane_id, MeshBase *displaced_mesh=nullptr) const

virtual Point	getCutPlaneNormal (unsigned int plane_id, MeshBase *displaced_mesh=nullptr) 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=nullptr) const

virtual const EFAElement *	getEFAElement () const

virtual unsigned int	numCutPlanes () const

virtual void	getIntersectionInfo (unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, MeshBase *displaced_mesh=nullptr) 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.

Real	getPhysicalFaceAreaFraction (unsigned int side) const
	Returns the surface area fraction of the element side.

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.

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.

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

std::vector< std::vector< Real > >	_new_face_weights
	face quadrature weights from surface area fraction

Private Member Functions
virtual Point	getNodeCoordinates (EFANode node, MeshBase displaced_mesh=nullptr) 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 24 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

elem	The element on which XFEMCutElem2D is built
CEMelem	The EFAFragment2D object that belongs to XFEMCutElem2D
n_qpoints	The number of quadrature points
n_sides	The number of sides which the element has

Definition at line 24 of file XFEMCutElem2D.C.

  : XFEMCutElem(elem, n_qpoints, n_sides), _efa_elem2d(CEMelem, true)
{
  computePhysicalVolumeFraction();
}

◆ ~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.

{
  // Purpose: calculate new weights via moment-fitting method
  std::vector<Point> elem_nodes(_n_nodes, Point(0.0, 0.0, 0.0));
  std::vector<std::vector<Real>> wsg;
 
  for (unsigned int i = 0; i < _n_nodes; ++i)
    elem_nodes[i] = (*_nodes[i]);
 
  if (_efa_elem2d.isPartial() && _n_qpoints <= 6) // ONLY work for <= 6 q_points
  {
    std::vector<std::vector<Real>> tsg;
    getPhysicalQuadraturePoints(tsg); // get tsg - QPs within partial element
    solveMomentFitting(
        _n_nodes, _n_qpoints, elem_nodes, tsg, wsg); // get wsg - QPs from moment-fitting
    _new_weights.resize(wsg.size(), 1.0);
    for (unsigned int i = 0; i < wsg.size(); ++i)
      _new_weights[i] = wsg[i][2]; // weight multiplier
  }
  else
    _new_weights.resize(_n_qpoints, _physical_volfrac);
}

◆ computePhysicalFaceAreaFraction()

void XFEMCutElem2D::computePhysicalFaceAreaFraction ( unsigned int side )

virtual

Computes the surface area fraction of the element side.

Parameters

side	The side of the element

find a fragment edge which is covered by element side

Implements XFEMCutElem.

Definition at line 80 of file XFEMCutElem2D.C.

{
  Real frag_surf = 0.0;
 
  EFAEdge * edge = _efa_elem2d.getEdge(side);
 
  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
  {
    EFANode * node_1 = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
    EFANode * node_2 = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
 
    if (edge->containsNode(node_1) && edge->containsNode(node_2))
    {
      Point edge_p1 = getNodeCoordinates(node_1);
      Point edge_p2 = getNodeCoordinates(node_2);
      frag_surf = (edge_p1 - edge_p2).norm();
      _physical_areafrac[side] = frag_surf / _elem_side_area[side];
      return;
    }
  }
  _physical_areafrac[side] = 1.0;
}

◆ computePhysicalVolumeFraction()

void XFEMCutElem2D::computePhysicalVolumeFraction ( )

virtual

Computes the volume fraction of the element fragment.

Implements XFEMCutElem.

Definition at line 63 of file XFEMCutElem2D.C.

{
  Real frag_vol = 0.0;
 
  // Calculate area of entire element and fragment using the formula:
  // A = 1/2 sum_{i=0}^{n-1} (x_i y_{i+1} - x_{i+1} y{i})
 
  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
  {
    Point edge_p1 = getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(0));
    Point edge_p2 = getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(1));
    frag_vol += 0.5 * (edge_p1(0) - edge_p2(0)) * (edge_p1(1) + edge_p2(1));
  }
  _physical_volfrac = frag_vol / _elem_volume;
}

Referenced by XFEMCutElem2D().

◆ 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

qrule	The standard MOOSE face quadrature rule
xfem_qrule	The integration scheme for the cut element (We use surface area fraction only)
q_points	The quadrature points for the element side
side	The side of the element

Definition at line 92 of file XFEMCutElem.C.

{
  _new_face_weights[side].clear();
  _new_face_weights[side].resize(qrule->n_points(), 1.0);
 
  computePhysicalFaceAreaFraction(side);
  Real surffrac = getPhysicalFaceAreaFraction(side);
  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
    _new_face_weights[side][qp] = surffrac;
 
  _have_face_weights[side] = true;
}

Referenced by XFEMCutElem::getFaceWeightMultipliers().

◆ 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

qrule	The standard MOOSE quadrature rule
xfem_qrule	The integration scheme for the cut element
q_points	The quadrature points for the element

Definition at line 109 of file XFEMCutElem.C.

{
  _new_weights.clear();
 
  _new_weights.resize(qrule->n_points(), 1.0);
 
  switch (xfem_qrule)
  {
    case Xfem::VOLFRAC:
    {
      computePhysicalVolumeFraction();
      Real volfrac = getPhysicalVolumeFraction();
      for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
        _new_weights[qp] = volfrac;
      break;
    }
    case Xfem::MOMENT_FITTING:
    {
      // These are the coordinates in parametric coordinates
      _qp_points = qrule->get_points();
      _qp_weights = qrule->get_weights();
 
      computeMomentFittingWeights();
 
      // Blend weights from moment fitting and volume fraction to avoid negative weights
      Real alpha = 1.0;
      if (*std::min_element(_new_weights.begin(), _new_weights.end()) < 0.0)
      {
        // One or more of the weights computed by moment fitting is negative.
        // Blend moment and volume fraction weights to keep them nonnegative.
        // Find the largest value of alpha that will keep all of the weights nonnegative.
        for (unsigned int i = 0; i < _n_qpoints; ++i)
        {
          const Real denominator = _physical_volfrac - _new_weights[i];
          if (denominator > 0.0) // Negative values would give a negative value for alpha, which
                                 // must be between 0 and 1.
          {
            const Real alpha_i = _physical_volfrac / denominator;
            if (alpha_i < alpha)
              alpha = alpha_i;
          }
        }
        for (unsigned int i = 0; i < _n_qpoints; ++i)
        {
          _new_weights[i] = alpha * _new_weights[i] + (1.0 - alpha) * _physical_volfrac;
          if (_new_weights[i] < 0.0) // We can end up with small (roundoff) negative weights
            _new_weights[i] = 0.0;
        }
      }
      break;
    }
    case Xfem::DIRECT: // remove q-points outside the partial element's physical domain
    {
      bool nonzero = false;
      for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
      {
        // q_points contains quadrature point locations in physical coordinates
        if (isPointPhysical(q_points[qp]))
        {
          nonzero = true;
          _new_weights[qp] = 1.0;
        }
        else
          _new_weights[qp] = 0.0;
      }
      if (!nonzero)
      {
        // Set the weights to a small value (1e-3) to avoid having DOFs
        // with zeros on the diagonal, which occurs for nodes that are
        // connected to no physical material.
        for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
          _new_weights[qp] = 1e-3;
      }
      break;
    }
    default:
      mooseError("Undefined option for XFEM_QRULE");
  }
  _have_weights = true;
}

Referenced by XFEMCutElem::getWeightMultipliers().

◆ getCrackTipOriginAndDirection()

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

virtual

Implements XFEMCutElem.

Definition at line 180 of file XFEMCutElem2D.C.

{
  // TODO: two cut plane case is not working
  std::vector<EFANode *> cut_line_nodes;
  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
  {
    if (_efa_elem2d.getFragment(0)->isEdgeInterior(i))
    {
      std::vector<EFANode *> node_line(2, nullptr);
      node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
      node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
      if (node_line[1]->id() == tip_id)
      {
        cut_line_nodes.push_back(node_line[0]);
        cut_line_nodes.push_back(node_line[1]);
      }
      else if (node_line[0]->id() == tip_id)
      {
        node_line[1] = node_line[0];
        node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
        cut_line_nodes.push_back(node_line[0]);
        cut_line_nodes.push_back(node_line[1]);
      }
    }
  }
  if (cut_line_nodes.size() == 0)
    mooseError("no cut line found in this element");
 
  Point cut_line_p1 = getNodeCoordinates(cut_line_nodes[0]);
  Point cut_line_p2 = getNodeCoordinates(cut_line_nodes[1]);
  Point cut_line = cut_line_p2 - cut_line_p1;
  Real len = std::sqrt(cut_line.norm_sq());
  cut_line /= len;
  origin = cut_line_p2;
  direction = Point(cut_line(0), cut_line(1), 0.0);
}

◆ getCutPlaneNormal()

Point XFEMCutElem2D::getCutPlaneNormal	(	unsigned int	plane_id,
		MeshBase *	displaced_mesh = `nullptr`
	)		const

virtual

Implements XFEMCutElem.

Definition at line 151 of file XFEMCutElem2D.C.

{
  Point normal(0.0, 0.0, 0.0);
  std::vector<std::vector<EFANode *>> cut_line_nodes;
  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
  {
    if (_efa_elem2d.getFragment(0)->isEdgeInterior(i))
    {
      std::vector<EFANode *> node_line(2, nullptr);
      node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
      node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
      cut_line_nodes.push_back(node_line);
    }
  }
  if (cut_line_nodes.size() == 0)
    mooseError("no cut line found in this element");
  if (plane_id < cut_line_nodes.size()) // valid plane_id
  {
    Point cut_line_p1 = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
    Point cut_line_p2 = getNodeCoordinates(cut_line_nodes[plane_id][1], displaced_mesh);
    Point cut_line = cut_line_p2 - cut_line_p1;
    Real len = std::sqrt(cut_line.norm_sq());
    cut_line /= len;
    normal = Point(cut_line(1), -cut_line(0), 0.0);
  }
  return normal;
}

Referenced by getIntersectionInfo().

◆ getCutPlaneOrigin()

Point XFEMCutElem2D::getCutPlaneOrigin	(	unsigned int	plane_id,
		MeshBase *	displaced_mesh = `nullptr`
	)		const

virtual

Implements XFEMCutElem.

Definition at line 129 of file XFEMCutElem2D.C.

{
  Point orig(0.0, 0.0, 0.0);
  std::vector<std::vector<EFANode *>> cut_line_nodes;
  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
  {
    if (_efa_elem2d.getFragment(0)->isEdgeInterior(i))
    {
      std::vector<EFANode *> node_line(2, nullptr);
      node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
      node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
      cut_line_nodes.push_back(node_line);
    }
  }
  if (cut_line_nodes.size() == 0)
    mooseError("no cut line found in this element");
  if (plane_id < cut_line_nodes.size()) // valid plane_id
    orig = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
  return orig;
}

◆ getEFAElement()

const EFAElement * XFEMCutElem2D::getEFAElement ( ) const

virtual

Implements XFEMCutElem.

Definition at line 236 of file XFEMCutElem2D.C.

{
  return &_efa_elem2d;
}

◆ 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 77 of file XFEMCutElem.C.

{
  if (!_have_face_weights[side])
    computeXFEMFaceWeights(qrule, xfem_qrule, q_points, side);
 
  face_weights.resize(_new_face_weights[side].size());
  for (unsigned int qp = 0; qp < _new_face_weights[side].size(); ++qp)
    face_weights[qp] = _new_face_weights[side][qp];
}

Referenced by XFEM::getXFEMFaceWeights().

◆ getFragmentFaces()

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

virtual

Implements XFEMCutElem.

Definition at line 220 of file XFEMCutElem2D.C.

{
  frag_faces.clear();
  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
  {
    std::vector<Point> edge_points(2, Point(0.0, 0.0, 0.0));
    edge_points[0] =
        getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(0), displaced_mesh);
    edge_points[1] =
        getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(1), displaced_mesh);
    frag_faces.push_back(edge_points);
  }
}

◆ getIntersectionInfo()

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

virtual

Implements XFEMCutElem.

Definition at line 441 of file XFEMCutElem2D.C.

{
  intersectionPoints.resize(2); // 2D: the intersection line is straight and can be represented by
                                // two ending points.(may have issues with double cuts case!)
 
  std::vector<std::vector<EFANode *>> cut_line_nodes;
  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
  {
    if (_efa_elem2d.getFragment(0)->isEdgeInterior(i))
    {
      std::vector<EFANode *> node_line(2, nullptr);
      node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
      node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
      cut_line_nodes.push_back(node_line);
    }
  }
  if (cut_line_nodes.size() == 0)
    mooseError("No cut line found in this element");
 
  if (plane_id < cut_line_nodes.size()) // valid plane_id
  {
    intersectionPoints[0] = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
    intersectionPoints[1] = getNodeCoordinates(cut_line_nodes[plane_id][1], displaced_mesh);
  }
 
  normal = getCutPlaneNormal(plane_id, displaced_mesh);
}

◆ getMomentFittingWeight()

Real XFEMCutElem::getMomentFittingWeight ( unsigned int i_qp ) const

inherited

◆ getNodeCoordinates()

Point XFEMCutElem2D::getNodeCoordinates	(	EFANode *	node,
		MeshBase *	displaced_mesh = `nullptr`
	)		const

privatevirtual

Implements XFEMCutElem.

Definition at line 36 of file XFEMCutElem2D.C.

{
  Point node_coor(0.0, 0.0, 0.0);
  std::vector<EFANode *> master_nodes;
  std::vector<Point> master_points;
  std::vector<double> master_weights;
 
  _efa_elem2d.getMasterInfo(CEMnode, master_nodes, master_weights);
  for (unsigned int i = 0; i < master_nodes.size(); ++i)
  {
    if (master_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
    {
      Node * node = _nodes[master_nodes[i]->id()];
      if (displaced_mesh)
        node = displaced_mesh->node_ptr(node->id());
      Point node_p((*node)(0), (*node)(1), 0.0);
      master_points.push_back(node_p);
    }
    else
      mooseError("master nodes must be local");
  }
  for (unsigned int i = 0; i < master_nodes.size(); ++i)
    node_coor += master_weights[i] * master_points[i];
  return node_coor;
}

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

◆ getPhysicalFaceAreaFraction()

Real XFEMCutElem::getPhysicalFaceAreaFraction ( unsigned int side ) const

inherited

Returns the surface area fraction of the element side.

Parameters

side	The side of the element

Definition at line 57 of file XFEMCutElem.C.

{
  return _physical_areafrac[side];
}

Referenced by XFEMCutElem::computeXFEMFaceWeights().

◆ getPhysicalQuadraturePoints()

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

private

Definition at line 252 of file XFEMCutElem2D.C.

{
  // Get the coords for parial element nodes
  EFAFragment2D * frag = _efa_elem2d.getFragment(0);
  unsigned int nnd_pe = frag->numEdges();
  std::vector<Point> frag_points(nnd_pe, Point(0.0, 0.0, 0.0)); // nodal coord of partial elem
  Real jac = 0.0;
 
  for (unsigned int j = 0; j < nnd_pe; ++j)
    frag_points[j] = getNodeCoordinates(frag->getEdge(j)->getNode(0));
 
  // Get centroid coords for partial elements
  Point xcrd(0.0, 0.0, 0.0);
  for (unsigned int j = 0; j < nnd_pe; ++j)
    xcrd += frag_points[j];
  xcrd *= (1.0 / nnd_pe);
 
  // Get tsg - the physical coords of Gaussian Q-points for partial elements
  if ((nnd_pe == 3) || (nnd_pe == 4)) // partial element is a triangle or quad
  {
    std::vector<std::vector<Real>> sg2;
    std::vector<std::vector<Real>> shape(nnd_pe, std::vector<Real>(3, 0.0));
    Xfem::stdQuadr2D(nnd_pe, 2, sg2);
    for (unsigned int l = 0; l < sg2.size(); ++l)
    {
      Xfem::shapeFunc2D(nnd_pe, sg2[l], frag_points, shape, jac, true); // Get shape
      std::vector<Real> tsg_line(3, 0.0);
      for (unsigned int k = 0; k < nnd_pe; ++k)
      {
        tsg_line[0] += shape[k][2] * frag_points[k](0);
        tsg_line[1] += shape[k][2] * frag_points[k](1);
      }
      if (nnd_pe == 3)                 // tri partial elem
        tsg_line[2] = sg2[l][3] * jac; // total weights
      else                             // quad partial elem
        tsg_line[2] = sg2[l][2] * jac; // total weights
      tsg.push_back(tsg_line);
    }
  }
  else if (nnd_pe >= 5) // partial element is a polygon
  {
    for (unsigned int j = 0; j < nnd_pe; ++j) // loop all sub-tris
    {
      std::vector<std::vector<Real>> shape(3, std::vector<Real>(3, 0.0));
      std::vector<Point> subtri_points(3, Point(0.0, 0.0, 0.0)); // sub-tri nodal coords
 
      int jplus1 = j < nnd_pe - 1 ? j + 1 : 0;
      subtri_points[0] = xcrd;
      subtri_points[1] = frag_points[j];
      subtri_points[2] = frag_points[jplus1];
 
      std::vector<std::vector<Real>> sg2;
      Xfem::stdQuadr2D(3, 2, sg2);                  // get sg2
      for (unsigned int l = 0; l < sg2.size(); ++l) // loop all int pts on a sub-tri
      {
        Xfem::shapeFunc2D(3, sg2[l], subtri_points, shape, jac, true); // Get shape
        std::vector<Real> tsg_line(3, 0.0);
        for (unsigned int k = 0; k < 3; ++k) // loop sub-tri nodes
        {
          tsg_line[0] += shape[k][2] * subtri_points[k](0);
          tsg_line[1] += shape[k][2] * subtri_points[k](1);
        }
        tsg_line[2] = sg2[l][3] * jac;
        tsg.push_back(tsg_line);
      }
    }
  }
  else
    mooseError("Invalid partial element!");
}

Referenced by computeMomentFittingWeights().

◆ getPhysicalVolumeFraction()

Real XFEMCutElem::getPhysicalVolumeFraction ( ) const

inherited

Returns the volume fraction of the element fragment.

Definition at line 51 of file XFEMCutElem.C.

{
  return _physical_volfrac;
}

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

◆ getWeightMultipliers()

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

inherited

Definition at line 63 of file XFEMCutElem.C.

{
  if (!_have_weights)
    computeXFEMWeights(qrule, xfem_qrule, q_points);
 
  weights.resize(_new_weights.size());
  for (unsigned int qp = 0; qp < _new_weights.size(); ++qp)
    weights[qp] = _new_weights[qp];
}

Referenced by XFEM::getXFEMWeights().

◆ isPointPhysical()

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

inherited

Definition at line 193 of file XFEMCutElem.C.

{
  // determine whether the point is inside the physical domain of a partial element
  bool physical_flag = true;
  unsigned int n_cut_planes = numCutPlanes();
  for (unsigned int plane_id = 0; plane_id < n_cut_planes; ++plane_id)
  {
    Point origin = getCutPlaneOrigin(plane_id);
    Point normal = getCutPlaneNormal(plane_id);
    Point origin2qp = p - origin;
    if (origin2qp * normal > 0.0)
    {
      physical_flag = false; // Point outside pysical domain
      break;
    }
  }
  return physical_flag;
}

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

◆ numCutPlanes()

unsigned int XFEMCutElem2D::numCutPlanes ( ) const

virtual

Implements XFEMCutElem.

Definition at line 242 of file XFEMCutElem2D.C.

{
  unsigned int counter = 0;
  for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
    if (_efa_elem2d.getFragment(0)->isEdgeInterior(i))
      counter += 1;
  return counter;
}

◆ 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.

{
  // Get physical coords for the new six-point rule
  std::vector<std::vector<Real>> shape(nen, std::vector<Real>(3, 0.0));
  std::vector<std::vector<Real>> wss;
 
  if (nen == 4)
  {
    wss.resize(_qp_points.size());
    for (unsigned int i = 0; i < _qp_points.size(); i++)
    {
      wss[i].resize(3);
      wss[i][0] = _qp_points[i](0);
      wss[i][1] = _qp_points[i](1);
      wss[i][2] = _qp_weights[i];
    }
  }
  else if (nen == 3)
  {
    wss.resize(_qp_points.size());
    for (unsigned int i = 0; i < _qp_points.size(); i++)
    {
      wss[i].resize(4);
      wss[i][0] = _qp_points[i](0);
      wss[i][1] = _qp_points[i](1);
      wss[i][2] = 1.0 - _qp_points[i](0) - _qp_points[i](1);
      wss[i][3] = _qp_weights[i];
    }
  }
  else
    mooseError("Invalid element");
 
  wsg.resize(wss.size());
  for (unsigned int i = 0; i < wsg.size(); ++i)
    wsg[i].resize(3, 0.0);
  Real jac = 0.0;
  std::vector<Real> old_weights(wss.size(), 0.0);
 
  for (unsigned int l = 0; l < wsg.size(); ++l)
  {
    Xfem::shapeFunc2D(nen, wss[l], elem_nodes, shape, jac, true); // Get shape
    if (nen == 4)                                                 // 2D quad elem
      old_weights[l] = wss[l][2] * jac;                           // weights for total element
    else if (nen == 3)                                            // 2D triangle elem
      old_weights[l] = wss[l][3] * jac;
    else
      mooseError("Invalid element!");
    for (unsigned int k = 0; k < nen; ++k) // physical coords of Q-pts
    {
      wsg[l][0] += shape[k][2] * elem_nodes[k](0);
      wsg[l][1] += shape[k][2] * elem_nodes[k](1);
    }
  }
 
  // Compute weights via moment fitting
  Real * A;
  A = new Real[wsg.size() * wsg.size()];
  unsigned ind = 0;
  for (unsigned int i = 0; i < wsg.size(); ++i)
  {
    A[ind] = 1.0; // const
    if (nqp > 1)
      A[1 + ind] = wsg[i][0]; // x
    if (nqp > 2)
      A[2 + ind] = wsg[i][1]; // y
    if (nqp > 3)
      A[3 + ind] = wsg[i][0] * wsg[i][1]; // x*y
    if (nqp > 4)
      A[4 + ind] = wsg[i][0] * wsg[i][0]; // x^2
    if (nqp > 5)
      A[5 + ind] = wsg[i][1] * wsg[i][1]; // y^2
    if (nqp > 6)
      mooseError("Q-points of more than 6 are not allowed now!");
    ind = ind + nqp;
  }
 
  Real * b;
  b = new Real[wsg.size()];
  for (unsigned int i = 0; i < wsg.size(); ++i)
    b[i] = 0.0;
  for (unsigned int i = 0; i < tsg.size(); ++i)
  {
    b[0] += tsg[i][2];
    if (nqp > 1)
      b[1] += tsg[i][2] * tsg[i][0];
    if (nqp > 2)
      b[2] += tsg[i][2] * tsg[i][1];
    if (nqp > 3)
      b[3] += tsg[i][2] * tsg[i][0] * tsg[i][1];
    if (nqp > 4)
      b[4] += tsg[i][2] * tsg[i][0] * tsg[i][0];
    if (nqp > 5)
      b[5] += tsg[i][2] * tsg[i][1] * tsg[i][1];
    if (nqp > 6)
      mooseError("Q-points of more than 6 are not allowed now!");
  }
 
  PetscBLASInt nrhs = 1;
  PetscBLASInt info;
  PetscBLASInt n = wsg.size();
  std::vector<PetscBLASInt> ipiv(n);
 
  LAPACKgesv_(&n, &nrhs, A, &n, &ipiv[0], b, &n, &info);
 
  for (unsigned int i = 0; i < wsg.size(); ++i)
    wsg[i][2] = b[i] / old_weights[i]; // get the multiplier
 
  // delete arrays
  delete[] A;
  delete[] b;
}

Referenced by computeMomentFittingWeights().

Member Data Documentation

◆ _efa_elem2d

EFAElement2D XFEMCutElem2D::_efa_elem2d

private

Definition at line 41 of file XFEMCutElem2D.h.

Referenced by computeMomentFittingWeights(), computePhysicalFaceAreaFraction(), computePhysicalVolumeFraction(), getCrackTipOriginAndDirection(), getCutPlaneNormal(), getCutPlaneOrigin(), getEFAElement(), getFragmentFaces(), getIntersectionInfo(), getNodeCoordinates(), getPhysicalQuadraturePoints(), and numCutPlanes().

◆ _elem_side_area

std::vector<Real> XFEMCutElem::_elem_side_area

protectedinherited

Definition at line 49 of file XFEMCutElem.h.

Referenced by computePhysicalFaceAreaFraction(), XFEMCutElem3D::computePhysicalFaceAreaFraction(), and XFEMCutElem::XFEMCutElem().

◆ _elem_volume

Real XFEMCutElem::_elem_volume

protectedinherited

Definition at line 48 of file XFEMCutElem.h.

Referenced by computePhysicalVolumeFraction(), XFEMCutElem3D::computePhysicalVolumeFraction(), and XFEMCutElem::XFEMCutElem().

◆ _have_face_weights

std::vector<bool> XFEMCutElem::_have_face_weights

protectedinherited

Definition at line 53 of file XFEMCutElem.h.

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

◆ _have_weights

bool XFEMCutElem::_have_weights

protectedinherited

Definition at line 52 of file XFEMCutElem.h.

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

◆ _n_nodes

unsigned int XFEMCutElem::_n_nodes

protectedinherited

Definition at line 42 of file XFEMCutElem.h.

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

◆ _n_qpoints

unsigned int XFEMCutElem::_n_qpoints

protectedinherited

Definition at line 43 of file XFEMCutElem.h.

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

◆ _n_sides

unsigned int XFEMCutElem::_n_sides

protectedinherited

Definition at line 44 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 57 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 55 of file XFEMCutElem.h.

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

◆ _nodes

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

protectedinherited

Definition at line 45 of file XFEMCutElem.h.

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

◆ _physical_areafrac

std::vector<Real> XFEMCutElem::_physical_areafrac

protectedinherited

Definition at line 51 of file XFEMCutElem.h.

Referenced by computePhysicalFaceAreaFraction(), XFEMCutElem3D::computePhysicalFaceAreaFraction(), XFEMCutElem::getPhysicalFaceAreaFraction(), and XFEMCutElem::XFEMCutElem().

◆ _physical_volfrac

Real XFEMCutElem::_physical_volfrac

protectedinherited

Definition at line 50 of file XFEMCutElem.h.

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

◆ _qp_points

std::vector<Point> XFEMCutElem::_qp_points

protectedinherited

Definition at line 46 of file XFEMCutElem.h.

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

◆ _qp_weights

std::vector<Real> XFEMCutElem::_qp_weights

protectedinherited

Definition at line 47 of file XFEMCutElem.h.

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

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

xfem/include/base/XFEMCutElem2D.h
xfem/src/base/XFEMCutElem2D.C

Public Member Functions

Protected Attributes

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ XFEMCutElem2D()

◆ ~XFEMCutElem2D()

Member Function Documentation

◆ computeMomentFittingWeights()

◆ computePhysicalFaceAreaFraction()

◆ computePhysicalVolumeFraction()

◆ computeXFEMFaceWeights()

◆ computeXFEMWeights()

◆ getCrackTipOriginAndDirection()

◆ getCutPlaneNormal()

◆ getCutPlaneOrigin()

◆ getEFAElement()

◆ getFaceWeightMultipliers()

◆ getFragmentFaces()

◆ getIntersectionInfo()

◆ getMomentFittingWeight()

◆ getNodeCoordinates()

◆ getPhysicalFaceAreaFraction()

◆ getPhysicalQuadraturePoints()

◆ getPhysicalVolumeFraction()

◆ getWeightMultipliers()

◆ isPointPhysical()

◆ numCutPlanes()

◆ setQuadraturePointsAndWeights()

◆ solveMomentFitting()

Member Data Documentation

◆ _efa_elem2d

◆ _elem_side_area

◆ _elem_volume

◆ _have_face_weights

◆ _have_weights

◆ _n_nodes

◆ _n_qpoints

◆ _n_sides

◆ _new_face_weights

◆ _new_weights

◆ _nodes

◆ _physical_areafrac

◆ _physical_volfrac

◆ _qp_points

◆ _qp_weights