EshelbianCore Struct Reference

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

Inheritance diagram for EshelbianCore:

Collaboration diagram for 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

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 12 of file EshelbianCore.hpp.

Constructor & Destructor Documentation

EshelbianCore()

EshelbianCore::EshelbianCore

(

MoFEM::Interface &

m_field

)

Examples

EshelbianPlasticity.cpp.

Definition at line 916 of file EshelbianPlasticity.cpp.

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

~EshelbianCore()

EshelbianCore::~EshelbianCore ( )

virtualdefault

Examples

EshelbianPlasticity.cpp.

Member Function Documentation

addBoundaryFiniteElement()

MoFEMErrorCode EshelbianCore::addBoundaryFiniteElement

(

const EntityHandle

meshset = 0

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 1583 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

)

Examples

EshelbianPlasticity.cpp.

Definition at line 5998 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
	)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 1783 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

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 1157 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
	)

Examples

ep.cpp.

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
	)

Examples

ep.cpp.

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
	)

Examples

ep.cpp.

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
	)

Examples

ep.cpp.

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

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 1520 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

)

Examples

EshelbianPlasticity.cpp.

Definition at line 7092 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

Examples

EshelbianPlasticity.cpp.

Definition at line 4329 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);
 
      // hybrid disp, evalated on face first
      EshelbianPlasticity::AddHOOps<2, 2, 3>::add(
          fe_ptr->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
      fe_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          hybridSpatialDisp, dataAtPts->getHybridDispAtPts()));
      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);
      fe_ptr->getOpPtrVector().push_back(op_loop_domain_side);
      fe_ptr->getOpPtrVector().push_back(new OpFaceMaterialForce(dataAtPts));
 
      // evaluated in side domain, that is op_loop_domain_side
      op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>();
 
      auto inv_jac = boost::make_shared<MatrixDouble>();
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
          materialH1Positions, frontAdjEdges, nullptr, nullptr, inv_jac);
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpSetHOInvJacVectorBase(HDIV, inv_jac));
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
              piolaStress, dataAtPts->getApproxPAtPts()));
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateHTensorTensorField<SPACE_DIM, SPACE_DIM>(
              bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>(
              rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
      if (noStretch) {
        op_loop_domain_side->getOpPtrVector().push_back(
            physicalEquations->returnOpCalculateStretchFromStress(
                dataAtPts, physicalEquations));
      } else {
        op_loop_domain_side->getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
                stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
      }
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpCalculateHVecVectorHessian<SPACE_DIM, SPACE_DIM>(
              hybridSpatialDisp, dataAtPts->getGradPAtPts()()));
      op_loop_domain_side->getOpPtrVector().push_back(
          new OpFaceSideMaterialForce(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();
                }
 
                // 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);
                }
 
                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.};
                  FTensor::Tensor1<double, SPACE_DIM> t_node_force_tmp;
                  t_node_force_tmp(I) = t_node_force(I);
                  t_node_force_tmp.normalize();
                  t_area_dir(I) = -t_node_force_tmp(I);
                  t_area_dir(I) *= l / 2;
                  return t_area_dir;
                }
 
                // 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 adj_faces = get_adj_range(adj_tets, 2);
            auto t_area_dir = calculate_crack_area_growth_direction(
                n, t_node_force);
            if (nbJIntegralLevels) {
              auto adj_edges = intersect(get_adj_range(adj_tets, 1),
                                        unite(*frontEdges, body_edges));
              auto seed_n = get_conn_range(adj_edges);
              auto skin_adj_edges = get_skin(mField, adj_edges);
              skin_adj_edges = subtract(skin_adj_edges, body_skin_conn);
              seed_n = subtract(seed_n, skin_adj_edges);
              t_area_dir(I) = 0;
              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);
              }
              for (auto sn : skin_adj_edges) {
                auto t_area_dir_sn = calculate_crack_area_growth_direction(
                    sn, t_node_force);
                t_area_dir(I) += t_area_dir_sn(I) / 2;
              }
            }
 
            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");
        }
 
        // iterate over edges, and calculate average edge material force
        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);
        };
 
        // iterate over edges, and calculate average edge griffith energy
        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 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

Examples

EshelbianPlasticity.cpp.

Definition at line 4929 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

(

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 6755 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

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 7060 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 EshelbianCore::d_f_linear ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 159 of file EshelbianCore.hpp.

{ return 1; }

d_f_log()

static double EshelbianCore::d_f_log ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 140 of file EshelbianCore.hpp.

                                        {
    return pow(exponentBase, v) * log(EshelbianCore::exponentBase);
  }

d_f_log_e()

static double EshelbianCore::d_f_log_e ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 95 of file EshelbianCore.hpp.

                                          {
    if constexpr (use_quadratic_exp) {
      return d_f_log_e_quadratic(v);
    } else {
      if (v > v_max)
        return std::exp(v_max);
      else 
        return std::exp(v);
    }
  }

d_f_log_e_quadratic()

static double EshelbianCore::d_f_log_e_quadratic ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 64 of file EshelbianCore.hpp.

                                                    {
    if (v > v_max) {
      double e = static_cast<double>(std::exp(v_max));
      double dv = v - v_max;
      return e * dv + e;
    } else {
      return static_cast<double>(std::exp(v));
    } 
  }

dd_f_linear()

static double EshelbianCore::dd_f_linear ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 160 of file EshelbianCore.hpp.

{ return 0; }

dd_f_log()

static double EshelbianCore::dd_f_log ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 143 of file EshelbianCore.hpp.

                                         {
    return pow(EshelbianCore::exponentBase, v) *
           log(EshelbianCore::exponentBase) * log(EshelbianCore::exponentBase);
  }

dd_f_log_e()

static double EshelbianCore::dd_f_log_e ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 105 of file EshelbianCore.hpp.

                                           {
    if constexpr (use_quadratic_exp) {
      return dd_f_log_e_quadratic(v);
    } else {
      if (v > v_max)
        return 0.;
      else 
        return std::exp(v);
    }
  }

dd_f_log_e_quadratic()

static double EshelbianCore::dd_f_log_e_quadratic ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 74 of file EshelbianCore.hpp.

                                                     {
    if (v > v_max) {
      return static_cast<double>(std::exp(v_max));
    } else {
      return static_cast<double>(std::exp(v));
    } 
  }

f_linear()

static double EshelbianCore::f_linear ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 158 of file EshelbianCore.hpp.

{ return v + 1; }

f_log()

static double EshelbianCore::f_log ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 137 of file EshelbianCore.hpp.

                                      {
    return pow(EshelbianCore::exponentBase, v);
  }

f_log_e()

static double EshelbianCore::f_log_e ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 82 of file EshelbianCore.hpp.

                                        {
    if constexpr(use_quadratic_exp) {
      return f_log_e_quadratic(v);
    } else {
      if (v > v_max)
        // y = exp(v_max) * v + exp(v_max) * (1 - v_max);
        // y/exp(v_max) = v + (1 - v_max);
        // y/exp(v_max) - (1 - v_max) = v;
        return std::exp(v_max) * v + std::exp(v_max) * (1 - v_max);
      else 
        return std::exp(v);
    }
  }

f_log_e_quadratic()

static double EshelbianCore::f_log_e_quadratic ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 54 of file EshelbianCore.hpp.

                                                  {
    if (v > v_max) {
      double e = static_cast<double>(std::exp(v_max));
      double dv = v - v_max;
      return 0.5 * e * dv * dv + e * dv + e;
    } else {
      return static_cast<double>(std::exp(v));
    } 
  }

getBc()

template<typename BC >

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

inline

Examples

EshelbianPlasticity.cpp.

Definition at line 243 of file EshelbianCore.hpp.

                                                                            {
    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());
      }
      Range faces;
      CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                                 true);
      bc_vec_ptr->emplace_back(it->getName(), block_attributes, faces);
    }
    MoFEMFunctionReturn(0);
  }

getExternalStrain()

MoFEMErrorCode EshelbianCore::getExternalStrain

(

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 7013 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

(

)

Examples

EshelbianPlasticity.cpp.

Definition at line 922 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_skeleton"};
  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 EshelbianCore::getSpatialDispBc

(

)

[Getting norms]

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 6835 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 EshelbianCore::getSpatialRotationBc ( )

inline

Examples

ep.cpp.

Definition at line 271 of file EshelbianCore.hpp.

                                               {
    MoFEMFunctionBegin;
    bcSpatialRotationVecPtr = boost::make_shared<BcRotVec>();
    CHKERR getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_BC", 4);
    CHKERR getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_AXIS_BC", 7);
 
    auto ts_rotation =
        boost::make_shared<DynamicRelaxationTimeScale>("rotation_history.txt");
    for (auto &bc : *bcSpatialRotationVecPtr) {
      timeScaleMap[bc.blockName] =
          GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
              ts_rotation, "rotation_history", ".txt", bc.blockName);
    }
 
    MoFEMFunctionReturn(0);
  }

getSpatialTractionBc()

MoFEMErrorCode EshelbianCore::getSpatialTractionBc

(

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 6929 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 EshelbianCore::getSpatialTractionFreeBc ( const EntityHandle  meshset = 0 )

inline

Examples

ep.cpp.

Definition at line 305 of file EshelbianCore.hpp.

                                                           {
    bcSpatialFreeTractionVecPtr =
        boost::shared_ptr<TractionFreeBc>(new TractionFreeBc());
    return getTractionFreeBc(meshset, bcSpatialFreeTractionVecPtr, "CONTACT");
  }

gettingNorms()

MoFEMErrorCode EshelbianCore::gettingNorms

(

)

[Getting norms]

Examples

EshelbianPlasticity.cpp.

Definition at line 6767 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

Examples

EshelbianPlasticity.cpp.

Definition at line 2073 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 EshelbianCore::inv_d_f_linear ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 163 of file EshelbianCore.hpp.

{ return 0; }

inv_d_f_log()

static double EshelbianCore::inv_d_f_log ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 151 of file EshelbianCore.hpp.

                                            {
    return (1. / v) / log(EshelbianCore::exponentBase);
  }

inv_d_f_log_e()

static double EshelbianCore::inv_d_f_log_e ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 124 of file EshelbianCore.hpp.

                                              {
    if (v > std::exp(v_min))
      return 1. / v;
    else
      return 1. / exp(v_min);
  }

inv_dd_f_linear()

static double EshelbianCore::inv_dd_f_linear ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 164 of file EshelbianCore.hpp.

{ return 0; }

inv_dd_f_log()

static double EshelbianCore::inv_dd_f_log ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 154 of file EshelbianCore.hpp.

                                             {
    return -(1. / (v * v)) / log(EshelbianCore::exponentBase);
  }

inv_dd_f_log_e()

static double EshelbianCore::inv_dd_f_log_e ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 130 of file EshelbianCore.hpp.

                                               {
    if (v > std::exp(v_min))
      return -1. / (v * v);
    else
      return 0.;
  }

inv_f_linear()

static double EshelbianCore::inv_f_linear ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 162 of file EshelbianCore.hpp.

{ return v - 1; }

inv_f_log()

static double EshelbianCore::inv_f_log ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 148 of file EshelbianCore.hpp.

                                          {
    return log(v) / log(EshelbianCore::exponentBase);
  }

inv_f_log_e()

static double EshelbianCore::inv_f_log_e ( const double  v )

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 115 of file EshelbianCore.hpp.

                                            {
    if (v > std::exp(v_min))
      return std::log(v);
    else
      // y = exp(v_min) * v + exp(v_min) * (1 - v_min);
      // y/exp(v_min) = v + (1 - v_min);
      // y/exp(v_min) - (1 - v_min) = v;
      return v / exp(v_min) - (1 - v_min);
  }

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 = {}
	)

Examples

EshelbianPlasticity.cpp.

Definition at line 3713 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 = {}
	)

Examples

EshelbianPlasticity.cpp.

Definition at line 4207 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 - 1, 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));
  post_proc_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);
  post_proc_ptr->getOpPtrVector().push_back(op_loop_domain_side);
 
  // evaluated in side domain, that is op_loop_domain_side
  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<SPACE_DIM, SPACE_DIM>(
          piolaStress, dataAtPts->getApproxPAtPts()));
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateHTensorTensorField<SPACE_DIM, SPACE_DIM>(
          bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  op_loop_domain_side->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>(
          rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  if (noStretch) {
    op_loop_domain_side->getOpPtrVector().push_back(
        physicalEquations->returnOpCalculateStretchFromStress(
            dataAtPts, physicalEquations));
  } else {
    op_loop_domain_side->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
            stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
  }
 
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
  OpPPMap::DataMapMat vec_fields;
  vec_fields["HybridDisplacement"] = hybrid_disp;
  vec_fields["Omega"] = dataAtPts->getRotAxisAtPts();
  OpPPMap::DataMapMat mat_fields;
  mat_fields["PiolaStress"] = dataAtPts->getApproxPAtPts();
  mat_fields["HybridDisplacementGradient"] =
      dataAtPts->getGradHybridDispAtPts();
  OpPPMap::DataMapMat mat_fields_symm;
  mat_fields_symm["LogSpatialStretch"] = dataAtPts->getLogStretchTensorAtPts();
 
  post_proc_ptr->getOpPtrVector().push_back(
 
      new OpPPMap(
 
          post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
 
          {},
 
          vec_fields,
 
          mat_fields,
 
          mat_fields_symm
 
          )
 
  );
 
  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

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 1406 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_with_both_nodes_not_at_front =
          subtract(crack_front_edges, 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 898 of file EshelbianPlasticity.cpp.

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

saveOrgCoords()

MoFEMErrorCode EshelbianCore::saveOrgCoords

(

)

Examples

EshelbianPlasticity.cpp.

Definition at line 6735 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 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
	)

Examples

EshelbianPlasticity.cpp.

Definition at line 2307 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

(

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 3652 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

Examples

EshelbianPlasticity.cpp.

Definition at line 3079 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 EshelbianCore::setElasticElementOps

(

const int

tag

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 3178 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

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 3201 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 EshelbianCore::setFaceElementOps	(	const bool	add_elastic,
		const bool	add_material,
		boost::shared_ptr< FaceElementForcesAndSourcesCore > &	fe_rhs,
		boost::shared_ptr< FaceElementForcesAndSourcesCore > &	fe_lhs
	)

Examples

EshelbianPlasticity.cpp.

Definition at line 2920 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 EshelbianCore::setNewFrontCoordinates

(

)

Examples

EshelbianPlasticity.cpp.

Definition at line 5959 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

Examples

EshelbianPlasticity.cpp.

Definition at line 2428 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
	)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 3561 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
	)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 3431 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 362 of file EshelbianCore.hpp.

Member Data Documentation

a00FieldList

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

Definition at line 446 of file EshelbianCore.hpp.

a00RangeList

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

Definition at line 448 of file EshelbianCore.hpp.

addCrackMeshsetId

int EshelbianCore::addCrackMeshsetId = 1000

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 31 of file EshelbianCore.hpp.

alphaOmega

double EshelbianCore::alphaOmega = 0

Examples

EshelbianPlasticity.cpp.

Definition at line 221 of file EshelbianCore.hpp.

alphaRho

double EshelbianCore::alphaRho = 0

Examples

EshelbianPlasticity.cpp.

Definition at line 222 of file EshelbianCore.hpp.

alphaU

double EshelbianCore::alphaU = 0

Examples

EshelbianPlasticity.cpp.

Definition at line 219 of file EshelbianCore.hpp.

alphaW

double EshelbianCore::alphaW = 0

Examples

EshelbianPlasticity.cpp.

Definition at line 220 of file EshelbianCore.hpp.

AnalyticalExprPythonPtr

boost::shared_ptr<AnalyticalExprPython> EshelbianCore::AnalyticalExprPythonPtr

Examples

EshelbianPlasticity.cpp.

Definition at line 179 of file EshelbianCore.hpp.

aoS

AO EshelbianCore::aoS = PETSC_NULLPTR

Definition at line 444 of file EshelbianCore.hpp.

bcSpatialAnalyticalDisplacementVecPtr

boost::shared_ptr<AnalyticalDisplacementBcVec> EshelbianCore::bcSpatialAnalyticalDisplacementVecPtr

Examples

EshelbianPlasticity.cpp.

Definition at line 235 of file EshelbianCore.hpp.

bcSpatialAnalyticalTractionVecPtr

boost::shared_ptr<AnalyticalTractionBcVec> EshelbianCore::bcSpatialAnalyticalTractionVecPtr

Examples

EshelbianPlasticity.cpp.

Definition at line 236 of file EshelbianCore.hpp.

bcSpatialDispVecPtr

boost::shared_ptr<BcDispVec> EshelbianCore::bcSpatialDispVecPtr

Examples

EshelbianPlasticity.cpp.

Definition at line 229 of file EshelbianCore.hpp.

bcSpatialFreeTractionVecPtr

boost::shared_ptr<TractionFreeBc> EshelbianCore::bcSpatialFreeTractionVecPtr

Examples

EshelbianPlasticity.cpp.

Definition at line 232 of file EshelbianCore.hpp.

bcSpatialNormalDisplacementVecPtr

boost::shared_ptr<NormalDisplacementBcVec> EshelbianCore::bcSpatialNormalDisplacementVecPtr

Examples

EshelbianPlasticity.cpp.

Definition at line 233 of file EshelbianCore.hpp.

bcSpatialPressureVecPtr

boost::shared_ptr<PressureBcVec> EshelbianCore::bcSpatialPressureVecPtr

Examples

EshelbianPlasticity.cpp.

Definition at line 237 of file EshelbianCore.hpp.

bcSpatialRotationVecPtr

boost::shared_ptr<BcRotVec> EshelbianCore::bcSpatialRotationVecPtr

Examples

EshelbianPlasticity.cpp.

Definition at line 230 of file EshelbianCore.hpp.

bcSpatialTractionVecPtr

boost::shared_ptr<TractionBcVec> EshelbianCore::bcSpatialTractionVecPtr

Examples

EshelbianPlasticity.cpp.

Definition at line 231 of file EshelbianCore.hpp.

bitAdjEnt

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

bit ref level for parent

Examples

EshelbianPlasticity.cpp.

Definition at line 429 of file EshelbianCore.hpp.

bitAdjEntMask

BitRefLevel EshelbianCore::bitAdjEntMask

Initial value:

=
      BitRefLevel().set()

bit ref level for parent parent

Examples

EshelbianPlasticity.cpp.

Definition at line 430 of file EshelbianCore.hpp.

bitAdjParent

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

bit ref level for parent

Examples

EshelbianPlasticity.cpp.

Definition at line 426 of file EshelbianCore.hpp.

bitAdjParentMask

BitRefLevel EshelbianCore::bitAdjParentMask

Initial value:

=
      BitRefLevel().set()

bit ref level for parent parent

Examples

EshelbianPlasticity.cpp.

Definition at line 427 of file EshelbianCore.hpp.

bubbleField

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

Examples

EshelbianPlasticity.cpp.

Definition at line 203 of file EshelbianCore.hpp.

contactDisp

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

Examples

EshelbianPlasticity.cpp.

Definition at line 200 of file EshelbianCore.hpp.

contactElement

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

Examples

EshelbianPlasticity.cpp.

Definition at line 209 of file EshelbianCore.hpp.

contactFaces

boost::shared_ptr<Range> EshelbianCore::contactFaces

Examples

EshelbianPlasticity.cpp.

Definition at line 414 of file EshelbianCore.hpp.

contactRefinementLevels

int EshelbianCore::contactRefinementLevels = 1

Examples

EshelbianPlasticity.cpp.

Definition at line 224 of file EshelbianCore.hpp.

contactTreeRhs

boost::shared_ptr<ContactTree> EshelbianCore::contactTreeRhs

Make a contact tree.

Examples

EshelbianPlasticity.cpp.

Definition at line 185 of file EshelbianCore.hpp.

crackFaces

boost::shared_ptr<Range> EshelbianCore::crackFaces

Examples

EshelbianPlasticity.cpp.

Definition at line 415 of file EshelbianCore.hpp.

crackHybridIs

SmartPetscObj<IS> EshelbianCore::crackHybridIs

Definition at line 445 of file EshelbianCore.hpp.

crackingOn

PetscBool EshelbianCore::crackingOn = PETSC_FALSE

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 22 of file EshelbianCore.hpp.

crackingStartTime

double EshelbianCore::crackingStartTime = 0

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 23 of file EshelbianCore.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 44 of file EshelbianCore.hpp.

dataAtPts

boost::shared_ptr<DataAtIntegrationPts> EshelbianCore::dataAtPts

Examples

EshelbianPlasticity.cpp.

Definition at line 177 of file EshelbianCore.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 45 of file EshelbianCore.hpp.

dM

SmartPetscObj<DM> EshelbianCore::dM

Coupled problem all fields.

Examples

EshelbianPlasticity.cpp.

Definition at line 187 of file EshelbianCore.hpp.

dmElastic

SmartPetscObj<DM> EshelbianCore::dmElastic

Elastic problem.

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 188 of file EshelbianCore.hpp.

dmMaterial

SmartPetscObj<DM> EshelbianCore::dmMaterial

Material problem.

Definition at line 189 of file EshelbianCore.hpp.

dmPrjSpatial

SmartPetscObj<DM> EshelbianCore::dmPrjSpatial

Projection spatial displacement.

Examples

EshelbianPlasticity.cpp.

Definition at line 190 of file EshelbianCore.hpp.

dynamicRelaxation

PetscBool EshelbianCore::dynamicRelaxation

inlinestatic

Initial value:

=
      PETSC_FALSE

Examples

EshelbianOperators.cpp, EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 20 of file EshelbianCore.hpp.

dynamicStep

int EshelbianCore::dynamicStep

inlinestatic

Initial value:

=
      0

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 28 of file EshelbianCore.hpp.

dynamicTime

double EshelbianCore::dynamicTime

inlinestatic

Initial value:

=
      0

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 26 of file EshelbianCore.hpp.

edgeExchange

CommInterface::EntitiesPetscVector EshelbianCore::edgeExchange

Examples

EshelbianPlasticity.cpp.

Definition at line 437 of file EshelbianCore.hpp.

elasticBcLhs

boost::shared_ptr<FaceElementForcesAndSourcesCore> EshelbianCore::elasticBcLhs

Examples

EshelbianPlasticity.cpp.

Definition at line 183 of file EshelbianCore.hpp.

elasticBcRhs

boost::shared_ptr<FaceElementForcesAndSourcesCore> EshelbianCore::elasticBcRhs

Examples

EshelbianPlasticity.cpp.

Definition at line 184 of file EshelbianCore.hpp.

elasticFeLhs

boost::shared_ptr<VolumeElementForcesAndSourcesCore> EshelbianCore::elasticFeLhs

Examples

EshelbianPlasticity.cpp.

Definition at line 182 of file EshelbianCore.hpp.

elasticFeRhs

boost::shared_ptr<VolumeElementForcesAndSourcesCore> EshelbianCore::elasticFeRhs

Examples

EshelbianPlasticity.cpp.

Definition at line 181 of file EshelbianCore.hpp.

elementVolumeName

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

Examples

EshelbianPlasticity.cpp.

Definition at line 205 of file EshelbianCore.hpp.

energyReleaseSelector

enum EnergyReleaseSelector EshelbianCore::energyReleaseSelector

inlinestatic

Initial value:

=
      GRIFFITH_SKELETON

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 33 of file EshelbianCore.hpp.

exponentBase

double EshelbianCore::exponentBase = exp(1)

static

Examples

EshelbianPlasticity.cpp.

Definition at line 42 of file EshelbianCore.hpp.

externalStrainVecPtr

boost::shared_ptr<ExternalStrainVec> EshelbianCore::externalStrainVecPtr

Examples

EshelbianPlasticity.cpp.

Definition at line 238 of file EshelbianCore.hpp.

f

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

static

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 43 of file EshelbianCore.hpp.

faceExchange

CommInterface::EntitiesPetscVector EshelbianCore::faceExchange

Examples

EshelbianPlasticity.cpp.

Definition at line 436 of file EshelbianCore.hpp.

frontAdjEdges

boost::shared_ptr<Range> EshelbianCore::frontAdjEdges

Examples

EshelbianPlasticity.cpp.

Definition at line 417 of file EshelbianCore.hpp.

frontEdges

boost::shared_ptr<Range> EshelbianCore::frontEdges

Examples

EshelbianPlasticity.cpp.

Definition at line 416 of file EshelbianCore.hpp.

frontLayers

int EshelbianCore::frontLayers = 3

Definition at line 225 of file EshelbianCore.hpp.

frontVertices

boost::shared_ptr<Range> EshelbianCore::frontVertices

Examples

EshelbianPlasticity.cpp.

Definition at line 418 of file EshelbianCore.hpp.

gradApproximator

enum RotSelector EshelbianCore::gradApproximator = LARGE_ROT

inlinestatic

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 16 of file EshelbianCore.hpp.

griffithEnergy

double EshelbianCore::griffithEnergy = 1

inlinestatic

Griffith energy.

Examples

EshelbianPlasticity.cpp.

Definition at line 32 of file EshelbianCore.hpp.

hybridSpatialDisp

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

Examples

EshelbianPlasticity.cpp.

Definition at line 198 of file EshelbianCore.hpp.

internalStressInterpOrder

int EshelbianCore::internalStressInterpOrder

inlinestatic

Initial value:

=
      1

Examples

EshelbianPlasticity.cpp.

Definition at line 37 of file EshelbianCore.hpp.

internalStressTagName

std::string EshelbianCore::internalStressTagName

inlinestatic

Initial value:

=
      ""

Examples

EshelbianPlasticity.cpp.

Definition at line 35 of file EshelbianCore.hpp.

internalStressVoigt

PetscBool EshelbianCore::internalStressVoigt

inlinestatic

Initial value:

=
      PETSC_FALSE

Examples

EshelbianPlasticity.cpp.

Definition at line 39 of file EshelbianCore.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 47 of file EshelbianCore.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 48 of file EshelbianCore.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 46 of file EshelbianCore.hpp.

l2UserBaseScale

PetscBool EshelbianCore::l2UserBaseScale = PETSC_TRUE

inlinestatic

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 30 of file EshelbianCore.hpp.

listTagsToTransfer

std::vector<Tag> EshelbianCore::listTagsToTransfer

list of tags to transfer to postprocessor

Examples

EshelbianPlasticity.cpp.

Definition at line 441 of file EshelbianCore.hpp.

materialH1Positions

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

Examples

EshelbianPlasticity.cpp.

Definition at line 197 of file EshelbianCore.hpp.

materialOrder

int EshelbianCore::materialOrder = 1

Examples

EshelbianPlasticity.cpp.

Definition at line 218 of file EshelbianCore.hpp.

materialSkeletonFaces

boost::shared_ptr<Range> EshelbianCore::materialSkeletonFaces

Definition at line 420 of file EshelbianCore.hpp.

maxMovedFaces

boost::shared_ptr<Range> EshelbianCore::maxMovedFaces

Examples

EshelbianPlasticity.cpp.

Definition at line 421 of file EshelbianCore.hpp.

mField

MoFEM::Interface& EshelbianCore::mField

Examples

EshelbianPlasticity.cpp.

Definition at line 175 of file EshelbianCore.hpp.

naturalBcElement

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

Examples

EshelbianPlasticity.cpp.

Definition at line 206 of file EshelbianCore.hpp.

nbCrackFaces

int EshelbianCore::nbCrackFaces = 0

Examples

EshelbianPlasticity.cpp.

Definition at line 450 of file EshelbianCore.hpp.

nbJIntegralLevels

int EshelbianCore::nbJIntegralLevels

inlinestatic

Initial value:

=
      0

Examples

EshelbianPlasticity.cpp.

Definition at line 24 of file EshelbianCore.hpp.

noStretch

PetscBool EshelbianCore::noStretch = PETSC_FALSE

inlinestatic

Examples

EshelbianPlasticity.cpp.

Definition at line 18 of file EshelbianCore.hpp.

parentAdjSkeletonFunctionDim2

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

Examples

EshelbianPlasticity.cpp.

Definition at line 424 of file EshelbianCore.hpp.

physicalEquations

boost::shared_ptr<PhysicalEquations> EshelbianCore::physicalEquations

Examples

EshelbianPlasticity.cpp.

Definition at line 178 of file EshelbianCore.hpp.

piolaStress

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

Examples

EshelbianPlasticity.cpp.

Definition at line 192 of file EshelbianCore.hpp.

rotAxis

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

Examples

EshelbianPlasticity.cpp.

Definition at line 202 of file EshelbianCore.hpp.

rotSelector

enum RotSelector EshelbianCore::rotSelector = LARGE_ROT

inlinestatic

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 15 of file EshelbianCore.hpp.

S

Mat EshelbianCore::S = PETSC_NULLPTR

Definition at line 443 of file EshelbianCore.hpp.

setSingularity

PetscBool EshelbianCore::setSingularity = PETSC_TRUE

inlinestatic

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 19 of file EshelbianCore.hpp.

skeletonElement

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

Examples

EshelbianPlasticity.cpp.

Definition at line 208 of file EshelbianCore.hpp.

skeletonFaces

boost::shared_ptr<Range> EshelbianCore::skeletonFaces

Examples

EshelbianPlasticity.cpp.

Definition at line 419 of file EshelbianCore.hpp.

skinElement

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

Examples

EshelbianPlasticity.cpp.

Definition at line 207 of file EshelbianCore.hpp.

solTSStep

SmartPetscObj<Vec> EshelbianCore::solTSStep

Examples

EshelbianPlasticity.cpp.

Definition at line 433 of file EshelbianCore.hpp.

spaceH1Order

int EshelbianCore::spaceH1Order = -1

Examples

EshelbianPlasticity.cpp.

Definition at line 217 of file EshelbianCore.hpp.

spaceOrder

int EshelbianCore::spaceOrder = 2

Examples

EshelbianPlasticity.cpp.

Definition at line 216 of file EshelbianCore.hpp.

spatialH1Disp

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

Examples

EshelbianPlasticity.cpp.

Definition at line 196 of file EshelbianCore.hpp.

spatialL2Disp

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

Examples

EshelbianPlasticity.cpp.

Definition at line 194 of file EshelbianCore.hpp.

stretchSelector

enum StretchSelector EshelbianCore::stretchSelector = LOG

inlinestatic

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 17 of file EshelbianCore.hpp.

stretchTensor

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

Examples

EshelbianPlasticity.cpp.

Definition at line 201 of file EshelbianCore.hpp.

symmetrySelector

enum SymmetrySelector EshelbianCore::symmetrySelector = NOT_SYMMETRIC

inlinestatic

Examples

EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 14 of file EshelbianCore.hpp.

timeScaleMap

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

Examples

EshelbianPlasticity.cpp.

Definition at line 240 of file EshelbianCore.hpp.

use_quadratic_exp

constexpr bool EshelbianCore::use_quadratic_exp = true

inlinestaticconstexpr

Definition at line 50 of file EshelbianCore.hpp.

v_max

constexpr double EshelbianCore::v_max = 12

inlinestaticconstexpr

Definition at line 51 of file EshelbianCore.hpp.

v_min

constexpr double EshelbianCore::v_min = -v_max

inlinestaticconstexpr

Definition at line 52 of file EshelbianCore.hpp.

vertexExchange

CommInterface::EntitiesPetscVector EshelbianCore::vertexExchange

Examples

EshelbianPlasticity.cpp.

Definition at line 438 of file EshelbianCore.hpp.

volumeExchange

CommInterface::EntitiesPetscVector EshelbianCore::volumeExchange

Examples

EshelbianPlasticity.cpp.

Definition at line 435 of file EshelbianCore.hpp.

---

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

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