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

This is the XFEM class. More...

#include <XFEM.h>

Inheritance diagram for XFEM:
[legend]

Public Member Functions

 XFEM (const InputParameters &params)
 
 ~XFEM ()
 
void addGeometricCut (const 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 ()
 
virtual bool update (Real time, NonlinearSystemBase &nl, AuxiliarySystem &aux)
 Method to update the mesh due to modified cut planes. More...
 
virtual void initSolution (NonlinearSystemBase &nl, AuxiliarySystem &aux)
 Initialize the solution on newly created nodes. More...
 
void buildEFAMesh ()
 
bool markCuts (Real time)
 
bool markCutEdgesByGeometry (Real time)
 
bool markCutEdgesByState (Real time)
 
bool markCutFacesByGeometry (Real time)
 
bool markCutFacesByState ()
 
bool initCutIntersectionEdge (Point cut_origin, RealVectorValue cut_normal, Point &edge_p1, Point &edge_p2, Real &dist)
 
bool cutMeshWithEFA (NonlinearSystemBase &nl, AuxiliarySystem &aux)
 
Point getEFANodeCoords (EFANode *CEMnode, EFAElement *CEMElem, const Elem *elem, MeshBase *displaced_mesh=NULL) const
 
Real getPhysicalVolumeFraction (const Elem *elem) const
 Get the volume fraction of an element that is physical. 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)
 
std::vector< Real > & getXFEMCutData ()
 Set and get xfem cut data and type. More...
 
void setXFEMCutData (std::vector< Real > &cut_data)
 
std::string & getXFEMCutType ()
 
void setXFEMCutType (std::string &cut_type)
 
Xfem::XFEM_QRULEgetXFEMQRule ()
 
void setXFEMQRule (std::string &xfem_qrule)
 
void setCrackGrowthMethod (bool use_crack_growth_increment, Real crack_growth_increment)
 
virtual bool getXFEMWeights (MooseArray< Real > &weights, const Elem *elem, QBase *qrule, const MooseArray< Point > &q_points)
 
virtual const ElementPairLocator::ElementPairList * getXFEMCutElemPairs () const
 
virtual const ElementPairLocator::ElementPairList * getXFEMDisplacedCutElemPairs () const
 
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 ()
 

Private Member Functions

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

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
 
ElementPairLocator::ElementPairList _sibling_elems
 
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
 
ElementFragmentAlgorithm _efa_mesh
 
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...
 

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 60 of file XFEM.h.

Constructor & Destructor Documentation

XFEM::XFEM ( const InputParameters &  params)
explicit

Definition at line 29 of file XFEM.C.

29  : XFEMInterface(params), _efa_mesh(Moose::out)
30 {
31 #ifndef LIBMESH_ENABLE_UNIQUE_ID
32  mooseError("MOOSE requires unique ids to be enabled in libmesh (configure with "
33  "--enable-unique-id) to use XFEM!");
34 #endif
35  _has_secondary_cut = false;
36 }
ElementFragmentAlgorithm _efa_mesh
Definition: XFEM.h:199
bool _has_secondary_cut
Definition: XFEM.h:177
XFEM::~XFEM ( )

Definition at line 38 of file XFEM.C.

39 {
40  for (std::map<unique_id_type, XFEMCutElem *>::iterator cemit = _cut_elem_map.begin();
41  cemit != _cut_elem_map.end();
42  ++cemit)
43  delete cemit->second;
44 }
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:186

Member Function Documentation

void XFEM::addGeometricCut ( const GeometricCutUserObject geometric_cut)

Definition at line 47 of file XFEM.C.

48 {
49  _geometric_cuts.push_back(geometric_cut);
50 }
std::vector< const GeometricCutUserObject * > _geometric_cuts
Definition: XFEM.h:184
void XFEM::addStateMarkedElem ( unsigned int  elem_id,
RealVectorValue &  normal 
)

Definition at line 101 of file XFEM.C.

Referenced by addStateMarkedElem(), and addStateMarkedFrag().

102 {
103  Elem * elem = _mesh->elem(elem_id);
104  std::map<const Elem *, RealVectorValue>::iterator mit;
105  mit = _state_marked_elems.find(elem);
106  if (mit != _state_marked_elems.end())
107  mooseError(" ERROR: element ", elem->id(), " already marked for crack growth.");
108  _state_marked_elems[elem] = normal;
109 }
std::map< const Elem *, RealVectorValue > _state_marked_elems
Definition: XFEM.h:195
void XFEM::addStateMarkedElem ( unsigned int  elem_id,
RealVectorValue &  normal,
unsigned int  marked_side 
)

Definition at line 112 of file XFEM.C.

113 {
114  addStateMarkedElem(elem_id, normal);
115  Elem * elem = _mesh->elem(elem_id);
116  std::map<const Elem *, unsigned int>::iterator mit;
117  mit = _state_marked_elem_sides.find(elem);
118  if (mit != _state_marked_elem_sides.end())
119  {
120  mooseError(" ERROR: side of element ", elem->id(), " already marked for crack initiation.");
121  exit(1);
122  }
123  _state_marked_elem_sides[elem] = marked_side;
124 }
void addStateMarkedElem(unsigned int elem_id, RealVectorValue &normal)
Definition: XFEM.C:101
std::map< const Elem *, unsigned int > _state_marked_elem_sides
Definition: XFEM.h:197
void XFEM::addStateMarkedFrag ( unsigned int  elem_id,
RealVectorValue &  normal 
)

Definition at line 127 of file XFEM.C.

128 {
129  addStateMarkedElem(elem_id, normal);
130  Elem * elem = _mesh->elem(elem_id);
131  std::set<const Elem *>::iterator mit;
132  mit = _state_marked_frags.find(elem);
133  if (mit != _state_marked_frags.end())
134  {
135  mooseError(
136  " ERROR: element ", elem->id(), " already marked for fragment-secondary crack initiation.");
137  exit(1);
138  }
139  _state_marked_frags.insert(elem);
140 }
std::set< const Elem * > _state_marked_frags
Definition: XFEM.h:196
void addStateMarkedElem(unsigned int elem_id, RealVectorValue &normal)
Definition: XFEM.C:101
void XFEM::buildEFAMesh ( )

Definition at line 265 of file XFEM.C.

Referenced by update().

266 {
267  _efa_mesh.reset();
268 
269  MeshBase::element_iterator elem_it = _mesh->elements_begin();
270  const MeshBase::element_iterator elem_end = _mesh->elements_end();
271 
272  // Load all existing elements in to EFA mesh
273  for (elem_it = _mesh->elements_begin(); elem_it != elem_end; ++elem_it)
274  {
275  Elem * elem = *elem_it;
276  std::vector<unsigned int> quad;
277  for (unsigned int i = 0; i < elem->n_nodes(); ++i)
278  quad.push_back(elem->node(i));
279  if (_mesh->mesh_dimension() == 2)
280  _efa_mesh.add2DElement(quad, elem->id());
281  else if (_mesh->mesh_dimension() == 3)
282  _efa_mesh.add3DElement(quad, elem->id());
283  else
284  mooseError("XFEM only works for 2D and 3D");
285  }
286 
287  // Restore fragment information for elements that have been previously cut
288  for (elem_it = _mesh->elements_begin(); elem_it != elem_end; ++elem_it)
289  {
290  Elem * elem = *elem_it;
291  std::map<unique_id_type, XFEMCutElem *>::iterator cemit = _cut_elem_map.find(elem->unique_id());
292  if (cemit != _cut_elem_map.end())
293  {
294  XFEMCutElem * xfce = cemit->second;
295  EFAElement * CEMElem = _efa_mesh.getElemByID(elem->id());
296  _efa_mesh.restoreFragmentInfo(CEMElem, xfce->getEFAElement());
297  }
298  }
299 
300  // Must update edge neighbors before restore edge intersections. Otherwise, when we
301  // add edge intersections, we do not have neighbor information to use.
302  // Correction: no need to use neighbor info now
305 }
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:186
ElementFragmentAlgorithm _efa_mesh
Definition: XFEM.h:199
void restoreFragmentInfo(EFAElement *const elem, const EFAElement *const from_elem)
EFAElement * add2DElement(std::vector< unsigned int > quad, unsigned int id)
EFAElement * getElemByID(unsigned int id)
virtual const EFAElement * getEFAElement() const =0
EFAElement * add3DElement(std::vector< unsigned int > quad, unsigned int id)
void XFEM::clearStateMarkedElems ( )

Definition at line 143 of file XFEM.C.

Referenced by update().

144 {
145  _state_marked_elems.clear();
146  _state_marked_frags.clear();
147  _state_marked_elem_sides.clear();
148 }
std::set< const Elem * > _state_marked_frags
Definition: XFEM.h:196
std::map< const Elem *, unsigned int > _state_marked_elem_sides
Definition: XFEM.h:197
std::map< const Elem *, RealVectorValue > _state_marked_elems
Definition: XFEM.h:195
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 407 of file XFEM.C.

Referenced by markCutEdgesByState().

416 {
417  std::vector<Point> edge_ends(2, Point(0.0, 0.0, 0.0));
418  Point edge1(0.0, 0.0, 0.0);
419  Point edge2(0.0, 0.0, 0.0);
420  Point left_angle(0.0, 0.0, 0.0);
421  Point right_angle(0.0, 0.0, 0.0);
422  Point left_angle_normal(0.0, 0.0, 0.0);
423  Point right_angle_normal(0.0, 0.0, 0.0);
424  Point crack_direction_normal(0.0, 0.0, 0.0);
425  Point edge1_to_tip(0.0, 0.0, 0.0);
426  Point edge2_to_tip(0.0, 0.0, 0.0);
427  Point edge1_to_tip_normal(0.0, 0.0, 0.0);
428  Point edge2_to_tip_normal(0.0, 0.0, 0.0);
429 
430  Real cos_45 = std::cos(45.0 / 180.0 * 3.14159);
431  Real sin_45 = std::sin(45.0 / 180.0 * 3.14159);
432 
433  left_angle(0) = cos_45 * crack_tip_direction(0) - sin_45 * crack_tip_direction(1);
434  left_angle(1) = sin_45 * crack_tip_direction(0) + cos_45 * crack_tip_direction(1);
435 
436  right_angle(0) = cos_45 * crack_tip_direction(0) + sin_45 * crack_tip_direction(1);
437  right_angle(1) = -sin_45 * crack_tip_direction(0) + cos_45 * crack_tip_direction(1);
438 
439  left_angle_normal(0) = -left_angle(1);
440  left_angle_normal(1) = left_angle(0);
441 
442  right_angle_normal(0) = -right_angle(1);
443  right_angle_normal(1) = right_angle(0);
444 
445  crack_direction_normal(0) = -crack_tip_direction(1);
446  crack_direction_normal(1) = crack_tip_direction(0);
447 
448  Real angle_min = 0.0;
449  Real distance = 0.0;
450  unsigned int nsides = CEMElem->numEdges();
451 
452  for (unsigned int i = 0; i < nsides; ++i)
453  {
454  if (!orig_edge->isPartialOverlap(*CEMElem->getEdge(i)))
455  {
456  edge_ends[0] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
457  edge_ends[1] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
458 
459  edge1_to_tip = (edge_ends[0] * 0.95 + edge_ends[1] * 0.05) - crack_tip_origin;
460  edge2_to_tip = (edge_ends[0] * 0.05 + edge_ends[1] * 0.95) - crack_tip_origin;
461 
462  edge1_to_tip /= pow(edge1_to_tip.norm_sq(), 0.5);
463  edge2_to_tip /= pow(edge2_to_tip.norm_sq(), 0.5);
464 
465  edge1_to_tip_normal(0) = -edge1_to_tip(1);
466  edge1_to_tip_normal(1) = edge1_to_tip(0);
467 
468  edge2_to_tip_normal(0) = -edge2_to_tip(1);
469  edge2_to_tip_normal(1) = edge2_to_tip(0);
470 
471  Real angle_edge1_normal = edge1_to_tip_normal * normal;
472  Real angle_edge2_normal = edge2_to_tip_normal * normal;
473 
474  if (std::abs(angle_edge1_normal) > std::abs(angle_min) &&
475  (edge1_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
476  {
477  edge_id_keep = i;
478  distance_keep = 0.05;
479  normal_keep = edge1_to_tip_normal;
480  angle_min = angle_edge1_normal;
481  }
482  else if (std::abs(angle_edge2_normal) > std::abs(angle_min) &&
483  (edge2_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
484  {
485  edge_id_keep = i;
486  distance_keep = 0.95;
487  normal_keep = edge2_to_tip_normal;
488  angle_min = angle_edge2_normal;
489  }
490 
492  crack_tip_origin, left_angle_normal, edge_ends[0], edge_ends[1], distance) &&
493  (!CEMElem->isEdgePhantom(i)))
494  {
495  if (std::abs(left_angle_normal * normal) > std::abs(angle_min) &&
496  (edge1_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
497  {
498  edge_id_keep = i;
499  distance_keep = distance;
500  normal_keep = left_angle_normal;
501  angle_min = left_angle_normal * normal;
502  }
503  }
504  else if (initCutIntersectionEdge(
505  crack_tip_origin, right_angle_normal, edge_ends[0], edge_ends[1], distance) &&
506  (!CEMElem->isEdgePhantom(i)))
507  {
508  if (std::abs(right_angle_normal * normal) > std::abs(angle_min) &&
509  (edge2_to_tip * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
510  {
511  edge_id_keep = i;
512  distance_keep = distance;
513  normal_keep = right_angle_normal;
514  angle_min = right_angle_normal * normal;
515  }
516  }
517  else if (initCutIntersectionEdge(crack_tip_origin,
518  crack_direction_normal,
519  edge_ends[0],
520  edge_ends[1],
521  distance) &&
522  (!CEMElem->isEdgePhantom(i)))
523  {
524  if (std::abs(crack_direction_normal * normal) > std::abs(angle_min) &&
525  (crack_tip_direction * crack_tip_direction) > std::cos(45.0 / 180.0 * 3.14159))
526  {
527  edge_id_keep = i;
528  distance_keep = distance;
529  normal_keep = crack_direction_normal;
530  angle_min = crack_direction_normal * normal;
531  }
532  }
533  }
534  }
535 
536  // avoid small volume fraction cut
537  if ((distance_keep - 0.05) < 0.0)
538  {
539  distance_keep = 0.05;
540  }
541  else if ((distance_keep - 0.95) > 0.0)
542  {
543  distance_keep = 0.95;
544  }
545 }
EFAEdge * getEdge(unsigned int edge_id) const
bool initCutIntersectionEdge(Point cut_origin, RealVectorValue cut_normal, Point &edge_p1, Point &edge_p2, Real &dist)
Definition: XFEM.C:905
bool isEdgePhantom(unsigned int edge_id) const
Point getEFANodeCoords(EFANode *CEMnode, EFAElement *CEMElem, const Elem *elem, MeshBase *displaced_mesh=NULL) const
Definition: XFEM.C:1296
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
bool isPartialOverlap(const EFAEdge &other) const
Definition: EFAEdge.C:55
unsigned int numEdges() const
bool XFEM::cutMeshWithEFA ( NonlinearSystemBase &  nl,
AuxiliarySystem &  aux 
)

Definition at line 924 of file XFEM.C.

Referenced by update().

925 {
926  std::map<unsigned int, Node *> efa_id_to_new_node;
927  std::map<unsigned int, Node *> efa_id_to_new_node2;
928  std::map<unsigned int, Elem *> efa_id_to_new_elem;
929  _cached_solution.clear();
930  _cached_aux_solution.clear();
931 
933  // DEBUG
934  //_efa_mesh.printMesh();
935 
937  // DEBUG
938  //_efa_mesh.printMesh();
939 
940  const std::vector<EFANode *> new_nodes = _efa_mesh.getNewNodes();
941  const std::vector<EFAElement *> new_elements = _efa_mesh.getChildElements();
942  const std::vector<EFAElement *> delete_elements = _efa_mesh.getParentElements();
943 
944  bool mesh_changed = (new_nodes.size() + new_elements.size() + delete_elements.size() > 0);
945 
946  // Prepare to cache solution on DOFs modified by XFEM
947  if (mesh_changed)
948  {
949  nl.serializeSolution();
950  aux.serializeSolution();
951  }
952  NumericVector<Number> & current_solution = *nl.system().current_local_solution;
953  NumericVector<Number> & old_solution = nl.solutionOld();
954  NumericVector<Number> & older_solution = nl.solutionOlder();
955  NumericVector<Number> & current_aux_solution = *aux.system().current_local_solution;
956  NumericVector<Number> & old_aux_solution = aux.solutionOld();
957  NumericVector<Number> & older_aux_solution = aux.solutionOlder();
958 
959  std::map<Node *, Node *> new_nodes_to_parents;
960 
961  // Add new nodes
962  for (unsigned int i = 0; i < new_nodes.size(); ++i)
963  {
964  unsigned int new_node_id = new_nodes[i]->id();
965  unsigned int parent_id = new_nodes[i]->parent()->id();
966 
967  Node * parent_node = _mesh->node_ptr(parent_id);
968  Node * new_node = Node::build(*parent_node, _mesh->n_nodes()).release();
969  new_node->processor_id() = parent_node->processor_id();
970  _mesh->add_node(new_node);
971 
972  new_nodes_to_parents[new_node] = parent_node;
973 
974  new_node->set_n_systems(parent_node->n_systems());
975  efa_id_to_new_node.insert(std::make_pair(new_node_id, new_node));
976  _console << "XFEM added new node: " << new_node->id() << "\n";
977  if (_displaced_mesh)
978  {
979  const Node * parent_node2 = _displaced_mesh->node_ptr(parent_id);
980  Node * new_node2 = Node::build(*parent_node2, _displaced_mesh->n_nodes()).release();
981  new_node2->processor_id() = parent_node2->processor_id();
982  _displaced_mesh->add_node(new_node2);
983 
984  new_node2->set_n_systems(parent_node2->n_systems());
985  efa_id_to_new_node2.insert(std::make_pair(new_node_id, new_node2));
986  }
987  }
988 
989  // Add new elements
990  std::map<unsigned int, std::vector<const Elem *>> temporary_parent_children_map;
991 
992  for (unsigned int i = 0; i < new_elements.size(); ++i)
993  {
994  unsigned int parent_id = new_elements[i]->getParent()->id();
995  unsigned int efa_child_id = new_elements[i]->id();
996 
997  Elem * parent_elem = _mesh->elem(parent_id);
998  Elem * libmesh_elem = Elem::build(parent_elem->type()).release();
999 
1000  for (ElementPairLocator::ElementPairList::iterator it = _sibling_elems.begin();
1001  it != _sibling_elems.end();
1002  ++it)
1003  {
1004  if (parent_elem == it->first)
1005  it->first = libmesh_elem;
1006  else if (parent_elem == it->second)
1007  it->second = libmesh_elem;
1008  }
1009 
1010  // parent has at least two children
1011  if (new_elements[i]->getParent()->numChildren() > 1)
1012  temporary_parent_children_map[parent_id].push_back(libmesh_elem);
1013 
1014  Elem * parent_elem2 = NULL;
1015  Elem * libmesh_elem2 = NULL;
1016  if (_displaced_mesh)
1017  {
1018  parent_elem2 = _displaced_mesh->elem(parent_id);
1019  libmesh_elem2 = Elem::build(parent_elem2->type()).release();
1020 
1021  for (ElementPairLocator::ElementPairList::iterator it = _sibling_displaced_elems.begin();
1022  it != _sibling_displaced_elems.end();
1023  ++it)
1024  {
1025  if (parent_elem2 == it->first)
1026  it->first = libmesh_elem2;
1027  else if (parent_elem2 == it->second)
1028  it->second = libmesh_elem2;
1029  }
1030  }
1031 
1032  for (unsigned int j = 0; j < new_elements[i]->numNodes(); ++j)
1033  {
1034  unsigned int node_id = new_elements[i]->getNode(j)->id();
1035  Node * libmesh_node;
1036 
1037  std::map<unsigned int, Node *>::iterator nit = efa_id_to_new_node.find(node_id);
1038  if (nit != efa_id_to_new_node.end())
1039  libmesh_node = nit->second;
1040  else
1041  libmesh_node = _mesh->node_ptr(node_id);
1042 
1043  libmesh_elem->set_node(j) = libmesh_node;
1044 
1045  // Store solution for all nodes affected by XFEM (even existing nodes)
1046  if (parent_elem->is_semilocal(_mesh->processor_id()))
1047  {
1048  Node * solution_node = libmesh_node; // Node from which to store solution
1049  if (new_nodes_to_parents.find(libmesh_node) != new_nodes_to_parents.end())
1050  solution_node = new_nodes_to_parents[libmesh_node];
1051 
1052  if ((_moose_mesh->isSemiLocal(solution_node)) ||
1053  (solution_node->processor_id() == _mesh->processor_id()))
1054  {
1055  storeSolutionForNode(libmesh_node,
1056  solution_node,
1057  nl,
1059  current_solution,
1060  old_solution,
1061  older_solution);
1062  storeSolutionForNode(libmesh_node,
1063  solution_node,
1064  aux,
1066  current_aux_solution,
1067  old_aux_solution,
1068  older_aux_solution);
1069  }
1070  }
1071 
1072  Node * parent_node = parent_elem->get_node(j);
1073  std::vector<boundary_id_type> parent_node_boundary_ids =
1074  _mesh->boundary_info->boundary_ids(parent_node);
1075  _mesh->boundary_info->add_node(libmesh_node, parent_node_boundary_ids);
1076 
1077  if (_displaced_mesh)
1078  {
1079  std::map<unsigned int, Node *>::iterator nit2 = efa_id_to_new_node2.find(node_id);
1080  if (nit2 != efa_id_to_new_node2.end())
1081  libmesh_node = nit2->second;
1082  else
1083  libmesh_node = _displaced_mesh->node_ptr(node_id);
1084 
1085  libmesh_elem2->set_node(j) = libmesh_node;
1086 
1087  parent_node = parent_elem2->get_node(j);
1088  parent_node_boundary_ids.clear();
1089  parent_node_boundary_ids = _displaced_mesh->boundary_info->boundary_ids(parent_node);
1090  _displaced_mesh->boundary_info->add_node(libmesh_node, parent_node_boundary_ids);
1091  }
1092  }
1093 
1094  libmesh_elem->set_p_level(parent_elem->p_level());
1095  libmesh_elem->set_p_refinement_flag(parent_elem->p_refinement_flag());
1096  _mesh->add_elem(libmesh_elem);
1097  libmesh_elem->set_n_systems(parent_elem->n_systems());
1098  libmesh_elem->subdomain_id() = parent_elem->subdomain_id();
1099  libmesh_elem->processor_id() = parent_elem->processor_id();
1100 
1101  // TODO: The 0 here is the thread ID. Need to sort out how to do this correctly
1102  // TODO: Also need to copy neighbor material data
1103  if (parent_elem->processor_id() == _mesh->processor_id())
1104  {
1105  (*_material_data)[0]->copy(*libmesh_elem, *parent_elem, 0);
1106  for (unsigned int side = 0; side < parent_elem->n_sides(); ++side)
1107  {
1108  std::vector<boundary_id_type> parent_elem_boundary_ids =
1109  _mesh->boundary_info->boundary_ids(parent_elem, side);
1110  std::vector<boundary_id_type>::iterator it_bd = parent_elem_boundary_ids.begin();
1111  for (; it_bd != parent_elem_boundary_ids.end(); ++it_bd)
1112  {
1113  if (_fe_problem->needMaterialOnSide(*it_bd, 0))
1114  (*_bnd_material_data)[0]->copy(*libmesh_elem, *parent_elem, side);
1115  }
1116  }
1117 
1118  // Store solution for all elements affected by XFEM
1119  storeSolutionForElement(libmesh_elem,
1120  parent_elem,
1121  nl,
1123  current_solution,
1124  old_solution,
1125  older_solution);
1126  storeSolutionForElement(libmesh_elem,
1127  parent_elem,
1128  aux,
1130  current_aux_solution,
1131  old_aux_solution,
1132  older_aux_solution);
1133  }
1134 
1135  // The crack tip origin map is stored before cut, thus the elem should be updated with new
1136  // element.
1137  std::map<const Elem *, std::vector<Point>>::iterator mit =
1138  _elem_crack_origin_direction_map.find(parent_elem);
1139  if (mit != _elem_crack_origin_direction_map.end())
1140  {
1141  std::vector<Point> crack_data = _elem_crack_origin_direction_map[parent_elem];
1143  _elem_crack_origin_direction_map[libmesh_elem] = crack_data;
1144  }
1145 
1146  _console << "XFEM added new element: " << libmesh_elem->id() << "\n";
1147 
1148  XFEMCutElem * xfce = NULL;
1149  if (_mesh->mesh_dimension() == 2)
1150  {
1151  EFAElement2D * new_efa_elem2d = dynamic_cast<EFAElement2D *>(new_elements[i]);
1152  if (!new_efa_elem2d)
1153  mooseError("EFAelem is not of EFAelement2D type");
1154  xfce = new XFEMCutElem2D(libmesh_elem, new_efa_elem2d, (*_material_data)[0]->nQPoints());
1155  }
1156  else if (_mesh->mesh_dimension() == 3)
1157  {
1158  EFAElement3D * new_efa_elem3d = dynamic_cast<EFAElement3D *>(new_elements[i]);
1159  if (!new_efa_elem3d)
1160  mooseError("EFAelem is not of EFAelement3D type");
1161  xfce = new XFEMCutElem3D(libmesh_elem, new_efa_elem3d, (*_material_data)[0]->nQPoints());
1162  }
1163  _cut_elem_map.insert(std::pair<unique_id_type, XFEMCutElem *>(libmesh_elem->unique_id(), xfce));
1164  efa_id_to_new_elem.insert(std::make_pair(efa_child_id, libmesh_elem));
1165 
1166  if (_displaced_mesh)
1167  {
1168  libmesh_elem2->set_p_level(parent_elem2->p_level());
1169  libmesh_elem2->set_p_refinement_flag(parent_elem2->p_refinement_flag());
1170  _displaced_mesh->add_elem(libmesh_elem2);
1171  libmesh_elem2->set_n_systems(parent_elem2->n_systems());
1172  libmesh_elem2->subdomain_id() = parent_elem2->subdomain_id();
1173  libmesh_elem2->processor_id() = parent_elem2->processor_id();
1174  }
1175 
1176  unsigned int n_sides = parent_elem->n_sides();
1177  for (unsigned int side = 0; side < n_sides; ++side)
1178  {
1179  std::vector<boundary_id_type> parent_elem_boundary_ids =
1180  _mesh->boundary_info->boundary_ids(parent_elem, side);
1181  _mesh->boundary_info->add_side(libmesh_elem, side, parent_elem_boundary_ids);
1182  }
1183  if (_displaced_mesh)
1184  {
1185  n_sides = parent_elem2->n_sides();
1186  for (unsigned int side = 0; side < n_sides; ++side)
1187  {
1188  std::vector<boundary_id_type> parent_elem_boundary_ids =
1189  _displaced_mesh->boundary_info->boundary_ids(parent_elem2, side);
1190  _displaced_mesh->boundary_info->add_side(libmesh_elem2, side, parent_elem_boundary_ids);
1191  }
1192  }
1193 
1194  unsigned int n_edges = parent_elem->n_edges();
1195  for (unsigned int edge = 0; edge < n_edges; ++edge)
1196  {
1197  std::vector<boundary_id_type> parent_elem_boundary_ids =
1198  _mesh->boundary_info->edge_boundary_ids(parent_elem, edge);
1199  _mesh->boundary_info->add_edge(libmesh_elem, edge, parent_elem_boundary_ids);
1200  }
1201  if (_displaced_mesh)
1202  {
1203  n_edges = parent_elem2->n_edges();
1204  for (unsigned int edge = 0; edge < n_edges; ++edge)
1205  {
1206  std::vector<boundary_id_type> parent_elem_boundary_ids =
1207  _displaced_mesh->boundary_info->edge_boundary_ids(parent_elem2, edge);
1208  _displaced_mesh->boundary_info->add_edge(libmesh_elem2, edge, parent_elem_boundary_ids);
1209  }
1210  }
1211  }
1212 
1213  // delete elements
1214  for (std::size_t i = 0; i < delete_elements.size(); ++i)
1215  {
1216  Elem * elem_to_delete = _mesh->elem(delete_elements[i]->id());
1217 
1218  // delete the XFEMCutElem object for any elements that are to be deleted
1219  std::map<unique_id_type, XFEMCutElem *>::iterator cemit =
1220  _cut_elem_map.find(elem_to_delete->unique_id());
1221  if (cemit != _cut_elem_map.end())
1222  {
1223  delete cemit->second;
1224  _cut_elem_map.erase(cemit);
1225  }
1226 
1227  elem_to_delete->nullify_neighbors();
1228  _mesh->boundary_info->remove(elem_to_delete);
1229  unsigned int deleted_elem_id = elem_to_delete->id();
1230  _mesh->delete_elem(elem_to_delete);
1231  _console << "XFEM deleted element: " << deleted_elem_id << "\n";
1232 
1233  if (_displaced_mesh)
1234  {
1235  Elem * elem_to_delete2 = _displaced_mesh->elem(delete_elements[i]->id());
1236  elem_to_delete2->nullify_neighbors();
1237  _displaced_mesh->boundary_info->remove(elem_to_delete2);
1238  _displaced_mesh->delete_elem(elem_to_delete2);
1239  }
1240  }
1241 
1242  for (std::map<unsigned int, std::vector<const Elem *>>::iterator it =
1243  temporary_parent_children_map.begin();
1244  it != temporary_parent_children_map.end();
1245  ++it)
1246  {
1247  std::vector<const Elem *> & sibling_elem_vec = it->second;
1248  // TODO: for cut-node case, how to find the sibling elements?
1249  // if (sibling_elem_vec.size() != 2)
1250  // mooseError("Must have exactly 2 sibling elements");
1251  _sibling_elems.push_back(std::make_pair(sibling_elem_vec[0], sibling_elem_vec[1]));
1252  }
1253 
1254  // add sibling elems on displaced mesh
1255  if (_displaced_mesh)
1256  {
1257  for (ElementPairLocator::ElementPairList::iterator it = _sibling_elems.begin();
1258  it != _sibling_elems.end();
1259  ++it)
1260  {
1261  Elem * elem = _displaced_mesh->elem(it->first->id());
1262  Elem * elem_pair = _displaced_mesh->elem(it->second->id());
1263  _sibling_displaced_elems.push_back(std::make_pair(elem, elem_pair));
1264  }
1265  }
1266 
1267  // clear the temporary map
1268  temporary_parent_children_map.clear();
1269 
1270  // Store information about crack tip elements
1271  if (mesh_changed)
1272  {
1273  _crack_tip_elems.clear();
1274  const std::set<EFAElement *> CrackTipElements = _efa_mesh.getCrackTipElements();
1275  std::set<EFAElement *>::const_iterator sit;
1276  for (sit = CrackTipElements.begin(); sit != CrackTipElements.end(); ++sit)
1277  {
1278  unsigned int eid = (*sit)->id();
1279  Elem * crack_tip_elem;
1280  std::map<unsigned int, Elem *>::iterator eit = efa_id_to_new_elem.find(eid);
1281  if (eit != efa_id_to_new_elem.end())
1282  crack_tip_elem = eit->second;
1283  else
1284  crack_tip_elem = _mesh->elem(eid);
1285  _crack_tip_elems.insert(crack_tip_elem);
1286  }
1287  }
1288  _console << std::flush;
1289 
1290  // store virtual nodes
1291  // store cut edge info
1292  return mesh_changed;
1293 }
ElementPairLocator::ElementPairList _sibling_elems
Definition: XFEM.h:188
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:186
const std::vector< EFAElement * > & getChildElements()
std::set< const Elem * > _crack_tip_elems
Definition: XFEM.h:187
ElementFragmentAlgorithm _efa_mesh
Definition: XFEM.h:199
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/elem...
Definition: XFEM.h:217
void updateTopology(bool mergeUncutVirtualEdges=true)
const std::vector< EFAElement * > & getParentElements()
const std::vector< EFANode * > & getNewNodes()
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.
Definition: XFEM.C:1637
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.
Definition: XFEM.C:1603
std::map< const Elem *, std::vector< Point > > _elem_crack_origin_direction_map
Definition: XFEM.h:191
ElementPairLocator::ElementPairList _sibling_displaced_elems
Definition: XFEM.h:189
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/elem...
Definition: XFEM.h:208
const std::set< EFAElement * > & getCrackTipElements()
void XFEM::getCrackTipOrigin ( std::map< unsigned int, const Elem * > &  elem_id_crack_tip,
std::vector< Point > &  crack_front_points 
)

Definition at line 53 of file XFEM.C.

55 {
56  elem_id_crack_tip.clear();
57  crack_front_points.clear();
58  crack_front_points.resize(_elem_crack_origin_direction_map.size());
59 
60  unsigned int crack_tip_index = 0;
61  // This map is used to sort the order in _elem_crack_origin_direction_map such that every process
62  // has same order
63  std::map<unsigned int, const Elem *> elem_id_map;
64 
65  int m = -1;
66  for (std::map<const Elem *, std::vector<Point>>::iterator mit1 =
69  ++mit1)
70  {
71  unsigned int elem_id = mit1->first->id();
72  if (elem_id > 999999)
73  {
74  elem_id_map[m] = mit1->first;
75  m--;
76  }
77  else
78  {
79  elem_id_map[elem_id] = mit1->first;
80  }
81  }
82 
83  for (std::map<unsigned int, const Elem *>::iterator mit1 = elem_id_map.begin();
84  mit1 != elem_id_map.end();
85  mit1++)
86  {
87  const Elem * elem = mit1->second;
88  std::map<const Elem *, std::vector<Point>>::iterator mit2 =
90  if (mit2 != _elem_crack_origin_direction_map.end())
91  {
92  elem_id_crack_tip[crack_tip_index] = mit2->first;
93  crack_front_points[crack_tip_index] =
94  (mit2->second)[0]; // [0] stores origin coordinates and [1] stores direction
95  crack_tip_index++;
96  }
97  }
98 }
std::map< const Elem *, std::vector< Point > > _elem_crack_origin_direction_map
Definition: XFEM.h:191
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 1348 of file XFEM.C.

1351 {
1352  Real comp = 0.0;
1353  Point planedata(0.0, 0.0, 0.0);
1354  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1355  it = _cut_elem_map.find(elem->unique_id());
1356  if (it != _cut_elem_map.end())
1357  {
1358  const XFEMCutElem * xfce = it->second;
1359  const EFAElement * EFAelem = xfce->getEFAElement();
1360  if (EFAelem->isPartial()) // exclude the full crack tip elements
1361  {
1362  if ((unsigned int)quantity < 3)
1363  {
1364  unsigned int index = (unsigned int)quantity;
1365  planedata = xfce->getCutPlaneOrigin(plane_id, _displaced_mesh);
1366  comp = planedata(index);
1367  }
1368  else if ((unsigned int)quantity < 6)
1369  {
1370  unsigned int index = (unsigned int)quantity - 3;
1371  planedata = xfce->getCutPlaneNormal(plane_id, _displaced_mesh);
1372  comp = planedata(index);
1373  }
1374  else
1375  mooseError("In get_cut_plane index out of range");
1376  }
1377  }
1378  return comp;
1379 }
virtual bool isPartial() const =0
virtual Point getCutPlaneOrigin(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const =0
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:186
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const =0
virtual const EFAElement * getEFAElement() const =0
Point XFEM::getEFANodeCoords ( EFANode CEMnode,
EFAElement CEMElem,
const Elem *  elem,
MeshBase *  displaced_mesh = NULL 
) const

Definition at line 1296 of file XFEM.C.

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

1300 {
1301  Point node_coor(0.0, 0.0, 0.0);
1302  std::vector<EFANode *> master_nodes;
1303  std::vector<Point> master_points;
1304  std::vector<double> master_weights;
1305 
1306  CEMElem->getMasterInfo(CEMnode, master_nodes, master_weights);
1307  for (std::size_t i = 0; i < master_nodes.size(); ++i)
1308  {
1309  if (master_nodes[i]->category() == EFANode::N_CATEGORY_PERMANENT)
1310  {
1311  unsigned int local_node_id = CEMElem->getLocalNodeIndex(master_nodes[i]);
1312  Node * node = elem->get_node(local_node_id);
1313  if (displaced_mesh)
1314  node = displaced_mesh->node_ptr(node->id());
1315  Point node_p((*node)(0), (*node)(1), (*node)(2));
1316  master_points.push_back(node_p);
1317  }
1318  else
1319  mooseError("master nodes must be permanent");
1320  }
1321  for (std::size_t i = 0; i < master_nodes.size(); ++i)
1322  node_coor += master_weights[i] * master_points[i];
1323 
1324  return node_coor;
1325 }
unsigned int getLocalNodeIndex(EFANode *node) const
Definition: EFAElement.C:109
virtual void getMasterInfo(EFANode *node, std::vector< EFANode * > &master_nodes, std::vector< double > &master_weights) const =0
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
elemElement for which dof indices are found
sysSystem for which the dof indices are found

Definition at line 1736 of file XFEM.C.

Referenced by setSolution(), and storeSolutionForElement().

1737 {
1738  SubdomainID sid = elem->subdomain_id();
1739  const std::vector<MooseVariable *> & vars = sys.getVariables(0);
1740  std::vector<dof_id_type> solution_dofs;
1741  solution_dofs.reserve(vars.size()); // just an approximation
1742  for (auto var : vars)
1743  {
1744  if (!var->isNodal())
1745  {
1746  const std::set<SubdomainID> & var_subdomains = sys.getSubdomainsForVar(var->number());
1747  if (var_subdomains.empty() || var_subdomains.find(sid) != var_subdomains.end())
1748  {
1749  unsigned int n_comp = elem->n_comp(sys.number(), var->number());
1750  for (unsigned int icomp = 0; icomp < n_comp; ++icomp)
1751  {
1752  dof_id_type elem_dof = elem->dof_number(sys.number(), var->number(), icomp);
1753  solution_dofs.push_back(elem_dof);
1754  }
1755  }
1756  }
1757  }
1758  return solution_dofs;
1759 }
void XFEM::getFragmentEdges ( const Elem *  elem,
EFAElement2D CEMElem,
std::vector< std::vector< Point >> &  frag_edges 
) const
private

Definition at line 1429 of file XFEM.C.

Referenced by markCutEdgesByGeometry().

1432 {
1433  // N.B. CEMElem here has global EFAnode
1434  frag_edges.clear();
1435  if (CEMElem->numFragments() > 0)
1436  {
1437  if (CEMElem->numFragments() > 1)
1438  mooseError("element ", elem->id(), " has more than one fragments at this point");
1439  for (unsigned int i = 0; i < CEMElem->getFragment(0)->numEdges(); ++i)
1440  {
1441  std::vector<Point> p_line(2, Point(0.0, 0.0, 0.0));
1442  p_line[0] = getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
1443  p_line[1] = getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
1444  frag_edges.push_back(p_line);
1445  }
1446  }
1447 }
unsigned int numEdges() const
EFAFragment2D * getFragment(unsigned int frag_id) const
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
Point getEFANodeCoords(EFANode *CEMnode, EFAElement *CEMElem, const Elem *elem, MeshBase *displaced_mesh=NULL) const
Definition: XFEM.C:1296
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
virtual unsigned int numFragments() const
Definition: EFAElement2D.C:200
void XFEM::getFragmentFaces ( const Elem *  elem,
std::vector< std::vector< Point >> &  frag_faces,
bool  displaced_mesh = false 
) const

Definition at line 1412 of file XFEM.C.

Referenced by markCutFacesByGeometry().

1415 {
1416  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1417  it = _cut_elem_map.find(elem->unique_id());
1418  if (it != _cut_elem_map.end())
1419  {
1420  const XFEMCutElem * xfce = it->second;
1421  if (displaced_mesh)
1422  xfce->getFragmentFaces(frag_faces, _displaced_mesh);
1423  else
1424  xfce->getFragmentFaces(frag_faces);
1425  }
1426 }
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:186
virtual void getFragmentFaces(std::vector< std::vector< Point >> &frag_faces, MeshBase *displaced_mesh=NULL) const =0
void XFEM::getFragmentFaces ( const Elem *  elem,
EFAElement3D CEMElem,
std::vector< std::vector< Point >> &  frag_faces 
) const
private

Definition at line 1450 of file XFEM.C.

1453 {
1454  // N.B. CEMElem here has global EFAnode
1455  frag_faces.clear();
1456  if (CEMElem->numFragments() > 0)
1457  {
1458  if (CEMElem->numFragments() > 1)
1459  mooseError("element ", elem->id(), " has more than one fragments at this point");
1460  for (unsigned int i = 0; i < CEMElem->getFragment(0)->numFaces(); ++i)
1461  {
1462  unsigned int num_face_nodes = CEMElem->getFragmentFace(0, i)->numNodes();
1463  std::vector<Point> p_line(num_face_nodes, Point(0.0, 0.0, 0.0));
1464  for (unsigned int j = 0; j < num_face_nodes; ++j)
1465  p_line[j] = getEFANodeCoords(CEMElem->getFragmentFace(0, i)->getNode(j), CEMElem, elem);
1466  frag_faces.push_back(p_line);
1467  }
1468  }
1469 }
unsigned int numNodes() const
Definition: EFAFace.C:85
EFAFragment3D * getFragment(unsigned int frag_id) const
EFAFace * getFragmentFace(unsigned int frag_id, unsigned int face_id) const
Point getEFANodeCoords(EFANode *CEMnode, EFAElement *CEMElem, const Elem *elem, MeshBase *displaced_mesh=NULL) const
Definition: XFEM.C:1296
virtual unsigned int numFragments() const
Definition: EFAElement3D.C:263
unsigned int numFaces() const
EFANode * getNode(unsigned int node_id) const
Definition: EFAFace.C:97
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
nodeNode for which dof indices are found
sysSystem for which the dof indices are found

Definition at line 1762 of file XFEM.C.

Referenced by setSolution(), and storeSolutionForNode().

1763 {
1764  const std::set<SubdomainID> & sids = _moose_mesh->getNodeBlockIds(*node);
1765  const std::vector<MooseVariable *> & vars = sys.getVariables(0);
1766  std::vector<dof_id_type> solution_dofs;
1767  solution_dofs.reserve(vars.size()); // just an approximation
1768  for (auto var : vars)
1769  {
1770  if (var->isNodal())
1771  {
1772  const std::set<SubdomainID> & var_subdomains = sys.getSubdomainsForVar(var->number());
1773  std::set<SubdomainID> intersect;
1774  set_intersection(var_subdomains.begin(),
1775  var_subdomains.end(),
1776  sids.begin(),
1777  sids.end(),
1778  std::inserter(intersect, intersect.begin()));
1779  if (var_subdomains.empty() || !intersect.empty())
1780  {
1781  unsigned int n_comp = node->n_comp(sys.number(), var->number());
1782  for (unsigned int icomp = 0; icomp < n_comp; ++icomp)
1783  {
1784  dof_id_type node_dof = node->dof_number(sys.number(), var->number(), icomp);
1785  solution_dofs.push_back(node_dof);
1786  }
1787  }
1788  }
1789  }
1790  return solution_dofs;
1791 }
Real XFEM::getPhysicalVolumeFraction ( const Elem *  elem) const

Get the volume fraction of an element that is physical.

Definition at line 1328 of file XFEM.C.

Referenced by markCutEdgesByState().

1329 {
1330  Real phys_volfrac = 1.0;
1331  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1332  it = _cut_elem_map.find(elem->unique_id());
1333  if (it != _cut_elem_map.end())
1334  {
1335  XFEMCutElem * xfce = it->second;
1336  const EFAElement * EFAelem = xfce->getEFAElement();
1337  if (EFAelem->isPartial())
1338  { // exclude the full crack tip elements
1340  phys_volfrac = xfce->getPhysicalVolumeFraction();
1341  }
1342  }
1343 
1344  return phys_volfrac;
1345 }
virtual bool isPartial() const =0
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:186
Real getPhysicalVolumeFraction() const
Definition: XFEMCutElem.C:31
virtual void computePhysicalVolumeFraction()=0
virtual const EFAElement * getEFAElement() const =0
std::vector<Real>& XFEM::getXFEMCutData ( )

Set and get xfem cut data and type.

virtual const ElementPairLocator::ElementPairList* XFEM::getXFEMCutElemPairs ( ) const
inlinevirtual

Definition at line 147 of file XFEM.h.

148  {
149  return &_sibling_elems;
150  }
ElementPairLocator::ElementPairList _sibling_elems
Definition: XFEM.h:188
std::string& XFEM::getXFEMCutType ( )
virtual const ElementPairLocator::ElementPairList* XFEM::getXFEMDisplacedCutElemPairs ( ) const
inlinevirtual

Definition at line 151 of file XFEM.h.

152  {
153  return &_sibling_displaced_elems;
154  }
ElementPairLocator::ElementPairList _sibling_displaced_elems
Definition: XFEM.h:189
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 1513 of file XFEM.C.

1518 {
1519  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1520  it = _cut_elem_map.find(elem->unique_id());
1521  if (it != _cut_elem_map.end())
1522  {
1523  const XFEMCutElem * xfce = it->second;
1524  if (displaced_mesh)
1525  xfce->getIntersectionInfo(plane_id, normal, intersectionPoints, _displaced_mesh);
1526  else
1527  xfce->getIntersectionInfo(plane_id, normal, intersectionPoints);
1528  }
1529 }
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:186
virtual void getIntersectionInfo(unsigned int plane_id, Point &normal, std::vector< Point > &intersectionPoints, MeshBase *displaced_mesh=NULL) const =0
Xfem::XFEM_QRULE & XFEM::getXFEMQRule ( )

Definition at line 1472 of file XFEM.C.

Referenced by getXFEMWeights().

1473 {
1474  return _XFEM_qrule;
1475 }
Xfem::XFEM_QRULE _XFEM_qrule
Definition: XFEM.h:179
void XFEM::getXFEMqRuleOnLine ( std::vector< Point > &  intersection_points,
std::vector< Point > &  quad_pts,
std::vector< Real > &  quad_wts 
) const
virtual

Definition at line 1532 of file XFEM.C.

1535 {
1536  Point p1 = intersection_points[0];
1537  Point p2 = intersection_points[1];
1538 
1539  // number of quadrature points
1540  std::size_t num_qpoints = 2;
1541 
1542  // quadrature coordinates
1543  Real xi0 = -std::sqrt(1.0 / 3.0);
1544  Real xi1 = std::sqrt(1.0 / 3.0);
1545 
1546  quad_wts.resize(num_qpoints);
1547  quad_pts.resize(num_qpoints);
1548 
1549  Real integ_jacobian = pow((p1 - p2).size_sq(), 0.5) * 0.5;
1550 
1551  quad_wts[0] = 1.0 * integ_jacobian;
1552  quad_wts[1] = 1.0 * integ_jacobian;
1553 
1554  quad_pts[0] = (1.0 - xi0) / 2.0 * p1 + (1.0 + xi0) / 2.0 * p2;
1555  quad_pts[1] = (1.0 - xi1) / 2.0 * p1 + (1.0 + xi1) / 2.0 * p2;
1556 }
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
void XFEM::getXFEMqRuleOnSurface ( std::vector< Point > &  intersection_points,
std::vector< Point > &  quad_pts,
std::vector< Real > &  quad_wts 
) const
virtual

Definition at line 1559 of file XFEM.C.

1562 {
1563  std::size_t nnd_pe = intersection_points.size();
1564  Point xcrd(0.0, 0.0, 0.0);
1565  for (std::size_t i = 0; i < nnd_pe; ++i)
1566  xcrd += intersection_points[i];
1567  xcrd /= nnd_pe;
1568 
1569  quad_pts.resize(nnd_pe);
1570  quad_wts.resize(nnd_pe);
1571 
1572  Real jac = 0.0;
1573 
1574  for (std::size_t j = 0; j < nnd_pe; ++j) // loop all sub-tris
1575  {
1576  std::vector<std::vector<Real>> shape(3, std::vector<Real>(3, 0.0));
1577  std::vector<Point> subtrig_points(3, Point(0.0, 0.0, 0.0)); // sub-trig nodal coords
1578 
1579  int jplus1 = j < nnd_pe - 1 ? j + 1 : 0;
1580  subtrig_points[0] = xcrd;
1581  subtrig_points[1] = intersection_points[j];
1582  subtrig_points[2] = intersection_points[jplus1];
1583 
1584  std::vector<std::vector<Real>> sg2;
1585  Xfem::stdQuadr2D(3, 1, sg2); // get sg2
1586  for (std::size_t l = 0; l < sg2.size(); ++l) // loop all int pts on a sub-trig
1587  {
1588  Xfem::shapeFunc2D(3, sg2[l], subtrig_points, shape, jac, true); // Get shape
1589  std::vector<Real> tsg_line(3, 0.0);
1590  for (std::size_t k = 0; k < 3; ++k) // loop sub-trig nodes
1591  {
1592  tsg_line[0] += shape[k][2] * subtrig_points[k](0);
1593  tsg_line[1] += shape[k][2] * subtrig_points[k](1);
1594  tsg_line[2] += shape[k][2] * subtrig_points[k](2);
1595  }
1596  quad_pts[j + l] = Point(tsg_line[0], tsg_line[1], tsg_line[2]);
1597  quad_wts[j + l] = sg2[l][3] * jac;
1598  }
1599  }
1600 }
void shapeFunc2D(unsigned int nen, std::vector< Real > &ss, std::vector< Point > &xl, std::vector< std::vector< Real >> &shp, Real &xsj, bool natl_flg)
Definition: XFEMFuncs.C:274
void stdQuadr2D(unsigned int nen, unsigned int iord, std::vector< std::vector< Real >> &sg2)
Definition: XFEMFuncs.C:101
bool XFEM::getXFEMWeights ( MooseArray< Real > &  weights,
const Elem *  elem,
QBase *  qrule,
const MooseArray< Point > &  q_points 
)
virtual

Definition at line 1496 of file XFEM.C.

1500 {
1501  bool have_weights = false;
1502  XFEMCutElem * xfce = NULL;
1503  if (isElemCut(elem, xfce))
1504  {
1505  mooseAssert(xfce != NULL, "Must have valid XFEMCutElem object here");
1506  xfce->getWeightMultipliers(weights, qrule, getXFEMQRule(), q_points);
1507  have_weights = true;
1508  }
1509  return have_weights;
1510 }
bool isElemCut(const Elem *elem, XFEMCutElem *&xfce) const
Definition: XFEM.C:1388
Xfem::XFEM_QRULE & getXFEMQRule()
Definition: XFEM.C:1472
void getWeightMultipliers(MooseArray< Real > &weights, QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
Definition: XFEMCutElem.C:37
bool XFEM::has_secondary_cut ( )
inline

Definition at line 166 of file XFEM.h.

166 { return _has_secondary_cut; }
bool _has_secondary_cut
Definition: XFEM.h:177
bool XFEM::initCutIntersectionEdge ( Point  cut_origin,
RealVectorValue  cut_normal,
Point &  edge_p1,
Point &  edge_p2,
Real &  dist 
)

Definition at line 905 of file XFEM.C.

Referenced by correctCrackExtensionDirection(), and markCutEdgesByState().

907 {
908  dist = 0.0;
909  bool does_intersect = false;
910  Point origin2p1 = edge_p1 - cut_origin;
911  Real plane2p1 = cut_normal(0) * origin2p1(0) + cut_normal(1) * origin2p1(1);
912  Point origin2p2 = edge_p2 - cut_origin;
913  Real plane2p2 = cut_normal(0) * origin2p2(0) + cut_normal(1) * origin2p2(1);
914 
915  if (plane2p1 * plane2p2 < 0.0)
916  {
917  dist = -plane2p1 / (plane2p2 - plane2p1);
918  does_intersect = true;
919  }
920  return does_intersect;
921 }
void XFEM::initSolution ( NonlinearSystemBase &  nl,
AuxiliarySystem &  aux 
)
virtual

Initialize the solution on newly created nodes.

Definition at line 238 of file XFEM.C.

239 {
240  nl.serializeSolution();
241  aux.serializeSolution();
242  NumericVector<Number> & current_solution = *nl.system().current_local_solution;
243  NumericVector<Number> & old_solution = nl.solutionOld();
244  NumericVector<Number> & older_solution = nl.solutionOlder();
245  NumericVector<Number> & current_aux_solution = *aux.system().current_local_solution;
246  NumericVector<Number> & old_aux_solution = aux.solutionOld();
247  NumericVector<Number> & older_aux_solution = aux.solutionOlder();
248 
249  setSolution(nl, _cached_solution, current_solution, old_solution, older_solution);
250  setSolution(
251  aux, _cached_aux_solution, current_aux_solution, old_aux_solution, older_aux_solution);
252 
253  current_solution.close();
254  old_solution.close();
255  older_solution.close();
256  current_aux_solution.close();
257  old_aux_solution.close();
258  older_aux_solution.close();
259 
260  _cached_solution.clear();
261  _cached_aux_solution.clear();
262 }
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.
Definition: XFEM.C:1671
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/elem...
Definition: XFEM.h:217
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/elem...
Definition: XFEM.h:208
bool XFEM::isElemAtCrackTip ( const Elem *  elem) const

Definition at line 1382 of file XFEM.C.

Referenced by markCutEdgesByGeometry(), and markCutEdgesByState().

1383 {
1384  return (_crack_tip_elems.find(elem) != _crack_tip_elems.end());
1385 }
std::set< const Elem * > _crack_tip_elems
Definition: XFEM.h:187
bool XFEM::isElemCut ( const Elem *  elem,
XFEMCutElem *&  xfce 
) const

Definition at line 1388 of file XFEM.C.

Referenced by getXFEMWeights(), and isElemCut().

1389 {
1390  xfce = NULL;
1391  bool is_cut = false;
1392  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
1393  it = _cut_elem_map.find(elem->unique_id());
1394  if (it != _cut_elem_map.end())
1395  {
1396  xfce = it->second;
1397  const EFAElement * EFAelem = xfce->getEFAElement();
1398  if (EFAelem->isPartial()) // exclude the full crack tip elements
1399  is_cut = true;
1400  }
1401  return is_cut;
1402 }
virtual bool isPartial() const =0
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:186
virtual const EFAElement * getEFAElement() const =0
bool XFEM::isElemCut ( const Elem *  elem) const

Definition at line 1405 of file XFEM.C.

1406 {
1407  XFEMCutElem * xfce = NULL;
1408  return isElemCut(elem, xfce);
1409 }
bool isElemCut(const Elem *elem, XFEMCutElem *&xfce) const
Definition: XFEM.C:1388
bool XFEM::markCutEdgesByGeometry ( Real  time)

Definition at line 325 of file XFEM.C.

Referenced by markCuts().

326 {
327  bool marked_edges = false;
328  bool marked_nodes = false;
329 
330  std::vector<const GeometricCutUserObject *> active_geometric_cuts;
331  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
332  if (_geometric_cuts[i]->active(time))
333  {
334  active_geometric_cuts.push_back(_geometric_cuts[i]);
335  }
336 
337  if (active_geometric_cuts.size() > 0)
338  {
339  for (MeshBase::element_iterator elem_it = _mesh->elements_begin();
340  elem_it != _mesh->elements_end();
341  ++elem_it)
342  {
343  const Elem * elem = *elem_it;
344  std::vector<CutEdge> elem_cut_edges;
345  std::vector<CutNode> elem_cut_nodes;
346  std::vector<CutEdge> frag_cut_edges;
347  std::vector<std::vector<Point>> frag_edges;
348  EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
349  EFAElement2D * CEMElem = dynamic_cast<EFAElement2D *>(EFAelem);
350 
351  if (!CEMElem)
352  mooseError("EFAelem is not of EFAelement2D type");
353 
354  // continue if elem has been already cut twice - IMPORTANT
355  if (CEMElem->isFinalCut())
356  continue;
357 
358  // get fragment edges
359  getFragmentEdges(elem, CEMElem, frag_edges);
360 
361  // mark cut edges for the element and its fragment
362  for (unsigned int i = 0; i < active_geometric_cuts.size(); ++i)
363  {
364  active_geometric_cuts[i]->cutElementByGeometry(elem, elem_cut_edges, elem_cut_nodes, time);
365  if (CEMElem->numFragments() > 0)
366  active_geometric_cuts[i]->cutFragmentByGeometry(frag_edges, frag_cut_edges, time);
367  }
368 
369  for (unsigned int i = 0; i < elem_cut_edges.size(); ++i) // mark element edges
370  {
371  if (!CEMElem->isEdgePhantom(elem_cut_edges[i].host_side_id)) // must not be phantom edge
372  {
374  elem->id(), elem_cut_edges[i].host_side_id, elem_cut_edges[i].distance);
375  marked_edges = true;
376  }
377  }
378 
379  for (unsigned int i = 0; i < elem_cut_nodes.size(); ++i) // mark element edges
380  {
381  _efa_mesh.addElemNodeIntersection(elem->id(), elem_cut_nodes[i].host_id);
382  marked_nodes = true;
383  }
384 
385  for (unsigned int i = 0; i < frag_cut_edges.size();
386  ++i) // MUST DO THIS AFTER MARKING ELEMENT EDGES
387  {
388  if (!CEMElem->getFragment(0)->isSecondaryInteriorEdge(frag_cut_edges[i].host_side_id))
389  {
391  elem->id(), frag_cut_edges[i].host_side_id, frag_cut_edges[i].distance))
392  {
393  marked_edges = true;
394 
395  if (!isElemAtCrackTip(elem))
396  _has_secondary_cut = true;
397  }
398  }
399  }
400  }
401  }
402 
403  return marked_edges || marked_nodes;
404 }
EFAFragment2D * getFragment(unsigned int frag_id) const
void addElemNodeIntersection(unsigned int elemid, unsigned int nodeid)
void addElemEdgeIntersection(unsigned int elemid, unsigned int edgeid, double position)
std::vector< const GeometricCutUserObject * > _geometric_cuts
Definition: XFEM.h:184
ElementFragmentAlgorithm _efa_mesh
Definition: XFEM.h:199
virtual bool isFinalCut() const
Definition: EFAElement2D.C:789
bool _has_secondary_cut
Definition: XFEM.h:177
bool isSecondaryInteriorEdge(unsigned int edge_id) const
bool addFragEdgeIntersection(unsigned int elemid, unsigned int frag_edge_id, double position)
bool isEdgePhantom(unsigned int edge_id) const
void getFragmentEdges(const Elem *elem, EFAElement2D *CEMElem, std::vector< std::vector< Point >> &frag_edges) const
Definition: XFEM.C:1429
virtual unsigned int numFragments() const
Definition: EFAElement2D.C:200
EFAElement * getElemByID(unsigned int id)
bool isElemAtCrackTip(const Elem *elem) const
Definition: XFEM.C:1382
bool XFEM::markCutEdgesByState ( Real  time)

Definition at line 548 of file XFEM.C.

Referenced by markCuts().

549 {
550  bool marked_edges = false;
551  for (std::map<const Elem *, RealVectorValue>::iterator pmeit = _state_marked_elems.begin();
552  pmeit != _state_marked_elems.end();
553  ++pmeit)
554  {
555  const Elem * elem = pmeit->first;
556  RealVectorValue normal = pmeit->second;
557  EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
558  EFAElement2D * CEMElem = dynamic_cast<EFAElement2D *>(EFAelem);
559 
560  Real volfrac_elem = getPhysicalVolumeFraction(elem);
561  if (volfrac_elem < 0.25)
562  continue;
563 
564  if (!CEMElem)
565  mooseError("EFAelem is not of EFAelement2D type");
566 
567  // continue if elem is already cut twice - IMPORTANT
568  if (CEMElem->isFinalCut())
569  continue;
570 
571  // find the first cut edge
572  unsigned int nsides = CEMElem->numEdges();
573  unsigned int orig_cut_side_id = 999999;
574  Real orig_cut_distance = -1.0;
575  EFANode * orig_node = NULL;
576  EFAEdge * orig_edge = NULL;
577 
578  // crack tip origin coordinates and direction
579  Point crack_tip_origin(0, 0, 0);
580  Point crack_tip_direction(0, 0, 0);
581 
582  if (isElemAtCrackTip(elem)) // crack tip element's crack intiation
583  {
584  orig_cut_side_id = CEMElem->getTipEdgeID();
585  if (orig_cut_side_id < nsides) // valid crack-tip edge found
586  {
587  orig_edge = CEMElem->getEdge(orig_cut_side_id);
588  orig_node = CEMElem->getTipEmbeddedNode();
589  }
590  else
591  mooseError("element ", elem->id(), " has no valid crack-tip edge");
592 
593  // obtain the crack tip origin coordinates and direction.
594  std::map<const Elem *, std::vector<Point>>::iterator ecodm =
596  if (ecodm != _elem_crack_origin_direction_map.end())
597  {
598  crack_tip_origin = (ecodm->second)[0];
599  crack_tip_direction = (ecodm->second)[1];
600  }
601  else
602  mooseError("element ", elem->id(), " cannot find its crack tip origin and direction.");
603  }
604  else
605  {
606  std::map<const Elem *, unsigned int>::iterator mit1;
607  mit1 = _state_marked_elem_sides.find(elem);
608  std::set<const Elem *>::iterator mit2;
609  mit2 = _state_marked_frags.find(elem);
610 
611  if (mit1 != _state_marked_elem_sides.end()) // specified boundary crack initiation
612  {
613  orig_cut_side_id = mit1->second;
614  if (!CEMElem->isEdgePhantom(orig_cut_side_id) &&
615  !CEMElem->getEdge(orig_cut_side_id)->hasIntersection())
616  {
617  orig_cut_distance = 0.5;
618  _efa_mesh.addElemEdgeIntersection(elem->id(), orig_cut_side_id, orig_cut_distance);
619  orig_edge = CEMElem->getEdge(orig_cut_side_id);
620  orig_node = orig_edge->getEmbeddedNode(0);
621  // get a virtual crack tip direction
622  Point elem_center(0.0, 0.0, 0.0);
623  Point edge_center;
624  for (unsigned int i = 0; i < nsides; ++i)
625  {
626  elem_center += getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
627  elem_center += getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
628  }
629  elem_center /= nsides * 2.0;
630  edge_center = getEFANodeCoords(orig_edge->getNode(0), CEMElem, elem) +
631  getEFANodeCoords(orig_edge->getNode(1), CEMElem, elem);
632  edge_center /= 2.0;
633  crack_tip_origin = edge_center;
634  crack_tip_direction = elem_center - edge_center;
635  crack_tip_direction /= pow(crack_tip_direction.norm_sq(), 0.5);
636  }
637  else
638  continue; // skip this elem if specified boundary edge is phantom
639  }
640  else if (mit2 != _state_marked_frags.end()) // cut-surface secondary crack initiation
641  {
642  if (CEMElem->numFragments() != 1)
643  mooseError("element ",
644  elem->id(),
645  " flagged for a secondary crack, but has ",
646  CEMElem->numFragments(),
647  " fragments");
648  std::vector<unsigned int> interior_edge_id = CEMElem->getFragment(0)->getInteriorEdgeID();
649  if (interior_edge_id.size() == 1)
650  orig_cut_side_id = interior_edge_id[0];
651  else
652  continue; // skip this elem if more than one interior edges found (i.e. elem's been cut
653  // twice)
654  orig_cut_distance = 0.5;
655  _efa_mesh.addFragEdgeIntersection(elem->id(), orig_cut_side_id, orig_cut_distance);
656  orig_edge = CEMElem->getFragmentEdge(0, orig_cut_side_id);
657  orig_node = orig_edge->getEmbeddedNode(0); // must be an interior embedded node
658  Point elem_center(0.0, 0.0, 0.0);
659  Point edge_center;
660  unsigned int nsides_frag = CEMElem->getFragment(0)->numEdges();
661  for (unsigned int i = 0; i < nsides_frag; ++i)
662  {
663  elem_center +=
664  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
665  elem_center +=
666  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
667  }
668  elem_center /= nsides_frag * 2.0;
669  edge_center = getEFANodeCoords(orig_edge->getNode(0), CEMElem, elem) +
670  getEFANodeCoords(orig_edge->getNode(1), CEMElem, elem);
671  edge_center /= 2.0;
672  crack_tip_origin = edge_center;
673  crack_tip_direction = elem_center - edge_center;
674  crack_tip_direction /= pow(crack_tip_direction.norm_sq(), 0.5);
675  }
676  else
677  mooseError("element ",
678  elem->id(),
679  " flagged for state-based growth, but has no edge intersections");
680  }
681 
682  Point cut_origin(0.0, 0.0, 0.0);
683  if (orig_node)
684  cut_origin = getEFANodeCoords(orig_node, CEMElem, elem); // cutting plane origin's coords
685  else
686  mooseError("element ", elem->id(), " does not have valid orig_node");
687 
688  // loop through element edges to add possible second cut points
689  std::vector<Point> edge_ends(2, Point(0.0, 0.0, 0.0));
690  Point edge1(0.0, 0.0, 0.0);
691  Point edge2(0.0, 0.0, 0.0);
692  Point cut_edge_point(0.0, 0.0, 0.0);
693  bool find_compatible_direction = false;
694  unsigned int edge_id_keep = 0;
695  Real distance_keep = 0.0;
696  Point normal_keep(0.0, 0.0, 0.0);
697  Real distance = 0.0;
698  bool edge_cut = false;
699 
700  for (unsigned int i = 0; i < nsides; ++i)
701  {
702  if (!orig_edge->isPartialOverlap(*CEMElem->getEdge(i)))
703  {
704  edge_ends[0] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(0), CEMElem, elem);
705  edge_ends[1] = getEFANodeCoords(CEMElem->getEdge(i)->getNode(1), CEMElem, elem);
707  crack_tip_origin, normal, edge_ends[0], edge_ends[1], distance) &&
708  (!CEMElem->isEdgePhantom(i))))
709  {
710  cut_edge_point = distance * edge_ends[1] + (1.0 - distance) * edge_ends[0];
711  distance_keep = distance;
712  edge_id_keep = i;
713  normal_keep = normal;
714  edge_cut = true;
715  break;
716  }
717  }
718  }
719 
720  Point between_two_cuts = (cut_edge_point - crack_tip_origin);
721  between_two_cuts /= pow(between_two_cuts.norm_sq(), 0.5);
722  Real angle_between_two_cuts = between_two_cuts * crack_tip_direction;
723 
724  if (angle_between_two_cuts > std::cos(45.0 / 180.0 * 3.14159)) // original cut direction is good
725  find_compatible_direction = true;
726 
727  if (!find_compatible_direction && edge_cut)
729  CEMElem,
730  orig_edge,
731  normal,
732  crack_tip_origin,
733  crack_tip_direction,
734  distance_keep,
735  edge_id_keep,
736  normal_keep);
737 
738  if (edge_cut)
739  {
741  _efa_mesh.addElemEdgeIntersection(elem->id(), edge_id_keep, distance_keep);
742  else
743  {
744  MeshBase::element_iterator elem_it = _mesh->elements_begin();
745  const MeshBase::element_iterator elem_end = _mesh->elements_end();
746 
747  Point growth_direction(0.0, 0.0, 0.0);
748 
749  growth_direction(0) = -normal_keep(1);
750  growth_direction(1) = normal_keep(0);
751 
752  if (growth_direction * crack_tip_direction < 1.0e-10)
753  growth_direction *= (-1.0);
754 
755  Real x0 = crack_tip_origin(0);
756  Real y0 = crack_tip_origin(1);
757  Real x1 = x0 + _crack_growth_increment * growth_direction(0);
758  Real y1 = y0 + _crack_growth_increment * growth_direction(1);
759 
760  XFEMCrackGrowthIncrement2DCut geometric_cut(x0, y0, x1, y1, time * 0.9, time * 0.9);
761 
762  for (; elem_it != elem_end; ++elem_it)
763  {
764  const Elem * elem = *elem_it;
765  std::vector<CutEdgeForCrackGrowthIncr> elem_cut_edges;
766  EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
767  EFAElement2D * CEMElem = dynamic_cast<EFAElement2D *>(EFAelem);
768 
769  if (!CEMElem)
770  mooseError("EFAelem is not of EFAelement2D type");
771 
772  // continue if elem has been already cut twice - IMPORTANT
773  if (CEMElem->isFinalCut())
774  continue;
775 
776  // mark cut edges for the element and its fragment
777  geometric_cut.cutElementByCrackGrowthIncrement(elem, elem_cut_edges, time);
778 
779  for (unsigned int i = 0; i < elem_cut_edges.size(); ++i) // mark element edges
780  {
781  if (!CEMElem->isEdgePhantom(elem_cut_edges[i].host_side_id)) // must not be phantom edge
782  {
784  elem->id(), elem_cut_edges[i].host_side_id, elem_cut_edges[i].distance);
785  }
786  }
787  }
788  }
789  }
790  // loop though framgent boundary edges to add possible second cut points
791  // N.B. must do this after marking element edges
792  if (CEMElem->numFragments() > 0 && !edge_cut)
793  {
794  for (unsigned int i = 0; i < CEMElem->getFragment(0)->numEdges(); ++i)
795  {
796  if (!orig_edge->isPartialOverlap(*CEMElem->getFragmentEdge(0, i)))
797  {
798  edge_ends[0] =
799  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(0), CEMElem, elem);
800  edge_ends[1] =
801  getEFANodeCoords(CEMElem->getFragmentEdge(0, i)->getNode(1), CEMElem, elem);
803  crack_tip_origin, normal, edge_ends[0], edge_ends[1], distance) &&
804  (!CEMElem->getFragment(0)->isSecondaryInteriorEdge(i)))
805  {
806  if (_efa_mesh.addFragEdgeIntersection(elem->id(), edge_id_keep, distance_keep))
807  if (!isElemAtCrackTip(elem))
808  _has_secondary_cut = true;
809  break;
810  }
811  }
812  }
813  }
814 
815  marked_edges = true;
816 
817  } // loop over all state_marked_elems
818 
819  return marked_edges;
820 }
std::vector< unsigned int > getInteriorEdgeID() const
EFANode * getEmbeddedNode(unsigned int index) const
Definition: EFAEdge.C:328
unsigned int numEdges() const
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)
Definition: XFEM.C:407
bool _use_crack_growth_increment
Definition: XFEM.h:181
EFAFragment2D * getFragment(unsigned int frag_id) const
void addElemEdgeIntersection(unsigned int elemid, unsigned int edgeid, double position)
Real _crack_growth_increment
Definition: XFEM.h:182
ElementFragmentAlgorithm _efa_mesh
Definition: XFEM.h:199
virtual bool isFinalCut() const
Definition: EFAElement2D.C:789
bool _has_secondary_cut
Definition: XFEM.h:177
EFAEdge * getEdge(unsigned int edge_id) const
EFAEdge * getFragmentEdge(unsigned int frag_id, unsigned int edge_id) const
EFANode * getTipEmbeddedNode() const
bool initCutIntersectionEdge(Point cut_origin, RealVectorValue cut_normal, Point &edge_p1, Point &edge_p2, Real &dist)
Definition: XFEM.C:905
bool isSecondaryInteriorEdge(unsigned int edge_id) const
bool addFragEdgeIntersection(unsigned int elemid, unsigned int frag_edge_id, double position)
bool isEdgePhantom(unsigned int edge_id) const
Point getEFANodeCoords(EFANode *CEMnode, EFAElement *CEMElem, const Elem *elem, MeshBase *displaced_mesh=NULL) const
Definition: XFEM.C:1296
std::map< const Elem *, std::vector< Point > > _elem_crack_origin_direction_map
Definition: XFEM.h:191
EFANode * getNode(unsigned int index) const
Definition: EFAEdge.C:177
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
std::set< const Elem * > _state_marked_frags
Definition: XFEM.h:196
bool hasIntersection() const
Definition: EFAEdge.C:196
bool isPartialOverlap(const EFAEdge &other) const
Definition: EFAEdge.C:55
virtual unsigned int numFragments() const
Definition: EFAElement2D.C:200
EFAElement * getElemByID(unsigned int id)
bool isElemAtCrackTip(const Elem *elem) const
Definition: XFEM.C:1382
Real getPhysicalVolumeFraction(const Elem *elem) const
Get the volume fraction of an element that is physical.
Definition: XFEM.C:1328
unsigned int numEdges() const
std::map< const Elem *, unsigned int > _state_marked_elem_sides
Definition: XFEM.h:197
unsigned int getTipEdgeID() const
std::map< const Elem *, RealVectorValue > _state_marked_elems
Definition: XFEM.h:195
bool XFEM::markCutFacesByGeometry ( Real  time)

Definition at line 823 of file XFEM.C.

Referenced by markCuts().

824 {
825  bool marked_faces = false;
826 
827  MeshBase::element_iterator elem_it = _mesh->elements_begin();
828  const MeshBase::element_iterator elem_end = _mesh->elements_end();
829 
830  std::vector<const GeometricCutUserObject *> active_geometric_cuts;
831  for (unsigned int i = 0; i < _geometric_cuts.size(); ++i)
832  if (_geometric_cuts[i]->active(time))
833  active_geometric_cuts.push_back(_geometric_cuts[i]);
834 
835  if (active_geometric_cuts.size() > 0)
836  {
837  for (MeshBase::element_iterator elem_it = _mesh->elements_begin();
838  elem_it != _mesh->elements_end();
839  ++elem_it)
840  {
841  const Elem * elem = *elem_it;
842  std::vector<CutFace> elem_cut_faces;
843  std::vector<CutFace> frag_cut_faces;
844  std::vector<std::vector<Point>> frag_faces;
845  EFAElement * EFAelem = _efa_mesh.getElemByID(elem->id());
846  EFAElement3D * CEMElem = dynamic_cast<EFAElement3D *>(EFAelem);
847  if (!CEMElem)
848  mooseError("EFAelem is not of EFAelement3D type");
849 
850  // continue if elem has been already cut twice - IMPORTANT
851  if (CEMElem->isFinalCut())
852  continue;
853 
854  // get fragment faces
855  getFragmentFaces(elem, CEMElem, frag_faces);
856 
857  // mark cut faces for the element and its fragment
858  for (unsigned int i = 0; i < active_geometric_cuts.size(); ++i)
859  {
860  active_geometric_cuts[i]->cutElementByGeometry(elem, elem_cut_faces, time);
861  // TODO: This would be done for branching, which is not yet supported in 3D
862  // if (CEMElem->numFragments() > 0)
863  // active_geometric_cuts[i]->cutFragmentByGeometry(frag_faces, frag_cut_faces, time);
864  }
865 
866  for (unsigned int i = 0; i < elem_cut_faces.size(); ++i) // mark element faces
867  {
868  if (!CEMElem->isFacePhantom(elem_cut_faces[i].face_id)) // must not be phantom face
869  {
871  elem_cut_faces[i].face_id,
872  elem_cut_faces[i].face_edge,
873  elem_cut_faces[i].position);
874  marked_faces = true;
875  }
876  }
877 
878  for (unsigned int i = 0; i < frag_cut_faces.size();
879  ++i) // MUST DO THIS AFTER MARKING ELEMENT EDGES
880  {
881  if (!CEMElem->getFragment(0)->isThirdInteriorFace(frag_cut_faces[i].face_id))
882  {
884  frag_cut_faces[i].face_id,
885  frag_cut_faces[i].face_edge,
886  frag_cut_faces[i].position);
887  marked_faces = true;
888  }
889  }
890  }
891  }
892 
893  return marked_faces;
894 }
std::vector< const GeometricCutUserObject * > _geometric_cuts
Definition: XFEM.h:184
void getFragmentFaces(const Elem *elem, std::vector< std::vector< Point >> &frag_faces, bool displaced_mesh=false) const
Definition: XFEM.C:1412
EFAFragment3D * getFragment(unsigned int frag_id) const
bool isFacePhantom(unsigned int face_id) const
ElementFragmentAlgorithm _efa_mesh
Definition: XFEM.h:199
bool isThirdInteriorFace(unsigned int face_id) const
void addFragFaceIntersection(unsigned int ElemID, unsigned int FragFaceID, std::vector< unsigned int > FragFaceEdgeID, std::vector< double > position)
void addElemFaceIntersection(unsigned int elemid, unsigned int faceid, std::vector< unsigned int > edgeid, std::vector< double > position)
virtual bool isFinalCut() const
Definition: EFAElement3D.C:777
EFAElement * getElemByID(unsigned int id)
bool XFEM::markCutFacesByState ( )

Definition at line 897 of file XFEM.C.

Referenced by markCuts().

898 {
899  bool marked_faces = false;
900  // TODO: need to finish this for 3D problems
901  return marked_faces;
902 }
bool XFEM::markCuts ( Real  time)

Definition at line 308 of file XFEM.C.

Referenced by update().

309 {
310  bool marked_sides = false;
311  if (_mesh->mesh_dimension() == 2)
312  {
313  marked_sides = markCutEdgesByGeometry(time);
314  marked_sides |= markCutEdgesByState(time);
315  }
316  else if (_mesh->mesh_dimension() == 3)
317  {
318  marked_sides = markCutFacesByGeometry(time);
319  marked_sides |= markCutFacesByState();
320  }
321  return marked_sides;
322 }
bool markCutEdgesByState(Real time)
Definition: XFEM.C:548
bool markCutFacesByGeometry(Real time)
Definition: XFEM.C:823
bool markCutEdgesByGeometry(Real time)
Definition: XFEM.C:325
bool markCutFacesByState()
Definition: XFEM.C:897
void XFEM::setCrackGrowthMethod ( bool  use_crack_growth_increment,
Real  crack_growth_increment 
)

Definition at line 1489 of file XFEM.C.

1490 {
1491  _use_crack_growth_increment = use_crack_growth_increment;
1492  _crack_growth_increment = crack_growth_increment;
1493 }
bool _use_crack_growth_increment
Definition: XFEM.h:181
Real _crack_growth_increment
Definition: XFEM.h:182
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
sysSystem for which the solution is set
stored_solutionData structure that the stored solution is obtained from
current_solutionCurrent solution vector that will be set
old_solutionOld solution vector that will be set
older_solutionOlder solution vector that will be set

Definition at line 1671 of file XFEM.C.

Referenced by initSolution().

1676 {
1677  const auto nodes_end = _mesh->local_nodes_end();
1678  for (auto node_it = _mesh->local_nodes_begin(); node_it != nodes_end; ++node_it)
1679  {
1680  Node * node = *node_it;
1681  auto mit = stored_solution.find(node->unique_id());
1682  if (mit != stored_solution.end())
1683  {
1684  const std::vector<Real> & stored_node_solution = mit->second;
1685  std::vector<dof_id_type> stored_solution_dofs = getNodeSolutionDofs(node, sys);
1686  setSolutionForDOFs(stored_node_solution,
1687  stored_solution_dofs,
1688  current_solution,
1689  old_solution,
1690  older_solution);
1691  }
1692  }
1693 
1694  const auto elems_end = _mesh->local_elements_end();
1695  for (auto elem_it = _mesh->local_elements_begin(); elem_it != elems_end; ++elem_it)
1696  {
1697  Elem * elem = *elem_it;
1698  auto mit = stored_solution.find(elem->unique_id());
1699  if (mit != stored_solution.end())
1700  {
1701  const std::vector<Real> & stored_elem_solution = mit->second;
1702  std::vector<dof_id_type> stored_solution_dofs = getElementSolutionDofs(elem, sys);
1703  setSolutionForDOFs(stored_elem_solution,
1704  stored_solution_dofs,
1705  current_solution,
1706  old_solution,
1707  older_solution);
1708  }
1709  }
1710 }
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.
Definition: XFEM.C:1713
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...
Definition: XFEM.C:1762
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...
Definition: XFEM.C:1736
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_solutionStored solution values to set the solution to
stored_solution_dofsDof indices for the entries in stored_solution
current_solutionCurrent solution vector that will be set
old_solutionOld solution vector that will be set
older_solutionOlder solution vector that will be set

Definition at line 1713 of file XFEM.C.

Referenced by setSolution().

1718 {
1719  // Solution vector is stored first for current, then old and older solutions.
1720  // These are the offsets to the beginning of the old and older solutions in the vector.
1721  const auto old_solution_offset = stored_solution_dofs.size();
1722  const auto older_solution_offset = old_solution_offset * 2;
1723 
1724  for (std::size_t i = 0; i < stored_solution_dofs.size(); ++i)
1725  {
1726  current_solution.set(stored_solution_dofs[i], stored_solution[i]);
1727  if (_fe_problem->isTransient())
1728  {
1729  old_solution.set(stored_solution_dofs[i], stored_solution[old_solution_offset + i]);
1730  older_solution.set(stored_solution_dofs[i], stored_solution[older_solution_offset + i]);
1731  }
1732  }
1733 }
void XFEM::setXFEMCutData ( std::vector< Real > &  cut_data)
void XFEM::setXFEMCutType ( std::string &  cut_type)
void XFEM::setXFEMQRule ( std::string &  xfem_qrule)

Definition at line 1478 of file XFEM.C.

1479 {
1480  if (xfem_qrule == "volfrac")
1482  else if (xfem_qrule == "moment_fitting")
1484  else if (xfem_qrule == "direct")
1486 }
Xfem::XFEM_QRULE _XFEM_qrule
Definition: XFEM.h:179
void XFEM::storeCrackTipOriginAndDirection ( )

Definition at line 151 of file XFEM.C.

Referenced by update().

152 {
154  std::set<EFAElement *> CrackTipElements = _efa_mesh.getCrackTipElements();
155  std::set<EFAElement *>::iterator sit;
156  for (sit = CrackTipElements.begin(); sit != CrackTipElements.end(); ++sit)
157  {
158  if (_mesh->mesh_dimension() == 2)
159  {
160  EFAElement2D * CEMElem = dynamic_cast<EFAElement2D *>(*sit);
161  EFANode * tip_node = CEMElem->getTipEmbeddedNode();
162  unsigned int cts_id = CEMElem->getCrackTipSplitElementID();
163 
164  Point origin(0, 0, 0);
165  Point direction(0, 0, 0);
166 
167  std::map<unique_id_type, XFEMCutElem *>::const_iterator it;
168  it = _cut_elem_map.find(_mesh->elem(cts_id)->unique_id());
169  if (it != _cut_elem_map.end())
170  {
171  const XFEMCutElem * xfce = it->second;
172  const EFAElement * EFAelem = xfce->getEFAElement();
173  if (EFAelem->isPartial()) // exclude the full crack tip elements
174  {
175  xfce->getCrackTipOriginAndDirection(tip_node->id(), origin, direction);
176  }
177  }
178 
179  std::vector<Point> tip_data;
180  tip_data.push_back(origin);
181  tip_data.push_back(direction);
182  const Elem * elem = _mesh->elem((*sit)->id());
184  std::pair<const Elem *, std::vector<Point>>(elem, tip_data));
185  }
186  }
187 }
virtual void getCrackTipOriginAndDirection(unsigned tip_id, Point &origin, Point &direction) const =0
virtual bool isPartial() const =0
std::map< unique_id_type, XFEMCutElem * > _cut_elem_map
Definition: XFEM.h:186
ElementFragmentAlgorithm _efa_mesh
Definition: XFEM.h:199
EFANode * getTipEmbeddedNode() const
unsigned int getCrackTipSplitElementID() const
Definition: EFAElement2D.C:592
std::map< const Elem *, std::vector< Point > > _elem_crack_origin_direction_map
Definition: XFEM.h:191
unsigned int id() const
Definition: EFANode.C:34
virtual const EFAElement * getEFAElement() const =0
const std::set< EFAElement * > & getCrackTipElements()
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_toElement for which the solution will be stored
elem_to_store_fromElement from which the solution to be stored is obtained
sysSystem from which the solution is stored
stored_solutionData structure that the stored solution is saved to
current_solutionCurrent solution vector that the solution is obtained from
old_solutionOld solution vector that the solution is obtained from
older_solutionOlder solution vector that the solution is obtained from

Definition at line 1637 of file XFEM.C.

Referenced by cutMeshWithEFA().

1644 {
1645  std::vector<dof_id_type> stored_solution_dofs = getElementSolutionDofs(elem_to_store_from, sys);
1646  std::vector<Real> stored_solution_scratch;
1647  // Size for current solution, as well as for old, and older solution only for transient case
1648  std::size_t stored_solution_size =
1649  (_fe_problem->isTransient() ? stored_solution_dofs.size() * 3 : stored_solution_dofs.size());
1650  stored_solution_scratch.reserve(stored_solution_size);
1651 
1652  // Store in the order defined in stored_solution_dofs first for the current, then for old and
1653  // older if applicable
1654  for (auto dof : stored_solution_dofs)
1655  stored_solution_scratch.push_back(current_solution(dof));
1656 
1657  if (_fe_problem->isTransient())
1658  {
1659  for (auto dof : stored_solution_dofs)
1660  stored_solution_scratch.push_back(old_solution(dof));
1661 
1662  for (auto dof : stored_solution_dofs)
1663  stored_solution_scratch.push_back(older_solution(dof));
1664  }
1665 
1666  if (stored_solution_scratch.size() > 0)
1667  stored_solution[elem_to_store_to->unique_id()] = stored_solution_scratch;
1668 }
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...
Definition: XFEM.C:1736
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_toNode for which the solution will be stored
node_to_store_fromNode from which the solution to be stored is obtained
sysSystem from which the solution is stored
stored_solutionData structure that the stored solution is saved to
current_solutionCurrent solution vector that the solution is obtained from
old_solutionOld solution vector that the solution is obtained from
older_solutionOlder solution vector that the solution is obtained from

Definition at line 1603 of file XFEM.C.

Referenced by cutMeshWithEFA().

1610 {
1611  std::vector<dof_id_type> stored_solution_dofs = getNodeSolutionDofs(node_to_store_from, sys);
1612  std::vector<Real> stored_solution_scratch;
1613  // Size for current solution, as well as for old, and older solution only for transient case
1614  std::size_t stored_solution_size =
1615  (_fe_problem->isTransient() ? stored_solution_dofs.size() * 3 : stored_solution_dofs.size());
1616  stored_solution_scratch.reserve(stored_solution_size);
1617 
1618  // Store in the order defined in stored_solution_dofs first for the current, then for old and
1619  // older if applicable
1620  for (auto dof : stored_solution_dofs)
1621  stored_solution_scratch.push_back(current_solution(dof));
1622 
1623  if (_fe_problem->isTransient())
1624  {
1625  for (auto dof : stored_solution_dofs)
1626  stored_solution_scratch.push_back(old_solution(dof));
1627 
1628  for (auto dof : stored_solution_dofs)
1629  stored_solution_scratch.push_back(older_solution(dof));
1630  }
1631 
1632  if (stored_solution_scratch.size() > 0)
1633  stored_solution[node_to_store_to->unique_id()] = stored_solution_scratch;
1634 }
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...
Definition: XFEM.C:1762
bool XFEM::update ( Real  time,
NonlinearSystemBase &  nl,
AuxiliarySystem &  aux 
)
virtual

Method to update the mesh due to modified cut planes.

Definition at line 190 of file XFEM.C.

191 {
192  bool mesh_changed = false;
193 
194  buildEFAMesh();
195 
197 
198  if (markCuts(time))
199  mesh_changed = cutMeshWithEFA(nl, aux);
200 
201  if (mesh_changed)
202  {
203  buildEFAMesh();
205  }
206 
207  if (mesh_changed)
208  {
209  _mesh->update_parallel_id_counts();
210  MeshCommunication().make_elems_parallel_consistent(*_mesh);
211  MeshCommunication().make_nodes_parallel_consistent(*_mesh);
212  // _mesh->find_neighbors();
213  // _mesh->contract();
214  _mesh->allow_renumbering(false);
215  _mesh->skip_partitioning(true);
216  _mesh->prepare_for_use();
217  // _mesh->prepare_for_use(true,true); //doing this preserves the numbering, but generates
218  // warning
219 
220  if (_displaced_mesh)
221  {
222  _displaced_mesh->update_parallel_id_counts();
223  MeshCommunication().make_elems_parallel_consistent(*_displaced_mesh);
224  MeshCommunication().make_nodes_parallel_consistent(*_displaced_mesh);
225  _displaced_mesh->allow_renumbering(false);
226  _displaced_mesh->skip_partitioning(true);
227  _displaced_mesh->prepare_for_use();
228  // _displaced_mesh->prepare_for_use(true,true);
229  }
230  }
231 
233 
234  return mesh_changed;
235 }
void storeCrackTipOriginAndDirection()
Definition: XFEM.C:151
bool markCuts(Real time)
Definition: XFEM.C:308
bool cutMeshWithEFA(NonlinearSystemBase &nl, AuxiliarySystem &aux)
Definition: XFEM.C:924
void clearStateMarkedElems()
Definition: XFEM.C:143
void buildEFAMesh()
Definition: XFEM.C:265

Member Data Documentation

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 217 of file XFEM.h.

Referenced by cutMeshWithEFA(), and initSolution().

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 208 of file XFEM.h.

Referenced by cutMeshWithEFA(), and initSolution().

Real XFEM::_crack_growth_increment
private

Definition at line 182 of file XFEM.h.

Referenced by markCutEdgesByState(), and setCrackGrowthMethod().

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

Definition at line 187 of file XFEM.h.

Referenced by cutMeshWithEFA(), and isElemAtCrackTip().

std::map<unique_id_type, XFEMCutElem *> XFEM::_cut_elem_map
private
ElementFragmentAlgorithm XFEM::_efa_mesh
private
std::map<const Elem *, std::vector<Point> > XFEM::_elem_crack_origin_direction_map
private
std::vector<const GeometricCutUserObject *> XFEM::_geometric_cuts
private

Definition at line 184 of file XFEM.h.

Referenced by addGeometricCut(), markCutEdgesByGeometry(), and markCutFacesByGeometry().

bool XFEM::_has_secondary_cut
private

Definition at line 177 of file XFEM.h.

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

ElementPairLocator::ElementPairList XFEM::_sibling_displaced_elems
private

Definition at line 189 of file XFEM.h.

Referenced by cutMeshWithEFA().

ElementPairLocator::ElementPairList XFEM::_sibling_elems
private

Definition at line 188 of file XFEM.h.

Referenced by cutMeshWithEFA().

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

Definition at line 197 of file XFEM.h.

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

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

Definition at line 195 of file XFEM.h.

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

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

Definition at line 196 of file XFEM.h.

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

bool XFEM::_use_crack_growth_increment
private

Definition at line 181 of file XFEM.h.

Referenced by markCutEdgesByState(), and setCrackGrowthMethod().

Xfem::XFEM_QRULE XFEM::_XFEM_qrule
private

Definition at line 179 of file XFEM.h.

Referenced by getXFEMQRule(), and setXFEMQRule().


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