Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "ElementPairLocator.h"
13 : #include "ElementFragmentAlgorithm.h"
14 : #include "XFEMInterface.h"
15 : #include "XFEMCrackGrowthIncrement2DCut.h"
16 :
17 : #include "libmesh/vector_value.h"
18 : #include "libmesh/quadrature.h"
19 :
20 : #include "GeometricCutUserObject.h"
21 :
22 : // Forward declarations
23 : class SystemBase;
24 :
25 : namespace Xfem
26 : {
27 : enum XFEM_CUTPLANE_QUANTITY
28 : {
29 : ORIGIN_X,
30 : ORIGIN_Y,
31 : ORIGIN_Z,
32 : NORMAL_X,
33 : NORMAL_Y,
34 : NORMAL_Z
35 : };
36 :
37 : enum XFEM_QRULE
38 : {
39 : VOLFRAC,
40 : MOMENT_FITTING,
41 : DIRECT
42 : };
43 :
44 : /**
45 : * Convenient typedef for local storage of stateful material properties. The first (component 0)
46 : * entry in the CachedMaterialProperties is a map for old material properties. The second
47 : * (component 1) entry is a map for older material properties. The second entry will be empty if
48 : * the material storage doesn't have older material properties.
49 : */
50 : typedef std::array<std::unordered_map<unsigned int, std::string>, 2> CachedMaterialProperties;
51 :
52 : /**
53 : * Information about a cut element. This is a tuple of (0) the parent
54 : * element, (1) the geometric cut userobject that cuts the element, (2) the cut
55 : * subdomain ID, and (3) the stateful material properties.
56 : */
57 13909 : struct CutElemInfo
58 : {
59 : const Elem * _parent_elem;
60 : const GeometricCutUserObject * _geometric_cut;
61 : CutSubdomainID _cut_subdomain_id;
62 : CachedMaterialProperties _elem_material_properties;
63 : CachedMaterialProperties _bnd_material_properties;
64 : // TODO: add neighbor material properties
65 : // CachedMaterialProperties _neighbor_material_properties;
66 :
67 0 : CutElemInfo()
68 0 : : _parent_elem(nullptr),
69 0 : _geometric_cut(nullptr),
70 0 : _cut_subdomain_id(std::numeric_limits<CutSubdomainID>::max())
71 : {
72 0 : }
73 :
74 : CutElemInfo(const Elem * parent_elem,
75 : const GeometricCutUserObject * geometric_cut,
76 : CutSubdomainID cut_subdomain_id)
77 1980 : : _parent_elem(parent_elem), _geometric_cut(geometric_cut), _cut_subdomain_id(cut_subdomain_id)
78 : {
79 : }
80 :
81 : // A new child element is said to be previously healed if an entry in the old CutElemInfo has the
82 : // same parent element ID, the same cut, AND the same cut subdomain ID.
83 : bool match(const CutElemInfo & rhs)
84 : {
85 9949 : return _parent_elem == rhs._parent_elem && _geometric_cut == rhs._geometric_cut &&
86 1209 : _cut_subdomain_id == rhs._cut_subdomain_id;
87 : }
88 : };
89 : } // namespace Xfem
90 :
91 : class XFEMCutElem;
92 : class XFEMCrackGrowthIncrement2DCut;
93 : class EFANode;
94 : class EFAEdge;
95 : class EFAElement;
96 : class EFAElement2D;
97 : class EFAElement3D;
98 :
99 : /**
100 : * This is the \p XFEM class. This class implements
101 : * algorithms for dynamic mesh modification in support of
102 : * a phantom node approach for XFEM
103 : */
104 :
105 : // ------------------------------------------------------------
106 : // XFEM class definition
107 : class XFEM : public XFEMInterface
108 : {
109 : public:
110 : explicit XFEM(const InputParameters & params);
111 :
112 : ~XFEM();
113 :
114 : void addGeometricCut(GeometricCutUserObject * geometric_cut);
115 :
116 : void addStateMarkedElem(unsigned int elem_id, RealVectorValue & normal);
117 : void addStateMarkedElem(unsigned int elem_id, RealVectorValue & normal, unsigned int marked_side);
118 : void addStateMarkedFrag(unsigned int elem_id, RealVectorValue & normal);
119 :
120 : void clearStateMarkedElems();
121 :
122 : /**
123 : * Add information about a new cut to be performed on a specific 2d element
124 : * @param elem_id The id of the element to be cut
125 : * @param geom_info The object containing information about the cut to be performed
126 : * @param interface_id The ID of the interface
127 : */
128 : void addGeomMarkedElem2D(const unsigned int elem_id,
129 : const Xfem::GeomMarkedElemInfo2D geom_info,
130 : const unsigned int interface_id);
131 :
132 : /**
133 : * Add information about a new cut to be performed on a specific 3d element
134 : * @param elem_id The id of the element to be cut
135 : * @param geom_info The object containing information about the cut to be performed
136 : * @param interface_id The ID of the interface
137 : */
138 : void addGeomMarkedElem3D(const unsigned int elem_id,
139 : const Xfem::GeomMarkedElemInfo3D geom_info,
140 : const unsigned int interface_id);
141 :
142 : /**
143 : * Clear out the list of elements to be marked for cutting. Called after cutting is done.
144 : */
145 : void clearGeomMarkedElems();
146 :
147 : virtual bool update(Real time,
148 : const std::vector<std::shared_ptr<NonlinearSystemBase>> & nl,
149 : AuxiliarySystem & aux) override;
150 :
151 : virtual void initSolution(const std::vector<std::shared_ptr<NonlinearSystemBase>> & nl,
152 : AuxiliarySystem & aux) override;
153 :
154 : void buildEFAMesh();
155 : bool markCuts(Real time);
156 : bool markCutEdgesByGeometry();
157 : bool markCutEdgesByState(Real time);
158 : bool markCutFacesByGeometry();
159 : bool markCutFacesByState();
160 : bool initCutIntersectionEdge(
161 : Point cut_origin, RealVectorValue cut_normal, Point & edge_p1, Point & edge_p2, Real & dist);
162 : bool cutMeshWithEFA(const std::vector<std::shared_ptr<NonlinearSystemBase>> & nl,
163 : AuxiliarySystem & aux);
164 :
165 : /**
166 : * Potentially heal the mesh by merging some of the pairs
167 : * of partial elements cut by XFEM back into single elements
168 : * if indicated by the cutting objects.
169 : * @return true if the mesh has been modified due to healing
170 : **/
171 : bool healMesh();
172 :
173 : virtual bool updateHeal() override;
174 : Point getEFANodeCoords(EFANode * CEMnode,
175 : EFAElement * CEMElem,
176 : const Elem * elem,
177 : MeshBase * displaced_mesh = nullptr) const;
178 :
179 : /**
180 : * Get the volume fraction of an element that is physical
181 : */
182 : Real getPhysicalVolumeFraction(const Elem * elem) const;
183 :
184 : /**
185 : * Return true if the point is inside the element physical domain
186 : * Note: if this element is not cut, return true too
187 : */
188 : bool isPointInsidePhysicalDomain(const Elem * elem, const Point & point) const;
189 :
190 : /**
191 : * Get specified component of normal or origin for cut plane for a given element
192 : */
193 : Real getCutPlane(const Elem * elem,
194 : const Xfem::XFEM_CUTPLANE_QUANTITY quantity,
195 : unsigned int plane_id) const;
196 :
197 : bool isElemAtCrackTip(const Elem * elem) const;
198 : bool isElemCut(const Elem * elem, XFEMCutElem *& xfce) const;
199 : bool isElemCut(const Elem * elem) const;
200 : void getFragmentFaces(const Elem * elem,
201 : std::vector<std::vector<Point>> & frag_faces,
202 : bool displaced_mesh = false) const;
203 : void storeCrackTipOriginAndDirection();
204 :
205 : void correctCrackExtensionDirection(const Elem * elem,
206 : EFAElement2D * CEMElem,
207 : EFAEdge * orig_edge,
208 : Point normal,
209 : Point crack_tip_origin,
210 : Point crack_tip_direction,
211 : Real & distance_keep,
212 : unsigned int & edge_id_keep,
213 : Point & normal_keep);
214 :
215 : void getCrackTipOrigin(std::map<unsigned int, const Elem *> & elem_id_crack_tip,
216 : std::vector<Point> & crack_front_points);
217 :
218 : Xfem::XFEM_QRULE & getXFEMQRule();
219 : void setXFEMQRule(std::string & xfem_qrule);
220 : void setCrackGrowthMethod(bool use_crack_growth_increment, Real crack_growth_increment);
221 :
222 : /**
223 : * Controls amount of debugging information output
224 : * @param debug_output_level How much information to output (see description of
225 : * _debug_output_level)
226 : */
227 : void setDebugOutputLevel(unsigned int debug_output_level);
228 :
229 : virtual bool getXFEMWeights(MooseArray<Real> & weights,
230 : const Elem * elem,
231 : QBase * qrule,
232 : const MooseArray<Point> & q_points) override;
233 : virtual bool getXFEMFaceWeights(MooseArray<Real> & weights,
234 : const Elem * elem,
235 : QBase * qrule,
236 : const MooseArray<Point> & q_points,
237 : unsigned int side) override;
238 :
239 : /**
240 : * Get the list of cut element pairs corresponding to a given
241 : * interface ID.
242 : * @param interface_id The ID of the interface
243 : * @return the list of elements cut by that interface
244 : **/
245 541 : virtual const ElementPairLocator::ElementPairList * getXFEMCutElemPairs(unsigned int interface_id)
246 : {
247 541 : return &_sibling_elems[interface_id];
248 : }
249 :
250 : /**
251 : * Get the list of cut element pairs on the displaced mesh
252 : * corresponding to a given interface ID.
253 : * @param interface_id The ID of the interface
254 : * @return the list of elements cut by that interface
255 : **/
256 : virtual const ElementPairLocator::ElementPairList *
257 225 : getXFEMDisplacedCutElemPairs(unsigned int interface_id)
258 : {
259 225 : return &_sibling_displaced_elems[interface_id];
260 : }
261 :
262 : /**
263 : * Get the interface ID corresponding to a given GeometricCutUserObject.
264 : * @param gcu pointer to the GeometricCutUserObject
265 : * @return the interface ID
266 : **/
267 268 : virtual unsigned int getGeometricCutID(const GeometricCutUserObject * gcu)
268 : {
269 268 : return _geom_marker_id_map[gcu];
270 : };
271 :
272 : virtual void getXFEMIntersectionInfo(const Elem * elem,
273 : unsigned int plane_id,
274 : Point & normal,
275 : std::vector<Point> & intersectionPoints,
276 : bool displaced_mesh = false) const;
277 : virtual void getXFEMqRuleOnLine(std::vector<Point> & intersection_points,
278 : std::vector<Point> & quad_pts,
279 : std::vector<Real> & quad_wts) const;
280 : virtual void getXFEMqRuleOnSurface(std::vector<Point> & intersection_points,
281 : std::vector<Point> & quad_pts,
282 : std::vector<Real> & quad_wts) const;
283 36812 : bool has_secondary_cut() { return _has_secondary_cut; }
284 :
285 : /**
286 : * Get the EFAElement2D object for a specified libMesh element
287 : * @param elem Pointer to the libMesh element for which the object is requested
288 : */
289 : EFAElement2D * getEFAElem2D(const Elem * elem);
290 :
291 : /**
292 : * Get the EFAElement3D object for a specified libMesh element
293 : * @param elem Pointer to the libMesh element for which the object is requested
294 : */
295 : EFAElement3D * getEFAElem3D(const Elem * elem);
296 :
297 : void getFragmentEdges(const Elem * elem,
298 : EFAElement2D * CEMElem,
299 : std::vector<std::vector<Point>> & frag_edges) const;
300 : void getFragmentFaces(const Elem * elem,
301 : EFAElement3D * CEMElem,
302 : std::vector<std::vector<Point>> & frag_faces) const;
303 :
304 : const std::map<const Elem *, std::vector<Point>> & getCrackTipOriginMap() const
305 : {
306 : return _elem_crack_origin_direction_map;
307 : }
308 :
309 : /**
310 : * Determine which cut subdomain the element belongs to relative to the cut
311 : * @param gcuo The GeometricCutUserObject for the cut
312 : * @param cut_elem The element being cut
313 : * @param parent_elem The parent element
314 : */
315 : CutSubdomainID getCutSubdomainID(const GeometricCutUserObject * gcuo,
316 : const Elem * cut_elem,
317 : const Elem * parent_elem = nullptr) const;
318 :
319 : private:
320 : bool _has_secondary_cut;
321 :
322 : Xfem::XFEM_QRULE _XFEM_qrule;
323 :
324 : bool _use_crack_growth_increment;
325 : Real _crack_growth_increment;
326 :
327 : std::vector<const GeometricCutUserObject *> _geometric_cuts;
328 :
329 : std::map<unique_id_type, XFEMCutElem *> _cut_elem_map;
330 : std::set<const Elem *> _crack_tip_elems;
331 : std::set<const Elem *> _crack_tip_elems_to_be_healed;
332 : std::map<unsigned int, ElementPairLocator::ElementPairList> _sibling_elems;
333 : std::map<unsigned int, ElementPairLocator::ElementPairList> _sibling_displaced_elems;
334 :
335 : std::map<const Elem *, std::vector<Point>> _elem_crack_origin_direction_map;
336 :
337 : // std::map<const Elem*, Point> _crack_propagation_direction_map;
338 :
339 : std::map<const Elem *, RealVectorValue> _state_marked_elems;
340 : std::set<const Elem *> _state_marked_frags;
341 : std::map<const Elem *, unsigned int> _state_marked_elem_sides;
342 :
343 : /// Data structure for storing information about all 2D elements to be cut by geometry
344 : std::map<const Elem *, std::vector<Xfem::GeomMarkedElemInfo2D>> _geom_marked_elems_2d;
345 :
346 : /// Data structure for storing information about all 3D elements to be cut by geometry
347 : std::map<const Elem *, std::vector<Xfem::GeomMarkedElemInfo3D>> _geom_marked_elems_3d;
348 :
349 : /// Data structure for storing the elements cut by specific geometric cutters
350 : std::map<unsigned int, std::set<unsigned int>> _geom_marker_id_elems;
351 :
352 : /// Data structure for storing the GeommetricCutUserObjects and their corresponding id
353 : std::map<const GeometricCutUserObject *, unsigned int> _geom_marker_id_map;
354 :
355 : ElementFragmentAlgorithm _efa_mesh;
356 :
357 : /// Controls amount of debugging output information
358 : /// 0: None
359 : /// 1: Summary
360 : /// 2: Details on modifications to mesh
361 : /// 3: Full dump of element fragment algorithm mesh
362 : unsigned int _debug_output_level;
363 :
364 : /**
365 : * Data structure to store the nonlinear solution for nodes/elements affected by XFEM
366 : * For each node/element, this is stored as a vector that contains all components
367 : * of all applicable variables in an order defined by getElementSolutionDofs() or
368 : * getNodeSolutionDofs(). This vector first contains the current solution in that
369 : * order, followed by the old and older solutions, also in that same order.
370 : */
371 : std::map<unique_id_type, std::vector<Real>> _cached_solution;
372 :
373 : /**
374 : * Data structure to store the auxiliary solution for nodes/elements affected by XFEM
375 : * For each node/element, this is stored as a vector that contains all components
376 : * of all applicable variables in an order defined by getElementSolutionDofs() or
377 : * getNodeSolutionDofs(). This vector first contains the current solution in that
378 : * order, followed by the old and older solutions, also in that same order.
379 : */
380 : std::map<unique_id_type, std::vector<Real>> _cached_aux_solution;
381 :
382 : /**
383 : * All geometrically cut elements and their CutElemInfo during the current execution of
384 : * XFEM_MARK. This data structure is updated everytime a new cut element is created.
385 : */
386 : std::unordered_map<const Elem *, Xfem::CutElemInfo> _geom_cut_elems;
387 :
388 : /**
389 : * All geometrically cut elements and their CutElemInfo before the current execution of
390 : * XFEM_MARK.
391 : */
392 : std::unordered_map<const Elem *, Xfem::CutElemInfo> _old_geom_cut_elems;
393 :
394 : /**
395 : * Store the solution in stored_solution for a given node
396 : * @param node_to_store_to Node for which the solution will be stored
397 : * @param node_to_store_from Node from which the solution to be stored is obtained
398 : * @param sys System from which the solution is stored
399 : * @param stored_solution Data structure that the stored solution is saved to
400 : * @param current_solution Current solution vector that the solution is obtained from
401 : * @param old_solution Old solution vector that the solution is obtained from
402 : * @param older_solution Older solution vector that the solution is obtained from
403 : */
404 : void storeSolutionForNode(const Node * node_to_store_to,
405 : const Node * node_to_store_from,
406 : SystemBase & sys,
407 : std::map<unique_id_type, std::vector<Real>> & stored_solution,
408 : const NumericVector<Number> & current_solution,
409 : const NumericVector<Number> & old_solution,
410 : const NumericVector<Number> & older_solution);
411 :
412 : /**
413 : * Store the solution in stored_solution for a given element
414 : * @param elem_to_store_to Element for which the solution will be stored
415 : * @param elem_to_store_from Element from which the solution to be stored is obtained
416 : * @param sys System from which the solution is stored
417 : * @param stored_solution Data structure that the stored solution is saved to
418 : * @param current_solution Current solution vector that the solution is obtained from
419 : * @param old_solution Old solution vector that the solution is obtained from
420 : * @param older_solution Older solution vector that the solution is obtained from
421 : */
422 : void storeSolutionForElement(const Elem * elem_to_store_to,
423 : const Elem * elem_to_store_from,
424 : SystemBase & sys,
425 : std::map<unique_id_type, std::vector<Real>> & stored_solution,
426 : const NumericVector<Number> & current_solution,
427 : const NumericVector<Number> & old_solution,
428 : const NumericVector<Number> & older_solution);
429 :
430 : /**
431 : * Set the solution for all locally-owned nodes/elements that have stored values
432 : * @param sys System for which the solution is set
433 : * @param stored_solution Data structure that the stored solution is obtained from
434 : * @param current_solution Current solution vector that will be set
435 : * @param old_solution Old solution vector that will be set
436 : * @param older_solution Older solution vector that will be set
437 : */
438 : void setSolution(SystemBase & sys,
439 : const std::map<unique_id_type, std::vector<Real>> & stored_solution,
440 : NumericVector<Number> & current_solution,
441 : NumericVector<Number> & old_solution,
442 : NumericVector<Number> & older_solution);
443 :
444 : /**
445 : * Set the solution for a set of DOFs
446 : * @param stored_solution Stored solution values to set the solution to
447 : * @param stored_solution_dofs Dof indices for the entries in stored_solution
448 : * @param current_solution Current solution vector that will be set
449 : * @param old_solution Old solution vector that will be set
450 : * @param older_solution Older solution vector that will be set
451 : */
452 : void setSolutionForDOFs(const std::vector<Real> & stored_solution,
453 : const std::vector<dof_id_type> & stored_solution_dofs,
454 : NumericVector<Number> & current_solution,
455 : NumericVector<Number> & old_solution,
456 : NumericVector<Number> & older_solution);
457 :
458 : /**
459 : * Get a vector of the dof indices for all components of all variables
460 : * associated with an element
461 : * @param elem Element for which dof indices are found
462 : * @param sys System for which the dof indices are found
463 : */
464 : std::vector<dof_id_type> getElementSolutionDofs(const Elem * elem, SystemBase & sys) const;
465 :
466 : /**
467 : * Get a vector of the dof indices for all components of all variables
468 : * associated with a node
469 : * @param node Node for which dof indices are found
470 : * @param sys System for which the dof indices are found
471 : */
472 : std::vector<dof_id_type> getNodeSolutionDofs(const Node * node, SystemBase & sys) const;
473 :
474 : /**
475 : * Get the GeometricCutUserObject associated with an element
476 : * @param elem The element
477 : * @return A constant pointer to the GeometricCutUserObject, nullptr if nothing found
478 : */
479 : const GeometricCutUserObject * getGeometricCutForElem(const Elem * elem) const;
480 :
481 : void storeMaterialPropertiesForElementHelper(const Elem * elem,
482 : MaterialPropertyStorage & storage);
483 :
484 : /**
485 : * Helper function to store the material properties of a healed element
486 : * @param parent_elem The parent element
487 : * @param elem1 The first child element
488 : * @param elem2 The second child element
489 : */
490 : void storeMaterialPropertiesForElement(const Elem * parent_elem, const Elem * child_elem);
491 :
492 : /**
493 : * Load the material properties
494 : * @param props_deserialized The material properties
495 : * @param props_serialized The serialized material properties
496 : *
497 : * This does very dirty things and writes back to MOOSE's stateful properties. It should
498 : * not do this in the future.
499 : */
500 : void loadMaterialPropertiesForElementHelper(const Elem * elem,
501 : const Xfem::CachedMaterialProperties & cached_props,
502 : MaterialPropertyStorage & storage) const;
503 :
504 : /**
505 : * Helper function to store the material properties of a healed element
506 : * @param elem The cut element to restore material properties to.
507 : * @param elem_from The element to copy material properties from.
508 : * @param cached_cei The material properties cache to use.
509 : */
510 : void loadMaterialPropertiesForElement(
511 : const Elem * elem,
512 : const Elem * elem_from,
513 : std::unordered_map<const Elem *, Xfem::CutElemInfo> & cached_cei) const;
514 :
515 : /**
516 : * Return the first node in the provided element that is found to be in the physical domain
517 : * @param e Constant pointer to the child element
518 : * @param e0 Constant pointer to the parent element whose nodes will be querried
519 : * @return A constant pointer to the node
520 : */
521 : const Node * pickFirstPhysicalNode(const Elem * e, const Elem * e0) const;
522 : };
|
