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

MeshCut3DUserObject: (1) reads in a mesh describing the crack surface, (2) uses the mesh to do initial cutting of 3D elements, and (3) grows the mesh based on prescribed growth functions. More...

#include <MeshCut3DUserObject.h>

Inheritance diagram for MeshCut3DUserObject:
[legend]

Public Member Functions

 MeshCut3DUserObject (const InputParameters &parameters)
 
virtual void initialize () override
 
virtual const std::vector< Point > getCrackFrontPoints (unsigned int num_crack_front_points) const override
 get a set of points along a crack front from a XFEM GeometricCutUserObject More...
 
virtual bool cutElementByGeometry (const Elem *elem, std::vector< Xfem::CutEdge > &cut_edges, std::vector< Xfem::CutNode > &cut_nodes, Real time) const override
 Check to see whether a specified 2D element should be cut based on geometric conditions. More...
 
virtual bool cutElementByGeometry (const Elem *elem, std::vector< Xfem::CutFace > &cut_faces, Real time) const override
 Check to see whether a specified 3D element should be cut based on geometric conditions. More...
 
virtual bool cutFragmentByGeometry (std::vector< std::vector< Point >> &frag_edges, std::vector< Xfem::CutEdge > &cut_edges, Real time) const override
 Check to see whether a fragment of a 2D element should be cut based on geometric conditions. More...
 
virtual bool cutFragmentByGeometry (std::vector< std::vector< Point >> &frag_faces, std::vector< Xfem::CutFace > &cut_faces, Real time) const override
 Check to see whether a fragment of a 3D element should be cut based on geometric conditions. More...
 
virtual void execute () override
 
virtual void threadJoin (const UserObject &y) override
 
virtual void finalize () override
 
unsigned int getInterfaceID () const
 Get the interface ID for this cutting object. More...
 
void setInterfaceID (unsigned int interface_id)
 Set the interface ID for this cutting object. More...
 
bool shouldHealMesh () const
 Should the elements cut by this cutting object be healed in the current time step? More...
 

Protected Member Functions

virtual bool intersectWithEdge (const Point &p1, const Point &p2, const std::vector< Point > &_vertices, Point &pint) const
 Check if a line intersects with an element. More...
 
bool findIntersection (const Point &p1, const Point &p2, const std::vector< Point > &vertices, Point &pint) const
 Find directional intersection along the positive extension of the vector from p1 to p2. More...
 
bool isInsideEdge (const Point &p1, const Point &p2, const Point &p) const
 Check if point p is inside the edge p1-p2. More...
 
Real getRelativePosition (const Point &p1, const Point &p2, const Point &p) const
 Get the relative position of p from p1. More...
 
bool isInsideCutPlane (const std::vector< Point > &_vertices, const Point &p) const
 Check if point p is inside a plane. More...
 
void findBoundaryNodes ()
 Find boundary nodes of the cutter mesh This is a simple algorithm simply based on the added angle = 360 degrees Works fine for planar cutting surface for curved cutting surface, need to re-work this subroutine to make it more general. More...
 
void findBoundaryEdges ()
 Find boundary edges of the cutter mesh. More...
 
void sortBoundaryNodes ()
 Sort boundary nodes to be in the right order along the boundary. More...
 
Real findDistance (dof_id_type node1, dof_id_type node2)
 Find distance between two nodes. More...
 
void refineBoundary ()
 If boundary nodes are too sparse, add nodes in between. More...
 
void findActiveBoundaryNodes ()
 Find all active boundary nodes in the cutter mesh Find boundary nodes that will grow; nodes outside of the structural mesh are inactive. More...
 
void findActiveBoundaryDirection ()
 Find growth direction at each active node. More...
 
void growFront ()
 Grow the cutter mesh. More...
 
void sortFrontNodes ()
 Sort the front nodes. More...
 
void findFrontIntersection ()
 Find front-structure intersections. More...
 
void refineFront ()
 Refine the mesh at the front. More...
 
void triangulation ()
 Create tri3 elements between the new front and the old front. More...
 
void joinBoundary ()
 Join active boundaries and inactive boundaries to be the new boundary. More...
 
void serialize (std::string &serialized_buffer)
 Methods to pack/unpack the _marked_elems_2d and _marked_elems_3d data into a structure suitable for parallel communication. More...
 
void deserialize (std::vector< std::string > &serialized_buffers)
 

Protected Attributes

std::unique_ptr< MeshBase > _cut_mesh
 The cutter mesh. More...
 
const unsigned int _cut_elem_nnode = 3
 The cutter mesh has triangluar elements only. More...
 
const unsigned int _cut_elem_dim = 2
 
MooseMesh & _mesh
 The structural mesh. More...
 
const unsigned int _elem_dim = 3
 The structural mesh must be 3D only. More...
 
const Real _const_intersection = 0.01
 Used to define intersection points. More...
 
Real _size_control
 Used for cutter mesh refinement and front advancement. More...
 
unsigned int _n_step_growth
 Number of steps to grow the mesh. More...
 
bool _stop
 variables to help control the work flow More...
 
bool _grow
 
std::vector< dof_id_type > _boundary
 Boundary nodes of the cutter mesh. More...
 
std::set< Xfem::CutEdge_boundary_edges
 Edges at the boundary. More...
 
std::map< dof_id_type, std::vector< dof_id_type > > _boundary_map
 A map of boundary nodes and their neighbors. More...
 
std::vector< std::vector< dof_id_type > > _active_boundary
 Active boundary nodes where growth is allowed. More...
 
std::vector< std::vector< Point > > _active_direction
 Growth direction for active boundaries. More...
 
std::vector< unsigned int > _inactive_boundary_pos
 Inactive boundary. More...
 
std::vector< std::vector< dof_id_type > > _front
 New boundary after growth. More...
 
Function & _func_x
 Parsed functions of front growth. More...
 
Function & _func_y
 
Function & _func_z
 
std::shared_ptr< XFEM_xfem
 Pointer to the XFEM controller object. More...
 
unsigned int _interface_id
 Associated interface id. More...
 
bool _heal_always
 Heal the mesh. More...
 
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo2D > > _marked_elems_2d
 Containers with information about all 2D and 3D elements marked for cutting by this object. More...
 
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo3D > > _marked_elems_3d
 

Detailed Description

MeshCut3DUserObject: (1) reads in a mesh describing the crack surface, (2) uses the mesh to do initial cutting of 3D elements, and (3) grows the mesh based on prescribed growth functions.

Definition at line 29 of file MeshCut3DUserObject.h.

Constructor & Destructor Documentation

◆ MeshCut3DUserObject()

MeshCut3DUserObject::MeshCut3DUserObject ( const InputParameters &  parameters)

Definition at line 44 of file MeshCut3DUserObject.C.

45  : GeometricCutUserObject(parameters),
46  _mesh(_subproblem.mesh()),
47  _n_step_growth(getParam<unsigned int>("n_step_growth")),
48  _func_x(getFunction("function_x")),
49  _func_y(getFunction("function_y")),
50  _func_z(getFunction("function_z"))
51 {
52  _grow = (_n_step_growth == 0 ? 0 : 1);
53 
54  if (_grow)
55  {
56  if (!isParamValid("size_control"))
57  mooseError("Crack growth needs size control");
58 
59  _size_control = getParam<Real>("size_control");
60  }
61 
62  // only the xda type is currently supported
63  MeshFileName xfem_cut_mesh_file = getParam<MeshFileName>("mesh_file");
64  _cut_mesh = libmesh_make_unique<ReplicatedMesh>(_communicator);
65  _cut_mesh->read(xfem_cut_mesh_file);
66 
67  // test element type; only tri3 elements are allowed
68  for (const auto & cut_elem : _cut_mesh->element_ptr_range())
69  {
70  if (cut_elem->n_nodes() != _cut_elem_nnode)
71  mooseError("The input cut mesh should include tri elements only!");
72  if (cut_elem->dim() != _cut_elem_dim)
73  mooseError("The input cut mesh should have 2D elements only!");
74  }
75 }
const unsigned int _cut_elem_nnode
The cutter mesh has triangluar elements only.
unsigned int _n_step_growth
Number of steps to grow the mesh.
Real _size_control
Used for cutter mesh refinement and front advancement.
MooseMesh & _mesh
The structural mesh.
GeometricCutUserObject(const InputParameters &parameters)
Factory constructor, takes parameters so that all derived classes can be built using the same constru...
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
Function & _func_x
Parsed functions of front growth.
const unsigned int _cut_elem_dim

Member Function Documentation

◆ cutElementByGeometry() [1/2]

bool MeshCut3DUserObject::cutElementByGeometry ( const Elem *  elem,
std::vector< Xfem::CutEdge > &  cut_edges,
std::vector< Xfem::CutNode > &  cut_nodes,
Real  time 
) const
overridevirtual

Check to see whether a specified 2D element should be cut based on geometric conditions.

Parameters
elemPointer to the libMesh element to be considered for cutting
cut_edgesData structure filled with information about edges to be cut
cut_nodesData structure filled with information about nodes to be cut
timeCurrent simulation time
Returns
bool true if element is to be cut

Implements GeometricCutUserObject.

Definition at line 112 of file MeshCut3DUserObject.C.

116 {
117  mooseError("invalid method for 3D mesh cutting");
118  return false;
119 }

◆ cutElementByGeometry() [2/2]

bool MeshCut3DUserObject::cutElementByGeometry ( const Elem *  elem,
std::vector< Xfem::CutFace > &  cut_faces,
Real  time 
) const
overridevirtual

Check to see whether a specified 3D element should be cut based on geometric conditions.

Parameters
elemPointer to the libMesh element to be considered for cutting
cut_facesData structure filled with information about edges to be cut
timeCurrent simulation time
Returns
bool true if element is to be cut

Implements GeometricCutUserObject.

Definition at line 122 of file MeshCut3DUserObject.C.

128 {
129  bool elem_cut = false;
130 
131  if (elem->dim() != _elem_dim)
132  mooseError("The structural mesh to be cut by a surface mesh must be 3D!");
133 
134  for (unsigned int i = 0; i < elem->n_sides(); ++i)
135  {
136  // This returns the lowest-order type of side.
137  std::unique_ptr<Elem> curr_side = elem->side(i);
138  if (curr_side->dim() != 2)
139  mooseError("In cutElementByGeometry dimension of side must be 2, but it is ",
140  curr_side->dim());
141  unsigned int n_edges = curr_side->n_sides();
142 
143  std::vector<unsigned int> cut_edges;
144  std::vector<Real> cut_pos;
145 
146  for (unsigned int j = 0; j < n_edges; j++)
147  {
148  // This returns the lowest-order type of side.
149  std::unique_ptr<Elem> curr_edge = curr_side->side(j);
150  if (curr_edge->type() != EDGE2)
151  mooseError("In cutElementByGeometry face edge must be EDGE2, but type is: ",
152  libMesh::Utility::enum_to_string(curr_edge->type()),
153  " base element type is: ",
154  libMesh::Utility::enum_to_string(elem->type()));
155  Node * node1 = curr_edge->get_node(0);
156  Node * node2 = curr_edge->get_node(1);
157 
158  for (const auto & cut_elem : _cut_mesh->element_ptr_range())
159  {
160  std::vector<Point> vertices;
161 
162  for (auto & node : cut_elem->node_ref_range())
163  {
164  Point & this_point = node;
165  vertices.push_back(this_point);
166  }
167 
168  Point intersection;
169  if (intersectWithEdge(*node1, *node2, vertices, intersection))
170  {
171  cut_edges.push_back(j);
172  cut_pos.emplace_back(getRelativePosition(*node1, *node2, intersection));
173  }
174  }
175  }
176 
177  // if two edges of an element are cut, it is considered as an element being cut
178  if (cut_edges.size() == 2)
179  {
180  elem_cut = true;
181  Xfem::CutFace mycut;
182  mycut._face_id = i;
183  mycut._face_edge.push_back(cut_edges[0]);
184  mycut._face_edge.push_back(cut_edges[1]);
185  mycut._position.push_back(cut_pos[0]);
186  mycut._position.push_back(cut_pos[1]);
187  cut_faces.push_back(mycut);
188  }
189  }
190  return elem_cut;
191 }
virtual bool intersectWithEdge(const Point &p1, const Point &p2, const std::vector< Point > &_vertices, Point &pint) const
Check if a line intersects with an element.
const unsigned int _elem_dim
The structural mesh must be 3D only.
std::vector< Real > _position
Fractional distance along the cut edges where the cut is located.
Data structure defining a cut through a face.
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
Real getRelativePosition(const Point &p1, const Point &p2, const Point &p) const
Get the relative position of p from p1.
std::vector< unsigned int > _face_edge
IDs of all cut faces.
unsigned int _face_id
ID of the cut face.

◆ cutFragmentByGeometry() [1/2]

bool MeshCut3DUserObject::cutFragmentByGeometry ( std::vector< std::vector< Point >> &  frag_edges,
std::vector< Xfem::CutEdge > &  cut_edges,
Real  time 
) const
overridevirtual

Check to see whether a fragment of a 2D element should be cut based on geometric conditions.

Parameters
frag_edgesData structure defining the current fragment to be considered
cut_edgesData structure filled with information about fragment edges to be cut
timeCurrent simulation time
Returns
bool true if fragment is to be cut

Implements GeometricCutUserObject.

Definition at line 194 of file MeshCut3DUserObject.C.

197 {
198  mooseError("invalid method for 3D mesh cutting");
199  return false;
200 }

◆ cutFragmentByGeometry() [2/2]

bool MeshCut3DUserObject::cutFragmentByGeometry ( std::vector< std::vector< Point >> &  frag_faces,
std::vector< Xfem::CutFace > &  cut_faces,
Real  time 
) const
overridevirtual

Check to see whether a fragment of a 3D element should be cut based on geometric conditions.

Parameters
frag_facesData structure defining the current fragment to be considered
cut_facesData structure filled with information about fragment faces to be cut
timeCurrent simulation time
Returns
bool true if fragment is to be cut

Implements GeometricCutUserObject.

Definition at line 203 of file MeshCut3DUserObject.C.

206 {
207  // TODO: Need this for branching in 3D
208  mooseError("cutFragmentByGeometry not yet implemented for 3D mesh cutting");
209  return false;
210 }

◆ deserialize()

void GeometricCutUserObject::deserialize ( std::vector< std::string > &  serialized_buffers)
protectedinherited

Definition at line 222 of file GeometricCutUserObject.C.

Referenced by GeometricCutUserObject::finalize().

223 {
224  mooseAssert(serialized_buffers.size() == _app.n_processors(),
225  "Unexpected size of serialized_buffers: " << serialized_buffers.size());
226 
227  // The input string stream used for deserialization
228  std::istringstream iss;
229 
230  // Loop over all datastructures for all processors to perfrom the gather operation
231  for (unsigned int rank = 0; rank < serialized_buffers.size(); ++rank)
232  {
233  // skip the current processor (its data is already in the structures)
234  if (rank == processor_id())
235  continue;
236 
237  // populate the stream with a new buffer and reset stream state
238  iss.clear();
239  iss.str(serialized_buffers[rank]);
240 
241  // Load the communicated data into temporary structures
242  std::map<unsigned int, std::vector<Xfem::GeomMarkedElemInfo2D>> other_marked_elems_2d;
243  std::map<unsigned int, std::vector<Xfem::GeomMarkedElemInfo3D>> other_marked_elems_3d;
244  dataLoad(iss, other_marked_elems_2d, this);
245  dataLoad(iss, other_marked_elems_3d, this);
246 
247  // merge the data in with the current processor's data
248  _marked_elems_2d.insert(other_marked_elems_2d.begin(), other_marked_elems_2d.end());
249  _marked_elems_3d.insert(other_marked_elems_3d.begin(), other_marked_elems_3d.end());
250  }
251 }
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo2D > > _marked_elems_2d
Containers with information about all 2D and 3D elements marked for cutting by this object...
void dataLoad(std::istream &stream, Xfem::CutFace &cf, void *context)
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo3D > > _marked_elems_3d

◆ execute()

void GeometricCutUserObject::execute ( )
overridevirtualinherited

Reimplemented in MovingLineSegmentCutSetUserObject.

Definition at line 67 of file GeometricCutUserObject.C.

Referenced by MovingLineSegmentCutSetUserObject::execute().

68 {
69  if (_current_elem->dim() == 2)
70  {
71  std::vector<Xfem::CutEdge> elem_cut_edges;
72  std::vector<Xfem::CutNode> elem_cut_nodes;
73  std::vector<Xfem::CutEdge> frag_cut_edges;
74  std::vector<std::vector<Point>> frag_edges;
75 
76  EFAElement2D * EFAElem = _xfem->getEFAElem2D(_current_elem);
77 
78  // Don't cut again if elem has been already cut twice
79  if (!EFAElem->isFinalCut())
80  {
81  // get fragment edges
82  _xfem->getFragmentEdges(_current_elem, EFAElem, frag_edges);
83 
84  // mark cut edges for the element and its fragment
85  bool cut = cutElementByGeometry(_current_elem, elem_cut_edges, elem_cut_nodes, _t);
86  if (EFAElem->numFragments() > 0)
87  cut |= cutFragmentByGeometry(frag_edges, frag_cut_edges, _t);
88 
89  if (cut)
90  {
92  gmei2d._elem_cut_edges = elem_cut_edges;
93  gmei2d._elem_cut_nodes = elem_cut_nodes;
94  gmei2d._frag_cut_edges = frag_cut_edges;
95  gmei2d._frag_edges = frag_edges;
96  _marked_elems_2d[_current_elem->id()].push_back(gmei2d);
97  }
98  }
99  }
100  else if (_current_elem->dim() == 3)
101  {
102  std::vector<Xfem::CutFace> elem_cut_faces;
103  std::vector<Xfem::CutFace> frag_cut_faces;
104  std::vector<std::vector<Point>> frag_faces;
105 
106  EFAElement3D * EFAElem = _xfem->getEFAElem3D(_current_elem);
107 
108  // Don't cut again if elem has been already cut twice
109  if (!EFAElem->isFinalCut())
110  {
111  // get fragment edges
112  _xfem->getFragmentFaces(_current_elem, EFAElem, frag_faces);
113 
114  // mark cut faces for the element and its fragment
115  bool cut = cutElementByGeometry(_current_elem, elem_cut_faces, _t);
116  // TODO: This would be done for branching, which is not yet supported in 3D
117  // if (EFAElem->numFragments() > 0)
118  // cut |= cutFragmentByGeometry(frag_faces, frag_cut_faces, _t);
119 
120  if (cut)
121  {
123  gmei3d._elem_cut_faces = elem_cut_faces;
124  gmei3d._frag_cut_faces = frag_cut_faces;
125  gmei3d._frag_faces = frag_faces;
126  _marked_elems_3d[_current_elem->id()].push_back(gmei3d);
127  }
128  }
129  }
130 }
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo2D > > _marked_elems_2d
Containers with information about all 2D and 3D elements marked for cutting by this object...
Data structure describing geometrically described cut through 3D element.
virtual bool cutFragmentByGeometry(std::vector< std::vector< Point >> &frag_edges, std::vector< Xfem::CutEdge > &cut_edges, Real time) const =0
Check to see whether a fragment of a 2D element should be cut based on geometric conditions.
virtual unsigned int numFragments() const
Definition: EFAElement2D.C:202
virtual bool cutElementByGeometry(const Elem *elem, std::vector< Xfem::CutEdge > &cut_edges, std::vector< Xfem::CutNode > &cut_nodes, Real time) const =0
Check to see whether a specified 2D element should be cut based on geometric conditions.
std::vector< CutNode > _elem_cut_nodes
Container for data about all cut nodes in this element.
std::vector< CutEdge > _elem_cut_edges
Container for data about all cut edges in this element.
std::vector< std::vector< Point > > _frag_edges
Container for data about all cut edges in cut fragments in this element.
virtual bool isFinalCut() const
Definition: EFAElement2D.C:791
std::vector< CutEdge > _frag_cut_edges
Container for data about all cut fragments in this element.
Data structure describing geometrically described cut through 2D element.
std::shared_ptr< XFEM > _xfem
Pointer to the XFEM controller object.
std::vector< CutFace > _frag_cut_faces
Container for data about all faces this element&#39;s fragment.
std::vector< std::vector< Point > > _frag_faces
Container for data about all cut faces in cut fragments in this element.
std::vector< CutFace > _elem_cut_faces
Container for data about all cut faces in this element.
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo3D > > _marked_elems_3d
virtual bool isFinalCut() const
Definition: EFAElement3D.C:838

◆ finalize()

void GeometricCutUserObject::finalize ( )
overridevirtualinherited

Reimplemented in MovingLineSegmentCutSetUserObject.

Definition at line 254 of file GeometricCutUserObject.C.

Referenced by MovingLineSegmentCutSetUserObject::finalize().

255 {
256  // for single processor runs we do not need to do anything here
257  if (_app.n_processors() > 1)
258  {
259  // create send buffer
260  std::string send_buffer;
261 
262  // create byte buffers for the streams received from all processors
263  std::vector<std::string> recv_buffers;
264 
265  // pack the complex datastructures into the string stream
266  serialize(send_buffer);
267 
268  // broadcast serialized data to and receive from all processors
269  _communicator.allgather(send_buffer, recv_buffers);
270 
271  // unpack the received data and merge it into the local data structures
272  deserialize(recv_buffers);
273  }
274 
275  for (const auto & it : _marked_elems_2d)
276  for (const auto & gmei : it.second)
277  _xfem->addGeomMarkedElem2D(it.first, gmei, _interface_id);
278 
279  for (const auto & it : _marked_elems_3d)
280  for (const auto & gmei : it.second)
281  _xfem->addGeomMarkedElem3D(it.first, gmei, _interface_id);
282 
283  _marked_elems_2d.clear();
284  _marked_elems_3d.clear();
285 }
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo2D > > _marked_elems_2d
Containers with information about all 2D and 3D elements marked for cutting by this object...
unsigned int _interface_id
Associated interface id.
void deserialize(std::vector< std::string > &serialized_buffers)
void serialize(std::string &serialized_buffer)
Methods to pack/unpack the _marked_elems_2d and _marked_elems_3d data into a structure suitable for p...
std::shared_ptr< XFEM > _xfem
Pointer to the XFEM controller object.
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo3D > > _marked_elems_3d

◆ findActiveBoundaryDirection()

void MeshCut3DUserObject::findActiveBoundaryDirection ( )
protected

Find growth direction at each active node.

Definition at line 666 of file MeshCut3DUserObject.C.

Referenced by initialize().

670 {
671  _active_direction.clear();
672 
673  for (unsigned int i = 0; i < _active_boundary.size(); ++i)
674  {
675  std::vector<Point> temp;
676  Point dir;
677 
678  if (_inactive_boundary_pos.size() != 0)
679  {
680  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
681  dir(j) = 0;
682  temp.push_back(dir);
683  }
684 
685  unsigned int i1 = 1;
686  unsigned int i2 = _active_boundary[i].size() - 1;
687  if (_inactive_boundary_pos.size() == 0)
688  {
689  i1 = 0;
690  i2 = _active_boundary[i].size();
691  }
692 
693  for (unsigned int j = i1; j < i2; ++j)
694  {
695  Node * this_node = _cut_mesh->node_ptr(_active_boundary[i][j]);
696  mooseAssert(this_node, "Node is NULL");
697  Point & this_point = *this_node;
698 
699  dir(0) = _func_x.value(0, this_point);
700  dir(1) = _func_y.value(0, this_point);
701  dir(2) = _func_z.value(0, this_point);
702 
703  temp.push_back(dir);
704  }
705 
706  if (_inactive_boundary_pos.size() != 0)
707  {
708  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
709  dir(j) = 0;
710  temp.push_back(dir);
711  }
712 
713  _active_direction.push_back(temp);
714  }
715 
716  Real maxl = 0;
717 
718  for (unsigned int i = 0; i < _active_direction.size(); ++i)
719  for (unsigned int j = 0; j < _active_direction[i].size(); ++j)
720  {
721  Point pt = _active_direction[i][j];
722  Real length = std::sqrt(pt * pt);
723  if (length > maxl)
724  maxl = length;
725  }
726 
727  for (unsigned int i = 0; i < _active_direction.size(); ++i)
728  for (unsigned int j = 0; j < _active_direction[i].size(); ++j)
729  _active_direction[i][j] /= maxl;
730 }
std::vector< std::vector< Point > > _active_direction
Growth direction for active boundaries.
std::vector< std::vector< dof_id_type > > _active_boundary
Active boundary nodes where growth is allowed.
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
Function & _func_x
Parsed functions of front growth.
std::vector< unsigned int > _inactive_boundary_pos
Inactive boundary.

◆ findActiveBoundaryNodes()

void MeshCut3DUserObject::findActiveBoundaryNodes ( )
protected

Find all active boundary nodes in the cutter mesh Find boundary nodes that will grow; nodes outside of the structural mesh are inactive.

Definition at line 606 of file MeshCut3DUserObject.C.

Referenced by initialize().

607 {
608  _active_boundary.clear();
609  _inactive_boundary_pos.clear();
610 
611  std::unique_ptr<PointLocatorBase> pl = _mesh.getPointLocator();
612  pl->enable_out_of_mesh_mode();
613 
614  unsigned int n_boundary = _boundary.size();
615 
616  // if the node is outside of the structural model, store its position in _boundary to
617  // _inactive_boundary_pos
618  for (unsigned int j = 0; j < n_boundary; ++j)
619  {
620  Node * this_node = _cut_mesh->node_ptr(_boundary[j]);
621  mooseAssert(this_node, "Node is NULL");
622  Point & this_point = *this_node;
623 
624  const Elem * elem = (*pl)(this_point);
625  if (elem == NULL)
626  _inactive_boundary_pos.push_back(j);
627  }
628 
629  unsigned int n_inactive_boundary = _inactive_boundary_pos.size();
630 
631  // all nodes are inactive, stop
632  if (n_inactive_boundary == n_boundary)
633  _stop = 1;
634 
635  // find and store active boundary segments in "_active_boundary"
636  if (n_inactive_boundary == 0)
637  _active_boundary.push_back(_boundary);
638  else
639  {
640  for (unsigned int i = 0; i < n_inactive_boundary - 1; ++i)
641  {
642  if (_inactive_boundary_pos[i + 1] - _inactive_boundary_pos[i] != 1)
643  {
644  std::vector<dof_id_type> temp;
645  for (unsigned int j = _inactive_boundary_pos[i]; j <= _inactive_boundary_pos[i + 1]; ++j)
646  {
647  temp.push_back(_boundary[j]);
648  }
649  _active_boundary.push_back(temp);
650  }
651  }
652  if (_inactive_boundary_pos[n_inactive_boundary - 1] - _inactive_boundary_pos[0] <
653  n_boundary - 1)
654  {
655  std::vector<dof_id_type> temp;
656  for (unsigned int j = _inactive_boundary_pos[n_inactive_boundary - 1]; j < n_boundary; ++j)
657  temp.push_back(_boundary[j]);
658  for (unsigned int j = 0; j <= _inactive_boundary_pos[0]; ++j)
659  temp.push_back(_boundary[j]);
660  _active_boundary.push_back(temp);
661  }
662  }
663 }
bool _stop
variables to help control the work flow
MooseMesh & _mesh
The structural mesh.
std::vector< std::vector< dof_id_type > > _active_boundary
Active boundary nodes where growth is allowed.
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
std::vector< dof_id_type > _boundary
Boundary nodes of the cutter mesh.
std::vector< unsigned int > _inactive_boundary_pos
Inactive boundary.

◆ findBoundaryEdges()

void MeshCut3DUserObject::findBoundaryEdges ( )
protected

Find boundary edges of the cutter mesh.

Definition at line 361 of file MeshCut3DUserObject.C.

Referenced by initialize().

362 {
363  _boundary_edges.clear();
364 
365  std::vector<dof_id_type> corner_elem_id;
366  unsigned int counter = 0;
367 
368  std::vector<dof_id_type> node_id(_cut_elem_nnode);
369  std::vector<bool> is_node_on_boundary(_cut_elem_nnode);
370 
371  for (const auto & cut_elem : _cut_mesh->element_ptr_range())
372  {
373  for (unsigned int i = 0; i < _cut_elem_nnode; ++i)
374  {
375  node_id[i] = cut_elem->get_node(i)->id();
376  is_node_on_boundary[i] = (_boundary_map.find(node_id[i]) != _boundary_map.end());
377  }
378 
379  if (is_node_on_boundary[0] && is_node_on_boundary[1] && is_node_on_boundary[2])
380  {
381  // this is an element at the corner; all nodes are on the boundary but not all edges are on
382  // the boundary
383  corner_elem_id.push_back(counter);
384  }
385  else
386  {
387  // for other elements, find and store boundary edges
388  for (unsigned int i = 0; i < _cut_elem_nnode; ++i)
389  {
390  // if both nodes on an edge are on the boundary, it is a boundary edge.
391  if (is_node_on_boundary[i] && is_node_on_boundary[(i + 1 <= 2) ? i + 1 : 0])
392  {
393  dof_id_type node1 = node_id[i];
394  dof_id_type node2 = node_id[(i + 1 <= 2) ? i + 1 : 0];
395  if (node1 > node2)
396  std::swap(node1, node2);
397 
398  Xfem::CutEdge ce;
399 
400  if (node1 > node2)
401  std::swap(node1, node2);
402  ce._id1 = node1;
403  ce._id2 = node2;
404 
405  _boundary_edges.insert(ce);
406  }
407  }
408  }
409  ++counter;
410  }
411 
412  // loop over edges in corner elements
413  // if an edge is shared by two elements, it is not an boundary edge (is_edge_inside = 1)
414  for (unsigned int i = 0; i < corner_elem_id.size(); ++i)
415  {
416  auto elem_it = _cut_mesh->elements_begin();
417 
418  for (dof_id_type j = 0; j < corner_elem_id[i]; ++j)
419  ++elem_it;
420  Elem * cut_elem = *elem_it;
421 
422  for (unsigned int j = 0; j < _cut_elem_nnode; ++j)
423  {
424  bool is_edge_inside = 0;
425 
426  dof_id_type node1 = cut_elem->get_node(j)->id();
427  dof_id_type node2 = cut_elem->get_node((j + 1 <= 2) ? j + 1 : 0)->id();
428  if (node1 > node2)
429  std::swap(node1, node2);
430 
431  unsigned int counter = 0;
432  for (const auto & cut_elem2 : _cut_mesh->element_ptr_range())
433  {
434  if (counter != corner_elem_id[i])
435  {
436  for (unsigned int k = 0; k < _cut_elem_nnode; ++k)
437  {
438  dof_id_type node3 = cut_elem2->get_node(k)->id();
439  dof_id_type node4 = cut_elem2->get_node((k + 1 <= 2) ? k + 1 : 0)->id();
440  if (node3 > node4)
441  std::swap(node3, node4);
442 
443  if (node1 == node3 && node2 == node4)
444  {
445  is_edge_inside = 1;
446  goto endloop;
447  }
448  }
449  }
450  ++counter;
451  }
452  endloop:
453  if (is_edge_inside == 0)
454  {
455  // store boundary edges
456  Xfem::CutEdge ce;
457 
458  if (node1 > node2)
459  std::swap(node1, node2);
460  ce._id1 = node1;
461  ce._id2 = node2;
462 
463  _boundary_edges.insert(ce);
464  }
465  else
466  {
467  // this is not a boundary edge; remove it from existing edge list
468  for (auto it = _boundary_edges.begin(); it != _boundary_edges.end();)
469  {
470  if ((*it)._id1 == node1 && (*it)._id2 == node2)
471  it = _boundary_edges.erase(it);
472  else
473  ++it;
474  }
475  }
476  }
477  }
478 }
const unsigned int _cut_elem_nnode
The cutter mesh has triangluar elements only.
std::map< dof_id_type, std::vector< dof_id_type > > _boundary_map
A map of boundary nodes and their neighbors.
Data structure defining a cut on an element edge.
std::set< Xfem::CutEdge > _boundary_edges
Edges at the boundary.
unsigned int _id1
ID of the first node on the edge.
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
unsigned int _id2
ID of the second node on the edge.
static unsigned int counter

◆ findBoundaryNodes()

void MeshCut3DUserObject::findBoundaryNodes ( )
protected

Find boundary nodes of the cutter mesh This is a simple algorithm simply based on the added angle = 360 degrees Works fine for planar cutting surface for curved cutting surface, need to re-work this subroutine to make it more general.

Definition at line 319 of file MeshCut3DUserObject.C.

Referenced by initialize().

320 {
321  unsigned int n_nodes = _cut_mesh->n_nodes();
322  std::vector<Real> angle(n_nodes, 0); // this assumes that the cutter mesh has compressed node id
323  std::vector<dof_id_type> node_id(_cut_elem_nnode);
324 
325  std::vector<Point> vertices(_cut_elem_nnode);
326 
327  for (const auto & cut_elem : _cut_mesh->element_ptr_range())
328  {
329  for (unsigned int i = 0; i < _cut_elem_nnode; ++i)
330  {
331  Node * this_node = cut_elem->get_node(i);
332  Point & this_point = *this_node;
333  vertices[i] = this_point;
334  node_id[i] = this_node->id();
335  }
336 
337  for (unsigned int i = 0; i < _cut_elem_nnode; ++i)
338  mooseAssert(node_id[i] < n_nodes, "Node ID is out of range");
339 
340  angle[node_id[0]] += Xfem::angle_rad_3d(&vertices[2](0), &vertices[0](0), &vertices[1](0));
341  angle[node_id[1]] += Xfem::angle_rad_3d(&vertices[0](0), &vertices[1](0), &vertices[2](0));
342  angle[node_id[2]] += Xfem::angle_rad_3d(&vertices[1](0), &vertices[2](0), &vertices[0](0));
343  }
344 
345  // In each element, angles at three vertices are calculated. Angles associated with all nodes are
346  // evaluated.
347  // Interior nodes will have a total angle = 2*pi; otherwise, it is a boundary node
348  // This assumes the cutter surface is flat.
349  for (const auto & node : _cut_mesh->node_ptr_range())
350  {
351  dof_id_type id = node->id();
352  if (!MooseUtils::relativeFuzzyEqual(angle[id], libMesh::pi * 2))
353  {
354  std::vector<dof_id_type> neighbors;
355  _boundary_map[id] = neighbors;
356  }
357  }
358 }
const unsigned int _cut_elem_nnode
The cutter mesh has triangluar elements only.
std::map< dof_id_type, std::vector< dof_id_type > > _boundary_map
A map of boundary nodes and their neighbors.
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
double angle_rad_3d(double p1[3], double p2[3], double p3[3])
Definition: XFEMFuncs.C:692

◆ findDistance()

Real MeshCut3DUserObject::findDistance ( dof_id_type  node1,
dof_id_type  node2 
)
protected

Find distance between two nodes.

Definition at line 541 of file MeshCut3DUserObject.C.

Referenced by refineBoundary(), refineFront(), and triangulation().

542 {
543  Node * n1 = _cut_mesh->node_ptr(node1);
544  mooseAssert(n1 != nullptr, "Node is NULL");
545  Node * n2 = _cut_mesh->node_ptr(node2);
546  mooseAssert(n2 != nullptr, "Node is NULL");
547  Real distance = (*n1 - *n2).norm();
548  return distance;
549 }
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.

◆ findFrontIntersection()

void MeshCut3DUserObject::findFrontIntersection ( )
protected

Find front-structure intersections.

Definition at line 778 of file MeshCut3DUserObject.C.

Referenced by initialize().

779 {
780  ConstBndElemRange & range = *_mesh.getBoundaryElementRange();
781 
782  for (unsigned int i = 0; i < _front.size(); ++i)
783  {
784  if (_front[i].size() >= 2)
785  {
786  std::vector<Point> pint1;
787  std::vector<Point> pint2;
788  std::vector<Real> length1;
789  std::vector<Real> length2;
790 
791  Real node_id = _front[i][0];
792  Node * this_node = _cut_mesh->node_ptr(node_id);
793  mooseAssert(this_node, "Node is NULL");
794  Point & p2 = *this_node;
795 
796  node_id = _front[i][1];
797  this_node = _cut_mesh->node_ptr(node_id);
798  mooseAssert(this_node, "Node is NULL");
799  Point & p1 = *this_node;
800 
801  node_id = _front[i].back();
802  this_node = _cut_mesh->node_ptr(node_id);
803  mooseAssert(this_node, "Node is NULL");
804  Point & p4 = *this_node;
805 
806  node_id = _front[i][_front[i].size() - 2];
807  this_node = _cut_mesh->node_ptr(node_id);
808  mooseAssert(this_node, "Node is NULL");
809  Point & p3 = *this_node;
810 
811  bool do_inter1 = 1;
812  bool do_inter2 = 1;
813 
814  std::unique_ptr<PointLocatorBase> pl = _mesh.getPointLocator();
815  pl->enable_out_of_mesh_mode();
816  const Elem * elem = (*pl)(p1);
817  if (elem == NULL)
818  do_inter1 = 0;
819  elem = (*pl)(p4);
820  if (elem == NULL)
821  do_inter2 = 0;
822 
823  for (const auto & belem : range)
824  {
825  Point pt;
826  std::vector<Point> vertices;
827 
828  elem = belem->_elem;
829  std::unique_ptr<Elem> curr_side = elem->side(belem->_side);
830  for (unsigned int j = 0; j < curr_side->n_nodes(); ++j)
831  {
832  Node * node = curr_side->get_node(j);
833  Point & this_point = *node;
834  vertices.push_back(this_point);
835  }
836 
837  if (findIntersection(p1, p2, vertices, pt))
838  {
839  pint1.push_back(pt);
840  length1.push_back((pt - p1) * (pt - p1));
841  }
842  if (findIntersection(p3, p4, vertices, pt))
843  {
844  pint2.push_back(pt);
845  length2.push_back((pt - p3) * (pt - p3));
846  }
847  }
848 
849  if (length1.size() != 0 && do_inter1)
850  {
851  auto it1 = std::min_element(length1.begin(), length1.end());
852  Point inter1 = pint1[std::distance(length1.begin(), it1)];
853  inter1 += (inter1 - p1) * _const_intersection;
854 
855  Node * this_node = Node::build(inter1, _cut_mesh->n_nodes()).release();
856  _cut_mesh->add_node(this_node);
857 
858  mooseAssert(_cut_mesh->n_nodes() - 1 > 0, "The cut mesh should have at least one element.");
859  unsigned int n = _cut_mesh->n_nodes() - 1;
860 
861  auto it = _front[i].begin();
862  _front[i].insert(it, n);
863  }
864 
865  if (length2.size() != 0 && do_inter2)
866  {
867  auto it2 = std::min_element(length2.begin(), length2.end());
868  Point inter2 = pint2[std::distance(length2.begin(), it2)];
869  inter2 += (inter2 - p2) * _const_intersection;
870 
871  Node * this_node = Node::build(inter2, _cut_mesh->n_nodes()).release();
872  _cut_mesh->add_node(this_node);
873 
874  dof_id_type n = _cut_mesh->n_nodes() - 1;
875 
876  auto it = _front[i].begin();
877  unsigned int m = _front[i].size();
878  _front[i].insert(it + m, n);
879  }
880  }
881  }
882 }
bool findIntersection(const Point &p1, const Point &p2, const std::vector< Point > &vertices, Point &pint) const
Find directional intersection along the positive extension of the vector from p1 to p2...
const Real _const_intersection
Used to define intersection points.
MooseMesh & _mesh
The structural mesh.
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
std::vector< std::vector< dof_id_type > > _front
New boundary after growth.

◆ findIntersection()

bool MeshCut3DUserObject::findIntersection ( const Point &  p1,
const Point &  p2,
const std::vector< Point > &  vertices,
Point &  pint 
) const
protected

Find directional intersection along the positive extension of the vector from p1 to p2.

Definition at line 244 of file MeshCut3DUserObject.C.

Referenced by findFrontIntersection().

248 {
249  bool has_intersection = false;
250 
251  Plane elem_plane(vertices[0], vertices[1], vertices[2]);
252  Point point = vertices[0];
253  Point normal = elem_plane.unit_normal(point);
254 
255  std::array<Real, 3> plane_point = {{point(0), point(1), point(2)}};
256  std::array<Real, 3> planenormal = {{normal(0), normal(1), normal(2)}};
257  std::array<Real, 3> p_begin = {{p1(0), p1(1), p1(2)}};
258  std::array<Real, 3> p_end = {{p2(0), p2(1), p2(2)}};
259  std::array<Real, 3> cut_point = {{0.0, 0.0, 0.0}};
260 
262  &plane_point[0], &planenormal[0], &p_begin[0], &p_end[0], &cut_point[0]) == 1)
263  {
264  Point p(cut_point[0], cut_point[1], cut_point[2]);
265  Real dotp = ((p - p1) * (p2 - p1)) / ((p2 - p1) * (p2 - p1));
266  if (isInsideCutPlane(vertices, p) && dotp > 1)
267  {
268  pint = p;
269  has_intersection = true;
270  }
271  }
272  return has_intersection;
273 }
int plane_normal_line_exp_int_3d(double pp[3], double normal[3], double p1[3], double p2[3], double pint[3])
Definition: XFEMFuncs.C:403
bool isInsideCutPlane(const std::vector< Point > &_vertices, const Point &p) const
Check if point p is inside a plane.

◆ getCrackFrontPoints()

const std::vector< Point > MeshCut3DUserObject::getCrackFrontPoints ( unsigned int  int) const
overridevirtual

get a set of points along a crack front from a XFEM GeometricCutUserObject

Returns
A vector which contains all crack front points

Implements CrackFrontPointsProvider.

Definition at line 1073 of file MeshCut3DUserObject.C.

1074 {
1075  mooseError("getCrackFrontPoints() is not implemented for this object.");
1076 }

◆ getInterfaceID()

unsigned int GeometricCutUserObject::getInterfaceID ( ) const
inlineinherited

Get the interface ID for this cutting object.

Returns
the interface ID

Definition at line 173 of file GeometricCutUserObject.h.

173 { return _interface_id; };
unsigned int _interface_id
Associated interface id.

◆ getRelativePosition()

Real MeshCut3DUserObject::getRelativePosition ( const Point &  p1,
const Point &  p2,
const Point &  p 
) const
protected

Get the relative position of p from p1.

Definition at line 284 of file MeshCut3DUserObject.C.

285 {
286  Real full_len = (p2 - p1).norm();
287  Real len_p1_p = (p - p1).norm();
288  return len_p1_p / full_len;
289 }

◆ growFront()

void MeshCut3DUserObject::growFront ( )
protected

Grow the cutter mesh.

Definition at line 733 of file MeshCut3DUserObject.C.

Referenced by initialize().

734 {
735  _front.clear();
736 
737  for (unsigned int i = 0; i < _active_boundary.size(); ++i)
738  {
739  std::vector<dof_id_type> temp;
740 
741  unsigned int i1 = 1;
742  unsigned int i2 = _active_boundary[i].size() - 1;
743  if (_inactive_boundary_pos.size() == 0)
744  {
745  i1 = 0;
746  i2 = _active_boundary[i].size();
747  }
748 
749  for (unsigned int j = i1; j < i2; ++j)
750  {
751  Node * this_node = _cut_mesh->node_ptr(_active_boundary[i][j]);
752  mooseAssert(this_node, "Node is NULL");
753  Point & this_point = *this_node;
754  Point dir = _active_direction[i][j];
755 
756  Point x;
757  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
758  x(k) = this_point(k) + dir(k) * _size_control;
759 
760  this_node = Node::build(x, _cut_mesh->n_nodes()).release();
761  _cut_mesh->add_node(this_node);
762 
763  dof_id_type id = _cut_mesh->n_nodes() - 1;
764  temp.push_back(id);
765  }
766 
767  _front.push_back(temp);
768  }
769 }
std::vector< std::vector< Point > > _active_direction
Growth direction for active boundaries.
Real _size_control
Used for cutter mesh refinement and front advancement.
std::vector< std::vector< dof_id_type > > _active_boundary
Active boundary nodes where growth is allowed.
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
std::vector< std::vector< dof_id_type > > _front
New boundary after growth.
std::vector< unsigned int > _inactive_boundary_pos
Inactive boundary.

◆ initialize()

void MeshCut3DUserObject::initialize ( )
overridevirtual

Reimplemented from GeometricCutUserObject.

Definition at line 78 of file MeshCut3DUserObject.C.

79 {
80  // this is a test of all methods relevant to front growth
81  if (_grow)
82  {
83  _stop = 0;
84 
89 
90  for (unsigned int i = 0; i < _n_step_growth; ++i)
91  {
93 
94  if (_stop != 1)
95  {
97  growFront();
99 
100  if (_inactive_boundary_pos.size() != 0)
102 
103  refineFront();
104  triangulation();
105  joinBoundary();
106  }
107  }
108  }
109 }
void findBoundaryNodes()
Find boundary nodes of the cutter mesh This is a simple algorithm simply based on the added angle = 3...
unsigned int _n_step_growth
Number of steps to grow the mesh.
void sortFrontNodes()
Sort the front nodes.
void findActiveBoundaryDirection()
Find growth direction at each active node.
bool _stop
variables to help control the work flow
void sortBoundaryNodes()
Sort boundary nodes to be in the right order along the boundary.
void findActiveBoundaryNodes()
Find all active boundary nodes in the cutter mesh Find boundary nodes that will grow; nodes outside o...
void joinBoundary()
Join active boundaries and inactive boundaries to be the new boundary.
void refineFront()
Refine the mesh at the front.
void findBoundaryEdges()
Find boundary edges of the cutter mesh.
void findFrontIntersection()
Find front-structure intersections.
void triangulation()
Create tri3 elements between the new front and the old front.
void growFront()
Grow the cutter mesh.
void refineBoundary()
If boundary nodes are too sparse, add nodes in between.
std::vector< unsigned int > _inactive_boundary_pos
Inactive boundary.

◆ intersectWithEdge()

bool MeshCut3DUserObject::intersectWithEdge ( const Point &  p1,
const Point &  p2,
const std::vector< Point > &  _vertices,
Point &  pint 
) const
protectedvirtual

Check if a line intersects with an element.

Definition at line 213 of file MeshCut3DUserObject.C.

217 {
218  bool has_intersection = false;
219 
220  Plane elem_plane(vertices[0], vertices[1], vertices[2]);
221  Point point = vertices[0];
222  Point normal = elem_plane.unit_normal(point);
223 
224  std::array<Real, 3> plane_point = {{point(0), point(1), point(2)}};
225  std::array<Real, 3> planenormal = {{normal(0), normal(1), normal(2)}};
226  std::array<Real, 3> edge_point1 = {{p1(0), p1(1), p1(2)}};
227  std::array<Real, 3> edge_point2 = {{p2(0), p2(1), p2(2)}};
228  std::array<Real, 3> cut_point = {{0.0, 0.0, 0.0}};
229 
231  &plane_point[0], &planenormal[0], &edge_point1[0], &edge_point2[0], &cut_point[0]) == 1)
232  {
233  Point temp_p(cut_point[0], cut_point[1], cut_point[2]);
234  if (isInsideCutPlane(vertices, temp_p) && isInsideEdge(p1, p2, temp_p))
235  {
236  pint = temp_p;
237  has_intersection = true;
238  }
239  }
240  return has_intersection;
241 }
int plane_normal_line_exp_int_3d(double pp[3], double normal[3], double p1[3], double p2[3], double pint[3])
Definition: XFEMFuncs.C:403
bool isInsideCutPlane(const std::vector< Point > &_vertices, const Point &p) const
Check if point p is inside a plane.
bool isInsideEdge(const Point &p1, const Point &p2, const Point &p) const
Check if point p is inside the edge p1-p2.

◆ isInsideCutPlane()

bool MeshCut3DUserObject::isInsideCutPlane ( const std::vector< Point > &  _vertices,
const Point &  p 
) const
protected

Check if point p is inside a plane.

Definition at line 292 of file MeshCut3DUserObject.C.

Referenced by findIntersection(), and intersectWithEdge().

293 {
294  unsigned int n_node = vertices.size();
295 
296  Plane elem_plane(vertices[0], vertices[1], vertices[2]);
297  Point normal = elem_plane.unit_normal(vertices[0]);
298 
299  bool inside = false;
300  unsigned int counter = 0;
301 
302  for (unsigned int i = 0; i < n_node; ++i)
303  {
304  unsigned int iplus1 = (i < n_node - 1 ? i + 1 : 0);
305  Point middle2p = p - 0.5 * (vertices[i] + vertices[iplus1]);
306  const Point side_tang = vertices[iplus1] - vertices[i];
307  Point side_norm = side_tang.cross(normal);
308  Xfem::normalizePoint(middle2p);
309  Xfem::normalizePoint(side_norm);
310  if (middle2p * side_norm <= 0.0)
311  counter += 1;
312  }
313  if (counter == n_node)
314  inside = true;
315  return inside;
316 }
void normalizePoint(Point &p)
Definition: XFEMFuncs.C:621
static unsigned int counter

◆ isInsideEdge()

bool MeshCut3DUserObject::isInsideEdge ( const Point &  p1,
const Point &  p2,
const Point &  p 
) const
protected

Check if point p is inside the edge p1-p2.

Definition at line 276 of file MeshCut3DUserObject.C.

Referenced by intersectWithEdge().

277 {
278  Real dotp1 = (p1 - p) * (p2 - p1);
279  Real dotp2 = (p2 - p) * (p2 - p1);
280  return (dotp1 * dotp2 <= 0.0);
281 }

◆ joinBoundary()

void MeshCut3DUserObject::joinBoundary ( )
protected

Join active boundaries and inactive boundaries to be the new boundary.

Definition at line 1033 of file MeshCut3DUserObject.C.

Referenced by initialize().

1034 {
1035  if (_inactive_boundary_pos.size() == 0)
1036  {
1037  _boundary = _front[0];
1038  _boundary.pop_back();
1039  return;
1040  }
1041 
1042  std::vector<dof_id_type> full_front;
1043 
1044  unsigned int size1 = _active_boundary.size();
1045 
1046  for (unsigned int i = 0; i < size1; ++i)
1047  {
1048  unsigned int size2 = _active_boundary[i].size();
1049 
1050  dof_id_type bd1 = _active_boundary[i][size2 - 1];
1051  dof_id_type bd2 = _active_boundary[i + 1 < size1 ? i + 1 : 0][0];
1052 
1053  full_front.insert(full_front.end(), _front[i].begin(), _front[i].end());
1054 
1055  auto it1 = std::find(_boundary.begin(), _boundary.end(), bd1);
1056  unsigned int pos1 = std::distance(_boundary.begin(), it1);
1057  auto it2 = std::find(_boundary.begin(), _boundary.end(), bd2);
1058  unsigned int pos2 = std::distance(_boundary.begin(), it2);
1059 
1060  if (pos1 <= pos2)
1061  full_front.insert(full_front.end(), _boundary.begin() + pos1, _boundary.begin() + pos2 + 1);
1062  else
1063  {
1064  full_front.insert(full_front.end(), _boundary.begin() + pos1, _boundary.end());
1065  full_front.insert(full_front.end(), _boundary.begin(), _boundary.begin() + pos2 + 1);
1066  }
1067  }
1068 
1069  _boundary = full_front;
1070 }
std::vector< std::vector< dof_id_type > > _active_boundary
Active boundary nodes where growth is allowed.
std::vector< dof_id_type > _boundary
Boundary nodes of the cutter mesh.
std::vector< std::vector< dof_id_type > > _front
New boundary after growth.
std::vector< unsigned int > _inactive_boundary_pos
Inactive boundary.

◆ refineBoundary()

void MeshCut3DUserObject::refineBoundary ( )
protected

If boundary nodes are too sparse, add nodes in between.

Definition at line 552 of file MeshCut3DUserObject.C.

Referenced by initialize().

553 {
554  std::vector<dof_id_type> new_boundary_order(_boundary.begin(), _boundary.end());
555 
556  mooseAssert(_boundary.size() >= 2, "Boundary should have at least two nodes");
557 
558  for (unsigned int i = _boundary.size() - 1; i >= 1; --i)
559  {
560  dof_id_type node1 = _boundary[i - 1];
561  dof_id_type node2 = _boundary[i];
562 
563  Real distance = findDistance(node1, node2);
564 
565  if (distance > _size_control)
566  {
567  unsigned int n = static_cast<unsigned int>(distance / _size_control);
568  std::array<Real, LIBMESH_DIM> x1;
569  std::array<Real, LIBMESH_DIM> x2;
570 
571  Node * n1 = _cut_mesh->node_ptr(node1);
572  mooseAssert(n1 != nullptr, "Node is NULL");
573  Point & p1 = *n1;
574  Node * n2 = _cut_mesh->node_ptr(node2);
575  mooseAssert(n2 != nullptr, "Node is NULL");
576  Point & p2 = *n2;
577 
578  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
579  {
580  x1[j] = p1(j);
581  x2[j] = p2(j);
582  }
583 
584  for (unsigned int j = 0; j < n; ++j)
585  {
586  Point x;
587  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
588  x(k) = x2[k] - (x2[k] - x1[k]) * (j + 1) / (n + 1);
589 
590  Node * this_node = Node::build(x, _cut_mesh->n_nodes()).release();
591  _cut_mesh->add_node(this_node);
592 
593  dof_id_type id = _cut_mesh->n_nodes() - 1;
594  auto it = new_boundary_order.begin();
595  new_boundary_order.insert(it + i, id);
596  }
597  }
598  }
599 
600  _boundary = new_boundary_order;
601  mooseAssert(_boundary.size() > 0, "Boundary should not have zero size");
602  _boundary.pop_back();
603 }
Real _size_control
Used for cutter mesh refinement and front advancement.
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
std::vector< dof_id_type > _boundary
Boundary nodes of the cutter mesh.
Real findDistance(dof_id_type node1, dof_id_type node2)
Find distance between two nodes.

◆ refineFront()

void MeshCut3DUserObject::refineFront ( )
protected

Refine the mesh at the front.

Definition at line 885 of file MeshCut3DUserObject.C.

Referenced by initialize().

886 {
887  std::vector<std::vector<dof_id_type>> new_front(_front.begin(), _front.end());
888 
889  for (unsigned int ifront = 0; ifront < _front.size(); ++ifront)
890  {
891  unsigned int i1 = _front[ifront].size() - 1;
892  if (_inactive_boundary_pos.size() == 0)
893  i1 = _front[ifront].size();
894 
895  for (unsigned int i = i1; i >= 1; --i)
896  {
897  unsigned int i2 = i;
898  if (_inactive_boundary_pos.size() == 0)
899  i2 = (i <= _front[ifront].size() - 1 ? i : 0);
900 
901  dof_id_type node1 = _front[ifront][i - 1];
902  dof_id_type node2 = _front[ifront][i2];
903  Real distance = findDistance(node1, node2);
904 
905  if (distance > _size_control)
906  {
907  unsigned int n = static_cast<int>(distance / _size_control);
908  std::array<Real, LIBMESH_DIM> x1;
909  std::array<Real, LIBMESH_DIM> x2;
910 
911  Node * this_node = _cut_mesh->node_ptr(node1);
912  mooseAssert(this_node, "Node is NULL");
913  Point & p1 = *this_node;
914  this_node = _cut_mesh->node_ptr(node2);
915  mooseAssert(this_node, "Node is NULL");
916  Point & p2 = *this_node;
917 
918  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
919  {
920  x1[j] = p1(j);
921  x2[j] = p2(j);
922  }
923 
924  for (unsigned int j = 0; j < n; ++j)
925  {
926  Point x;
927  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
928  x(k) = x2[k] - (x2[k] - x1[k]) * (j + 1) / (n + 1);
929 
930  Node * this_node = Node::build(x, _cut_mesh->n_nodes()).release();
931  _cut_mesh->add_node(this_node);
932 
933  dof_id_type id = _cut_mesh->n_nodes() - 1;
934 
935  auto it = new_front[ifront].begin();
936  new_front[ifront].insert(it + i, id);
937  }
938  }
939  }
940  }
941 
942  _front = new_front;
943 }
Real _size_control
Used for cutter mesh refinement and front advancement.
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
std::vector< std::vector< dof_id_type > > _front
New boundary after growth.
Real findDistance(dof_id_type node1, dof_id_type node2)
Find distance between two nodes.
std::vector< unsigned int > _inactive_boundary_pos
Inactive boundary.

◆ serialize()

void GeometricCutUserObject::serialize ( std::string &  serialized_buffer)
protectedinherited

Methods to pack/unpack the _marked_elems_2d and _marked_elems_3d data into a structure suitable for parallel communication.

Definition at line 209 of file GeometricCutUserObject.C.

Referenced by GeometricCutUserObject::finalize().

210 {
211  // stream for serializing the _marked_elems_2d and _marked_elems_3d data structures to a byte
212  // stream
213  std::ostringstream oss;
214  dataStore(oss, _marked_elems_2d, this);
215  dataStore(oss, _marked_elems_3d, this);
216 
217  // Populate the passed in string pointer with the string stream's buffer contents
218  serialized_buffer.assign(oss.str());
219 }
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo2D > > _marked_elems_2d
Containers with information about all 2D and 3D elements marked for cutting by this object...
void dataStore(std::ostream &stream, Xfem::CutFace &cf, void *context)
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo3D > > _marked_elems_3d

◆ setInterfaceID()

void GeometricCutUserObject::setInterfaceID ( unsigned int  interface_id)
inlineinherited

Set the interface ID for this cutting object.

Parameters
theinterface ID

Definition at line 179 of file GeometricCutUserObject.h.

Referenced by XFEM::addGeometricCut().

179 { _interface_id = interface_id; };
unsigned int _interface_id
Associated interface id.

◆ shouldHealMesh()

bool GeometricCutUserObject::shouldHealMesh ( ) const
inlineinherited

Should the elements cut by this cutting object be healed in the current time step?

Returns
true if the cut element should be healed

Definition at line 186 of file GeometricCutUserObject.h.

186 { return _heal_always; };
bool _heal_always
Heal the mesh.

◆ sortBoundaryNodes()

void MeshCut3DUserObject::sortBoundaryNodes ( )
protected

Sort boundary nodes to be in the right order along the boundary.

Definition at line 481 of file MeshCut3DUserObject.C.

Referenced by initialize().

482 {
483  _boundary.clear();
484 
485  for (auto it = _boundary_edges.begin(); it != _boundary_edges.end(); ++it)
486  {
487  dof_id_type node1 = (*it)._id1;
488  dof_id_type node2 = (*it)._id2;
489 
490  mooseAssert(_boundary_map.find(node1) != _boundary_map.end(),
491  "_boundary_map does not have this key");
492  mooseAssert(_boundary_map.find(node2) != _boundary_map.end(),
493  "_boundary_map does not have this key");
494 
495  _boundary_map.find(node1)->second.push_back(node2);
496  _boundary_map.find(node2)->second.push_back(node1);
497  }
498 
499  auto it = _boundary_map.begin();
500  while (it != _boundary_map.end())
501  {
502  if (it->second.size() != 2)
503  mooseError(
504  "Boundary nodes in the cutter mesh must have exactly two neighbors; this one has: ",
505  it->second.size());
506  ++it;
507  }
508 
509  auto it2 = _boundary_edges.begin();
510  dof_id_type node1 = (*it2)._id1;
511  dof_id_type node2 = (*it2)._id2;
512  _boundary.push_back(node1);
513  _boundary.push_back(node2);
514 
515  for (unsigned int i = 0; i < _boundary_edges.size() - 1; ++i)
516  {
517  mooseAssert(_boundary_map.find(node2) != _boundary_map.end(),
518  "_boundary_map does not have this key");
519 
520  dof_id_type node3 = _boundary_map.find(node2)->second[0];
521  dof_id_type node4 = _boundary_map.find(node2)->second[1];
522 
523  if (node3 == node1)
524  {
525  _boundary.push_back(node4);
526  node1 = node2;
527  node2 = node4;
528  }
529  else if (node4 == node1)
530  {
531  _boundary.push_back(node3);
532  node1 = node2;
533  node2 = node3;
534  }
535  else
536  mooseError("Discontinuity in cutter boundary");
537  }
538 }
std::map< dof_id_type, std::vector< dof_id_type > > _boundary_map
A map of boundary nodes and their neighbors.
std::set< Xfem::CutEdge > _boundary_edges
Edges at the boundary.
std::vector< dof_id_type > _boundary
Boundary nodes of the cutter mesh.

◆ sortFrontNodes()

void MeshCut3DUserObject::sortFrontNodes ( )
protected

Sort the front nodes.

Definition at line 772 of file MeshCut3DUserObject.C.

Referenced by initialize().

774 {
775 }

◆ threadJoin()

void GeometricCutUserObject::threadJoin ( const UserObject &  y)
overridevirtualinherited

Definition at line 133 of file GeometricCutUserObject.C.

134 {
135  const GeometricCutUserObject & gcuo = dynamic_cast<const GeometricCutUserObject &>(y);
136 
137  for (const auto & it : gcuo._marked_elems_2d)
138  {
139  mooseAssert(_marked_elems_2d.find(it.first) == _marked_elems_2d.end(),
140  "Element already inserted in map from a different thread");
141  _marked_elems_2d[it.first] = it.second;
142  }
143  for (const auto & it : gcuo._marked_elems_3d)
144  {
145  mooseAssert(_marked_elems_3d.find(it.first) == _marked_elems_3d.end(),
146  "Element already inserted in map from a different thread");
147  _marked_elems_3d[it.first] = it.second;
148  }
149 }
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo2D > > _marked_elems_2d
Containers with information about all 2D and 3D elements marked for cutting by this object...
std::map< unsigned int, std::vector< Xfem::GeomMarkedElemInfo3D > > _marked_elems_3d

◆ triangulation()

void MeshCut3DUserObject::triangulation ( )
protected

Create tri3 elements between the new front and the old front.

Definition at line 946 of file MeshCut3DUserObject.C.

Referenced by initialize().

947 {
948 
949  mooseAssert(_active_boundary.size() == _front.size(),
950  "_active_boundary and _front should have the same size!");
951 
952  if (_inactive_boundary_pos.size() == 0)
953  {
954  _active_boundary[0].push_back(_active_boundary[0][0]);
955  _front[0].push_back(_front[0][0]);
956  }
957 
958  // loop over active segments
959  for (unsigned int k = 0; k < _front.size(); ++k)
960  {
961  unsigned int n1 = _active_boundary[k].size();
962  unsigned int n2 = _front[k].size();
963 
964  unsigned int i1 = 0;
965  unsigned int i2 = 0;
966 
967  // stop when all nodes are associated with an element
968  while (!(i1 == n1 - 1 && i2 == n2 - 1))
969  {
970  std::vector<dof_id_type> elem;
971 
972  dof_id_type p1 = _active_boundary[k][i1]; // node in the old front
973  dof_id_type p2 = _front[k][i2]; // node in the new front
974 
975  if (i1 != n1 - 1 && i2 != n2 - 1)
976  {
977  dof_id_type p3 = _active_boundary[k][i1 + 1]; // next node in the old front
978  dof_id_type p4 = _front[k][i2 + 1]; // next node in the new front
979 
980  elem.push_back(p1);
981  elem.push_back(p2);
982 
983  Real d1 = findDistance(p1, p4);
984  Real d2 = findDistance(p3, p2);
985 
986  if (d1 < d2)
987  {
988  elem.push_back(p4);
989  i2++;
990  }
991 
992  else
993  {
994  elem.push_back(p3);
995  i1++;
996  }
997  }
998 
999  else if (i1 == n1 - 1)
1000  {
1001  dof_id_type p4 = _front[k][i2 + 1]; // next node in the new front
1002 
1003  elem.push_back(p1);
1004  elem.push_back(p2);
1005  elem.push_back(p4);
1006  i2++;
1007  }
1008 
1009  else if (i2 == n2 - 1)
1010  {
1011  dof_id_type p3 = _active_boundary[k][i1 + 1]; // next node in the old front
1012 
1013  elem.push_back(p1);
1014  elem.push_back(p2);
1015  elem.push_back(p3);
1016  i1++;
1017  }
1018 
1019  Elem * new_elem = Elem::build(TRI3).release();
1020 
1021  for (unsigned int i = 0; i < _cut_elem_nnode; ++i)
1022  {
1023  mooseAssert(_cut_mesh->node_ptr(elem[i]) != nullptr, "Node is NULL");
1024  new_elem->set_node(i) = _cut_mesh->node_ptr(elem[i]);
1025  }
1026 
1027  _cut_mesh->add_elem(new_elem);
1028  }
1029  }
1030 }
const unsigned int _cut_elem_nnode
The cutter mesh has triangluar elements only.
std::vector< std::vector< dof_id_type > > _active_boundary
Active boundary nodes where growth is allowed.
std::unique_ptr< MeshBase > _cut_mesh
The cutter mesh.
std::vector< std::vector< dof_id_type > > _front
New boundary after growth.
Real findDistance(dof_id_type node1, dof_id_type node2)
Find distance between two nodes.
std::vector< unsigned int > _inactive_boundary_pos
Inactive boundary.

Member Data Documentation

◆ _active_boundary

std::vector<std::vector<dof_id_type> > MeshCut3DUserObject::_active_boundary
protected

Active boundary nodes where growth is allowed.

Definition at line 89 of file MeshCut3DUserObject.h.

Referenced by findActiveBoundaryDirection(), findActiveBoundaryNodes(), growFront(), joinBoundary(), and triangulation().

◆ _active_direction

std::vector<std::vector<Point> > MeshCut3DUserObject::_active_direction
protected

Growth direction for active boundaries.

Definition at line 92 of file MeshCut3DUserObject.h.

Referenced by findActiveBoundaryDirection(), and growFront().

◆ _boundary

std::vector<dof_id_type> MeshCut3DUserObject::_boundary
protected

Boundary nodes of the cutter mesh.

Definition at line 80 of file MeshCut3DUserObject.h.

Referenced by findActiveBoundaryNodes(), joinBoundary(), refineBoundary(), and sortBoundaryNodes().

◆ _boundary_edges

std::set<Xfem::CutEdge> MeshCut3DUserObject::_boundary_edges
protected

Edges at the boundary.

Definition at line 83 of file MeshCut3DUserObject.h.

Referenced by findBoundaryEdges(), and sortBoundaryNodes().

◆ _boundary_map

std::map<dof_id_type, std::vector<dof_id_type> > MeshCut3DUserObject::_boundary_map
protected

A map of boundary nodes and their neighbors.

Definition at line 86 of file MeshCut3DUserObject.h.

Referenced by findBoundaryEdges(), findBoundaryNodes(), and sortBoundaryNodes().

◆ _const_intersection

const Real MeshCut3DUserObject::_const_intersection = 0.01
protected

Used to define intersection points.

Definition at line 67 of file MeshCut3DUserObject.h.

Referenced by findFrontIntersection().

◆ _cut_elem_dim

const unsigned int MeshCut3DUserObject::_cut_elem_dim = 2
protected

Definition at line 58 of file MeshCut3DUserObject.h.

Referenced by MeshCut3DUserObject().

◆ _cut_elem_nnode

const unsigned int MeshCut3DUserObject::_cut_elem_nnode = 3
protected

The cutter mesh has triangluar elements only.

Definition at line 57 of file MeshCut3DUserObject.h.

Referenced by findBoundaryEdges(), findBoundaryNodes(), MeshCut3DUserObject(), and triangulation().

◆ _cut_mesh

std::unique_ptr<MeshBase> MeshCut3DUserObject::_cut_mesh
protected

◆ _elem_dim

const unsigned int MeshCut3DUserObject::_elem_dim = 3
protected

The structural mesh must be 3D only.

Definition at line 64 of file MeshCut3DUserObject.h.

◆ _front

std::vector<std::vector<dof_id_type> > MeshCut3DUserObject::_front
protected

New boundary after growth.

Definition at line 98 of file MeshCut3DUserObject.h.

Referenced by findFrontIntersection(), growFront(), joinBoundary(), refineFront(), and triangulation().

◆ _func_x

Function& MeshCut3DUserObject::_func_x
protected

Parsed functions of front growth.

Definition at line 203 of file MeshCut3DUserObject.h.

Referenced by findActiveBoundaryDirection().

◆ _func_y

Function& MeshCut3DUserObject::_func_y
protected

Definition at line 204 of file MeshCut3DUserObject.h.

Referenced by findActiveBoundaryDirection().

◆ _func_z

Function& MeshCut3DUserObject::_func_z
protected

Definition at line 205 of file MeshCut3DUserObject.h.

Referenced by findActiveBoundaryDirection().

◆ _grow

bool MeshCut3DUserObject::_grow
protected

Definition at line 77 of file MeshCut3DUserObject.h.

Referenced by initialize(), and MeshCut3DUserObject().

◆ _heal_always

bool GeometricCutUserObject::_heal_always
protectedinherited

Heal the mesh.

Definition at line 196 of file GeometricCutUserObject.h.

◆ _inactive_boundary_pos

std::vector<unsigned int> MeshCut3DUserObject::_inactive_boundary_pos
protected

◆ _interface_id

unsigned int GeometricCutUserObject::_interface_id
protectedinherited

Associated interface id.

Definition at line 193 of file GeometricCutUserObject.h.

Referenced by GeometricCutUserObject::finalize(), and GeometricCutUserObject::GeometricCutUserObject().

◆ _marked_elems_2d

std::map<unsigned int, std::vector<Xfem::GeomMarkedElemInfo2D> > GeometricCutUserObject::_marked_elems_2d
protectedinherited

◆ _marked_elems_3d

std::map<unsigned int, std::vector<Xfem::GeomMarkedElemInfo3D> > GeometricCutUserObject::_marked_elems_3d
protectedinherited

◆ _mesh

MooseMesh& MeshCut3DUserObject::_mesh
protected

The structural mesh.

Definition at line 61 of file MeshCut3DUserObject.h.

Referenced by findActiveBoundaryNodes(), and findFrontIntersection().

◆ _n_step_growth

unsigned int MeshCut3DUserObject::_n_step_growth
protected

Number of steps to grow the mesh.

Definition at line 73 of file MeshCut3DUserObject.h.

Referenced by initialize(), and MeshCut3DUserObject().

◆ _size_control

Real MeshCut3DUserObject::_size_control
protected

Used for cutter mesh refinement and front advancement.

Definition at line 70 of file MeshCut3DUserObject.h.

Referenced by growFront(), MeshCut3DUserObject(), refineBoundary(), and refineFront().

◆ _stop

bool MeshCut3DUserObject::_stop
protected

variables to help control the work flow

Definition at line 76 of file MeshCut3DUserObject.h.

Referenced by findActiveBoundaryNodes(), and initialize().

◆ _xfem

std::shared_ptr<XFEM> GeometricCutUserObject::_xfem
protectedinherited

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