XFEM.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#pragma once
 
#include "ElementPairLocator.h"
#include "ElementFragmentAlgorithm.h"
#include "XFEMInterface.h"
#include "XFEMCrackGrowthIncrement2DCut.h"
 
#include "libmesh/vector_value.h"
#include "libmesh/quadrature.h"
 
#include "GeometricCutUserObject.h"
 
// Forward declarations
class SystemBase;
 
namespace Xfem
{
enum XFEM_CUTPLANE_QUANTITY
{
  ORIGIN_X,
  ORIGIN_Y,
  ORIGIN_Z,
  NORMAL_X,
  NORMAL_Y,
  NORMAL_Z
};
 
enum XFEM_QRULE
{
  VOLFRAC,
  MOMENT_FITTING,
  DIRECT
};
 
typedef std::array<std::unordered_map<unsigned int, std::string>, 2> CachedMaterialProperties;
 
struct CutElemInfo
{
  const Elem * _parent_elem;
  const GeometricCutUserObject * _geometric_cut;
  CutSubdomainID _cut_subdomain_id;
  CachedMaterialProperties _elem_material_properties;
  CachedMaterialProperties _bnd_material_properties;
  // TODO: add neighbor material properties
  // CachedMaterialProperties _neighbor_material_properties;
 
  CutElemInfo()
    : _parent_elem(nullptr),
      _geometric_cut(nullptr),
      _cut_subdomain_id(std::numeric_limits<CutSubdomainID>::max())
  {
  }
 
  CutElemInfo(const Elem * parent_elem,
              const GeometricCutUserObject * geometric_cut,
              CutSubdomainID cut_subdomain_id)
    : _parent_elem(parent_elem), _geometric_cut(geometric_cut), _cut_subdomain_id(cut_subdomain_id)
  {
  }
 
  // A new child element is said to be previously healed if an entry in the old CutElemInfo has the
  // same parent element ID, the same cut, AND the same cut subdomain ID.
  bool match(const CutElemInfo & rhs)
  {
    return _parent_elem == rhs._parent_elem && _geometric_cut == rhs._geometric_cut &&
           _cut_subdomain_id == rhs._cut_subdomain_id;
  }
};
} // namespace Xfem
 
class XFEMCutElem;
class XFEMCrackGrowthIncrement2DCut;
class EFANode;
class EFAEdge;
class EFAElement;
class EFAElement2D;
class EFAElement3D;
 
// ------------------------------------------------------------
// XFEM class definition
class XFEM : public XFEMInterface
{
public:
  explicit XFEM(const InputParameters & params);
 
  ~XFEM();
 
  void addGeometricCut(GeometricCutUserObject * geometric_cut);
 
  void addStateMarkedElem(unsigned int elem_id, RealVectorValue & normal);
  void addStateMarkedElem(unsigned int elem_id, RealVectorValue & normal, unsigned int marked_side);
  void addStateMarkedFrag(unsigned int elem_id, RealVectorValue & normal);
 
  void clearStateMarkedElems();
 
  void addGeomMarkedElem2D(const unsigned int elem_id,
                           const Xfem::GeomMarkedElemInfo2D geom_info,
                           const unsigned int interface_id);
 
  void addGeomMarkedElem3D(const unsigned int elem_id,
                           const Xfem::GeomMarkedElemInfo3D geom_info,
                           const unsigned int interface_id);
 
  void clearGeomMarkedElems();
 
  virtual bool update(Real time,
                      const std::vector<std::shared_ptr<NonlinearSystemBase>> & nl,
                      AuxiliarySystem & aux) override;
 
  virtual void initSolution(const std::vector<std::shared_ptr<NonlinearSystemBase>> & nl,
                            AuxiliarySystem & aux) override;
 
  void buildEFAMesh();
  bool markCuts(Real time);
  bool markCutEdgesByGeometry();
  bool markCutEdgesByState(Real time);
  bool markCutFacesByGeometry();
  bool markCutFacesByState();
  bool initCutIntersectionEdge(
      Point cut_origin, RealVectorValue cut_normal, Point & edge_p1, Point & edge_p2, Real & dist);
  bool cutMeshWithEFA(const std::vector<std::shared_ptr<NonlinearSystemBase>> & nl,
                      AuxiliarySystem & aux);
 
  bool healMesh();
 
  virtual bool updateHeal() override;
  Point getEFANodeCoords(EFANode * CEMnode,
                         EFAElement * CEMElem,
                         const Elem * elem,
                         MeshBase * displaced_mesh = nullptr) const;
 
  Real getPhysicalVolumeFraction(const Elem * elem) const;
 
  bool isPointInsidePhysicalDomain(const Elem * elem, const Point & point) const;
 
  Real getCutPlane(const Elem * elem,
                   const Xfem::XFEM_CUTPLANE_QUANTITY quantity,
                   unsigned int plane_id) const;
 
  bool isElemAtCrackTip(const Elem * elem) const;
  bool isElemCut(const Elem * elem, XFEMCutElem *& xfce) const;
  bool isElemCut(const Elem * elem) const;
  void getFragmentFaces(const Elem * elem,
                        std::vector<std::vector<Point>> & frag_faces,
                        bool displaced_mesh = false) const;
  void storeCrackTipOriginAndDirection();
 
  void correctCrackExtensionDirection(const Elem * elem,
                                      EFAElement2D * CEMElem,
                                      EFAEdge * orig_edge,
                                      Point normal,
                                      Point crack_tip_origin,
                                      Point crack_tip_direction,
                                      Real & distance_keep,
                                      unsigned int & edge_id_keep,
                                      Point & normal_keep);
 
  void getCrackTipOrigin(std::map<unsigned int, const Elem *> & elem_id_crack_tip,
                         std::vector<Point> & crack_front_points);
 
  Xfem::XFEM_QRULE & getXFEMQRule();
  void setXFEMQRule(std::string & xfem_qrule);
  void setCrackGrowthMethod(bool use_crack_growth_increment, Real crack_growth_increment);
 
  void setDebugOutputLevel(unsigned int debug_output_level);
 
  void setMinWeightMultiplier(Real min_weight_multiplier);
 
  virtual bool getXFEMWeights(MooseArray<Real> & weights,
                              const Elem * elem,
                              QBase * qrule,
                              const MooseArray<Point> & q_points) override;
  virtual bool getXFEMFaceWeights(MooseArray<Real> & weights,
                                  const Elem * elem,
                                  QBase * qrule,
                                  const MooseArray<Point> & q_points,
                                  unsigned int side) override;
 
  virtual const ElementPairLocator::ElementPairList * getXFEMCutElemPairs(unsigned int interface_id)
  {
    return &_sibling_elems[interface_id];
  }
 
  virtual const ElementPairLocator::ElementPairList *
  getXFEMDisplacedCutElemPairs(unsigned int interface_id)
  {
    return &_sibling_displaced_elems[interface_id];
  }
 
  virtual unsigned int getGeometricCutID(const GeometricCutUserObject * gcu)
  {
    return _geom_marker_id_map[gcu];
  };
 
  virtual void getXFEMIntersectionInfo(const Elem * elem,
                                       unsigned int plane_id,
                                       Point & normal,
                                       std::vector<Point> & intersectionPoints,
                                       bool displaced_mesh = false) const;
  virtual void getXFEMqRuleOnLine(std::vector<Point> & intersection_points,
                                  std::vector<Point> & quad_pts,
                                  std::vector<Real> & quad_wts) const;
  virtual void getXFEMqRuleOnSurface(std::vector<Point> & intersection_points,
                                     std::vector<Point> & quad_pts,
                                     std::vector<Real> & quad_wts) const;
  bool has_secondary_cut() { return _has_secondary_cut; }
 
  EFAElement2D * getEFAElem2D(const Elem * elem);
 
  EFAElement3D * getEFAElem3D(const Elem * elem);
 
  void getFragmentEdges(const Elem * elem,
                        EFAElement2D * CEMElem,
                        std::vector<std::vector<Point>> & frag_edges) const;
  void getFragmentFaces(const Elem * elem,
                        EFAElement3D * CEMElem,
                        std::vector<std::vector<Point>> & frag_faces) const;
 
  const std::map<const Elem *, std::vector<Point>> & getCrackTipOriginMap() const
  {
    return _elem_crack_origin_direction_map;
  }
 
  CutSubdomainID getCutSubdomainID(const GeometricCutUserObject * gcuo,
                                   const Elem * cut_elem,
                                   const Elem * parent_elem = nullptr) const;
 
private:
  bool _has_secondary_cut;
 
  Xfem::XFEM_QRULE _XFEM_qrule;
 
  bool _use_crack_growth_increment;
  Real _crack_growth_increment;
 
  std::vector<const GeometricCutUserObject *> _geometric_cuts;
 
  std::map<unique_id_type, XFEMCutElem *> _cut_elem_map;
  std::set<const Elem *> _crack_tip_elems;
  std::set<const Elem *> _crack_tip_elems_to_be_healed;
  std::map<unsigned int, ElementPairLocator::ElementPairList> _sibling_elems;
  std::map<unsigned int, ElementPairLocator::ElementPairList> _sibling_displaced_elems;
 
  std::map<const Elem *, std::vector<Point>> _elem_crack_origin_direction_map;
 
  // std::map<const Elem*, Point> _crack_propagation_direction_map;
 
  std::map<const Elem *, RealVectorValue> _state_marked_elems;
  std::set<const Elem *> _state_marked_frags;
  std::map<const Elem *, unsigned int> _state_marked_elem_sides;
 
  int _state_marked_t_step;
 
  std::map<const Elem *, std::vector<Xfem::GeomMarkedElemInfo2D>> _geom_marked_elems_2d;
 
  std::map<const Elem *, std::vector<Xfem::GeomMarkedElemInfo3D>> _geom_marked_elems_3d;
 
  std::map<unsigned int, std::set<unsigned int>> _geom_marker_id_elems;
 
  std::map<const GeometricCutUserObject *, unsigned int> _geom_marker_id_map;
 
  ElementFragmentAlgorithm _efa_mesh;
 
  unsigned int _debug_output_level;
 
  Real _min_weight_multiplier;
 
  std::map<unique_id_type, std::vector<Real>> _cached_solution;
 
  std::map<unique_id_type, std::vector<Real>> _cached_aux_solution;
 
  std::unordered_map<const Elem *, Xfem::CutElemInfo> _geom_cut_elems;
 
  std::unordered_map<const Elem *, Xfem::CutElemInfo> _old_geom_cut_elems;
 
  void storeSolutionForNode(const Node * node_to_store_to,
                            const Node * node_to_store_from,
                            SystemBase & sys,
                            std::map<unique_id_type, std::vector<Real>> & stored_solution,
                            const NumericVector<Number> & current_solution,
                            const NumericVector<Number> & old_solution,
                            const NumericVector<Number> & older_solution);
 
  void storeSolutionForElement(const Elem * elem_to_store_to,
                               const Elem * elem_to_store_from,
                               SystemBase & sys,
                               std::map<unique_id_type, std::vector<Real>> & stored_solution,
                               const NumericVector<Number> & current_solution,
                               const NumericVector<Number> & old_solution,
                               const NumericVector<Number> & older_solution);
 
  void setSolution(SystemBase & sys,
                   const std::map<unique_id_type, std::vector<Real>> & stored_solution,
                   NumericVector<Number> & current_solution,
                   NumericVector<Number> & old_solution,
                   NumericVector<Number> & older_solution);
 
  void setSolutionForDOFs(const std::vector<Real> & stored_solution,
                          const std::vector<dof_id_type> & stored_solution_dofs,
                          NumericVector<Number> & current_solution,
                          NumericVector<Number> & old_solution,
                          NumericVector<Number> & older_solution);
 
  std::vector<dof_id_type> getElementSolutionDofs(const Elem * elem, SystemBase & sys) const;
 
  std::vector<dof_id_type> getNodeSolutionDofs(const Node * node, SystemBase & sys) const;
 
  const GeometricCutUserObject * getGeometricCutForElem(const Elem * elem) const;
 
  void storeMaterialPropertiesForElementHelper(const Elem * elem,
                                               MaterialPropertyStorage & storage);
 
  void storeMaterialPropertiesForElement(const Elem * parent_elem, const Elem * child_elem);
 
  void loadMaterialPropertiesForElementHelper(const Elem * elem,
                                              const Xfem::CachedMaterialProperties & cached_props,
                                              MaterialPropertyStorage & storage) const;
 
  void loadMaterialPropertiesForElement(
      const Elem * elem,
      const Elem * elem_from,
      std::unordered_map<const Elem *, Xfem::CutElemInfo> & cached_cei) const;
 
  const Node * pickFirstPhysicalNode(const Elem * e, const Elem * e0) const;
};