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 : /**
230 : * Controls the minimum average weight multiplier for each element
231 : * @param min_weight_multiplier The minimum average weight multiplier applied
232 : * by XFEM to the standard quadrature weights
233 : */
234 : void setMinWeightMultiplier(Real min_weight_multiplier);
235 :
236 : virtual bool getXFEMWeights(MooseArray<Real> & weights,
237 : const Elem * elem,
238 : QBase * qrule,
239 : const MooseArray<Point> & q_points) override;
240 : virtual bool getXFEMFaceWeights(MooseArray<Real> & weights,
241 : const Elem * elem,
242 : QBase * qrule,
243 : const MooseArray<Point> & q_points,
244 : unsigned int side) override;
245 :
246 : /**
247 : * Get the list of cut element pairs corresponding to a given
248 : * interface ID.
249 : * @param interface_id The ID of the interface
250 : * @return the list of elements cut by that interface
251 : **/
252 561 : virtual const ElementPairLocator::ElementPairList * getXFEMCutElemPairs(unsigned int interface_id)
253 : {
254 561 : return &_sibling_elems[interface_id];
255 : }
256 :
257 : /**
258 : * Get the list of cut element pairs on the displaced mesh
259 : * corresponding to a given interface ID.
260 : * @param interface_id The ID of the interface
261 : * @return the list of elements cut by that interface
262 : **/
263 : virtual const ElementPairLocator::ElementPairList *
264 245 : getXFEMDisplacedCutElemPairs(unsigned int interface_id)
265 : {
266 245 : return &_sibling_displaced_elems[interface_id];
267 : }
268 :
269 : /**
270 : * Get the interface ID corresponding to a given GeometricCutUserObject.
271 : * @param gcu pointer to the GeometricCutUserObject
272 : * @return the interface ID
273 : **/
274 268 : virtual unsigned int getGeometricCutID(const GeometricCutUserObject * gcu)
275 : {
276 268 : return _geom_marker_id_map[gcu];
277 : };
278 :
279 : virtual void getXFEMIntersectionInfo(const Elem * elem,
280 : unsigned int plane_id,
281 : Point & normal,
282 : std::vector<Point> & intersectionPoints,
283 : bool displaced_mesh = false) const;
284 : virtual void getXFEMqRuleOnLine(std::vector<Point> & intersection_points,
285 : std::vector<Point> & quad_pts,
286 : std::vector<Real> & quad_wts) const;
287 : virtual void getXFEMqRuleOnSurface(std::vector<Point> & intersection_points,
288 : std::vector<Point> & quad_pts,
289 : std::vector<Real> & quad_wts) const;
290 39468 : bool has_secondary_cut() { return _has_secondary_cut; }
291 :
292 : /**
293 : * Get the EFAElement2D object for a specified libMesh element
294 : * @param elem Pointer to the libMesh element for which the object is requested
295 : */
296 : EFAElement2D * getEFAElem2D(const Elem * elem);
297 :
298 : /**
299 : * Get the EFAElement3D object for a specified libMesh element
300 : * @param elem Pointer to the libMesh element for which the object is requested
301 : */
302 : EFAElement3D * getEFAElem3D(const Elem * elem);
303 :
304 : void getFragmentEdges(const Elem * elem,
305 : EFAElement2D * CEMElem,
306 : std::vector<std::vector<Point>> & frag_edges) const;
307 : void getFragmentFaces(const Elem * elem,
308 : EFAElement3D * CEMElem,
309 : std::vector<std::vector<Point>> & frag_faces) const;
310 :
311 : const std::map<const Elem *, std::vector<Point>> & getCrackTipOriginMap() const
312 : {
313 : return _elem_crack_origin_direction_map;
314 : }
315 :
316 : /**
317 : * Determine which cut subdomain the element belongs to relative to the cut
318 : * @param gcuo The GeometricCutUserObject for the cut
319 : * @param cut_elem The element being cut
320 : * @param parent_elem The parent element
321 : */
322 : CutSubdomainID getCutSubdomainID(const GeometricCutUserObject * gcuo,
323 : const Elem * cut_elem,
324 : const Elem * parent_elem = nullptr) const;
325 :
326 : /**
327 : * Check whether the near-tip enrichment changed during the most recent call to
328 : * XFEM::update().
329 : * @return true if it changed
330 : */
331 : virtual bool didNearTipEnrichmentChange() override;
332 :
333 : private:
334 : bool _has_secondary_cut;
335 :
336 : Xfem::XFEM_QRULE _XFEM_qrule;
337 :
338 : bool _use_crack_growth_increment;
339 : Real _crack_growth_increment;
340 :
341 : std::vector<const GeometricCutUserObject *> _geometric_cuts;
342 :
343 : std::map<unique_id_type, XFEMCutElem *> _cut_elem_map;
344 : std::set<const Elem *> _crack_tip_elems;
345 : std::set<const Elem *> _crack_tip_elems_to_be_healed;
346 : std::map<unsigned int, ElementPairLocator::ElementPairList> _sibling_elems;
347 : std::map<unsigned int, ElementPairLocator::ElementPairList> _sibling_displaced_elems;
348 :
349 : std::map<const Elem *, std::vector<Point>> _elem_crack_origin_direction_map;
350 :
351 : // std::map<const Elem*, Point> _crack_propagation_direction_map;
352 :
353 : std::map<const Elem *, RealVectorValue> _state_marked_elems;
354 : std::set<const Elem *> _state_marked_frags;
355 : std::map<const Elem *, unsigned int> _state_marked_elem_sides;
356 :
357 : /// Data structure for storing information about all 2D elements to be cut by geometry
358 : std::map<const Elem *, std::vector<Xfem::GeomMarkedElemInfo2D>> _geom_marked_elems_2d;
359 :
360 : /// Data structure for storing information about all 3D elements to be cut by geometry
361 : std::map<const Elem *, std::vector<Xfem::GeomMarkedElemInfo3D>> _geom_marked_elems_3d;
362 :
363 : /// Data structure for storing the elements cut by specific geometric cutters
364 : std::map<unsigned int, std::set<unsigned int>> _geom_marker_id_elems;
365 :
366 : /// Data structure for storing the GeommetricCutUserObjects and their corresponding id
367 : std::map<const GeometricCutUserObject *, unsigned int> _geom_marker_id_map;
368 :
369 : ElementFragmentAlgorithm _efa_mesh;
370 :
371 : /// Controls amount of debugging output information
372 : /// 0: None
373 : /// 1: Summary
374 : /// 2: Details on modifications to mesh
375 : /// 3: Full dump of element fragment algorithm mesh
376 : unsigned int _debug_output_level;
377 :
378 : /// The minimum average multiplier applied by XFEM to the standard quadrature weights
379 : /// to integrate partial elements
380 : Real _min_weight_multiplier;
381 :
382 : /**
383 : * Data structure to store the nonlinear solution for nodes/elements affected by XFEM
384 : * For each node/element, this is stored as a vector that contains all components
385 : * of all applicable variables in an order defined by getElementSolutionDofs() or
386 : * getNodeSolutionDofs(). This vector first contains the current solution in that
387 : * order, followed by the old and older solutions, also in that same order.
388 : */
389 : std::map<unique_id_type, std::vector<Real>> _cached_solution;
390 :
391 : /**
392 : * Data structure to store the auxiliary solution for nodes/elements affected by XFEM
393 : * For each node/element, this is stored as a vector that contains all components
394 : * of all applicable variables in an order defined by getElementSolutionDofs() or
395 : * getNodeSolutionDofs(). This vector first contains the current solution in that
396 : * order, followed by the old and older solutions, also in that same order.
397 : */
398 : std::map<unique_id_type, std::vector<Real>> _cached_aux_solution;
399 :
400 : /**
401 : * All geometrically cut elements and their CutElemInfo during the current execution of
402 : * XFEM_MARK. This data structure is updated everytime a new cut element is created.
403 : */
404 : std::unordered_map<const Elem *, Xfem::CutElemInfo> _geom_cut_elems;
405 :
406 : /**
407 : * All geometrically cut elements and their CutElemInfo before the current execution of
408 : * XFEM_MARK.
409 : */
410 : std::unordered_map<const Elem *, Xfem::CutElemInfo> _old_geom_cut_elems;
411 :
412 : /**
413 : * Store the solution in stored_solution for a given node
414 : * @param node_to_store_to Node for which the solution will be stored
415 : * @param node_to_store_from Node 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 storeSolutionForNode(const Node * node_to_store_to,
423 : const Node * node_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 : * Store the solution in stored_solution for a given element
432 : * @param elem_to_store_to Element for which the solution will be stored
433 : * @param elem_to_store_from Element from which the solution to be stored is obtained
434 : * @param sys System from which the solution is stored
435 : * @param stored_solution Data structure that the stored solution is saved to
436 : * @param current_solution Current solution vector that the solution is obtained from
437 : * @param old_solution Old solution vector that the solution is obtained from
438 : * @param older_solution Older solution vector that the solution is obtained from
439 : */
440 : void storeSolutionForElement(const Elem * elem_to_store_to,
441 : const Elem * elem_to_store_from,
442 : SystemBase & sys,
443 : std::map<unique_id_type, std::vector<Real>> & stored_solution,
444 : const NumericVector<Number> & current_solution,
445 : const NumericVector<Number> & old_solution,
446 : const NumericVector<Number> & older_solution);
447 :
448 : /**
449 : * Set the solution for all locally-owned nodes/elements that have stored values
450 : * @param sys System for which the solution is set
451 : * @param stored_solution Data structure that the stored solution is obtained from
452 : * @param current_solution Current solution vector that will be set
453 : * @param old_solution Old solution vector that will be set
454 : * @param older_solution Older solution vector that will be set
455 : */
456 : void setSolution(SystemBase & sys,
457 : const std::map<unique_id_type, std::vector<Real>> & stored_solution,
458 : NumericVector<Number> & current_solution,
459 : NumericVector<Number> & old_solution,
460 : NumericVector<Number> & older_solution);
461 :
462 : /**
463 : * Set the solution for a set of DOFs
464 : * @param stored_solution Stored solution values to set the solution to
465 : * @param stored_solution_dofs Dof indices for the entries in stored_solution
466 : * @param current_solution Current solution vector that will be set
467 : * @param old_solution Old solution vector that will be set
468 : * @param older_solution Older solution vector that will be set
469 : */
470 : void setSolutionForDOFs(const std::vector<Real> & stored_solution,
471 : const std::vector<dof_id_type> & stored_solution_dofs,
472 : NumericVector<Number> & current_solution,
473 : NumericVector<Number> & old_solution,
474 : NumericVector<Number> & older_solution);
475 :
476 : /**
477 : * Get a vector of the dof indices for all components of all variables
478 : * associated with an element
479 : * @param elem Element for which dof indices are found
480 : * @param sys System for which the dof indices are found
481 : */
482 : std::vector<dof_id_type> getElementSolutionDofs(const Elem * elem, SystemBase & sys) const;
483 :
484 : /**
485 : * Get a vector of the dof indices for all components of all variables
486 : * associated with a node
487 : * @param node Node for which dof indices are found
488 : * @param sys System for which the dof indices are found
489 : */
490 : std::vector<dof_id_type> getNodeSolutionDofs(const Node * node, SystemBase & sys) const;
491 :
492 : /**
493 : * Get the GeometricCutUserObject associated with an element
494 : * @param elem The element
495 : * @return A constant pointer to the GeometricCutUserObject, nullptr if nothing found
496 : */
497 : const GeometricCutUserObject * getGeometricCutForElem(const Elem * elem) const;
498 :
499 : void storeMaterialPropertiesForElementHelper(const Elem * elem,
500 : MaterialPropertyStorage & storage);
501 :
502 : /**
503 : * Helper function to store the material properties of a healed element
504 : * @param parent_elem The parent element
505 : * @param elem1 The first child element
506 : * @param elem2 The second child element
507 : */
508 : void storeMaterialPropertiesForElement(const Elem * parent_elem, const Elem * child_elem);
509 :
510 : /**
511 : * Load the material properties
512 : * @param props_deserialized The material properties
513 : * @param props_serialized The serialized material properties
514 : *
515 : * This does very dirty things and writes back to MOOSE's stateful properties. It should
516 : * not do this in the future.
517 : */
518 : void loadMaterialPropertiesForElementHelper(const Elem * elem,
519 : const Xfem::CachedMaterialProperties & cached_props,
520 : MaterialPropertyStorage & storage) const;
521 :
522 : /**
523 : * Helper function to store the material properties of a healed element
524 : * @param elem The cut element to restore material properties to.
525 : * @param elem_from The element to copy material properties from.
526 : * @param cached_cei The material properties cache to use.
527 : */
528 : void loadMaterialPropertiesForElement(
529 : const Elem * elem,
530 : const Elem * elem_from,
531 : std::unordered_map<const Elem *, Xfem::CutElemInfo> & cached_cei) const;
532 :
533 : /**
534 : * Return the first node in the provided element that is found to be in the physical domain
535 : * @param e Constant pointer to the child element
536 : * @param e0 Constant pointer to the parent element whose nodes will be querried
537 : * @return A constant pointer to the node
538 : */
539 : const Node * pickFirstPhysicalNode(const Elem * e, const Elem * e0) const;
540 : };
|
