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;

  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;
};