This is the XFEM class. More...

#include <XFEM.h>

Inheritance diagram for XFEM:

Public Member Functions
	XFEM (const InputParameters &params)

	~XFEM ()

void	addGeometricCut (GeometricCutUserObject *geometric_cut)

void	addStateMarkedElem (unsigned int elem_id, RealVectorValue &normal)

void	addStateMarkedElem (unsigned int elem_id, RealVectorValue &normal, unsigned int marked_side)

void	addStateMarkedFrag (unsigned int elem_id, RealVectorValue &normal)

void	clearStateMarkedElems ()

void	addGeomMarkedElem2D (const unsigned int elem_id, const Xfem::GeomMarkedElemInfo2D geom_info, const unsigned int interface_id)
	Add information about a new cut to be performed on a specific 2d element. More...

void	addGeomMarkedElem3D (const unsigned int elem_id, const Xfem::GeomMarkedElemInfo3D geom_info, const unsigned int interface_id)
	Add information about a new cut to be performed on a specific 3d element. More...

void	clearGeomMarkedElems ()
	Clear out the list of elements to be marked for cutting. More...

virtual bool	update (Real time, const std::vector< std::shared_ptr< NonlinearSystemBase >> &nl, AuxiliarySystem &aux) override

virtual void	initSolution (const std::vector< std::shared_ptr< NonlinearSystemBase >> &nl, AuxiliarySystem &aux) override

void	buildEFAMesh ()

bool	markCuts (Real time)

bool	markCutEdgesByGeometry ()

bool	markCutEdgesByState (Real time)

bool	markCutFacesByGeometry ()

bool	markCutFacesByState ()

bool	initCutIntersectionEdge (Point cut_origin, RealVectorValue cut_normal, Point &edge_p1, Point &edge_p2, Real &dist)

bool	cutMeshWithEFA (const std::vector< std::shared_ptr< NonlinearSystemBase >> &nl, AuxiliarySystem &aux)

bool	healMesh ()
	Potentially heal the mesh by merging some of the pairs of partial elements cut by XFEM back into single elements if indicated by the cutting objects. More...

virtual bool	updateHeal () override

Point	getEFANodeCoords (EFANode CEMnode, EFAElement CEMElem, const Elem elem, MeshBase displaced_mesh=nullptr) const

Real	getPhysicalVolumeFraction (const Elem *elem) const
	Get the volume fraction of an element that is physical. More...

bool	isPointInsidePhysicalDomain (const Elem *elem, const Point &point) const
	Return true if the point is inside the element physical domain Note: if this element is not cut, return true too. More...

Real	getCutPlane (const Elem *elem, const Xfem::XFEM_CUTPLANE_QUANTITY quantity, unsigned int plane_id) const
	Get specified component of normal or origin for cut plane for a given element. More...

bool	isElemAtCrackTip (const Elem *elem) const

bool	isElemCut (const Elem elem, XFEMCutElem &xfce) const

bool	isElemCut (const Elem *elem) const

void	getFragmentFaces (const Elem *elem, std::vector< std::vector< Point >> &frag_faces, bool displaced_mesh=false) const

void	storeCrackTipOriginAndDirection ()

void	correctCrackExtensionDirection (const Elem elem, EFAElement2D CEMElem, EFAEdge *orig_edge, Point normal, Point crack_tip_origin, Point crack_tip_direction, Real &distance_keep, unsigned int &edge_id_keep, Point &normal_keep)

void	getCrackTipOrigin (std::map< unsigned int, const Elem *> &elem_id_crack_tip, std::vector< Point > &crack_front_points)

Xfem::XFEM_QRULE &	getXFEMQRule ()

void	setXFEMQRule (std::string &xfem_qrule)

void	setCrackGrowthMethod (bool use_crack_growth_increment, Real crack_growth_increment)

void	setDebugOutputLevel (unsigned int debug_output_level)
	Controls amount of debugging information output. More...

void	setMinWeightMultiplier (Real min_weight_multiplier)
	Controls the minimum average weight multiplier for each element. More...

virtual bool	getXFEMWeights (MooseArray< Real > &weights, const Elem elem, QBase qrule, const MooseArray< Point > &q_points) override

virtual bool	getXFEMFaceWeights (MooseArray< Real > &weights, const Elem elem, QBase qrule, const MooseArray< Point > &q_points, unsigned int side) override

virtual const ElementPairLocator::ElementPairList *	getXFEMCutElemPairs (unsigned int interface_id)
	Get the list of cut element pairs corresponding to a given interface ID. More...

virtual const ElementPairLocator::ElementPairList *	getXFEMDisplacedCutElemPairs (unsigned int interface_id)
	Get the list of cut element pairs on the displaced mesh corresponding to a given interface ID. More...

virtual unsigned int	getGeometricCutID (const GeometricCutUserObject *gcu)
	Get the interface ID corresponding to a given GeometricCutUserObject. More...

virtual void	getXFEMIntersectionInfo (const Elem *elem, unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, bool displaced_mesh=false) const

virtual void	getXFEMqRuleOnLine (std::vector< Point > &intersection_points, std::vector< Point > &quad_pts, std::vector< Real > &quad_wts) const

virtual void	getXFEMqRuleOnSurface (std::vector< Point > &intersection_points, std::vector< Point > &quad_pts, std::vector< Real > &quad_wts) const

bool	has_secondary_cut ()

EFAElement2D *	getEFAElem2D (const Elem *elem)
	Get the EFAElement2D object for a specified libMesh element. More...

EFAElement3D *	getEFAElem3D (const Elem *elem)
	Get the EFAElement3D object for a specified libMesh element. More...

void	getFragmentEdges (const Elem elem, EFAElement2D CEMElem, std::vector< std::vector< Point >> &frag_edges) const

void	getFragmentFaces (const Elem elem, EFAElement3D CEMElem, std::vector< std::vector< Point >> &frag_faces) const

const std::map< const Elem *, std::vector< Point > > &	getCrackTipOriginMap () const

CutSubdomainID	getCutSubdomainID (const GeometricCutUserObject gcuo, const Elem cut_elem, const Elem *parent_elem=nullptr) const
	Determine which cut subdomain the element belongs to relative to the cut. More...

void	setMesh (MooseMesh *mesh)

void	setDisplacedMesh (MooseMesh *displaced_mesh)

void	setMaterialData (const std::vector< MaterialData * > &data)

void	setBoundaryMaterialData (const std::vector< MaterialData * > &data)

virtual bool	getXFEMWeights (MooseArray< Real > &weights, const Elem elem, libMesh::QBase qrule, const MooseArray< Point > &q_points)=0

virtual bool	getXFEMFaceWeights (MooseArray< Real > &weights, const Elem elem, libMesh::QBase qrule, const MooseArray< Point > &q_points, unsigned int side)=0

Public Attributes
const ConsoleStream	_console

Protected Attributes
FEProblemBase *	_fe_problem

std::vector< MaterialData *>	_material_data

std::vector< MaterialData *>	_bnd_material_data

MooseMesh *	_moose_mesh

MooseMesh *	_moose_displaced_mesh

MeshBase *	_mesh

MeshBase *	_displaced_mesh

Private Member Functions
void	storeSolutionForNode (const Node node_to_store_to, const Node node_to_store_from, SystemBase &sys, std::map< unique_id_type, std::vector< Real >> &stored_solution, const NumericVector< Number > &current_solution, const NumericVector< Number > &old_solution, const NumericVector< Number > &older_solution)
	Store the solution in stored_solution for a given node. More...

void	storeSolutionForElement (const Elem elem_to_store_to, const Elem elem_to_store_from, SystemBase &sys, std::map< unique_id_type, std::vector< Real >> &stored_solution, const NumericVector< Number > &current_solution, const NumericVector< Number > &old_solution, const NumericVector< Number > &older_solution)
	Store the solution in stored_solution for a given element. More...

void	setSolution (SystemBase &sys, const std::map< unique_id_type, std::vector< Real >> &stored_solution, NumericVector< Number > &current_solution, NumericVector< Number > &old_solution, NumericVector< Number > &older_solution)
	Set the solution for all locally-owned nodes/elements that have stored values. More...

void	setSolutionForDOFs (const std::vector< Real > &stored_solution, const std::vector< dof_id_type > &stored_solution_dofs, NumericVector< Number > &current_solution, NumericVector< Number > &old_solution, NumericVector< Number > &older_solution)
	Set the solution for a set of DOFs. More...

std::vector< dof_id_type >	getElementSolutionDofs (const Elem *elem, SystemBase &sys) const
	Get a vector of the dof indices for all components of all variables associated with an element. More...

std::vector< dof_id_type >	getNodeSolutionDofs (const Node *node, SystemBase &sys) const
	Get a vector of the dof indices for all components of all variables associated with a node. More...

const GeometricCutUserObject *	getGeometricCutForElem (const Elem *elem) const
	Get the GeometricCutUserObject associated with an element. More...

void	storeMaterialPropertiesForElementHelper (const Elem *elem, MaterialPropertyStorage &storage)

void	storeMaterialPropertiesForElement (const Elem parent_elem, const Elem child_elem)
	Helper function to store the material properties of a healed element. More...

void	loadMaterialPropertiesForElementHelper (const Elem *elem, const Xfem::CachedMaterialProperties &cached_props, MaterialPropertyStorage &storage) const
	Load the material properties. More...

void	loadMaterialPropertiesForElement (const Elem elem, const Elem elem_from, std::unordered_map< const Elem *, Xfem::CutElemInfo > &cached_cei) const
	Helper function to store the material properties of a healed element. More...

const Node *	pickFirstPhysicalNode (const Elem e, const Elem e0) const
	Return the first node in the provided element that is found to be in the physical domain. More...

Private Attributes
bool	_has_secondary_cut

Xfem::XFEM_QRULE	_XFEM_qrule

bool	_use_crack_growth_increment

Real	_crack_growth_increment

std::vector< const GeometricCutUserObject * >	_geometric_cuts

std::map< unique_id_type, XFEMCutElem * >	_cut_elem_map

std::set< const Elem * >	_crack_tip_elems

std::set< const Elem * >	_crack_tip_elems_to_be_healed

std::map< unsigned int, ElementPairLocator::ElementPairList >	_sibling_elems

std::map< unsigned int, ElementPairLocator::ElementPairList >	_sibling_displaced_elems

std::map< const Elem *, std::vector< Point > >	_elem_crack_origin_direction_map

std::map< const Elem *, RealVectorValue >	_state_marked_elems

std::set< const Elem * >	_state_marked_frags

std::map< const Elem *, unsigned int >	_state_marked_elem_sides

std::map< const Elem *, std::vector< Xfem::GeomMarkedElemInfo2D > >	_geom_marked_elems_2d
	Data structure for storing information about all 2D elements to be cut by geometry. More...

std::map< const Elem *, std::vector< Xfem::GeomMarkedElemInfo3D > >	_geom_marked_elems_3d
	Data structure for storing information about all 3D elements to be cut by geometry. More...

std::map< unsigned int, std::set< unsigned int > >	_geom_marker_id_elems
	Data structure for storing the elements cut by specific geometric cutters. More...

std::map< const GeometricCutUserObject *, unsigned int >	_geom_marker_id_map
	Data structure for storing the GeommetricCutUserObjects and their corresponding id. More...

ElementFragmentAlgorithm	_efa_mesh

unsigned int	_debug_output_level
	Controls amount of debugging output information 0: None 1: Summary 2: Details on modifications to mesh 3: Full dump of element fragment algorithm mesh. More...

Real	_min_weight_multiplier
	The minimum average multiplier applied by XFEM to the standard quadrature weights to integrate partial elements. More...

std::map< unique_id_type, std::vector< Real > >	_cached_solution
	Data structure to store the nonlinear solution for nodes/elements affected by XFEM For each node/element, this is stored as a vector that contains all components of all applicable variables in an order defined by getElementSolutionDofs() or getNodeSolutionDofs(). More...

std::map< unique_id_type, std::vector< Real > >	_cached_aux_solution
	Data structure to store the auxiliary solution for nodes/elements affected by XFEM For each node/element, this is stored as a vector that contains all components of all applicable variables in an order defined by getElementSolutionDofs() or getNodeSolutionDofs(). More...

std::unordered_map< const Elem *, Xfem::CutElemInfo >	_geom_cut_elems
	All geometrically cut elements and their CutElemInfo during the current execution of XFEM_MARK. More...

std::unordered_map< const Elem *, Xfem::CutElemInfo >	_old_geom_cut_elems
	All geometrically cut elements and their CutElemInfo before the current execution of XFEM_MARK. More...

Detailed Description

This is the XFEM class.

This class implements algorithms for dynamic mesh modification in support of a phantom node approach for XFEM

Definition at line 107 of file XFEM.h.

Constructor & Destructor Documentation

◆ XFEM()

XFEM::XFEM ( const InputParameters & params )

explicit

Definition at line 36 of file XFEM.C.

   : XFEMInterface(params),
     _efa_mesh(Moose::out),
     _debug_output_level(1),
     _min_weight_multiplier(0.0)
 {
 #ifndef LIBMESH_ENABLE_UNIQUE_ID
   mooseError("MOOSE requires unique ids to be enabled in libmesh (configure with "
              "--enable-unique-id) to use XFEM!");
 #endif
   _has_secondary_cut = false;
 }

◆ ~XFEM()

XFEM::~XFEM ( )

Definition at line 49 of file XFEM.C.

 {
   for (std::map<unique_id_type, XFEMCutElem *>::iterator cemit = _cut_elem_map.begin();
        cemit != _cut_elem_map.end();
        ++cemit)
     delete cemit->second;
 }

Member Function Documentation

◆ addGeometricCut()

void XFEM::addGeometricCut ( GeometricCutUserObject * geometric_cut )

Definition at line 58 of file XFEM.C.

 {
   _geometric_cuts.push_back(geometric_cut);
 
   geometric_cut->setInterfaceID(_geometric_cuts.size() - 1);
 
   _geom_marker_id_map[geometric_cut] = _geometric_cuts.size() - 1;
 }

◆ addGeomMarkedElem2D()

void XFEM::addGeomMarkedElem2D	(	const unsigned int	elem_id,
		const Xfem::GeomMarkedElemInfo2D	geom_info,
		const unsigned int	interface_id
	)

Add information about a new cut to be performed on a specific 2d element.

Parameters

elem_id	The id of the element to be cut
geom_info	The object containing information about the cut to be performed
interface_id	The ID of the interface

Definition at line 160 of file XFEM.C.

 {
   Elem * elem = _mesh->elem_ptr(elem_id);
   _geom_marked_elems_2d[elem].push_back(geom_info);
   _geom_marker_id_elems[interface_id].insert(elem_id);
 }

◆ addGeomMarkedElem3D()

void XFEM::addGeomMarkedElem3D	(	const unsigned int	elem_id,
		const Xfem::GeomMarkedElemInfo3D	geom_info,
		const unsigned int	interface_id
	)

Add information about a new cut to be performed on a specific 3d element.

Parameters

elem_id	The id of the element to be cut
geom_info	The object containing information about the cut to be performed
interface_id	The ID of the interface

Definition at line 170 of file XFEM.C.

 {
   Elem * elem = _mesh->elem_ptr(elem_id);
   _geom_marked_elems_3d[elem].push_back(geom_info);
   _geom_marker_id_elems[interface_id].insert(elem_id);
 }

◆ addStateMarkedElem() [1/2]

void XFEM::addStateMarkedElem	(	unsigned int	elem_id,
		RealVectorValue &	normal
	)

Definition at line 114 of file XFEM.C.

Referenced by addStateMarkedElem(), and addStateMarkedFrag().

 {
   Elem * elem = _mesh->elem_ptr(elem_id);
   std::map<const Elem *, RealVectorValue>::iterator mit;
   mit = _state_marked_elems.find(elem);
   if (mit != _state_marked_elems.end())
     mooseError(" ERROR: element ", elem->id(), " already marked for crack growth.");
   _state_marked_elems[elem] = normal;
 }

◆ addStateMarkedElem() [2/2]

void XFEM::addStateMarkedElem	(	unsigned int	elem_id,
		RealVectorValue &	normal,
		unsigned int	marked_side
	)

Definition at line 125 of file XFEM.C.

 {
   addStateMarkedElem(elem_id, normal);
   Elem * elem = _mesh->elem_ptr(elem_id);
   std::map<const Elem *, unsigned int>::iterator mit;
   mit = _state_marked_elem_sides.find(elem);
   if (mit != _state_marked_elem_sides.end())
     mooseError(" ERROR: side of element ", elem->id(), " already marked for crack initiation.");
 
   _state_marked_elem_sides[elem] = marked_side;
 }

◆ addStateMarkedFrag()

void XFEM::addStateMarkedFrag	(	unsigned int	elem_id,
		RealVectorValue &	normal
	)

Definition at line 138 of file XFEM.C.

 {
   addStateMarkedElem(elem_id, normal);
   Elem * elem = _mesh->elem_ptr(elem_id);
   std::set<const Elem *>::iterator mit;
   mit = _state_marked_frags.find(elem);
   if (mit != _state_marked_frags.end())
     mooseError(
         " ERROR: element ", elem->id(), " already marked for fragment-secondary crack initiation.");
 
   _state_marked_frags.insert(elem);
 }

◆ buildEFAMesh()

void XFEM::buildEFAMesh ( )

Definition at line 339 of file XFEM.C.

Referenced by update(), and updateHeal().

 {
   _efa_mesh.reset();
 
   // Load all existing elements in to EFA mesh
   for (auto & elem : _mesh->element_ptr_range())
   {
     std::vector<unsigned int> quad;
     for (unsigned int i = 0; i < elem->n_nodes(); ++i)
       quad.push_back(elem->node_id(i));
 
     if (_mesh->mesh_dimension() == 2)
       _efa_mesh.add2DElement(quad, elem->id());
     else if (_mesh->mesh_dimension() == 3)
       _efa_mesh.add3DElement(quad, elem->id());
     else
       mooseError("XFEM only works for 2D and 3D");
   }
 
   // Restore fragment information for elements that have been previously cut
   for (auto & elem : _mesh->element_ptr_range())
   {
     std::map<unique_id_type, XFEMCutElem *>::iterator cemit = _cut_elem_map.find(elem->unique_id());
     if (cemit != _cut_elem_map.end())
     {
       XFEMCutElem * xfce = cemit->second;
       EFAElement * CEMElem = _efa_mesh.getElemByID(elem->id());
       _efa_mesh.restoreFragmentInfo(CEMElem, xfce->getEFAElement());
     }
   }
 
   // Must update edge neighbors before restore edge intersections. Otherwise, when we
   // add edge intersections, we do not have neighbor information to use.
   // Correction: no need to use neighbor info now
   _efa_mesh.updateEdgeNeighbors();
   _efa_mesh.initCrackTipTopology();
 }

◆ clearGeomMarkedElems()

void XFEM::clearGeomMarkedElems ( )

Clear out the list of elements to be marked for cutting.

Called after cutting is done.

Definition at line 180 of file XFEM.C.

Referenced by update().

 {
   _geom_marked_elems_2d.clear();
   _geom_marked_elems_3d.clear();
 }

◆ clearStateMarkedElems()

void XFEM::clearStateMarkedElems ( )

Definition at line 152 of file XFEM.C.

Referenced by update().

 {
   _state_marked_elems.clear();
   _state_marked_frags.clear();
   _state_marked_elem_sides.clear();
 }

◆ correctCrackExtensionDirection()

void XFEM::correctCrackExtensionDirection	(	const Elem *	elem,
		EFAElement2D *	CEMElem,
		EFAEdge *	orig_edge,
		Point	normal,
		Point	crack_tip_origin,
		Point	crack_tip_direction,
		Real &	distance_keep,
		unsigned int &	edge_id_keep,
		Point &	normal_keep
	)

Definition at line 447 of file XFEM.C.

Referenced by markCutEdgesByState().

 {
   std::vector<Point> edge_ends(2, Point(0.0, 0.0, 0.0));
   Point edge1(0.0, 0.0, 0.0);
   Point edge2(0.0, 0.0, 0.0);
   Point left_angle(0.0, 0.0, 0.0);
   Point right_angle(0.0, 0.0, 0.0);
   Point left_angle_normal(0.0, 0.0, 0.0);
   Point right_angle_normal(0.0, 0.0, 0.0);
   Point crack_direction_normal(0.0, 0.0, 0.0);
   Point edge1_to_tip(0.0, 0.0, 0.0);
   Point edge2_to_tip(0.0, 0.0, 0.0);
   Point edge1_to_tip_normal(0.0, 0.0, 0.0);
   Point edge2_to_tip_normal(0.0, 0.0, 0.0);
 
   Real cos_45 = std::cos(45.0 / 180.0 * 3.14159);
   Real sin_45 = std::sin(45.0 / 180.0 * 3.14159);
 
   left_angle(0) = cos_45 * crack_tip_direction(0) - sin_45 * crack_tip_direction(1);
   left_angle(1) = sin_45 * crack_tip_direction(0) + cos_45 * crack_tip_direction(1);
 
   right_angle(0) = cos_45 * crack_tip_direction(0) + sin_45 * crack_tip_direction(1);
   right_angle(1) = -sin_45 * crack_tip_direction(0) + cos_45 * crack_tip_direction(1);
 
   left_angle_normal(0) = -left_angle(1);
   left_angle_normal(1) = left_angle(0);
 
   right_angle_normal(0) = -right_angle(1);
   right_angle_normal(1) = right_angle(0);
 
   crack_direction_normal(0) = -crack_tip_direction(1);
   crack_direction_normal(1) = crack_tip_direction(0);
 
   Real angle_min = 0.0;
   Real distance = 0.0;
   unsigned int nsides = CEMElem->numEdges();
 
   for (unsigned int i = 0; i < nsides; ++i)
   {
     if (!orig_edge->isPartialOverlap(*CEMElem->getEdge(i)))
     {
       edge_ends[0] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
       edge_ends[1] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
 
       edge1_to_tip = (edge_ends[0] * 0.95 + edge_ends[1] * 0.05) - crack_tip_origin;
       edge2_to_tip = (edge_ends[0] * 0.05 + edge_ends[1] * 0.95) - crack_tip_origin;
 
       edge1_to_tip /= pow(edge1_to_tip.norm_sq(), 0.5);
       edge2_to_tip /= pow(edge2_to_tip.norm_sq(), 0.5);
 
       edge1_to_tip_normal(0) = -edge1_to_tip(1);
       edge1_to_tip_normal(1) = edge1_to_tip(0);
 
       edge2_to_tip_normal(0) = -edge2_to_tip(1);
       edge2_to_tip_normal(1) = edge2_to_tip(0);
 
       Real angle_edge1_normal = edge1_to_tip_normal * normal;
       Real angle_edge2_normal = edge2_to_tip_normal * normal;
 
       if (std::abs(angle_edge1_normal) > std::abs(angle_min) &&
           (edge1_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
       {
         edge_id_keep = i;
         distance_keep = 0.05;
         normal_keep = edge1_to_tip_normal;
         angle_min = angle_edge1_normal;
       }
       else if (std::abs(angle_edge2_normal) > std::abs(angle_min) &&
                (edge2_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
       {
         edge_id_keep = i;
         distance_keep = 0.95;
         normal_keep = edge2_to_tip_normal;
         angle_min = angle_edge2_normal;
       }
 
       if (initCutIntersectionEdge(
               crack_tip_origin, left_angle_normal, edge_ends[0], edge_ends[1], distance) &&
           (!CEMElem->isEdgePhantom(i)))
       {
         if (std::abs(left_angle_normal * normal) > std::abs(angle_min) &&
             (edge1_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
         {
           edge_id_keep = i;
           distance_keep = distance;
           normal_keep = left_angle_normal;
           angle_min = left_angle_normal * normal;
         }
       }
       else if (initCutIntersectionEdge(
                    crack_tip_origin, right_angle_normal, edge_ends[0], edge_ends[1], distance) &&
                (!CEMElem->isEdgePhantom(i)))
       {
         if (std::abs(right_angle_normal * normal) > std::abs(angle_min) &&
             (edge2_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
         {
           edge_id_keep = i;
           distance_keep = distance;
           normal_keep = right_angle_normal;
           angle_min = right_angle_normal * normal;
         }
       }
       else if (initCutIntersectionEdge(crack_tip_origin,
                                        crack_direction_normal,
                                        edge_ends[0],
                                        edge_ends[1],
                                        distance) &&
                (!CEMElem->isEdgePhantom(i)))
       {
         if (std::abs(crack_direction_normal * normal) > std::abs(angle_min) &&
             (crack_tip_direction * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
         {
           edge_id_keep = i;
           distance_keep = distance;
           normal_keep = crack_direction_normal;
           angle_min = crack_direction_normal * normal;
         }
       }
     }
   }
 
   // avoid small volume fraction cut
   if ((distance_keep - 0.05) < 0.0)
   {
     distance_keep = 0.05;
   }
   else if ((distance_keep - 0.95) > 0.0)
   {
     distance_keep = 0.95;
   }
 }

◆ cutMeshWithEFA()

bool XFEM::cutMeshWithEFA	(	const std::vector< std::shared_ptr< NonlinearSystemBase >> &	nl,
		AuxiliarySystem &	aux
	)

Definition at line 1104 of file XFEM.C.

Referenced by update().

 {
   if (nls.size() != 1)
     mooseError("XFEM does not currently support multiple nonlinear systems");
 
   std::map<unsigned int, Node *> efa_id_to_new_node;
   std::map<unsigned int, Node *> efa_id_to_new_node2;
   std::map<unsigned int, Elem *> efa_id_to_new_elem;
   _cached_solution.clear();
   _cached_aux_solution.clear();
 
   // Copy the current geometric cut element info (from last time) into the
   // _old_geom_cut_elems.
   _old_geom_cut_elems.swap(_geom_cut_elems);
   _geom_cut_elems.clear();
 
   _efa_mesh.updatePhysicalLinksAndFragments();
 
   if (_debug_output_level > 2)
   {
     _console << "\nXFEM Element fragment algorithm mesh prior to cutting:\n";
     _console << std::flush;
     _efa_mesh.printMesh();
   }
 
   _efa_mesh.updateTopology();
 
   if (_debug_output_level > 2)
   {
     _console << "\nXFEM Element fragment algorithm mesh after cutting:\n";
     _console << std::flush;
     _efa_mesh.printMesh();
   }
 
   const std::vector<EFANode *> new_nodes = _efa_mesh.getNewNodes();
   const std::vector<EFAElement *> new_elements = _efa_mesh.getChildElements();
   const std::vector<EFAElement *> delete_elements = _efa_mesh.getParentElements();
 
   bool mesh_changed = (new_nodes.size() + new_elements.size() + delete_elements.size() > 0);
 
   // Prepare to cache solution on DOFs modified by XFEM
   if (mesh_changed)
   {
     nls[0]->serializeSolution();
     aux.serializeSolution();
     if (_debug_output_level > 1)
       _console << "\n";
   }
   NumericVector<Number> & current_solution = *nls[0]->system().current_local_solution;
   NumericVector<Number> & old_solution = nls[0]->solutionOld();
   NumericVector<Number> & older_solution = nls[0]->solutionOlder();
   NumericVector<Number> & current_aux_solution = *aux.system().current_local_solution;
   NumericVector<Number> & old_aux_solution = aux.solutionOld();
   NumericVector<Number> & older_aux_solution = aux.solutionOlder();
 
   std::map<Node *, Node *> new_nodes_to_parents;
 
   // Add new nodes
   for (unsigned int i = 0; i < new_nodes.size(); ++i)
   {
     unsigned int new_node_id = new_nodes[i]->id();
     unsigned int parent_id = new_nodes[i]->parent()->id();
 
     Node * parent_node = _mesh->node_ptr(parent_id);
     Node * new_node = Node::build(*parent_node, _mesh->max_node_id()).release();
     _mesh->add_node(new_node);
 
     new_nodes_to_parents[new_node] = parent_node;
 
     new_node->set_n_systems(parent_node->n_systems());
     efa_id_to_new_node.insert(std::make_pair(new_node_id, new_node));
     if (_debug_output_level > 1)
       _console << "XFEM added new node: " << new_node->id() << std::endl;
     if (_displaced_mesh)
     {
       const Node * parent_node2 = _displaced_mesh->node_ptr(parent_id);
       Node * new_node2 = Node::build(*parent_node2, _displaced_mesh->max_node_id()).release();
       _displaced_mesh->add_node(new_node2);
 
       new_node2->set_n_systems(parent_node2->n_systems());
       efa_id_to_new_node2.insert(std::make_pair(new_node_id, new_node2));
     }
   }
 
   // Add new elements
   std::map<unsigned int, std::vector<const Elem *>> temporary_parent_children_map;
 
   std::vector<boundary_id_type> parent_boundary_ids;
 
   for (unsigned int i = 0; i < new_elements.size(); ++i)
   {
     unsigned int parent_id = new_elements[i]->getParent()->id();
     unsigned int efa_child_id = new_elements[i]->id();
 
     Elem * parent_elem = _mesh->elem_ptr(parent_id);
     Elem * libmesh_elem = Elem::build(parent_elem->type()).release();
 
     for (unsigned int m = 0; m < _geometric_cuts.size(); ++m)
     {
       for (auto & it : _sibling_elems[_geometric_cuts[m]->getInterfaceID()])
       {
         if (parent_elem == it.first)
           it.first = libmesh_elem;
         else if (parent_elem == it.second)
           it.second = libmesh_elem;
       }
     }
 
     // parent has at least two children
     if (new_elements[i]->getParent()->numChildren() > 1)
       temporary_parent_children_map[parent_elem->id()].push_back(libmesh_elem);
 
     Elem * parent_elem2 = nullptr;
     Elem * libmesh_elem2 = nullptr;
     if (_displaced_mesh)
     {
       parent_elem2 = _displaced_mesh->elem_ptr(parent_id);
       libmesh_elem2 = Elem::build(parent_elem2->type()).release();
 
       for (unsigned int m = 0; m < _geometric_cuts.size(); ++m)
       {
         for (auto & it : _sibling_displaced_elems[_geometric_cuts[m]->getInterfaceID()])
         {
           if (parent_elem2 == it.first)
             it.first = libmesh_elem2;
           else if (parent_elem2 == it.second)
             it.second = libmesh_elem2;
         }
       }
     }
 
     for (unsigned int j = 0; j < new_elements[i]->numNodes(); ++j)
     {
       unsigned int node_id = new_elements[i]->getNode(j)->id();
       Node * libmesh_node;
 
       std::map<unsigned int, Node *>::iterator nit = efa_id_to_new_node.find(node_id);
       if (nit != efa_id_to_new_node.end())
         libmesh_node = nit->second;
       else
         libmesh_node = _mesh->node_ptr(node_id);
 
       if (libmesh_node->processor_id() == DofObject::invalid_processor_id)
         libmesh_node->processor_id() = parent_elem->processor_id();
 
       libmesh_elem->set_node(j, libmesh_node);
 
       // Store solution for all nodes affected by XFEM (even existing nodes)
       if (parent_elem->is_semilocal(_mesh->processor_id()))
       {
         Node * solution_node = libmesh_node; // Node from which to store solution
         if (new_nodes_to_parents.find(libmesh_node) != new_nodes_to_parents.end())
           solution_node = new_nodes_to_parents[libmesh_node];
 
         if ((_moose_mesh->isSemiLocal(solution_node)) ||
             (libmesh_node->processor_id() == _mesh->processor_id()))
         {
           storeSolutionForNode(libmesh_node,
                                solution_node,
                                *nls[0],
                                _cached_solution,
                                current_solution,
                                old_solution,
                                older_solution);
           storeSolutionForNode(libmesh_node,
                                solution_node,
                                aux,
                                _cached_aux_solution,
                                current_aux_solution,
                                old_aux_solution,
                                older_aux_solution);
         }
       }
 
       Node * parent_node = parent_elem->node_ptr(j);
       _mesh->get_boundary_info().boundary_ids(parent_node, parent_boundary_ids);
       _mesh->get_boundary_info().add_node(libmesh_node, parent_boundary_ids);
 
       if (_displaced_mesh)
       {
         std::map<unsigned int, Node *>::iterator nit2 = efa_id_to_new_node2.find(node_id);
         if (nit2 != efa_id_to_new_node2.end())
           libmesh_node = nit2->second;
         else
           libmesh_node = _displaced_mesh->node_ptr(node_id);
 
         if (libmesh_node->processor_id() == DofObject::invalid_processor_id)
           libmesh_node->processor_id() = parent_elem2->processor_id();
 
         libmesh_elem2->set_node(j, libmesh_node);
 
         parent_node = parent_elem2->node_ptr(j);
         _displaced_mesh->get_boundary_info().boundary_ids(parent_node, parent_boundary_ids);
         _displaced_mesh->get_boundary_info().add_node(libmesh_node, parent_boundary_ids);
       }
     }
 
     libmesh_elem->set_p_level(parent_elem->p_level());
     libmesh_elem->set_p_refinement_flag(parent_elem->p_refinement_flag());
     _mesh->add_elem(libmesh_elem);
     libmesh_elem->set_n_systems(parent_elem->n_systems());
     libmesh_elem->subdomain_id() = parent_elem->subdomain_id();
     libmesh_elem->processor_id() = parent_elem->processor_id();
 
     // The crack tip origin map is stored before cut, thus the elem should be updated with new
     // element.
     std::map<const Elem *, std::vector<Point>>::iterator mit =
         _elem_crack_origin_direction_map.find(parent_elem);
     if (mit != _elem_crack_origin_direction_map.end())
     {
       std::vector<Point> crack_data = _elem_crack_origin_direction_map[parent_elem];
       _elem_crack_origin_direction_map.erase(mit);
       _elem_crack_origin_direction_map[libmesh_elem] = crack_data;
     }
 
     if (_debug_output_level > 1)
       _console << "XFEM added new element: " << libmesh_elem->id() << std::endl;
 
     XFEMCutElem * xfce = nullptr;
     if (_mesh->mesh_dimension() == 2)
     {
       EFAElement2D * new_efa_elem2d = dynamic_cast<EFAElement2D *>(new_elements[i]);
       if (!new_efa_elem2d)
         mooseError("EFAelem is not of EFAelement2D type");
       xfce = new XFEMCutElem2D(libmesh_elem,
                                new_efa_elem2d,
                                _fe_problem->assembly(0, /*nl_sys_num=*/0).qRule()->n_points(),
                                libmesh_elem->n_sides());
     }
     else if (_mesh->mesh_dimension() == 3)
     {
       EFAElement3D * new_efa_elem3d = dynamic_cast<EFAElement3D *>(new_elements[i]);
       if (!new_efa_elem3d)
         mooseError("EFAelem is not of EFAelement3D type");
       xfce = new XFEMCutElem3D(libmesh_elem,
                                new_efa_elem3d,
                                _fe_problem->assembly(0, /*nl_sys_num=*/0).qRule()->n_points(),
                                libmesh_elem->n_sides());
     }
     _cut_elem_map.insert(std::pair<unique_id_type, XFEMCutElem *>(libmesh_elem->unique_id(), xfce));
     efa_id_to_new_elem.insert(std::make_pair(efa_child_id, libmesh_elem));
 
     if (_displaced_mesh)
     {
       libmesh_elem2->set_p_level(parent_elem2->p_level());
       libmesh_elem2->set_p_refinement_flag(parent_elem2->p_refinement_flag());
       _displaced_mesh->add_elem(libmesh_elem2);
       libmesh_elem2->set_n_systems(parent_elem2->n_systems());
       libmesh_elem2->subdomain_id() = parent_elem2->subdomain_id();
       libmesh_elem2->processor_id() = parent_elem2->processor_id();
     }
 
     unsigned int n_sides = parent_elem->n_sides();
     for (unsigned int side = 0; side < n_sides; ++side)
     {
       _mesh->get_boundary_info().boundary_ids(parent_elem, side, parent_boundary_ids);
       _mesh->get_boundary_info().add_side(libmesh_elem, side, parent_boundary_ids);
     }
     if (_displaced_mesh)
     {
       n_sides = parent_elem2->n_sides();
       for (unsigned int side = 0; side < n_sides; ++side)
       {
         _displaced_mesh->get_boundary_info().boundary_ids(parent_elem2, side, parent_boundary_ids);
         _displaced_mesh->get_boundary_info().add_side(libmesh_elem2, side, parent_boundary_ids);
       }
     }
 
     unsigned int n_edges = parent_elem->n_edges();
     for (unsigned int edge = 0; edge < n_edges; ++edge)
     {
       _mesh->get_boundary_info().edge_boundary_ids(parent_elem, edge, parent_boundary_ids);
       _mesh->get_boundary_info().add_edge(libmesh_elem, edge, parent_boundary_ids);
     }
     if (_displaced_mesh)
     {
       n_edges = parent_elem2->n_edges();
       for (unsigned int edge = 0; edge < n_edges; ++edge)
       {
         _displaced_mesh->get_boundary_info().edge_boundary_ids(
             parent_elem2, edge, parent_boundary_ids);
         _displaced_mesh->get_boundary_info().add_edge(libmesh_elem2, edge, parent_boundary_ids);
       }
     }
 
     // TODO: Also need to copy neighbor material data
     if (parent_elem->processor_id() == _mesh->processor_id())
     {
       if (_material_data[0]->getMaterialPropertyStorage().hasStatefulProperties())
         _material_data[0]->copy(*libmesh_elem, *parent_elem, 0);
 
       if (_bnd_material_data[0]->getMaterialPropertyStorage().hasStatefulProperties())
         for (unsigned int side = 0; side < parent_elem->n_sides(); ++side)
         {
           _mesh->get_boundary_info().boundary_ids(parent_elem, side, parent_boundary_ids);
           std::vector<boundary_id_type>::iterator it_bd = parent_boundary_ids.begin();
           for (; it_bd != parent_boundary_ids.end(); ++it_bd)
           {
             if (_fe_problem->needBoundaryMaterialOnSide(*it_bd, 0))
               _bnd_material_data[0]->copy(*libmesh_elem, *parent_elem, side);
           }
         }
 
       // Store the current information about the geometrically cut element, and load cached material
       // properties into the new child element, if any.
       const GeometricCutUserObject * gcuo = getGeometricCutForElem(parent_elem);
       if (gcuo && gcuo->shouldHealMesh())
       {
         CutSubdomainID gcsid = getCutSubdomainID(gcuo, libmesh_elem, parent_elem);
         Xfem::CutElemInfo cei(parent_elem, gcuo, gcsid);
         _geom_cut_elems.emplace(libmesh_elem, cei);
         // Find the element to copy data from.
         // Iterate through the old geometrically cut elements, if its parent element AND the
         // geometric cut user object AND the cut subdomain ID are the same as the
         // current element, then that must be it.
         for (auto old_cei : _old_geom_cut_elems)
           if (cei.match(old_cei.second))
           {
             loadMaterialPropertiesForElement(libmesh_elem, old_cei.first, _old_geom_cut_elems);
             if (_debug_output_level > 1)
               _console << "XFEM set material properties for element: " << libmesh_elem->id()
                        << "\n";
             break;
           }
       }
 
       // Store solution for all elements affected by XFEM
       storeSolutionForElement(libmesh_elem,
                               parent_elem,
                               *nls[0],
                               _cached_solution,
                               current_solution,
                               old_solution,
                               older_solution);
       storeSolutionForElement(libmesh_elem,
                               parent_elem,
                               aux,
                               _cached_aux_solution,
                               current_aux_solution,
                               old_aux_solution,
                               older_aux_solution);
     }
   }
 
   // delete elements
   for (std::size_t i = 0; i < delete_elements.size(); ++i)
   {
     Elem * elem_to_delete = _mesh->elem_ptr(delete_elements[i]->id());
 
     // delete the XFEMCutElem object for any elements that are to be deleted
     std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
         _cut_elem_map.find(elem_to_delete->unique_id());
     if (cemit != _cut_elem_map.end())
     {
       delete cemit->second;
       _cut_elem_map.erase(cemit);
     }
 
     // remove the property storage of deleted element/side
     _material_data[0]->eraseProperty(elem_to_delete);
     _bnd_material_data[0]->eraseProperty(elem_to_delete);
 
     elem_to_delete->nullify_neighbors();
     _mesh->get_boundary_info().remove(elem_to_delete);
     unsigned int deleted_elem_id = elem_to_delete->id();
     _mesh->delete_elem(elem_to_delete);
     if (_debug_output_level > 1)
       _console << "XFEM deleted element: " << deleted_elem_id << std::endl;
 
     if (_displaced_mesh)
     {
       Elem * elem_to_delete2 = _displaced_mesh->elem_ptr(delete_elements[i]->id());
       elem_to_delete2->nullify_neighbors();
       _displaced_mesh->get_boundary_info().remove(elem_to_delete2);
       _displaced_mesh->delete_elem(elem_to_delete2);
     }
   }
 
   for (std::map<unsigned int, std::vector<const Elem *>>::iterator it =
            temporary_parent_children_map.begin();
        it != temporary_parent_children_map.end();
        ++it)
   {
     std::vector<const Elem *> & sibling_elem_vec = it->second;
     // TODO: for cut-node case, how to find the sibling elements?
     // if (sibling_elem_vec.size() != 2)
     // mooseError("Must have exactly 2 sibling elements");
 
     for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
       for (auto const & elem_id : _geom_marker_id_elems[_geometric_cuts[i]->getInterfaceID()])
         if (it->first == elem_id)
           _sibling_elems[_geometric_cuts[i]->getInterfaceID()].push_back(
               std::make_pair(sibling_elem_vec[0], sibling_elem_vec[1]));
   }
 
   // add sibling elems on displaced mesh
   if (_displaced_mesh)
   {
     for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
     {
       for (auto & se : _sibling_elems[_geometric_cuts[i]->getInterfaceID()])
       {
         Elem * elem = _displaced_mesh->elem_ptr(se.first->id());
         Elem * elem_pair = _displaced_mesh->elem_ptr(se.second->id());
         _sibling_displaced_elems[_geometric_cuts[i]->getInterfaceID()].push_back(
             std::make_pair(elem, elem_pair));
       }
     }
   }
 
   // clear the temporary map
   temporary_parent_children_map.clear();
 
   // Store information about crack tip elements
   if (mesh_changed)
   {
     _crack_tip_elems.clear();
     _crack_tip_elems_to_be_healed.clear();
     const std::set<EFAElement *> CrackTipElements = _efa_mesh.getCrackTipElements();
     std::set<EFAElement *>::const_iterator sit;
     for (sit = CrackTipElements.begin(); sit != CrackTipElements.end(); ++sit)
     {
       unsigned int eid = (*sit)->id();
       Elem * crack_tip_elem;
       std::map<unsigned int, Elem *>::iterator eit = efa_id_to_new_elem.find(eid);
       if (eit != efa_id_to_new_elem.end())
         crack_tip_elem = eit->second;
       else
         crack_tip_elem = _mesh->elem_ptr(eid);
       _crack_tip_elems.insert(crack_tip_elem);
 
       // Store the crack tip elements which are going to be healed
       for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
       {
         if (_geometric_cuts[i]->shouldHealMesh())
         {
           for (auto const & mie : _geom_marker_id_elems[_geometric_cuts[i]->getInterfaceID()])
             if ((*sit)->getParent() != nullptr)
             {
               if (_mesh->mesh_dimension() == 2)
               {
                 EFAElement2D * efa_elem2d = dynamic_cast<EFAElement2D *>((*sit)->getParent());
                 if (!efa_elem2d)
                   mooseError("EFAelem is not of EFAelement2D type");
 
                 for (unsigned int edge_id = 0; edge_id < efa_elem2d->numEdges(); ++edge_id)
                 {
                   for (unsigned int en_iter = 0; en_iter < efa_elem2d->numEdgeNeighbors(edge_id);
                        ++en_iter)
                   {
                     EFAElement2D * edge_neighbor = efa_elem2d->getEdgeNeighbor(edge_id, en_iter);
                     if (edge_neighbor != nullptr && edge_neighbor->id() == mie)
                       _crack_tip_elems_to_be_healed.insert(crack_tip_elem);
                   }
                 }
               }
               else if (_mesh->mesh_dimension() == 3)
               {
                 EFAElement3D * efa_elem3d = dynamic_cast<EFAElement3D *>((*sit)->getParent());
                 if (!efa_elem3d)
                   mooseError("EFAelem is not of EFAelement3D type");
 
                 for (unsigned int face_id = 0; face_id < efa_elem3d->numFaces(); ++face_id)
                 {
                   for (unsigned int fn_iter = 0; fn_iter < efa_elem3d->numFaceNeighbors(face_id);
                        ++fn_iter)
                   {
                     EFAElement3D * face_neighbor = efa_elem3d->getFaceNeighbor(face_id, fn_iter);
                     if (face_neighbor != nullptr && face_neighbor->id() == mie)
                       _crack_tip_elems_to_be_healed.insert(crack_tip_elem);
                   }
                 }
               }
             }
         }
       }
     }
   }
 
   if (_debug_output_level > 0)
   {
     _console << "\nXFEM mesh cutting with element fragment algorithm complete\n";
     _console << "# new nodes:        " << new_nodes.size() << "\n";
     _console << "# new elements:     " << new_elements.size() << "\n";
     _console << "# deleted elements: " << delete_elements.size() << "\n";
     _console << std::flush;
   }
 
   // store virtual nodes
   // store cut edge info
   return mesh_changed;
 }

◆ getCrackTipOrigin()

void XFEM::getCrackTipOrigin	(	std::map< unsigned int, const Elem *> &	elem_id_crack_tip,
		std::vector< Point > &	crack_front_points
	)

Definition at line 68 of file XFEM.C.

 {
   elem_id_crack_tip.clear();
   crack_front_points.clear();
   crack_front_points.resize(_elem_crack_origin_direction_map.size());
 
   unsigned int crack_tip_index = 0;
   // This map is used to sort the order in _elem_crack_origin_direction_map such that every process
   // has same order
   std::map<unsigned int, const Elem *> elem_id_map;
 
   int m = -1;
   for (std::map<const Elem *, std::vector<Point>>::iterator mit1 =
            _elem_crack_origin_direction_map.begin();
        mit1 != _elem_crack_origin_direction_map.end();
        ++mit1)
   {
     unsigned int elem_id = mit1->first->id();
     if (elem_id == std::numeric_limits<unsigned int>::max())
     {
       elem_id_map[m] = mit1->first;
       m--;
     }
     else
       elem_id_map[elem_id] = mit1->first;
   }
 
   for (std::map<unsigned int, const Elem *>::iterator mit1 = elem_id_map.begin();
        mit1 != elem_id_map.end();
        mit1++)
   {
     const Elem * elem = mit1->second;
     std::map<const Elem *, std::vector<Point>>::iterator mit2 =
         _elem_crack_origin_direction_map.find(elem);
     if (mit2 != _elem_crack_origin_direction_map.end())
     {
       elem_id_crack_tip[crack_tip_index] = mit2->first;
       crack_front_points[crack_tip_index] =
           (mit2->second)[0]; // [0] stores origin coordinates and [1] stores direction
       crack_tip_index++;
     }
   }
 }

◆ getCrackTipOriginMap()

const std::map<const Elem *, std::vector<Point> >& XFEM::getCrackTipOriginMap ( ) const

inline

Definition at line 311 of file XFEM.h.

   {
     return _elem_crack_origin_direction_map;
   }

◆ getCutPlane()

Real XFEM::getCutPlane	(	const Elem *	elem,
		const Xfem::XFEM_CUTPLANE_QUANTITY	quantity,
		unsigned int	plane_id
	)		const

Get specified component of normal or origin for cut plane for a given element.

Definition at line 1669 of file XFEM.C.

 {
   Real comp = 0.0;
   Point planedata(0.0, 0.0, 0.0);
   std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
   it = _cut_elem_map.find(elem->unique_id());
   if (it != _cut_elem_map.end())
   {
     const XFEMCutElem * xfce = it->second;
     const EFAElement * EFAelem = xfce->getEFAElement();
     if (EFAelem->isPartial()) // exclude the full crack tip elements
     {
       if ((unsigned int)quantity < 3)
       {
         unsigned int index = (unsigned int)quantity;
         planedata = xfce->getCutPlaneOrigin(plane_id, _displaced_mesh);
         comp = planedata(index);
       }
       else if ((unsigned int)quantity < 6)
       {
         unsigned int index = (unsigned int)quantity - 3;
         planedata = xfce->getCutPlaneNormal(plane_id, _displaced_mesh);
         comp = planedata(index);
       }
       else
         mooseError("In get_cut_plane index out of range");
     }
   }
   return comp;
 }

◆ getCutSubdomainID()

CutSubdomainID XFEM::getCutSubdomainID	(	const GeometricCutUserObject *	gcuo,
		const Elem *	cut_elem,
		const Elem *	parent_elem = `nullptr`
	)		const

Determine which cut subdomain the element belongs to relative to the cut.

Parameters

gcuo	The GeometricCutUserObject for the cut
cut_elem	The element being cut
parent_elem	The parent element

Definition at line 2292 of file XFEM.C.

Referenced by cutMeshWithEFA().

 {
   if (!parent_elem)
     parent_elem = cut_elem;
   // Pick any node from the parent element that is inside the physical domain and return its
   // CutSubdomainID.
   const Node * node = pickFirstPhysicalNode(cut_elem, parent_elem);
   return gcuo->getCutSubdomainID(node);
 }

◆ getEFAElem2D()

EFAElement2D * XFEM::getEFAElem2D ( const Elem * elem )

Get the EFAElement2D object for a specified libMesh element.

Parameters

elem	Pointer to the libMesh element for which the object is requested

Definition at line 1749 of file XFEM.C.

Referenced by markCutEdgesByGeometry(), and markCutEdgesByState().

 {
   EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
   EFAElement2D * EFAelem2D = dynamic_cast<EFAElement2D *>(EFAelem);
 
   if (!EFAelem2D)
     mooseError("EFAelem is not of EFAelement2D type");
 
   return EFAelem2D;
 }

◆ getEFAElem3D()

EFAElement3D * XFEM::getEFAElem3D ( const Elem * elem )

Get the EFAElement3D object for a specified libMesh element.

Parameters

elem	Pointer to the libMesh element for which the object is requested

Definition at line 1761 of file XFEM.C.

Referenced by markCutFacesByGeometry().

 {
   EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
   EFAElement3D * EFAelem3D = dynamic_cast<EFAElement3D *>(EFAelem);
 
   if (!EFAelem3D)
     mooseError("EFAelem is not of EFAelement3D type");
 
   return EFAelem3D;
 }

◆ getEFANodeCoords()

Point XFEM::getEFANodeCoords	(	EFANode *	CEMnode,
		EFAElement *	CEMElem,
		const Elem *	elem,
		MeshBase *	displaced_mesh = `nullptr`
	)		const

Definition at line 1599 of file XFEM.C.

Referenced by correctCrackExtensionDirection(), getFragmentEdges(), getFragmentFaces(), and markCutEdgesByState().

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

◆ getElementSolutionDofs()

std::vector< dof_id_type > XFEM::getElementSolutionDofs	(	const Elem *	elem,
		SystemBase &	sys
	)		const

private

Get a vector of the dof indices for all components of all variables associated with an element.

Parameters

elem	Element for which dof indices are found
sys	System for which the dof indices are found

Definition at line 2118 of file XFEM.C.

Referenced by setSolution(), and storeSolutionForElement().

 {
   SubdomainID sid = elem->subdomain_id();
   const std::vector<MooseVariableFEBase *> & vars = sys.getVariables(0);
   std::vector<dof_id_type> solution_dofs;
   solution_dofs.reserve(vars.size()); // just an approximation
   for (auto var : vars)
   {
     if (!var->isNodal())
     {
       const std::set<SubdomainID> & var_subdomains = sys.getSubdomainsForVar(var->number());
       if (var_subdomains.empty() || var_subdomains.find(sid) != var_subdomains.end())
       {
         unsigned int n_comp = elem->n_comp(sys.number(), var->number());
         for (unsigned int icomp = 0; icomp < n_comp; ++icomp)
         {
           dof_id_type elem_dof = elem->dof_number(sys.number(), var->number(), icomp);
           solution_dofs.push_back(elem_dof);
         }
       }
     }
   }
   return solution_dofs;
 }

◆ getFragmentEdges()

void XFEM::getFragmentEdges	(	const Elem *	elem,
		EFAElement2D *	CEMElem,
		std::vector< std::vector< Point >> &	frag_edges
	)		const

Definition at line 1773 of file XFEM.C.

 {
   // N.B. CEMElem here has global EFAnode
   frag_edges.clear();
   if (CEMElem->numFragments() > 0)
   {
     if (CEMElem->numFragments() > 1)
       mooseError("element ", elem->id(), " has more than one fragment at this point");
     for (unsigned int i = 0; i < CEMElem->getFragment(0)->numEdges(); ++i)
     {
       std::vector<Point> p_line(2, Point(0.0, 0.0, 0.0));
       p_line[0] = getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
       p_line[1] = getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
       frag_edges.push_back(p_line);
     }
   }
 }

◆ getFragmentFaces() [1/2]

void XFEM::getFragmentFaces	(	const Elem *	elem,
		std::vector< std::vector< Point >> &	frag_faces,
		bool	displaced_mesh = `false`
	)		const

Definition at line 1732 of file XFEM.C.

 {
   std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
   it = _cut_elem_map.find(elem->unique_id());
   if (it != _cut_elem_map.end())
   {
     const XFEMCutElem * xfce = it->second;
     if (displaced_mesh)
       xfce->getFragmentFaces(frag_faces, _displaced_mesh);
     else
       xfce->getFragmentFaces(frag_faces);
   }
 }

◆ getFragmentFaces() [2/2]

void XFEM::getFragmentFaces	(	const Elem *	elem,
		EFAElement3D *	CEMElem,
		std::vector< std::vector< Point >> &	frag_faces
	)		const

Definition at line 1794 of file XFEM.C.

 {
   // N.B. CEMElem here has global EFAnode
   frag_faces.clear();
   if (CEMElem->numFragments() > 0)
   {
     if (CEMElem->numFragments() > 1)
       mooseError("element ", elem->id(), " has more than one fragment at this point");
     for (unsigned int i = 0; i < CEMElem->getFragment(0)->numFaces(); ++i)
     {
       unsigned int num_face_nodes = CEMElem->getFragmentFace(0, i)->numNodes();
       std::vector<Point> p_line(num_face_nodes, Point(0.0, 0.0, 0.0));
       for (unsigned int j = 0; j < num_face_nodes; ++j)
         p_line[j] = getEFANodeCoords(CEMElem->getFragmentFace(0, i)->getNode(j), CEMElem, elem);
       frag_faces.push_back(p_line);
     }
   }
 }

◆ getGeometricCutForElem()

const GeometricCutUserObject * XFEM::getGeometricCutForElem ( const Elem * elem ) const

private

Get the GeometricCutUserObject associated with an element.

Parameters

elem	The element

Returns: A constant pointer to the GeometricCutUserObject, nullptr if nothing found

Definition at line 2176 of file XFEM.C.

Referenced by cutMeshWithEFA().

 {
   for (auto gcmit : _geom_marker_id_elems)
   {
     std::set<unsigned int> elems = gcmit.second;
     if (elems.find(elem->id()) != elems.end())
       return _geometric_cuts[gcmit.first];
   }
   return nullptr;
 }

◆ getGeometricCutID()

virtual unsigned int XFEM::getGeometricCutID ( const GeometricCutUserObject * gcu )

inlinevirtual

Get the interface ID corresponding to a given GeometricCutUserObject.

Parameters

gcu	pointer to the GeometricCutUserObject

Returns: the interface ID

Definition at line 274 of file XFEM.h.

   {
     return _geom_marker_id_map[gcu];
   };

◆ getNodeSolutionDofs()

std::vector< dof_id_type > XFEM::getNodeSolutionDofs	(	const Node *	node,
		SystemBase &	sys
	)		const

private

Get a vector of the dof indices for all components of all variables associated with a node.

Parameters

node	Node for which dof indices are found
sys	System for which the dof indices are found

Definition at line 2144 of file XFEM.C.

Referenced by setSolution(), and storeSolutionForNode().

 {
   const std::set<SubdomainID> & sids = _moose_mesh->getNodeBlockIds(*node);
   const std::vector<MooseVariableFEBase *> & vars = sys.getVariables(0);
   std::vector<dof_id_type> solution_dofs;
   solution_dofs.reserve(vars.size()); // just an approximation
   for (auto var : vars)
   {
     if (var->isNodal())
     {
       const std::set<SubdomainID> & var_subdomains = sys.getSubdomainsForVar(var->number());
       std::set<SubdomainID> intersect;
       set_intersection(var_subdomains.begin(),
                        var_subdomains.end(),
                        sids.begin(),
                        sids.end(),
                        std::inserter(intersect, intersect.begin()));
       if (var_subdomains.empty() || !intersect.empty())
       {
         unsigned int n_comp = node->n_comp(sys.number(), var->number());
         for (unsigned int icomp = 0; icomp < n_comp; ++icomp)
         {
           dof_id_type node_dof = node->dof_number(sys.number(), var->number(), icomp);
           solution_dofs.push_back(node_dof);
         }
       }
     }
   }
   return solution_dofs;
 }

◆ getPhysicalVolumeFraction()

Real XFEM::getPhysicalVolumeFraction ( const Elem * elem ) const

Get the volume fraction of an element that is physical.

Definition at line 1631 of file XFEM.C.

Referenced by markCutEdgesByState().

 {
   Real phys_volfrac = 1.0;
   std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
   it = _cut_elem_map.find(elem->unique_id());
   if (it != _cut_elem_map.end())
   {
     XFEMCutElem * xfce = it->second;
     const EFAElement * EFAelem = xfce->getEFAElement();
     if (EFAelem->isPartial())
     { // exclude the full crack tip elements
       xfce->computePhysicalVolumeFraction();
       phys_volfrac = xfce->getPhysicalVolumeFraction();
     }
   }
 
   return phys_volfrac;
 }

◆ getXFEMCutElemPairs()

virtual const ElementPairLocator::ElementPairList* XFEM::getXFEMCutElemPairs ( unsigned int interface_id )

inlinevirtual

Get the list of cut element pairs corresponding to a given interface ID.

Parameters

interface_id The ID of the interface

Returns: the list of elements cut by that interface

Definition at line 252 of file XFEM.h.

   {
     return &_sibling_elems[interface_id];
   }

◆ getXFEMDisplacedCutElemPairs()

virtual const ElementPairLocator::ElementPairList* XFEM::getXFEMDisplacedCutElemPairs ( unsigned int interface_id )

inlinevirtual

Get the list of cut element pairs on the displaced mesh corresponding to a given interface ID.

Parameters

interface_id The ID of the interface

Returns: the list of elements cut by that interface

Definition at line 264 of file XFEM.h.

   {
     return &_sibling_displaced_elems[interface_id];
   }

◆ getXFEMFaceWeights()

bool XFEM::getXFEMFaceWeights	(	MooseArray< Real > &	weights,
		const Elem *	elem,
		QBase *	qrule,
		const MooseArray< Point > &	q_points,
		unsigned int	side
	)

overridevirtual

Definition at line 1881 of file XFEM.C.

 {
   bool have_weights = false;
   XFEMCutElem * xfce = nullptr;
   if (isElemCut(elem, xfce))
   {
     mooseAssert(xfce != nullptr, "Must have valid XFEMCutElem object here");
     xfce->getFaceWeightMultipliers(weights, qrule, getXFEMQRule(), q_points, side);
     have_weights = true;
   }
   return have_weights;
 }

◆ getXFEMIntersectionInfo()

void XFEM::getXFEMIntersectionInfo	(	const Elem *	elem,
		unsigned int	plane_id,
		Point &	normal,
		std::vector< Point > &	intersectionPoints,
		bool	displaced_mesh = `false`
	)		const

virtual

Definition at line 1899 of file XFEM.C.

 {
   std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
   it = _cut_elem_map.find(elem->unique_id());
   if (it != _cut_elem_map.end())
   {
     const XFEMCutElem * xfce = it->second;
     if (displaced_mesh)
       xfce->getIntersectionInfo(plane_id, normal, intersectionPoints, _displaced_mesh);
     else
       xfce->getIntersectionInfo(plane_id, normal, intersectionPoints);
   }
 }

◆ getXFEMQRule()

Xfem::XFEM_QRULE & XFEM::getXFEMQRule ( )

Definition at line 1816 of file XFEM.C.

Referenced by getXFEMFaceWeights(), and getXFEMWeights().

 {
   return _XFEM_qrule;
 }

◆ getXFEMqRuleOnLine()

void XFEM::getXFEMqRuleOnLine	(	std::vector< Point > &	intersection_points,
		std::vector< Point > &	quad_pts,
		std::vector< Real > &	quad_wts
	)		const

virtual

Definition at line 1918 of file XFEM.C.

 {
   Point p1 = intersection_points[0];
   Point p2 = intersection_points[1];
 
   // number of quadrature points
   std::size_t num_qpoints = 2;
 
   // quadrature coordinates
   Real xi0 = -std::sqrt(1.0 / 3.0);
   Real xi1 = std::sqrt(1.0 / 3.0);
 
   quad_wts.resize(num_qpoints);
   quad_pts.resize(num_qpoints);
 
   Real integ_jacobian = pow((p1 - p2).norm_sq(), 0.5) * 0.5;
 
   quad_wts[0] = 1.0 * integ_jacobian;
   quad_wts[1] = 1.0 * integ_jacobian;
 
   quad_pts[0] = (1.0 - xi0) / 2.0 * p1 + (1.0 + xi0) / 2.0 * p2;
   quad_pts[1] = (1.0 - xi1) / 2.0 * p1 + (1.0 + xi1) / 2.0 * p2;
 }

◆ getXFEMqRuleOnSurface()

void XFEM::getXFEMqRuleOnSurface	(	std::vector< Point > &	intersection_points,
		std::vector< Point > &	quad_pts,
		std::vector< Real > &	quad_wts
	)		const

virtual

Definition at line 1945 of file XFEM.C.

 {
   std::size_t nnd_pe = intersection_points.size();
   Point xcrd(0.0, 0.0, 0.0);
   for (std::size_t i = 0; i < nnd_pe; ++i)
     xcrd += intersection_points[i];
   xcrd /= nnd_pe;
 
   quad_pts.resize(nnd_pe);
   quad_wts.resize(nnd_pe);
 
   Real jac = 0.0;
 
   for (std::size_t 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> subtrig_points(3, Point(0.0, 0.0, 0.0)); // sub-trig nodal coords
 
     int jplus1 = j < nnd_pe - 1 ? j + 1 : 0;
     subtrig_points[0] = xcrd;
     subtrig_points[1] = intersection_points[j];
     subtrig_points[2] = intersection_points[jplus1];
 
     std::vector<std::vector<Real>> sg2;
     Xfem::stdQuadr2D(3, 1, sg2);                 // get sg2
     for (std::size_t l = 0; l < sg2.size(); ++l) // loop all int pts on a sub-trig
     {
       Xfem::shapeFunc2D(3, sg2[l], subtrig_points, shape, jac, true); // Get shape
       std::vector<Real> tsg_line(3, 0.0);
       for (std::size_t k = 0; k < 3; ++k) // loop sub-trig nodes
       {
         tsg_line[0] += shape[k][2] * subtrig_points[k](0);
         tsg_line[1] += shape[k][2] * subtrig_points[k](1);
         tsg_line[2] += shape[k][2] * subtrig_points[k](2);
       }
       quad_pts[j + l] = Point(tsg_line[0], tsg_line[1], tsg_line[2]);
       quad_wts[j + l] = sg2[l][3] * jac;
     }
   }
 }

◆ getXFEMWeights()

bool XFEM::getXFEMWeights	(	MooseArray< Real > &	weights,
		const Elem *	elem,
		QBase *	qrule,
		const MooseArray< Point > &	q_points
	)

overridevirtual

Definition at line 1852 of file XFEM.C.

 {
   bool have_weights = false;
   XFEMCutElem * xfce = nullptr;
   if (isElemCut(elem, xfce))
   {
     mooseAssert(xfce != nullptr, "Must have valid XFEMCutElem object here");
     xfce->getWeightMultipliers(weights, qrule, getXFEMQRule(), q_points);
     have_weights = true;
 
     Real ave_weight_multiplier = 0;
     for (unsigned int i = 0; i < weights.size(); ++i)
       ave_weight_multiplier += weights[i];
     ave_weight_multiplier /= weights.size();
 
     if (ave_weight_multiplier < _min_weight_multiplier)
     {
       const Real amount_to_add = _min_weight_multiplier - ave_weight_multiplier;
       for (unsigned int i = 0; i < weights.size(); ++i)
         weights[i] += amount_to_add;
     }
   }
   return have_weights;
 }

◆ has_secondary_cut()

bool XFEM::has_secondary_cut ( )

inline

Definition at line 290 of file XFEM.h.

290 { return _has_secondary_cut; }

XFEM::_has_secondary_cut

bool _has_secondary_cut

Definition: XFEM.h:327

◆ healMesh()

bool XFEM::healMesh ( )

Potentially heal the mesh by merging some of the pairs of partial elements cut by XFEM back into single elements if indicated by the cutting objects.

Returns: true if the mesh has been modified due to healing

Definition at line 920 of file XFEM.C.

Referenced by updateHeal().

 {
   bool mesh_changed = false;
 
   std::set<Node *> nodes_to_delete;
   std::set<Node *> nodes_to_delete_displaced;
   std::set<unsigned int> cutelems_to_delete;
   unsigned int deleted_elem_count = 0;
   std::vector<std::string> healed_geometric_cuts;
 
   for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
   {
     if (_geometric_cuts[i]->shouldHealMesh())
     {
       healed_geometric_cuts.push_back(_geometric_cuts[i]->name());
       for (auto & it : _sibling_elems[_geometric_cuts[i]->getInterfaceID()])
       {
         Elem * elem1 = const_cast<Elem *>(it.first);
         Elem * elem2 = const_cast<Elem *>(it.second);
 
         std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
             _cut_elem_map.find(elem1->unique_id());
         if (cemit != _cut_elem_map.end())
         {
           const XFEMCutElem * xfce = cemit->second;
 
           cutelems_to_delete.insert(elem1->unique_id());
 
           for (unsigned int in = 0; in < elem1->n_nodes(); ++in)
           {
             Node * e1node = elem1->node_ptr(in);
             Node * e2node = elem2->node_ptr(in);
             if (!xfce->isPointPhysical(*e1node) &&
                 e1node != e2node) // This would happen at the crack tip
             {
               elem1->set_node(in, e2node);
               nodes_to_delete.insert(e1node);
             }
             else if (e1node != e2node)
               nodes_to_delete.insert(e2node);
           }
         }
         else
           mooseError("Could not find XFEMCutElem for element to be kept in healing");
 
         // Store the material properties of the elements to be healed. So that if the element is
         // immediately re-cut, we can restore the material properties (especially those stateful
         // ones).
         std::vector<const Elem *> healed_elems = {elem1, elem2};
 
         if (_geometric_cuts[i]->shouldHealMesh())
           // If the parent element will not be re-cut, then all of its nodes must have the same
           // CutSubdomainID. Therefore, just query the first node in this parent element to
           // get its CutSubdomainID.
           for (auto e : healed_elems)
             if (elem1->processor_id() == _mesh->processor_id() &&
                 e->processor_id() == _mesh->processor_id())
             {
               storeMaterialPropertiesForElement(/*parent_elem = */ elem1, /*child_elem = */ e);
               // In case the healed element is not re-cut, copy the corresponding material
               // properties to the parent element now than later.
               CutSubdomainID parent_gcsid =
                   _geometric_cuts[i]->getCutSubdomainID(elem1->node_ptr(0));
               CutSubdomainID gcsid = _geom_cut_elems[e]._cut_subdomain_id;
               if (parent_gcsid == gcsid)
                 loadMaterialPropertiesForElement(elem1, e, _geom_cut_elems);
             }
 
         if (_displaced_mesh)
         {
           Elem * elem1_displaced = _displaced_mesh->elem_ptr(it.first->id());
           Elem * elem2_displaced = _displaced_mesh->elem_ptr(it.second->id());
 
           std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
               _cut_elem_map.find(elem1_displaced->unique_id());
           if (cemit != _cut_elem_map.end())
           {
             const XFEMCutElem * xfce = cemit->second;
 
             for (unsigned int in = 0; in < elem1_displaced->n_nodes(); ++in)
             {
               Node * e1node_displaced = elem1_displaced->node_ptr(in);
               Node * e2node_displaced = elem2_displaced->node_ptr(in);
               if (!xfce->isPointPhysical(*elem1->node_ptr(in)) &&
                   e1node_displaced != e2node_displaced)
               {
                 elem1_displaced->set_node(in, e2node_displaced);
                 nodes_to_delete_displaced.insert(e1node_displaced);
               }
               else if (e1node_displaced != e2node_displaced)
                 nodes_to_delete_displaced.insert(e2node_displaced);
             }
           }
           else
             mooseError("Could not find XFEMCutElem for element to be kept in healing");
 
           elem2_displaced->nullify_neighbors();
           _displaced_mesh->get_boundary_info().remove(elem2_displaced);
           _displaced_mesh->delete_elem(elem2_displaced);
         }
 
         // remove the property storage of deleted element/side
         _material_data[0]->eraseProperty(elem2);
         _bnd_material_data[0]->eraseProperty(elem2);
 
         cutelems_to_delete.insert(elem2->unique_id());
         elem2->nullify_neighbors();
         _mesh->get_boundary_info().remove(elem2);
         unsigned int deleted_elem_id = elem2->id();
         _mesh->delete_elem(elem2);
         if (_debug_output_level > 1)
         {
           if (deleted_elem_count == 0)
             _console << "\n";
           _console << "XFEM healing deleted element: " << deleted_elem_id << std::endl;
         }
         ++deleted_elem_count;
         mesh_changed = true;
       }
     }
   }
 
   for (auto & sit : nodes_to_delete)
   {
     Node * node_to_delete = sit;
     dof_id_type deleted_node_id = node_to_delete->id();
     _mesh->get_boundary_info().remove(node_to_delete);
     _mesh->delete_node(node_to_delete);
     if (_debug_output_level > 1)
       _console << "XFEM healing deleted node: " << deleted_node_id << std::endl;
   }
 
   if (_displaced_mesh)
   {
     for (auto & sit : nodes_to_delete_displaced)
     {
       Node * node_to_delete_displaced = sit;
       _displaced_mesh->get_boundary_info().remove(node_to_delete_displaced);
       _displaced_mesh->delete_node(node_to_delete_displaced);
     }
   }
 
   for (auto & ced : cutelems_to_delete)
     if (_cut_elem_map.find(ced) != _cut_elem_map.end())
     {
       delete _cut_elem_map.find(ced)->second;
       _cut_elem_map.erase(ced);
     }
 
   for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
     if (_geometric_cuts[i]->shouldHealMesh())
       _sibling_elems[_geometric_cuts[i]->getInterfaceID()].clear();
 
   if (_displaced_mesh)
   {
     for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
       if (_geometric_cuts[i]->shouldHealMesh())
         _sibling_displaced_elems[_geometric_cuts[i]->getInterfaceID()].clear();
   }
 
   for (auto & ceh : _crack_tip_elems_to_be_healed)
   {
     _crack_tip_elems.erase(ceh);
     _elem_crack_origin_direction_map.erase(ceh);
     delete _cut_elem_map.find(ceh->unique_id())->second;
     _cut_elem_map.erase(ceh->unique_id());
   }
 
   if (!healed_geometric_cuts.empty() && _debug_output_level > 0)
   {
     _console << "\nXFEM mesh healing complete\n";
     _console << "Names of healed geometric cut objects: ";
     for (auto geomcut : healed_geometric_cuts)
       _console << geomcut << " ";
     _console << "\n";
     _console << "# deleted nodes:    " << nodes_to_delete.size() << "\n";
     _console << "# deleted elements: " << deleted_elem_count << "\n";
     _console << std::flush;
   }
 
   return mesh_changed;
 }

◆ initCutIntersectionEdge()

bool XFEM::initCutIntersectionEdge	(	Point	cut_origin,
		RealVectorValue	cut_normal,
		Point &	edge_p1,
		Point &	edge_p2,
		Real &	dist
	)

Definition at line 901 of file XFEM.C.

Referenced by correctCrackExtensionDirection(), and markCutEdgesByState().

 {
   dist = 0.0;
   bool does_intersect = false;
   Point origin2p1 = edge_p1 - cut_origin;
   Real plane2p1 = cut_normal(0) * origin2p1(0) + cut_normal(1) * origin2p1(1);
   Point origin2p2 = edge_p2 - cut_origin;
   Real plane2p2 = cut_normal(0) * origin2p2(0) + cut_normal(1) * origin2p2(1);
 
   if (plane2p1 * plane2p2 < 0.0)
   {
     dist = -plane2p1 / (plane2p2 - plane2p1);
     does_intersect = true;
   }
   return does_intersect;
 }

◆ initSolution()

void XFEM::initSolution	(	const std::vector< std::shared_ptr< NonlinearSystemBase >> &	nl,
		AuxiliarySystem &	aux
	)

overridevirtual

Implements XFEMInterface.

Definition at line 308 of file XFEM.C.

 {
   if (nls.size() != 1)
     mooseError("XFEM does not currently support multiple nonlinear systems");
 
   nls[0]->serializeSolution();
   aux.serializeSolution();
   NumericVector<Number> & current_solution = *nls[0]->system().current_local_solution;
   NumericVector<Number> & old_solution = nls[0]->solutionOld();
   NumericVector<Number> & older_solution = nls[0]->solutionOlder();
   NumericVector<Number> & current_aux_solution = *aux.system().current_local_solution;
   NumericVector<Number> & old_aux_solution = aux.solutionOld();
   NumericVector<Number> & older_aux_solution = aux.solutionOlder();
 
   setSolution(*nls[0], _cached_solution, current_solution, old_solution, older_solution);
   setSolution(
       aux, _cached_aux_solution, current_aux_solution, old_aux_solution, older_aux_solution);
 
   current_solution.close();
   old_solution.close();
   older_solution.close();
   current_aux_solution.close();
   old_aux_solution.close();
   older_aux_solution.close();
 
   _cached_solution.clear();
   _cached_aux_solution.clear();
 }

◆ isElemAtCrackTip()

bool XFEM::isElemAtCrackTip ( const Elem * elem ) const

Definition at line 1703 of file XFEM.C.

Referenced by markCutEdgesByGeometry(), and markCutEdgesByState().

 {
   return (_crack_tip_elems.find(elem) != _crack_tip_elems.end());
 }

◆ isElemCut() [1/2]

bool XFEM::isElemCut	(	const Elem *	elem,
		XFEMCutElem *&	xfce
	)		const

Definition at line 1709 of file XFEM.C.

Referenced by getXFEMFaceWeights(), getXFEMWeights(), and isElemCut().

 {
   const auto it = _cut_elem_map.find(elem->unique_id());
   if (it != _cut_elem_map.end())
   {
     xfce = it->second;
     const EFAElement * EFAelem = xfce->getEFAElement();
     if (EFAelem->isPartial()) // exclude the full crack tip elements
       return true;
   }
 
   xfce = nullptr;
   return false;
 }

◆ isElemCut() [2/2]

bool XFEM::isElemCut ( const Elem * elem ) const

Definition at line 1725 of file XFEM.C.

 {
   XFEMCutElem * xfce;
   return isElemCut(elem, xfce);
 }

◆ isPointInsidePhysicalDomain()

bool XFEM::isPointInsidePhysicalDomain	(	const Elem *	elem,
		const Point &	point
	)		const

Return true if the point is inside the element physical domain Note: if this element is not cut, return true too.

Definition at line 1651 of file XFEM.C.

Referenced by pickFirstPhysicalNode().

 {
   std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
   it = _cut_elem_map.find(elem->unique_id());
   if (it != _cut_elem_map.end())
   {
     XFEMCutElem * xfce = it->second;
 
     if (xfce->isPointPhysical(point))
       return true;
   }
   else
     return true;
 
   return false;
 }

◆ loadMaterialPropertiesForElement()

void XFEM::loadMaterialPropertiesForElement	(	const Elem *	elem,
		const Elem *	elem_from,
		std::unordered_map< const Elem *, Xfem::CutElemInfo > &	cached_cei
	)		const

private

Helper function to store the material properties of a healed element.

Parameters

elem	The cut element to restore material properties to.
elem_from	The element to copy material properties from.
cached_cei	The material properties cache to use.

Definition at line 2258 of file XFEM.C.

Referenced by cutMeshWithEFA(), and healMesh().

 {
   // Restore the element material properties
   mooseAssert(cached_cei.count(elem_from) > 0, "XFEM: Unable to find cached material properties.");
   Xfem::CutElemInfo & cei = cached_cei[elem_from];
 
   // Load element material properties from cached properties
   loadMaterialPropertiesForElementHelper(elem,
                                          cei._elem_material_properties,
                                          _material_data[0]->getMaterialPropertyStorageForXFEM({}));
 
   // Check if any of the element side need material properties
   bool need_boundary_materials = false;
   for (unsigned int side = 0; side < elem->n_sides(); ++side)
   {
     std::vector<boundary_id_type> elem_boundary_ids;
     _mesh->get_boundary_info().boundary_ids(elem, side, elem_boundary_ids);
     for (auto bdid : elem_boundary_ids)
       if (_fe_problem->needBoundaryMaterialOnSide(bdid, 0))
         need_boundary_materials = true;
   }
 
   // Load boundary material properties from cached properties
   if (need_boundary_materials)
     loadMaterialPropertiesForElementHelper(
         elem,
         cei._bnd_material_properties,
         _bnd_material_data[0]->getMaterialPropertyStorageForXFEM({}));
 }

◆ loadMaterialPropertiesForElementHelper()

void XFEM::loadMaterialPropertiesForElementHelper	(	const Elem *	elem,
		const Xfem::CachedMaterialProperties &	cached_props,
		MaterialPropertyStorage &	storage
	)		const

private

Load the material properties.

Parameters

props_deserialized	The material properties
props_serialized	The serialized material properties

This does very dirty things and writes back to MOOSE's stateful properties. It should not do this in the future.

Definition at line 2234 of file XFEM.C.

Referenced by loadMaterialPropertiesForElement().

 {
   if (!storage.hasStatefulProperties())
     return;
 
   for (const auto state : storage.statefulIndexRange())
   {
     const auto & serialized_props = cached_props[state - 1];
     for (const auto & [side, serialized_side_props] : serialized_props)
     {
       std::istringstream iss;
       iss.str(serialized_side_props);
       iss.clear();
 
       // This is very dirty. We should not write to MOOSE's stateful properties.
       // Please remove me :(
       dataLoad(iss, storage.setProps(elem, side, state), nullptr);
     }
   }
 }

◆ markCutEdgesByGeometry()

bool XFEM::markCutEdgesByGeometry ( )

Definition at line 395 of file XFEM.C.

Referenced by markCuts().

 {
   bool marked_edges = false;
   bool marked_nodes = false;
 
   for (const auto & gme : _geom_marked_elems_2d)
   {
     for (const auto & gmei : gme.second)
     {
       EFAElement2D * EFAElem = getEFAElem2D(gme.first);
 
       for (unsigned int i = 0; i < gmei._elem_cut_edges.size(); ++i) // mark element edges
       {
         if (!EFAElem->isEdgePhantom(
                 gmei._elem_cut_edges[i]._host_side_id)) // must not be phantom edge
         {
           _efa_mesh.addElemEdgeIntersection(gme.first->id(),
                                             gmei._elem_cut_edges[i]._host_side_id,
                                             gmei._elem_cut_edges[i]._distance);
           marked_edges = true;
         }
       }
 
       for (unsigned int i = 0; i < gmei._elem_cut_nodes.size(); ++i) // mark element edges
       {
         _efa_mesh.addElemNodeIntersection(gme.first->id(), gmei._elem_cut_nodes[i]._host_id);
         marked_nodes = true;
       }
 
       for (unsigned int i = 0; i < gmei._frag_cut_edges.size();
            ++i) // MUST DO THIS AFTER MARKING ELEMENT EDGES
       {
         if (!EFAElem->getFragment(0)->isSecondaryInteriorEdge(
                 gmei._frag_cut_edges[i]._host_side_id))
         {
           if (_efa_mesh.addFragEdgeIntersection(gme.first->id(),
                                                 gmei._frag_cut_edges[i]._host_side_id,
                                                 gmei._frag_cut_edges[i]._distance))
           {
             marked_edges = true;
             if (!isElemAtCrackTip(gme.first))
               _has_secondary_cut = true;
           }
         }
       }
     }
   }
 
   return marked_edges || marked_nodes;
 }

◆ markCutEdgesByState()

bool XFEM::markCutEdgesByState ( Real time )

Definition at line 588 of file XFEM.C.

Referenced by markCuts().

 {
   bool marked_edges = false;
   for (std::map<const Elem *, RealVectorValue>::iterator pmeit = _state_marked_elems.begin();
        pmeit != _state_marked_elems.end();
        ++pmeit)
   {
     const Elem * elem = pmeit->first;
     RealVectorValue normal = pmeit->second;
     EFAElement2D * CEMElem = getEFAElem2D(elem);
 
     Real volfrac_elem = getPhysicalVolumeFraction(elem);
     if (volfrac_elem < 0.25)
       continue;
 
     // continue if elem is already cut twice - IMPORTANT
     if (CEMElem->isFinalCut())
       continue;
 
     // find the first cut edge
     unsigned int nsides = CEMElem->numEdges();
     unsigned int orig_cut_side_id = std::numeric_limits<unsigned int>::max();
     Real orig_cut_distance = -1.0;
     EFANode * orig_node = nullptr;
     EFAEdge * orig_edge = nullptr;
 
     // crack tip origin coordinates and direction
     Point crack_tip_origin(0, 0, 0);
     Point crack_tip_direction(0, 0, 0);
 
     if (isElemAtCrackTip(elem)) // crack tip element's crack intiation
     {
       orig_cut_side_id = CEMElem->getTipEdgeID();
       if (orig_cut_side_id < nsides) // valid crack-tip edge found
       {
         orig_edge = CEMElem->getEdge(orig_cut_side_id);
         orig_node = CEMElem->getTipEmbeddedNode();
       }
       else
         mooseError("element ", elem->id(), " has no valid crack-tip edge");
 
       // obtain the crack tip origin coordinates and direction.
       std::map<const Elem *, std::vector<Point>>::iterator ecodm =
           _elem_crack_origin_direction_map.find(elem);
       if (ecodm != _elem_crack_origin_direction_map.end())
       {
         crack_tip_origin = (ecodm->second)[0];
         crack_tip_direction = (ecodm->second)[1];
       }
       else
         mooseError("element ", elem->id(), " cannot find its crack tip origin and direction.");
     }
     else
     {
       std::map<const Elem *, unsigned int>::iterator mit1;
       mit1 = _state_marked_elem_sides.find(elem);
       std::set<const Elem *>::iterator mit2;
       mit2 = _state_marked_frags.find(elem);
 
       if (mit1 != _state_marked_elem_sides.end()) // specified boundary crack initiation
       {
         orig_cut_side_id = mit1->second;
         if (!CEMElem->isEdgePhantom(orig_cut_side_id) &&
             !CEMElem->getEdge(orig_cut_side_id)->hasIntersection())
         {
           orig_cut_distance = 0.5;
           _efa_mesh.addElemEdgeIntersection(elem->id(), orig_cut_side_id, orig_cut_distance);
           orig_edge = CEMElem->getEdge(orig_cut_side_id);
           orig_node = orig_edge->getEmbeddedNode(0);
           // get a virtual crack tip direction
           Point elem_center(0.0, 0.0, 0.0);
           Point edge_center;
           for (unsigned int i = 0; i < nsides; ++i)
           {
             elem_center += getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
             elem_center += getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
           }
           elem_center /= nsides * 2.0;
           edge_center = getEFANodeCoords(orig_edge->getNode(0), CEMElem, elem) +
                         getEFANodeCoords(orig_edge->getNode(1), CEMElem, elem);
           edge_center /= 2.0;
           crack_tip_origin = edge_center;
           crack_tip_direction = elem_center - edge_center;
           crack_tip_direction /= pow(crack_tip_direction.norm_sq(), 0.5);
         }
         else
           continue; // skip this elem if specified boundary edge is phantom
       }
       else if (mit2 != _state_marked_frags.end()) // cut-surface secondary crack initiation
       {
         if (CEMElem->numFragments() != 1)
           mooseError("element ",
                      elem->id(),
                      " flagged for a secondary crack, but has ",
                      CEMElem->numFragments(),
                      " fragments");
         std::vector<unsigned int> interior_edge_id = CEMElem->getFragment(0)->getInteriorEdgeID();
         if (interior_edge_id.size() == 1)
           orig_cut_side_id = interior_edge_id[0];
         else
           continue; // skip this elem if more than one interior edges found (i.e. elem's been cut
                     // twice)
         orig_cut_distance = 0.5;
         _efa_mesh.addFragEdgeIntersection(elem->id(), orig_cut_side_id, orig_cut_distance);
         orig_edge = CEMElem->getFragmentEdge(0, orig_cut_side_id);
         orig_node = orig_edge->getEmbeddedNode(0); // must be an interior embedded node
         Point elem_center(0.0, 0.0, 0.0);
         Point edge_center;
         unsigned int nsides_frag = CEMElem->getFragment(0)->numEdges();
         for (unsigned int i = 0; i < nsides_frag; ++i)
         {
           elem_center +=
               getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
           elem_center +=
               getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
         }
         elem_center /= nsides_frag * 2.0;
         edge_center = getEFANodeCoords(orig_edge->getNode(0), CEMElem, elem) +
                       getEFANodeCoords(orig_edge->getNode(1), CEMElem, elem);
         edge_center /= 2.0;
         crack_tip_origin = edge_center;
         crack_tip_direction = elem_center - edge_center;
         crack_tip_direction /= pow(crack_tip_direction.norm_sq(), 0.5);
       }
       else
         mooseError("element ",
                    elem->id(),
                    " flagged for state-based growth, but has no edge intersections");
     }
 
     Point cut_origin(0.0, 0.0, 0.0);
     if (orig_node)
       cut_origin = getEFANodeCoords(orig_node, CEMElem, elem); // cutting plane origin's coords
     else
       mooseError("element ", elem->id(), " does not have valid orig_node");
 
     // loop through element edges to add possible second cut points
     std::vector<Point> edge_ends(2, Point(0.0, 0.0, 0.0));
     Point edge1(0.0, 0.0, 0.0);
     Point edge2(0.0, 0.0, 0.0);
     Point cut_edge_point(0.0, 0.0, 0.0);
     bool find_compatible_direction = false;
     unsigned int edge_id_keep = 0;
     Real distance_keep = 0.0;
     Point normal_keep(0.0, 0.0, 0.0);
     Real distance = 0.0;
     bool edge_cut = false;
 
     for (unsigned int i = 0; i < nsides; ++i)
     {
       if (!orig_edge->isPartialOverlap(*CEMElem->getEdge(i)))
       {
         edge_ends[0] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
         edge_ends[1] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
         if ((initCutIntersectionEdge(
                  crack_tip_origin, normal, edge_ends[0], edge_ends[1], distance) &&
              (!CEMElem->isEdgePhantom(i))))
         {
           cut_edge_point = distance * edge_ends[1] + (1.0 - distance) * edge_ends[0];
           distance_keep = distance;
           edge_id_keep = i;
           normal_keep = normal;
           edge_cut = true;
           break;
         }
       }
     }
 
     Point between_two_cuts = (cut_edge_point - crack_tip_origin);
     between_two_cuts /= pow(between_two_cuts.norm_sq(), 0.5);
     Real angle_between_two_cuts = between_two_cuts * crack_tip_direction;
 
     if (angle_between_two_cuts > std::cos(45.0 / 180.0 * 3.14159)) // original cut direction is good
       find_compatible_direction = true;
 
     if (!find_compatible_direction && edge_cut)
       correctCrackExtensionDirection(elem,
                                      CEMElem,
                                      orig_edge,
                                      normal,
                                      crack_tip_origin,
                                      crack_tip_direction,
                                      distance_keep,
                                      edge_id_keep,
                                      normal_keep);
 
     if (edge_cut)
     {
       if (!_use_crack_growth_increment)
         _efa_mesh.addElemEdgeIntersection(elem->id(), edge_id_keep, distance_keep);
       else
       {
         Point growth_direction(0.0, 0.0, 0.0);
 
         growth_direction(0) = -normal_keep(1);
         growth_direction(1) = normal_keep(0);
 
         if (growth_direction * crack_tip_direction < 1.0e-10)
           growth_direction *= (-1.0);
 
         Real x0 = crack_tip_origin(0);
         Real y0 = crack_tip_origin(1);
         Real x1 = x0 + _crack_growth_increment * growth_direction(0);
         Real y1 = y0 + _crack_growth_increment * growth_direction(1);
 
         XFEMCrackGrowthIncrement2DCut geometric_cut(x0, y0, x1, y1, time * 0.9, time * 0.9);
 
         for (const auto & elem : _mesh->element_ptr_range())
         {
           std::vector<CutEdgeForCrackGrowthIncr> elem_cut_edges;
           EFAElement2D * CEMElem = getEFAElem2D(elem);
 
           // continue if elem has been already cut twice - IMPORTANT
           if (CEMElem->isFinalCut())
             continue;
 
           // mark cut edges for the element and its fragment
           geometric_cut.cutElementByCrackGrowthIncrement(elem, elem_cut_edges, time);
 
           for (unsigned int i = 0; i < elem_cut_edges.size(); ++i) // mark element edges
           {
             if (!CEMElem->isEdgePhantom(
                     elem_cut_edges[i]._host_side_id)) // must not be phantom edge
             {
               _efa_mesh.addElemEdgeIntersection(
                   elem->id(), elem_cut_edges[i]._host_side_id, elem_cut_edges[i]._distance);
             }
           }
         }
       }
     }
     // loop though framgent boundary edges to add possible second cut points
     // N.B. must do this after marking element edges
     if (CEMElem->numFragments() > 0 && !edge_cut)
     {
       for (unsigned int i = 0; i < CEMElem->getFragment(0)->numEdges(); ++i)
       {
         if (!orig_edge->isPartialOverlap(*CEMElem->getFragmentEdge(0, i)))
         {
           edge_ends[0] =
               getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
           edge_ends[1] =
               getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
           if (initCutIntersectionEdge(
                   crack_tip_origin, normal, edge_ends[0], edge_ends[1], distance) &&
               (!CEMElem->getFragment(0)->isSecondaryInteriorEdge(i)))
           {
             if (_efa_mesh.addFragEdgeIntersection(elem->id(), edge_id_keep, distance_keep))
               if (!isElemAtCrackTip(elem))
                 _has_secondary_cut = true;
             break;
           }
         }
       }
     }
 
     marked_edges = true;
 
   } // loop over all state_marked_elems
 
   return marked_edges;
 }

◆ markCutFacesByGeometry()

bool XFEM::markCutFacesByGeometry ( )

Definition at line 852 of file XFEM.C.

Referenced by markCuts().

 {
   bool marked_faces = false;
 
   for (const auto & gme : _geom_marked_elems_3d)
   {
     for (const auto & gmei : gme.second)
     {
       EFAElement3D * EFAElem = getEFAElem3D(gme.first);
 
       for (unsigned int i = 0; i < gmei._elem_cut_faces.size(); ++i) // mark element faces
       {
         if (!EFAElem->isFacePhantom(gmei._elem_cut_faces[i]._face_id)) // must not be phantom face
         {
           _efa_mesh.addElemFaceIntersection(gme.first->id(),
                                             gmei._elem_cut_faces[i]._face_id,
                                             gmei._elem_cut_faces[i]._face_edge,
                                             gmei._elem_cut_faces[i]._position);
           marked_faces = true;
         }
       }
 
       for (unsigned int i = 0; i < gmei._frag_cut_faces.size();
            ++i) // MUST DO THIS AFTER MARKING ELEMENT EDGES
       {
         if (!EFAElem->getFragment(0)->isThirdInteriorFace(gmei._frag_cut_faces[i]._face_id))
         {
           _efa_mesh.addFragFaceIntersection(gme.first->id(),
                                             gmei._frag_cut_faces[i]._face_id,
                                             gmei._frag_cut_faces[i]._face_edge,
                                             gmei._frag_cut_faces[i]._position);
           marked_faces = true;
         }
       }
     }
   }
 
   return marked_faces;
 }

◆ markCutFacesByState()

bool XFEM::markCutFacesByState ( )

Definition at line 893 of file XFEM.C.

Referenced by markCuts().

 {
   bool marked_faces = false;
   // TODO: need to finish this for 3D problems
   return marked_faces;
 }

◆ markCuts()

bool XFEM::markCuts ( Real time )

Definition at line 378 of file XFEM.C.

Referenced by update().

 {
   bool marked_sides = false;
   if (_mesh->mesh_dimension() == 2)
   {
     marked_sides = markCutEdgesByGeometry();
     marked_sides |= markCutEdgesByState(time);
   }
   else if (_mesh->mesh_dimension() == 3)
   {
     marked_sides = markCutFacesByGeometry();
     marked_sides |= markCutFacesByState();
   }
   return marked_sides;
 }

◆ pickFirstPhysicalNode()

const Node * XFEM::pickFirstPhysicalNode	(	const Elem *	e,
		const Elem *	e0
	)		const

private

Return the first node in the provided element that is found to be in the physical domain.

Parameters

e	Constant pointer to the child element
e0	Constant pointer to the parent element whose nodes will be querried

Returns: A constant pointer to the node

Definition at line 2305 of file XFEM.C.

Referenced by getCutSubdomainID().

 {
   for (auto i : e0->node_index_range())
     if (isPointInsidePhysicalDomain(e, e0->node_ref(i)))
       return e0->node_ptr(i);
   mooseError("cannot find a physical node in the current element");
   return nullptr;
 }

◆ setCrackGrowthMethod()

void XFEM::setCrackGrowthMethod	(	bool	use_crack_growth_increment,
		Real	crack_growth_increment
	)

Definition at line 1833 of file XFEM.C.

 {
   _use_crack_growth_increment = use_crack_growth_increment;
   _crack_growth_increment = crack_growth_increment;
 }

◆ setDebugOutputLevel()

void XFEM::setDebugOutputLevel ( unsigned int debug_output_level )

Controls amount of debugging information output.

Parameters

debug_output_level How much information to output (see description of _debug_output_level)

Definition at line 1840 of file XFEM.C.

 {
   _debug_output_level = debug_output_level;
 }

◆ setMinWeightMultiplier()

void XFEM::setMinWeightMultiplier ( Real min_weight_multiplier )

Controls the minimum average weight multiplier for each element.

Parameters

min_weight_multiplier The minimum average weight multiplier applied by XFEM to the standard quadrature weights

Definition at line 1846 of file XFEM.C.

 {
   _min_weight_multiplier = min_weight_multiplier;
 }

◆ setSolution()

void XFEM::setSolution	(	SystemBase &	sys,
		const std::map< unique_id_type, std::vector< Real >> &	stored_solution,
		NumericVector< Number > &	current_solution,
		NumericVector< Number > &	old_solution,
		NumericVector< Number > &	older_solution
	)

private

Set the solution for all locally-owned nodes/elements that have stored values.

Parameters

sys	System for which the solution is set
stored_solution	Data structure that the stored solution is obtained from
current_solution	Current solution vector that will be set
old_solution	Old solution vector that will be set
older_solution	Older solution vector that will be set

Definition at line 2057 of file XFEM.C.

Referenced by initSolution().

 {
   for (auto & node : _mesh->local_node_ptr_range())
   {
     auto mit = stored_solution.find(node->unique_id());
     if (mit != stored_solution.end())
     {
       const std::vector<Real> & stored_node_solution = mit->second;
       std::vector<dof_id_type> stored_solution_dofs = getNodeSolutionDofs(node, sys);
       setSolutionForDOFs(stored_node_solution,
                          stored_solution_dofs,
                          current_solution,
                          old_solution,
                          older_solution);
     }
   }
 
   for (auto & elem : as_range(_mesh->local_elements_begin(), _mesh->local_elements_end()))
   {
     auto mit = stored_solution.find(elem->unique_id());
     if (mit != stored_solution.end())
     {
       const std::vector<Real> & stored_elem_solution = mit->second;
       std::vector<dof_id_type> stored_solution_dofs = getElementSolutionDofs(elem, sys);
       setSolutionForDOFs(stored_elem_solution,
                          stored_solution_dofs,
                          current_solution,
                          old_solution,
                          older_solution);
     }
   }
 }

◆ setSolutionForDOFs()

void XFEM::setSolutionForDOFs	(	const std::vector< Real > &	stored_solution,
		const std::vector< dof_id_type > &	stored_solution_dofs,
		NumericVector< Number > &	current_solution,
		NumericVector< Number > &	old_solution,
		NumericVector< Number > &	older_solution
	)

private

Set the solution for a set of DOFs.

Parameters

stored_solution	Stored solution values to set the solution to
stored_solution_dofs	Dof indices for the entries in stored_solution
current_solution	Current solution vector that will be set
old_solution	Old solution vector that will be set
older_solution	Older solution vector that will be set

Definition at line 2095 of file XFEM.C.

Referenced by setSolution().

 {
   // Solution vector is stored first for current, then old and older solutions.
   // These are the offsets to the beginning of the old and older solutions in the vector.
   const auto old_solution_offset = stored_solution_dofs.size();
   const auto older_solution_offset = old_solution_offset * 2;
 
   for (std::size_t i = 0; i < stored_solution_dofs.size(); ++i)
   {
     current_solution.set(stored_solution_dofs[i], stored_solution[i]);
     if (_fe_problem->isTransient())
     {
       old_solution.set(stored_solution_dofs[i], stored_solution[old_solution_offset + i]);
       older_solution.set(stored_solution_dofs[i], stored_solution[older_solution_offset + i]);
     }
   }
 }

◆ setXFEMQRule()

void XFEM::setXFEMQRule ( std::string & xfem_qrule )

Definition at line 1822 of file XFEM.C.

 {
   if (xfem_qrule == "volfrac")
     _XFEM_qrule = Xfem::VOLFRAC;
   else if (xfem_qrule == "moment_fitting")
     _XFEM_qrule = Xfem::MOMENT_FITTING;
   else if (xfem_qrule == "direct")
     _XFEM_qrule = Xfem::DIRECT;
 }

◆ storeCrackTipOriginAndDirection()

void XFEM::storeCrackTipOriginAndDirection ( )

Definition at line 187 of file XFEM.C.

Referenced by update().

 {
   _elem_crack_origin_direction_map.clear();
   std::set<EFAElement *> CrackTipElements = _efa_mesh.getCrackTipElements();
   std::set<EFAElement *>::iterator sit;
   for (sit = CrackTipElements.begin(); sit != CrackTipElements.end(); ++sit)
   {
     if (_mesh->mesh_dimension() == 2)
     {
       EFAElement2D * CEMElem = dynamic_cast<EFAElement2D *>(*sit);
       EFANode * tip_node = CEMElem->getTipEmbeddedNode();
       unsigned int cts_id = CEMElem->getCrackTipSplitElementID();
 
       Point origin(0, 0, 0);
       Point direction(0, 0, 0);
 
       std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
       it = _cut_elem_map.find(_mesh->elem_ptr(cts_id)->unique_id());
       if (it != _cut_elem_map.end())
       {
         const XFEMCutElem * xfce = it->second;
         const EFAElement * EFAelem = xfce->getEFAElement();
         if (EFAelem->isPartial()) // exclude the full crack tip elements
         {
           xfce->getCrackTipOriginAndDirection(tip_node->id(), origin, direction);
         }
       }
 
       std::vector<Point> tip_data;
       tip_data.push_back(origin);
       tip_data.push_back(direction);
       const Elem * elem = _mesh->elem_ptr((*sit)->id());
       _elem_crack_origin_direction_map.insert(
           std::pair<const Elem *, std::vector<Point>>(elem, tip_data));
     }
   }
 }

◆ storeMaterialPropertiesForElement()

void XFEM::storeMaterialPropertiesForElement	(	const Elem *	parent_elem,
		const Elem *	child_elem
	)

private

Helper function to store the material properties of a healed element.

Parameters

parent_elem	The parent element
elem1	The first child element
elem2	The second child element

Definition at line 2206 of file XFEM.C.

Referenced by healMesh().

 {
   // Set the parent element so that it is consistent post-healing
   _geom_cut_elems[child_elem]._parent_elem = parent_elem;
 
   // Locally store the element material properties
   storeMaterialPropertiesForElementHelper(child_elem,
                                           _material_data[0]->getMaterialPropertyStorageForXFEM({}));
 
   // Locally store the boundary material properties
   // First check if any of the side need material properties
   bool need_boundary_materials = false;
   for (unsigned int side = 0; side < child_elem->n_sides(); ++side)
   {
     std::vector<boundary_id_type> elem_boundary_ids;
     _mesh->get_boundary_info().boundary_ids(child_elem, side, elem_boundary_ids);
     for (auto bdid : elem_boundary_ids)
       if (_fe_problem->needBoundaryMaterialOnSide(bdid, 0))
         need_boundary_materials = true;
   }
 
   // If boundary material properties are needed for this element, then store them.
   if (need_boundary_materials)
     storeMaterialPropertiesForElementHelper(
         child_elem, _bnd_material_data[0]->getMaterialPropertyStorageForXFEM({}));
 }

◆ storeMaterialPropertiesForElementHelper()

void XFEM::storeMaterialPropertiesForElementHelper	(	const Elem *	elem,
		MaterialPropertyStorage &	storage
	)

private

Definition at line 2188 of file XFEM.C.

Referenced by storeMaterialPropertiesForElement().

 {
   for (const auto state : storage.statefulIndexRange())
   {
     const auto & elem_props = storage.props(state).at(elem);
     auto & serialized_props = _geom_cut_elems[elem]._elem_material_properties[state - 1];
     serialized_props.clear();
     for (const auto & side_props_pair : elem_props)
     {
       const auto side = side_props_pair.first;
       std::ostringstream oss;
       dataStore(oss, storage.setProps(elem, side, state), nullptr);
       serialized_props[side].assign(oss.str());
     }
   }
 }

◆ storeSolutionForElement()

void XFEM::storeSolutionForElement	(	const Elem *	elem_to_store_to,
		const Elem *	elem_to_store_from,
		SystemBase &	sys,
		std::map< unique_id_type, std::vector< Real >> &	stored_solution,
		const NumericVector< Number > &	current_solution,
		const NumericVector< Number > &	old_solution,
		const NumericVector< Number > &	older_solution
	)

private

Store the solution in stored_solution for a given element.

Parameters

elem_to_store_to	Element for which the solution will be stored
elem_to_store_from	Element from which the solution to be stored is obtained
sys	System from which the solution is stored
stored_solution	Data structure that the stored solution is saved to
current_solution	Current solution vector that the solution is obtained from
old_solution	Old solution vector that the solution is obtained from
older_solution	Older solution vector that the solution is obtained from

Definition at line 2023 of file XFEM.C.

Referenced by cutMeshWithEFA().

 {
   std::vector<dof_id_type> stored_solution_dofs = getElementSolutionDofs(elem_to_store_from, sys);
   std::vector<Real> stored_solution_scratch;
   // Size for current solution, as well as for old, and older solution only for transient case
   std::size_t stored_solution_size =
       (_fe_problem->isTransient() ? stored_solution_dofs.size() * 3 : stored_solution_dofs.size());
   stored_solution_scratch.reserve(stored_solution_size);
 
   // Store in the order defined in stored_solution_dofs first for the current, then for old and
   // older if applicable
   for (auto dof : stored_solution_dofs)
     stored_solution_scratch.push_back(current_solution(dof));
 
   if (_fe_problem->isTransient())
   {
     for (auto dof : stored_solution_dofs)
       stored_solution_scratch.push_back(old_solution(dof));
 
     for (auto dof : stored_solution_dofs)
       stored_solution_scratch.push_back(older_solution(dof));
   }
 
   if (stored_solution_scratch.size() > 0)
     stored_solution[elem_to_store_to->unique_id()] = stored_solution_scratch;
 }

◆ storeSolutionForNode()

void XFEM::storeSolutionForNode	(	const Node *	node_to_store_to,
		const Node *	node_to_store_from,
		SystemBase &	sys,
		std::map< unique_id_type, std::vector< Real >> &	stored_solution,
		const NumericVector< Number > &	current_solution,
		const NumericVector< Number > &	old_solution,
		const NumericVector< Number > &	older_solution
	)

private

Store the solution in stored_solution for a given node.

Parameters

node_to_store_to	Node for which the solution will be stored
node_to_store_from	Node from which the solution to be stored is obtained
sys	System from which the solution is stored
stored_solution	Data structure that the stored solution is saved to
current_solution	Current solution vector that the solution is obtained from
old_solution	Old solution vector that the solution is obtained from
older_solution	Older solution vector that the solution is obtained from

Definition at line 1989 of file XFEM.C.

Referenced by cutMeshWithEFA().

 {
   std::vector<dof_id_type> stored_solution_dofs = getNodeSolutionDofs(node_to_store_from, sys);
   std::vector<Real> stored_solution_scratch;
   // Size for current solution, as well as for old, and older solution only for transient case
   std::size_t stored_solution_size =
       (_fe_problem->isTransient() ? stored_solution_dofs.size() * 3 : stored_solution_dofs.size());
   stored_solution_scratch.reserve(stored_solution_size);
 
   // Store in the order defined in stored_solution_dofs first for the current, then for old and
   // older if applicable
   for (auto dof : stored_solution_dofs)
     stored_solution_scratch.push_back(current_solution(dof));
 
   if (_fe_problem->isTransient())
   {
     for (auto dof : stored_solution_dofs)
       stored_solution_scratch.push_back(old_solution(dof));
 
     for (auto dof : stored_solution_dofs)
       stored_solution_scratch.push_back(older_solution(dof));
   }
 
   if (stored_solution_scratch.size() > 0)
     stored_solution[node_to_store_to->unique_id()] = stored_solution_scratch;
 }

◆ update()

bool XFEM::update	(	Real	time,
		const std::vector< std::shared_ptr< NonlinearSystemBase >> &	nl,
		AuxiliarySystem &	aux
	)

overridevirtual

Implements XFEMInterface.

Definition at line 263 of file XFEM.C.

 {
   if (_moose_mesh->isDistributedMesh())
     mooseError("Use of XFEM with distributed mesh is not yet supported");
 
   bool mesh_changed = false;
 
   buildEFAMesh();
 
   _fe_problem->execute(EXEC_XFEM_MARK);
 
   storeCrackTipOriginAndDirection();
 
   if (markCuts(time))
     mesh_changed = cutMeshWithEFA(nl, aux);
 
   if (mesh_changed)
   {
     buildEFAMesh();
     storeCrackTipOriginAndDirection();
   }
 
   if (mesh_changed)
   {
     _mesh->allow_renumbering(false);
     _mesh->skip_partitioning(true);
     _mesh->prepare_for_use();
 
     if (_displaced_mesh)
     {
       _displaced_mesh->allow_renumbering(false);
       _displaced_mesh->skip_partitioning(true);
       _displaced_mesh->prepare_for_use();
     }
   }
 
   clearStateMarkedElems();
   clearGeomMarkedElems();
 
   return mesh_changed;
 }

◆ updateHeal()

bool XFEM::updateHeal ( )

overridevirtual

Implements XFEMInterface.

Definition at line 226 of file XFEM.C.

 {
   bool mesh_changed = false;
 
   mesh_changed = healMesh();
 
   if (mesh_changed)
     buildEFAMesh();
 
   if (mesh_changed)
   {
     _mesh->update_parallel_id_counts();
     MeshCommunication().make_elems_parallel_consistent(*_mesh);
     MeshCommunication().make_nodes_parallel_consistent(*_mesh);
     //    _mesh->find_neighbors();
     //    _mesh->contract();
     _mesh->allow_renumbering(false);
     _mesh->skip_partitioning(true);
     _mesh->prepare_for_use();
 
     if (_displaced_mesh)
     {
       _displaced_mesh->update_parallel_id_counts();
       MeshCommunication().make_elems_parallel_consistent(*_displaced_mesh);
       MeshCommunication().make_nodes_parallel_consistent(*_displaced_mesh);
       _displaced_mesh->allow_renumbering(false);
       _displaced_mesh->skip_partitioning(true);
       _displaced_mesh->prepare_for_use();
     }
   }
 
   _geom_marker_id_elems.clear();
 
   return mesh_changed;
 }

Member Data Documentation

◆ _cached_aux_solution

std::map<unique_id_type, std::vector<Real> > XFEM::_cached_aux_solution

private

Data structure to store the auxiliary solution for nodes/elements affected by XFEM For each node/element, this is stored as a vector that contains all components of all applicable variables in an order defined by getElementSolutionDofs() or getNodeSolutionDofs().

This vector first contains the current solution in that order, followed by the old and older solutions, also in that same order.

Definition at line 391 of file XFEM.h.

Referenced by cutMeshWithEFA(), and initSolution().

◆ _cached_solution

std::map<unique_id_type, std::vector<Real> > XFEM::_cached_solution

private

Data structure to store the nonlinear solution for nodes/elements affected by XFEM For each node/element, this is stored as a vector that contains all components of all applicable variables in an order defined by getElementSolutionDofs() or getNodeSolutionDofs().

This vector first contains the current solution in that order, followed by the old and older solutions, also in that same order.

Definition at line 382 of file XFEM.h.

Referenced by cutMeshWithEFA(), and initSolution().

◆ _crack_growth_increment

Real XFEM::_crack_growth_increment

private

Definition at line 332 of file XFEM.h.

Referenced by markCutEdgesByState(), and setCrackGrowthMethod().

◆ _crack_tip_elems

std::set<const Elem *> XFEM::_crack_tip_elems

private

Definition at line 337 of file XFEM.h.

Referenced by cutMeshWithEFA(), healMesh(), and isElemAtCrackTip().

◆ _crack_tip_elems_to_be_healed

std::set<const Elem *> XFEM::_crack_tip_elems_to_be_healed

private

Definition at line 338 of file XFEM.h.

Referenced by cutMeshWithEFA(), and healMesh().

◆ _cut_elem_map

std::map<unique_id_type, XFEMCutElem *> XFEM::_cut_elem_map

private

Definition at line 336 of file XFEM.h.

Referenced by buildEFAMesh(), cutMeshWithEFA(), getCutPlane(), getFragmentFaces(), getPhysicalVolumeFraction(), getXFEMIntersectionInfo(), healMesh(), isElemCut(), isPointInsidePhysicalDomain(), storeCrackTipOriginAndDirection(), and ~XFEM().

◆ _debug_output_level

unsigned int XFEM::_debug_output_level

private

Controls amount of debugging output information 0: None 1: Summary 2: Details on modifications to mesh 3: Full dump of element fragment algorithm mesh.

Definition at line 369 of file XFEM.h.

Referenced by cutMeshWithEFA(), healMesh(), and setDebugOutputLevel().

◆ _efa_mesh

ElementFragmentAlgorithm XFEM::_efa_mesh

private

Definition at line 362 of file XFEM.h.

Referenced by buildEFAMesh(), cutMeshWithEFA(), getEFAElem2D(), getEFAElem3D(), markCutEdgesByGeometry(), markCutEdgesByState(), markCutFacesByGeometry(), and storeCrackTipOriginAndDirection().

◆ _elem_crack_origin_direction_map

std::map<const Elem *, std::vector<Point> > XFEM::_elem_crack_origin_direction_map

private

Definition at line 342 of file XFEM.h.

Referenced by cutMeshWithEFA(), getCrackTipOrigin(), getCrackTipOriginMap(), healMesh(), markCutEdgesByState(), and storeCrackTipOriginAndDirection().

◆ _geom_cut_elems

std::unordered_map<const Elem *, Xfem::CutElemInfo> XFEM::_geom_cut_elems

private

All geometrically cut elements and their CutElemInfo during the current execution of XFEM_MARK.

This data structure is updated everytime a new cut element is created.

Definition at line 397 of file XFEM.h.

Referenced by cutMeshWithEFA(), healMesh(), storeMaterialPropertiesForElement(), and storeMaterialPropertiesForElementHelper().

◆ _geom_marked_elems_2d

std::map<const Elem *, std::vector<Xfem::GeomMarkedElemInfo2D> > XFEM::_geom_marked_elems_2d

private

Data structure for storing information about all 2D elements to be cut by geometry.

Definition at line 351 of file XFEM.h.

Referenced by addGeomMarkedElem2D(), clearGeomMarkedElems(), and markCutEdgesByGeometry().

◆ _geom_marked_elems_3d

std::map<const Elem *, std::vector<Xfem::GeomMarkedElemInfo3D> > XFEM::_geom_marked_elems_3d

private

Data structure for storing information about all 3D elements to be cut by geometry.

Definition at line 354 of file XFEM.h.

Referenced by addGeomMarkedElem3D(), clearGeomMarkedElems(), and markCutFacesByGeometry().

◆ _geom_marker_id_elems

std::map<unsigned int, std::set<unsigned int> > XFEM::_geom_marker_id_elems

private

Data structure for storing the elements cut by specific geometric cutters.

Definition at line 357 of file XFEM.h.

Referenced by addGeomMarkedElem2D(), addGeomMarkedElem3D(), cutMeshWithEFA(), getGeometricCutForElem(), and updateHeal().

◆ _geom_marker_id_map

std::map<const GeometricCutUserObject *, unsigned int> XFEM::_geom_marker_id_map

private

Data structure for storing the GeommetricCutUserObjects and their corresponding id.

Definition at line 360 of file XFEM.h.

Referenced by addGeometricCut(), and getGeometricCutID().

◆ _geometric_cuts

std::vector<const GeometricCutUserObject *> XFEM::_geometric_cuts

private

Definition at line 334 of file XFEM.h.

Referenced by addGeometricCut(), cutMeshWithEFA(), getGeometricCutForElem(), and healMesh().

◆ _has_secondary_cut

bool XFEM::_has_secondary_cut

private

Definition at line 327 of file XFEM.h.

Referenced by has_secondary_cut(), markCutEdgesByGeometry(), markCutEdgesByState(), and XFEM().

◆ _min_weight_multiplier

Real XFEM::_min_weight_multiplier

private

The minimum average multiplier applied by XFEM to the standard quadrature weights to integrate partial elements.

Definition at line 373 of file XFEM.h.

Referenced by getXFEMWeights(), and setMinWeightMultiplier().

◆ _old_geom_cut_elems

std::unordered_map<const Elem *, Xfem::CutElemInfo> XFEM::_old_geom_cut_elems

private

All geometrically cut elements and their CutElemInfo before the current execution of XFEM_MARK.

Definition at line 403 of file XFEM.h.

Referenced by cutMeshWithEFA().

◆ _sibling_displaced_elems

std::map<unsigned int, ElementPairLocator::ElementPairList> XFEM::_sibling_displaced_elems

private

Definition at line 340 of file XFEM.h.

Referenced by cutMeshWithEFA(), getXFEMDisplacedCutElemPairs(), and healMesh().

◆ _sibling_elems

std::map<unsigned int, ElementPairLocator::ElementPairList> XFEM::_sibling_elems

private

Definition at line 339 of file XFEM.h.

Referenced by cutMeshWithEFA(), getXFEMCutElemPairs(), and healMesh().

◆ _state_marked_elem_sides

std::map<const Elem *, unsigned int> XFEM::_state_marked_elem_sides

private

Definition at line 348 of file XFEM.h.

Referenced by addStateMarkedElem(), clearStateMarkedElems(), and markCutEdgesByState().

◆ _state_marked_elems

std::map<const Elem *, RealVectorValue> XFEM::_state_marked_elems

private

Definition at line 346 of file XFEM.h.

Referenced by addStateMarkedElem(), clearStateMarkedElems(), and markCutEdgesByState().

◆ _state_marked_frags

std::set<const Elem *> XFEM::_state_marked_frags

private

Definition at line 347 of file XFEM.h.

Referenced by addStateMarkedFrag(), clearStateMarkedElems(), and markCutEdgesByState().

◆ _use_crack_growth_increment

bool XFEM::_use_crack_growth_increment

private

Definition at line 331 of file XFEM.h.

Referenced by markCutEdgesByState(), and setCrackGrowthMethod().

◆ _XFEM_qrule

Xfem::XFEM_QRULE XFEM::_XFEM_qrule

private

Definition at line 329 of file XFEM.h.

Referenced by getXFEMQRule(), and setXFEMQRule().

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

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

Public Member Functions

Public Attributes

Protected Attributes

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ XFEM()

◆ ~XFEM()

Member Function Documentation

◆ addGeometricCut()

◆ addGeomMarkedElem2D()

◆ addGeomMarkedElem3D()

◆ addStateMarkedElem() [1/2]

◆ addStateMarkedElem() [2/2]

◆ addStateMarkedFrag()

◆ buildEFAMesh()

◆ clearGeomMarkedElems()

◆ clearStateMarkedElems()

◆ correctCrackExtensionDirection()

◆ cutMeshWithEFA()

◆ getCrackTipOrigin()

◆ getCrackTipOriginMap()

◆ getCutPlane()

◆ getCutSubdomainID()

◆ getEFAElem2D()

◆ getEFAElem3D()

◆ getEFANodeCoords()

◆ getElementSolutionDofs()

◆ getFragmentEdges()

◆ getFragmentFaces() [1/2]

◆ getFragmentFaces() [2/2]

◆ getGeometricCutForElem()

◆ getGeometricCutID()

◆ getNodeSolutionDofs()

◆ getPhysicalVolumeFraction()

◆ getXFEMCutElemPairs()

◆ getXFEMDisplacedCutElemPairs()

◆ getXFEMFaceWeights()

◆ getXFEMIntersectionInfo()

◆ getXFEMQRule()

◆ getXFEMqRuleOnLine()

◆ getXFEMqRuleOnSurface()

◆ getXFEMWeights()

◆ has_secondary_cut()

◆ healMesh()

◆ initCutIntersectionEdge()

◆ initSolution()

◆ isElemAtCrackTip()

◆ isElemCut() [1/2]

◆ isElemCut() [2/2]

◆ isPointInsidePhysicalDomain()

◆ loadMaterialPropertiesForElement()

◆ loadMaterialPropertiesForElementHelper()

◆ markCutEdgesByGeometry()

◆ markCutEdgesByState()

◆ markCutFacesByGeometry()

◆ markCutFacesByState()

◆ markCuts()

◆ pickFirstPhysicalNode()

◆ setCrackGrowthMethod()

◆ setDebugOutputLevel()

◆ setMinWeightMultiplier()

◆ setSolution()

◆ setSolutionForDOFs()

◆ setXFEMQRule()

◆ storeCrackTipOriginAndDirection()

◆ storeMaterialPropertiesForElement()

◆ storeMaterialPropertiesForElementHelper()

◆ storeSolutionForElement()

◆ storeSolutionForNode()

◆ update()

◆ updateHeal()

Member Data Documentation

◆ _cached_aux_solution

◆ _cached_solution

◆ _crack_growth_increment

◆ _crack_tip_elems

◆ _crack_tip_elems_to_be_healed

◆ _cut_elem_map