EshelbianPlasticity::EshelbianCore Struct Reference

#include "users_modules/eshelbian_plasticity/src/EshelbianPlasticity.hpp"

Inheritance diagram for EshelbianPlasticity::EshelbianCore:

Collaboration diagram for EshelbianPlasticity::EshelbianCore:

Classes
struct	DynamicRelaxationTimeScale
struct	SetUpSchur

Public Member Functions
MoFEMErrorCode	query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
	Getting interface of core database.
	EshelbianCore (MoFEM::Interface &m_field)
virtual	~EshelbianCore ()
MoFEMErrorCode	getOptions ()
template<typename BC >
MoFEMErrorCode	getBc (boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
MoFEMErrorCode	getSpatialDispBc ()
	[Getting norms]
MoFEMErrorCode	getSpatialRotationBc ()
MoFEMErrorCode	getSpatialTractionBc ()
MoFEMErrorCode	getTractionFreeBc (const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string contact_set_name)
	Remove all, but entities where kinematic constrains are applied.
MoFEMErrorCode	getSpatialTractionFreeBc (const EntityHandle meshset=0)
MoFEMErrorCode	getExternalStrain ()
MoFEMErrorCode	createExchangeVectors (Sev sev)
MoFEMErrorCode	addFields (const EntityHandle meshset=0)
MoFEMErrorCode	projectGeometry (const EntityHandle meshset=0)
MoFEMErrorCode	addVolumeFiniteElement (const EntityHandle meshset=0)
MoFEMErrorCode	addBoundaryFiniteElement (const EntityHandle meshset=0)
MoFEMErrorCode	addDMs (const BitRefLevel bit=BitRefLevel().set(0), const EntityHandle meshset=0)
MoFEMErrorCode	addMaterial_HMHHStVenantKirchhoff (const int tape, const double lambda, const double mu, const double sigma_y)
MoFEMErrorCode	addMaterial_HMHMooneyRivlin (const int tape, const double alpha, const double beta, const double lambda, const double sigma_y)
MoFEMErrorCode	addMaterial_HMHNeohookean (const int tape, const double c10, const double K)
MoFEMErrorCode	addMaterial_Hencky (double E, double nu)
MoFEMErrorCode	setBaseVolumeElementOps (const int tag, const bool do_rhs, const bool do_lhs, const bool calc_rates, SmartPetscObj< Vec > ver_vec, boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe)
MoFEMErrorCode	setVolumeElementOps (const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_lhs)
MoFEMErrorCode	setFaceElementOps (const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
MoFEMErrorCode	setContactElementRhsOps (boost::shared_ptr< ContactTree > &fe_contact_tree)
MoFEMErrorCode	setElasticElementOps (const int tag)
MoFEMErrorCode	setElasticElementToTs (DM dm)
MoFEMErrorCode	solveElastic (TS ts, Vec x)
MoFEMErrorCode	solveDynamicRelaxation (TS ts, Vec x)
MoFEMErrorCode	setBlockTagsOnSkin ()
MoFEMErrorCode	postProcessResults (const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, Vec var_vec=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
MoFEMErrorCode	postProcessSkeletonResults (const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
MoFEMErrorCode	calculateCrackArea (boost::shared_ptr< double > area_ptr)
MoFEMErrorCode	gettingNorms ()
	[Getting norms]
MoFEMErrorCode	calculateFaceMaterialForce (const int tag, TS ts)
MoFEMErrorCode	calculateOrientation (const int tag, bool set_orientation)
MoFEMErrorCode	setNewFrontCoordinates ()
MoFEMErrorCode	addCrackSurfaces (const bool debug=false)
MoFEMErrorCode	saveOrgCoords ()
MoFEMErrorCode	createCrackSurfaceMeshset ()
Public Member Functions inherited from MoFEM::UnknownInterface
virtual MoFEMErrorCode	query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
template<class IFACE >
MoFEMErrorCode	registerInterface (bool error_if_registration_failed=true)
	Register interface.
template<class IFACE >
MoFEMErrorCode	getInterface (IFACE *&iface) const
	Get interface reference to pointer of interface.
template<class IFACE >
MoFEMErrorCode	getInterface (IFACE **const iface) const
	Get interface pointer to pointer of interface.
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE	getInterface () const
	Get interface pointer to pointer of interface.
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE	getInterface () const
	Get reference to interface.
template<class IFACE >
IFACE *	getInterface () const
	Function returning pointer to interface.
virtual	~UnknownInterface ()=default

Static Public Member Functions
static double	f_log_e_quadratic (const double v)
static double	d_f_log_e_quadratic (const double v)
static double	dd_f_log_e_quadratic (const double v)
static double	f_log_e (const double v)
static double	d_f_log_e (const double v)
static double	dd_f_log_e (const double v)
static double	inv_f_log_e (const double v)
static double	inv_d_f_log_e (const double v)
static double	inv_dd_f_log_e (const double v)
static double	f_log (const double v)
static double	d_f_log (const double v)
static double	dd_f_log (const double v)
static double	inv_f_log (const double v)
static double	inv_d_f_log (const double v)
static double	inv_dd_f_log (const double v)
static double	f_linear (const double v)
static double	d_f_linear (const double v)
static double	dd_f_linear (const double v)
static double	inv_f_linear (const double v)
static double	inv_d_f_linear (const double v)
static double	inv_dd_f_linear (const double v)
Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode	getLibVersion (Version &version)
	Get library version.
static MoFEMErrorCode	getFileVersion (moab::Interface &moab, Version &version)
	Get database major version.
static MoFEMErrorCode	setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
	Get database major version.
static MoFEMErrorCode	getInterfaceVersion (Version &version)
	Get database major version.

Public Attributes
MoFEM::Interface &	mField
boost::shared_ptr< DataAtIntegrationPts >	dataAtPts
boost::shared_ptr< PhysicalEquations >	physicalEquations
boost::shared_ptr< AnalyticalExprPython >	AnalyticalExprPythonPtr
boost::shared_ptr< VolumeElementForcesAndSourcesCore >	elasticFeRhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore >	elasticFeLhs
boost::shared_ptr< FaceElementForcesAndSourcesCore >	elasticBcLhs
boost::shared_ptr< FaceElementForcesAndSourcesCore >	elasticBcRhs
boost::shared_ptr< ContactTree >	contactTreeRhs
	Make a contact tree.
SmartPetscObj< DM >	dM
	Coupled problem all fields.
SmartPetscObj< DM >	dmElastic
	Elastic problem.
SmartPetscObj< DM >	dmMaterial
	Material problem.
SmartPetscObj< DM >	dmPrjSpatial
	Projection spatial displacement.
const std::string	piolaStress = "P"
const std::string	spatialL2Disp = "wL2"
const std::string	spatialH1Disp = "wH1"
const std::string	materialH1Positions = "XH1"
const std::string	hybridSpatialDisp = "hybridSpatialDisp"
const std::string	contactDisp = "contactDisp"
const std::string	stretchTensor = "u"
const std::string	rotAxis = "omega"
const std::string	bubbleField = "bubble"
const std::string	elementVolumeName = "EP"
const std::string	naturalBcElement = "NATURAL_BC"
const std::string	skinElement = "SKIN"
const std::string	skeletonElement = "SKELETON"
const std::string	contactElement = "CONTACT"
int	spaceOrder = 2
int	spaceH1Order = -1
int	materialOrder = 1
double	alphaU = 0
double	alphaW = 0
double	alphaOmega = 0
double	alphaRho = 0
int	contactRefinementLevels = 1
int	frontLayers = 3
boost::shared_ptr< BcDispVec >	bcSpatialDispVecPtr
boost::shared_ptr< BcRotVec >	bcSpatialRotationVecPtr
boost::shared_ptr< TractionBcVec >	bcSpatialTractionVecPtr
boost::shared_ptr< TractionFreeBc >	bcSpatialFreeTractionVecPtr
boost::shared_ptr< NormalDisplacementBcVec >	bcSpatialNormalDisplacementVecPtr
boost::shared_ptr< AnalyticalDisplacementBcVec >	bcSpatialAnalyticalDisplacementVecPtr
boost::shared_ptr< AnalyticalTractionBcVec >	bcSpatialAnalyticalTractionVecPtr
boost::shared_ptr< PressureBcVec >	bcSpatialPressureVecPtr
boost::shared_ptr< ExternalStrainVec >	externalStrainVecPtr
std::map< std::string, boost::shared_ptr< ScalingMethod > >	timeScaleMap
boost::shared_ptr< Range >	contactFaces
boost::shared_ptr< Range >	crackFaces
boost::shared_ptr< Range >	frontEdges
boost::shared_ptr< Range >	frontAdjEdges
boost::shared_ptr< Range >	frontVertices
boost::shared_ptr< Range >	skeletonFaces
boost::shared_ptr< Range >	materialSkeletonFaces
boost::shared_ptr< Range >	maxMovedFaces
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > >	parentAdjSkeletonFunctionDim2
BitRefLevel	bitAdjParent = BitRefLevel().set()
	bit ref level for parent
BitRefLevel	bitAdjParentMask
	bit ref level for parent parent
BitRefLevel	bitAdjEnt = BitRefLevel().set()
	bit ref level for parent
BitRefLevel	bitAdjEntMask
	bit ref level for parent parent
SmartPetscObj< Vec >	solTSStep
CommInterface::EntitiesPetscVector	volumeExchange
CommInterface::EntitiesPetscVector	faceExchange
CommInterface::EntitiesPetscVector	edgeExchange
CommInterface::EntitiesPetscVector	vertexExchange
std::vector< Tag >	listTagsToTransfer
	list of tags to transfer to postprocessor
Mat	S = PETSC_NULLPTR
AO	aoS = PETSC_NULLPTR
SmartPetscObj< IS >	crackHybridIs
std::vector< std::string >	a00FieldList
std::vector< boost::shared_ptr< Range > >	a00RangeList
int	nbCrackFaces = 0

Static Public Attributes
static enum SymmetrySelector	symmetrySelector = NOT_SYMMETRIC
static enum RotSelector	rotSelector = LARGE_ROT
static enum RotSelector	gradApproximator = LARGE_ROT
static enum StretchSelector	stretchSelector = LOG
static PetscBool	noStretch = PETSC_FALSE
static PetscBool	setSingularity = PETSC_TRUE
static PetscBool	dynamicRelaxation
static PetscBool	crackingOn = PETSC_FALSE
static double	crackingStartTime = 0
static int	nbJIntegralLevels
static double	dynamicTime
static int	dynamicStep
static PetscBool	l2UserBaseScale = PETSC_TRUE
static int	addCrackMeshsetId = 1000
static double	griffithEnergy = 1
	Griffith energy.
static enum EnergyReleaseSelector	energyReleaseSelector
static std::string	internalStressTagName
static int	internalStressInterpOrder
static PetscBool	internalStressVoigt
static double	exponentBase = exp(1)
static boost::function< double(const double)>	f = EshelbianCore::f_log_e
static boost::function< double(const double)>	d_f
static boost::function< double(const double)>	dd_f
static boost::function< double(const double)>	inv_f
static boost::function< double(const double)>	inv_d_f
static boost::function< double(const double)>	inv_dd_f
static constexpr bool	use_quadratic_exp = true
static constexpr double	v_max = 12
static constexpr double	v_min = -v_max

Friends
struct	solve_elastic_set_up

Detailed Description

Definition at line 13 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

EshelbianCore()

EshelbianCore::EshelbianCore

(

MoFEM::Interface &

m_field

)

Examples

EshelbianPlasticity.cpp.

Definition at line 913 of file EshelbianPlasticity.cpp.

                                                    : mField(m_field) {
  CHK_THROW_MESSAGE(getOptions(), "getOptions failed");
}

~EshelbianCore()

EshelbianCore::~EshelbianCore ( )

virtualdefault

Member Function Documentation

addBoundaryFiniteElement()

MoFEMErrorCode EshelbianCore::addBoundaryFiniteElement

(

const EntityHandle

meshset = 0

)

Definition at line 1593 of file EshelbianPlasticity.cpp.

                                                                  {
  MoFEMFunctionBegin;
 
  Range meshset_ents;
  CHKERR mField.get_moab().get_entities_by_handle(meshset, meshset_ents);
 
  auto set_fe_adjacency = [&](auto fe_name) {
    MoFEMFunctionBegin;
    parentAdjSkeletonFunctionDim2 =
        boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
            bitAdjParent, bitAdjParentMask, bitAdjEnt, bitAdjEntMask);
    CHKERR mField.modify_finite_element_adjacency_table(
        fe_name, MBTRI, *parentAdjSkeletonFunctionDim2);
    MoFEMFunctionReturn(0);
  };
 
  // set finite element fields
  auto add_field_to_fe = [this](const std::string fe,
                                const std::string field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, materialH1Positions);
    MoFEMFunctionReturn(0);
  };
 
  if (!mField.check_finite_element(naturalBcElement)) {
 
    Range natural_bc_elements;
    if (bcSpatialDispVecPtr) {
      for (auto &v : *bcSpatialDispVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialRotationVecPtr) {
      for (auto &v : *bcSpatialRotationVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialNormalDisplacementVecPtr) {
      for (auto &v : *bcSpatialNormalDisplacementVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialAnalyticalDisplacementVecPtr) {
      for (auto &v : *bcSpatialAnalyticalDisplacementVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialTractionVecPtr) {
      for (auto &v : *bcSpatialTractionVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialAnalyticalTractionVecPtr) {
      for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    if (bcSpatialPressureVecPtr) {
      for (auto &v : *bcSpatialPressureVecPtr) {
        natural_bc_elements.merge(v.faces);
      }
    }
    natural_bc_elements = intersect(natural_bc_elements, meshset_ents);
 
    CHKERR mField.add_finite_element(naturalBcElement, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
                                                     naturalBcElement);
    CHKERR add_field_to_fe(naturalBcElement, piolaStress);
    CHKERR add_field_to_fe(naturalBcElement, hybridSpatialDisp);
    CHKERR set_fe_adjacency(naturalBcElement);
    CHKERR mField.build_finite_elements(naturalBcElement);
  }
 
  auto get_skin = [&](auto &body_ents) {
    Skinner skin(&mField.get_moab());
    Range skin_ents;
    CHKERR skin.find_skin(0, body_ents, false, skin_ents);
    return skin_ents;
  };
 
  auto filter_true_skin = [&](auto &&skin) {
    Range boundary_ents;
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &boundary_ents);
    return boundary_ents;
  };
 
  if (!mField.check_finite_element(skinElement)) {
 
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM,
                                                       body_ents);
    auto skin = filter_true_skin(get_skin(body_ents));
 
    CHKERR mField.add_finite_element(skinElement, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       spatialH1Disp);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       materialH1Positions);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       contactDisp);
    CHKERR mField.modify_finite_element_add_field_data(skinElement,
                                                       hybridSpatialDisp);
    // CHKERR add_field_to_fe(skinElement, hybridSpatialDisp);
    // CHKERR add_field_to_fe(skinElement, hybridMaterialDisp);
 
    CHKERR mField.build_finite_elements(skinElement);
  }
 
  if (!mField.check_finite_element(contactElement)) {
    if (contactFaces) {
      MOFEM_LOG("EP", Sev::inform)
          << "Contact elements " << contactFaces->size();
      CHKERR mField.add_finite_element(contactElement, MF_ZERO);
      CHKERR mField.add_ents_to_finite_element_by_type(*contactFaces, MBTRI,
                                                       contactElement);
      CHKERR add_field_to_fe(contactElement, piolaStress);
      CHKERR add_field_to_fe(contactElement, hybridSpatialDisp);
      CHKERR add_field_to_fe(contactElement, contactDisp);
      CHKERR add_field_to_fe(contactElement, spatialH1Disp);
      CHKERR set_fe_adjacency(contactElement);
      CHKERR mField.build_finite_elements(contactElement);
    }
  }
 
  if (!mField.check_finite_element(skeletonElement)) {
    if (!skeletonFaces)
      SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
    MOFEM_LOG("EP", Sev::inform)
        << "Skeleton elements " << skeletonFaces->size();
    CHKERR mField.add_finite_element(skeletonElement, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(*skeletonFaces, MBTRI,
                                                     skeletonElement);
    CHKERR add_field_to_fe(skeletonElement, piolaStress);
    CHKERR add_field_to_fe(skeletonElement, hybridSpatialDisp);
    CHKERR add_field_to_fe(skeletonElement, contactDisp);
    CHKERR add_field_to_fe(skeletonElement, spatialH1Disp);
    CHKERR set_fe_adjacency(skeletonElement);
    CHKERR mField.build_finite_elements(skeletonElement);
  }
 
  //   if (!mField.check_finite_element(materialSkeletonElement)) {
 
  //     Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM -
  //     2);
 
  //     Range front_elements;
  //     for (auto l = 0; l != frontLayers; ++l) {
  //       Range front_elements_layer;
  //       CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM,
  //       true,
  //                                                front_elements_layer,
  //                                                moab::Interface::UNION);
  //       front_elements.merge(front_elements_layer);
  //       front_edges.clear();
  //       CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
  //                                                SPACE_DIM - 2, true,
  //                                                front_edges,
  //                                                moab::Interface::UNION);
  //     }
  //     Range body_ents;
  //     CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
  //     body_ents); Range front_elements_faces; CHKERR
  //     mField.get_moab().get_adjacencies(front_elements, SPACE_DIM - 1,
  //                                              true, front_elements_faces,
  //                                              moab::Interface::UNION);
  //     auto body_skin = filter_true_skin(get_skin(body_ents));
  //     auto front_elements_skin = filter_true_skin(get_skin(front_elements));
  //     Range material_skeleton_faces =
  //         subtract(front_elements_faces, front_elements_skin);
  //     material_skeleton_faces.merge(intersect(front_elements_skin,
  //     body_skin));
 
  // #ifndef NDEBUG
  //     CHKERR save_range(mField.get_moab(), "front_elements.vtk",
  //     front_elements); CHKERR save_range(mField.get_moab(),
  //     "front_elements_skin.vtk",
  //                       front_elements_skin);
  //     CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
  //                       material_skeleton_faces);
  // #endif
 
  //     CHKERR mField.add_finite_element(materialSkeletonElement, MF_ZERO);
  //     CHKERR mField.add_ents_to_finite_element_by_type(
  //         material_skeleton_faces, MBTRI, materialSkeletonElement);
  //     // CHKERR add_field_to_fe(materialSkeletonElement, eshelbyStress);
  //     // CHKERR add_field_to_fe(materialSkeletonElement, hybridMaterialDisp);
  //     CHKERR set_fe_adjacency(materialSkeletonElement);
  //     CHKERR mField.build_finite_elements(materialSkeletonElement);
  //   }
 
  MoFEMFunctionReturn(0);
}

addCrackSurfaces()

MoFEMErrorCode EshelbianCore::addCrackSurfaces

(

const bool

debug = false

)

Definition at line 5941 of file EshelbianPlasticity.cpp.

                                                               {
  MoFEMFunctionBegin;
 
  constexpr bool potential_crack_debug = false;
  if constexpr (potential_crack_debug) {
 
    auto add_ents = get_range_from_block(mField, "POTENTIAL", SPACE_DIM - 1);
    Range crack_front_verts;
    CHKERR mField.get_moab().get_connectivity(*frontEdges, crack_front_verts,
                                              true);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        crack_front_verts);
    Range crack_front_faces;
    CHKERR mField.get_moab().get_adjacencies(crack_front_verts, SPACE_DIM - 1,
                                             true, crack_front_faces,
                                             moab::Interface::UNION);
    crack_front_faces = intersect(crack_front_faces, add_ents);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        crack_front_faces);
    CHKERR mField.getInterface<MeshsetsManager>()->addEntitiesToMeshset(
        BLOCKSET, addCrackMeshsetId, crack_front_faces);
  }
 
  auto get_crack_faces = [&]() {
    if (maxMovedFaces) {
      return unite(*crackFaces, *maxMovedFaces);
    } else {
      return *crackFaces;
    }
  };
 
  auto get_extended_crack_faces = [&]() {
    auto get_faces_of_crack_front_verts = [&](auto crack_faces_org) {
      ParallelComm *pcomm =
          ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
 
      Range crack_faces;
 
      if (!pcomm->rank()) {
 
        auto get_nodes = [&](auto &&e) {
          Range nodes;
          CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
                         "get connectivity");
          return nodes;
        };
 
        auto get_adj = [&](auto &&e, auto dim,
                           auto t = moab::Interface::UNION) {
          Range adj;
          CHK_MOAB_THROW(
              mField.get_moab().get_adjacencies(e, dim, true, adj, t),
              "get adj");
          return adj;
        };
 
        Range body_ents;
        CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                           body_ents);
        auto body_skin = get_skin(mField, body_ents);
        auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
        auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
        auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
        auto front_block_nodes = get_nodes(front_block_edges);
 
        size_t s;
        do {
          s = crack_faces.size();
 
          auto crack_face_nodes = get_nodes(crack_faces_org);
          auto crack_faces_edges =
              get_adj(crack_faces_org, 1, moab::Interface::UNION);
 
          auto crack_skin = get_skin(mField, crack_faces_org);
          front_block_edges = subtract(front_block_edges, crack_skin);
          auto crack_skin_nodes = get_nodes(crack_skin);
          crack_skin_nodes.merge(front_block_nodes);
 
          auto crack_skin_faces =
              get_adj(crack_skin, 2, moab::Interface::UNION);
          crack_skin_faces =
              subtract(subtract(crack_skin_faces, crack_faces_org), body_skin);
 
          crack_faces = crack_faces_org;
          for (auto f : crack_skin_faces) {
            auto edges = intersect(
                get_adj(Range(f, f), 1, moab::Interface::UNION), crack_skin);
 
            // if other edge is part of body skin, e.g. crack punching through
            // body surface
            if (edges.size() == 2) {
              edges.merge(
                  intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
                            body_skin_edges));
            }
 
            if (edges.size() == 2) {
              auto edge_conn = get_nodes(Range(edges));
              auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
                                     crack_faces_org);
              if (faces.size() == 2) {
                auto edge0_conn = get_nodes(Range(edges[0], edges[0]));
                auto edge1_conn = get_nodes(Range(edges[1], edges[1]));
                auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
                                            crack_skin_nodes); // node at apex
                if (edges_conn.size() == 1) {
 
                  auto node_edges =
                      subtract(intersect(get_adj(edges_conn, 1,
                                                 moab::Interface::INTERSECT),
                                         crack_faces_edges),
                               crack_skin); // nodes on crack surface, but not
                                            // at the skin
 
                  if (node_edges.size()) {
                    FTENSOR_INDEX(SPACE_DIM, i);
                    FTensor::Tensor1<double, SPACE_DIM> t_v0;
                    CHKERR mField.get_moab().get_coords(edges_conn, &t_v0(0));
 
                    auto get_t_dir = [&](auto e_conn) {
                      auto other_node = subtract(e_conn, edges_conn);
                      FTensor::Tensor1<double, SPACE_DIM> t_dir;
                      CHKERR mField.get_moab().get_coords(other_node,
                                                          &t_dir(0));
                      t_dir(i) -= t_v0(i);
                      return t_dir;
                    };
 
                    FTensor::Tensor1<double, SPACE_DIM> t_ave_dir;
                    t_ave_dir(i) =
                        get_t_dir(edge0_conn)(i) + get_t_dir(edge1_conn)(i);
 
                    FTensor::Tensor1<double, SPACE_DIM> t_crack_surface_ave_dir;
                    t_crack_surface_ave_dir(i) = 0;
                    for (auto e : node_edges) {
                      auto e_conn = get_nodes(Range(e, e));
                      auto t_dir = get_t_dir(e_conn);
                      t_crack_surface_ave_dir(i) += t_dir(i);
                    }
 
                    auto dot = t_ave_dir(i) * t_crack_surface_ave_dir(i);
                    // ave edges is in opposite direction to crack surface, so
                    // thus crack is not turning back
                    if (dot < -std::numeric_limits<double>::epsilon()) {
                      crack_faces.insert(f);
                    }
                  } else {
                    crack_faces.insert(f);
                  }
                }
              }
            } else if (edges.size() == 3) {
              crack_faces.insert(f);
            }
 
            // if other edge is part of geometry edge, e.g. keyway
            if (edges.size() == 1) {
              edges.merge(
                  intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
                            geometry_edges));
              edges.merge(
                  intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
                            front_block_edges));
              if (edges.size() == 2) {
                crack_faces.insert(f);
                continue;
              }
            }
          }
 
          crack_faces_org = crack_faces;
 
        } while (s != crack_faces.size());
      };
 
      return crack_faces; // send_type(mField, crack_faces, MBTRI);
    };
 
    return get_faces_of_crack_front_verts(get_crack_faces());
  };
 
  if (debug) {
    CHKERR save_range(mField.get_moab(), "new_crack_surface_debug.vtk",
                      get_extended_crack_faces());
  }
 
  auto reconstruct_crack_faces = [&](auto crack_faces) {
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
 
    auto impl = [&]() {
      MoFEMFunctionBegin;
 
      Range new_crack_faces;
      if (!pcomm->rank()) {
 
        auto get_nodes = [&](auto &&e) {
          Range nodes;
          CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
                         "get connectivity");
          return nodes;
        };
 
        auto get_adj = [&](auto &&e, auto dim,
                           auto t = moab::Interface::UNION) {
          Range adj;
          CHK_MOAB_THROW(
              mField.get_moab().get_adjacencies(e, dim, true, adj, t),
              "get adj");
          return adj;
        };
 
        auto get_test_on_crack_surface = [&]() {
          auto crack_faces_nodes =
              get_nodes(crack_faces); // nodes on crac faces
          auto crack_faces_tets =
              get_adj(crack_faces_nodes, 3,
                      moab::Interface::UNION); // adjacent
                                               // tets to
                                               // crack
                                               // faces throug nodes
          auto crack_faces_tets_nodes =
              get_nodes(crack_faces_tets); // nodes on crack faces tets
          crack_faces_tets_nodes =
              subtract(crack_faces_tets_nodes, crack_faces_nodes);
          crack_faces_tets =
              subtract(crack_faces_tets, get_adj(crack_faces_tets_nodes, 3,
                                                 moab::Interface::UNION));
          new_crack_faces =
              get_adj(crack_faces_tets, 2,
                      moab::Interface::UNION); // adjacency faces to crack
                                               // faces through tets
          new_crack_faces.merge(crack_faces);  // add original crack faces
 
          return std::make_tuple(new_crack_faces, crack_faces_tets);
        };
 
        auto carck_faces_test_edges = [&](auto faces, auto tets) {
          auto adj_tets_faces = get_adj(tets, 2, moab::Interface::UNION);
          auto adj_faces_edges = get_adj(subtract(faces, adj_tets_faces), 1,
                                         moab::Interface::UNION);
          auto adj_tets_edges = get_adj(tets, 1, moab::Interface::UNION);
          auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
          auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
          adj_faces_edges.merge(geometry_edges);    // geometry edges
          adj_faces_edges.merge(front_block_edges); // front block edges
 
          auto boundary_tets_edges = intersect(adj_tets_edges, adj_faces_edges);
          auto boundary_test_nodes = get_nodes(boundary_tets_edges);
          auto boundary_test_nodes_edges =
              get_adj(boundary_test_nodes, 1, moab::Interface::UNION);
          auto boundary_test_nodes_edges_nodes = subtract(
              get_nodes(boundary_test_nodes_edges), boundary_test_nodes);
 
          boundary_tets_edges =
              subtract(boundary_test_nodes_edges,
                       get_adj(boundary_test_nodes_edges_nodes, 1,
                               moab::Interface::UNION));
 
          Range body_ents;
          CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                             body_ents);
          auto body_skin = get_skin(mField, body_ents);
 
          auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
          body_skin_edges = intersect(get_adj(tets, 1, moab::Interface::UNION),
                                      body_skin_edges);
          body_skin = intersect(body_skin, adj_tets_faces);
          body_skin_edges = subtract(
              body_skin_edges, get_adj(body_skin, 1, moab::Interface::UNION));
 
          save_range(mField.get_moab(), "body_skin_edges.vtk", body_skin_edges);
          for (auto e : body_skin_edges) {
            auto adj_tet = intersect(
                get_adj(Range(e, e), 3, moab::Interface::INTERSECT), tets);
            if (adj_tet.size() == 1) {
              boundary_tets_edges.insert(e);
            }
          }
 
          return boundary_tets_edges;
        };
 
        auto p = get_test_on_crack_surface();
        auto &[new_crack_faces, crack_faces_tets] = p;
 
        if (debug) {
          CHKERR save_range(mField.get_moab(), "hole_crack_faces_debug.vtk",
                            crack_faces);
          CHKERR save_range(mField.get_moab(), "new_crack_faces_debug.vtk",
                            new_crack_faces);
          CHKERR save_range(mField.get_moab(), "new_crack_tets_debug.vtk",
                            crack_faces_tets);
        }
 
        auto boundary_tets_edges =
            carck_faces_test_edges(new_crack_faces, crack_faces_tets);
        CHKERR save_range(mField.get_moab(), "boundary_tets_edges.vtk",
                          boundary_tets_edges);
 
        auto resolve_surface = [&](auto boundary_tets_edges,
                                   auto crack_faces_tets) {
          auto boundary_tets_edges_nodes = get_nodes(boundary_tets_edges);
          auto crack_faces_tets_faces =
              get_adj(crack_faces_tets, 2, moab::Interface::UNION);
 
          Range all_removed_faces;
          Range all_removed_tets;
          int counter = 0;
 
          int size = 0;
          while (size != crack_faces_tets.size()) {
            auto tets_faces =
                get_adj(crack_faces_tets, 2, moab::Interface::UNION);
            auto skin_tets = get_skin(mField, crack_faces_tets);
            auto skin_skin =
                get_skin(mField, subtract(crack_faces_tets_faces, tets_faces));
            auto skin_skin_nodes = get_nodes(skin_skin);
 
            size = crack_faces_tets.size();
            MOFEM_LOG("SELF", Sev::inform)
                << "Crack faces tets size " << crack_faces_tets.size()
                << " crack faces size " << crack_faces_tets_faces.size();
            auto skin_tets_nodes = subtract(
                get_nodes(skin_tets),
                boundary_tets_edges_nodes); // not remove tets which are
                                            // adjagasent to crack faces nodes
            skin_tets_nodes = subtract(skin_tets_nodes, skin_skin_nodes);
 
            Range removed_nodes;
            Range tets_to_remove;
            Range faces_to_remove;
            for (auto n : skin_tets_nodes) {
              auto tets =
                  intersect(get_adj(Range(n, n), 3, moab::Interface::INTERSECT),
                            crack_faces_tets);
              if (tets.size() == 0) {
                continue;
              }
 
              auto hole_detetction = [&]() {
                auto adj_tets =
                    get_adj(Range(n, n), 3, moab::Interface::INTERSECT);
                adj_tets =
                    subtract(adj_tets,
                             crack_faces_tets); // tetst adjacent to the node
                                                // but not part of crack surface
                if (adj_tets.size() == 0) {
                  return std::make_pair(
                      intersect(
                          get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
                          tets_faces),
                      tets);
                }
 
                std::vector<Range> tets_groups;
                auto test_adj_tets = adj_tets;
                while (test_adj_tets.size()) {
                  auto seed_size = 0;
                  Range seed = Range(test_adj_tets[0], test_adj_tets[0]);
                  while (seed.size() != seed_size) {
                    auto adj_faces =
                        subtract(get_adj(seed, 2, moab::Interface::UNION),
                                 tets_faces); // edges which are not
                                              // part of the node
                    seed_size = seed.size();
                    seed.merge(
                        intersect(get_adj(adj_faces, 3, moab::Interface::UNION),
                                  test_adj_tets ));
                  }
                  tets_groups.push_back(seed);
                  test_adj_tets = subtract(test_adj_tets, seed);
                }
                if (tets_groups.size() == 1) {
 
                  return std::make_pair(
                      intersect(
                          get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
                          tets_faces),
                      tets);
 
                }
 
                Range tets_to_remove;
                Range faces_to_remove;
                for (auto &r : tets_groups) {
                  auto f = get_adj(r, 2, moab::Interface::UNION);
                  auto t = intersect(get_adj(f, 3, moab::Interface::UNION),
                                     crack_faces_tets); // tets
 
                  if (f.size() > faces_to_remove.size() ||
                      faces_to_remove.size() == 0) {
                    faces_to_remove = f;
                    tets_to_remove = t; // largest group of tets
                  }
                }
                MOFEM_LOG("EPSELF", Sev::inform)
                    << "Hole detection: faces to remove "
                    << faces_to_remove.size() << " tets to remove "
                    << tets_to_remove.size();
                return std::make_pair(faces_to_remove, tets_to_remove);
              };
 
              if (tets.size() < tets_to_remove.size() ||
                  tets_to_remove.size() == 0) {
                removed_nodes = Range(n, n);
                auto [h_faces_to_remove, h_tets_to_remove] =
                    hole_detetction(); // find faces and tets to remove
                faces_to_remove = h_faces_to_remove;
                tets_to_remove = h_tets_to_remove;
 
                // intersect(
                //     get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
                //     tets_faces);
 
              } // find tets which is largest adjacencty size, so that it is
                // removed first, and then faces are removed
              all_removed_faces.merge(faces_to_remove);
              all_removed_tets.merge(tets_to_remove);
            }
 
            crack_faces_tets = subtract(crack_faces_tets, tets_to_remove);
            crack_faces_tets_faces =
                subtract(crack_faces_tets_faces, faces_to_remove);
 
            if (debug) {
              save_range(mField.get_moab(),
                         "removed_nodes_" +
                             boost::lexical_cast<std::string>(counter) + ".vtk",
                         removed_nodes);
              save_range(mField.get_moab(),
                         "faces_to_remove_" +
                             boost::lexical_cast<std::string>(counter) + ".vtk",
                         faces_to_remove);
              save_range(mField.get_moab(),
                         "tets_to_remove_" +
                             boost::lexical_cast<std::string>(counter) + ".vtk",
                         tets_to_remove);
              save_range(mField.get_moab(),
                         "crack_faces_tets_faces_" +
                             boost::lexical_cast<std::string>(counter) + ".vtk",
                         crack_faces_tets_faces);
              save_range(mField.get_moab(),
                         "crack_faces_tets_" +
                             boost::lexical_cast<std::string>(counter) + ".vtk",
                         crack_faces_tets);
            }
            counter++;
          }
 
          auto cese_internal_faces = [&]() {
            MoFEMFunctionBegin;
            auto skin_tets = get_skin(mField, crack_faces_tets);
            auto adj_faces = get_adj(skin_tets, 2, moab::Interface::UNION);
            adj_faces =
                subtract(adj_faces, skin_tets); // remove skin tets faces
            auto adj_tets = get_adj(adj_faces, 3,
                                    moab::Interface::UNION); // tets which are
                                                             // adjacent to skin
            crack_faces_tets =
                subtract(crack_faces_tets,
                         adj_tets); // remove tets which are adjacent to
                                    // skin, so that they are not removed
            crack_faces_tets_faces =
                subtract(crack_faces_tets_faces, adj_faces);
 
            all_removed_faces.merge(adj_faces);
            all_removed_tets.merge(adj_tets);
            
 
            MOFEM_LOG("EPSELF", Sev::inform)
                << "Remove internal faces size " << adj_faces.size()
                << " tets size " << adj_tets.size();
            MoFEMFunctionReturn(0);
          };
 
          auto case_only_one_free_edge = [&]() {
            MoFEMFunctionBegin;
 
            for (auto t : Range(crack_faces_tets)) {
 
              auto adj_faces = get_adj(
                  Range(t, t), 2,
                  moab::Interface::UNION); // faces of tet which can be removed
              auto crack_surface_edges =
                  get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
                                   adj_faces),
                          1,
                          moab::Interface::UNION); // edges not on the tet but
                                                   // on crack surface
              auto adj_edges =
                  subtract(get_adj(Range(t, t), 1, moab::Interface::INTERSECT),
                           crack_surface_edges); // free edges
              adj_edges = subtract(
                  adj_edges,
                  boundary_tets_edges); // edges which are not part of gemetry
 
              if (adj_edges.size() == 1) {
                crack_faces_tets =
                    subtract(crack_faces_tets,
                             Range(t, t)); // remove tets which are adjacent to
                                           // skin, so that they are not removed
 
                auto faces_to_remove =
                    get_adj(adj_edges, 2, moab::Interface::UNION); // faces
                                                                   // which can
                                                                   // be removed
                crack_faces_tets_faces =
                    subtract(crack_faces_tets_faces, faces_to_remove);
 
                all_removed_faces.merge(faces_to_remove);
                all_removed_tets.merge(Range(t, t));
 
                MOFEM_LOG("EPSELF", Sev::inform) << "Remove free one edges ";
              }
            }
 
            crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
            crack_faces_tets_faces =
                subtract(crack_faces_tets_faces, all_removed_faces);
 
            MoFEMFunctionReturn(0);
          };
 
          auto cese_flat_tet = [&](auto max_adj_edges) {
            MoFEMFunctionBegin;
 
            Range body_ents;
            CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                               body_ents);
            auto body_skin = get_skin(mField, body_ents);
            auto body_skin_edges =
                get_adj(body_skin, 1, moab::Interface::UNION);
 
            for (auto t : Range(crack_faces_tets)) {
 
              auto adj_faces = get_adj(
                  Range(t, t), 2,
                  moab::Interface::UNION); // faces of tet which can be removed
              auto crack_surface_edges =
                  get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
                                   adj_faces),
                          1,
                          moab::Interface::UNION); // edges not on the tet but
                                                   // on crack surface
              auto adj_edges =
                  subtract(get_adj(Range(t, t), 1, moab::Interface::INTERSECT),
                           crack_surface_edges); // free edges
              adj_edges = subtract(adj_edges, body_skin_edges);
 
              auto tet_edges = get_adj(Range(t, t), 1,
                                          moab::Interface::UNION); // edges of
                                                                    // tet
              tet_edges = subtract(tet_edges, adj_edges);
 
              for (auto e : tet_edges) {
                constexpr int opposite_edge[] = {5, 3, 4, 1, 2, 0};
                auto get_side = [&](auto e) {
                  int side, sense, offset;
                  CHK_MOAB_THROW(
                      mField.get_moab().side_number(t, e, side, sense, offset),
                      "get side number failed");
                  return side;
                };
                auto get_side_ent = [&](auto side) {
                  EntityHandle side_edge;
                  CHK_MOAB_THROW(
                      mField.get_moab().side_element(t, 1, side, side_edge),
                      "get side failed");
                  return side_edge;
                };
                adj_edges.erase(get_side_ent(opposite_edge[get_side(e)]));
              }
 
              if (adj_edges.size() <= max_adj_edges) {
 
                double dot = 1;
                Range faces_to_remove;
                for (auto e : adj_edges) {
                  auto edge_adj_faces =
                      get_adj(Range(e, e), 2, moab::Interface::UNION);
                  edge_adj_faces = intersect(edge_adj_faces, adj_faces);
                  if (edge_adj_faces.size() != 2) {
                    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                                      "Adj faces size is not 2 for edge " +
                                          boost::lexical_cast<std::string>(e));
                  }
 
                  auto get_normal = [&](auto f) {
                    FTensor::Tensor1<double, SPACE_DIM> t_n;
                    CHK_THROW_MESSAGE(
                        mField.getInterface<Tools>()->getTriNormal(f, &t_n(0)),
                        "get tri normal failed");
                    return t_n;
                  };
                  auto t_n0 = get_normal(edge_adj_faces[0]);
                  auto t_n1 = get_normal(edge_adj_faces[1]);
                  auto get_sense = [&](auto f) {
                    int side, sense, offset;
                    CHK_MOAB_THROW(mField.get_moab().side_number(t, f, side,
                                                                 sense, offset),
                                   "get side number failed");
                    return sense;
                  };
                  auto sense0 = get_sense(edge_adj_faces[0]);
                  auto sense1 = get_sense(edge_adj_faces[1]);
                  t_n0.normalize();
                  t_n1.normalize();
 
                  FTENSOR_INDEX(SPACE_DIM, i);
                  auto dot_e = (sense0 * sense1) * t_n0(i) * t_n1(i);
                  if (dot_e < dot || e == adj_edges[0]) {
                    dot = dot_e;
                    faces_to_remove = edge_adj_faces;
                  }
                }
 
                all_removed_faces.merge(faces_to_remove);
                all_removed_tets.merge(Range(t, t));
 
                MOFEM_LOG("EPSELF", Sev::inform)
                    << "Remove free edges on flat tet, with considered nb. of "
                       "edges "
                    << adj_edges.size();
              }
            }
 
            crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
            crack_faces_tets_faces =
                subtract(crack_faces_tets_faces, all_removed_faces);
 
            MoFEMFunctionReturn(0);
          };
 
          CHK_THROW_MESSAGE(case_only_one_free_edge(),
                            "Case only one free edge failed");
          for (auto max_adj_edges : {0, 1, 2, 3}) {
            CHK_THROW_MESSAGE(cese_flat_tet(max_adj_edges),
                              "Case only one free edge failed");
          }
          CHK_THROW_MESSAGE(cese_internal_faces(),
                            "Case internal faces failed");
 
          if (debug) {
            save_range(mField.get_moab(),
                       "crack_faces_tets_faces_" +
                           boost::lexical_cast<std::string>(counter) + ".vtk",
                       crack_faces_tets_faces);
            save_range(mField.get_moab(),
                       "crack_faces_tets_" +
                           boost::lexical_cast<std::string>(counter) + ".vtk",
                       crack_faces_tets);
          }
 
          return std::make_tuple(crack_faces_tets_faces, crack_faces_tets,
                                 all_removed_faces, all_removed_tets);
        };
 
        auto [resolved_faces, resolved_tets, all_removed_faces,
              all_removed_tets] =
            resolve_surface(boundary_tets_edges, crack_faces_tets);
        resolved_faces.merge(subtract(crack_faces, all_removed_faces));
        if (debug) {
          CHKERR save_range(mField.get_moab(), "resolved_faces.vtk",
                            resolved_faces);
          CHKERR save_range(mField.get_moab(), "resolved_tets.vtk",
                            resolved_tets);
        }
 
        crack_faces = resolved_faces;
      }
 
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(impl(), "resolve new crack surfaces");
 
    return crack_faces; // send_type(mField, crack_faces, MBTRI);
  };
 
 
  auto resolve_consisten_crack_extension = [&]() {
    MoFEMFunctionBegin;
    auto crack_meshset =
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
            addCrackMeshsetId, BLOCKSET);
    auto meshset = crack_meshset->getMeshset();
 
 
    if (!mField.get_comm_rank()) {
      Range old_crack_faces;
      CHKERR mField.get_moab().get_entities_by_type(meshset, MBTRI,
                                                    old_crack_faces);
      auto extendeded_crack_faces = get_extended_crack_faces();
      auto reconstructed_crack_faces =
          subtract(reconstruct_crack_faces(extendeded_crack_faces),
                   subtract(*crackFaces, old_crack_faces));
      if (nbCrackFaces >= reconstructed_crack_faces.size()) {
        MOFEM_LOG("EPSELF", Sev::warning)
            << "No new crack faces to add, skipping adding to meshset";
        extendeded_crack_faces = subtract(
            extendeded_crack_faces, subtract(*crackFaces, old_crack_faces));
        MOFEM_LOG("EPSELF", Sev::inform)
            << "Number crack faces size (extended) "
            << extendeded_crack_faces.size();
        CHKERR mField.get_moab().clear_meshset(&meshset, 1);
        CHKERR mField.get_moab().add_entities(meshset, extendeded_crack_faces);
      } else {
        CHKERR mField.get_moab().clear_meshset(&meshset, 1);
        CHKERR mField.get_moab().add_entities(meshset,
                                              reconstructed_crack_faces);
        MOFEM_LOG("EPSELF", Sev::inform)
            << "Number crack faces size (reconstructed) "
            << reconstructed_crack_faces.size();
        nbCrackFaces = reconstructed_crack_faces.size();
      }
    }
 
    Range crack_faces;
    if (!mField.get_comm_rank()) {
      CHKERR mField.get_moab().get_entities_by_type(meshset, MBTRI,
                                                    crack_faces);
    }
    crack_faces = send_type(mField, crack_faces, MBTRI);
    if (mField.get_comm_rank()) {
      CHKERR mField.get_moab().clear_meshset(&meshset, 1);
      CHKERR mField.get_moab().add_entities(meshset, crack_faces);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR resolve_consisten_crack_extension();
 
  MoFEMFunctionReturn(0);
};

addDMs()

MoFEMErrorCode EshelbianCore::addDMs	(	const BitRefLevel	bit = BitRefLevel().set(0),
		const EntityHandle	meshset = 0 )

Definition at line 1793 of file EshelbianPlasticity.cpp.

                                                                 {
  MoFEMFunctionBegin;
 
  // find adjacencies between finite elements and dofs
  CHKERR mField.build_adjacencies(bit, QUIET);
 
  // Create coupled problem
  dM = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
                            BitRefLevel().set());
  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
  CHKERR DMMoFEMAddElement(dM, elementVolumeName);
  CHKERR DMMoFEMAddElement(dM, naturalBcElement);
  CHKERR DMMoFEMAddElement(dM, skeletonElement);
  CHKERR DMMoFEMAddElement(dM, contactElement);
  CHKERR DMMoFEMAddElement(dM, skinElement);
 
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
  CHKERR DMSetUp(dM);
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
 
  auto remove_dofs_on_broken_skin = [&](const std::string prb_name) {
    MoFEMFunctionBegin;
    for (int d : {0, 1, 2}) {
      std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
      CHKERR mField.getInterface<ProblemsManager>()
          ->getSideDofsOnBrokenSpaceEntities(
              dofs_to_remove, prb_name, ROW, piolaStress,
              (*bcSpatialFreeTractionVecPtr)[d], SPACE_DIM, d, d);
      // remove piola dofs, i.e. traction free boundary
      CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, ROW,
                                                                dofs_to_remove);
      CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, COL,
                                                                dofs_to_remove);
    }
    MoFEMFunctionReturn(0);
  };
  CHKERR remove_dofs_on_broken_skin("ESHELBY_PLASTICITY");
 
  // Create elastic sub-problem
  dmElastic = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
  CHKERR DMMoFEMSetDestroyProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, piolaStress);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis);
  if (!noStretch) {
    CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor);
  }
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, hybridSpatialDisp);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, contactDisp);
  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
  CHKERR DMMoFEMAddElement(dmElastic, skeletonElement);
  CHKERR DMMoFEMAddElement(dmElastic, contactElement);
  CHKERR DMMoFEMAddElement(dmElastic, skinElement);
  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
  CHKERR DMSetUp(dmElastic);
 
  // dmMaterial = createDM(mField.get_comm(), "DMMOFEM");
  // CHKERR DMMoFEMCreateSubDM(dmMaterial, dM, "MATERIAL_PROBLEM");
  // CHKERR DMMoFEMSetDestroyProblem(dmMaterial, PETSC_TRUE);
  // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, eshelbyStress);
  // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, materialL2Disp);
  // CHKERR DMMoFEMAddElement(dmMaterh elementVolumeName);
  // CHKERR DMMoFEMAddElement(dmMaterial, naturalBcElement);
  // CHKERR DMMoFEMAddElement(dmMaterial, skinElement);
  // CHKERR DMMoFEMSetSquareProblem(dmMaterial, PETSC_TRUE);
  // CHKERR DMMoFEMSetIsPartitioned(dmMaterial, PETSC_TRUE);
  // CHKERR DMSetUp(dmMaterial);
 
  auto set_zero_block = [&]() {
    MoFEMFunctionBegin;
    if (!noStretch) {
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
    }
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
    CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
        "ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
    if (!noStretch) {
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", bubbleField, bubbleField);
      CHKERR
      mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", piolaStress, piolaStress);
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", bubbleField, piolaStress);
      CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
          "ELASTIC_PROBLEM", piolaStress, bubbleField);
    }
 
    solTSStep = createDMVector(dmElastic);
    MoFEMFunctionReturn(0);
  };
 
  auto set_section = [&]() {
    MoFEMFunctionBegin;
    PetscSection section;
    CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
                                                           &section);
    CHKERR DMSetSection(dmElastic, section);
    CHKERR DMSetGlobalSection(dmElastic, section);
    CHKERR PetscSectionDestroy(&section);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_zero_block();
  CHKERR set_section();
 
  dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
  CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp);
  CHKERR DMMoFEMAddElement(dmPrjSpatial, elementVolumeName);
  CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMSetUp(dmPrjSpatial);
 
  // CHKERR mField.getInterface<BcManager>()
  //     ->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
  //         "PROJECT_SPATIAL", spatialH1Disp, true, false);
 
  MoFEMFunctionReturn(0);
}

addFields()

MoFEMErrorCode EshelbianCore::addFields

(

const EntityHandle

meshset = 0

)

Definition at line 1154 of file EshelbianPlasticity.cpp.

                                                                  {
  MoFEMFunctionBegin;
 
  auto get_tets = [&]() {
    Range tets;
    CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
    return tets;
  };
 
  auto get_tets_skin = [&]() {
    Range tets_skin_part;
    Skinner skin(&mField.get_moab());
    CHKERR skin.find_skin(0, get_tets(), false, tets_skin_part);
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    Range tets_skin;
    CHKERR pcomm->filter_pstatus(tets_skin_part,
                                 PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &tets_skin);
    return tets_skin;
  };
 
  auto subtract_boundary_conditions = [&](auto &&tets_skin) {
    // That mean, that hybrid field on all faces on which traction is applied,
    // on other faces, or enforcing displacements as
    // natural boundary condition.
    if (bcSpatialTractionVecPtr)
      for (auto &v : *bcSpatialTractionVecPtr) {
        tets_skin = subtract(tets_skin, v.faces);
      }
    if (bcSpatialAnalyticalTractionVecPtr)
      for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
        tets_skin = subtract(tets_skin, v.faces);
      }      
 
    if (bcSpatialPressureVecPtr)
      for (auto &v : *bcSpatialPressureVecPtr) {
        tets_skin = subtract(tets_skin, v.faces);
      }
 
    return tets_skin;
  };
 
  auto add_blockset = [&](auto block_name, auto &&tets_skin) {
    auto crack_faces =
        get_range_from_block(mField, "block_name", SPACE_DIM - 1);
    tets_skin.merge(crack_faces);
    return tets_skin;
  };
 
  auto subtract_blockset = [&](auto block_name, auto &&tets_skin) {
    auto contact_range =
        get_range_from_block(mField, block_name, SPACE_DIM - 1);
    tets_skin = subtract(tets_skin, contact_range);
    return tets_skin;
  };
 
  auto get_stress_trace_faces = [&](auto &&tets_skin) {
    Range faces;
    CHKERR mField.get_moab().get_adjacencies(get_tets(), SPACE_DIM - 1, true,
                                             faces, moab::Interface::UNION);
    Range trace_faces = subtract(faces, tets_skin);
    return trace_faces;
  };
 
  auto tets = get_tets();
 
  // remove also contact faces, i.e. that is also kind of hybrid field but
  // named but used to enforce contact conditions
  auto trace_faces = get_stress_trace_faces(
 
      subtract_blockset("CONTACT",
                        subtract_boundary_conditions(get_tets_skin()))
 
  );
 
  contactFaces = boost::make_shared<Range>(intersect(
      trace_faces, get_range_from_block(mField, "CONTACT", SPACE_DIM - 1)));
  skeletonFaces =
      boost::make_shared<Range>(subtract(trace_faces, *contactFaces));
  // materialSkeletonFaces =
  //     boost::make_shared<Range>(get_stress_trace_faces(Range()));
 
#ifndef NDEBUG
  if (contactFaces->size())
    CHKERR save_range(mField.get_moab(),
                      "contact_faces_" +
                          std::to_string(mField.get_comm_rank()) + ".vtk",
                      *contactFaces);
  if (skeletonFaces->size())
    CHKERR save_range(mField.get_moab(),
                      "skeleton_faces_" +
                          std::to_string(mField.get_comm_rank()) + ".vtk",
                      *skeletonFaces);
    // if (skeletonFaces->size())
    //   CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
    //                     *materialSkeletonFaces);
#endif
 
  auto add_broken_hdiv_field = [this, meshset](const std::string field_name,
                                               const int order) {
    MoFEMFunctionBegin;
 
    constexpr FieldApproximationBase base = DEMKOWICZ_JACOBI_BASE;
 
    auto get_side_map_hdiv = [&]() {
      return std::vector<
 
          std::pair<EntityType,
                    std::function<MoFEMErrorCode(BaseFunction::DofsSideMap &)>
 
                    >>{
 
          {MBTET,
           [&](BaseFunction::DofsSideMap &dofs_side_map) -> MoFEMErrorCode {
             return TetPolynomialBase::setDofsSideMap(HDIV, DISCONTINUOUS, base,
                                                      dofs_side_map);
           }}
 
      };
    };
 
    CHKERR mField.add_broken_field(field_name, HDIV, base, SPACE_DIM,
                                   get_side_map_hdiv(), MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_l2_field = [this, meshset](const std::string field_name,
                                      const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_h1_field = [this, meshset](const std::string field_name,
                                      const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
    CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_l2_field_by_range = [this](const std::string field_name,
                                      const int order, const int dim,
                                      const int field_dim, Range &&r) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, field_dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(r);
    CHKERR mField.add_ents_to_field_by_dim(r, dim, field_name);
    CHKERR mField.set_field_order(r, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_bubble_field = [this, meshset](const std::string field_name,
                                          const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
                            MF_ZERO);
    // Modify field
    auto field_ptr = mField.get_field_structure(field_name);
    auto field_order_table =
        const_cast<Field *>(field_ptr)->getFieldOrderTable();
    auto get_cgg_bubble_order_zero = [](int p) { return 0; };
    auto get_cgg_bubble_order_tet = [](int p) {
      return NBVOLUMETET_CCG_BUBBLE(p);
    };
    field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
    field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
    field_order_table[MBTRI] = get_cgg_bubble_order_zero;
    field_order_table[MBTET] = get_cgg_bubble_order_tet;
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_user_l2_field = [this, meshset](const std::string field_name,
                                           const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, USER_BASE, dim, MB_TAG_DENSE,
                            MF_ZERO);
    // Modify field
    auto field_ptr = mField.get_field_structure(field_name);
    auto field_order_table =
        const_cast<Field *>(field_ptr)->getFieldOrderTable();
    auto zero_dofs = [](int p) { return 0; };
    auto dof_l2_tet = [](int p) { return NBVOLUMETET_L2(p); };
    field_order_table[MBVERTEX] = zero_dofs;
    field_order_table[MBEDGE] = zero_dofs;
    field_order_table[MBTRI] = zero_dofs;
    field_order_table[MBTET] = dof_l2_tet;
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  // spatial fields
  CHKERR add_broken_hdiv_field(piolaStress, spaceOrder);
  CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
  CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
  CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
  CHKERR add_user_l2_field(stretchTensor, noStretch ? -1 : spaceOrder, 6);
 
  if (!skeletonFaces)
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
  if (!contactFaces)
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No contact faces");
 
  auto get_hybridised_disp = [&]() {
    auto faces = *skeletonFaces;
    auto skin = subtract_boundary_conditions(get_tets_skin());
    for (auto &bc : *bcSpatialNormalDisplacementVecPtr) {
      faces.merge(intersect(bc.faces, skin));
    }
    return faces;
  };
 
  CHKERR add_l2_field_by_range(hybridSpatialDisp, spaceOrder - 1, 2, 3,
                               get_hybridised_disp());
  CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
                               Range(*contactFaces));
 
  // spatial displacement
  CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
  // material positions
  CHKERR add_h1_field(materialH1Positions, 2, 3);
 
  // Eshelby stress
  // CHKERR add_broken_hdiv_field(eshelbyStress, spaceOrder);
  // CHKERR add_l2_field(materialL2Disp, spaceOrder - 1, 3);
  // CHKERR add_l2_field_by_range(hybridMaterialDisp, spaceOrder - 1, 2, 3,
  //                              Range(*materialSkeletonFaces));
 
  CHKERR mField.build_fields();
 
  MoFEMFunctionReturn(0);
}

addMaterial_Hencky()

MoFEMErrorCode EshelbianCore::addMaterial_Hencky	(	double	E,
		double	nu )

Definition at line 540 of file EshelbianADOL-C.cpp.

                                                                    {
  MoFEMFunctionBegin;
  physicalEquations = boost::make_shared<HMHHencky>(mField, E, nu);
  MoFEMFunctionReturn(0);
}

addMaterial_HMHHStVenantKirchhoff()

MoFEMErrorCode EshelbianCore::addMaterial_HMHHStVenantKirchhoff	(	const int	tape,
		const double	lambda,
		const double	mu,
		const double	sigma_y )

Definition at line 513 of file EshelbianADOL-C.cpp.

                          {
  MoFEMFunctionBegin;
  physicalEquations = boost::make_shared<HMHStVenantKirchhoff>(lambda, mu);
  CHKERR physicalEquations->recordTape(tape, nullptr);
  MoFEMFunctionReturn(0);
}

addMaterial_HMHMooneyRivlin()

MoFEMErrorCode EshelbianCore::addMaterial_HMHMooneyRivlin	(	const int	tape,
		const double	alpha,
		const double	beta,
		const double	lambda,
		const double	sigma_y )

Definition at line 522 of file EshelbianADOL-C.cpp.

                          {
  MoFEMFunctionBegin;
  physicalEquations =
      boost::make_shared<HMHPMooneyRivlinWriggersEq63>(alpha, beta, lambda);
  CHKERR physicalEquations->recordTape(tape, nullptr);
  MoFEMFunctionReturn(0);
}

addMaterial_HMHNeohookean()

MoFEMErrorCode EshelbianCore::addMaterial_HMHNeohookean	(	const int	tape,
		const double	c10,
		const double	K )

Definition at line 531 of file EshelbianADOL-C.cpp.

                                                                        {
  MoFEMFunctionBegin;
  physicalEquations = boost::make_shared<HMHNeohookean>(mField, c10, K);
  CHKERR physicalEquations->recordTape(tape, nullptr);
  MoFEMFunctionReturn(0);
}

addVolumeFiniteElement()

MoFEMErrorCode EshelbianCore::addVolumeFiniteElement

(

const EntityHandle

meshset = 0

)

Definition at line 1530 of file EshelbianPlasticity.cpp.

                                                                {
  MoFEMFunctionBegin;
 
  // set finite element fields
  auto add_field_to_fe = [this](const std::string fe,
                                const std::string field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
    MoFEMFunctionReturn(0);
  };
 
  if (!mField.check_finite_element(elementVolumeName)) {
    CHKERR mField.add_finite_element(elementVolumeName, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(meshset, MBTET,
                                                     elementVolumeName);
 
    CHKERR add_field_to_fe(elementVolumeName, piolaStress);
    CHKERR add_field_to_fe(elementVolumeName, bubbleField);
    if (!noStretch)
      CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
    CHKERR add_field_to_fe(elementVolumeName, rotAxis);
    CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
    CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
    CHKERR add_field_to_fe(elementVolumeName, contactDisp);
    CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
                                                       materialH1Positions);
 
    // build finite elements data structures
    CHKERR mField.build_finite_elements(elementVolumeName);
  }
 
  // if (!mField.check_finite_element(materialVolumeElement)) {
 
  //   Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM - 2);
 
  //   Range front_elements;
  //   for (auto l = 0; l != frontLayers; ++l) {
  //     Range front_elements_layer;
  //     CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM, true,
  //                                              front_elements_layer,
  //                                              moab::Interface::UNION);
  //     front_elements.merge(front_elements_layer);
  //     front_edges.clear();
  //     CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
  //                                              SPACE_DIM - 2, true,
  //                                              front_edges,
  //                                              moab::Interface::UNION);
  //   }
 
  //   CHKERR mField.add_finite_element(materialVolumeElement, MF_ZERO);
  //   CHKERR mField.add_ents_to_finite_element_by_type(front_elements, MBTET,
  //                                                    materialVolumeElement);
  //   // CHKERR add_field_to_fe(materialVolumeElement, eshelbyStress);
  //   // CHKERR add_field_to_fe(materialVolumeElement, materialL2Disp);
  //   CHKERR mField.build_finite_elements(materialVolumeElement);
  // }
 
  MoFEMFunctionReturn(0);
}

calculateCrackArea()

MoFEMErrorCode EshelbianCore::calculateCrackArea

(

boost::shared_ptr< double >

area_ptr

)

Definition at line 7035 of file EshelbianPlasticity.cpp.

                                                                  {
  MoFEMFunctionBegin;
 
  if (!area_ptr) {
    CHK_THROW_MESSAGE(MOFEM_INVALID_DATA, "area_ptr is null");
  }
 
  int success;
  *area_ptr = 0;
  if (mField.get_comm_rank() == 0) {
    MOFEM_LOG("EP", Sev::inform) << "Calculate crack area";
    auto crack_faces = get_range_from_block(mField, "CRACK", SPACE_DIM - 1);
    for (auto f : crack_faces) {
      *area_ptr += mField.getInterface<Tools>()->getTriArea(f);
    }
    success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0, mField.get_comm());
  } else {
    success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0, mField.get_comm());
  }
  if (success != MPI_SUCCESS) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "bcast failed");
  }
  MoFEMFunctionReturn(0);
}

calculateFaceMaterialForce()

MoFEMErrorCode EshelbianCore::calculateFaceMaterialForce	(	const int	tag,
		TS	ts )

Create element to integration faces energies

Definition at line 4299 of file EshelbianPlasticity.cpp.

                                                                             {
  MoFEMFunctionBegin;
 
  constexpr bool debug = false;
 
  auto get_tags_vec = [&](std::vector<std::pair<std::string, int>> names) {
    std::vector<Tag> tags;
    tags.reserve(names.size());
    auto create_and_clean = [&]() {
      MoFEMFunctionBegin;
      for (auto n : names) {
        tags.push_back(Tag());
        auto &tag = tags.back();
        auto &moab = mField.get_moab();
        auto rval = moab.tag_get_handle(n.first.c_str(), tag);
        if (rval == MB_SUCCESS) {
          moab.tag_delete(tag);
        }
        double def_val[] = {0., 0., 0.};
        CHKERR moab.tag_get_handle(n.first.c_str(), n.second, MB_TYPE_DOUBLE,
                                   tag, MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
      }
      MoFEMFunctionReturn(0);
    };
    CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
    return tags;
  };
 
  enum ExhangeTags { MATERIALFORCE, AREAGROWTH, GRIFFITHFORCE, FACEPRESSURE };
 
  auto tags = get_tags_vec({{"MaterialForce", 3},
                            {"AreaGrowth", 3},
                            {"GriffithForce", 1},
                            {"FacePressure", 1}});
 
  auto calculate_material_forces = [&]() {
    MoFEMFunctionBegin;
 
    /**
     * @brief Create element to integration faces energies
     */
    auto get_face_material_force_fe = [&]() {
      using FaceEle = FaceElementForcesAndSourcesCore;
      auto fe_ptr = boost::make_shared<FaceEle>(mField);
      fe_ptr->getRuleHook = [](int, int, int) { return -1; };
      fe_ptr->setRuleHook =
          SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
      EshelbianPlasticity::AddHOOps<2, 2, 3>::add(
          fe_ptr->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
      fe_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
              hybridSpatialDisp, dataAtPts->getGradHybridDispAtPts()));
 
      auto op_loop_domain_side =
          new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
              mField, elementVolumeName, SPACE_DIM, Sev::noisy);
      op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>();
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
          materialH1Positions, frontAdjEdges);
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<3, 3>(piolaStress,
                                               dataAtPts->getApproxPAtPts()));
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateHTensorTensorField<3, 3>(
              bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
      if (noStretch) {
        op_loop_domain_side->getOpPtrVector().push_back(
            physicalEquations->returnOpCalculateStretchFromStress(
                dataAtPts, physicalEquations));
      } 
 
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpFaceSideMaterialForce(dataAtPts));
      fe_ptr->getOpPtrVector().push_back(op_loop_domain_side);
      fe_ptr->getOpPtrVector().push_back(new OpFaceMaterialForce(dataAtPts));
 
      return fe_ptr;
    };
 
    auto integrate_face_material_force_fe = [&](auto &&face_energy_fe) {
      MoFEMFunctionBegin;
      CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
          dM, skeletonElement, face_energy_fe, 0, mField.get_comm_size());
 
      auto face_exchange = CommInterface::createEntitiesPetscVector(
          mField.get_comm(), mField.get_moab(), 2, 3, Sev::inform);
 
      auto print_loc_size = [this](auto v, auto str, auto sev) {
        MoFEMFunctionBegin;
        int size;
        CHKERR VecGetLocalSize(v.second, &size);
        int low, high;
        CHKERR VecGetOwnershipRange(v.second, &low, &high);
        MOFEM_LOG("EPSYNC", sev) << str << " local size " << size << " ( "
                                 << low << " " << high << " ) ";
        MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), sev);
        MoFEMFunctionReturn(0);
      };
      CHKERR print_loc_size(face_exchange, "material face_exchange",
                            Sev::verbose);
 
      CHKERR CommInterface::updateEntitiesPetscVector(
          mField.get_moab(), face_exchange, tags[ExhangeTags::MATERIALFORCE]);
      CHKERR CommInterface::updateEntitiesPetscVector(
          mField.get_moab(), faceExchange, tags[ExhangeTags::FACEPRESSURE]);
 
      // #ifndef NDEBUG
      if (debug) {
        CHKERR save_range(mField.get_moab(),
                          "front_skin_faces_material_force_" +
                              std::to_string(mField.get_comm_rank()) + ".vtk",
                          *skeletonFaces);
      }
      // #endif
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR integrate_face_material_force_fe(get_face_material_force_fe());
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_front_material_force = [&]() {
    MoFEMFunctionBegin;
    FTENSOR_INDEX(SPACE_DIM, I);
 
    auto get_conn = [&](auto e) {
      Range conn;
      CHK_MOAB_THROW(mField.get_moab().get_connectivity(&e, 1, conn, true),
                     "get connectivity");
      return conn;
    };
 
    auto get_conn_range = [&](auto e) {
      Range conn;
      CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
                     "get connectivity");
      return conn;
    };
 
    auto get_adj = [&](auto e, auto dim) {
      Range adj;
      CHK_MOAB_THROW(mField.get_moab().get_adjacencies(&e, 1, dim, true, adj),
                     "get adj");
      return adj;
    };
 
    auto get_adj_range = [&](auto e, auto dim) {
      Range adj;
      CHK_MOAB_THROW(mField.get_moab().get_adjacencies(e, dim, true, adj,
                                                       moab::Interface::UNION),
                     "get adj");
      return adj;
    };
 
    auto get_material_force = [&](auto r, auto th) {
      MatrixDouble material_forces(r.size(), 3, false);
      CHK_MOAB_THROW(
          mField.get_moab().tag_get_data(th, r, material_forces.data().data()),
          "get data");
      return material_forces;
    };
 
    if (mField.get_comm_rank() == 0) {
 
      auto crack_edges = get_adj_range(*crackFaces, 1);
      auto front_nodes = get_conn_range(*frontEdges);
      auto body_edges = get_range_from_block(mField, "EDGES", 1);
      // auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
      // front_block_edges = subtract(front_block_edges, crack_edges);
      // auto front_block_edges_conn = get_conn_range(front_block_edges);
 
      // #ifndef NDEBUG
      Range all_skin_faces;
      Range all_front_faces;
      // #endif
 
      auto calculate_edge_direction = [&](auto e) {
        const EntityHandle *conn;
        int num_nodes;
        CHK_MOAB_THROW(
            mField.get_moab().get_connectivity(e, conn, num_nodes, true),
            "get connectivity");
        std::array<double, 6> coords;
        CHK_MOAB_THROW(
            mField.get_moab().get_coords(conn, num_nodes, coords.data()),
            "get coords");
        FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p0{
            &coords[0], &coords[1], &coords[2]};
        FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p1{
            &coords[3], &coords[4], &coords[5]};
        FTensor::Tensor1<double, SPACE_DIM> t_dir;
        FTENSOR_INDEX(SPACE_DIM, i);
        t_dir(i) = t_p1(i) - t_p0(i);
        return t_dir;
      };
 
      // take bubble tets at node, and then avarage over the edges
      auto calculate_force_through_node = [&]() {
        MoFEMFunctionBegin;
 
        FTENSOR_INDEX(SPACE_DIM, I);
        FTENSOR_INDEX(SPACE_DIM, J);
        FTENSOR_INDEX(SPACE_DIM, K);
        FTENSOR_INDEX(SPACE_DIM, L);
 
        Range body_ents;
        CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                           body_ents);
        auto body_skin = get_skin(mField, body_ents);
        auto body_skin_conn = get_conn_range(body_skin);
 
        // calculate nodal material force
        for (auto n : front_nodes) {
          auto adj_tets = get_adj(n, 3);
          for (int ll = 0; ll < nbJIntegralLevels; ++ll) {
            auto conn = get_conn_range(adj_tets);
            adj_tets = get_adj_range(conn, 3);
          }
          auto skin_faces = get_skin(mField, adj_tets);
          auto material_forces = get_material_force(skin_faces, tags[0]);
 
          // #ifndef NDEBUG
          if (debug) {
            all_skin_faces.merge(skin_faces);
          }
          // #endif
 
          auto calculate_node_material_force = [&]() {
            auto t_face_T =
                getFTensor1FromPtr<SPACE_DIM>(material_forces.data().data());
            FTensor::Tensor1<double, SPACE_DIM> t_node_force{0., 0., 0.};
            for (auto face : skin_faces) {
 
              FTensor::Tensor1<double, SPACE_DIM> t_face_force_tmp{0., 0., 0.};
              t_face_force_tmp(I) = t_face_T(I);
              ++t_face_T;
 
              auto face_tets = intersect(get_adj(face, 3), adj_tets);
 
              if (face_tets.empty()) {
                continue;
              }
 
              if (face_tets.size() != 1) {
                CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                                  "face_tets.size() != 1");
              }
 
              int side_number, sense, offset;
              CHK_MOAB_THROW(mField.get_moab().side_number(face_tets[0], face,
                                                           side_number, sense,
                                                           offset),
                             "moab side number");
              t_face_force_tmp(I) *= sense;
              t_node_force(I) += t_face_force_tmp(I);
            }
 
            return t_node_force;
          };
 
          auto calculate_crack_area_growth_direction = [&](auto n,
                                                           auto &t_node_force) {
            // if skin is on body surface, project the direction on it
            FTensor::Tensor1<double, SPACE_DIM> t_project{0., 0., 0.};
            auto boundary_node = intersect(Range(n, n), body_skin_conn);
            if (boundary_node.size()) {
              auto faces = intersect(get_adj(n, 2), body_skin);
              for (auto f : faces) {
                FTensor::Tensor1<double, 3> t_normal_face;
                CHKERR mField.getInterface<Tools>()->getTriNormal(
                    f, &t_normal_face(0));
                t_project(I) += t_normal_face(I);
              }
              t_project.normalize();
            }
 
            auto adj_faces = intersect(get_adj(n, 2), *crackFaces);
            if (adj_faces.empty()) {
              auto adj_edges =
                  intersect(get_adj(n, 1), unite(*frontEdges, body_edges));
              double l = 0;
              for (auto e : adj_edges) {
                auto t_dir = calculate_edge_direction(e);
                l += t_dir.l2();
              }
              l /= 2;
              FTensor::Tensor1<double, SPACE_DIM> t_area_dir{0., 0., 0.};
              t_area_dir(I) = -t_node_force(I) / t_node_force.l2();
              t_area_dir(I) *= l / 2;
              return t_area_dir;
            }
 
            // calculate surface projection matrix
            FTensor::Tensor2<double, 3, 3> t_Q;
            auto t_kd = FTensor::Kronecker_Delta<double>();
            t_Q(I, J) = t_kd(I, J);
            if (boundary_node.size()) {
              t_Q(I, J) -= t_project(I) * t_project(J);
            }
 
            // calculate direction
            auto front_edges = get_adj(n, 1);
            FTensor::Tensor1<double, 3> t_area_dir{0., 0., 0.};
            for (auto f : adj_faces) {
              int num_nodes;
              const EntityHandle *conn;
              CHKERR mField.get_moab().get_connectivity(f, conn, num_nodes,
                                                        true);
              std::array<double, 9> coords;
              CHKERR mField.get_moab().get_coords(conn, num_nodes,
                                                  coords.data());
              FTensor::Tensor1<double, 3> t_face_normal;
              FTensor::Tensor3<double, 3, 3, 3> t_d_normal;
              CHKERR mField.getInterface<Tools>()->getTriNormal(
                  coords.data(), &t_face_normal(0), &t_d_normal(0, 0, 0));
              auto n_it = std::find(conn, conn + num_nodes, n);
              auto n_index = std::distance(conn, n_it);
 
              FTensor::Tensor2<double, 3, 3> t_face_hessian{
                  t_d_normal(0, n_index, 0), t_d_normal(0, n_index, 1),
                  t_d_normal(0, n_index, 2),
 
                  t_d_normal(1, n_index, 0), t_d_normal(1, n_index, 1),
                  t_d_normal(1, n_index, 2),
 
                  t_d_normal(2, n_index, 0), t_d_normal(2, n_index, 1),
                  t_d_normal(2, n_index, 2)};
 
              FTensor::Tensor2<double, 3, 3> t_projected_hessian;
              t_projected_hessian(I, J) =
                  t_Q(I, K) * (t_face_hessian(K, L) * t_Q(L, J));
              t_face_normal.normalize();
              t_area_dir(K) +=
                  t_face_normal(I) * t_projected_hessian(I, K) / 2.;
            }
 
            return t_area_dir;
          };
 
          auto t_node_force = calculate_node_material_force();
          t_node_force(I) /= griffithEnergy; // scale all by griffith energy
          CHK_MOAB_THROW(
              mField.get_moab().tag_set_data(tags[ExhangeTags::MATERIALFORCE],
                                             &n, 1, &t_node_force(0)),
              "set data");
 
          auto get_area_dir = [&]() {
            auto t_area_dir =
                calculate_crack_area_growth_direction(n, t_node_force);
            if (nbJIntegralLevels) {
              Range seed_n = Range(n, n);
              Range adj_edges;
              for (int l = 0; l < nbJIntegralLevels; ++l) {
                adj_edges.merge(intersect(get_adj_range(seed_n, 1),
                                          unite(*frontEdges, body_edges)));
                seed_n = get_conn_range(adj_edges);
              }
              auto skin_adj_edges = get_skin(mField, adj_edges);
              skin_adj_edges.merge(Range(n, n));
              seed_n = subtract(seed_n, skin_adj_edges);
              for (auto sn : seed_n) {
                auto t_area_dir_sn =
                    calculate_crack_area_growth_direction(sn, t_node_force);
                t_area_dir(I) += t_area_dir_sn(I);
              }
            }
            return t_area_dir;
          };
 
          auto t_area_dir = get_area_dir();
 
          CHK_MOAB_THROW(
              mField.get_moab().tag_set_data(tags[ExhangeTags::AREAGROWTH], &n,
                                             1, &t_area_dir(0)),
              "set data");
          auto griffith = -t_node_force(I) * t_area_dir(I) /
                          (t_area_dir(K) * t_area_dir(K));
          CHK_MOAB_THROW(
              mField.get_moab().tag_set_data(tags[ExhangeTags::GRIFFITHFORCE],
                                             &n, 1, &griffith),
              "set data");
        }
 
        auto ave_node_force = [&](auto th) {
          MoFEMFunctionBegin;
 
          for (auto e : *frontEdges) {
 
            auto conn = get_conn(e);
            auto data = get_material_force(conn, th);
            auto t_node = getFTensor1FromPtr<SPACE_DIM>(data.data().data());
            FTensor::Tensor1<double, SPACE_DIM> t_edge{0., 0., 0.};
            for (auto n : conn) {
              t_edge(I) += t_node(I);
              ++t_node;
            }
            t_edge(I) /= conn.size();
 
            FTensor::Tensor1<double, SPACE_DIM> t_edge_direction =
                calculate_edge_direction(e);
            t_edge_direction.normalize();
 
            FTENSOR_INDEX(SPACE_DIM, I);
            FTENSOR_INDEX(SPACE_DIM, J);
            FTENSOR_INDEX(SPACE_DIM, K);
            FTensor::Tensor1<double, SPACE_DIM> t_cross;
            t_cross(K) =
                FTensor::levi_civita(I, J, K) * t_edge_direction(I) * t_edge(J);
            t_edge(K) = FTensor::levi_civita(I, J, K) * t_edge_direction(J) *
                        t_cross(I);
 
            CHKERR mField.get_moab().tag_set_data(th, &e, 1, &t_edge(0));
          }
          MoFEMFunctionReturn(0);
        };
 
        auto ave_node_griffith_energy = [&](auto th) {
          MoFEMFunctionBegin;
          for (auto e : *frontEdges) {
            FTensor::Tensor1<double, SPACE_DIM> t_edge_force;
            CHKERR mField.get_moab().tag_get_data(
                tags[ExhangeTags::MATERIALFORCE], &e, 1, &t_edge_force(0));
            FTensor::Tensor1<double, SPACE_DIM> t_edge_area_dir;
            CHKERR mField.get_moab().tag_get_data(tags[ExhangeTags::AREAGROWTH],
                                                  &e, 1, &t_edge_area_dir(0));
            double griffith_energy = -t_edge_force(I) * t_edge_area_dir(I) /
                                     (t_edge_area_dir(K) * t_edge_area_dir(K));
            CHKERR mField.get_moab().tag_set_data(th, &e, 1, &griffith_energy);
          }
          MoFEMFunctionReturn(0);
        };
 
        CHKERR ave_node_force(tags[ExhangeTags::MATERIALFORCE]);
        CHKERR ave_node_force(tags[ExhangeTags::AREAGROWTH]);
        CHKERR ave_node_griffith_energy(tags[ExhangeTags::GRIFFITHFORCE]);
 
        MoFEMFunctionReturn(0);
      };
 
      CHKERR calculate_force_through_node();
 
      // calculate face cross
      for (auto e : *frontEdges) {
        auto adj_faces = get_adj(e, 2);
        auto crack_face = intersect(get_adj(e, 2), *crackFaces);
 
        // #ifndef NDEBUG
        if (debug) {
          all_front_faces.merge(adj_faces);
        }
        // #endif
 
        FTensor::Tensor1<double, SPACE_DIM> t_edge_force;
        CHKERR mField.get_moab().tag_get_data(tags[ExhangeTags::MATERIALFORCE],
                                              &e, 1, &t_edge_force(0));
        FTensor::Tensor1<double, SPACE_DIM> t_edge_direction =
            calculate_edge_direction(e);
        t_edge_direction.normalize();
        FTENSOR_INDEX(SPACE_DIM, I);
        FTENSOR_INDEX(SPACE_DIM, J);
        FTENSOR_INDEX(SPACE_DIM, K);
        FTensor::Tensor1<double, SPACE_DIM> t_cross;
        t_cross(K) = FTensor::levi_civita(I, J, K) * t_edge_direction(I) *
                     t_edge_force(J);
 
        for (auto f : adj_faces) {
          FTensor::Tensor1<double, SPACE_DIM> t_normal;
          CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
          t_normal.normalize();
          int side_number, sense, offset;
          CHKERR mField.get_moab().side_number(f, e, side_number, sense,
                                               offset);
          auto dot = -sense * t_cross(I) * t_normal(I);
          CHK_MOAB_THROW(mField.get_moab().tag_set_data(
                             tags[ExhangeTags::GRIFFITHFORCE], &f, 1, &dot),
                         "set data");
        }
      }
 
      // #ifndef NDEBUG
      if (debug) {
        int ts_step;
        CHKERR TSGetStepNumber(ts, &ts_step);
        CHKERR save_range(mField.get_moab(),
                          "front_edges_material_force_" +
                              std::to_string(ts_step) + ".vtk",
                          *frontEdges);
        CHKERR save_range(mField.get_moab(),
                          "front_skin_faces_material_force_" +
                              std::to_string(ts_step) + ".vtk",
                          all_skin_faces);
        CHKERR save_range(mField.get_moab(),
                          "front_faces_material_force_" +
                              std::to_string(ts_step) + ".vtk",
                          all_front_faces);
      }
      // #endif
    }
 
    auto edge_exchange = CommInterface::createEntitiesPetscVector(
        mField.get_comm(), mField.get_moab(), 1, 3, Sev::inform);
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), edge_exchange, tags[ExhangeTags::MATERIALFORCE]);
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), edge_exchange, tags[ExhangeTags::AREAGROWTH]);
    CHKERR CommInterface::updateEntitiesPetscVector(
        mField.get_moab(), edgeExchange, tags[ExhangeTags::GRIFFITHFORCE]);
 
    MoFEMFunctionReturn(0);
  };
 
  auto print_results = [&]() {
    MoFEMFunctionBegin;
 
    auto get_conn_range = [&](auto e) {
      Range conn;
      CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
                     "get connectivity");
      return conn;
    };
 
    auto get_tag_data = [&](auto &ents, auto tag, auto dim) {
      std::vector<double> data(ents.size() * dim);
      CHK_MOAB_THROW(mField.get_moab().tag_get_data(tag, ents, data.data()),
                     "get data");
      return data;
    };
 
    if (mField.get_comm_rank() == 0) {
      auto at_nodes = [&]() {
        MoFEMFunctionBegin;
        auto conn = get_conn_range(*frontEdges);
        auto material_force =
            get_tag_data(conn, tags[ExhangeTags::MATERIALFORCE], 3);
        auto area_growth = get_tag_data(conn, tags[ExhangeTags::AREAGROWTH], 3);
        auto griffith_force =
            get_tag_data(conn, tags[ExhangeTags::GRIFFITHFORCE], 1);
        std::vector<double> coords(conn.size() * 3);
        CHK_MOAB_THROW(mField.get_moab().get_coords(conn, coords.data()),
                       "get coords");
        if (conn.size())
          MOFEM_LOG("EPSELF", Sev::inform) << "Material force at nodes";
        for (size_t i = 0; i < conn.size(); ++i) {
          MOFEM_LOG("EPSELF", Sev::inform)
              << "Node " << conn[i] << " coords "
              << coords[i * 3 + 0] << " " << coords[i * 3 + 1] << " "
              << coords[i * 3 + 2] << " material force "
              << material_force[i * 3 + 0] << " "
              << material_force[i * 3 + 1] << " "
              << material_force[i * 3 + 2] << " area growth "
              << area_growth[i * 3 + 0] << " " << area_growth[i * 3 + 1]
              << " " << area_growth[i * 3 + 2] << " griffith force "
              << griffith_force[i];
        }
        MoFEMFunctionReturn(0);
      };
 
      at_nodes();
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR calculate_material_forces();
  CHKERR calculate_front_material_force();
  CHKERR print_results();
 
  MoFEMFunctionReturn(0);
}

calculateOrientation()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::calculateOrientation	(	const int	tag,
		bool	set_orientation )

Iterate over front edges, get adjacent faces, find maximal face energy. Maximal face energy is stored in the edge. Maximal face energy is magnitude of edge Griffith force.

For each front edge, find maximal face energy and orientation. This is by finding angle between edge material force and maximal face normal

Definition at line 4872 of file EshelbianPlasticity.cpp.

                                                                         {
  MoFEMFunctionBegin;
 
  constexpr bool debug = false;
  constexpr auto sev = Sev::verbose;
 
  Range body_ents;
  CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
  auto body_skin = get_skin(mField, body_ents);
  Range body_skin_edges;
  CHKERR mField.get_moab().get_adjacencies(body_skin, 1, false, body_skin_edges,
                                           moab::Interface::UNION);
  Range boundary_skin_verts;
  CHKERR mField.get_moab().get_connectivity(body_skin_edges,
                                            boundary_skin_verts, true);
 
  auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
  Range geometry_edges_verts;
  CHKERR mField.get_moab().get_connectivity(geometry_edges,
                                            geometry_edges_verts, true);
  Range crack_faces_verts;
  CHKERR mField.get_moab().get_connectivity(*crackFaces, crack_faces_verts,
                                            true);
  Range crack_faces_edges;
  CHKERR mField.get_moab().get_adjacencies(
      *crackFaces, 1, true, crack_faces_edges, moab::Interface::UNION);
  Range crack_faces_tets;
  CHKERR mField.get_moab().get_adjacencies(
      *crackFaces, 3, true, crack_faces_tets, moab::Interface::UNION);
 
  Range front_verts;
  CHKERR mField.get_moab().get_connectivity(*frontEdges, front_verts, true);
  Range front_faces;
  CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2, true, front_faces,
                                           moab::Interface::UNION);
  Range front_verts_edges;
  CHKERR mField.get_moab().get_adjacencies(
      front_verts, 1, true, front_verts_edges, moab::Interface::UNION);
 
  auto get_tags_vec = [&](auto tag_name, int dim) {
    std::vector<Tag> tags(1);
 
    if (dim > 3)
      CHK_THROW_MESSAGE(MOFEM_ATOM_TEST_INVALID, "dim>3");
 
    auto create_and_clean = [&]() {
      MoFEMFunctionBegin;
      auto &moab = mField.get_moab();
      auto rval = moab.tag_get_handle(tag_name, tags[0]);
      if (rval == MB_SUCCESS) {
        moab.tag_delete(tags[0]);
      }
      double def_val[] = {0., 0., 0.};
      CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
                                 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
 
    return tags;
  };
 
  auto get_adj_front = [&](bool subtract_crack) {
    Range adj_front;
    CHKERR mField.get_moab().get_adjacencies(*frontEdges, SPACE_DIM - 1, true,
                                             adj_front, moab::Interface::UNION);
    if (subtract_crack)
      adj_front = subtract(adj_front, *crackFaces);
    return adj_front;
  };
 
  MOFEM_LOG_CHANNEL("SELF");
 
  auto th_front_position = get_tags_vec("FrontPosition", 3);
  auto th_max_face_energy = get_tags_vec("MaxFaceEnergy", 1);
 
  if (mField.get_comm_rank() == 0) {
 
    auto get_crack_adj_tets = [&](auto r) {
      Range crack_faces_conn;
      CHKERR mField.get_moab().get_connectivity(r, crack_faces_conn);
      Range crack_faces_conn_tets;
      CHKERR mField.get_moab().get_adjacencies(crack_faces_conn, SPACE_DIM,
                                               true, crack_faces_conn_tets,
                                               moab::Interface::UNION);
      return crack_faces_conn_tets;
    };
 
    auto get_layers_for_sides = [&](auto &side) {
      std::vector<Range> layers;
      auto get = [&]() {
        MoFEMFunctionBegin;
 
        auto get_adj = [&](auto &r, int dim) {
          Range adj;
          CHKERR mField.get_moab().get_adjacencies(r, dim, true, adj,
                                                   moab::Interface::UNION);
          return adj;
        };
 
        auto get_tets = [&](auto r) { return get_adj(r, SPACE_DIM); };
 
        Range front_nodes;
        CHKERR mField.get_moab().get_connectivity(*frontEdges, front_nodes,
                                                  true);
        Range front_faces = get_adj(front_nodes, 2);
        front_faces = subtract(front_faces, *crackFaces);
        auto front_tets = get_tets(front_nodes);
        auto front_side = intersect(side, front_tets);
        layers.push_back(front_side);
        for (;;) {
          auto adj_faces = get_skin(mField, layers.back());
          adj_faces = intersect(adj_faces, front_faces);
          auto adj_faces_tets = get_tets(adj_faces);
          adj_faces_tets = intersect(adj_faces_tets, front_tets);
          layers.push_back(unite(layers.back(), adj_faces_tets));
          if (layers.back().size() == layers[layers.size() - 2].size()) {
            break;
          }
        }
        MoFEMFunctionReturn(0);
      };
      CHK_THROW_MESSAGE(get(), "get_layers_for_sides");
      return layers;
    };
 
    auto sides_pair = get_two_sides_of_crack_surface(mField, *crackFaces);
    auto layers_top = get_layers_for_sides(sides_pair.first);
    auto layers_bottom = get_layers_for_sides(sides_pair.second);
 
#ifndef NDEBUG
    if (debug) {
      CHKERR save_range(
          mField.get_moab(),
          "crack_tets_" +
              boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
          get_crack_adj_tets(*crackFaces));
      CHKERR save_range(mField.get_moab(), "sides_first.vtk", sides_pair.first);
      CHKERR save_range(mField.get_moab(), "sides_second.vtk",
                        sides_pair.second);
      MOFEM_LOG("EP", sev) << "Nb. layers " << layers_top.size();
      int l = 0;
      for (auto &r : layers_top) {
        MOFEM_LOG("EP", sev) << "Layer " << l << " size " << r.size();
        CHKERR save_range(
            mField.get_moab(),
            "layers_top_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
        ++l;
      }
 
      l = 0;
      for (auto &r : layers_bottom) {
        MOFEM_LOG("EP", sev) << "Layer " << l << " size " << r.size();
        CHKERR save_range(
            mField.get_moab(),
            "layers_bottom_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
        ++l;
      }
    }
#endif
 
    auto get_cross = [&](auto t_dir, auto f) {
      FTensor::Tensor1<double, SPACE_DIM> t_normal;
      CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
      t_normal.normalize();
      FTensor::Tensor1<double, SPACE_DIM> t_cross;
      FTENSOR_INDEX(SPACE_DIM, i);
      FTENSOR_INDEX(SPACE_DIM, j);
      FTENSOR_INDEX(SPACE_DIM, k);
      t_cross(i) = FTensor::levi_civita(i, j, k) * t_normal(j) * t_dir(k);
      return t_cross;
    };
 
    auto get_sense = [&](auto f, auto e) {
      int side, sense, offset;
      CHK_MOAB_THROW(mField.get_moab().side_number(f, e, side, sense, offset),
                     "get sense");
      return std::make_tuple(side, sense, offset);
    };
 
    auto calculate_edge_direction = [&](auto e, auto normalize = true) {
      const EntityHandle *conn;
      int num_nodes;
      CHKERR mField.get_moab().get_connectivity(e, conn, num_nodes, true);
      std::array<double, 6> coords;
      CHKERR mField.get_moab().get_coords(conn, num_nodes, coords.data());
      FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p0{
          &coords[0], &coords[1], &coords[2]};
      FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p1{
          &coords[3], &coords[4], &coords[5]};
      FTensor::Tensor1<double, SPACE_DIM> t_dir;
      FTENSOR_INDEX(SPACE_DIM, i);
      t_dir(i) = t_p1(i) - t_p0(i);
      if (normalize)
        t_dir.normalize();
      return t_dir;
    };
 
    auto evaluate_face_energy_and_set_orientation = [&](auto front_edges,
                                                        auto front_faces,
                                                        auto &sides_pair,
                                                        auto th_position) {
      MoFEMFunctionBegin;
 
      Tag th_face_energy;
      Tag th_material_force;
      switch (energyReleaseSelector) {
      case GRIFFITH_FORCE:
      case GRIFFITH_SKELETON:
        CHKERR mField.get_moab().tag_get_handle("GriffithForce",
                                                th_face_energy);
        CHKERR mField.get_moab().tag_get_handle("MaterialForce",
                                                th_material_force);
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
                "Unknown energy release selector");
      };
 
      /**
       * Iterate over front edges, get adjacent faces, find maximal face energy.
       * Maximal face energy is stored in the edge. Maximal face energy is
       * magnitude of edge Griffith force.
       */
      auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
                                          auto &edge_face_max_energy_map) {
        MoFEMFunctionBegin;
 
        Range body_ents;
        CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
        auto body_skin = get_skin(mField, body_ents);
 
        Range max_faces;
 
        for (auto e : front_edges) {
 
          double griffith_force;
          CHKERR mField.get_moab().tag_get_data(th_face_energy, &e, 1,
                                                &griffith_force);
 
          Range faces;
          CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
          faces = subtract(intersect(faces, front_faces), body_skin);
          std::vector<double> face_energy(faces.size());
          CHKERR mField.get_moab().tag_get_data(th_face_energy, faces,
                                                face_energy.data());
          auto max_energy_it =
              std::max_element(face_energy.begin(), face_energy.end());
          double max_energy =
              max_energy_it != face_energy.end() ? *max_energy_it : 0;
 
          edge_face_max_energy_map[e] =
              std::make_tuple(faces[max_energy_it - face_energy.begin()],
                              griffith_force, static_cast<double>(0));
          MOFEM_LOG("EP", Sev::inform)
              << "Edge " << e << " griffith force " << griffith_force
              << " max face energy " << max_energy << " factor "
              << max_energy / griffith_force;
 
          max_faces.insert(faces[max_energy_it - face_energy.begin()]);
        }
 
#ifndef NDEBUG
        if (debug) {
          CHKERR save_range(
              mField.get_moab(),
              "max_faces_" +
                  boost::lexical_cast<std::string>(mField.get_comm_rank()) +
                  ".vtk",
              max_faces);
        }
#endif
 
        MoFEMFunctionReturn(0);
      };
 
      /**
       * For each front edge, find maximal face energy and orientation. This is
       * by finding angle between edge material force and maximal face normal
       *
       */
      auto calculate_face_orientation = [&](auto &edge_face_max_energy_map) {
        MoFEMFunctionBegin;
 
        auto up_down_face = [&](
 
                                auto &face_angle_map_up,
                                auto &face_angle_map_down
 
                            ) {
          MoFEMFunctionBegin;
 
          for (auto &m : edge_face_max_energy_map) {
            auto e = m.first;
            auto [max_face, energy, opt_angle] = m.second;
 
            Range faces;
            CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
            faces = intersect(faces, front_faces);
            Range adj_tets; // tetrahedrons adjacent to the face
            CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
                                                     false, adj_tets,
                                                     moab::Interface::UNION);
            if (adj_tets.size()) {
 
              Range adj_tets; // tetrahedrons adjacent to the face
              CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
                                                       false, adj_tets,
                                                       moab::Interface::UNION);
              if (adj_tets.size()) {
 
                Range adj_tets_faces;
                // get faces
                CHKERR mField.get_moab().get_adjacencies(
                    adj_tets, SPACE_DIM - 1, false, adj_tets_faces,
                    moab::Interface::UNION);
                adj_tets_faces = intersect(adj_tets_faces, faces);
                FTENSOR_INDEX(SPACE_DIM, i);
 
                // cross product of face normal and edge direction
                auto t_cross_max =
                    get_cross(calculate_edge_direction(e, true), max_face);
                auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
                t_cross_max(i) *= sense_max;
 
                for (auto t : adj_tets) {
                  Range adj_tets_faces;
                  CHKERR mField.get_moab().get_adjacencies(
                      &t, 1, SPACE_DIM - 1, false, adj_tets_faces);
                  adj_tets_faces = intersect(adj_tets_faces, faces);
                  adj_tets_faces =
                      subtract(adj_tets_faces, Range(max_face, max_face));
 
                  if (adj_tets_faces.size() == 1) {
 
                    // cross product of adjacent face normal and edge
                    // direction
                    auto t_cross = get_cross(calculate_edge_direction(e, true),
                                             adj_tets_faces[0]);
                    auto [side, sense, offset] =
                        get_sense(adj_tets_faces[0], e);
                    t_cross(i) *= sense;
                    double dot = t_cross(i) * t_cross_max(i);
                    auto angle = std::acos(dot);
 
                    double face_energy;
                    CHKERR mField.get_moab().tag_get_data(
                        th_face_energy, adj_tets_faces, &face_energy);
 
                    auto [side_face, sense_face, offset_face] =
                        get_sense(t, max_face);
 
                    if (sense_face > 0) {
                      face_angle_map_up[e] = std::make_tuple(face_energy, angle,
                                                             adj_tets_faces[0]);
 
                    } else {
                      face_angle_map_down[e] = std::make_tuple(
                          face_energy, -angle, adj_tets_faces[0]);
                    }
                  }
                }
              }
            }
          }
 
          MoFEMFunctionReturn(0);
        };
 
        auto calc_optimal_angle = [&](
 
                                      auto &face_angle_map_up,
                                      auto &face_angle_map_down
 
                                  ) {
          MoFEMFunctionBegin;
 
          for (auto &m : edge_face_max_energy_map) {
            auto e = m.first;
            auto &[max_face, e0, a0] = m.second;
 
            if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
 
              if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
                  face_angle_map_down.find(e) == face_angle_map_down.end()) {
                // Do nothing
              } else {
 
                switch (energyReleaseSelector) {
                case GRIFFITH_FORCE:
                case GRIFFITH_SKELETON: {
 
                  Tag th_material_force;
                  CHKERR mField.get_moab().tag_get_handle("MaterialForce",
                                                          th_material_force);
                  FTensor::Tensor1<double, SPACE_DIM> t_material_force;
                  CHKERR mField.get_moab().tag_get_data(
                      th_material_force, &e, 1, &t_material_force(0));
                  auto material_force_magnitude = t_material_force.l2();
                  if (material_force_magnitude <
                      std::numeric_limits<double>::epsilon()) {
                    a0 = 0;
 
                  } else {
 
                    auto t_edge_dir = calculate_edge_direction(e, true);
                    auto t_cross_max = get_cross(t_edge_dir, max_face);
                    auto [side, sense, offset] = get_sense(max_face, e);
                    t_cross_max(sense) *= sense;
 
                    FTENSOR_INDEX(SPACE_DIM, I);
                    FTENSOR_INDEX(SPACE_DIM, J);
                    FTENSOR_INDEX(SPACE_DIM, K);
 
                    t_material_force.normalize();
                    t_cross_max.normalize();
                    FTensor::Tensor1<double, SPACE_DIM> t_cross;
                    t_cross(I) = FTensor::levi_civita(I, J, K) *
                                 t_material_force(J) * t_cross_max(K);
                    a0 = -std::asin(t_cross(I) * t_edge_dir(I));
 
                    MOFEM_LOG("EP", sev)
                        << "Optimal angle " << a0 << " energy " << e0;
                  }
                  break;
                }
                default: {
 
                  SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
                          "Unknown energy release selector");
                }
                }
              }
            }
          }
 
          MoFEMFunctionReturn(0);
        };
 
        std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
            face_angle_map_up;
        std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
            face_angle_map_down;
        CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
        CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
 
#ifndef NDEBUG
        if (debug) {
          auto th_angle = get_tags_vec("Angle", 1);
          Range up;
          for (auto &m : face_angle_map_up) {
            auto [e, a, face] = m.second;
            up.insert(face);
            CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
          }
          Range down;
          for (auto &m : face_angle_map_down) {
            auto [e, a, face] = m.second;
            down.insert(face);
            CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
          }
 
          Range max_energy_faces;
          for (auto &m : edge_face_max_energy_map) {
            auto [face, e, angle] = m.second;
            max_energy_faces.insert(face);
            CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1,
                                                  &angle);
          }
          if (mField.get_comm_rank() == 0) {
            CHKERR save_range(mField.get_moab(), "up_faces.vtk", up);
            CHKERR save_range(mField.get_moab(), "down_faces.vtk", down);
            CHKERR save_range(mField.get_moab(), "max_energy_faces.vtk",
                              max_energy_faces);
          }
        }
#endif // NDEBUG
 
        MoFEMFunctionReturn(0);
      };
 
      auto get_conn = [&](auto e) {
        Range conn;
        CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
                       "get conn");
        return conn;
      };
 
      auto get_adj = [&](auto e, auto dim) {
        Range adj;
        CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
                           e, dim, false, adj, moab::Interface::UNION),
                       "get adj");
        return adj;
      };
 
      auto get_coords = [&](auto v) {
        FTensor::Tensor1<double, 3> t_coords;
        CHK_MOAB_THROW(mField.get_moab().get_coords(v, &t_coords(0)),
                       "get coords");
        return t_coords;
      };
 
      // calulate normal of the max energy face
      auto get_rotated_normal = [&](auto e, auto f, auto angle) {
        FTENSOR_INDEX(SPACE_DIM, i);
        FTENSOR_INDEX(SPACE_DIM, j);
        auto t_edge_dir = calculate_edge_direction(e, true);
        auto [side, sense, offset] = get_sense(f, e);
        t_edge_dir(i) *= sense;
        t_edge_dir.normalize();
        t_edge_dir(i) *= angle;
        auto t_R = LieGroups::SO3::exp(t_edge_dir, angle);
        FTensor::Tensor1<double, SPACE_DIM> t_normal;
        mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
        FTensor::Tensor1<double, SPACE_DIM> t_rotated_normal;
        t_rotated_normal(i) = t_R(i, j) * t_normal(j);
        return std::make_tuple(t_normal, t_rotated_normal);
      };
 
      auto set_coord = [&](auto v, auto &adj_vertex_tets_verts, auto &coords,
                           auto &t_move, auto gamma) {
        auto index = adj_vertex_tets_verts.index(v);
        if (index >= 0) {
          for (auto ii : {0, 1, 2}) {
            coords[3 * index + ii] += gamma * t_move(ii);
          }
          return true;
        }
        return false;
      };
 
      auto tets_quality = [&](auto quality, auto &adj_vertex_tets_verts,
                              auto &adj_vertex_tets, auto &coords) {
        for (auto t : adj_vertex_tets) {
          const EntityHandle *conn;
          int num_nodes;
          CHKERR mField.get_moab().get_connectivity(t, conn, num_nodes, true);
          std::array<double, 12> tet_coords;
          for (auto n = 0; n != 4; ++n) {
            auto index = adj_vertex_tets_verts.index(conn[n]);
            if (index < 0) {
              CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong index");
            }
            for (auto ii = 0; ii != 3; ++ii) {
              tet_coords[3 * n + ii] = coords[3 * index + ii];
            }
          }
          double q = Tools::volumeLengthQuality(tet_coords.data());
          if (!std::isnormal(q))
            q = -2;
          quality = std::min(quality, q);
        };
 
        return quality;
      };
 
      auto calculate_free_face_node_displacement =
          [&](auto &edge_face_max_energy_map) {
            // get edges adjacent to vertex along which nodes are moving
            auto get_vertex_edges = [&](auto vertex) {
              Range vertex_edges; // edges adjacent to vertex
 
              auto impl = [&]() {
                MoFEMFunctionBegin;
                CHKERR mField.get_moab().get_adjacencies(vertex, 1, false,
                                                         vertex_edges);
                vertex_edges = subtract(vertex_edges, front_verts_edges);
 
                if (boundary_skin_verts.size() &&
                    boundary_skin_verts.find(vertex[0]) !=
                        boundary_skin_verts.end()) {
                  MOFEM_LOG("EP", sev) << "Boundary vertex";
                  vertex_edges = intersect(vertex_edges, body_skin_edges);
                }
                if (geometry_edges_verts.size() &&
                    geometry_edges_verts.find(vertex[0]) !=
                        geometry_edges_verts.end()) {
                  MOFEM_LOG("EP", sev) << "Geometry edge vertex";
                  vertex_edges = intersect(vertex_edges, geometry_edges);
                }
                if (crack_faces_verts.size() &&
                    crack_faces_verts.find(vertex[0]) !=
                        crack_faces_verts.end()) {
                  MOFEM_LOG("EP", sev) << "Crack face vertex";
                  vertex_edges = intersect(vertex_edges, crack_faces_edges);
                }
                MoFEMFunctionReturn(0);
              };
 
              CHK_THROW_MESSAGE(impl(), "get_vertex_edges");
 
              return vertex_edges;
            };
 
            // vector of rotated faces, edge along node is moved, moved edge,
            // moved displacement, quality, cardinality, gamma
            using Bundle = std::vector<
 
                std::tuple<EntityHandle, EntityHandle, EntityHandle,
                           FTensor::Tensor1<double, 3>, double, double, double>
 
                >;
            std::map<EntityHandle, Bundle> edge_bundle_map;
 
            for (auto &m : edge_face_max_energy_map) {
 
              auto edge = m.first;
              auto &[max_face, energy, opt_angle] = m.second;
 
              // calculate rotation of max energy face
              auto [t_normal, t_rotated_normal] =
                  get_rotated_normal(edge, max_face, opt_angle);
 
              auto front_vertex = get_conn(Range(m.first, m.first));
              auto adj_tets = get_adj(Range(max_face, max_face), 3);
              auto adj_tets_faces = get_adj(adj_tets, 2);
              auto adj_front_faces = subtract(
                  intersect(get_adj(Range(edge, edge), 2), adj_tets_faces),
                  *crackFaces);
              if (adj_front_faces.size() > 3)
                CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                                  "adj_front_faces.size()>3");
 
              FTensor::Tensor1<double, SPACE_DIM> t_material_force;
              CHKERR mField.get_moab().tag_get_data(th_material_force, &edge, 1,
                                                    &t_material_force(0));
              std::vector<double> griffith_energy(adj_front_faces.size());
              CHKERR mField.get_moab().tag_get_data(
                  th_face_energy, adj_front_faces, griffith_energy.data());
 
 
              auto set_edge_bundle = [&](auto min_gamma) {
                for (auto rotated_f : adj_front_faces) {
 
                  double rotated_face_energy =
                      griffith_energy[adj_front_faces.index(rotated_f)];
 
                  auto vertex = subtract(get_conn(Range(rotated_f, rotated_f)),
                                         front_vertex);
                  if (vertex.size() != 1) {
                    CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                                      "Wrong number of vertex to move");
                  }
                  auto front_vertex_edges_vertex = get_conn(
                      intersect(get_adj(front_vertex, 1), crack_faces_edges));
                  vertex = subtract(
                      vertex, front_vertex_edges_vertex); // vertex free to move
                  if (vertex.empty()) {
                    continue;
                  }
 
                  auto face_cardinality = [&](auto f, auto &seen_front_edges) {
                    auto whole_front =
                        unite(*frontEdges,
                              subtract(body_skin_edges, crack_faces_edges));
                    auto faces = Range(f, f);
                    int c = 0;
                    for (; c < 10; ++c) {
                      auto front_edges =
                          subtract(get_adj(faces, 1), seen_front_edges);
                      if (front_edges.size() == 0) {
                        return 0;
                      }
                      auto front_connected_edges =
                          intersect(front_edges, whole_front);
                      if (front_connected_edges.size()) {
                        seen_front_edges.merge(front_connected_edges);
                        return c;
                      }
                      faces.merge(get_adj(front_edges, 2));
                      ++c;
                    }
                    return c;
                  };
 
                  Range seen_edges = Range(edge, edge);
                  double rotated_face_cardinality = face_cardinality(
                      rotated_f,
                      seen_edges); // add cardinality of max energy
                                   // face to rotated face cardinality
                  // rotated_face_cardinality +=
                  //     face_cardinality(max_face, seen_edges);
                  rotated_face_cardinality = std::max(rotated_face_cardinality,
                                                      1.); // at least one edge
 
                  auto t_vertex_coords = get_coords(vertex);
                  auto vertex_edges = get_vertex_edges(vertex);
 
                  EntityHandle f0 = front_vertex[0];
                  EntityHandle f1 = front_vertex[1];
                  FTensor::Tensor1<double, 3> t_v_e0, t_v_e1;
                  CHKERR mField.get_moab().get_coords(&f0, 1, &t_v_e0(0));
                  CHKERR mField.get_moab().get_coords(&f1, 1, &t_v_e1(0));
 
                  FTENSOR_INDEX(SPACE_DIM, i);
                  for (auto e_used_to_move_detection : vertex_edges) {
                    auto edge_conn = get_conn(Range(e_used_to_move_detection,
                                                    e_used_to_move_detection));
                    edge_conn = subtract(edge_conn, vertex);
                    // Find displacement of the edge such that dot porduct with
                    // normal is zero.
                    //
                    // { (t_v0 - t_vertex_coords) + gamma * (t_v3 -
                    // t_vertex_coords) } * n = 0
                    // where t_v0 is the edge vertex, t_v3 is the edge end
                    // point, n is the rotated normal of the face gamma is the
                    // factor by which the edge is moved
                    FTensor::Tensor1<double, 3> t_v0;
                    t_v0(i) = (t_v_e0(i) + t_v_e1(i)) / 2;
                    FTensor::Tensor1<double, 3> t_v3;
                    CHKERR mField.get_moab().get_coords(edge_conn, &t_v3(0));
                    auto a =
                        (t_v0(i) - t_vertex_coords(i)) * t_rotated_normal(i);
                    auto b =
                        (t_v3(i) - t_vertex_coords(i)) * t_rotated_normal(i);
                    auto gamma = a / b;
 
                    constexpr double eps =
                        std::numeric_limits<double>::epsilon();
                    if (std::isnormal(gamma) && gamma < 1.0 - eps &&
                        gamma > -0.1) {
                      FTensor::Tensor1<double, 3> t_move;
                      t_move(i) = gamma * (t_v3(i) - t_vertex_coords(i));
 
                      auto check_rotated_face_directoon = [&]() {
                        FTensor::Tensor1<double, SPACE_DIM> t_delta;
                        t_delta(i) = t_vertex_coords(i) + t_move(i) - t_v0(i);
                        t_delta.normalize();
                        auto dot =
                            (t_material_force(i) / t_material_force.l2()) *
                            t_delta(i);
                        return -dot > 0 ? true : false;
                      };
 
                      if (check_rotated_face_directoon()) {
 
                        MOFEM_LOG("EP", Sev::inform)
                            << "Crack edge " << edge << " moved face "
                            << rotated_f
                            << " edge: " << e_used_to_move_detection
                            << " face direction/energy " << rotated_face_energy
                            << " face cardinality " << rotated_face_cardinality
                            << " gamma: " << gamma;
 
                        auto &bundle = edge_bundle_map[edge];
                        bundle.emplace_back(rotated_f, e_used_to_move_detection,
                                            vertex[0], t_move, 1,
                                            rotated_face_cardinality, gamma);
                      }
                    }
                  }
                }
              };
 
              set_edge_bundle(std::numeric_limits<double>::epsilon());
              if (edge_bundle_map[edge].empty()) {
                set_edge_bundle(-1.);
              }
            }
 
            return edge_bundle_map;
          };
 
      auto get_sort_by_energy = [&](auto &edge_face_max_energy_map) {
        std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
            sort_by_energy;
 
        for (auto &m : edge_face_max_energy_map) {
          auto e = m.first;
          auto &[max_face, energy, opt_angle] = m.second;
            sort_by_energy[energy] = std::make_tuple(e, max_face, opt_angle);
        }
 
        return sort_by_energy;
      };
 
      auto set_tag = [&](auto &&adj_edges_map, auto &&sort_by_energy) {
        MoFEMFunctionBegin;
 
        Tag th_face_pressure;
        CHK_MOAB_THROW(
            mField.get_moab().tag_get_handle("FacePressure", th_face_pressure),
            "get tag");
        auto get_face_pressure = [&](auto face) {
          double pressure;
          CHK_MOAB_THROW(mField.get_moab().tag_get_data(th_face_pressure, &face,
                                                        1, &pressure),
                         "get rag data");
          return pressure;
        };
 
        MOFEM_LOG("EPSELF", Sev::inform)
            << "Number of edges to check " << sort_by_energy.size();
 
        enum face_energy { POSITIVE, NEGATIVE };
        constexpr bool skip_negative = true;
 
        for (auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
 
          // iterate edges wih maximal energy, and make them seed. Such edges,
          // will most likely will have also smallest node displacement
          for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
               ++it) {
 
            auto energy = it->first;
            auto [max_edge, max_face, opt_angle] = it->second;
 
            auto face_pressure = get_face_pressure(max_face);
            if (skip_negative) {
              if (fe == face_energy::POSITIVE) {
                if (face_pressure < 0) {
                  MOFEM_LOG("EPSELF", Sev::inform)
                      << "Skip negative face " << max_face << " with energy "
                      << energy << " and pressure " << face_pressure;
                  continue;
                }
              }
            }
 
            MOFEM_LOG("EPSELF", Sev::inform)
                << "Check face " << max_face << " edge " << max_edge
                << " energy " << energy << " optimal angle " << opt_angle
                << " face pressure " << face_pressure;
 
            auto jt = adj_edges_map.find(max_edge);
            if (jt == adj_edges_map.end()) {
              MOFEM_LOG("EPSELF", Sev::warning)
                  << "Edge " << max_edge << " not found in adj_edges_map";
              continue;
            }
            auto &bundle = jt->second;
 
            auto find_max_in_bundle_impl = [&](auto edge, auto &bundle,
                                               auto gamma) {
              FTENSOR_INDEX(SPACE_DIM, i);
 
              EntityHandle vertex_max = 0;
              EntityHandle face_max = 0;
              EntityHandle move_edge_max = 0;
              double max_quality = -2;
              double max_quality_evaluated = -2;
              double min_cardinality = std::numeric_limits<double>::max();
 
              FTensor::Tensor1<double, SPACE_DIM> t_move_last{0., 0., 0.};
 
              for (auto &b : bundle) {
                auto &[face, move_edge, vertex, t_move, quality, cardinality,
                       edge_gamma] = b;
 
                auto adj_vertex_tets = get_adj(Range(vertex, vertex), 3);
                auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
                std::vector<double> coords(3 * adj_vertex_tets_verts.size());
                CHK_MOAB_THROW(mField.get_moab().get_coords(
                                   adj_vertex_tets_verts, coords.data()),
                               "get coords");
 
                set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
                quality = tets_quality(quality, adj_vertex_tets_verts,
                                       adj_vertex_tets, coords);
 
                auto eval_quality = [](auto q, auto c, auto edge_gamma) {
                  if (q < 0) {
                    return q;
                  } else {
                    return ((edge_gamma < 0) ? (q / 2) : q) / pow(c, 2);
                  }
                };
 
                if (eval_quality(quality, cardinality, edge_gamma) >=
                    max_quality_evaluated) {
                  max_quality = quality;
                  min_cardinality = cardinality;
                  vertex_max = vertex;
                  face_max = face;
                  move_edge_max = move_edge;
                  t_move_last(i) = t_move(i);
                  max_quality_evaluated =
                      eval_quality(max_quality, min_cardinality, edge_gamma);
                }
              }
 
              return std::make_tuple(vertex_max, face_max, t_move_last,
                                     max_quality, min_cardinality);
            };
 
            auto find_max_in_bundle = [&](auto edge, auto &bundle) {
              auto b_org_bundle = bundle;
              auto r = find_max_in_bundle_impl(edge, bundle, 1.);
              auto &[vertex_max, face_max, t_move_last, max_quality,
                     cardinality] = r;
              if (max_quality < 0) {
                for (double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
                  bundle = b_org_bundle;
                  r = find_max_in_bundle_impl(edge, bundle, gamma);
                  auto &[vertex_max, face_max, t_move_last, max_quality,
                         cardinality] = r;
                  MOFEM_LOG("EPSELF", Sev::warning)
                      << "Back tracking: gamma " << gamma << " edge " << edge
                      << " quality " << max_quality << " cardinality "
                      << cardinality;
                  if (max_quality > 0.01) {
                    FTENSOR_INDEX(SPACE_DIM, I);
                    t_move_last(I) *= gamma;
                    return r;
                  }
                }
                FTENSOR_INDEX(SPACE_DIM, I);
                t_move_last(I) = 0;
              }
              return r;
            };
 
            // set tags with displacement of node and face energy
            auto set_tag_to_vertex_and_face = [&](auto &&r, auto &quality) {
              MoFEMFunctionBegin;
              auto &[v, f, t_move, q, cardinality] = r;
 
              if ((q > 0 && std::isnormal(q)) && energy > 0) {
 
                MOFEM_LOG("EPSELF", Sev::inform)
                    << "Set tag: vertex " << v << " face " << f << " "
                    << max_edge << " move " << t_move << " energy " << energy
                    << " quality " << q << " cardinality " << cardinality;
                CHKERR mField.get_moab().tag_set_data(th_position[0], &v, 1,
                                                      &t_move(0));
                CHKERR mField.get_moab().tag_set_data(th_max_face_energy[0], &f,
                                                      1, &energy);
              }
 
              quality = q;
              MoFEMFunctionReturn(0);
            };
 
            double quality = -2;
            CHKERR set_tag_to_vertex_and_face(
 
                find_max_in_bundle(max_edge, bundle),
 
                quality
 
            );
 
            if (quality > 0 && std::isnormal(quality) && energy > 0) {
              MOFEM_LOG("EPSELF", Sev::inform)
                  << "Crack face set with quality: " << quality;
              MoFEMFunctionReturnHot(0);
            }
          }
 
          if (!skip_negative)
            break;
        }
 
        MoFEMFunctionReturn(0);
      };
 
      // map: {edge, {face, energy, optimal_angle}}
      MOFEM_LOG("EP", sev) << "Calculate orientation";
      std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
          edge_face_max_energy_map;
      CHKERR find_maximal_face_energy(front_edges, front_faces,
                                      edge_face_max_energy_map);
      CHKERR calculate_face_orientation(edge_face_max_energy_map);
 
      MOFEM_LOG("EP", sev) << "Calculate node positions";
      CHKERR set_tag(
 
          calculate_free_face_node_displacement(edge_face_max_energy_map),
          get_sort_by_energy(edge_face_max_energy_map)
 
      );
 
      MoFEMFunctionReturn(0);
    };
 
    MOFEM_LOG("EP", sev) << "Front edges " << frontEdges->size();
    CHKERR evaluate_face_energy_and_set_orientation(
        *frontEdges, get_adj_front(false), sides_pair, th_front_position);
  }
 
  // exchange positions and energies from processor zero to all other
  CHKERR VecZeroEntries(vertexExchange.second);
  CHKERR VecGhostUpdateBegin(vertexExchange.second, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(vertexExchange.second, INSERT_VALUES,
                           SCATTER_FORWARD);
  CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
      mField.get_moab(), vertexExchange, th_front_position[0]);
  CHKERR VecZeroEntries(faceExchange.second);
  CHKERR VecGhostUpdateBegin(faceExchange.second, INSERT_VALUES,
                             SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(faceExchange.second, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
      mField.get_moab(), faceExchange, th_max_face_energy[0]);
 
  auto get_max_moved_faces = [&]() {
    Range max_moved_faces;
    auto adj_front = get_adj_front(false);
    std::vector<double> face_energy(adj_front.size());
    CHKERR mField.get_moab().tag_get_data(th_max_face_energy[0], adj_front,
                                          face_energy.data());
    for (int i = 0; i != adj_front.size(); ++i) {
      if (face_energy[i] > std::numeric_limits<double>::epsilon()) {
        max_moved_faces.insert(adj_front[i]);
      }
    }
 
    return boost::make_shared<Range>(max_moved_faces);
  };
 
  // move all faces with energy larger than 0
  maxMovedFaces = get_max_moved_faces();
  MOFEM_LOG("EP", sev) << "Number of of moved faces: " << maxMovedFaces->size();
 
#ifndef NDEBUG
  if (debug) {
    CHKERR save_range(
        mField.get_moab(),
        "max_moved_faces_" +
            boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
        *maxMovedFaces);
  }
#endif
 
  MoFEMFunctionReturn(0);
}

createCrackSurfaceMeshset()

MoFEMErrorCode EshelbianCore::createCrackSurfaceMeshset

(

)

Definition at line 6698 of file EshelbianPlasticity.cpp.

                                                        {
  MoFEMFunctionBegin;
  auto meshset_mng = mField.getInterface<MeshsetsManager>();
  while (meshset_mng->checkMeshset(addCrackMeshsetId, BLOCKSET))
    ++addCrackMeshsetId;
  MOFEM_LOG("EP", Sev::inform)
      << "Crack added surface meshset " << addCrackMeshsetId;
  CHKERR meshset_mng->addMeshset(BLOCKSET, addCrackMeshsetId, "CRACK_COMPUTED");
  MoFEMFunctionReturn(0);
};

createExchangeVectors()

MoFEMErrorCode EshelbianCore::createExchangeVectors

(

Sev

sev

)

Definition at line 7003 of file EshelbianPlasticity.cpp.

                                                           {
  MoFEMFunctionBegin;
 
  auto print_loc_size = [this](auto v, auto str, auto sev) {
    MoFEMFunctionBegin;
    int size;
    CHKERR VecGetLocalSize(v.second, &size);
    int low, high;
    CHKERR VecGetOwnershipRange(v.second, &low, &high);
    MOFEM_LOG("EPSYNC", sev) << str << " local size " << size << " ( " << low
                             << " " << high << " ) ";
    MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), sev);
    MoFEMFunctionReturn(0);
  };
 
  volumeExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 3, 1, sev);
  CHKERR print_loc_size(volumeExchange, "volumeExchange", sev);
  faceExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 2, 1, Sev::inform);
  CHKERR print_loc_size(faceExchange, "faceExchange", sev);
  edgeExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 1, 1, Sev::inform);
  CHKERR print_loc_size(edgeExchange, "edgeExchange", sev);
  vertexExchange = CommInterface::createEntitiesPetscVector(
      mField.get_comm(), mField.get_moab(), 0, 3, Sev::inform);
  CHKERR print_loc_size(vertexExchange, "vertexExchange", sev);
 
  MoFEMFunctionReturn(0);
}

d_f_linear()

static double EshelbianPlasticity::EshelbianCore::d_f_linear ( const double v )

inlinestatic

Definition at line 160 of file EshelbianPlasticity.hpp.

d_f_log()

static double EshelbianPlasticity::EshelbianCore::d_f_log ( const double v )

inlinestatic

Definition at line 141 of file EshelbianPlasticity.hpp.

d_f_log_e()

static double EshelbianPlasticity::EshelbianCore::d_f_log_e ( const double v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 96 of file EshelbianPlasticity.hpp.

d_f_log_e_quadratic()

static double EshelbianPlasticity::EshelbianCore::d_f_log_e_quadratic ( const double v )

inlinestatic

Definition at line 65 of file EshelbianPlasticity.hpp.

                           { GRIFFITH_FORCE, GRIFFITH_SKELETON };
 
using MatrixPtr = boost::shared_ptr<MatrixDouble>;
using VectorPtr = boost::shared_ptr<VectorDouble>;
 
using EntData = EntitiesFieldData::EntData;
using UserDataOperator = ForcesAndSourcesCore::UserDataOperator;
using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator;

dd_f_linear()

static double EshelbianPlasticity::EshelbianCore::dd_f_linear ( const double v )

inlinestatic

Definition at line 161 of file EshelbianPlasticity.hpp.

dd_f_log()

static double EshelbianPlasticity::EshelbianCore::dd_f_log ( const double v )

inlinestatic

Definition at line 144 of file EshelbianPlasticity.hpp.

dd_f_log_e()

static double EshelbianPlasticity::EshelbianCore::dd_f_log_e ( const double v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 106 of file EshelbianPlasticity.hpp.

dd_f_log_e_quadratic()

static double EshelbianPlasticity::EshelbianCore::dd_f_log_e_quadratic ( const double v )

inlinestatic

Definition at line 75 of file EshelbianPlasticity.hpp.

f_linear()

static double EshelbianPlasticity::EshelbianCore::f_linear ( const double v )

inlinestatic

Definition at line 159 of file EshelbianPlasticity.hpp.

f_log()

static double EshelbianPlasticity::EshelbianCore::f_log ( const double v )

inlinestatic

Definition at line 138 of file EshelbianPlasticity.hpp.

f_log_e()

static double EshelbianPlasticity::EshelbianCore::f_log_e ( const double v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 83 of file EshelbianPlasticity.hpp.

    : public boost::enable_shared_from_this<DataAtIntegrationPts> {
 
  MatrixDouble approxPAtPts;
  MatrixDouble approxSigmaAtPts;
  MatrixDouble divPAtPts;
  MatrixDouble divSigmaAtPts;
  MatrixDouble wL2AtPts;
  MatrixDouble wL2DotAtPts;
  MatrixDouble wL2DotDotAtPts;
  MatrixDouble logStretchTensorAtPts;
  MatrixDouble stretchTensorAtPts;
  VectorDouble jacobianAtPts;
  MatrixDouble tractionAtPts;

f_log_e_quadratic()

static double EshelbianPlasticity::EshelbianCore::f_log_e_quadratic ( const double v )

inlinestatic

Definition at line 55 of file EshelbianPlasticity.hpp.

                      { SYMMETRIC, NOT_SYMMETRIC, ANIT_SYMMETRIC };
enum RotSelector { SMALL_ROT, MODERATE_ROT, LARGE_ROT, NO_H1_CONFIGURATION };
enum StretchSelector {
  LINEAR,
  LOG /*linar extenion*/,
  LOG_QUADRATIC /*quadratic extension*/,

getBc()

template<typename BC >

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getBc ( boost::shared_ptr< BC > & bc_vec_ptr, const std::string block_name, const int nb_attributes )

inline

Definition at line 244 of file EshelbianPlasticity.hpp.

                                {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matD);
  }
 
  inline MatrixPtr getMatInvDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matInvD);
  }
 
  inline MatrixPtr getMatAxiatorDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matAxiatorD);
  }
 
  inline MatrixPtr getMatDeviatorDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matDeviatorD);
  }
 
  inline MatrixPtr getSmallWH1AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wH1AtPts);
  }
 
  inline MatrixPtr getLargeXH1AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &XH1AtPts);
  }
 
  inline MatrixPtr getContactL2AtPts() {

getExternalStrain()

MoFEMErrorCode EshelbianCore::getExternalStrain

(

)

Definition at line 6956 of file EshelbianPlasticity.cpp.

                                                {
  MoFEMFunctionBegin;
 
  auto getExternalStrain = [&](boost::shared_ptr<ExternalStrainVec> &ext_strain_vec_ptr,
                             const std::string block_name,
                             const int nb_attributes) {
    MoFEMFunctionBegin;
    for (auto it : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
             std::regex((boost::format("%s(.*)") % block_name).str()))
    ) {
      std::vector<double> block_attributes;
      CHKERR it->getAttributes(block_attributes);
      if (block_attributes.size() < nb_attributes) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "In block %s expected %d attributes, but given %ld",
                 it->getName().c_str(), nb_attributes, block_attributes.size());
      }
 
      auto get_block_ents = [&]() {
        Range ents;
        CHKERR mField.get_moab().get_entities_by_handle(it->meshset, ents,
                                                        true);
        return ents;
      };
      auto Ents = get_block_ents();
      ext_strain_vec_ptr->emplace_back(it->getName(), block_attributes,
                               get_block_ents());
    }
    MoFEMFunctionReturn(0);
  };
 
  externalStrainVecPtr = boost::make_shared<ExternalStrainVec>();
  CHKERR getExternalStrain(externalStrainVecPtr, "EXTERNALSTRAIN", 2);
 
  auto ts_pre_stretch =
      boost::make_shared<DynamicRelaxationTimeScale>("externalstrain_history.txt");
  for (auto &ext_strain_block: *externalStrainVecPtr) {
    MOFEM_LOG("EP", Sev::noisy)
        << "Add time scaling external strain: " << ext_strain_block.blockName;
    timeScaleMap[ext_strain_block.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_pre_stretch, "externalstrain_history", ".txt", ext_strain_block.blockName);
  }
 
  MoFEMFunctionReturn(0);
}

getOptions()

MoFEMErrorCode EshelbianCore::getOptions

(

)

Definition at line 919 of file EshelbianPlasticity.cpp.

                                         {
  MoFEMFunctionBegin;
  const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
  const char *list_symm[] = {"symm", "not_symm"};
  const char *list_release[] = {"griffith_force", "griffith_skelton"};
  const char *list_stretches[] = {"linear", "log", "log_quadratic"};
  PetscInt choice_rot = EshelbianCore::rotSelector;
  PetscInt choice_grad = EshelbianCore::gradApproximator;
  PetscInt choice_symm = EshelbianCore::symmetrySelector;
  PetscInt choice_release = EshelbianCore::energyReleaseSelector;
  PetscInt choice_stretch = StretchSelector::LOG;
  char analytical_expr_file_name[255] = "analytical_expr.py";
 
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
                           "none");
  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
                         spaceOrder, &spaceOrder, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
                         spaceH1Order, &spaceH1Order, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
                         "", materialOrder, &materialOrder, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
                            &alphaU, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
                            &alphaW, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-viscosity_alpha_omega", "rot viscosity", "",
                            alphaOmega, &alphaOmega, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
                            &alphaRho, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
                           LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
                           PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-grad", "gradient of defamation approximate", "",
                           list_rots, NO_H1_CONFIGURATION + 1,
                           list_rots[choice_grad], &choice_grad, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-symm", "symmetric variant", "", list_symm, 2,
                           list_symm[choice_symm], &choice_symm, PETSC_NULLPTR);
 
  CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
                            &EshelbianCore::exponentBase, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-stretches", "stretches", "", list_stretches,
                           StretchSelector::STRETCH_SELECTOR_LAST,
                           list_stretches[choice_stretch], &choice_stretch,
                           PETSC_NULLPTR);
 
  CHKERR PetscOptionsBool("-no_stretch", "do not solve for stretch", "",
                          noStretch, &noStretch, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-set_singularity", "set singularity", "",
                          setSingularity, &setSingularity, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-l2_user_base_scale", "streach scale", "",
                          l2UserBaseScale, &l2UserBaseScale, PETSC_NULLPTR);
  
  // dynamic relaxation
  CHKERR PetscOptionsBool("-dynamic_relaxation", "dynamic time relaxation", "",
                          dynamicRelaxation, &dynamicRelaxation, PETSC_NULLPTR);
 
  // contact parameters
  CHKERR PetscOptionsInt("-contact_max_post_proc_ref_level", "refinement level",
                         "", contactRefinementLevels, &contactRefinementLevels,
                         PETSC_NULLPTR);
 
  // cracking parameters
  CHKERR PetscOptionsBool("-cracking_on", "cracking ON", "", crackingOn,
                          &crackingOn, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-cracking_start_time", "cracking start time", "",
                            crackingStartTime, &crackingStartTime, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-griffith_energy", "Griffith energy", "",
                            griffithEnergy, &griffithEnergy, PETSC_NULLPTR);
  CHKERR PetscOptionsEList("-energy_release_variant", "energy release variant",
                           "", list_release, 2, list_release[choice_release],
                           &choice_release, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-nb_J_integral_levels", "Number of J integarl levels",
                         "", nbJIntegralLevels, &nbJIntegralLevels, PETSC_NULLPTR);
 
  // internal stress
  char tag_name[255] = "";
  CHKERR PetscOptionsString("-internal_stress_tag_name",
                            "internal stress tag name", "", "", tag_name, 255,
                            PETSC_NULLPTR);
  internalStressTagName = string(tag_name);
  CHKERR PetscOptionsInt(
      "-internal_stress_interp_order", "internal stress interpolation order",
      "", internalStressInterpOrder, &internalStressInterpOrder, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-internal_stress_voigt", "Voigt index notation", "",
                          internalStressVoigt, &internalStressVoigt,
                          PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-analytical_expr_file",
                               analytical_expr_file_name, 255, PETSC_NULLPTR);
 
  PetscOptionsEnd();
 
  if (internalStressInterpOrder != 0 && internalStressInterpOrder != 1) {
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
             "Unsupported internal stress interpolation order %d",
             internalStressInterpOrder);
  }
 
  if (setSingularity) {
    l2UserBaseScale = PETSC_TRUE;
  } 
 
  EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
  EshelbianCore::gradApproximator = static_cast<RotSelector>(choice_grad);
  EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
  EshelbianCore::symmetrySelector = static_cast<SymmetrySelector>(choice_symm);
  EshelbianCore::energyReleaseSelector =
      static_cast<EnergyReleaseSelector>(choice_release);
 
  switch (EshelbianCore::stretchSelector) {
  case StretchSelector::LINEAR:
    EshelbianCore::f = f_linear;
    EshelbianCore::d_f = d_f_linear;
    EshelbianCore::dd_f = dd_f_linear;
    EshelbianCore::inv_f = inv_f_linear;
    EshelbianCore::inv_d_f = inv_d_f_linear;
    EshelbianCore::inv_dd_f = inv_dd_f_linear;
    break;
  case StretchSelector::LOG:
    if (std::fabs(EshelbianCore::exponentBase - exp(1)) >
        std::numeric_limits<float>::epsilon()) {
      EshelbianCore::f = f_log;
      EshelbianCore::d_f = d_f_log;
      EshelbianCore::dd_f = dd_f_log;
      EshelbianCore::inv_f = inv_f_log;
      EshelbianCore::inv_d_f = inv_d_f_log;
      EshelbianCore::inv_dd_f = inv_dd_f_log;
    } else {
      EshelbianCore::f = f_log_e;
      EshelbianCore::d_f = d_f_log_e;
      EshelbianCore::dd_f = dd_f_log_e;
      EshelbianCore::inv_f = inv_f_log;
      EshelbianCore::inv_d_f = inv_d_f_log;
      EshelbianCore::inv_dd_f = inv_dd_f_log;
    }
    break;
  case StretchSelector::LOG_QUADRATIC:
    EshelbianCore::f = f_log_e_quadratic;
    EshelbianCore::d_f = d_f_log_e_quadratic;
    EshelbianCore::dd_f = dd_f_log_e_quadratic;
    EshelbianCore::inv_f = [](const double x) {
      CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
                        "No logarithmic quadratic stretch for this case");
      return 0;
    };
    EshelbianCore::inv_d_f = EshelbianCore::inv_f;
    EshelbianCore::inv_dd_f = EshelbianCore::inv_f;
    break; // no stretch, do not use stretch functions
  default:
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
    break;
  };
 
  MOFEM_LOG("EP", Sev::inform) << "spaceOrder: -space_order " << spaceOrder;
  MOFEM_LOG("EP", Sev::inform)
      << "spaceH1Order: -space_h1_order " << spaceH1Order;
  MOFEM_LOG("EP", Sev::inform)
      << "materialOrder: -material_order " << materialOrder;
  MOFEM_LOG("EP", Sev::inform) << "alphaU: -viscosity_alpha_u " << alphaU;
  MOFEM_LOG("EP", Sev::inform) << "alphaW: -viscosity_alpha_w " << alphaW;
  MOFEM_LOG("EP", Sev::inform)
      << "alphaOmega: -viscosity_alpha_omega " << alphaOmega;
  MOFEM_LOG("EP", Sev::inform) << "alphaRho: -density_alpha_rho " << alphaRho;
  MOFEM_LOG("EP", Sev::inform)
      << "Rotations: -rotations " << list_rots[EshelbianCore::rotSelector];
  MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
                               << list_rots[EshelbianCore::gradApproximator];
  MOFEM_LOG("EP", Sev::inform)
      << "Symmetry: -symm " << list_symm[EshelbianCore::symmetrySelector];
  if (exponentBase != exp(1))
    MOFEM_LOG("EP", Sev::inform)
        << "Base exponent: -exponent_base " << EshelbianCore::exponentBase;
  else
    MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
  MOFEM_LOG("EP", Sev::inform)
      << "Stretch: -stretches " << list_stretches[choice_stretch];
  MOFEM_LOG("EP", Sev::inform) << "No stretch: -no_stretch " << (noStretch)
      ? "yes"
      : "no";
 
  MOFEM_LOG("EP", Sev::inform)
          << "Dynamic relaxation: -dynamic_relaxation " << (dynamicRelaxation)
      ? "yes"
      : "no";
  MOFEM_LOG("EP", Sev::inform)
          << "Singularity: -set_singularity " << (setSingularity)
      ? "yes"
      : "no";
  MOFEM_LOG("EP", Sev::inform)
          << "L2 user base scale: -l2_user_base_scale " << (l2UserBaseScale)
      ? "yes"
      : "no";
 
  MOFEM_LOG("EP", Sev::inform) << "Cracking on: -cracking_on " << (crackingOn)
      ? "yes"
      : "no";
  MOFEM_LOG("EP", Sev::inform)
      << "Cracking start time: -cracking_start_time " << crackingStartTime;
  MOFEM_LOG("EP", Sev::inform)
      << "Griffith energy: -griffith_energy " << griffithEnergy;
  MOFEM_LOG("EP", Sev::inform)
      << "Energy release variant: -energy_release_variant "
      << list_release[EshelbianCore::energyReleaseSelector];
  MOFEM_LOG("EP", Sev::inform)
      << "Number of J integral levels: -nb_J_integral_levels "
      << nbJIntegralLevels;
 
#ifdef ENABLE_PYTHON_BINDING
  auto file_exists = [](std::string myfile) {
    std::ifstream file(myfile.c_str());
    if (file) {
      return true;
    }
    return false;
  };
 
  if (file_exists(analytical_expr_file_name)) {
    MOFEM_LOG("EP", Sev::inform) << analytical_expr_file_name << " file found";
 
    AnalyticalExprPythonPtr = boost::make_shared<AnalyticalExprPython>();
    CHKERR AnalyticalExprPythonPtr->analyticalExprInit(
        analytical_expr_file_name);
    AnalyticalExprPythonWeakPtr = AnalyticalExprPythonPtr;
  } else {
    MOFEM_LOG("EP", Sev::warning)
        << analytical_expr_file_name << " file NOT found";
  }
#endif
 
  if (spaceH1Order == -1)
    spaceH1Order = spaceOrder;
 
  MoFEMFunctionReturn(0);
}

getSpatialDispBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialDispBc

(

)

[Getting norms]

Definition at line 6778 of file EshelbianPlasticity.cpp.

                                               {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      "", piolaStress, false, false);
 
  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto disp_bc = bc.second->dispBcPtr) {
 
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      MOFEM_LOG("EP", Sev::noisy) << "Displacement BC: " << *disp_bc;
 
      std::vector<double> block_attributes(6, 0.);
      if (disp_bc->data.flag1 == 1) {
        block_attributes[0] = disp_bc->data.value1;
        block_attributes[3] = 1;
      }
      if (disp_bc->data.flag2 == 1) {
        block_attributes[1] = disp_bc->data.value2;
        block_attributes[4] = 1;
      }
      if (disp_bc->data.flag3 == 1) {
        block_attributes[2] = disp_bc->data.value3;
        block_attributes[5] = 1;
      }
      auto faces = bc.second->bcEnts.subset_by_dimension(2);
      bcSpatialDispVecPtr->emplace_back(block_name, block_attributes, faces);
    }
  }
  // old way of naming blocksets for displacement BCs
  CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
 
  bcSpatialNormalDisplacementVecPtr =
      boost::make_shared<NormalDisplacementBcVec>();
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)NORMAL_DISPLACEMENT(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialNormalDisplacementVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
  bcSpatialAnalyticalDisplacementVecPtr =
      boost::make_shared<AnalyticalDisplacementBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)ANALYTICAL_DISPLACEMENT(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialAnalyticalDisplacementVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
  auto ts_displacement =
      boost::make_shared<DynamicRelaxationTimeScale>("disp_history.txt");
  for (auto &bc : *bcSpatialDispVecPtr) {
    MOFEM_LOG("EP", Sev::noisy)
        << "Add time scaling displacement BC: " << bc.blockName;
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_displacement, "disp_history", ".txt", bc.blockName);
  }
 
  auto ts_normal_displacement =
      boost::make_shared<DynamicRelaxationTimeScale>("normal_disp_history.txt");
  for (auto &bc : *bcSpatialNormalDisplacementVecPtr) {
    MOFEM_LOG("EP", Sev::noisy)
        << "Add time scaling normal displacement BC: " << bc.blockName;
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_normal_displacement, "normal_disp_history", ".txt",
            bc.blockName);
  }
 
  MoFEMFunctionReturn(0);
}

getSpatialRotationBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialRotationBc ( )

inline

Definition at line 272 of file EshelbianPlasticity.hpp.

                                          {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wGradH1AtPts);
  }
 
  inline VectorPtr getJacobianAtPts() {
    return boost::shared_ptr<VectorDouble>(shared_from_this(), &jacobianAtPts);
  }
 
  inline MatrixPtr getLeviKirchhoffAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &leviKirchhoffAtPts);
  };
 
  inline MatrixPtr getVarRotAxisPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varRotAxis);
  }

getSpatialTractionBc()

MoFEMErrorCode EshelbianCore::getSpatialTractionBc

(

)

Definition at line 6872 of file EshelbianPlasticity.cpp.

                                                   {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
                                                          false, false);
 
  bcSpatialTractionVecPtr = boost::make_shared<TractionBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto force_bc = bc.second->forceBcPtr) {
 
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      MOFEM_LOG("EP", Sev::noisy) << "Force BC: " << *force_bc;
 
      std::vector<double> block_attributes(6, 0.);
      block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
      block_attributes[3] = 1;
      block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
      block_attributes[4] = 1;
      block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
      block_attributes[5] = 1;
      auto faces = bc.second->bcEnts.subset_by_dimension(2);
      bcSpatialTractionVecPtr->emplace_back(block_name, block_attributes,
                                            faces);
    }
  }
  CHKERR getBc(bcSpatialTractionVecPtr, "SPATIAL_TRACTION_BC", 6);
 
  bcSpatialPressureVecPtr = boost::make_shared<PressureBcVec>();
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)PRESSURE(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
 
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialPressureVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
bcSpatialAnalyticalTractionVecPtr =
      boost::make_shared<AnalyticalTractionBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    auto block_name = "(.*)ANALYTICAL_TRACTION(.*)";
    std::regex reg_name(block_name);
    if (std::regex_match(bc.first, reg_name)) {
      auto [field_name, block_name] =
          BcManager::extractStringFromBlockId(bc.first, "");
      MOFEM_LOG("EP", Sev::inform)
          << "Field name: " << field_name << " Block name: " << block_name;
      bcSpatialAnalyticalTractionVecPtr->emplace_back(
          block_name, bc.second->bcAttributes,
          bc.second->bcEnts.subset_by_dimension(2));
    }
  }
 
  auto ts_traction =
      boost::make_shared<DynamicRelaxationTimeScale>("traction_history.txt");
  for (auto &bc : *bcSpatialTractionVecPtr) {
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_traction, "traction_history", ".txt", bc.blockName);
  }
 
  auto ts_pressure =
      boost::make_shared<DynamicRelaxationTimeScale>("pressure_history.txt");
  for (auto &bc : *bcSpatialPressureVecPtr) {
    timeScaleMap[bc.blockName] =
        GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
            ts_pressure, "pressure_history", ".txt", bc.blockName);
  }
 
  MoFEMFunctionReturn(0);
}

getSpatialTractionFreeBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialTractionFreeBc ( const EntityHandle meshset = 0 )

inline

Definition at line 306 of file EshelbianPlasticity.hpp.

                                       {

gettingNorms()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::gettingNorms

(

)

[Getting norms]

Definition at line 6710 of file EshelbianPlasticity.cpp.

                                           {
  MoFEMFunctionBegin;
 
  auto post_proc_norm_fe =
      boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
 
  auto post_proc_norm_rule_hook = [](int, int, int p) -> int {
    return 2 * (p);
  };
  post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
 
  post_proc_norm_fe->getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
 
  CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV}, materialH1Positions,
      frontAdjEdges);
 
  enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
  auto norms_vec =
      createVectorMPI(mField.get_comm(), LAST_NORM, PETSC_DETERMINE);
  CHKERR VecZeroEntries(norms_vec);
 
  auto u_l2_ptr = boost::make_shared<MatrixDouble>();
  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, u_l2_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(spatialH1Disp, u_h1_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_NORM_L2));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_h1_ptr, norms_vec, U_NORM_H1));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_ERROR_L2,
                                         u_h1_ptr));
 
  auto piola_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateHVecTensorField<3, 3>(piolaStress, piola_ptr));
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateHTensorTensorField<3, 3>(bubbleField, piola_ptr,
                                              MBMAXTYPE));
 
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor2<3, 3>(piola_ptr, norms_vec, PIOLA_NORM));
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
  CHKERR mField.loop_finite_elements("ELASTIC_PROBLEM", elementVolumeName,
                                     *post_proc_norm_fe);
  TetPolynomialBase::switchCacheBaseOff<HDIV>({post_proc_norm_fe.get()});
 
  CHKERR VecAssemblyBegin(norms_vec);
  CHKERR VecAssemblyEnd(norms_vec);
  const double *norms;
  CHKERR VecGetArrayRead(norms_vec, &norms);
  MOFEM_LOG("EP", Sev::inform) << "norm_u: " << std::sqrt(norms[U_NORM_L2]);
  MOFEM_LOG("EP", Sev::inform) << "norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
  MOFEM_LOG("EP", Sev::inform)
      << "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
  MOFEM_LOG("EP", Sev::inform)
      << "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
  CHKERR VecRestoreArrayRead(norms_vec, &norms);
 
  MoFEMFunctionReturn(0);
}

getTractionFreeBc()

MoFEMErrorCode EshelbianCore::getTractionFreeBc	(	const EntityHandle	meshset,
		boost::shared_ptr< TractionFreeBc > &	bc_ptr,
		const std::string	contact_set_name )

Remove all, but entities where kinematic constrains are applied.

Parameters

meshset
bc_ptr
disp_block_set_name
rot_block_set_name
contact_set_name

Returns

MoFEMErrorCode

Definition at line 2083 of file EshelbianPlasticity.cpp.

                                                                   {
  MoFEMFunctionBegin;
 
  // get skin from all tets
  Range tets;
  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
  Range tets_skin_part;
  Skinner skin(&mField.get_moab());
  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  Range tets_skin;
  CHKERR pcomm->filter_pstatus(tets_skin_part,
                               PSTATUS_SHARED | PSTATUS_MULTISHARED,
                               PSTATUS_NOT, -1, &tets_skin);
 
  bc_ptr->resize(3);
  for (int dd = 0; dd != 3; ++dd)
    (*bc_ptr)[dd] = tets_skin;
 
  // Do not remove dofs on which traction is applied
  if (bcSpatialDispVecPtr)
    for (auto &v : *bcSpatialDispVecPtr) {
      if (v.flags[0])
        (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      if (v.flags[1])
        (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      if (v.flags[2])
        (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  // Do not remove dofs on which rotation is applied
  if (bcSpatialRotationVecPtr)
    for (auto &v : *bcSpatialRotationVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialNormalDisplacementVecPtr)
    for (auto &v : *bcSpatialNormalDisplacementVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialAnalyticalDisplacementVecPtr)
    for (auto &v : *bcSpatialAnalyticalDisplacementVecPtr) {
      if (v.flags[0])
        (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      if (v.flags[1])
        (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      if (v.flags[2])
        (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialTractionVecPtr)
    for (auto &v : *bcSpatialTractionVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialAnalyticalTractionVecPtr)
    for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
        
  if (bcSpatialPressureVecPtr)
    for (auto &v : *bcSpatialPressureVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  // remove contact
  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
           std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
    Range faces;
    CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                              true);
    (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
    (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
    (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
  }
 
  MoFEMFunctionReturn(0);
}

inv_d_f_linear()

static double EshelbianPlasticity::EshelbianCore::inv_d_f_linear ( const double v )

inlinestatic

Definition at line 164 of file EshelbianPlasticity.hpp.

inv_d_f_log()

static double EshelbianPlasticity::EshelbianCore::inv_d_f_log ( const double v )

inlinestatic

Definition at line 152 of file EshelbianPlasticity.hpp.

inv_d_f_log_e()

static double EshelbianPlasticity::EshelbianCore::inv_d_f_log_e ( const double v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 125 of file EshelbianPlasticity.hpp.

inv_dd_f_linear()

static double EshelbianPlasticity::EshelbianCore::inv_dd_f_linear ( const double v )

inlinestatic

Definition at line 165 of file EshelbianPlasticity.hpp.

inv_dd_f_log()

static double EshelbianPlasticity::EshelbianCore::inv_dd_f_log ( const double v )

inlinestatic

Definition at line 155 of file EshelbianPlasticity.hpp.

inv_dd_f_log_e()

static double EshelbianPlasticity::EshelbianCore::inv_dd_f_log_e ( const double v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 131 of file EshelbianPlasticity.hpp.

inv_f_linear()

static double EshelbianPlasticity::EshelbianCore::inv_f_linear ( const double v )

inlinestatic

Definition at line 163 of file EshelbianPlasticity.hpp.

inv_f_log()

static double EshelbianPlasticity::EshelbianCore::inv_f_log ( const double v )

inlinestatic

Definition at line 149 of file EshelbianPlasticity.hpp.

inv_f_log_e()

static double EshelbianPlasticity::EshelbianCore::inv_f_log_e ( const double v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 116 of file EshelbianPlasticity.hpp.

postProcessResults()

MoFEMErrorCode EshelbianCore::postProcessResults	(	const int	tag,
		const std::string	file,
		Vec	f_residual = PETSC_NULLPTR,
		Vec	var_vec = PETSC_NULLPTR,
		std::vector< Tag >	tags_to_transfer = {} )

Definition at line 3723 of file EshelbianPlasticity.cpp.

                                                                   {
  MoFEMFunctionBegin;
 
  // mark crack surface
  if (crackingOn) {
    Tag th;
    rval = mField.get_moab().tag_get_handle("CRACK", th);
    if (rval == MB_SUCCESS) {
      CHKERR mField.get_moab().tag_delete(th);
    }
    int def_val[] = {0};
    CHKERR mField.get_moab().tag_get_handle(
        "CRACK", 1, MB_TYPE_INTEGER, th, MB_TAG_SPARSE | MB_TAG_CREAT, def_val);
    int mark[] = {1};
    CHKERR mField.get_moab().tag_clear_data(th, *crackFaces, mark);
    tags_to_transfer.push_back(th);
 
    auto get_tag = [&](auto name, auto dim) {
      auto &mob = mField.get_moab();
      Tag tag;
      double def_val[] = {0., 0., 0.};
      CHK_MOAB_THROW(mob.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
                                        MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
                     "create tag");
      return tag;
    };
 
    tags_to_transfer.push_back(get_tag("MaterialForce", 3));
    // tags_to_transfer.push_back(get_tag("GriffithForce", 1));
  }
 
  // add tags to transfer
  for (auto t : listTagsToTransfer) {
    tags_to_transfer.push_back(t);
  }
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
  }
 
  struct exclude_sdf {
    exclude_sdf(Range &&r) : map(r) {}
    bool operator()(FEMethod *fe_method_ptr) {
      auto ent = fe_method_ptr->getFEEntityHandle();
      if (map.find(ent) != map.end()) {
        return false;
      }
      return true;
    }
 
  private:
    Range map;
  };
 
  contactTreeRhs->exeTestHook =
      exclude_sdf(get_range_from_block(mField, "CONTACT_SDF", SPACE_DIM - 1));
 
  CHKERR DMoFEMLoopFiniteElements(dM, contactElement, contactTreeRhs);
 
  auto get_post_proc = [&](auto &post_proc_mesh, auto sense) {
    using FaceEle = FaceElementForcesAndSourcesCore;
    auto post_proc_ptr =
        boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
            mField, post_proc_mesh);
    EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
        frontAdjEdges);
 
    auto domain_ops = [&](auto &fe, int sense) {
      MoFEMFunctionBegin;
      fe.getUserPolynomialBase() =
          boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          fe.getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions,
          frontAdjEdges);
      auto piola_scale_ptr = boost::make_shared<double>(1.0);
      fe.getOpPtrVector().push_back(physicalEquations->returnOpSetScale(
          piola_scale_ptr, physicalEquations));
      fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
          piolaStress, dataAtPts->getApproxPAtPts(), piola_scale_ptr));
      fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
          bubbleField, dataAtPts->getApproxPAtPts(), piola_scale_ptr,
          SmartPetscObj<Vec>(), MBMAXTYPE));
      if (noStretch) {
        fe.getOpPtrVector().push_back(
            physicalEquations->returnOpCalculateStretchFromStress(
                dataAtPts, physicalEquations));
      } else {
        fe.getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<3>(
                stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
      }
      if (var_vector) {
        auto vec = SmartPetscObj<Vec>(var_vector, true);
        fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
            piolaStress, dataAtPts->getVarPiolaPts(),
            boost::make_shared<double>(1), vec));
        fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
            bubbleField, dataAtPts->getVarPiolaPts(),
            boost::make_shared<double>(1), vec, MBMAXTYPE));
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            rotAxis, dataAtPts->getVarRotAxisPts(), vec, MBTET));
        if (noStretch) {
          fe.getOpPtrVector().push_back(
              physicalEquations->returnOpCalculateVarStretchFromStress(
                  dataAtPts, physicalEquations));
        } else {
          fe.getOpPtrVector().push_back(
              new OpCalculateTensor2SymmetricFieldValues<3>(
                  stretchTensor, dataAtPts->getVarLogStreachPts(), vec, MBTET));
        }
      }
 
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
 
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
      fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
          spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
      // evaluate derived quantities
      fe.getOpPtrVector().push_back(
          new OpCalculateRotationAndSpatialGradient(dataAtPts));
 
      // evaluate integration points
      fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
          tag, true, false, dataAtPts, physicalEquations));
      if (auto op =
              physicalEquations->returnOpCalculateEnergy(dataAtPts, nullptr)) {
        fe.getOpPtrVector().push_back(op);
        fe.getOpPtrVector().push_back(new OpCalculateEshelbyStress(dataAtPts));
      }
 
      // // post-proc
      using OpPPMap = OpPostProcMapInMoab<
          SPACE_DIM, SPACE_DIM,
          VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator>;
 
      struct OpSidePPMap : public OpPPMap {
        OpSidePPMap(moab::Interface &post_proc_mesh,
                    std::vector<EntityHandle> &map_gauss_pts,
                    DataMapVec data_map_scalar, DataMapMat data_map_vec,
                    DataMapMat data_map_mat, DataMapMat data_symm_map_mat,
                    int sense)
            : OpPPMap(post_proc_mesh, map_gauss_pts, data_map_scalar,
                      data_map_vec, data_map_mat, data_symm_map_mat),
              tagSense(sense) {}
 
        MoFEMErrorCode doWork(int side, EntityType type,
                              EntitiesFieldData::EntData &data) {
          MoFEMFunctionBegin;
 
          if (tagSense != OpPPMap::getSkeletonSense())
            MoFEMFunctionReturnHot(0);
 
          CHKERR OpPPMap::doWork(side, type, data);
          MoFEMFunctionReturn(0);
        }
 
      private:
        int tagSense;
      };
 
      OpPPMap::DataMapMat vec_fields;
      vec_fields["SpatialDisplacementL2"] = dataAtPts->getSmallWL2AtPts();
      vec_fields["SpatialDisplacementH1"] = dataAtPts->getSmallWH1AtPts();
      vec_fields["Omega"] = dataAtPts->getRotAxisAtPts();
      vec_fields["ContactDisplacement"] = dataAtPts->getContactL2AtPts();
      vec_fields["AngularMomentum"] = dataAtPts->getLeviKirchhoffAtPts();
      vec_fields["X"] = dataAtPts->getLargeXH1AtPts();
      if (!noStretch) {
        vec_fields["EiegnLogStreach"] = dataAtPts->getEigenValsAtPts();
      }
      if (var_vector) {
        auto vec = SmartPetscObj<Vec>(var_vector, true);
        vec_fields["VarOmega"] = dataAtPts->getVarRotAxisPts();
        vec_fields["VarSpatialDisplacementL2"] =
            boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            spatialL2Disp, vec_fields["VarSpatialDisplacementL2"], vec, MBTET));
      }
      if (f_residual) {
        auto vec = SmartPetscObj<Vec>(f_residual, true);
        vec_fields["ResSpatialDisplacementL2"] =
            boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            spatialL2Disp, vec_fields["ResSpatialDisplacementL2"], vec, MBTET));
        vec_fields["ResOmega"] = boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            rotAxis, vec_fields["ResOmega"], vec, MBTET));
      }
 
      OpPPMap::DataMapMat mat_fields;
      mat_fields["PiolaStress"] = dataAtPts->getApproxPAtPts();
      if (var_vector) {
        mat_fields["VarPiolaStress"] = dataAtPts->getVarPiolaPts();
      }
      if (f_residual) {
        auto vec = SmartPetscObj<Vec>(f_residual, true);
        mat_fields["ResPiolaStress"] = boost::make_shared<MatrixDouble>();
        fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
            piolaStress, mat_fields["ResPiolaStress"],
            boost::make_shared<double>(1), vec));
        fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
            bubbleField, mat_fields["ResPiolaStress"],
            boost::make_shared<double>(1), vec, MBMAXTYPE));
      }
      if (!internalStressTagName.empty()) {
        mat_fields[internalStressTagName] = dataAtPts->getInternalStressAtPts();
        switch (internalStressInterpOrder) {
        case 0:
          fe.getOpPtrVector().push_back(
              new OpGetInternalStress<0>(dataAtPts, internalStressTagName));
          break;
        case 1:
          fe.getOpPtrVector().push_back(
              new OpGetInternalStress<1>(dataAtPts, internalStressTagName));
          break;
        default:
          SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
                   "Unsupported internal stress interpolation order %d",
                   internalStressInterpOrder);
        }
      }
 
      OpPPMap::DataMapMat mat_fields_symm;
      mat_fields_symm["LogSpatialStretch"] =
          dataAtPts->getLogStretchTensorAtPts();
      mat_fields_symm["SpatialStretch"] = dataAtPts->getStretchTensorAtPts();
      if (var_vector) {
        mat_fields_symm["VarLogSpatialStretch"] =
            dataAtPts->getVarLogStreachPts();
      }
      if (f_residual) {
        auto vec = SmartPetscObj<Vec>(f_residual, true);
        if (!noStretch) {
          mat_fields_symm["ResLogSpatialStretch"] =
              boost::make_shared<MatrixDouble>();
          fe.getOpPtrVector().push_back(
              new OpCalculateTensor2SymmetricFieldValues<3>(
                  stretchTensor, mat_fields_symm["ResLogSpatialStretch"], vec,
                  MBTET));
        }
      }
 
      fe.getOpPtrVector().push_back(
 
          new OpSidePPMap(
 
              post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
              {},
 
              vec_fields,
 
              mat_fields,
 
              mat_fields_symm,
 
              sense
 
              )
 
      );
 
      fe.getOpPtrVector().push_back(new OpPostProcDataStructure(
          post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
          dataAtPts, sense));
 
      MoFEMFunctionReturn(0);
    };
 
    post_proc_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(contactDisp,
                                            dataAtPts->getContactL2AtPts()));
 
    auto X_h1_ptr = boost::make_shared<MatrixDouble>();
    // H1 material positions
    post_proc_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(materialH1Positions,
                                            dataAtPts->getLargeXH1AtPts()));
 
    // domain
    auto op_loop_side = new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
        mField, elementVolumeName, SPACE_DIM);
    domain_ops(*(op_loop_side->getSideFEPtr()), sense);
    post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
 
    return post_proc_ptr;
  };
 
  // contact
  auto calcs_side_traction_and_displacements = [&](auto &post_proc_ptr,
                                                   auto &pip) {
    MoFEMFunctionBegin;
    auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
    // evaluate traction
    using EleOnSide =
        PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
    using SideEleOp = EleOnSide::UserDataOperator;
    auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
        mField, elementVolumeName, SPACE_DIM, Sev::noisy);
    op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
        boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
    EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
        materialH1Positions, frontAdjEdges);
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
            piolaStress, contact_common_data_ptr->contactTractionPtr(),
            boost::make_shared<double>(1.0)));
    pip.push_back(op_loop_domain_side);
    // evaluate contact displacement and contact conditions
    auto u_h1_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
        contactDisp, contact_common_data_ptr->contactDispPtr()));
    pip.push_back(new OpTreeSearch(
        contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
        get_range_from_block(mField, "CONTACT", SPACE_DIM - 1),
        &post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
 
    if (f_residual) {
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      auto op_this = new OpLoopThis<BoundaryEle>(mField, contactElement);
      pip.push_back(op_this);
      auto contact_residual = boost::make_shared<MatrixDouble>();
      op_this->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(
              contactDisp, contact_residual,
              SmartPetscObj<Vec>(f_residual, true)));
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
      op_this->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
              {},
 
              {{"res_contact", contact_residual}},
 
              {},
 
              {}
 
              )
 
      );
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto post_proc_mesh = boost::make_shared<moab::Core>();
  auto post_proc_ptr = get_post_proc(post_proc_mesh, 1);
  auto post_proc_negative_sense_ptr = get_post_proc(post_proc_mesh, -1);
  CHKERR calcs_side_traction_and_displacements(post_proc_ptr,
                                               post_proc_ptr->getOpPtrVector());
 
  auto own_tets =
      CommInterface::getPartEntities(mField.get_moab(), mField.get_comm_rank())
          .subset_by_dimension(SPACE_DIM);
  Range own_faces;
  CHKERR mField.get_moab().get_adjacencies(own_tets, SPACE_DIM - 1, true,
                                           own_faces, moab::Interface::UNION);
 
  auto get_post_negative = [&](auto &&ents) {
    auto crack_faces_pos = ents;
    auto crack_faces_neg = crack_faces_pos;
    auto skin = get_skin(mField, own_tets);
    auto crack_on_proc_skin = intersect(crack_faces_pos, skin);
    for (auto f : crack_on_proc_skin) {
      Range tet;
      CHKERR mField.get_moab().get_adjacencies(&f, 1, SPACE_DIM, false, tet);
      tet = intersect(tet, own_tets);
      int side_number, sense, offset;
      CHKERR mField.get_moab().side_number(tet[0], f, side_number, sense,
                                           offset);
      if (sense == 1) {
        crack_faces_neg.erase(f);
      } else {
        crack_faces_pos.erase(f);
      }
    }
    return std::make_pair(crack_faces_pos, crack_faces_neg);
  };
 
  auto get_crack_faces = [&](auto crack_faces) {
    auto get_adj = [&](auto e, auto dim) {
      Range adj;
      CHKERR mField.get_moab().get_adjacencies(e, dim, true, adj,
                                               moab::Interface::UNION);
      return adj;
    };
    auto tets = get_adj(crack_faces, 3);
    auto faces = subtract(get_adj(tets, 2), crack_faces);
    tets = subtract(tets, get_adj(faces, 3));
    return subtract(crack_faces, get_adj(tets, 2));
  };
 
  auto [crack_faces_pos, crack_faces_neg] =
      get_post_negative(intersect(own_faces, get_crack_faces(*crackFaces)));
 
  auto only_crack_faces = [&](FEMethod *fe_method_ptr) {
    auto ent = fe_method_ptr->getFEEntityHandle();
    if (crack_faces_pos.find(ent) == crack_faces_pos.end()) {
      return false;
    }
    return true;
  };
 
  auto only_negative_crack_faces = [&](FEMethod *fe_method_ptr) {
    auto ent = fe_method_ptr->getFEEntityHandle();
    if (crack_faces_neg.find(ent) == crack_faces_neg.end()) {
      return false;
    }
    return true;
  };
 
  auto get_adj_front = [&]() {
    auto skeleton_faces = *skeletonFaces;
    Range adj_front;
    CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2, true, adj_front,
                                             moab::Interface::UNION);
 
    adj_front = intersect(adj_front, skeleton_faces);
    adj_front = subtract(adj_front, *crackFaces);
    adj_front = intersect(own_faces, adj_front);
    return adj_front;
  };
 
  post_proc_ptr->setTagsToTransfer(tags_to_transfer);
  post_proc_negative_sense_ptr->setTagsToTransfer(tags_to_transfer);
 
  auto post_proc_begin =
      PostProcBrokenMeshInMoabBaseBegin(mField, post_proc_mesh);
  CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
  CHKERR DMoFEMLoopFiniteElements(dM, skinElement, post_proc_ptr);
  post_proc_ptr->exeTestHook = only_crack_faces;
  post_proc_negative_sense_ptr->exeTestHook = only_negative_crack_faces;
  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
      dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, skeletonElement,
                                              post_proc_negative_sense_ptr, 0,
                                              mField.get_comm_size());
 
  constexpr bool debug = false;
  if (debug) {
 
    auto [adj_front_pos, adj_front_neg] =
        get_post_negative(filter_owners(mField, get_adj_front()));
 
    auto only_front_faces_pos = [adj_front_pos](FEMethod *fe_method_ptr) {
      auto ent = fe_method_ptr->getFEEntityHandle();
      if (adj_front_pos.find(ent) == adj_front_pos.end()) {
        return false;
      }
      return true;
    };
 
    auto only_front_faces_neg = [adj_front_neg](FEMethod *fe_method_ptr) {
      auto ent = fe_method_ptr->getFEEntityHandle();
      if (adj_front_neg.find(ent) == adj_front_neg.end()) {
        return false;
      }
      return true;
    };
 
    post_proc_ptr->exeTestHook = only_front_faces_pos;
    CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
        dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
    post_proc_negative_sense_ptr->exeTestHook = only_front_faces_neg;
    CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, skeletonElement,
                                                post_proc_negative_sense_ptr, 0,
                                                mField.get_comm_size());
  }
  auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
  CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
 
  CHKERR post_proc_end.writeFile(file.c_str());
  MoFEMFunctionReturn(0);
}

postProcessSkeletonResults()

MoFEMErrorCode EshelbianCore::postProcessSkeletonResults	(	const int	tag,
		const std::string	file,
		Vec	f_residual = PETSC_NULLPTR,
		std::vector< Tag >	tags_to_transfer = {} )

Definition at line 4217 of file EshelbianPlasticity.cpp.

                                                                           {
  MoFEMFunctionBegin;
 
  using FaceEle = FaceElementForcesAndSourcesCore;
 
  auto post_proc_mesh = boost::make_shared<moab::Core>();
  auto post_proc_ptr =
      boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
          mField, post_proc_mesh);
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
      frontAdjEdges);
 
  auto hybrid_disp = boost::make_shared<MatrixDouble>();
  post_proc_ptr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>(hybridSpatialDisp, hybrid_disp));
 
  auto hybrid_res = boost::make_shared<MatrixDouble>();
  post_proc_ptr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>(hybridSpatialDisp, hybrid_disp));
 
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
  post_proc_ptr->getOpPtrVector().push_back(
 
      new OpPPMap(
 
          post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
          {},
 
          {{"hybrid_disp", hybrid_disp}},
 
          {},
 
          {}
 
          )
 
  );
 
  if (f_residual) {
    auto hybrid_res = boost::make_shared<MatrixDouble>();
    post_proc_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(
            hybridSpatialDisp, hybrid_res,
            SmartPetscObj<Vec>(f_residual, true)));
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
    post_proc_ptr->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
            {},
 
            {{"res_hybrid", hybrid_res}},
 
            {},
 
            {}
 
            )
 
    );
  }
 
  post_proc_ptr->setTagsToTransfer(tags_to_transfer);
 
  auto post_proc_begin =
      PostProcBrokenMeshInMoabBaseBegin(mField, post_proc_mesh);
  CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
  CHKERR DMoFEMLoopFiniteElements(dM, skeletonElement, post_proc_ptr);
  auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
  CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
 
  CHKERR post_proc_end.writeFile(file.c_str());
 
  MoFEMFunctionReturn(0);
}

projectGeometry()

MoFEMErrorCode EshelbianCore::projectGeometry

(

const EntityHandle

meshset = 0

)

Definition at line 1403 of file EshelbianPlasticity.cpp.

                                                                        {
  MoFEMFunctionBegin;
 
  Range meshset_ents;
  CHKERR mField.get_moab().get_entities_by_handle(meshset, meshset_ents);
 
  auto project_ho_geometry = [&](auto field) {
    Projection10NodeCoordsOnField ent_method(mField, materialH1Positions);
    return mField.loop_dofs(field, ent_method);
  };
  CHKERR project_ho_geometry(materialH1Positions);
 
  auto get_adj_front_edges = [&](auto &front_edges) {
    Range front_crack_nodes;
    Range crack_front_edges_with_both_nodes_not_at_front;
 
    if (mField.get_comm_rank() == 0) {
      auto &moab = mField.get_moab();
      CHKERR moab.get_connectivity(front_edges, front_crack_nodes, true);
      Range crack_front_edges;
      CHKERR moab.get_adjacencies(front_crack_nodes, SPACE_DIM - 2, false,
                                  crack_front_edges, moab::Interface::UNION);
      crack_front_edges =
          intersect(subtract(crack_front_edges, front_edges), meshset_ents);
      Range crack_front_edges_nodes;
      CHKERR moab.get_connectivity(crack_front_edges, crack_front_edges_nodes,
                                   true);
      crack_front_edges_nodes =
          subtract(crack_front_edges_nodes, front_crack_nodes);
      CHKERR moab.get_adjacencies(
          crack_front_edges_nodes, SPACE_DIM - 2, false,
          crack_front_edges_with_both_nodes_not_at_front,
          moab::Interface::UNION);
      crack_front_edges_with_both_nodes_not_at_front = intersect(
          crack_front_edges_with_both_nodes_not_at_front, meshset_ents);
      crack_front_edges_with_both_nodes_not_at_front = intersect(
          crack_front_edges_with_both_nodes_not_at_front, crack_front_edges);
    }
 
    front_crack_nodes = send_type(mField, front_crack_nodes, MBVERTEX);
    crack_front_edges_with_both_nodes_not_at_front = send_type(
        mField, crack_front_edges_with_both_nodes_not_at_front, MBEDGE);
 
    return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
                          boost::make_shared<Range>(
                              crack_front_edges_with_both_nodes_not_at_front));
  };
 
  crackFaces = boost::make_shared<Range>(
      get_range_from_block(mField, "CRACK", SPACE_DIM - 1));
  frontEdges =
      boost::make_shared<Range>(get_crack_front_edges(mField, *crackFaces));
  auto [front_vertices, front_adj_edges] = get_adj_front_edges(*frontEdges);
  frontVertices = front_vertices;
  frontAdjEdges = front_adj_edges;
 
#ifndef NDEBUG
  if (crackingOn) {
    auto rank = mField.get_comm_rank();
    // CHKERR save_range(mField.get_moab(),
    //                   (boost::format("meshset_ents_%d.vtk") % rank).str(),
    //                   meshset_ents);
    CHKERR save_range(mField.get_moab(),
                      (boost::format("crack_faces_%d.vtk") % rank).str(),
                      *crackFaces);
    CHKERR save_range(mField.get_moab(),
                      (boost::format("front_edges_%d.vtk") % rank).str(),
                      *frontEdges);
    // CHKERR save_range(mField.get_moab(),
    //                   (boost::format("front_vertices_%d.vtk") % rank).str(),
    //                   *frontVertices);
    // CHKERR save_range(mField.get_moab(),
    //                   (boost::format("front_adj_edges_%d.vtk") % rank).str(),
    //                   *frontAdjEdges);
  }
#endif // NDEBUG
 
  auto set_singular_dofs = [&](auto &front_adj_edges, auto &front_vertices) {
    MoFEMFunctionBegin;
    auto &moab = mField.get_moab();
 
    double eps = 1;
    double beta = 0;
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-singularity_eps", &beta,
                                 PETSC_NULLPTR);
    MOFEM_LOG("EP", Sev::inform) << "Singularity eps " << beta;
    eps -= beta;
 
    for (auto edge : front_adj_edges) {
      int num_nodes;
      const EntityHandle *conn;
      CHKERR moab.get_connectivity(edge, conn, num_nodes, false);
      double coords[6];
      CHKERR moab.get_coords(conn, num_nodes, coords);
      const double dir[3] = {coords[3] - coords[0], coords[4] - coords[1],
                             coords[5] - coords[2]};
      double dof[3] = {0, 0, 0};
      if (front_vertices.find(conn[0]) != front_vertices.end()) {
        for (int dd = 0; dd != 3; dd++) {
          dof[dd] = -dir[dd] * eps;
        }
      } else if (front_vertices.find(conn[1]) != front_vertices.end()) {
        for (int dd = 0; dd != 3; dd++) {
          dof[dd] = +dir[dd] * eps;
        }
      }
      for (_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(
               mField, materialH1Positions, edge, dit)) {
        const int idx = dit->get()->getEntDofIdx();
        if (idx > 2) {
          dit->get()->getFieldData() = 0;
        } else {
          dit->get()->getFieldData() = dof[idx];
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  if (setSingularity)
    CHKERR set_singular_dofs(*frontAdjEdges, *frontVertices);
 
  MoFEMFunctionReturn(0);
}

query_interface()

MoFEMErrorCode EshelbianCore::query_interface	(	boost::typeindex::type_index	type_index,
		UnknownInterface **	iface ) const

Getting interface of core database.

Parameters

uuid	unique ID of interface
iface	returned pointer to interface

Returns

error code

Examples

EshelbianPlasticity.cpp.

Definition at line 895 of file EshelbianPlasticity.cpp.

                                                               {
  *iface = const_cast<EshelbianCore *>(this);
  return 0;
}

saveOrgCoords()

MoFEMErrorCode EshelbianCore::saveOrgCoords

(

)

Definition at line 6678 of file EshelbianPlasticity.cpp.

                                            {
  MoFEMFunctionBegin;
  auto crack_faces =
      get_range_from_block(mField, "CRACK_COMPUTED", SPACE_DIM - 1);
  Range conn;
  CHKERR mField.get_moab().get_connectivity(crack_faces, conn, true);
  Range verts;
  CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, verts);
  verts = subtract(verts, conn);
  std::vector<double> coords(3 * verts.size());
  CHKERR mField.get_moab().get_coords(verts, coords.data());
  double def_coords[] = {0., 0., 0.};
  Tag th_org_coords;
  CHKERR mField.get_moab().tag_get_handle(
      "ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
      MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
  CHKERR mField.get_moab().tag_set_data(th_org_coords, verts, coords.data());
  MoFEMFunctionReturn(0);
}

setBaseVolumeElementOps()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setBaseVolumeElementOps	(	const int	tag,
		const bool	do_rhs,
		const bool	do_lhs,
		const bool	calc_rates,
		SmartPetscObj< Vec >	ver_vec,
		boost::shared_ptr< VolumeElementForcesAndSourcesCore >	fe )

Definition at line 2317 of file EshelbianPlasticity.cpp.

                                                           {
  MoFEMFunctionBegin;
 
  auto bubble_cache =
      boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0, MatrixDouble());
  fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      fe->getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions, frontAdjEdges);
 
  // set integration rule
  fe->getRuleHook = [](int, int, int) { return -1; };
  fe->setRuleHook = SetIntegrationAtFrontVolume(frontVertices, frontAdjEdges);
  // fe->getRuleHook = VolRule();
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
    dataAtPts->physicsPtr = physicalEquations;
  }
 
  // calculate fields values
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      piolaStress, dataAtPts->getApproxPAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
      bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
      piolaStress, dataAtPts->getDivPAtPts()));
 
  if (noStretch) {
    fe->getOpPtrVector().push_back(
        physicalEquations->returnOpCalculateStretchFromStress(
            dataAtPts, physicalEquations));
  } else {
    fe->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<3>(
            stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
  }
 
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
 
  // H1 displacements
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
      spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
 
  // velocities
  if (calc_rates) {
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
        spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
    if (noStretch) {
    } else {
      fe->getOpPtrVector().push_back(
          new OpCalculateTensor2SymmetricFieldValuesDot<3>(
              stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
      fe->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradientDot<6, 3>(
              stretchTensor, dataAtPts->getGradLogStretchDotTensorAtPts(),
              MBTET));
    }
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
        rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradientDot<3, 3>(
        rotAxis, dataAtPts->getRotAxisGradDotAtPts(), MBTET));
 
    // acceleration
    if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
      fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
          spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
    }
  }
 
  // variations
  if (var_vec.use_count()) {
    fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
        piolaStress, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
        var_vec));
    fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
        bubbleField, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
        var_vec, MBMAXTYPE));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        rotAxis, dataAtPts->getVarRotAxisPts(), var_vec, MBTET));
    fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
        piolaStress, dataAtPts->getDivVarPiolaPts(), var_vec));
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        spatialL2Disp, dataAtPts->getVarWL2Pts(), var_vec, MBTET));
 
    if (noStretch) {
      fe->getOpPtrVector().push_back(
          physicalEquations->returnOpCalculateVarStretchFromStress(
              dataAtPts, physicalEquations));
    } else {
      fe->getOpPtrVector().push_back(
          new OpCalculateTensor2SymmetricFieldValues<3>(
              stretchTensor, dataAtPts->getVarLogStreachPts(), var_vec, MBTET));
    }
  }
 
  // calculate other derived quantities
  fe->getOpPtrVector().push_back(new OpCalculateRotationAndSpatialGradient(
      dataAtPts, ((do_rhs || do_lhs) && calc_rates) ? alphaOmega : 0.));
 
  // evaluate integration points
  if (noStretch) {
  } else {
    fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
        tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
  }
 
  MoFEMFunctionReturn(0);
}

setBlockTagsOnSkin()

MoFEMErrorCode EshelbianCore::setBlockTagsOnSkin

(

)

Definition at line 3662 of file EshelbianPlasticity.cpp.

                                                 {
  MoFEMFunctionBegin;
 
  auto set_block = [&](auto name, int dim) {
    std::map<int, Range> map;
    auto set_tag_impl = [&](auto name) {
      MoFEMFunctionBegin;
      auto mesh_mng = mField.getInterface<MeshsetsManager>();
      auto bcs = mesh_mng->getCubitMeshsetPtr(
 
          std::regex((boost::format("%s(.*)") % name).str())
 
      );
      for (auto bc : bcs) {
        Range r;
        CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), dim, r,
                                                   true);
        map[bc->getMeshsetId()] = r;
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_tag_impl(name);
 
    return std::make_pair(name, map);
  };
 
  auto set_skin = [&](auto &&map) {
    for (auto &m : map.second) {
      auto s = filter_true_skin(mField, get_skin(mField, m.second));
      m.second.swap(s);
    }
    return map;
  };
 
  auto set_tag = [&](auto &&map) {
    Tag th;
    auto name = map.first;
    int def_val[] = {-1};
    CHK_MOAB_THROW(
        mField.get_moab().tag_get_handle(name, 1, MB_TYPE_INTEGER, th,
                                         MB_TAG_SPARSE | MB_TAG_CREAT, def_val),
        "create tag");
    for (auto &m : map.second) {
      int id = m.first;
      CHK_MOAB_THROW(mField.get_moab().tag_clear_data(th, m.second, &id),
                     "clear tag");
    }
    return th;
  };
 
  listTagsToTransfer.push_back(set_tag(set_skin(set_block("BODY", 3))));
  listTagsToTransfer.push_back(set_tag(set_skin(set_block("MAT_ELASTIC", 3))));
  listTagsToTransfer.push_back(
      set_tag(set_skin(set_block("MAT_NEOHOOKEAN", 3))));
  listTagsToTransfer.push_back(set_tag(set_block("CONTACT", 2)));
 
  MoFEMFunctionReturn(0);
}

setContactElementRhsOps()

MoFEMErrorCode EshelbianCore::setContactElementRhsOps

(

boost::shared_ptr< ContactTree > &

fe_contact_tree

)

Contact requires that body is marked

Definition at line 3089 of file EshelbianPlasticity.cpp.

  {
  MoFEMFunctionBegin;
 
  /** Contact requires that body is marked */
  auto get_body_range = [this](auto name, int dim, auto sev) {
    std::map<int, Range> map;
 
    for (auto m_ptr :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % name).str()
 
                 ))
 
    ) {
      Range ents;
      CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
                                                            dim, ents, true),
                     "by dim");
      map[m_ptr->getMeshsetId()] = ents;
      MOFEM_LOG("EPSYNC", sev) << "Meshset: " << m_ptr->getMeshsetId() << " "
                               << ents.size() << " entities";
    }
 
    MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), sev);
    return map;
  };
 
  auto get_map_skin = [this](auto &&map) {
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
 
    Skinner skin(&mField.get_moab());
    for (auto &m : map) {
      Range skin_faces;
      CHKERR skin.find_skin(0, m.second, false, skin_faces);
      CHK_MOAB_THROW(pcomm->filter_pstatus(skin_faces,
                                           PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                           PSTATUS_NOT, -1, nullptr),
                     "filter");
      m.second.swap(skin_faces);
    }
    return map;
  };
 
  /* The above code is written in C++ and it appears to be defining and using
  various operations on boundary elements and side elements. */
  using BoundaryEle = FaceElementForcesAndSourcesCore;
  using BoundaryEleOp = BoundaryEle::UserDataOperator;
 
  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
 
  auto calcs_side_traction = [&](auto &pip) {
    MoFEMFunctionBegin;
    using EleOnSide =
        PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
    using SideEleOp = EleOnSide::UserDataOperator;
    auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
        mField, elementVolumeName, SPACE_DIM, Sev::noisy);
    op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
        boost::make_shared<CGGUserPolynomialBase>();
    EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
        materialH1Positions, frontAdjEdges);
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
            piolaStress, contact_common_data_ptr->contactTractionPtr(),
            boost::make_shared<double>(1.0)));
    pip.push_back(op_loop_domain_side);
    MoFEMFunctionReturn(0);
  };
 
  auto add_contact_three = [&]() {
    MoFEMFunctionBegin;
    auto tree_moab_ptr = boost::make_shared<moab::Core>();
    fe_contact_tree = boost::make_shared<ContactTree>(
        mField, tree_moab_ptr, spaceOrder,
        get_body_range("CONTACT", SPACE_DIM - 1, Sev::inform));
    fe_contact_tree->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(
            contactDisp, contact_common_data_ptr->contactDispPtr()));
    CHKERR calcs_side_traction(fe_contact_tree->getOpPtrVector());
    auto u_h1_ptr = boost::make_shared<MatrixDouble>();
    fe_contact_tree->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
    fe_contact_tree->getOpPtrVector().push_back(
        new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_contact_three();
 
  MoFEMFunctionReturn(0);
}

setElasticElementOps()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setElasticElementOps

(

const int

tag

)

Definition at line 3188 of file EshelbianPlasticity.cpp.

                                                                {
  MoFEMFunctionBegin;
 
  // Add contact operators. Note that only for rhs. THe lhs is assembled with
  // volume element, to enable schur complement evaluation.
  CHKERR setContactElementRhsOps(contactTreeRhs);
 
  CHKERR setVolumeElementOps(tag, true, false, elasticFeRhs, elasticFeLhs);
  CHKERR setFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
 
  auto adj_cache =
      boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
 
  auto get_op_contact_bc = [&]() {
    using SideEle = FaceElementForcesAndSourcesCore;
    auto op_loop_side = new OpLoopSide<SideEle>(
        mField, contactElement, SPACE_DIM - 1, Sev::noisy, adj_cache);
    return op_loop_side;
  };
 
  MoFEMFunctionReturn(0);
}

setElasticElementToTs()

MoFEMErrorCode EshelbianCore::setElasticElementToTs

(

DM

dm

)

Definition at line 3211 of file EshelbianPlasticity.cpp.

                                                         {
  MoFEMFunctionBegin;
  boost::shared_ptr<FEMethod> null;
 
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
 
    CHKERR DMMoFEMTSSetI2Function(dm, elementVolumeName, elasticFeRhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Function(dm, naturalBcElement, elasticBcRhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
                                  null);
 
  } else {
    CHKERR DMMoFEMTSSetIFunction(dm, elementVolumeName, elasticFeRhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIFunction(dm, naturalBcElement, elasticBcRhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
                                 null);
  }
 
  MoFEMFunctionReturn(0);
}

setFaceElementOps()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setFaceElementOps	(	const bool	add_elastic,
		const bool	add_material,
		boost::shared_ptr< FaceElementForcesAndSourcesCore > &	fe_rhs,
		boost::shared_ptr< FaceElementForcesAndSourcesCore > &	fe_lhs )

Definition at line 2930 of file EshelbianPlasticity.cpp.

                                                              {
  MoFEMFunctionBegin;
 
  fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
  fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
 
  // set integration rule
  // fe_rhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
  // fe_lhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
  fe_rhs->getRuleHook = [](int, int, int) { return -1; };
  fe_lhs->getRuleHook = [](int, int, int) { return -1; };
  fe_rhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
  fe_lhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
 
  CHKERR
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe_rhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
  CHKERR
  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe_lhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
 
  if (add_elastic) {
 
    auto get_broken_op_side = [this](auto &pip) {
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      // Iterate over domain FEs adjacent to boundary.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
      // Note: EleOnSide, i.e. uses on domain projected skeleton rule
      auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
          mField, elementVolumeName, SPACE_DIM, Sev::noisy);
      op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>();
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
          materialH1Positions, frontAdjEdges);
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
      boost::shared_ptr<double> piola_scale_ptr(new double);
      op_loop_domain_side->getOpPtrVector().push_back(
          physicalEquations->returnOpSetScale(piola_scale_ptr,
                                              physicalEquations));
 
      auto piola_stress_mat_ptr = boost::make_shared<MatrixDouble>();
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<3, 3>(piolaStress,
                                               piola_stress_mat_ptr));
      pip.push_back(op_loop_domain_side);
      return std::make_tuple(broken_data_ptr, piola_scale_ptr,
                             piola_stress_mat_ptr);
    };
 
    auto set_rhs = [&]() {
      MoFEMFunctionBegin;
 
      auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
          get_broken_op_side(fe_rhs->getOpPtrVector());
 
      fe_rhs->getOpPtrVector().push_back(
          new OpDispBc(broken_data_ptr, bcSpatialDispVecPtr, timeScaleMap));
      fe_rhs->getOpPtrVector().push_back(new OpAnalyticalDispBc(
          broken_data_ptr, bcSpatialAnalyticalDisplacementVecPtr,
          timeScaleMap));
      fe_rhs->getOpPtrVector().push_back(new OpRotationBc(
          broken_data_ptr, bcSpatialRotationVecPtr, timeScaleMap));
 
      fe_rhs->getOpPtrVector().push_back(
          new OpBrokenTractionBc(hybridSpatialDisp, bcSpatialTractionVecPtr,
                                 piola_scale_ptr, timeScaleMap));
      auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
      // if you push gradient of L2 base to physical element, it will not work.
      fe_rhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM - 1>(
              hybridSpatialDisp, hybrid_grad_ptr));
      fe_rhs->getOpPtrVector().push_back(new OpBrokenPressureBc(
          hybridSpatialDisp, bcSpatialPressureVecPtr, piola_scale_ptr,
          hybrid_grad_ptr, timeScaleMap));
      fe_rhs->getOpPtrVector().push_back(new OpBrokenAnalyticalTractionBc(
          hybridSpatialDisp, bcSpatialAnalyticalTractionVecPtr, piola_scale_ptr,
          timeScaleMap));
 
      auto hybrid_ptr = boost::make_shared<MatrixDouble>();
      fe_rhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(hybridSpatialDisp,
                                                      hybrid_ptr));
      fe_rhs->getOpPtrVector().push_back(new OpNormalDispRhsBc(
          hybridSpatialDisp, hybrid_ptr, piola_stress_mat_ptr,
          bcSpatialNormalDisplacementVecPtr, timeScaleMap));
 
      auto get_normal_disp_bc_faces = [&]() {
        auto faces =
            get_range_from_block(mField, "NORMAL_DISPLACEMENT", SPACE_DIM - 1);
        return boost::make_shared<Range>(faces);
      };
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using BdyEleOp = BoundaryEle::UserDataOperator;
      using OpC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
          GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
      fe_rhs->getOpPtrVector().push_back(new OpC_dBroken(
          broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0),
          get_normal_disp_bc_faces()));
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_lhs = [&]() {
      MoFEMFunctionBegin;
 
      auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
          get_broken_op_side(fe_lhs->getOpPtrVector());
 
      fe_lhs->getOpPtrVector().push_back(new OpNormalDispLhsBc_dU(
          hybridSpatialDisp, bcSpatialNormalDisplacementVecPtr, timeScaleMap));
      fe_lhs->getOpPtrVector().push_back(new OpNormalDispLhsBc_dP(
          hybridSpatialDisp, broken_data_ptr, bcSpatialNormalDisplacementVecPtr,
          timeScaleMap));
 
      auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
      // if you push gradient of L2 base to physical element, it will not work.
      fe_rhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM - 1>(
              hybridSpatialDisp, hybrid_grad_ptr));
      fe_lhs->getOpPtrVector().push_back(new OpBrokenPressureBcLhs_dU(
          hybridSpatialDisp, bcSpatialPressureVecPtr, hybrid_grad_ptr,
          timeScaleMap));
 
      auto get_normal_disp_bc_faces = [&]() {
        auto faces =
            get_range_from_block(mField, "NORMAL_DISPLACEMENT", SPACE_DIM - 1);
        return boost::make_shared<Range>(faces);
      };
      
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using BdyEleOp = BoundaryEle::UserDataOperator;
      using OpC = FormsIntegrators<BdyEleOp>::Assembly<A>::BiLinearForm<
          GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
      fe_lhs->getOpPtrVector().push_back(new OpC(
          hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0),
          true, true, get_normal_disp_bc_faces()));
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_rhs();
    CHKERR set_lhs();
  }
 
  MoFEMFunctionReturn(0);
}

setNewFrontCoordinates()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setNewFrontCoordinates

(

)

Definition at line 5902 of file EshelbianPlasticity.cpp.

                                                     {
  MoFEMFunctionBegin;
 
  if (!maxMovedFaces)
    MoFEMFunctionReturnHot(0);
 
  Tag th_front_position;
  auto rval =
      mField.get_moab().tag_get_handle("FrontPosition", th_front_position);
  if (rval == MB_SUCCESS && maxMovedFaces) {
    Range verts;
    CHKERR mField.get_moab().get_connectivity(*maxMovedFaces, verts, true);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
    std::vector<double> coords(3 * verts.size());
    CHKERR mField.get_moab().get_coords(verts, coords.data());
    std::vector<double> pos(3 * verts.size());
    CHKERR mField.get_moab().tag_get_data(th_front_position, verts, pos.data());
    for (int i = 0; i != 3 * verts.size(); ++i) {
      coords[i] += pos[i];
    }
    CHKERR mField.get_moab().set_coords(verts, coords.data());
    double zero[] = {0., 0., 0.};
    CHKERR mField.get_moab().tag_clear_data(th_front_position, verts, zero);
  }
 
#ifndef NDEBUG
  constexpr bool debug = false;
  if (debug) {
 
    CHKERR save_range(
        mField.get_moab(),
        "set_coords_faces_" +
            boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
        *maxMovedFaces);
  }
#endif
  MoFEMFunctionReturn(0);
}

setVolumeElementOps()

MoFEMErrorCode EshelbianCore::setVolumeElementOps	(	const int	tag,
		const bool	add_elastic,
		const bool	add_material,
		boost::shared_ptr< VolumeElementForcesAndSourcesCore > &	fe_rhs,
		boost::shared_ptr< VolumeElementForcesAndSourcesCore > &	fe_lhs )

Contact requires that body is marked

Definition at line 2438 of file EshelbianPlasticity.cpp.

                                                                {
  MoFEMFunctionBegin;
 
  /** Contact requires that body is marked */
  auto get_body_range = [this](auto name, int dim) {
    std::map<int, Range> map;
 
    for (auto m_ptr :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % name).str()
 
                 ))
 
    ) {
      Range ents;
      CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
                                                            dim, ents, true),
                     "by dim");
      map[m_ptr->getMeshsetId()] = ents;
    }
 
    return map;
  };
 
  auto rule_contact = [](int, int, int o) { return -1; };
  auto refine = Tools::refineTriangle(contactRefinementLevels);
 
  auto set_rule_contact = [refine](
 
                              ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                              int order_col, int order_data
 
                          ) {
    MoFEMFunctionBegin;
    auto rule = 2 * order_data;
    fe_raw_ptr->gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
    MoFEMFunctionReturn(0);
  };
 
  // Right hand side
  fe_rhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
  CHKERR setBaseVolumeElementOps(tag, true, false, true, SmartPetscObj<Vec>(),
                                 fe_rhs);
 
  // elastic
  if (add_elastic) {
 
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialEquilibrium(spatialL2Disp, dataAtPts, alphaW, alphaRho));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialRotation(rotAxis, dataAtPts, alphaOmega));
    if (noStretch) {
      // do nothing - no stretch approximation
    } else {
      if (!internalStressTagName.empty()) {
        switch (internalStressInterpOrder) {
        case 0:
          fe_rhs->getOpPtrVector().push_back(
              new OpGetInternalStress<0>(dataAtPts, internalStressTagName));
          break;
        case 1:
          fe_rhs->getOpPtrVector().push_back(
              new OpGetInternalStress<1>(dataAtPts, internalStressTagName));
          break;
        default:
          SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
                   "Unsupported internal stress interpolation order %d",
                   internalStressInterpOrder);
        }
        if (internalStressVoigt) {
          fe_rhs->getOpPtrVector().push_back(
              new OpSpatialPhysicalInternalStress<true>(stretchTensor,
                                                        dataAtPts));
        } else {
          fe_rhs->getOpPtrVector().push_back(
              new OpSpatialPhysicalInternalStress<false>(stretchTensor,
                                                         dataAtPts));
        }
      }
      if (auto op = physicalEquations->returnOpSpatialPhysicalExternalStrain(
              stretchTensor, dataAtPts, externalStrainVecPtr, timeScaleMap)) {
        fe_rhs->getOpPtrVector().push_back(op);
      } else if (externalStrainVecPtr && !externalStrainVecPtr->empty()) {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
                "OpSpatialPhysicalExternalStrain not implemented for this "
                "material");
      }
 
      fe_rhs->getOpPtrVector().push_back(
          physicalEquations->returnOpSpatialPhysical(stretchTensor, dataAtPts,
                                                     alphaU));
    }
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyP(piolaStress, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
 
    auto set_hybridisation = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      using BdyEleOp = BoundaryEle::UserDataOperator;
 
      // First: Iterate over skeleton FEs adjacent to Domain FEs
      // Note:  BoundaryEle, i.e. uses skeleton interation rule
      auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
          mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
      // op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
      op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
        return -1;
      };
      op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
          SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
 
      CHKERR EshelbianPlasticity::
          AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              op_loop_skeleton_side->getOpPtrVector(), {L2},
              materialH1Positions, frontAdjEdges);
 
      // Second: Iterate over domain FEs adjacent to skelton, particularly one
      // domain element.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
      // Note: EleOnSide, i.e. uses on domain projected skeleton rule
      auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
          mField, elementVolumeName, SPACE_DIM, Sev::noisy);
      op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>();
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
          materialH1Positions, frontAdjEdges);
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
      auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(piolaStress,
                                                               flux_mat_ptr));
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
 
      // Assemble on skeleton
      op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
      using OpC_dHybrid = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
          GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
      using OpC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
          GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
      op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dHybrid(
          hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0)));
      auto hybrid_ptr = boost::make_shared<MatrixDouble>();
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(hybridSpatialDisp,
                                                      hybrid_ptr));
      op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dBroken(
          broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0)));
 
      // Add skeleton to domain pipeline
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_contact = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      using BdyEleOp = BoundaryEle::UserDataOperator;
 
      // First: Iterate over skeleton FEs adjacent to Domain FEs
      // Note:  BoundaryEle, i.e. uses skeleton interation rule
      auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
          mField, contactElement, SPACE_DIM - 1, Sev::noisy);
 
      op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
      op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
      CHKERR EshelbianPlasticity::
          AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              op_loop_skeleton_side->getOpPtrVector(), {L2},
              materialH1Positions, frontAdjEdges);
 
      // Second: Iterate over domain FEs adjacent to skelton, particularly
      // one domain element.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
 
      // Data storing contact fields
      auto contact_common_data_ptr =
          boost::make_shared<ContactOps::CommonData>();
 
      auto add_ops_domain_side = [&](auto &pip) {
        MoFEMFunctionBegin;
        // Note: EleOnSide, i.e. uses on domain projected skeleton rule
        auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
            mField, elementVolumeName, SPACE_DIM, Sev::noisy);
        op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
            boost::make_shared<CGGUserPolynomialBase>();
        CHKERR
        EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
            materialH1Positions, frontAdjEdges);
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpGetBrokenBaseSideData<SideEleOp>(piolaStress,
                                                   broken_data_ptr));
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
                piolaStress, contact_common_data_ptr->contactTractionPtr()));
        pip.push_back(op_loop_domain_side);
        MoFEMFunctionReturn(0);
      };
 
      auto add_ops_contact_rhs = [&](auto &pip) {
        MoFEMFunctionBegin;
        // get body id and SDF range
        auto contact_sfd_map_range_ptr =
            boost::make_shared<std::map<int, Range>>(
                get_body_range("CONTACT_SDF", SPACE_DIM - 1));
 
        pip.push_back(new OpCalculateVectorFieldValues<3>(
            contactDisp, contact_common_data_ptr->contactDispPtr()));
        auto u_h1_ptr = boost::make_shared<MatrixDouble>();
        pip.push_back(
            new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
        pip.push_back(new OpTreeSearch(
            contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
            get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
            nullptr));
        pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Rhs(
            contactDisp, contact_common_data_ptr, contactTreeRhs,
            contact_sfd_map_range_ptr));
        pip.push_back(
            new EshelbianPlasticity::OpConstrainBoundaryHDivRhs<A, GAUSS>(
                broken_data_ptr, contact_common_data_ptr, contactTreeRhs));
 
        MoFEMFunctionReturn(0);
      };
 
      // push ops to face/side pipeline
      CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
      CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
 
      // Add skeleton to domain pipeline
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
    CHKERR set_contact(fe_rhs->getOpPtrVector());
 
    // Body forces
    using BodyNaturalBC =
        NaturalBC<VolumeElementForcesAndSourcesCore::UserDataOperator>::
            Assembly<PETSC>::LinearForm<GAUSS>;
    using OpBodyForce =
        BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
 
    auto body_time_scale =
        boost::make_shared<DynamicRelaxationTimeScale>("body_force.txt");
    CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
        fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {body_time_scale},
        "BODY_FORCE", Sev::inform);
  }
 
  // Left hand side
  fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
  CHKERR setBaseVolumeElementOps(tag, true, true, true, SmartPetscObj<Vec>(),
                                 fe_lhs);
 
  // elastic
  if (add_elastic) {
 
    if (noStretch) {
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialConsistency_dP_dP(piolaStress, piolaStress, dataAtPts));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dP(
          bubbleField, piolaStress, dataAtPts));
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialConsistency_dBubble_dBubble(bubbleField, bubbleField,
                                                   dataAtPts));
    } else {
      fe_lhs->getOpPtrVector().push_back(
          physicalEquations->returnOpSpatialPhysical_du_du(
              stretchTensor, stretchTensor, dataAtPts, alphaU));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
          stretchTensor, piolaStress, dataAtPts, true));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
          stretchTensor, bubbleField, dataAtPts, true));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
          stretchTensor, rotAxis, dataAtPts,
          symmetrySelector == SYMMETRIC ? true : false));
    }
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
        spatialL2Disp, piolaStress, dataAtPts, true));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
        spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
        piolaStress, rotAxis, dataAtPts,
        symmetrySelector == SYMMETRIC ? true : false));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
        bubbleField, rotAxis, dataAtPts,
        symmetrySelector == SYMMETRIC ? true : false));
 
    if (symmetrySelector > SYMMETRIC) {
      if (!noStretch) {
        fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_du(
            rotAxis, stretchTensor, dataAtPts, false));
      }
      fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
          rotAxis, piolaStress, dataAtPts, false));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
          rotAxis, bubbleField, dataAtPts, false));
    }
    fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_domega(
        rotAxis, rotAxis, dataAtPts, alphaOmega));
    
    // stabilise 
    fe_lhs->getOpPtrVector().push_back(new OpAssembleVolumeStabilize());
 
 
    auto set_hybridisation = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      using BdyEleOp = BoundaryEle::UserDataOperator;
 
      // First: Iterate over skeleton FEs adjacent to Domain FEs
      // Note:  BoundaryEle, i.e. uses skeleton interation rule
      auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
          mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
      // op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
      op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
        return -1;
      };
      op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
          SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
      CHKERR EshelbianPlasticity::
          AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              op_loop_skeleton_side->getOpPtrVector(), {L2},
              materialH1Positions, frontAdjEdges);
 
      // Second: Iterate over domain FEs adjacent to skelton, particularly one
      // domain element.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
      // Note: EleOnSide, i.e. uses on domain projected skeleton rule
      auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
          mField, elementVolumeName, SPACE_DIM, Sev::noisy);
      op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>();
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
          materialH1Positions, frontAdjEdges);
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
 
      op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
      using OpC = FormsIntegrators<BdyEleOp>::Assembly<A>::BiLinearForm<
          GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
      op_loop_skeleton_side->getOpPtrVector().push_back(
          new OpC(hybridSpatialDisp, broken_data_ptr,
                  boost::make_shared<double>(1.0), true, false));
 
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_contact = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      using BoundaryEle =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
      using EleOnSide =
          PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
      using SideEleOp = EleOnSide::UserDataOperator;
      using BdyEleOp = BoundaryEle::UserDataOperator;
 
      // First: Iterate over skeleton FEs adjacent to Domain FEs
      // Note:  BoundaryEle, i.e. uses skeleton interation rule
      auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
          mField, contactElement, SPACE_DIM - 1, Sev::noisy);
 
      op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
      op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
      CHKERR EshelbianPlasticity::
          AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              op_loop_skeleton_side->getOpPtrVector(), {L2},
              materialH1Positions, frontAdjEdges);
 
      // Second: Iterate over domain FEs adjacent to skelton, particularly
      // one domain element.
      auto broken_data_ptr =
          boost::make_shared<std::vector<BrokenBaseSideData>>();
 
      // Data storing contact fields
      auto contact_common_data_ptr =
          boost::make_shared<ContactOps::CommonData>();
 
      auto add_ops_domain_side = [&](auto &pip) {
        MoFEMFunctionBegin;
        // Note: EleOnSide, i.e. uses on domain projected skeleton rule
        auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
            mField, elementVolumeName, SPACE_DIM, Sev::noisy);
        op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
            boost::make_shared<CGGUserPolynomialBase>();
        CHKERR
        EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
            materialH1Positions, frontAdjEdges);
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpGetBrokenBaseSideData<SideEleOp>(piolaStress,
                                                   broken_data_ptr));
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
                piolaStress, contact_common_data_ptr->contactTractionPtr()));
        pip.push_back(op_loop_domain_side);
        MoFEMFunctionReturn(0);
      };
 
      auto add_ops_contact_lhs = [&](auto &pip) {
        MoFEMFunctionBegin;
        pip.push_back(new OpCalculateVectorFieldValues<3>(
            contactDisp, contact_common_data_ptr->contactDispPtr()));
        auto u_h1_ptr = boost::make_shared<MatrixDouble>();
        pip.push_back(
            new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
        pip.push_back(new OpTreeSearch(
            contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
            get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
            nullptr));
 
        // get body id and SDF range
        auto contact_sfd_map_range_ptr =
            boost::make_shared<std::map<int, Range>>(
                get_body_range("CONTACT_SDF", SPACE_DIM - 1));
 
        pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU(
            contactDisp, contactDisp, contact_common_data_ptr, contactTreeRhs,
            contact_sfd_map_range_ptr));
        pip.push_back(
            new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP<A, GAUSS>(
                contactDisp, broken_data_ptr, contact_common_data_ptr,
                contactTreeRhs, contact_sfd_map_range_ptr));
        pip.push_back(
            new EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU<A, GAUSS>(
                broken_data_ptr, contactDisp, contact_common_data_ptr,
                contactTreeRhs));
 
        MoFEMFunctionReturn(0);
      };
 
      // push ops to face/side pipeline
      CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
      CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
 
      // Add skeleton to domain pipeline
      pip.push_back(op_loop_skeleton_side);
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
    CHKERR set_contact(fe_lhs->getOpPtrVector());
  }
 
  if (add_material) {
  }
 
  MoFEMFunctionReturn(0);
}

solveDynamicRelaxation()

MoFEMErrorCode EshelbianCore::solveDynamicRelaxation	(	TS	ts,
		Vec	x )

Definition at line 3571 of file EshelbianPlasticity.cpp.

                                                                 {
  MoFEMFunctionBegin;
 
  auto storage = solve_elastic_setup::setup(this, ts, x, false);
 
  double final_time = 1;
  double delta_time = 0.1;
  int max_it = 10;
  PetscBool ts_h1_update = PETSC_FALSE;
 
  PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options",
                           "none");
 
  CHKERR PetscOptionsScalar("-dynamic_final_time",
                            "dynamic relaxation final time", "", final_time,
                            &final_time, PETSC_NULLPTR);
  CHKERR PetscOptionsScalar("-dynamic_delta_time",
                            "dynamic relaxation final time", "", delta_time,
                            &delta_time, PETSC_NULLPTR);
  CHKERR PetscOptionsInt("-dynamic_max_it", "dynamic relaxation iterations", "",
                         max_it, &max_it, PETSC_NULLPTR);
  CHKERR PetscOptionsBool("-dynamic_h1_update", "update each ts step", "",
                          ts_h1_update, &ts_h1_update, PETSC_NULLPTR);
 
  PetscOptionsEnd();
 
  auto setup_ts_monitor = [&]() {
    auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
    return monitor_ptr;
  };
  auto monitor_ptr = setup_ts_monitor();
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
  CHKERR TSSetUp(ts);
  CHKERR TSElasticPostStep::postStepInitialise(this);
 
  double ts_delta_time;
  CHKERR TSGetTimeStep(ts, &ts_delta_time);
 
  if (ts_h1_update) {
    CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
    CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  }
 
  CHKERR TSElasticPostStep::preStepFun(ts);
  CHKERR TSElasticPostStep::postStepFun(ts);
 
  dynamicTime = 0;
  dynamicStep = 0;
  monitor_ptr->ts_u = PETSC_NULLPTR;
  monitor_ptr->ts_t = dynamicTime;
  monitor_ptr->ts_step = dynamicStep;
  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
  MOFEM_LOG("EP", Sev::inform) << "IT " << dynamicStep << " Time "
                               << dynamicTime << " delta time " << delta_time;
  dynamicTime += delta_time;
  ++dynamicStep;
 
  for (; dynamicTime < final_time; dynamicTime += delta_time) {
    MOFEM_LOG("EP", Sev::inform) << "IT " << dynamicStep << " Time "
                                 << dynamicTime << " delta time " << delta_time;
 
    CHKERR TSSetStepNumber(ts, 0);
    CHKERR TSSetTime(ts, 0);
    CHKERR TSSetTimeStep(ts, ts_delta_time);
    if (!ts_h1_update) {
      CHKERR TSElasticPostStep::preStepFun(ts);
    }
    CHKERR TSSolve(ts, PETSC_NULLPTR);
    if (!ts_h1_update) {
      CHKERR TSElasticPostStep::postStepFun(ts);
    }
 
    monitor_ptr->ts_u = PETSC_NULLPTR;
    monitor_ptr->ts_t = dynamicTime;
    monitor_ptr->ts_step = dynamicStep;
    CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
 
    ++dynamicStep;
    if (dynamicStep > max_it)
      break;
  }
 
  CHKERR TSElasticPostStep::postStepDestroy();
  TetPolynomialBase::switchCacheBaseOff<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
 
  MoFEMFunctionReturn(0);
}

solveElastic()

MoFEMErrorCode EshelbianCore::solveElastic	(	TS	ts,
		Vec	x )

Definition at line 3441 of file EshelbianPlasticity.cpp.

                                                       {
  MoFEMFunctionBegin;
 
  PetscBool debug_model = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-debug_model", &debug_model,
                             PETSC_NULLPTR);
 
  if (debug_model == PETSC_TRUE) {
    auto ts_ctx_ptr = getDMTsCtx(dmElastic);
    auto post_proc = [&](TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, Vec F,
                         void *ctx) {
      MoFEMFunctionBegin;
 
      SNES snes;
      CHKERR TSGetSNES(ts, &snes);
      int it;
      CHKERR SNESGetIterationNumber(snes, &it);
      std::string file_name = "snes_iteration_" + std::to_string(it) + ".h5m";
      CHKERR postProcessResults(1, file_name, F, u_t);
      std::string file_skel_name =
          "snes_iteration_skel_" + std::to_string(it) + ".h5m";
 
      auto get_material_force_tag = [&]() {
        auto &moab = mField.get_moab();
        Tag tag;
        CHK_MOAB_THROW(moab.tag_get_handle("MaterialForce", tag),
                       "can't get tag");
        return tag;
      };
 
      CHKERR calculateFaceMaterialForce(1, ts);
      CHKERR postProcessSkeletonResults(1, file_skel_name, F,
                                        {get_material_force_tag()});
 
      MoFEMFunctionReturn(0);
    };
    ts_ctx_ptr->tsDebugHook = post_proc;
  }
 
  auto storage = solve_elastic_setup::setup(this, ts, x, true);
 
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
    Vec xx;
    CHKERR VecDuplicate(x, &xx);
    CHKERR VecZeroEntries(xx);
    CHKERR TS2SetSolution(ts, x, xx);
    CHKERR VecDestroy(&xx);
  } else {
    CHKERR TSSetSolution(ts, x);
  }
 
  TetPolynomialBase::switchCacheBaseOn<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
  CHKERR TSSetUp(ts);
  CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
  CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  CHKERR TSElasticPostStep::postStepInitialise(this);
  CHKERR TSSolve(ts, PETSC_NULLPTR);
  CHKERR TSElasticPostStep::postStepDestroy();
  TetPolynomialBase::switchCacheBaseOff<HDIV>(
      {elasticFeLhs.get(), elasticFeRhs.get()});
 
  SNES snes;
  CHKERR TSGetSNES(ts, &snes);
  int lin_solver_iterations;
  CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
  MOFEM_LOG("EP", Sev::inform)
      << "Number of linear solver iterations " << lin_solver_iterations;
 
  PetscBool test_cook_flg = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_cook", &test_cook_flg,
                             PETSC_NULLPTR);
  if (test_cook_flg) {
    constexpr int expected_lin_solver_iterations = 11;
    if (lin_solver_iterations > expected_lin_solver_iterations)
      SETERRQ(
          PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
          "Expected number of iterations is different than expected %d > %d",
          lin_solver_iterations, expected_lin_solver_iterations);
  }
 
  PetscBool test_sslv116_flag = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_sslv116",
                             &test_sslv116_flag, PETSC_NULLPTR);
 
  if (test_sslv116_flag) {
    double max_val = 0.0;
    double min_val = 0.0;
    auto field_min_max = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
      MoFEMFunctionBeginHot;
      auto ent_type = ent_ptr->getEntType();
      if (ent_type == MBVERTEX) {
        max_val = std::max(ent_ptr->getEntFieldData()[SPACE_DIM - 1], max_val);
        min_val = std::min(ent_ptr->getEntFieldData()[SPACE_DIM - 1], min_val);
      }
      MoFEMFunctionReturnHot(0);
    };
    CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
        field_min_max, spatialH1Disp);
 
    double global_max_val = 0.0;
    double global_min_val = 0.0;
    MPI_Allreduce(&max_val, &global_max_val, 1, MPI_DOUBLE, MPI_MAX,
                  mField.get_comm());
    MPI_Allreduce(&min_val, &global_min_val, 1, MPI_DOUBLE, MPI_MIN,
                  mField.get_comm());
    MOFEM_LOG("EP", Sev::inform)
        << "Max " << spatialH1Disp << " value: " << global_max_val;
    MOFEM_LOG("EP", Sev::inform)
        << "Min " << spatialH1Disp << " value: " << global_min_val;
 
    double ref_max_val = 0.00767;
    double ref_min_val = -0.00329;
    if (std::abs(global_max_val - ref_max_val) > 1e-5) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "Incorrect max value of the displacement field: %f != %f",
               global_max_val, ref_max_val);
    }
    if (std::abs(global_min_val - ref_min_val) > 4e-5) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "Incorrect min value of the displacement field: %f != %f",
               global_min_val, ref_min_val);
    }
  }
 
  CHKERR gettingNorms();
 
  MoFEMFunctionReturn(0);
}

Friends And Related Symbol Documentation

solve_elastic_set_up

friend struct solve_elastic_set_up

friend

Definition at line 363 of file EshelbianPlasticity.hpp.

Member Data Documentation

a00FieldList

std::vector<std::string> EshelbianPlasticity::EshelbianCore::a00FieldList

Definition at line 447 of file EshelbianPlasticity.hpp.

a00RangeList

std::vector<boost::shared_ptr<Range> > EshelbianPlasticity::EshelbianCore::a00RangeList

Definition at line 449 of file EshelbianPlasticity.hpp.

addCrackMeshsetId

int EshelbianPlasticity::EshelbianCore::addCrackMeshsetId = 1000

inlinestatic

Definition at line 32 of file EshelbianPlasticity.hpp.

alphaOmega

double EshelbianPlasticity::EshelbianCore::alphaOmega = 0

Definition at line 222 of file EshelbianPlasticity.hpp.

alphaRho

double EshelbianPlasticity::EshelbianCore::alphaRho = 0

Definition at line 223 of file EshelbianPlasticity.hpp.

alphaU

double EshelbianPlasticity::EshelbianCore::alphaU = 0

Definition at line 220 of file EshelbianPlasticity.hpp.

alphaW

double EshelbianPlasticity::EshelbianCore::alphaW = 0

Definition at line 221 of file EshelbianPlasticity.hpp.

AnalyticalExprPythonPtr

boost::shared_ptr<AnalyticalExprPython> EshelbianPlasticity::EshelbianCore::AnalyticalExprPythonPtr

Definition at line 180 of file EshelbianPlasticity.hpp.

aoS

AO EshelbianPlasticity::EshelbianCore::aoS = PETSC_NULLPTR

Definition at line 445 of file EshelbianPlasticity.hpp.

bcSpatialAnalyticalDisplacementVecPtr

boost::shared_ptr<AnalyticalDisplacementBcVec> EshelbianPlasticity::EshelbianCore::bcSpatialAnalyticalDisplacementVecPtr

Definition at line 236 of file EshelbianPlasticity.hpp.

bcSpatialAnalyticalTractionVecPtr

boost::shared_ptr<AnalyticalTractionBcVec> EshelbianPlasticity::EshelbianCore::bcSpatialAnalyticalTractionVecPtr

Definition at line 237 of file EshelbianPlasticity.hpp.

bcSpatialDispVecPtr

boost::shared_ptr<BcDispVec> EshelbianPlasticity::EshelbianCore::bcSpatialDispVecPtr

Definition at line 230 of file EshelbianPlasticity.hpp.

bcSpatialFreeTractionVecPtr

boost::shared_ptr<TractionFreeBc> EshelbianPlasticity::EshelbianCore::bcSpatialFreeTractionVecPtr

Definition at line 233 of file EshelbianPlasticity.hpp.

bcSpatialNormalDisplacementVecPtr

boost::shared_ptr<NormalDisplacementBcVec> EshelbianPlasticity::EshelbianCore::bcSpatialNormalDisplacementVecPtr

Definition at line 234 of file EshelbianPlasticity.hpp.

bcSpatialPressureVecPtr

boost::shared_ptr<PressureBcVec> EshelbianPlasticity::EshelbianCore::bcSpatialPressureVecPtr

Definition at line 238 of file EshelbianPlasticity.hpp.

bcSpatialRotationVecPtr

boost::shared_ptr<BcRotVec> EshelbianPlasticity::EshelbianCore::bcSpatialRotationVecPtr

Definition at line 231 of file EshelbianPlasticity.hpp.

bcSpatialTractionVecPtr

boost::shared_ptr<TractionBcVec> EshelbianPlasticity::EshelbianCore::bcSpatialTractionVecPtr

Definition at line 232 of file EshelbianPlasticity.hpp.

bitAdjEnt

BitRefLevel EshelbianPlasticity::EshelbianCore::bitAdjEnt = BitRefLevel().set()

bit ref level for parent

Definition at line 430 of file EshelbianPlasticity.hpp.

bitAdjEntMask

BitRefLevel EshelbianPlasticity::EshelbianCore::bitAdjEntMask

Initial value:

=
      BitRefLevel().set()

bit ref level for parent parent

Definition at line 431 of file EshelbianPlasticity.hpp.

bitAdjParent

BitRefLevel EshelbianPlasticity::EshelbianCore::bitAdjParent = BitRefLevel().set()

bit ref level for parent

Definition at line 427 of file EshelbianPlasticity.hpp.

bitAdjParentMask

BitRefLevel EshelbianPlasticity::EshelbianCore::bitAdjParentMask

Initial value:

=
      BitRefLevel().set()

bit ref level for parent parent

Definition at line 428 of file EshelbianPlasticity.hpp.

bubbleField

const std::string EshelbianPlasticity::EshelbianCore::bubbleField = "bubble"

Definition at line 204 of file EshelbianPlasticity.hpp.

contactDisp

const std::string EshelbianPlasticity::EshelbianCore::contactDisp = "contactDisp"

Definition at line 201 of file EshelbianPlasticity.hpp.

contactElement

const std::string EshelbianPlasticity::EshelbianCore::contactElement = "CONTACT"

Definition at line 210 of file EshelbianPlasticity.hpp.

contactFaces

boost::shared_ptr<Range> EshelbianPlasticity::EshelbianCore::contactFaces

Definition at line 415 of file EshelbianPlasticity.hpp.

contactRefinementLevels

int EshelbianPlasticity::EshelbianCore::contactRefinementLevels = 1

Definition at line 225 of file EshelbianPlasticity.hpp.

contactTreeRhs

boost::shared_ptr<ContactTree> EshelbianPlasticity::EshelbianCore::contactTreeRhs

Make a contact tree.

Definition at line 186 of file EshelbianPlasticity.hpp.

crackFaces

boost::shared_ptr<Range> EshelbianPlasticity::EshelbianCore::crackFaces

Definition at line 416 of file EshelbianPlasticity.hpp.

crackHybridIs

SmartPetscObj<IS> EshelbianPlasticity::EshelbianCore::crackHybridIs

Definition at line 446 of file EshelbianPlasticity.hpp.

crackingOn

PetscBool EshelbianPlasticity::EshelbianCore::crackingOn = PETSC_FALSE

inlinestatic

Definition at line 23 of file EshelbianPlasticity.hpp.

crackingStartTime

double EshelbianPlasticity::EshelbianCore::crackingStartTime = 0

inlinestatic

Definition at line 24 of file EshelbianPlasticity.hpp.

d_f

boost::function< double(const double)> EshelbianCore::d_f

static

Initial value:

=
    EshelbianCore::d_f_log_e

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 45 of file EshelbianPlasticity.hpp.

dataAtPts

boost::shared_ptr<DataAtIntegrationPts> EshelbianPlasticity::EshelbianCore::dataAtPts

Definition at line 178 of file EshelbianPlasticity.hpp.

dd_f

boost::function< double(const double)> EshelbianCore::dd_f

static

Initial value:

=
    EshelbianCore::dd_f_log_e

Examples

EshelbianPlasticity.cpp.

Definition at line 46 of file EshelbianPlasticity.hpp.

dM

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dM

Coupled problem all fields.

Definition at line 188 of file EshelbianPlasticity.hpp.

dmElastic

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElastic

Elastic problem.

Definition at line 189 of file EshelbianPlasticity.hpp.

dmMaterial

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmMaterial

Material problem.

Definition at line 190 of file EshelbianPlasticity.hpp.

dmPrjSpatial

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmPrjSpatial

Projection spatial displacement.

Definition at line 191 of file EshelbianPlasticity.hpp.

dynamicRelaxation

PetscBool EshelbianPlasticity::EshelbianCore::dynamicRelaxation

inlinestatic

Initial value:

=
      PETSC_FALSE

Examples

EshelbianOperators.cpp.

Definition at line 21 of file EshelbianPlasticity.hpp.

dynamicStep

int EshelbianPlasticity::EshelbianCore::dynamicStep

inlinestatic

Initial value:

=
      0

Examples

EshelbianOperators.cpp.

Definition at line 29 of file EshelbianPlasticity.hpp.

dynamicTime

double EshelbianPlasticity::EshelbianCore::dynamicTime

inlinestatic

Initial value:

=
      0

Examples

EshelbianOperators.cpp.

Definition at line 27 of file EshelbianPlasticity.hpp.

edgeExchange

CommInterface::EntitiesPetscVector EshelbianPlasticity::EshelbianCore::edgeExchange

Definition at line 438 of file EshelbianPlasticity.hpp.

elasticBcLhs

boost::shared_ptr<FaceElementForcesAndSourcesCore> EshelbianPlasticity::EshelbianCore::elasticBcLhs

Definition at line 184 of file EshelbianPlasticity.hpp.

elasticBcRhs

boost::shared_ptr<FaceElementForcesAndSourcesCore> EshelbianPlasticity::EshelbianCore::elasticBcRhs

Definition at line 185 of file EshelbianPlasticity.hpp.

elasticFeLhs

boost::shared_ptr<VolumeElementForcesAndSourcesCore> EshelbianPlasticity::EshelbianCore::elasticFeLhs

Definition at line 183 of file EshelbianPlasticity.hpp.

elasticFeRhs

boost::shared_ptr<VolumeElementForcesAndSourcesCore> EshelbianPlasticity::EshelbianCore::elasticFeRhs

Definition at line 182 of file EshelbianPlasticity.hpp.

elementVolumeName

const std::string EshelbianPlasticity::EshelbianCore::elementVolumeName = "EP"

Definition at line 206 of file EshelbianPlasticity.hpp.

energyReleaseSelector

enum EnergyReleaseSelector EshelbianPlasticity::EshelbianCore::energyReleaseSelector

inlinestatic

Initial value:

=
      GRIFFITH_SKELETON

Examples

EshelbianOperators.cpp.

Definition at line 34 of file EshelbianPlasticity.hpp.

exponentBase

double EshelbianCore::exponentBase = exp(1)

static

Examples

EshelbianPlasticity.cpp.

Definition at line 43 of file EshelbianPlasticity.hpp.

externalStrainVecPtr

boost::shared_ptr<ExternalStrainVec> EshelbianPlasticity::EshelbianCore::externalStrainVecPtr

Definition at line 239 of file EshelbianPlasticity.hpp.

f

boost::function< double(const double)> EshelbianCore::f = EshelbianCore::f_log_e

static

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 44 of file EshelbianPlasticity.hpp.

faceExchange

CommInterface::EntitiesPetscVector EshelbianPlasticity::EshelbianCore::faceExchange

Definition at line 437 of file EshelbianPlasticity.hpp.

frontAdjEdges

boost::shared_ptr<Range> EshelbianPlasticity::EshelbianCore::frontAdjEdges

Definition at line 418 of file EshelbianPlasticity.hpp.

frontEdges

boost::shared_ptr<Range> EshelbianPlasticity::EshelbianCore::frontEdges

Definition at line 417 of file EshelbianPlasticity.hpp.

frontLayers

int EshelbianPlasticity::EshelbianCore::frontLayers = 3

Definition at line 226 of file EshelbianPlasticity.hpp.

frontVertices

boost::shared_ptr<Range> EshelbianPlasticity::EshelbianCore::frontVertices

Definition at line 419 of file EshelbianPlasticity.hpp.

gradApproximator

enum RotSelector EshelbianPlasticity::EshelbianCore::gradApproximator = LARGE_ROT

inlinestatic

Examples

EshelbianOperators.cpp.

Definition at line 17 of file EshelbianPlasticity.hpp.

griffithEnergy

double EshelbianPlasticity::EshelbianCore::griffithEnergy = 1

inlinestatic

Griffith energy.

Definition at line 33 of file EshelbianPlasticity.hpp.

hybridSpatialDisp

const std::string EshelbianPlasticity::EshelbianCore::hybridSpatialDisp = "hybridSpatialDisp"

Definition at line 199 of file EshelbianPlasticity.hpp.

internalStressInterpOrder

int EshelbianPlasticity::EshelbianCore::internalStressInterpOrder

inlinestatic

Initial value:

=
      1

Definition at line 38 of file EshelbianPlasticity.hpp.

internalStressTagName

std::string EshelbianPlasticity::EshelbianCore::internalStressTagName

inlinestatic

Initial value:

=
      ""

Definition at line 36 of file EshelbianPlasticity.hpp.

internalStressVoigt

PetscBool EshelbianPlasticity::EshelbianCore::internalStressVoigt

inlinestatic

Initial value:

=
      PETSC_FALSE

Definition at line 40 of file EshelbianPlasticity.hpp.

inv_d_f

boost::function< double(const double)> EshelbianCore::inv_d_f

static

Initial value:

=
    EshelbianCore::inv_d_f_log_e

Examples

EshelbianPlasticity.cpp.

Definition at line 48 of file EshelbianPlasticity.hpp.

inv_dd_f

boost::function< double(const double)> EshelbianCore::inv_dd_f

static

Initial value:

=
    EshelbianCore::inv_dd_f_log_e

Examples

EshelbianPlasticity.cpp.

Definition at line 49 of file EshelbianPlasticity.hpp.

inv_f

boost::function< double(const double)> EshelbianCore::inv_f

static

Initial value:

=
    EshelbianCore::inv_f_log_e

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 47 of file EshelbianPlasticity.hpp.

l2UserBaseScale

PetscBool EshelbianPlasticity::EshelbianCore::l2UserBaseScale = PETSC_TRUE

inlinestatic

Examples

EshelbianOperators.cpp.

Definition at line 31 of file EshelbianPlasticity.hpp.

listTagsToTransfer

std::vector<Tag> EshelbianPlasticity::EshelbianCore::listTagsToTransfer

list of tags to transfer to postprocessor

Definition at line 442 of file EshelbianPlasticity.hpp.

materialH1Positions

const std::string EshelbianPlasticity::EshelbianCore::materialH1Positions = "XH1"

Definition at line 198 of file EshelbianPlasticity.hpp.

materialOrder

int EshelbianPlasticity::EshelbianCore::materialOrder = 1

Definition at line 219 of file EshelbianPlasticity.hpp.

materialSkeletonFaces

boost::shared_ptr<Range> EshelbianPlasticity::EshelbianCore::materialSkeletonFaces

Definition at line 421 of file EshelbianPlasticity.hpp.

maxMovedFaces

boost::shared_ptr<Range> EshelbianPlasticity::EshelbianCore::maxMovedFaces

Definition at line 422 of file EshelbianPlasticity.hpp.

mField

MoFEM::Interface& EshelbianPlasticity::EshelbianCore::mField

Definition at line 176 of file EshelbianPlasticity.hpp.

naturalBcElement

const std::string EshelbianPlasticity::EshelbianCore::naturalBcElement = "NATURAL_BC"

Definition at line 207 of file EshelbianPlasticity.hpp.

nbCrackFaces

int EshelbianPlasticity::EshelbianCore::nbCrackFaces = 0

Definition at line 451 of file EshelbianPlasticity.hpp.

nbJIntegralLevels

int EshelbianPlasticity::EshelbianCore::nbJIntegralLevels

inlinestatic

Initial value:

=
      0

Definition at line 25 of file EshelbianPlasticity.hpp.

noStretch

PetscBool EshelbianPlasticity::EshelbianCore::noStretch = PETSC_FALSE

inlinestatic

Definition at line 19 of file EshelbianPlasticity.hpp.

parentAdjSkeletonFunctionDim2

boost::shared_ptr<ParentFiniteElementAdjacencyFunctionSkeleton<2> > EshelbianPlasticity::EshelbianCore::parentAdjSkeletonFunctionDim2

Definition at line 425 of file EshelbianPlasticity.hpp.

physicalEquations

boost::shared_ptr<PhysicalEquations> EshelbianPlasticity::EshelbianCore::physicalEquations

Definition at line 179 of file EshelbianPlasticity.hpp.

piolaStress

const std::string EshelbianPlasticity::EshelbianCore::piolaStress = "P"

Definition at line 193 of file EshelbianPlasticity.hpp.

rotAxis

const std::string EshelbianPlasticity::EshelbianCore::rotAxis = "omega"

Definition at line 203 of file EshelbianPlasticity.hpp.

rotSelector

enum RotSelector EshelbianPlasticity::EshelbianCore::rotSelector = LARGE_ROT

inlinestatic

Examples

EshelbianOperators.cpp.

Definition at line 16 of file EshelbianPlasticity.hpp.

S

Mat EshelbianPlasticity::EshelbianCore::S = PETSC_NULLPTR

Definition at line 444 of file EshelbianPlasticity.hpp.

setSingularity

PetscBool EshelbianPlasticity::EshelbianCore::setSingularity = PETSC_TRUE

inlinestatic

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 20 of file EshelbianPlasticity.hpp.

skeletonElement

const std::string EshelbianPlasticity::EshelbianCore::skeletonElement = "SKELETON"

Definition at line 209 of file EshelbianPlasticity.hpp.

skeletonFaces

boost::shared_ptr<Range> EshelbianPlasticity::EshelbianCore::skeletonFaces

Definition at line 420 of file EshelbianPlasticity.hpp.

skinElement

const std::string EshelbianPlasticity::EshelbianCore::skinElement = "SKIN"

Definition at line 208 of file EshelbianPlasticity.hpp.

solTSStep

SmartPetscObj<Vec> EshelbianPlasticity::EshelbianCore::solTSStep

Definition at line 434 of file EshelbianPlasticity.hpp.

spaceH1Order

int EshelbianPlasticity::EshelbianCore::spaceH1Order = -1

Definition at line 218 of file EshelbianPlasticity.hpp.

spaceOrder

int EshelbianPlasticity::EshelbianCore::spaceOrder = 2

Definition at line 217 of file EshelbianPlasticity.hpp.

spatialH1Disp

const std::string EshelbianPlasticity::EshelbianCore::spatialH1Disp = "wH1"

Definition at line 197 of file EshelbianPlasticity.hpp.

spatialL2Disp

const std::string EshelbianPlasticity::EshelbianCore::spatialL2Disp = "wL2"

Definition at line 195 of file EshelbianPlasticity.hpp.

stretchSelector

enum StretchSelector EshelbianPlasticity::EshelbianCore::stretchSelector = LOG

inlinestatic

Examples

EshelbianOperators.cpp.

Definition at line 18 of file EshelbianPlasticity.hpp.

stretchTensor

const std::string EshelbianPlasticity::EshelbianCore::stretchTensor = "u"

Definition at line 202 of file EshelbianPlasticity.hpp.

symmetrySelector

enum SymmetrySelector EshelbianPlasticity::EshelbianCore::symmetrySelector = NOT_SYMMETRIC

inlinestatic

Examples

EshelbianOperators.cpp.

Definition at line 15 of file EshelbianPlasticity.hpp.

timeScaleMap

std::map<std::string, boost::shared_ptr<ScalingMethod> > EshelbianPlasticity::EshelbianCore::timeScaleMap

Definition at line 241 of file EshelbianPlasticity.hpp.

use_quadratic_exp

bool EshelbianPlasticity::EshelbianCore::use_quadratic_exp = true

inlinestaticconstexpr

Definition at line 51 of file EshelbianPlasticity.hpp.

v_max

double EshelbianPlasticity::EshelbianCore::v_max = 12

inlinestaticconstexpr

Definition at line 52 of file EshelbianPlasticity.hpp.

v_min

double EshelbianPlasticity::EshelbianCore::v_min = -v_max

inlinestaticconstexpr

Definition at line 53 of file EshelbianPlasticity.hpp.

vertexExchange

CommInterface::EntitiesPetscVector EshelbianPlasticity::EshelbianCore::vertexExchange

Definition at line 439 of file EshelbianPlasticity.hpp.

volumeExchange

CommInterface::EntitiesPetscVector EshelbianPlasticity::EshelbianCore::volumeExchange

Definition at line 436 of file EshelbianPlasticity.hpp.

---

The documentation for this struct was generated from the following files:

• users_modules/eshelbian_plasticity/src/EshelbianPlasticity.hpp
• users_modules/eshelbian_plasticity/src/impl/EshelbianADOL-C.cpp
• users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp