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 553 : virtual const ElementPairLocator::ElementPairList * getXFEMCutElemPairs(unsigned int interface_id)
253 : {
254 553 : 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 237 : getXFEMDisplacedCutElemPairs(unsigned int interface_id)
265 : {
266 237 : 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 38746 : 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 : private:
327 : bool _has_secondary_cut;
328 :
329 : Xfem::XFEM_QRULE _XFEM_qrule;
330 :
331 : bool _use_crack_growth_increment;
332 : Real _crack_growth_increment;
333 :
334 : std::vector<const GeometricCutUserObject *> _geometric_cuts;
335 :
336 : std::map<unique_id_type, XFEMCutElem *> _cut_elem_map;
337 : std::set<const Elem *> _crack_tip_elems;
338 : std::set<const Elem *> _crack_tip_elems_to_be_healed;
339 : std::map<unsigned int, ElementPairLocator::ElementPairList> _sibling_elems;
340 : std::map<unsigned int, ElementPairLocator::ElementPairList> _sibling_displaced_elems;
341 :
342 : std::map<const Elem *, std::vector<Point>> _elem_crack_origin_direction_map;
343 :
344 : // std::map<const Elem*, Point> _crack_propagation_direction_map;
345 :
346 : std::map<const Elem *, RealVectorValue> _state_marked_elems;
347 : std::set<const Elem *> _state_marked_frags;
348 : std::map<const Elem *, unsigned int> _state_marked_elem_sides;
349 :
350 : /// Data structure for storing information about all 2D elements to be cut by geometry
351 : std::map<const Elem *, std::vector<Xfem::GeomMarkedElemInfo2D>> _geom_marked_elems_2d;
352 :
353 : /// Data structure for storing information about all 3D elements to be cut by geometry
354 : std::map<const Elem *, std::vector<Xfem::GeomMarkedElemInfo3D>> _geom_marked_elems_3d;
355 :
356 : /// Data structure for storing the elements cut by specific geometric cutters
357 : std::map<unsigned int, std::set<unsigned int>> _geom_marker_id_elems;
358 :
359 : /// Data structure for storing the GeommetricCutUserObjects and their corresponding id
360 : std::map<const GeometricCutUserObject *, unsigned int> _geom_marker_id_map;
361 :
362 : ElementFragmentAlgorithm _efa_mesh;
363 :
364 : /// Controls amount of debugging output information
365 : /// 0: None
366 : /// 1: Summary
367 : /// 2: Details on modifications to mesh
368 : /// 3: Full dump of element fragment algorithm mesh
369 : unsigned int _debug_output_level;
370 :
371 : /// The minimum average multiplier applied by XFEM to the standard quadrature weights
372 : /// to integrate partial elements
373 : Real _min_weight_multiplier;
374 :
375 : /**
376 : * Data structure to store the nonlinear solution for nodes/elements affected by XFEM
377 : * For each node/element, this is stored as a vector that contains all components
378 : * of all applicable variables in an order defined by getElementSolutionDofs() or
379 : * getNodeSolutionDofs(). This vector first contains the current solution in that
380 : * order, followed by the old and older solutions, also in that same order.
381 : */
382 : std::map<unique_id_type, std::vector<Real>> _cached_solution;
383 :
384 : /**
385 : * Data structure to store the auxiliary solution for nodes/elements affected by XFEM
386 : * For each node/element, this is stored as a vector that contains all components
387 : * of all applicable variables in an order defined by getElementSolutionDofs() or
388 : * getNodeSolutionDofs(). This vector first contains the current solution in that
389 : * order, followed by the old and older solutions, also in that same order.
390 : */
391 : std::map<unique_id_type, std::vector<Real>> _cached_aux_solution;
392 :
393 : /**
394 : * All geometrically cut elements and their CutElemInfo during the current execution of
395 : * XFEM_MARK. This data structure is updated everytime a new cut element is created.
396 : */
397 : std::unordered_map<const Elem *, Xfem::CutElemInfo> _geom_cut_elems;
398 :
399 : /**
400 : * All geometrically cut elements and their CutElemInfo before the current execution of
401 : * XFEM_MARK.
402 : */
403 : std::unordered_map<const Elem *, Xfem::CutElemInfo> _old_geom_cut_elems;
404 :
405 : /**
406 : * Store the solution in stored_solution for a given node
407 : * @param node_to_store_to Node for which the solution will be stored
408 : * @param node_to_store_from Node from which the solution to be stored is obtained
409 : * @param sys System from which the solution is stored
410 : * @param stored_solution Data structure that the stored solution is saved to
411 : * @param current_solution Current solution vector that the solution is obtained from
412 : * @param old_solution Old solution vector that the solution is obtained from
413 : * @param older_solution Older solution vector that the solution is obtained from
414 : */
415 : void storeSolutionForNode(const Node * node_to_store_to,
416 : const Node * node_to_store_from,
417 : SystemBase & sys,
418 : std::map<unique_id_type, std::vector<Real>> & stored_solution,
419 : const NumericVector<Number> & current_solution,
420 : const NumericVector<Number> & old_solution,
421 : const NumericVector<Number> & older_solution);
422 :
423 : /**
424 : * Store the solution in stored_solution for a given element
425 : * @param elem_to_store_to Element for which the solution will be stored
426 : * @param elem_to_store_from Element from which the solution to be stored is obtained
427 : * @param sys System from which the solution is stored
428 : * @param stored_solution Data structure that the stored solution is saved to
429 : * @param current_solution Current solution vector that the solution is obtained from
430 : * @param old_solution Old solution vector that the solution is obtained from
431 : * @param older_solution Older solution vector that the solution is obtained from
432 : */
433 : void storeSolutionForElement(const Elem * elem_to_store_to,
434 : const Elem * elem_to_store_from,
435 : SystemBase & sys,
436 : std::map<unique_id_type, std::vector<Real>> & stored_solution,
437 : const NumericVector<Number> & current_solution,
438 : const NumericVector<Number> & old_solution,
439 : const NumericVector<Number> & older_solution);
440 :
441 : /**
442 : * Set the solution for all locally-owned nodes/elements that have stored values
443 : * @param sys System for which the solution is set
444 : * @param stored_solution Data structure that the stored solution is obtained from
445 : * @param current_solution Current solution vector that will be set
446 : * @param old_solution Old solution vector that will be set
447 : * @param older_solution Older solution vector that will be set
448 : */
449 : void setSolution(SystemBase & sys,
450 : const std::map<unique_id_type, std::vector<Real>> & stored_solution,
451 : NumericVector<Number> & current_solution,
452 : NumericVector<Number> & old_solution,
453 : NumericVector<Number> & older_solution);
454 :
455 : /**
456 : * Set the solution for a set of DOFs
457 : * @param stored_solution Stored solution values to set the solution to
458 : * @param stored_solution_dofs Dof indices for the entries in stored_solution
459 : * @param current_solution Current solution vector that will be set
460 : * @param old_solution Old solution vector that will be set
461 : * @param older_solution Older solution vector that will be set
462 : */
463 : void setSolutionForDOFs(const std::vector<Real> & stored_solution,
464 : const std::vector<dof_id_type> & stored_solution_dofs,
465 : NumericVector<Number> & current_solution,
466 : NumericVector<Number> & old_solution,
467 : NumericVector<Number> & older_solution);
468 :
469 : /**
470 : * Get a vector of the dof indices for all components of all variables
471 : * associated with an element
472 : * @param elem Element for which dof indices are found
473 : * @param sys System for which the dof indices are found
474 : */
475 : std::vector<dof_id_type> getElementSolutionDofs(const Elem * elem, SystemBase & sys) const;
476 :
477 : /**
478 : * Get a vector of the dof indices for all components of all variables
479 : * associated with a node
480 : * @param node Node for which dof indices are found
481 : * @param sys System for which the dof indices are found
482 : */
483 : std::vector<dof_id_type> getNodeSolutionDofs(const Node * node, SystemBase & sys) const;
484 :
485 : /**
486 : * Get the GeometricCutUserObject associated with an element
487 : * @param elem The element
488 : * @return A constant pointer to the GeometricCutUserObject, nullptr if nothing found
489 : */
490 : const GeometricCutUserObject * getGeometricCutForElem(const Elem * elem) const;
491 :
492 : void storeMaterialPropertiesForElementHelper(const Elem * elem,
493 : MaterialPropertyStorage & storage);
494 :
495 : /**
496 : * Helper function to store the material properties of a healed element
497 : * @param parent_elem The parent element
498 : * @param elem1 The first child element
499 : * @param elem2 The second child element
500 : */
501 : void storeMaterialPropertiesForElement(const Elem * parent_elem, const Elem * child_elem);
502 :
503 : /**
504 : * Load the material properties
505 : * @param props_deserialized The material properties
506 : * @param props_serialized The serialized material properties
507 : *
508 : * This does very dirty things and writes back to MOOSE's stateful properties. It should
509 : * not do this in the future.
510 : */
511 : void loadMaterialPropertiesForElementHelper(const Elem * elem,
512 : const Xfem::CachedMaterialProperties & cached_props,
513 : MaterialPropertyStorage & storage) const;
514 :
515 : /**
516 : * Helper function to store the material properties of a healed element
517 : * @param elem The cut element to restore material properties to.
518 : * @param elem_from The element to copy material properties from.
519 : * @param cached_cei The material properties cache to use.
520 : */
521 : void loadMaterialPropertiesForElement(
522 : const Elem * elem,
523 : const Elem * elem_from,
524 : std::unordered_map<const Elem *, Xfem::CutElemInfo> & cached_cei) const;
525 :
526 : /**
527 : * Return the first node in the provided element that is found to be in the physical domain
528 : * @param e Constant pointer to the child element
529 : * @param e0 Constant pointer to the parent element whose nodes will be querried
530 : * @return A constant pointer to the node
531 : */
532 : const Node * pickFirstPhysicalNode(const Elem * e, const Elem * e0) const;
533 : };
|
