v0.15.0
Loading...
Searching...
No Matches
Classes | Public Member Functions | Static Public Member Functions | Public Attributes | Static Public Attributes | Friends | List of all members
EshelbianCore Struct Reference

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

Inheritance diagram for EshelbianCore:
[legend]
Collaboration diagram for EshelbianCore:
[legend]

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::InterfacemField
 
boost::shared_ptr< DataAtIntegrationPtsdataAtPts
 
boost::shared_ptr< PhysicalEquations > physicalEquations
 
boost::shared_ptr< AnalyticalExprPython > AnalyticalExprPythonPtr
 
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
 
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
 
boost::shared_ptr< FaceElementForcesAndSourcesCoreelasticBcLhs
 
boost::shared_ptr< FaceElementForcesAndSourcesCoreelasticBcRhs
 
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< RangecontactFaces
 
boost::shared_ptr< RangecrackFaces
 
boost::shared_ptr< RangefrontEdges
 
boost::shared_ptr< RangefrontAdjEdges
 
boost::shared_ptr< RangefrontVertices
 
boost::shared_ptr< RangeskeletonFaces
 
boost::shared_ptr< RangematerialSkeletonFaces
 
boost::shared_ptr< RangemaxMovedFaces
 
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< TaglistTagsToTransfer
 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
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 12 of file EshelbianCore.hpp.

Constructor & Destructor Documentation

◆ EshelbianCore()

EshelbianCore::EshelbianCore ( MoFEM::Interface m_field)
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 916 of file EshelbianPlasticity.cpp.

916 : mField(m_field) {
917 CHK_THROW_MESSAGE(getOptions(), "getOptions failed");
918}
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
MoFEM::Interface & mField
MoFEMErrorCode getOptions()

◆ ~EshelbianCore()

EshelbianCore::~EshelbianCore ( )
virtualdefault

Member Function Documentation

◆ addBoundaryFiniteElement()

MoFEMErrorCode EshelbianCore::addBoundaryFiniteElement ( const EntityHandle  meshset = 0)
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 1583 of file EshelbianPlasticity.cpp.

1583 {
1585
1586 Range meshset_ents;
1587 CHKERR mField.get_moab().get_entities_by_handle(meshset, meshset_ents);
1588
1589 auto set_fe_adjacency = [&](auto fe_name) {
1592 boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
1595 fe_name, MBTRI, *parentAdjSkeletonFunctionDim2);
1597 };
1598
1599 // set finite element fields
1600 auto add_field_to_fe = [this](const std::string fe,
1601 const std::string field_name) {
1608 };
1609
1611
1612 Range natural_bc_elements;
1613 if (bcSpatialDispVecPtr) {
1614 for (auto &v : *bcSpatialDispVecPtr) {
1615 natural_bc_elements.merge(v.faces);
1616 }
1617 }
1619 for (auto &v : *bcSpatialRotationVecPtr) {
1620 natural_bc_elements.merge(v.faces);
1621 }
1622 }
1624 for (auto &v : *bcSpatialNormalDisplacementVecPtr) {
1625 natural_bc_elements.merge(v.faces);
1626 }
1627 }
1630 natural_bc_elements.merge(v.faces);
1631 }
1632 }
1634 for (auto &v : *bcSpatialTractionVecPtr) {
1635 natural_bc_elements.merge(v.faces);
1636 }
1637 }
1639 for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
1640 natural_bc_elements.merge(v.faces);
1641 }
1642 }
1644 for (auto &v : *bcSpatialPressureVecPtr) {
1645 natural_bc_elements.merge(v.faces);
1646 }
1647 }
1648 natural_bc_elements = intersect(natural_bc_elements, meshset_ents);
1649
1651 CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
1653 CHKERR add_field_to_fe(naturalBcElement, piolaStress);
1654 CHKERR add_field_to_fe(naturalBcElement, hybridSpatialDisp);
1655 CHKERR set_fe_adjacency(naturalBcElement);
1657 }
1658
1659 auto get_skin = [&](auto &body_ents) {
1660 Skinner skin(&mField.get_moab());
1661 Range skin_ents;
1662 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
1663 return skin_ents;
1664 };
1665
1666 auto filter_true_skin = [&](auto &&skin) {
1667 Range boundary_ents;
1668 ParallelComm *pcomm =
1669 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1670 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
1671 PSTATUS_NOT, -1, &boundary_ents);
1672 return boundary_ents;
1673 };
1674
1676
1677 Range body_ents;
1678 CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM,
1679 body_ents);
1680 auto skin = filter_true_skin(get_skin(body_ents));
1681
1689 contactDisp);
1692 // CHKERR add_field_to_fe(skinElement, hybridSpatialDisp);
1693 // CHKERR add_field_to_fe(skinElement, hybridMaterialDisp);
1694
1696 }
1697
1699 if (contactFaces) {
1700 MOFEM_LOG("EP", Sev::inform)
1701 << "Contact elements " << contactFaces->size();
1705 CHKERR add_field_to_fe(contactElement, piolaStress);
1706 CHKERR add_field_to_fe(contactElement, hybridSpatialDisp);
1707 CHKERR add_field_to_fe(contactElement, contactDisp);
1708 CHKERR add_field_to_fe(contactElement, spatialH1Disp);
1709 CHKERR set_fe_adjacency(contactElement);
1711 }
1712 }
1713
1715 if (!skeletonFaces)
1716 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
1717 MOFEM_LOG("EP", Sev::inform)
1718 << "Skeleton elements " << skeletonFaces->size();
1722 CHKERR add_field_to_fe(skeletonElement, piolaStress);
1723 CHKERR add_field_to_fe(skeletonElement, hybridSpatialDisp);
1724 CHKERR add_field_to_fe(skeletonElement, contactDisp);
1725 CHKERR add_field_to_fe(skeletonElement, spatialH1Disp);
1726 CHKERR set_fe_adjacency(skeletonElement);
1728 }
1729
1730 // if (!mField.check_finite_element(materialSkeletonElement)) {
1731
1732 // Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM -
1733 // 2);
1734
1735 // Range front_elements;
1736 // for (auto l = 0; l != frontLayers; ++l) {
1737 // Range front_elements_layer;
1738 // CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM,
1739 // true,
1740 // front_elements_layer,
1741 // moab::Interface::UNION);
1742 // front_elements.merge(front_elements_layer);
1743 // front_edges.clear();
1744 // CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
1745 // SPACE_DIM - 2, true,
1746 // front_edges,
1747 // moab::Interface::UNION);
1748 // }
1749 // Range body_ents;
1750 // CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
1751 // body_ents); Range front_elements_faces; CHKERR
1752 // mField.get_moab().get_adjacencies(front_elements, SPACE_DIM - 1,
1753 // true, front_elements_faces,
1754 // moab::Interface::UNION);
1755 // auto body_skin = filter_true_skin(get_skin(body_ents));
1756 // auto front_elements_skin = filter_true_skin(get_skin(front_elements));
1757 // Range material_skeleton_faces =
1758 // subtract(front_elements_faces, front_elements_skin);
1759 // material_skeleton_faces.merge(intersect(front_elements_skin,
1760 // body_skin));
1761
1762 // #ifndef NDEBUG
1763 // CHKERR save_range(mField.get_moab(), "front_elements.vtk",
1764 // front_elements); CHKERR save_range(mField.get_moab(),
1765 // "front_elements_skin.vtk",
1766 // front_elements_skin);
1767 // CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
1768 // material_skeleton_faces);
1769 // #endif
1770
1771 // CHKERR mField.add_finite_element(materialSkeletonElement, MF_ZERO);
1772 // CHKERR mField.add_ents_to_finite_element_by_type(
1773 // material_skeleton_faces, MBTRI, materialSkeletonElement);
1774 // // CHKERR add_field_to_fe(materialSkeletonElement, eshelbyStress);
1775 // // CHKERR add_field_to_fe(materialSkeletonElement, hybridMaterialDisp);
1776 // CHKERR set_fe_adjacency(materialSkeletonElement);
1777 // CHKERR mField.build_finite_elements(materialSkeletonElement);
1778 // }
1779
1781}
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
static auto get_skin(MoFEM::Interface &m_field, Range body_ents)
constexpr int SPACE_DIM
@ MF_ZERO
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)=0
modify finite element table, only for advanced user
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
#define MOFEM_LOG(channel, severity)
Log.
const double v
phase velocity of light in medium (cm/ns)
constexpr auto field_name
const std::string skeletonElement
boost::shared_ptr< TractionBcVec > bcSpatialTractionVecPtr
boost::shared_ptr< Range > contactFaces
BitRefLevel bitAdjEnt
bit ref level for parent
const std::string materialH1Positions
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
const std::string spatialH1Disp
boost::shared_ptr< NormalDisplacementBcVec > bcSpatialNormalDisplacementVecPtr
const std::string piolaStress
boost::shared_ptr< AnalyticalDisplacementBcVec > bcSpatialAnalyticalDisplacementVecPtr
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > > parentAdjSkeletonFunctionDim2
boost::shared_ptr< Range > skeletonFaces
BitRefLevel bitAdjParentMask
bit ref level for parent parent
const std::string contactDisp
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
const std::string skinElement
boost::shared_ptr< AnalyticalTractionBcVec > bcSpatialAnalyticalTractionVecPtr
BitRefLevel bitAdjParent
bit ref level for parent
const std::string naturalBcElement
boost::shared_ptr< PressureBcVec > bcSpatialPressureVecPtr
const std::string hybridSpatialDisp
BitRefLevel bitAdjEntMask
bit ref level for parent parent
const std::string contactElement
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MPI_Comm & get_comm() const =0

◆ addCrackSurfaces()

MoFEMErrorCode EshelbianCore::addCrackSurfaces ( const bool  debug = false)
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 5881 of file EshelbianPlasticity.cpp.

5881 {
5883
5884 constexpr bool potential_crack_debug = false;
5885 if constexpr (potential_crack_debug) {
5886
5887 auto add_ents = get_range_from_block(mField, "POTENTIAL", SPACE_DIM - 1);
5888 Range crack_front_verts;
5889 CHKERR mField.get_moab().get_connectivity(*frontEdges, crack_front_verts,
5890 true);
5891 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
5892 crack_front_verts);
5893 Range crack_front_faces;
5894 CHKERR mField.get_moab().get_adjacencies(crack_front_verts, SPACE_DIM - 1,
5895 true, crack_front_faces,
5896 moab::Interface::UNION);
5897 crack_front_faces = intersect(crack_front_faces, add_ents);
5898 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
5899 crack_front_faces);
5900 CHKERR mField.getInterface<MeshsetsManager>()->addEntitiesToMeshset(
5901 BLOCKSET, addCrackMeshsetId, crack_front_faces);
5902 }
5903
5904 auto get_crack_faces = [&]() {
5905 if (maxMovedFaces) {
5906 return unite(*crackFaces, *maxMovedFaces);
5907 } else {
5908 return *crackFaces;
5909 }
5910 };
5911
5912 auto get_extended_crack_faces = [&]() {
5913 auto get_faces_of_crack_front_verts = [&](auto crack_faces_org) {
5914 ParallelComm *pcomm =
5915 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
5916
5917 Range crack_faces;
5918
5919 if (!pcomm->rank()) {
5920
5921 auto get_nodes = [&](auto &&e) {
5922 Range nodes;
5923 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
5924 "get connectivity");
5925 return nodes;
5926 };
5927
5928 auto get_adj = [&](auto &&e, auto dim,
5929 auto t = moab::Interface::UNION) {
5930 Range adj;
5932 mField.get_moab().get_adjacencies(e, dim, true, adj, t),
5933 "get adj");
5934 return adj;
5935 };
5936
5937 Range body_ents;
5938 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
5939 body_ents);
5940 auto body_skin = get_skin(mField, body_ents);
5941 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
5942 auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
5943 auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
5944 auto front_block_nodes = get_nodes(front_block_edges);
5945
5946 size_t s;
5947 do {
5948 s = crack_faces.size();
5949
5950 auto crack_face_nodes = get_nodes(crack_faces_org);
5951 auto crack_faces_edges =
5952 get_adj(crack_faces_org, 1, moab::Interface::UNION);
5953
5954 auto crack_skin = get_skin(mField, crack_faces_org);
5955 front_block_edges = subtract(front_block_edges, crack_skin);
5956 auto crack_skin_nodes = get_nodes(crack_skin);
5957 crack_skin_nodes.merge(front_block_nodes);
5958
5959 auto crack_skin_faces =
5960 get_adj(crack_skin, 2, moab::Interface::UNION);
5961 crack_skin_faces =
5962 subtract(subtract(crack_skin_faces, crack_faces_org), body_skin);
5963
5964 crack_faces = crack_faces_org;
5965 for (auto f : crack_skin_faces) {
5966 auto edges = intersect(
5967 get_adj(Range(f, f), 1, moab::Interface::UNION), crack_skin);
5968
5969 // if other edge is part of body skin, e.g. crack punching through
5970 // body surface
5971 if (edges.size() == 2) {
5972 edges.merge(
5973 intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
5974 body_skin_edges));
5975 }
5976
5977 if (edges.size() == 2) {
5978 auto edge_conn = get_nodes(Range(edges));
5979 auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
5980 crack_faces_org);
5981 if (faces.size() == 2) {
5982 auto edge0_conn = get_nodes(Range(edges[0], edges[0]));
5983 auto edge1_conn = get_nodes(Range(edges[1], edges[1]));
5984 auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
5985 crack_skin_nodes); // node at apex
5986 if (edges_conn.size() == 1) {
5987
5988 auto node_edges =
5989 subtract(intersect(get_adj(edges_conn, 1,
5990 moab::Interface::INTERSECT),
5991 crack_faces_edges),
5992 crack_skin); // nodes on crack surface, but not
5993 // at the skin
5994
5995 if (node_edges.size()) {
5998 CHKERR mField.get_moab().get_coords(edges_conn, &t_v0(0));
5999
6000 auto get_t_dir = [&](auto e_conn) {
6001 auto other_node = subtract(e_conn, edges_conn);
6003 CHKERR mField.get_moab().get_coords(other_node,
6004 &t_dir(0));
6005 t_dir(i) -= t_v0(i);
6006 return t_dir;
6007 };
6008
6010 t_ave_dir(i) =
6011 get_t_dir(edge0_conn)(i) + get_t_dir(edge1_conn)(i);
6012
6013 FTensor::Tensor1<double, SPACE_DIM> t_crack_surface_ave_dir;
6014 t_crack_surface_ave_dir(i) = 0;
6015 for (auto e : node_edges) {
6016 auto e_conn = get_nodes(Range(e, e));
6017 auto t_dir = get_t_dir(e_conn);
6018 t_crack_surface_ave_dir(i) += t_dir(i);
6019 }
6020
6021 auto dot = t_ave_dir(i) * t_crack_surface_ave_dir(i);
6022 // ave edges is in opposite direction to crack surface, so
6023 // thus crack is not turning back
6024 if (dot < -std::numeric_limits<double>::epsilon()) {
6025 crack_faces.insert(f);
6026 }
6027 } else {
6028 crack_faces.insert(f);
6029 }
6030 }
6031 }
6032 } else if (edges.size() == 3) {
6033 crack_faces.insert(f);
6034 }
6035
6036 // if other edge is part of geometry edge, e.g. keyway
6037 if (edges.size() == 1) {
6038 edges.merge(
6039 intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
6040 geometry_edges));
6041 edges.merge(
6042 intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
6043 front_block_edges));
6044 if (edges.size() == 2) {
6045 crack_faces.insert(f);
6046 continue;
6047 }
6048 }
6049 }
6050
6051 crack_faces_org = crack_faces;
6052
6053 } while (s != crack_faces.size());
6054 };
6055
6056 return crack_faces; // send_type(mField, crack_faces, MBTRI);
6057 };
6058
6059 return get_faces_of_crack_front_verts(get_crack_faces());
6060 };
6061
6062 if (debug) {
6063 CHKERR save_range(mField.get_moab(), "new_crack_surface_debug.vtk",
6064 get_extended_crack_faces());
6065 }
6066
6067 auto reconstruct_crack_faces = [&](auto crack_faces) {
6068 ParallelComm *pcomm =
6069 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
6070
6071 auto impl = [&]() {
6073
6074 Range new_crack_faces;
6075 if (!pcomm->rank()) {
6076
6077 auto get_nodes = [&](auto &&e) {
6078 Range nodes;
6079 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
6080 "get connectivity");
6081 return nodes;
6082 };
6083
6084 auto get_adj = [&](auto &&e, auto dim,
6085 auto t = moab::Interface::UNION) {
6086 Range adj;
6088 mField.get_moab().get_adjacencies(e, dim, true, adj, t),
6089 "get adj");
6090 return adj;
6091 };
6092
6093 auto get_test_on_crack_surface = [&]() {
6094 auto crack_faces_nodes =
6095 get_nodes(crack_faces); // nodes on crac faces
6096 auto crack_faces_tets =
6097 get_adj(crack_faces_nodes, 3,
6098 moab::Interface::UNION); // adjacent
6099 // tets to
6100 // crack
6101 // faces throug nodes
6102 auto crack_faces_tets_nodes =
6103 get_nodes(crack_faces_tets); // nodes on crack faces tets
6104 crack_faces_tets_nodes =
6105 subtract(crack_faces_tets_nodes, crack_faces_nodes);
6106 crack_faces_tets =
6107 subtract(crack_faces_tets, get_adj(crack_faces_tets_nodes, 3,
6108 moab::Interface::UNION));
6109 new_crack_faces =
6110 get_adj(crack_faces_tets, 2,
6111 moab::Interface::UNION); // adjacency faces to crack
6112 // faces through tets
6113 new_crack_faces.merge(crack_faces); // add original crack faces
6114
6115 return std::make_tuple(new_crack_faces, crack_faces_tets);
6116 };
6117
6118 auto carck_faces_test_edges = [&](auto faces, auto tets) {
6119 auto adj_tets_faces = get_adj(tets, 2, moab::Interface::UNION);
6120 auto adj_faces_edges = get_adj(subtract(faces, adj_tets_faces), 1,
6121 moab::Interface::UNION);
6122 auto adj_tets_edges = get_adj(tets, 1, moab::Interface::UNION);
6123 auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
6124 auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
6125 adj_faces_edges.merge(geometry_edges); // geometry edges
6126 adj_faces_edges.merge(front_block_edges); // front block edges
6127
6128 auto boundary_tets_edges = intersect(adj_tets_edges, adj_faces_edges);
6129 auto boundary_test_nodes = get_nodes(boundary_tets_edges);
6130 auto boundary_test_nodes_edges =
6131 get_adj(boundary_test_nodes, 1, moab::Interface::UNION);
6132 auto boundary_test_nodes_edges_nodes = subtract(
6133 get_nodes(boundary_test_nodes_edges), boundary_test_nodes);
6134
6135 boundary_tets_edges =
6136 subtract(boundary_test_nodes_edges,
6137 get_adj(boundary_test_nodes_edges_nodes, 1,
6138 moab::Interface::UNION));
6139
6140 Range body_ents;
6141 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
6142 body_ents);
6143 auto body_skin = get_skin(mField, body_ents);
6144
6145 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
6146 body_skin_edges = intersect(get_adj(tets, 1, moab::Interface::UNION),
6147 body_skin_edges);
6148 body_skin = intersect(body_skin, adj_tets_faces);
6149 body_skin_edges = subtract(
6150 body_skin_edges, get_adj(body_skin, 1, moab::Interface::UNION));
6151
6152 save_range(mField.get_moab(), "body_skin_edges.vtk", body_skin_edges);
6153 for (auto e : body_skin_edges) {
6154 auto adj_tet = intersect(
6155 get_adj(Range(e, e), 3, moab::Interface::INTERSECT), tets);
6156 if (adj_tet.size() == 1) {
6157 boundary_tets_edges.insert(e);
6158 }
6159 }
6160
6161 return boundary_tets_edges;
6162 };
6163
6164 auto p = get_test_on_crack_surface();
6165 auto &[new_crack_faces, crack_faces_tets] = p;
6166
6167 if (debug) {
6168 CHKERR save_range(mField.get_moab(), "hole_crack_faces_debug.vtk",
6169 crack_faces);
6170 CHKERR save_range(mField.get_moab(), "new_crack_faces_debug.vtk",
6171 new_crack_faces);
6172 CHKERR save_range(mField.get_moab(), "new_crack_tets_debug.vtk",
6173 crack_faces_tets);
6174 }
6175
6176 auto boundary_tets_edges =
6177 carck_faces_test_edges(new_crack_faces, crack_faces_tets);
6178 CHKERR save_range(mField.get_moab(), "boundary_tets_edges.vtk",
6179 boundary_tets_edges);
6180
6181 auto resolve_surface = [&](auto boundary_tets_edges,
6182 auto crack_faces_tets) {
6183 auto boundary_tets_edges_nodes = get_nodes(boundary_tets_edges);
6184 auto crack_faces_tets_faces =
6185 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
6186
6187 Range all_removed_faces;
6188 Range all_removed_tets;
6189 int counter = 0;
6190
6191 int size = 0;
6192 while (size != crack_faces_tets.size()) {
6193 auto tets_faces =
6194 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
6195 auto skin_tets = get_skin(mField, crack_faces_tets);
6196 auto skin_skin =
6197 get_skin(mField, subtract(crack_faces_tets_faces, tets_faces));
6198 auto skin_skin_nodes = get_nodes(skin_skin);
6199
6200 size = crack_faces_tets.size();
6201 MOFEM_LOG("SELF", Sev::inform)
6202 << "Crack faces tets size " << crack_faces_tets.size()
6203 << " crack faces size " << crack_faces_tets_faces.size();
6204 auto skin_tets_nodes = subtract(
6205 get_nodes(skin_tets),
6206 boundary_tets_edges_nodes); // not remove tets which are
6207 // adjagasent to crack faces nodes
6208 skin_tets_nodes = subtract(skin_tets_nodes, skin_skin_nodes);
6209
6210 Range removed_nodes;
6211 Range tets_to_remove;
6212 Range faces_to_remove;
6213 for (auto n : skin_tets_nodes) {
6214 auto tets =
6215 intersect(get_adj(Range(n, n), 3, moab::Interface::INTERSECT),
6216 crack_faces_tets);
6217 if (tets.size() == 0) {
6218 continue;
6219 }
6220
6221 auto hole_detetction = [&]() {
6222 auto adj_tets =
6223 get_adj(Range(n, n), 3, moab::Interface::INTERSECT);
6224 adj_tets =
6225 subtract(adj_tets,
6226 crack_faces_tets); // tetst adjacent to the node
6227 // but not part of crack surface
6228 if (adj_tets.size() == 0) {
6229 return std::make_pair(
6230 intersect(
6231 get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
6232 tets_faces),
6233 tets);
6234 }
6235
6236 std::vector<Range> tets_groups;
6237 auto test_adj_tets = adj_tets;
6238 while (test_adj_tets.size()) {
6239 auto seed_size = 0;
6240 Range seed = Range(test_adj_tets[0], test_adj_tets[0]);
6241 while (seed.size() != seed_size) {
6242 auto adj_faces =
6243 subtract(get_adj(seed, 2, moab::Interface::UNION),
6244 tets_faces); // edges which are not
6245 // part of the node
6246 seed_size = seed.size();
6247 seed.merge(
6248 intersect(get_adj(adj_faces, 3, moab::Interface::UNION),
6249 test_adj_tets ));
6250 }
6251 tets_groups.push_back(seed);
6252 test_adj_tets = subtract(test_adj_tets, seed);
6253 }
6254 if (tets_groups.size() == 1) {
6255
6256 return std::make_pair(
6257 intersect(
6258 get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
6259 tets_faces),
6260 tets);
6261
6262 }
6263
6264 Range tets_to_remove;
6265 Range faces_to_remove;
6266 for (auto &r : tets_groups) {
6267 auto f = get_adj(r, 2, moab::Interface::UNION);
6268 auto t = intersect(get_adj(f, 3, moab::Interface::UNION),
6269 crack_faces_tets); // tets
6270
6271 if (f.size() > faces_to_remove.size() ||
6272 faces_to_remove.size() == 0) {
6273 faces_to_remove = f;
6274 tets_to_remove = t; // largest group of tets
6275 }
6276 }
6277 MOFEM_LOG("EPSELF", Sev::inform)
6278 << "Hole detection: faces to remove "
6279 << faces_to_remove.size() << " tets to remove "
6280 << tets_to_remove.size();
6281 return std::make_pair(faces_to_remove, tets_to_remove);
6282 };
6283
6284 if (tets.size() < tets_to_remove.size() ||
6285 tets_to_remove.size() == 0) {
6286 removed_nodes = Range(n, n);
6287 auto [h_faces_to_remove, h_tets_to_remove] =
6288 hole_detetction(); // find faces and tets to remove
6289 faces_to_remove = h_faces_to_remove;
6290 tets_to_remove = h_tets_to_remove;
6291
6292 // intersect(
6293 // get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
6294 // tets_faces);
6295
6296 } // find tets which is largest adjacencty size, so that it is
6297 // removed first, and then faces are removed
6298 all_removed_faces.merge(faces_to_remove);
6299 all_removed_tets.merge(tets_to_remove);
6300 }
6301
6302 crack_faces_tets = subtract(crack_faces_tets, tets_to_remove);
6303 crack_faces_tets_faces =
6304 subtract(crack_faces_tets_faces, faces_to_remove);
6305
6306 if (debug) {
6308 "removed_nodes_" +
6309 boost::lexical_cast<std::string>(counter) + ".vtk",
6310 removed_nodes);
6312 "faces_to_remove_" +
6313 boost::lexical_cast<std::string>(counter) + ".vtk",
6314 faces_to_remove);
6316 "tets_to_remove_" +
6317 boost::lexical_cast<std::string>(counter) + ".vtk",
6318 tets_to_remove);
6320 "crack_faces_tets_faces_" +
6321 boost::lexical_cast<std::string>(counter) + ".vtk",
6322 crack_faces_tets_faces);
6324 "crack_faces_tets_" +
6325 boost::lexical_cast<std::string>(counter) + ".vtk",
6326 crack_faces_tets);
6327 }
6328 counter++;
6329 }
6330
6331 auto cese_internal_faces = [&]() {
6333 auto skin_tets = get_skin(mField, crack_faces_tets);
6334 auto adj_faces = get_adj(skin_tets, 2, moab::Interface::UNION);
6335 adj_faces =
6336 subtract(adj_faces, skin_tets); // remove skin tets faces
6337 auto adj_tets = get_adj(adj_faces, 3,
6338 moab::Interface::UNION); // tets which are
6339 // adjacent to skin
6340 crack_faces_tets =
6341 subtract(crack_faces_tets,
6342 adj_tets); // remove tets which are adjacent to
6343 // skin, so that they are not removed
6344 crack_faces_tets_faces =
6345 subtract(crack_faces_tets_faces, adj_faces);
6346
6347 all_removed_faces.merge(adj_faces);
6348 all_removed_tets.merge(adj_tets);
6349
6350
6351 MOFEM_LOG("EPSELF", Sev::inform)
6352 << "Remove internal faces size " << adj_faces.size()
6353 << " tets size " << adj_tets.size();
6355 };
6356
6357 auto case_only_one_free_edge = [&]() {
6359
6360 for (auto t : Range(crack_faces_tets)) {
6361
6362 auto adj_faces = get_adj(
6363 Range(t, t), 2,
6364 moab::Interface::UNION); // faces of tet which can be removed
6365 auto crack_surface_edges =
6366 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
6367 adj_faces),
6368 1,
6369 moab::Interface::UNION); // edges not on the tet but
6370 // on crack surface
6371 auto adj_edges =
6372 subtract(get_adj(Range(t, t), 1, moab::Interface::INTERSECT),
6373 crack_surface_edges); // free edges
6374 adj_edges = subtract(
6375 adj_edges,
6376 boundary_tets_edges); // edges which are not part of gemetry
6377
6378 if (adj_edges.size() == 1) {
6379 crack_faces_tets =
6380 subtract(crack_faces_tets,
6381 Range(t, t)); // remove tets which are adjacent to
6382 // skin, so that they are not removed
6383
6384 auto faces_to_remove =
6385 get_adj(adj_edges, 2, moab::Interface::UNION); // faces
6386 // which can
6387 // be removed
6388 crack_faces_tets_faces =
6389 subtract(crack_faces_tets_faces, faces_to_remove);
6390
6391 all_removed_faces.merge(faces_to_remove);
6392 all_removed_tets.merge(Range(t, t));
6393
6394 MOFEM_LOG("EPSELF", Sev::inform) << "Remove free one edges ";
6395 }
6396 }
6397
6398 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
6399 crack_faces_tets_faces =
6400 subtract(crack_faces_tets_faces, all_removed_faces);
6401
6403 };
6404
6405 auto cese_flat_tet = [&](auto max_adj_edges) {
6407
6408 Range body_ents;
6409 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
6410 body_ents);
6411 auto body_skin = get_skin(mField, body_ents);
6412 auto body_skin_edges =
6413 get_adj(body_skin, 1, moab::Interface::UNION);
6414
6415 for (auto t : Range(crack_faces_tets)) {
6416
6417 auto adj_faces = get_adj(
6418 Range(t, t), 2,
6419 moab::Interface::UNION); // faces of tet which can be removed
6420 auto crack_surface_edges =
6421 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
6422 adj_faces),
6423 1,
6424 moab::Interface::UNION); // edges not on the tet but
6425 // on crack surface
6426 auto adj_edges =
6427 subtract(get_adj(Range(t, t), 1, moab::Interface::INTERSECT),
6428 crack_surface_edges); // free edges
6429 adj_edges = subtract(adj_edges, body_skin_edges);
6430
6431 auto tet_edges = get_adj(Range(t, t), 1,
6432 moab::Interface::UNION); // edges of
6433 // tet
6434 tet_edges = subtract(tet_edges, adj_edges);
6435
6436 for (auto e : tet_edges) {
6437 constexpr int opposite_edge[] = {5, 3, 4, 1, 2, 0};
6438 auto get_side = [&](auto e) {
6439 int side, sense, offset;
6441 mField.get_moab().side_number(t, e, side, sense, offset),
6442 "get side number failed");
6443 return side;
6444 };
6445 auto get_side_ent = [&](auto side) {
6446 EntityHandle side_edge;
6448 mField.get_moab().side_element(t, 1, side, side_edge),
6449 "get side failed");
6450 return side_edge;
6451 };
6452 adj_edges.erase(get_side_ent(opposite_edge[get_side(e)]));
6453 }
6454
6455 if (adj_edges.size() <= max_adj_edges) {
6456
6457 double dot = 1;
6458 Range faces_to_remove;
6459 for (auto e : adj_edges) {
6460 auto edge_adj_faces =
6461 get_adj(Range(e, e), 2, moab::Interface::UNION);
6462 edge_adj_faces = intersect(edge_adj_faces, adj_faces);
6463 if (edge_adj_faces.size() != 2) {
6465 "Adj faces size is not 2 for edge " +
6466 boost::lexical_cast<std::string>(e));
6467 }
6468
6469 auto get_normal = [&](auto f) {
6472 mField.getInterface<Tools>()->getTriNormal(f, &t_n(0)),
6473 "get tri normal failed");
6474 return t_n;
6475 };
6476 auto t_n0 = get_normal(edge_adj_faces[0]);
6477 auto t_n1 = get_normal(edge_adj_faces[1]);
6478 auto get_sense = [&](auto f) {
6479 int side, sense, offset;
6480 CHK_MOAB_THROW(mField.get_moab().side_number(t, f, side,
6481 sense, offset),
6482 "get side number failed");
6483 return sense;
6484 };
6485 auto sense0 = get_sense(edge_adj_faces[0]);
6486 auto sense1 = get_sense(edge_adj_faces[1]);
6487 t_n0.normalize();
6488 t_n1.normalize();
6489
6491 auto dot_e = (sense0 * sense1) * t_n0(i) * t_n1(i);
6492 if (dot_e < dot || e == adj_edges[0]) {
6493 dot = dot_e;
6494 faces_to_remove = edge_adj_faces;
6495 }
6496 }
6497
6498 all_removed_faces.merge(faces_to_remove);
6499 all_removed_tets.merge(Range(t, t));
6500
6501 MOFEM_LOG("EPSELF", Sev::inform)
6502 << "Remove free edges on flat tet, with considered nb. of "
6503 "edges "
6504 << adj_edges.size();
6505 }
6506 }
6507
6508 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
6509 crack_faces_tets_faces =
6510 subtract(crack_faces_tets_faces, all_removed_faces);
6511
6513 };
6514
6515 CHK_THROW_MESSAGE(case_only_one_free_edge(),
6516 "Case only one free edge failed");
6517 for (auto max_adj_edges : {0, 1, 2, 3}) {
6518 CHK_THROW_MESSAGE(cese_flat_tet(max_adj_edges),
6519 "Case only one free edge failed");
6520 }
6521 CHK_THROW_MESSAGE(cese_internal_faces(),
6522 "Case internal faces failed");
6523
6524 if (debug) {
6526 "crack_faces_tets_faces_" +
6527 boost::lexical_cast<std::string>(counter) + ".vtk",
6528 crack_faces_tets_faces);
6530 "crack_faces_tets_" +
6531 boost::lexical_cast<std::string>(counter) + ".vtk",
6532 crack_faces_tets);
6533 }
6534
6535 return std::make_tuple(crack_faces_tets_faces, crack_faces_tets,
6536 all_removed_faces, all_removed_tets);
6537 };
6538
6539 auto [resolved_faces, resolved_tets, all_removed_faces,
6540 all_removed_tets] =
6541 resolve_surface(boundary_tets_edges, crack_faces_tets);
6542 resolved_faces.merge(subtract(crack_faces, all_removed_faces));
6543 if (debug) {
6544 CHKERR save_range(mField.get_moab(), "resolved_faces.vtk",
6545 resolved_faces);
6546 CHKERR save_range(mField.get_moab(), "resolved_tets.vtk",
6547 resolved_tets);
6548 }
6549
6550 crack_faces = resolved_faces;
6551 }
6552
6554 };
6555
6556 CHK_THROW_MESSAGE(impl(), "resolve new crack surfaces");
6557
6558 return crack_faces; // send_type(mField, crack_faces, MBTRI);
6559 };
6560
6561
6562 auto resolve_consisten_crack_extension = [&]() {
6564 auto crack_meshset =
6565 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
6567 auto meshset = crack_meshset->getMeshset();
6568
6569
6570 if (!mField.get_comm_rank()) {
6571 Range old_crack_faces;
6572 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTRI,
6573 old_crack_faces);
6574 auto extendeded_crack_faces = get_extended_crack_faces();
6575 auto reconstructed_crack_faces =
6576 subtract(reconstruct_crack_faces(extendeded_crack_faces),
6577 subtract(*crackFaces, old_crack_faces));
6578 if (nbCrackFaces >= reconstructed_crack_faces.size()) {
6579 MOFEM_LOG("EPSELF", Sev::warning)
6580 << "No new crack faces to add, skipping adding to meshset";
6581 extendeded_crack_faces = subtract(
6582 extendeded_crack_faces, subtract(*crackFaces, old_crack_faces));
6583 MOFEM_LOG("EPSELF", Sev::inform)
6584 << "Number crack faces size (extended) "
6585 << extendeded_crack_faces.size();
6586 CHKERR mField.get_moab().clear_meshset(&meshset, 1);
6587 CHKERR mField.get_moab().add_entities(meshset, extendeded_crack_faces);
6588 } else {
6589 CHKERR mField.get_moab().clear_meshset(&meshset, 1);
6590 CHKERR mField.get_moab().add_entities(meshset,
6591 reconstructed_crack_faces);
6592 MOFEM_LOG("EPSELF", Sev::inform)
6593 << "Number crack faces size (reconstructed) "
6594 << reconstructed_crack_faces.size();
6595 nbCrackFaces = reconstructed_crack_faces.size();
6596 }
6597 }
6598
6599 Range crack_faces;
6600 if (!mField.get_comm_rank()) {
6601 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTRI,
6602 crack_faces);
6603 }
6604 crack_faces = send_type(mField, crack_faces, MBTRI);
6605 if (mField.get_comm_rank()) {
6606 CHKERR mField.get_moab().clear_meshset(&meshset, 1);
6607 CHKERR mField.get_moab().add_entities(meshset, crack_faces);
6608 }
6609
6611 };
6612
6613 CHKERR resolve_consisten_crack_extension();
6614
6616};
static auto send_type(MoFEM::Interface &m_field, Range r, const EntityType type)
#define FTENSOR_INDEX(DIM, I)
Range get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
Definition adjoint.cpp:3290
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ BLOCKSET
static const bool debug
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
static constexpr int edges_conn[]
constexpr double t
plate stiffness
Definition plate.cpp:58
boost::shared_ptr< Range > maxMovedFaces
static boost::function< double(const double)> f
static int addCrackMeshsetId
boost::shared_ptr< Range > crackFaces
boost::shared_ptr< Range > frontEdges
Managing BitRefLevels.
virtual int get_comm_rank() const =0
Interface for managing meshsets containing materials and boundary conditions.
Auxiliary tools.
Definition Tools.hpp:19
static MoFEMErrorCode getTriNormal(const double *coords, double *normal, double *d_normal=nullptr)
Get the Tri Normal objectGet triangle normal.
Definition Tools.cpp:353
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
auto save_range

◆ addDMs()

MoFEMErrorCode EshelbianCore::addDMs ( const BitRefLevel  bit = BitRefLevel().set(0),
const EntityHandle  meshset = 0 
)
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 1783 of file EshelbianPlasticity.cpp.

1784 {
1786
1787 // find adjacencies between finite elements and dofs
1789
1790 // Create coupled problem
1791 dM = createDM(mField.get_comm(), "DMMOFEM");
1792 CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
1793 BitRefLevel().set());
1794 CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
1795 CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1801
1802 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
1803 CHKERR DMSetUp(dM);
1804 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
1805
1806 auto remove_dofs_on_broken_skin = [&](const std::string prb_name) {
1808 for (int d : {0, 1, 2}) {
1809 std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
1811 ->getSideDofsOnBrokenSpaceEntities(
1812 dofs_to_remove, prb_name, ROW, piolaStress,
1814 // remove piola dofs, i.e. traction free boundary
1815 CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, ROW,
1816 dofs_to_remove);
1817 CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, COL,
1818 dofs_to_remove);
1819 }
1821 };
1822 CHKERR remove_dofs_on_broken_skin("ESHELBY_PLASTICITY");
1823
1824 // Create elastic sub-problem
1825 dmElastic = createDM(mField.get_comm(), "DMMOFEM");
1826 CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
1832 if (!noStretch) {
1834 }
1844 CHKERR DMSetUp(dmElastic);
1845
1846 // dmMaterial = createDM(mField.get_comm(), "DMMOFEM");
1847 // CHKERR DMMoFEMCreateSubDM(dmMaterial, dM, "MATERIAL_PROBLEM");
1848 // CHKERR DMMoFEMSetDestroyProblem(dmMaterial, PETSC_TRUE);
1849 // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, eshelbyStress);
1850 // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, materialL2Disp);
1851 // CHKERR DMMoFEMAddElement(dmMaterh elementVolumeName);
1852 // CHKERR DMMoFEMAddElement(dmMaterial, naturalBcElement);
1853 // CHKERR DMMoFEMAddElement(dmMaterial, skinElement);
1854 // CHKERR DMMoFEMSetSquareProblem(dmMaterial, PETSC_TRUE);
1855 // CHKERR DMMoFEMSetIsPartitioned(dmMaterial, PETSC_TRUE);
1856 // CHKERR DMSetUp(dmMaterial);
1857
1858 auto set_zero_block = [&]() {
1860 if (!noStretch) {
1861 CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1862 "ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
1863 CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1864 "ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
1865 }
1866 CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1867 "ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
1868 CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1869 "ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
1870 CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1871 "ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
1872 CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1873 "ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
1874 if (!noStretch) {
1875 CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1876 "ELASTIC_PROBLEM", bubbleField, bubbleField);
1877 CHKERR
1878 mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1879 "ELASTIC_PROBLEM", piolaStress, piolaStress);
1880 CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1881 "ELASTIC_PROBLEM", bubbleField, piolaStress);
1882 CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1883 "ELASTIC_PROBLEM", piolaStress, bubbleField);
1884 }
1885
1888 };
1889
1890 auto set_section = [&]() {
1892 PetscSection section;
1893 CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
1894 &section);
1895 CHKERR DMSetSection(dmElastic, section);
1896 CHKERR DMSetGlobalSection(dmElastic, section);
1897 CHKERR PetscSectionDestroy(&section);
1899 };
1900
1901 CHKERR set_zero_block();
1902 CHKERR set_section();
1903
1904 dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
1905 CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
1911 CHKERR DMSetUp(dmPrjSpatial);
1912
1913 // CHKERR mField.getInterface<BcManager>()
1914 // ->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
1915 // "PROJECT_SPATIAL", spatialH1Disp, true, false);
1916
1918}
@ QUIET
@ COL
@ ROW
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1113
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:215
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition DMMoFEM.cpp:114
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1234
auto bit
set bit
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition DMMoFEM.cpp:434
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
const std::string spatialL2Disp
SmartPetscObj< DM > dM
Coupled problem all fields.
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTractionVecPtr
const std::string elementVolumeName
const std::string bubbleField
static PetscBool noStretch
const std::string rotAxis
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
SmartPetscObj< Vec > solTSStep
SmartPetscObj< DM > dmElastic
Elastic problem.
const std::string stretchTensor
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
Problem manager is used to build and partition problems.

◆ addFields()

MoFEMErrorCode EshelbianCore::addFields ( const EntityHandle  meshset = 0)
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 1157 of file EshelbianPlasticity.cpp.

1157 {
1159
1160 auto get_tets = [&]() {
1161 Range tets;
1162 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
1163 return tets;
1164 };
1165
1166 auto get_tets_skin = [&]() {
1167 Range tets_skin_part;
1168 Skinner skin(&mField.get_moab());
1169 CHKERR skin.find_skin(0, get_tets(), false, tets_skin_part);
1170 ParallelComm *pcomm =
1171 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1172 Range tets_skin;
1173 CHKERR pcomm->filter_pstatus(tets_skin_part,
1174 PSTATUS_SHARED | PSTATUS_MULTISHARED,
1175 PSTATUS_NOT, -1, &tets_skin);
1176 return tets_skin;
1177 };
1178
1179 auto subtract_boundary_conditions = [&](auto &&tets_skin) {
1180 // That mean, that hybrid field on all faces on which traction is applied,
1181 // on other faces, or enforcing displacements as
1182 // natural boundary condition.
1184 for (auto &v : *bcSpatialTractionVecPtr) {
1185 tets_skin = subtract(tets_skin, v.faces);
1186 }
1188 for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
1189 tets_skin = subtract(tets_skin, v.faces);
1190 }
1191
1193 for (auto &v : *bcSpatialPressureVecPtr) {
1194 tets_skin = subtract(tets_skin, v.faces);
1195 }
1196
1197 return tets_skin;
1198 };
1199
1200 auto add_blockset = [&](auto block_name, auto &&tets_skin) {
1201 auto crack_faces =
1202 get_range_from_block(mField, "block_name", SPACE_DIM - 1);
1203 tets_skin.merge(crack_faces);
1204 return tets_skin;
1205 };
1206
1207 auto subtract_blockset = [&](auto block_name, auto &&tets_skin) {
1208 auto contact_range =
1209 get_range_from_block(mField, block_name, SPACE_DIM - 1);
1210 tets_skin = subtract(tets_skin, contact_range);
1211 return tets_skin;
1212 };
1213
1214 auto get_stress_trace_faces = [&](auto &&tets_skin) {
1215 Range faces;
1216 CHKERR mField.get_moab().get_adjacencies(get_tets(), SPACE_DIM - 1, true,
1217 faces, moab::Interface::UNION);
1218 Range trace_faces = subtract(faces, tets_skin);
1219 return trace_faces;
1220 };
1221
1222 auto tets = get_tets();
1223
1224 // remove also contact faces, i.e. that is also kind of hybrid field but
1225 // named but used to enforce contact conditions
1226 auto trace_faces = get_stress_trace_faces(
1227
1228 subtract_blockset("CONTACT",
1229 subtract_boundary_conditions(get_tets_skin()))
1230
1231 );
1232
1233 contactFaces = boost::make_shared<Range>(intersect(
1234 trace_faces, get_range_from_block(mField, "CONTACT", SPACE_DIM - 1)));
1236 boost::make_shared<Range>(subtract(trace_faces, *contactFaces));
1237 // materialSkeletonFaces =
1238 // boost::make_shared<Range>(get_stress_trace_faces(Range()));
1239
1240#ifndef NDEBUG
1241 if (contactFaces->size())
1243 "contact_faces_" +
1244 std::to_string(mField.get_comm_rank()) + ".vtk",
1245 *contactFaces);
1246 if (skeletonFaces->size())
1248 "skeleton_faces_" +
1249 std::to_string(mField.get_comm_rank()) + ".vtk",
1250 *skeletonFaces);
1251 // if (skeletonFaces->size())
1252 // CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
1253 // *materialSkeletonFaces);
1254#endif
1255
1256 auto add_broken_hdiv_field = [this, meshset](const std::string field_name,
1257 const int order) {
1259
1261
1262 auto get_side_map_hdiv = [&]() {
1263 return std::vector<
1264
1265 std::pair<EntityType,
1267
1268 >>{
1269
1270 {MBTET,
1271 [&](BaseFunction::DofsSideMap &dofs_side_map) -> MoFEMErrorCode {
1272 return TetPolynomialBase::setDofsSideMap(HDIV, DISCONTINUOUS, base,
1273 dofs_side_map);
1274 }}
1275
1276 };
1277 };
1278
1280 get_side_map_hdiv(), MB_TAG_DENSE, MF_ZERO);
1282 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
1284 };
1285
1286 auto add_l2_field = [this, meshset](const std::string field_name,
1287 const int order, const int dim) {
1290 MB_TAG_DENSE, MF_ZERO);
1292 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
1294 };
1295
1296 auto add_h1_field = [this, meshset](const std::string field_name,
1297 const int order, const int dim) {
1300 MB_TAG_DENSE, MF_ZERO);
1302 CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
1303 CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
1304 CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
1305 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
1307 };
1308
1309 auto add_l2_field_by_range = [this](const std::string field_name,
1310 const int order, const int dim,
1311 const int field_dim, Range &&r) {
1314 MB_TAG_DENSE, MF_ZERO);
1315 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(r);
1319 };
1320
1321 auto add_bubble_field = [this, meshset](const std::string field_name,
1322 const int order, const int dim) {
1324 CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
1325 MF_ZERO);
1326 // Modify field
1327 auto field_ptr = mField.get_field_structure(field_name);
1328 auto field_order_table =
1329 const_cast<Field *>(field_ptr)->getFieldOrderTable();
1330 auto get_cgg_bubble_order_zero = [](int p) { return 0; };
1331 auto get_cgg_bubble_order_tet = [](int p) {
1332 return NBVOLUMETET_CCG_BUBBLE(p);
1333 };
1334 field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
1335 field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
1336 field_order_table[MBTRI] = get_cgg_bubble_order_zero;
1337 field_order_table[MBTET] = get_cgg_bubble_order_tet;
1339 CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
1340 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
1342 };
1343
1344 auto add_user_l2_field = [this, meshset](const std::string field_name,
1345 const int order, const int dim) {
1347 CHKERR mField.add_field(field_name, L2, USER_BASE, dim, MB_TAG_DENSE,
1348 MF_ZERO);
1349 // Modify field
1350 auto field_ptr = mField.get_field_structure(field_name);
1351 auto field_order_table =
1352 const_cast<Field *>(field_ptr)->getFieldOrderTable();
1353 auto zero_dofs = [](int p) { return 0; };
1354 auto dof_l2_tet = [](int p) { return NBVOLUMETET_L2(p); };
1355 field_order_table[MBVERTEX] = zero_dofs;
1356 field_order_table[MBEDGE] = zero_dofs;
1357 field_order_table[MBTRI] = zero_dofs;
1358 field_order_table[MBTET] = dof_l2_tet;
1360 CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
1362 };
1363
1364 // spatial fields
1365 CHKERR add_broken_hdiv_field(piolaStress, spaceOrder);
1366 CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
1367 CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
1368 CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
1369 CHKERR add_user_l2_field(stretchTensor, noStretch ? -1 : spaceOrder, 6);
1370
1371 if (!skeletonFaces)
1372 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
1373 if (!contactFaces)
1374 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No contact faces");
1375
1376 auto get_hybridised_disp = [&]() {
1377 auto faces = *skeletonFaces;
1378 auto skin = subtract_boundary_conditions(get_tets_skin());
1379 for (auto &bc : *bcSpatialNormalDisplacementVecPtr) {
1380 faces.merge(intersect(bc.faces, skin));
1381 }
1382 return faces;
1383 };
1384
1385 CHKERR add_l2_field_by_range(hybridSpatialDisp, spaceOrder - 1, 2, 3,
1386 get_hybridised_disp());
1387 CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
1389
1390 // spatial displacement
1391 CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
1392 // material positions
1393 CHKERR add_h1_field(materialH1Positions, 2, 3);
1394
1395 // Eshelby stress
1396 // CHKERR add_broken_hdiv_field(eshelbyStress, spaceOrder);
1397 // CHKERR add_l2_field(materialL2Disp, spaceOrder - 1, 3);
1398 // CHKERR add_l2_field_by_range(hybridMaterialDisp, spaceOrder - 1, 2, 3,
1399 // Range(*materialSkeletonFaces));
1400
1402
1404}
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ USER_BASE
user implemented approximation base
Definition definitions.h:68
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
@ L2
field with C-1 continuity
Definition definitions.h:88
@ H1
continuous field
Definition definitions.h:85
@ HDIV
field with continuous normal traction
Definition definitions.h:87
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
constexpr int order
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
int r
Definition sdf.py:8
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
virtual MoFEMErrorCode add_broken_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const std::vector< std::pair< EntityType, std::function< MoFEMErrorCode(BaseFunction::DofsSideMap &)> > > list_dof_side_map, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Field data structure for finite element approximation.

◆ addMaterial_Hencky()

MoFEMErrorCode EshelbianCore::addMaterial_Hencky ( double  E,
double  nu 
)
Examples
ep.cpp, and mofem/users_modules/eshelbian_plasticity/ep.cpp.

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

540 {
542 physicalEquations = boost::make_shared<HMHHencky>(mField, E, nu);
544}
boost::shared_ptr< PhysicalEquations > physicalEquations

◆ addMaterial_HMHHStVenantKirchhoff()

MoFEMErrorCode EshelbianCore::addMaterial_HMHHStVenantKirchhoff ( const int  tape,
const double  lambda,
const double  mu,
const double  sigma_y 
)
Examples
ep.cpp, and mofem/users_modules/eshelbian_plasticity/ep.cpp.

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

515 {
517 physicalEquations = boost::make_shared<HMHStVenantKirchhoff>(lambda, mu);
518 CHKERR physicalEquations->recordTape(tape, nullptr);
520}
static double lambda

◆ 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, and mofem/users_modules/eshelbian_plasticity/ep.cpp.

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

524 {
527 boost::make_shared<HMHPMooneyRivlinWriggersEq63>(alpha, beta, lambda);
528 CHKERR physicalEquations->recordTape(tape, nullptr);
530}

◆ addMaterial_HMHNeohookean()

MoFEMErrorCode EshelbianCore::addMaterial_HMHNeohookean ( const int  tape,
const double  c10,
const double  K 
)
Examples
ep.cpp, and mofem/users_modules/eshelbian_plasticity/ep.cpp.

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

533 {
535 physicalEquations = boost::make_shared<HMHNeohookean>(mField, c10, K);
536 CHKERR physicalEquations->recordTape(tape, nullptr);
538}

◆ addVolumeFiniteElement()

MoFEMErrorCode EshelbianCore::addVolumeFiniteElement ( const EntityHandle  meshset = 0)
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 1520 of file EshelbianPlasticity.cpp.

1520 {
1522
1523 // set finite element fields
1524 auto add_field_to_fe = [this](const std::string fe,
1525 const std::string field_name) {
1531 };
1532
1537
1538 CHKERR add_field_to_fe(elementVolumeName, piolaStress);
1539 CHKERR add_field_to_fe(elementVolumeName, bubbleField);
1540 if (!noStretch)
1541 CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
1542 CHKERR add_field_to_fe(elementVolumeName, rotAxis);
1543 CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
1544 CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
1545 CHKERR add_field_to_fe(elementVolumeName, contactDisp);
1548
1549 // build finite elements data structures
1551 }
1552
1553 // if (!mField.check_finite_element(materialVolumeElement)) {
1554
1555 // Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM - 2);
1556
1557 // Range front_elements;
1558 // for (auto l = 0; l != frontLayers; ++l) {
1559 // Range front_elements_layer;
1560 // CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM, true,
1561 // front_elements_layer,
1562 // moab::Interface::UNION);
1563 // front_elements.merge(front_elements_layer);
1564 // front_edges.clear();
1565 // CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
1566 // SPACE_DIM - 2, true,
1567 // front_edges,
1568 // moab::Interface::UNION);
1569 // }
1570
1571 // CHKERR mField.add_finite_element(materialVolumeElement, MF_ZERO);
1572 // CHKERR mField.add_ents_to_finite_element_by_type(front_elements, MBTET,
1573 // materialVolumeElement);
1574 // // CHKERR add_field_to_fe(materialVolumeElement, eshelbyStress);
1575 // // CHKERR add_field_to_fe(materialVolumeElement, materialL2Disp);
1576 // CHKERR mField.build_finite_elements(materialVolumeElement);
1577 // }
1578
1580}

◆ calculateCrackArea()

MoFEMErrorCode EshelbianCore::calculateCrackArea ( boost::shared_ptr< double area_ptr)
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 6975 of file EshelbianPlasticity.cpp.

6975 {
6977
6978 if (!area_ptr) {
6979 CHK_THROW_MESSAGE(MOFEM_INVALID_DATA, "area_ptr is null");
6980 }
6981
6982 int success;
6983 *area_ptr = 0;
6984 if (mField.get_comm_rank() == 0) {
6985 MOFEM_LOG("EP", Sev::inform) << "Calculate crack area";
6986 auto crack_faces = get_range_from_block(mField, "CRACK", SPACE_DIM - 1);
6987 for (auto f : crack_faces) {
6988 *area_ptr += mField.getInterface<Tools>()->getTriArea(f);
6989 }
6990 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0, mField.get_comm());
6991 } else {
6992 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0, mField.get_comm());
6993 }
6994 if (success != MPI_SUCCESS) {
6996 }
6998}
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_INVALID_DATA
Definition definitions.h:36

◆ calculateFaceMaterialForce()

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

Create element to integration faces energies

Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 4208 of file EshelbianPlasticity.cpp.

4208 {
4210
4211 constexpr bool debug = false;
4212
4213 auto get_tags_vec = [&](std::vector<std::pair<std::string, int>> names) {
4214 std::vector<Tag> tags;
4215 tags.reserve(names.size());
4216 auto create_and_clean = [&]() {
4218 for (auto n : names) {
4219 tags.push_back(Tag());
4220 auto &tag = tags.back();
4221 auto &moab = mField.get_moab();
4222 auto rval = moab.tag_get_handle(n.first.c_str(), tag);
4223 if (rval == MB_SUCCESS) {
4224 moab.tag_delete(tag);
4225 }
4226 double def_val[] = {0., 0., 0.};
4227 CHKERR moab.tag_get_handle(n.first.c_str(), n.second, MB_TYPE_DOUBLE,
4228 tag, MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
4229 }
4231 };
4232 CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
4233 return tags;
4234 };
4235
4236 enum ExhangeTags { MATERIALFORCE, AREAGROWTH, GRIFFITHFORCE, FACEPRESSURE };
4237
4238 auto tags = get_tags_vec({{"MaterialForce", 3},
4239 {"AreaGrowth", 3},
4240 {"GriffithForce", 1},
4241 {"FacePressure", 1}});
4242
4243 auto calculate_material_forces = [&]() {
4245
4246 /**
4247 * @brief Create element to integration faces energies
4248 */
4249 auto get_face_material_force_fe = [&]() {
4251 auto fe_ptr = boost::make_shared<FaceEle>(mField);
4252 fe_ptr->getRuleHook = [](int, int, int) { return -1; };
4253 fe_ptr->setRuleHook =
4254 SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
4255
4256 // hybrid disp, evalated on face first
4257 EshelbianPlasticity::AddHOOps<2, 2, 3>::add(
4258 fe_ptr->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
4259 fe_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4260 hybridSpatialDisp, dataAtPts->getHybridDispAtPts()));
4261 fe_ptr->getOpPtrVector().push_back(
4263 hybridSpatialDisp, dataAtPts->getGradHybridDispAtPts()));
4264 auto op_loop_domain_side =
4266 mField, elementVolumeName, SPACE_DIM, Sev::noisy);
4267 fe_ptr->getOpPtrVector().push_back(op_loop_domain_side);
4268 fe_ptr->getOpPtrVector().push_back(new OpFaceMaterialForce(dataAtPts));
4269
4270 // evaluated in side domain, that is op_loop_domain_side
4271 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
4272 boost::make_shared<CGGUserPolynomialBase>();
4273
4274 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
4275 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
4276 materialH1Positions, frontAdjEdges, nullptr, nullptr, nullptr);
4277 op_loop_domain_side->getOpPtrVector().push_back(
4279 piolaStress, dataAtPts->getApproxPAtPts()));
4280 op_loop_domain_side->getOpPtrVector().push_back(
4282 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
4283 // op_loop_domain_side->getOpPtrVector().push_back(
4284 // new OpZeroNoSideBrokenDofs(piolaStress));
4285 // op_loop_domain_side->getOpPtrVector().push_back(
4286 // new OpCalculateHVecTensorGradient<3, 9, 3>(
4287 // piolaStress, dataAtPts->getGradPAtPts()));
4288 // op_loop_domain_side->getOpPtrVector().push_back(
4289 // new OpCalculateHVecTensorGradient<9, 9, 3>(
4290 // bubbleField, dataAtPts->getGradPAtPts(), MBMAXTYPE));
4291
4292 op_loop_domain_side->getOpPtrVector().push_back(
4294 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
4295 if (noStretch) {
4296 op_loop_domain_side->getOpPtrVector().push_back(
4297 physicalEquations->returnOpCalculateStretchFromStress(
4299 } else {
4300 op_loop_domain_side->getOpPtrVector().push_back(
4302 stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
4303 }
4304
4305 op_loop_domain_side->getOpPtrVector().push_back(
4307
4308 return fe_ptr;
4309 };
4310
4311 auto integrate_face_material_force_fe = [&](auto &&face_energy_fe) {
4314 dM, skeletonElement, face_energy_fe, 0, mField.get_comm_size());
4315
4316 auto face_exchange = CommInterface::createEntitiesPetscVector(
4317 mField.get_comm(), mField.get_moab(), 2, 3, Sev::inform);
4318
4319 auto print_loc_size = [this](auto v, auto str, auto sev) {
4321 int size;
4322 CHKERR VecGetLocalSize(v.second, &size);
4323 int low, high;
4324 CHKERR VecGetOwnershipRange(v.second, &low, &high);
4325 MOFEM_LOG("EPSYNC", sev) << str << " local size " << size << " ( "
4326 << low << " " << high << " ) ";
4329 };
4330 CHKERR print_loc_size(face_exchange, "material face_exchange",
4331 Sev::verbose);
4332
4334 mField.get_moab(), face_exchange, tags[ExhangeTags::MATERIALFORCE]);
4336 mField.get_moab(), faceExchange, tags[ExhangeTags::FACEPRESSURE]);
4337
4338 // #ifndef NDEBUG
4339 if (debug) {
4341 "front_skin_faces_material_force_" +
4342 std::to_string(mField.get_comm_rank()) + ".vtk",
4343 *skeletonFaces);
4344 }
4345 // #endif
4346
4348 };
4349
4350 CHKERR integrate_face_material_force_fe(get_face_material_force_fe());
4351
4353 };
4354
4355 auto calculate_front_material_force = [&]() {
4358
4359 auto get_conn = [&](auto e) {
4360 Range conn;
4361 CHK_MOAB_THROW(mField.get_moab().get_connectivity(&e, 1, conn, true),
4362 "get connectivity");
4363 return conn;
4364 };
4365
4366 auto get_conn_range = [&](auto e) {
4367 Range conn;
4368 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
4369 "get connectivity");
4370 return conn;
4371 };
4372
4373 auto get_adj = [&](auto e, auto dim) {
4374 Range adj;
4375 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(&e, 1, dim, true, adj),
4376 "get adj");
4377 return adj;
4378 };
4379
4380 auto get_adj_range = [&](auto e, auto dim) {
4381 Range adj;
4382 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(e, dim, true, adj,
4383 moab::Interface::UNION),
4384 "get adj");
4385 return adj;
4386 };
4387
4388 auto get_material_force = [&](auto r, auto th) {
4389 MatrixDouble material_forces(r.size(), 3, false);
4391 mField.get_moab().tag_get_data(th, r, material_forces.data().data()),
4392 "get data");
4393 return material_forces;
4394 };
4395
4396 if (mField.get_comm_rank() == 0) {
4397
4398 auto crack_edges = get_adj_range(*crackFaces, 1);
4399 auto front_nodes = get_conn_range(*frontEdges);
4400 auto body_edges = get_range_from_block(mField, "EDGES", 1);
4401 // auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
4402 // front_block_edges = subtract(front_block_edges, crack_edges);
4403 // auto front_block_edges_conn = get_conn_range(front_block_edges);
4404
4405 // #ifndef NDEBUG
4406 Range all_skin_faces;
4407 Range all_front_faces;
4408 // #endif
4409
4410 auto calculate_edge_direction = [&](auto e) {
4411 const EntityHandle *conn;
4412 int num_nodes;
4414 mField.get_moab().get_connectivity(e, conn, num_nodes, true),
4415 "get connectivity");
4416 std::array<double, 6> coords;
4418 mField.get_moab().get_coords(conn, num_nodes, coords.data()),
4419 "get coords");
4421 &coords[0], &coords[1], &coords[2]};
4423 &coords[3], &coords[4], &coords[5]};
4426 t_dir(i) = t_p1(i) - t_p0(i);
4427 return t_dir;
4428 };
4429
4430 // take bubble tets at node, and then avarage over the edges
4431 auto calculate_force_through_node = [&]() {
4433
4438
4439 Range body_ents;
4440 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
4441 body_ents);
4442 auto body_skin = get_skin(mField, body_ents);
4443 auto body_skin_conn = get_conn_range(body_skin);
4444
4445 // calculate nodal material force
4446 for (auto n : front_nodes) {
4447 auto adj_tets = get_adj(n, 3);
4448 for (int ll = 0; ll < nbJIntegralLevels; ++ll) {
4449 auto conn = get_conn_range(adj_tets);
4450 adj_tets = get_adj_range(conn, 3);
4451 }
4452 auto skin_faces = get_skin(mField, adj_tets);
4453 auto material_forces = get_material_force(skin_faces, tags[0]);
4454
4455 // #ifndef NDEBUG
4456 if (debug) {
4457 all_skin_faces.merge(skin_faces);
4458 }
4459 // #endif
4460
4461 auto calculate_node_material_force = [&]() {
4462 auto t_face_T =
4463 getFTensor1FromPtr<SPACE_DIM>(material_forces.data().data());
4464 FTensor::Tensor1<double, SPACE_DIM> t_node_force{0., 0., 0.};
4465 for (auto face : skin_faces) {
4466
4467 FTensor::Tensor1<double, SPACE_DIM> t_face_force_tmp{0., 0., 0.};
4468 t_face_force_tmp(I) = t_face_T(I);
4469 ++t_face_T;
4470
4471 auto face_tets = intersect(get_adj(face, 3), adj_tets);
4472
4473 if (face_tets.empty()) {
4474 continue;
4475 }
4476
4477 if (face_tets.size() != 1) {
4479 "face_tets.size() != 1");
4480 }
4481
4482 int side_number, sense, offset;
4483 CHK_MOAB_THROW(mField.get_moab().side_number(face_tets[0], face,
4484 side_number, sense,
4485 offset),
4486 "moab side number");
4487 t_face_force_tmp(I) *= sense;
4488 t_node_force(I) += t_face_force_tmp(I);
4489 }
4490
4491 return t_node_force;
4492 };
4493
4494 auto calculate_crack_area_growth_direction =
4495 [&](auto n, auto &t_node_force) {
4496 // if skin is on body surface, project the direction on it
4497 FTensor::Tensor1<double, SPACE_DIM> t_project{0., 0., 0.};
4498 auto boundary_node = intersect(Range(n, n), body_skin_conn);
4499 if (boundary_node.size()) {
4500 auto faces = intersect(get_adj(n, 2), body_skin);
4501 for (auto f : faces) {
4502 FTensor::Tensor1<double, 3> t_normal_face;
4503 CHKERR mField.getInterface<Tools>()->getTriNormal(
4504 f, &t_normal_face(0));
4505 t_project(I) += t_normal_face(I);
4506 }
4507 t_project.normalize();
4508 }
4509
4510 // calculate surface projection matrix
4513 t_Q(I, J) = t_kd(I, J);
4514 if (boundary_node.size()) {
4515 t_Q(I, J) -= t_project(I) * t_project(J);
4516 }
4517
4518 auto adj_faces = intersect(get_adj(n, 2), *crackFaces);
4519 if (adj_faces.empty()) {
4520 auto adj_edges =
4521 intersect(get_adj(n, 1), unite(*frontEdges, body_edges));
4522 double l = 0;
4523 for (auto e : adj_edges) {
4524 auto t_dir = calculate_edge_direction(e);
4525 l += t_dir.l2();
4526 }
4527 l /= 2;
4528 FTensor::Tensor1<double, SPACE_DIM> t_area_dir{0., 0., 0.};
4529 FTensor::Tensor1<double, SPACE_DIM> t_node_force_tmp;
4530 t_node_force_tmp(I) = t_node_force(I);
4531 t_node_force_tmp.normalize();
4532 t_area_dir(I) = -t_node_force_tmp(I);
4533 t_area_dir(I) *= l / 2;
4534 return t_area_dir;
4535 }
4536
4537 // calculate direction
4538 auto front_edges = get_adj(n, 1);
4539 FTensor::Tensor1<double, 3> t_area_dir{0., 0., 0.};
4540 for (auto f : adj_faces) {
4541 int num_nodes;
4542 const EntityHandle *conn;
4543 CHKERR mField.get_moab().get_connectivity(f, conn, num_nodes,
4544 true);
4545 std::array<double, 9> coords;
4546 CHKERR mField.get_moab().get_coords(conn, num_nodes,
4547 coords.data());
4548 FTensor::Tensor1<double, 3> t_face_normal;
4550 CHKERR mField.getInterface<Tools>()->getTriNormal(
4551 coords.data(), &t_face_normal(0), &t_d_normal(0, 0, 0));
4552 auto n_it = std::find(conn, conn + num_nodes, n);
4553 auto n_index = std::distance(conn, n_it);
4554
4555 FTensor::Tensor2<double, 3, 3> t_face_hessian{
4556 t_d_normal(0, n_index, 0), t_d_normal(0, n_index, 1),
4557 t_d_normal(0, n_index, 2),
4558
4559 t_d_normal(1, n_index, 0), t_d_normal(1, n_index, 1),
4560 t_d_normal(1, n_index, 2),
4561
4562 t_d_normal(2, n_index, 0), t_d_normal(2, n_index, 1),
4563 t_d_normal(2, n_index, 2)};
4564
4565 FTensor::Tensor2<double, 3, 3> t_projected_hessian;
4566 t_projected_hessian(I, J) =
4567 t_Q(I, K) * (t_face_hessian(K, L) * t_Q(L, J));
4568 t_face_normal.normalize();
4569 t_area_dir(K) +=
4570 t_face_normal(I) * t_projected_hessian(I, K) / 2.;
4571 }
4572
4573 return t_area_dir;
4574 };
4575
4576 auto t_node_force = calculate_node_material_force();
4577 t_node_force(I) /= griffithEnergy; // scale all by griffith energy
4579 mField.get_moab().tag_set_data(tags[ExhangeTags::MATERIALFORCE],
4580 &n, 1, &t_node_force(0)),
4581 "set data");
4582
4583 auto get_area_dir = [&]() {
4584 auto adj_faces = get_adj_range(adj_tets, 2);
4585 auto t_area_dir = calculate_crack_area_growth_direction(
4586 n, t_node_force);
4587 if (nbJIntegralLevels) {
4588 auto adj_edges = intersect(get_adj_range(adj_tets, 1),
4589 unite(*frontEdges, body_edges));
4590 auto seed_n = get_conn_range(adj_edges);
4591 auto skin_adj_edges = get_skin(mField, adj_edges);
4592 skin_adj_edges = subtract(skin_adj_edges, body_skin_conn);
4593 seed_n = subtract(seed_n, skin_adj_edges);
4594 t_area_dir(I) = 0;
4595 for (auto sn : seed_n) {
4596 auto t_area_dir_sn = calculate_crack_area_growth_direction(
4597 sn, t_node_force);
4598 t_area_dir(I) += t_area_dir_sn(I);
4599 }
4600 for (auto sn : skin_adj_edges) {
4601 auto t_area_dir_sn = calculate_crack_area_growth_direction(
4602 sn, t_node_force);
4603 t_area_dir(I) += t_area_dir_sn(I) / 2;
4604 }
4605 }
4606
4607 return t_area_dir;
4608 };
4609
4610 auto t_area_dir = get_area_dir();
4611
4613 mField.get_moab().tag_set_data(tags[ExhangeTags::AREAGROWTH], &n,
4614 1, &t_area_dir(0)),
4615 "set data");
4616 auto griffith = -t_node_force(I) * t_area_dir(I) /
4617 (t_area_dir(K) * t_area_dir(K));
4619 mField.get_moab().tag_set_data(tags[ExhangeTags::GRIFFITHFORCE],
4620 &n, 1, &griffith),
4621 "set data");
4622 }
4623
4624 // iterate over edges, and calculate average edge material force
4625 auto ave_node_force = [&](auto th) {
4627
4628 for (auto e : *frontEdges) {
4629
4630 auto conn = get_conn(e);
4631 auto data = get_material_force(conn, th);
4632 auto t_node = getFTensor1FromPtr<SPACE_DIM>(data.data().data());
4633 FTensor::Tensor1<double, SPACE_DIM> t_edge{0., 0., 0.};
4634 for (auto n : conn) {
4635 t_edge(I) += t_node(I);
4636 ++t_node;
4637 }
4638 t_edge(I) /= conn.size();
4639
4640 FTensor::Tensor1<double, SPACE_DIM> t_edge_direction =
4641 calculate_edge_direction(e);
4642 t_edge_direction.normalize();
4643
4648 t_cross(K) =
4649 FTensor::levi_civita(I, J, K) * t_edge_direction(I) * t_edge(J);
4650 t_edge(K) = FTensor::levi_civita(I, J, K) * t_edge_direction(J) *
4651 t_cross(I);
4652
4653 CHKERR mField.get_moab().tag_set_data(th, &e, 1, &t_edge(0));
4654 }
4656 };
4657
4658 // iterate over edges, and calculate average edge griffith energy
4659 auto ave_node_griffith_energy = [&](auto th) {
4661 for (auto e : *frontEdges) {
4663 CHKERR mField.get_moab().tag_get_data(
4664 tags[ExhangeTags::MATERIALFORCE], &e, 1, &t_edge_force(0));
4666 CHKERR mField.get_moab().tag_get_data(tags[ExhangeTags::AREAGROWTH],
4667 &e, 1, &t_edge_area_dir(0));
4668 double griffith_energy = -t_edge_force(I) * t_edge_area_dir(I) /
4669 (t_edge_area_dir(K) * t_edge_area_dir(K));
4670 CHKERR mField.get_moab().tag_set_data(th, &e, 1, &griffith_energy);
4671 }
4673 };
4674
4675 CHKERR ave_node_force(tags[ExhangeTags::MATERIALFORCE]);
4676 CHKERR ave_node_force(tags[ExhangeTags::AREAGROWTH]);
4677 CHKERR ave_node_griffith_energy(tags[ExhangeTags::GRIFFITHFORCE]);
4678
4680 };
4681
4682 CHKERR calculate_force_through_node();
4683
4684 // calculate face cross
4685 for (auto e : *frontEdges) {
4686 auto adj_faces = get_adj(e, 2);
4687 auto crack_face = intersect(get_adj(e, 2), *crackFaces);
4688
4689 // #ifndef NDEBUG
4690 if (debug) {
4691 all_front_faces.merge(adj_faces);
4692 }
4693 // #endif
4694
4696 CHKERR mField.get_moab().tag_get_data(tags[ExhangeTags::MATERIALFORCE],
4697 &e, 1, &t_edge_force(0));
4698 FTensor::Tensor1<double, SPACE_DIM> t_edge_direction =
4699 calculate_edge_direction(e);
4700 t_edge_direction.normalize();
4705 t_cross(K) = FTensor::levi_civita(I, J, K) * t_edge_direction(I) *
4706 t_edge_force(J);
4707
4708 for (auto f : adj_faces) {
4710 CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
4711 t_normal.normalize();
4712 int side_number, sense, offset;
4713 CHKERR mField.get_moab().side_number(f, e, side_number, sense,
4714 offset);
4715 auto dot = -sense * t_cross(I) * t_normal(I);
4716 CHK_MOAB_THROW(mField.get_moab().tag_set_data(
4717 tags[ExhangeTags::GRIFFITHFORCE], &f, 1, &dot),
4718 "set data");
4719 }
4720 }
4721
4722 // #ifndef NDEBUG
4723 if (debug) {
4724 int ts_step;
4725 CHKERR TSGetStepNumber(ts, &ts_step);
4727 "front_edges_material_force_" +
4728 std::to_string(ts_step) + ".vtk",
4729 *frontEdges);
4731 "front_skin_faces_material_force_" +
4732 std::to_string(ts_step) + ".vtk",
4733 all_skin_faces);
4735 "front_faces_material_force_" +
4736 std::to_string(ts_step) + ".vtk",
4737 all_front_faces);
4738 }
4739 // #endif
4740 }
4741
4742 auto edge_exchange = CommInterface::createEntitiesPetscVector(
4743 mField.get_comm(), mField.get_moab(), 1, 3, Sev::inform);
4745 mField.get_moab(), edge_exchange, tags[ExhangeTags::MATERIALFORCE]);
4747 mField.get_moab(), edge_exchange, tags[ExhangeTags::AREAGROWTH]);
4749 mField.get_moab(), edgeExchange, tags[ExhangeTags::GRIFFITHFORCE]);
4750
4752 };
4753
4754 auto print_results = [&]() {
4756
4757 auto get_conn_range = [&](auto e) {
4758 Range conn;
4759 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
4760 "get connectivity");
4761 return conn;
4762 };
4763
4764 auto get_tag_data = [&](auto &ents, auto tag, auto dim) {
4765 std::vector<double> data(ents.size() * dim);
4766 CHK_MOAB_THROW(mField.get_moab().tag_get_data(tag, ents, data.data()),
4767 "get data");
4768 return data;
4769 };
4770
4771 if (mField.get_comm_rank() == 0) {
4772 auto at_nodes = [&]() {
4774 auto conn = get_conn_range(*frontEdges);
4775 auto material_force =
4776 get_tag_data(conn, tags[ExhangeTags::MATERIALFORCE], 3);
4777 auto area_growth = get_tag_data(conn, tags[ExhangeTags::AREAGROWTH], 3);
4778 auto griffith_force =
4779 get_tag_data(conn, tags[ExhangeTags::GRIFFITHFORCE], 1);
4780 std::vector<double> coords(conn.size() * 3);
4781 CHK_MOAB_THROW(mField.get_moab().get_coords(conn, coords.data()),
4782 "get coords");
4783 if (conn.size())
4784 MOFEM_LOG("EPSELF", Sev::inform) << "Material force at nodes";
4785 for (size_t i = 0; i < conn.size(); ++i) {
4786 MOFEM_LOG("EPSELF", Sev::inform)
4787 << "Node " << conn[i] << " coords "
4788 << coords[i * 3 + 0] << " " << coords[i * 3 + 1] << " "
4789 << coords[i * 3 + 2] << " material force "
4790 << material_force[i * 3 + 0] << " "
4791 << material_force[i * 3 + 1] << " "
4792 << material_force[i * 3 + 2] << " area growth "
4793 << area_growth[i * 3 + 0] << " " << area_growth[i * 3 + 1]
4794 << " " << area_growth[i * 3 + 2] << " griffith force "
4795 << griffith_force[i];
4796 }
4798 };
4799
4800 at_nodes();
4801 }
4803 };
4804
4805 CHKERR calculate_material_forces();
4806 CHKERR calculate_front_material_force();
4807 CHKERR print_results();
4808
4810}
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
constexpr auto t_kd
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:557
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
FTensor::Index< 'l', 3 > l
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
constexpr IntegrationType I
boost::shared_ptr< Range > frontAdjEdges
static double griffithEnergy
Griffith energy.
CommInterface::EntitiesPetscVector edgeExchange
static int nbJIntegralLevels
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< Range > frontVertices
CommInterface::EntitiesPetscVector faceExchange
static MoFEMErrorCode updateEntitiesPetscVector(moab::Interface &moab, EntitiesPetscVector &vec, Tag tag, UpdateGhosts update_gosts=defaultUpdateGhosts)
Exchange data between vector and data.
static EntitiesPetscVector createEntitiesPetscVector(MPI_Comm comm, moab::Interface &moab, int dim, const int nb_coeffs, Sev sev=Sev::verbose, int root_rank=0)
Create a ghost vector for exchanging data.
virtual int get_comm_size() const =0
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for double precision vector field values calculation.
Element used to execute operators on side of the element.

◆ 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
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 4812 of file EshelbianPlasticity.cpp.

4813 {
4815
4816 constexpr bool debug = false;
4817 constexpr auto sev = Sev::verbose;
4818
4819 Range body_ents;
4820 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
4821 auto body_skin = get_skin(mField, body_ents);
4822 Range body_skin_edges;
4823 CHKERR mField.get_moab().get_adjacencies(body_skin, 1, false, body_skin_edges,
4824 moab::Interface::UNION);
4825 Range boundary_skin_verts;
4826 CHKERR mField.get_moab().get_connectivity(body_skin_edges,
4827 boundary_skin_verts, true);
4828
4829 auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
4830 Range geometry_edges_verts;
4831 CHKERR mField.get_moab().get_connectivity(geometry_edges,
4832 geometry_edges_verts, true);
4833 Range crack_faces_verts;
4834 CHKERR mField.get_moab().get_connectivity(*crackFaces, crack_faces_verts,
4835 true);
4836 Range crack_faces_edges;
4837 CHKERR mField.get_moab().get_adjacencies(
4838 *crackFaces, 1, true, crack_faces_edges, moab::Interface::UNION);
4839 Range crack_faces_tets;
4840 CHKERR mField.get_moab().get_adjacencies(
4841 *crackFaces, 3, true, crack_faces_tets, moab::Interface::UNION);
4842
4843 Range front_verts;
4844 CHKERR mField.get_moab().get_connectivity(*frontEdges, front_verts, true);
4845 Range front_faces;
4846 CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2, true, front_faces,
4847 moab::Interface::UNION);
4848 Range front_verts_edges;
4849 CHKERR mField.get_moab().get_adjacencies(
4850 front_verts, 1, true, front_verts_edges, moab::Interface::UNION);
4851
4852 auto get_tags_vec = [&](auto tag_name, int dim) {
4853 std::vector<Tag> tags(1);
4854
4855 if (dim > 3)
4857
4858 auto create_and_clean = [&]() {
4860 auto &moab = mField.get_moab();
4861 auto rval = moab.tag_get_handle(tag_name, tags[0]);
4862 if (rval == MB_SUCCESS) {
4863 moab.tag_delete(tags[0]);
4864 }
4865 double def_val[] = {0., 0., 0.};
4866 CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
4867 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
4869 };
4870
4871 CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
4872
4873 return tags;
4874 };
4875
4876 auto get_adj_front = [&](bool subtract_crack) {
4877 Range adj_front;
4878 CHKERR mField.get_moab().get_adjacencies(*frontEdges, SPACE_DIM - 1, true,
4879 adj_front, moab::Interface::UNION);
4880 if (subtract_crack)
4881 adj_front = subtract(adj_front, *crackFaces);
4882 return adj_front;
4883 };
4884
4885 MOFEM_LOG_CHANNEL("SELF");
4886
4887 auto th_front_position = get_tags_vec("FrontPosition", 3);
4888 auto th_max_face_energy = get_tags_vec("MaxFaceEnergy", 1);
4889
4890 if (mField.get_comm_rank() == 0) {
4891
4892 auto get_crack_adj_tets = [&](auto r) {
4893 Range crack_faces_conn;
4894 CHKERR mField.get_moab().get_connectivity(r, crack_faces_conn);
4895 Range crack_faces_conn_tets;
4896 CHKERR mField.get_moab().get_adjacencies(crack_faces_conn, SPACE_DIM,
4897 true, crack_faces_conn_tets,
4898 moab::Interface::UNION);
4899 return crack_faces_conn_tets;
4900 };
4901
4902 auto get_layers_for_sides = [&](auto &side) {
4903 std::vector<Range> layers;
4904 auto get = [&]() {
4906
4907 auto get_adj = [&](auto &r, int dim) {
4908 Range adj;
4909 CHKERR mField.get_moab().get_adjacencies(r, dim, true, adj,
4910 moab::Interface::UNION);
4911 return adj;
4912 };
4913
4914 auto get_tets = [&](auto r) { return get_adj(r, SPACE_DIM); };
4915
4916 Range front_nodes;
4917 CHKERR mField.get_moab().get_connectivity(*frontEdges, front_nodes,
4918 true);
4919 Range front_faces = get_adj(front_nodes, 2);
4920 front_faces = subtract(front_faces, *crackFaces);
4921 auto front_tets = get_tets(front_nodes);
4922 auto front_side = intersect(side, front_tets);
4923 layers.push_back(front_side);
4924 for (;;) {
4925 auto adj_faces = get_skin(mField, layers.back());
4926 adj_faces = intersect(adj_faces, front_faces);
4927 auto adj_faces_tets = get_tets(adj_faces);
4928 adj_faces_tets = intersect(adj_faces_tets, front_tets);
4929 layers.push_back(unite(layers.back(), adj_faces_tets));
4930 if (layers.back().size() == layers[layers.size() - 2].size()) {
4931 break;
4932 }
4933 }
4935 };
4936 CHK_THROW_MESSAGE(get(), "get_layers_for_sides");
4937 return layers;
4938 };
4939
4941 auto layers_top = get_layers_for_sides(sides_pair.first);
4942 auto layers_bottom = get_layers_for_sides(sides_pair.second);
4943
4944#ifndef NDEBUG
4945 if (debug) {
4947 mField.get_moab(),
4948 "crack_tets_" +
4949 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
4950 get_crack_adj_tets(*crackFaces));
4951 CHKERR save_range(mField.get_moab(), "sides_first.vtk", sides_pair.first);
4952 CHKERR save_range(mField.get_moab(), "sides_second.vtk",
4953 sides_pair.second);
4954 MOFEM_LOG("EP", sev) << "Nb. layers " << layers_top.size();
4955 int l = 0;
4956 for (auto &r : layers_top) {
4957 MOFEM_LOG("EP", sev) << "Layer " << l << " size " << r.size();
4959 mField.get_moab(),
4960 "layers_top_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
4961 ++l;
4962 }
4963
4964 l = 0;
4965 for (auto &r : layers_bottom) {
4966 MOFEM_LOG("EP", sev) << "Layer " << l << " size " << r.size();
4968 mField.get_moab(),
4969 "layers_bottom_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
4970 ++l;
4971 }
4972 }
4973#endif
4974
4975 auto get_cross = [&](auto t_dir, auto f) {
4977 CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
4978 t_normal.normalize();
4983 t_cross(i) = FTensor::levi_civita(i, j, k) * t_normal(j) * t_dir(k);
4984 return t_cross;
4985 };
4986
4987 auto get_sense = [&](auto f, auto e) {
4988 int side, sense, offset;
4989 CHK_MOAB_THROW(mField.get_moab().side_number(f, e, side, sense, offset),
4990 "get sense");
4991 return std::make_tuple(side, sense, offset);
4992 };
4993
4994 auto calculate_edge_direction = [&](auto e, auto normalize = true) {
4995 const EntityHandle *conn;
4996 int num_nodes;
4997 CHKERR mField.get_moab().get_connectivity(e, conn, num_nodes, true);
4998 std::array<double, 6> coords;
4999 CHKERR mField.get_moab().get_coords(conn, num_nodes, coords.data());
5001 &coords[0], &coords[1], &coords[2]};
5003 &coords[3], &coords[4], &coords[5]};
5006 t_dir(i) = t_p1(i) - t_p0(i);
5007 if (normalize)
5008 t_dir.normalize();
5009 return t_dir;
5010 };
5011
5012 auto evaluate_face_energy_and_set_orientation = [&](auto front_edges,
5013 auto front_faces,
5014 auto &sides_pair,
5015 auto th_position) {
5017
5018 Tag th_face_energy;
5019 Tag th_material_force;
5020 switch (energyReleaseSelector) {
5021 case GRIFFITH_FORCE:
5022 case GRIFFITH_SKELETON:
5023 CHKERR mField.get_moab().tag_get_handle("GriffithForce",
5024 th_face_energy);
5025 CHKERR mField.get_moab().tag_get_handle("MaterialForce",
5026 th_material_force);
5027 break;
5028 default:
5029 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
5030 "Unknown energy release selector");
5031 };
5032
5033 /**
5034 * Iterate over front edges, get adjacent faces, find maximal face energy.
5035 * Maximal face energy is stored in the edge. Maximal face energy is
5036 * magnitude of edge Griffith force.
5037 */
5038 auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
5039 auto &edge_face_max_energy_map) {
5041
5042 Range body_ents;
5043 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
5044 auto body_skin = get_skin(mField, body_ents);
5045
5046 Range max_faces;
5047
5048 for (auto e : front_edges) {
5049
5050 double griffith_force;
5051 CHKERR mField.get_moab().tag_get_data(th_face_energy, &e, 1,
5052 &griffith_force);
5053
5054 Range faces;
5055 CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
5056 faces = subtract(intersect(faces, front_faces), body_skin);
5057 std::vector<double> face_energy(faces.size());
5058 CHKERR mField.get_moab().tag_get_data(th_face_energy, faces,
5059 face_energy.data());
5060 auto max_energy_it =
5061 std::max_element(face_energy.begin(), face_energy.end());
5062 double max_energy =
5063 max_energy_it != face_energy.end() ? *max_energy_it : 0;
5064
5065 edge_face_max_energy_map[e] =
5066 std::make_tuple(faces[max_energy_it - face_energy.begin()],
5067 griffith_force, static_cast<double>(0));
5068 MOFEM_LOG("EP", Sev::inform)
5069 << "Edge " << e << " griffith force " << griffith_force
5070 << " max face energy " << max_energy << " factor "
5071 << max_energy / griffith_force;
5072
5073 max_faces.insert(faces[max_energy_it - face_energy.begin()]);
5074 }
5075
5076#ifndef NDEBUG
5077 if (debug) {
5079 mField.get_moab(),
5080 "max_faces_" +
5081 boost::lexical_cast<std::string>(mField.get_comm_rank()) +
5082 ".vtk",
5083 max_faces);
5084 }
5085#endif
5086
5088 };
5089
5090 /**
5091 * For each front edge, find maximal face energy and orientation. This is
5092 * by finding angle between edge material force and maximal face normal
5093 *
5094 */
5095 auto calculate_face_orientation = [&](auto &edge_face_max_energy_map) {
5097
5098 auto up_down_face = [&](
5099
5100 auto &face_angle_map_up,
5101 auto &face_angle_map_down
5102
5103 ) {
5105
5106 for (auto &m : edge_face_max_energy_map) {
5107 auto e = m.first;
5108 auto [max_face, energy, opt_angle] = m.second;
5109
5110 Range faces;
5111 CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
5112 faces = intersect(faces, front_faces);
5113 Range adj_tets; // tetrahedrons adjacent to the face
5114 CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
5115 false, adj_tets,
5116 moab::Interface::UNION);
5117 if (adj_tets.size()) {
5118
5119 Range adj_tets; // tetrahedrons adjacent to the face
5120 CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
5121 false, adj_tets,
5122 moab::Interface::UNION);
5123 if (adj_tets.size()) {
5124
5125 Range adj_tets_faces;
5126 // get faces
5127 CHKERR mField.get_moab().get_adjacencies(
5128 adj_tets, SPACE_DIM - 1, false, adj_tets_faces,
5129 moab::Interface::UNION);
5130 adj_tets_faces = intersect(adj_tets_faces, faces);
5132
5133 // cross product of face normal and edge direction
5134 auto t_cross_max =
5135 get_cross(calculate_edge_direction(e, true), max_face);
5136 auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
5137 t_cross_max(i) *= sense_max;
5138
5139 for (auto t : adj_tets) {
5140 Range adj_tets_faces;
5141 CHKERR mField.get_moab().get_adjacencies(
5142 &t, 1, SPACE_DIM - 1, false, adj_tets_faces);
5143 adj_tets_faces = intersect(adj_tets_faces, faces);
5144 adj_tets_faces =
5145 subtract(adj_tets_faces, Range(max_face, max_face));
5146
5147 if (adj_tets_faces.size() == 1) {
5148
5149 // cross product of adjacent face normal and edge
5150 // direction
5151 auto t_cross = get_cross(calculate_edge_direction(e, true),
5152 adj_tets_faces[0]);
5153 auto [side, sense, offset] =
5154 get_sense(adj_tets_faces[0], e);
5155 t_cross(i) *= sense;
5156 double dot = t_cross(i) * t_cross_max(i);
5157 auto angle = std::acos(dot);
5158
5159 double face_energy;
5160 CHKERR mField.get_moab().tag_get_data(
5161 th_face_energy, adj_tets_faces, &face_energy);
5162
5163 auto [side_face, sense_face, offset_face] =
5164 get_sense(t, max_face);
5165
5166 if (sense_face > 0) {
5167 face_angle_map_up[e] = std::make_tuple(face_energy, angle,
5168 adj_tets_faces[0]);
5169
5170 } else {
5171 face_angle_map_down[e] = std::make_tuple(
5172 face_energy, -angle, adj_tets_faces[0]);
5173 }
5174 }
5175 }
5176 }
5177 }
5178 }
5179
5181 };
5182
5183 auto calc_optimal_angle = [&](
5184
5185 auto &face_angle_map_up,
5186 auto &face_angle_map_down
5187
5188 ) {
5190
5191 for (auto &m : edge_face_max_energy_map) {
5192 auto e = m.first;
5193 auto &[max_face, e0, a0] = m.second;
5194
5195 if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
5196
5197 if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
5198 face_angle_map_down.find(e) == face_angle_map_down.end()) {
5199 // Do nothing
5200 } else {
5201
5202 switch (energyReleaseSelector) {
5203 case GRIFFITH_FORCE:
5204 case GRIFFITH_SKELETON: {
5205
5206 Tag th_material_force;
5207 CHKERR mField.get_moab().tag_get_handle("MaterialForce",
5208 th_material_force);
5209 FTensor::Tensor1<double, SPACE_DIM> t_material_force;
5210 CHKERR mField.get_moab().tag_get_data(
5211 th_material_force, &e, 1, &t_material_force(0));
5212 auto material_force_magnitude = t_material_force.l2();
5213 if (material_force_magnitude <
5214 std::numeric_limits<double>::epsilon()) {
5215 a0 = 0;
5216
5217 } else {
5218
5219 auto t_edge_dir = calculate_edge_direction(e, true);
5220 auto t_cross_max = get_cross(t_edge_dir, max_face);
5221 auto [side, sense, offset] = get_sense(max_face, e);
5222 t_cross_max(sense) *= sense;
5223
5227
5228 t_material_force.normalize();
5229 t_cross_max.normalize();
5231 t_cross(I) = FTensor::levi_civita(I, J, K) *
5232 t_material_force(J) * t_cross_max(K);
5233 a0 = -std::asin(t_cross(I) * t_edge_dir(I));
5234
5235 MOFEM_LOG("EP", sev)
5236 << "Optimal angle " << a0 << " energy " << e0;
5237 }
5238 break;
5239 }
5240 default: {
5241
5242 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
5243 "Unknown energy release selector");
5244 }
5245 }
5246 }
5247 }
5248 }
5249
5251 };
5252
5253 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
5254 face_angle_map_up;
5255 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
5256 face_angle_map_down;
5257 CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
5258 CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
5259
5260#ifndef NDEBUG
5261 if (debug) {
5262 auto th_angle = get_tags_vec("Angle", 1);
5263 Range up;
5264 for (auto &m : face_angle_map_up) {
5265 auto [e, a, face] = m.second;
5266 up.insert(face);
5267 CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
5268 }
5269 Range down;
5270 for (auto &m : face_angle_map_down) {
5271 auto [e, a, face] = m.second;
5272 down.insert(face);
5273 CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
5274 }
5275
5276 Range max_energy_faces;
5277 for (auto &m : edge_face_max_energy_map) {
5278 auto [face, e, angle] = m.second;
5279 max_energy_faces.insert(face);
5280 CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1,
5281 &angle);
5282 }
5283 if (mField.get_comm_rank() == 0) {
5284 CHKERR save_range(mField.get_moab(), "up_faces.vtk", up);
5285 CHKERR save_range(mField.get_moab(), "down_faces.vtk", down);
5286 CHKERR save_range(mField.get_moab(), "max_energy_faces.vtk",
5287 max_energy_faces);
5288 }
5289 }
5290#endif // NDEBUG
5291
5293 };
5294
5295 auto get_conn = [&](auto e) {
5296 Range conn;
5297 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
5298 "get conn");
5299 return conn;
5300 };
5301
5302 auto get_adj = [&](auto e, auto dim) {
5303 Range adj;
5304 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
5305 e, dim, false, adj, moab::Interface::UNION),
5306 "get adj");
5307 return adj;
5308 };
5309
5310 auto get_coords = [&](auto v) {
5312 CHK_MOAB_THROW(mField.get_moab().get_coords(v, &t_coords(0)),
5313 "get coords");
5314 return t_coords;
5315 };
5316
5317 // calulate normal of the max energy face
5318 auto get_rotated_normal = [&](auto e, auto f, auto angle) {
5321 auto t_edge_dir = calculate_edge_direction(e, true);
5322 auto [side, sense, offset] = get_sense(f, e);
5323 t_edge_dir(i) *= sense;
5324 t_edge_dir.normalize();
5325 t_edge_dir(i) *= angle;
5326 auto t_R = LieGroups::SO3::exp(t_edge_dir, angle);
5328 mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
5329 FTensor::Tensor1<double, SPACE_DIM> t_rotated_normal;
5330 t_rotated_normal(i) = t_R(i, j) * t_normal(j);
5331 return std::make_tuple(t_normal, t_rotated_normal);
5332 };
5333
5334 auto set_coord = [&](auto v, auto &adj_vertex_tets_verts, auto &coords,
5335 auto &t_move, auto gamma) {
5336 auto index = adj_vertex_tets_verts.index(v);
5337 if (index >= 0) {
5338 for (auto ii : {0, 1, 2}) {
5339 coords[3 * index + ii] += gamma * t_move(ii);
5340 }
5341 return true;
5342 }
5343 return false;
5344 };
5345
5346 auto tets_quality = [&](auto quality, auto &adj_vertex_tets_verts,
5347 auto &adj_vertex_tets, auto &coords) {
5348 for (auto t : adj_vertex_tets) {
5349 const EntityHandle *conn;
5350 int num_nodes;
5351 CHKERR mField.get_moab().get_connectivity(t, conn, num_nodes, true);
5352 std::array<double, 12> tet_coords;
5353 for (auto n = 0; n != 4; ++n) {
5354 auto index = adj_vertex_tets_verts.index(conn[n]);
5355 if (index < 0) {
5357 }
5358 for (auto ii = 0; ii != 3; ++ii) {
5359 tet_coords[3 * n + ii] = coords[3 * index + ii];
5360 }
5361 }
5362 double q = Tools::volumeLengthQuality(tet_coords.data());
5363 if (!std::isnormal(q))
5364 q = -2;
5365 quality = std::min(quality, q);
5366 };
5367
5368 return quality;
5369 };
5370
5371 auto calculate_free_face_node_displacement =
5372 [&](auto &edge_face_max_energy_map) {
5373 // get edges adjacent to vertex along which nodes are moving
5374 auto get_vertex_edges = [&](auto vertex) {
5375 Range vertex_edges; // edges adjacent to vertex
5376
5377 auto impl = [&]() {
5379 CHKERR mField.get_moab().get_adjacencies(vertex, 1, false,
5380 vertex_edges);
5381 vertex_edges = subtract(vertex_edges, front_verts_edges);
5382
5383 if (boundary_skin_verts.size() &&
5384 boundary_skin_verts.find(vertex[0]) !=
5385 boundary_skin_verts.end()) {
5386 MOFEM_LOG("EP", sev) << "Boundary vertex";
5387 vertex_edges = intersect(vertex_edges, body_skin_edges);
5388 }
5389 if (geometry_edges_verts.size() &&
5390 geometry_edges_verts.find(vertex[0]) !=
5391 geometry_edges_verts.end()) {
5392 MOFEM_LOG("EP", sev) << "Geometry edge vertex";
5393 vertex_edges = intersect(vertex_edges, geometry_edges);
5394 }
5395 if (crack_faces_verts.size() &&
5396 crack_faces_verts.find(vertex[0]) !=
5397 crack_faces_verts.end()) {
5398 MOFEM_LOG("EP", sev) << "Crack face vertex";
5399 vertex_edges = intersect(vertex_edges, crack_faces_edges);
5400 }
5402 };
5403
5404 CHK_THROW_MESSAGE(impl(), "get_vertex_edges");
5405
5406 return vertex_edges;
5407 };
5408
5409 // vector of rotated faces, edge along node is moved, moved edge,
5410 // moved displacement, quality, cardinality, gamma
5411 using Bundle = std::vector<
5412
5415
5416 >;
5417 std::map<EntityHandle, Bundle> edge_bundle_map;
5418
5419 for (auto &m : edge_face_max_energy_map) {
5420
5421 auto edge = m.first;
5422 auto &[max_face, energy, opt_angle] = m.second;
5423
5424 // calculate rotation of max energy face
5425 auto [t_normal, t_rotated_normal] =
5426 get_rotated_normal(edge, max_face, opt_angle);
5427
5428 auto front_vertex = get_conn(Range(m.first, m.first));
5429 auto adj_tets = get_adj(Range(max_face, max_face), 3);
5430 auto adj_tets_faces = get_adj(adj_tets, 2);
5431 auto adj_front_faces = subtract(
5432 intersect(get_adj(Range(edge, edge), 2), adj_tets_faces),
5433 *crackFaces);
5434 if (adj_front_faces.size() > 3)
5436 "adj_front_faces.size()>3");
5437
5438 FTensor::Tensor1<double, SPACE_DIM> t_material_force;
5439 CHKERR mField.get_moab().tag_get_data(th_material_force, &edge, 1,
5440 &t_material_force(0));
5441 std::vector<double> griffith_energy(adj_front_faces.size());
5442 CHKERR mField.get_moab().tag_get_data(
5443 th_face_energy, adj_front_faces, griffith_energy.data());
5444
5445
5446 auto set_edge_bundle = [&](auto min_gamma) {
5447 for (auto rotated_f : adj_front_faces) {
5448
5449 double rotated_face_energy =
5450 griffith_energy[adj_front_faces.index(rotated_f)];
5451
5452 auto vertex = subtract(get_conn(Range(rotated_f, rotated_f)),
5453 front_vertex);
5454 if (vertex.size() != 1) {
5456 "Wrong number of vertex to move");
5457 }
5458 auto front_vertex_edges_vertex = get_conn(
5459 intersect(get_adj(front_vertex, 1), crack_faces_edges));
5460 vertex = subtract(
5461 vertex, front_vertex_edges_vertex); // vertex free to move
5462 if (vertex.empty()) {
5463 continue;
5464 }
5465
5466 auto face_cardinality = [&](auto f, auto &seen_front_edges) {
5467 auto whole_front =
5468 unite(*frontEdges,
5469 subtract(body_skin_edges, crack_faces_edges));
5470 auto faces = Range(f, f);
5471 int c = 0;
5472 for (; c < 10; ++c) {
5473 auto front_edges =
5474 subtract(get_adj(faces, 1), seen_front_edges);
5475 if (front_edges.size() == 0) {
5476 return 0;
5477 }
5478 auto front_connected_edges =
5479 intersect(front_edges, whole_front);
5480 if (front_connected_edges.size()) {
5481 seen_front_edges.merge(front_connected_edges);
5482 return c;
5483 }
5484 faces.merge(get_adj(front_edges, 2));
5485 ++c;
5486 }
5487 return c;
5488 };
5489
5490 Range seen_edges = Range(edge, edge);
5491 double rotated_face_cardinality = face_cardinality(
5492 rotated_f,
5493 seen_edges); // add cardinality of max energy
5494 // face to rotated face cardinality
5495 // rotated_face_cardinality +=
5496 // face_cardinality(max_face, seen_edges);
5497 rotated_face_cardinality = std::max(rotated_face_cardinality,
5498 1.); // at least one edge
5499
5500 auto t_vertex_coords = get_coords(vertex);
5501 auto vertex_edges = get_vertex_edges(vertex);
5502
5503 EntityHandle f0 = front_vertex[0];
5504 EntityHandle f1 = front_vertex[1];
5505 FTensor::Tensor1<double, 3> t_v_e0, t_v_e1;
5506 CHKERR mField.get_moab().get_coords(&f0, 1, &t_v_e0(0));
5507 CHKERR mField.get_moab().get_coords(&f1, 1, &t_v_e1(0));
5508
5510 for (auto e_used_to_move_detection : vertex_edges) {
5511 auto edge_conn = get_conn(Range(e_used_to_move_detection,
5512 e_used_to_move_detection));
5513 edge_conn = subtract(edge_conn, vertex);
5514 // Find displacement of the edge such that dot porduct with
5515 // normal is zero.
5516 //
5517 // { (t_v0 - t_vertex_coords) + gamma * (t_v3 -
5518 // t_vertex_coords) } * n = 0
5519 // where t_v0 is the edge vertex, t_v3 is the edge end
5520 // point, n is the rotated normal of the face gamma is the
5521 // factor by which the edge is moved
5523 t_v0(i) = (t_v_e0(i) + t_v_e1(i)) / 2;
5525 CHKERR mField.get_moab().get_coords(edge_conn, &t_v3(0));
5526 auto a =
5527 (t_v0(i) - t_vertex_coords(i)) * t_rotated_normal(i);
5528 auto b =
5529 (t_v3(i) - t_vertex_coords(i)) * t_rotated_normal(i);
5530 auto gamma = a / b;
5531
5532 constexpr double eps =
5533 std::numeric_limits<double>::epsilon();
5534 if (std::isnormal(gamma) && gamma < 1.0 - eps &&
5535 gamma > -0.1) {
5537 t_move(i) = gamma * (t_v3(i) - t_vertex_coords(i));
5538
5539 auto check_rotated_face_directoon = [&]() {
5541 t_delta(i) = t_vertex_coords(i) + t_move(i) - t_v0(i);
5542 t_delta.normalize();
5543 auto dot =
5544 (t_material_force(i) / t_material_force.l2()) *
5545 t_delta(i);
5546 return -dot > 0 ? true : false;
5547 };
5548
5549 if (check_rotated_face_directoon()) {
5550
5551 MOFEM_LOG("EP", Sev::inform)
5552 << "Crack edge " << edge << " moved face "
5553 << rotated_f
5554 << " edge: " << e_used_to_move_detection
5555 << " face direction/energy " << rotated_face_energy
5556 << " face cardinality " << rotated_face_cardinality
5557 << " gamma: " << gamma;
5558
5559 auto &bundle = edge_bundle_map[edge];
5560 bundle.emplace_back(rotated_f, e_used_to_move_detection,
5561 vertex[0], t_move, 1,
5562 rotated_face_cardinality, gamma);
5563 }
5564 }
5565 }
5566 }
5567 };
5568
5569 set_edge_bundle(std::numeric_limits<double>::epsilon());
5570 if (edge_bundle_map[edge].empty()) {
5571 set_edge_bundle(-1.);
5572 }
5573 }
5574
5575 return edge_bundle_map;
5576 };
5577
5578 auto get_sort_by_energy = [&](auto &edge_face_max_energy_map) {
5579 std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
5580 sort_by_energy;
5581
5582 for (auto &m : edge_face_max_energy_map) {
5583 auto e = m.first;
5584 auto &[max_face, energy, opt_angle] = m.second;
5585 sort_by_energy[energy] = std::make_tuple(e, max_face, opt_angle);
5586 }
5587
5588 return sort_by_energy;
5589 };
5590
5591 auto set_tag = [&](auto &&adj_edges_map, auto &&sort_by_energy) {
5593
5594 Tag th_face_pressure;
5596 mField.get_moab().tag_get_handle("FacePressure", th_face_pressure),
5597 "get tag");
5598 auto get_face_pressure = [&](auto face) {
5599 double pressure;
5600 CHK_MOAB_THROW(mField.get_moab().tag_get_data(th_face_pressure, &face,
5601 1, &pressure),
5602 "get rag data");
5603 return pressure;
5604 };
5605
5606 MOFEM_LOG("EPSELF", Sev::inform)
5607 << "Number of edges to check " << sort_by_energy.size();
5608
5609 enum face_energy { POSITIVE, NEGATIVE };
5610 constexpr bool skip_negative = true;
5611
5612 for (auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
5613
5614 // iterate edges wih maximal energy, and make them seed. Such edges,
5615 // will most likely will have also smallest node displacement
5616 for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
5617 ++it) {
5618
5619 auto energy = it->first;
5620 auto [max_edge, max_face, opt_angle] = it->second;
5621
5622 auto face_pressure = get_face_pressure(max_face);
5623 if (skip_negative) {
5624 if (fe == face_energy::POSITIVE) {
5625 if (face_pressure < 0) {
5626 MOFEM_LOG("EPSELF", Sev::inform)
5627 << "Skip negative face " << max_face << " with energy "
5628 << energy << " and pressure " << face_pressure;
5629 continue;
5630 }
5631 }
5632 }
5633
5634 MOFEM_LOG("EPSELF", Sev::inform)
5635 << "Check face " << max_face << " edge " << max_edge
5636 << " energy " << energy << " optimal angle " << opt_angle
5637 << " face pressure " << face_pressure;
5638
5639 auto jt = adj_edges_map.find(max_edge);
5640 if (jt == adj_edges_map.end()) {
5641 MOFEM_LOG("EPSELF", Sev::warning)
5642 << "Edge " << max_edge << " not found in adj_edges_map";
5643 continue;
5644 }
5645 auto &bundle = jt->second;
5646
5647 auto find_max_in_bundle_impl = [&](auto edge, auto &bundle,
5648 auto gamma) {
5650
5651 EntityHandle vertex_max = 0;
5652 EntityHandle face_max = 0;
5653 EntityHandle move_edge_max = 0;
5654 double max_quality = -2;
5655 double max_quality_evaluated = -2;
5656 double min_cardinality = std::numeric_limits<double>::max();
5657
5658 FTensor::Tensor1<double, SPACE_DIM> t_move_last{0., 0., 0.};
5659
5660 for (auto &b : bundle) {
5661 auto &[face, move_edge, vertex, t_move, quality, cardinality,
5662 edge_gamma] = b;
5663
5664 auto adj_vertex_tets = get_adj(Range(vertex, vertex), 3);
5665 auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
5666 std::vector<double> coords(3 * adj_vertex_tets_verts.size());
5667 CHK_MOAB_THROW(mField.get_moab().get_coords(
5668 adj_vertex_tets_verts, coords.data()),
5669 "get coords");
5670
5671 set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
5672 quality = tets_quality(quality, adj_vertex_tets_verts,
5673 adj_vertex_tets, coords);
5674
5675 auto eval_quality = [](auto q, auto c, auto edge_gamma) {
5676 if (q < 0) {
5677 return q;
5678 } else {
5679 return ((edge_gamma < 0) ? (q / 2) : q) / pow(c, 2);
5680 }
5681 };
5682
5683 if (eval_quality(quality, cardinality, edge_gamma) >=
5684 max_quality_evaluated) {
5685 max_quality = quality;
5686 min_cardinality = cardinality;
5687 vertex_max = vertex;
5688 face_max = face;
5689 move_edge_max = move_edge;
5690 t_move_last(i) = t_move(i);
5691 max_quality_evaluated =
5692 eval_quality(max_quality, min_cardinality, edge_gamma);
5693 }
5694 }
5695
5696 return std::make_tuple(vertex_max, face_max, t_move_last,
5697 max_quality, min_cardinality);
5698 };
5699
5700 auto find_max_in_bundle = [&](auto edge, auto &bundle) {
5701 auto b_org_bundle = bundle;
5702 auto r = find_max_in_bundle_impl(edge, bundle, 1.);
5703 auto &[vertex_max, face_max, t_move_last, max_quality,
5704 cardinality] = r;
5705 if (max_quality < 0) {
5706 for (double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
5707 bundle = b_org_bundle;
5708 r = find_max_in_bundle_impl(edge, bundle, gamma);
5709 auto &[vertex_max, face_max, t_move_last, max_quality,
5710 cardinality] = r;
5711 MOFEM_LOG("EPSELF", Sev::warning)
5712 << "Back tracking: gamma " << gamma << " edge " << edge
5713 << " quality " << max_quality << " cardinality "
5714 << cardinality;
5715 if (max_quality > 0.01) {
5717 t_move_last(I) *= gamma;
5718 return r;
5719 }
5720 }
5722 t_move_last(I) = 0;
5723 }
5724 return r;
5725 };
5726
5727 // set tags with displacement of node and face energy
5728 auto set_tag_to_vertex_and_face = [&](auto &&r, auto &quality) {
5730 auto &[v, f, t_move, q, cardinality] = r;
5731
5732 if ((q > 0 && std::isnormal(q)) && energy > 0) {
5733
5734 MOFEM_LOG("EPSELF", Sev::inform)
5735 << "Set tag: vertex " << v << " face " << f << " "
5736 << max_edge << " move " << t_move << " energy " << energy
5737 << " quality " << q << " cardinality " << cardinality;
5738 CHKERR mField.get_moab().tag_set_data(th_position[0], &v, 1,
5739 &t_move(0));
5740 CHKERR mField.get_moab().tag_set_data(th_max_face_energy[0], &f,
5741 1, &energy);
5742 }
5743
5744 quality = q;
5746 };
5747
5748 double quality = -2;
5749 CHKERR set_tag_to_vertex_and_face(
5750
5751 find_max_in_bundle(max_edge, bundle),
5752
5753 quality
5754
5755 );
5756
5757 if (quality > 0 && std::isnormal(quality) && energy > 0) {
5758 MOFEM_LOG("EPSELF", Sev::inform)
5759 << "Crack face set with quality: " << quality;
5761 }
5762 }
5763
5764 if (!skip_negative)
5765 break;
5766 }
5767
5769 };
5770
5771 // map: {edge, {face, energy, optimal_angle}}
5772 MOFEM_LOG("EP", sev) << "Calculate orientation";
5773 std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
5774 edge_face_max_energy_map;
5775 CHKERR find_maximal_face_energy(front_edges, front_faces,
5776 edge_face_max_energy_map);
5777 CHKERR calculate_face_orientation(edge_face_max_energy_map);
5778
5779 MOFEM_LOG("EP", sev) << "Calculate node positions";
5780 CHKERR set_tag(
5781
5782 calculate_free_face_node_displacement(edge_face_max_energy_map),
5783 get_sort_by_energy(edge_face_max_energy_map)
5784
5785 );
5786
5788 };
5789
5790 MOFEM_LOG("EP", sev) << "Front edges " << frontEdges->size();
5791 CHKERR evaluate_face_energy_and_set_orientation(
5792 *frontEdges, get_adj_front(false), sides_pair, th_front_position);
5793 }
5794
5795 // exchange positions and energies from processor zero to all other
5796 CHKERR VecZeroEntries(vertexExchange.second);
5797 CHKERR VecGhostUpdateBegin(vertexExchange.second, INSERT_VALUES,
5798 SCATTER_FORWARD);
5799 CHKERR VecGhostUpdateEnd(vertexExchange.second, INSERT_VALUES,
5800 SCATTER_FORWARD);
5801 CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
5802 mField.get_moab(), vertexExchange, th_front_position[0]);
5803 CHKERR VecZeroEntries(faceExchange.second);
5804 CHKERR VecGhostUpdateBegin(faceExchange.second, INSERT_VALUES,
5805 SCATTER_FORWARD);
5806 CHKERR VecGhostUpdateEnd(faceExchange.second, INSERT_VALUES, SCATTER_FORWARD);
5807 CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
5808 mField.get_moab(), faceExchange, th_max_face_energy[0]);
5809
5810 auto get_max_moved_faces = [&]() {
5811 Range max_moved_faces;
5812 auto adj_front = get_adj_front(false);
5813 std::vector<double> face_energy(adj_front.size());
5814 CHKERR mField.get_moab().tag_get_data(th_max_face_energy[0], adj_front,
5815 face_energy.data());
5816 for (int i = 0; i != adj_front.size(); ++i) {
5817 if (face_energy[i] > std::numeric_limits<double>::epsilon()) {
5818 max_moved_faces.insert(adj_front[i]);
5819 }
5820 }
5821
5822 return boost::make_shared<Range>(max_moved_faces);
5823 };
5824
5825 // move all faces with energy larger than 0
5826 maxMovedFaces = get_max_moved_faces();
5827 MOFEM_LOG("EP", sev) << "Number of of moved faces: " << maxMovedFaces->size();
5828
5829#ifndef NDEBUG
5830 if (debug) {
5832 mField.get_moab(),
5833 "max_moved_faces_" +
5834 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
5835 *maxMovedFaces);
5836 }
5837#endif
5838
5840}
static auto get_two_sides_of_crack_surface(MoFEM::Interface &m_field, Range crack_faces)
constexpr double a
static const double eps
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
constexpr double a0
const double c
speed of light (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Index< 'm', 3 > m
CommInterface::EntitiesPetscVector vertexExchange
static enum EnergyReleaseSelector energyReleaseSelector
static auto exp(A &&t_w_vee, B &&theta)
Definition Lie.hpp:48
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition Tools.cpp:15

◆ createCrackSurfaceMeshset()

MoFEMErrorCode EshelbianCore::createCrackSurfaceMeshset ( )
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 6638 of file EshelbianPlasticity.cpp.

6638 {
6640 auto meshset_mng = mField.getInterface<MeshsetsManager>();
6641 while (meshset_mng->checkMeshset(addCrackMeshsetId, BLOCKSET))
6643 MOFEM_LOG("EP", Sev::inform)
6644 << "Crack added surface meshset " << addCrackMeshsetId;
6645 CHKERR meshset_mng->addMeshset(BLOCKSET, addCrackMeshsetId, "CRACK_COMPUTED");
6647};
MoFEMErrorCode addMeshset(const CubitBCType cubit_bc_type, const int ms_id, const std::string name="")
Add CUBIT meshset to manager.

◆ createExchangeVectors()

MoFEMErrorCode EshelbianCore::createExchangeVectors ( Sev  sev)
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 6943 of file EshelbianPlasticity.cpp.

6943 {
6945
6946 auto print_loc_size = [this](auto v, auto str, auto sev) {
6948 int size;
6949 CHKERR VecGetLocalSize(v.second, &size);
6950 int low, high;
6951 CHKERR VecGetOwnershipRange(v.second, &low, &high);
6952 MOFEM_LOG("EPSYNC", sev) << str << " local size " << size << " ( " << low
6953 << " " << high << " ) ";
6956 };
6957
6959 mField.get_comm(), mField.get_moab(), 3, 1, sev);
6960 CHKERR print_loc_size(volumeExchange, "volumeExchange", sev);
6962 mField.get_comm(), mField.get_moab(), 2, 1, Sev::inform);
6963 CHKERR print_loc_size(faceExchange, "faceExchange", sev);
6965 mField.get_comm(), mField.get_moab(), 1, 1, Sev::inform);
6966 CHKERR print_loc_size(edgeExchange, "edgeExchange", sev);
6968 mField.get_comm(), mField.get_moab(), 0, 3, Sev::inform);
6969 CHKERR print_loc_size(vertexExchange, "vertexExchange", sev);
6970
6972}
CommInterface::EntitiesPetscVector volumeExchange

◆ d_f_linear()

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

◆ d_f_log()

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

◆ d_f_log_e()

static double EshelbianCore::d_f_log_e ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 95 of file EshelbianCore.hpp.

95 {
96 if constexpr (use_quadratic_exp) {
97 return d_f_log_e_quadratic(v);
98 } else {
99 if (v > v_max)
100 return std::exp(v_max);
101 else
102 return std::exp(v);
103 }
104 }
static constexpr bool use_quadratic_exp
static constexpr double v_max
static double d_f_log_e_quadratic(const double v)

◆ d_f_log_e_quadratic()

static double EshelbianCore::d_f_log_e_quadratic ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 64 of file EshelbianCore.hpp.

64 {
65 if (v > v_max) {
66 double e = static_cast<double>(std::exp(v_max));
67 double dv = v - v_max;
68 return e * dv + e;
69 } else {
70 return static_cast<double>(std::exp(v));
71 }
72 }

◆ dd_f_linear()

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

◆ dd_f_log()

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

◆ dd_f_log_e()

static double EshelbianCore::dd_f_log_e ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 105 of file EshelbianCore.hpp.

105 {
106 if constexpr (use_quadratic_exp) {
107 return dd_f_log_e_quadratic(v);
108 } else {
109 if (v > v_max)
110 return 0.;
111 else
112 return std::exp(v);
113 }
114 }
static double dd_f_log_e_quadratic(const double v)

◆ dd_f_log_e_quadratic()

static double EshelbianCore::dd_f_log_e_quadratic ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 74 of file EshelbianCore.hpp.

74 {
75 if (v > v_max) {
76 return static_cast<double>(std::exp(v_max));
77 } else {
78 return static_cast<double>(std::exp(v));
79 }
80 }

◆ f_linear()

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

◆ f_log()

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

◆ f_log_e()

static double EshelbianCore::f_log_e ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 82 of file EshelbianCore.hpp.

82 {
83 if constexpr(use_quadratic_exp) {
84 return f_log_e_quadratic(v);
85 } else {
86 if (v > v_max)
87 // y = exp(v_max) * v + exp(v_max) * (1 - v_max);
88 // y/exp(v_max) = v + (1 - v_max);
89 // y/exp(v_max) - (1 - v_max) = v;
90 return std::exp(v_max) * v + std::exp(v_max) * (1 - v_max);
91 else
92 return std::exp(v);
93 }
94 }
static double f_log_e_quadratic(const double v)

◆ f_log_e_quadratic()

static double EshelbianCore::f_log_e_quadratic ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 54 of file EshelbianCore.hpp.

54 {
55 if (v > v_max) {
56 double e = static_cast<double>(std::exp(v_max));
57 double dv = v - v_max;
58 return 0.5 * e * dv * dv + e * dv + e;
59 } else {
60 return static_cast<double>(std::exp(v));
61 }
62 }

◆ 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
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 243 of file EshelbianCore.hpp.

244 {
246 for (auto it :
247 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
248
249 (boost::format("%s(.*)") % block_name).str()
250
251 ))
252
253 ) {
254 std::vector<double> block_attributes;
255 CHKERR it->getAttributes(block_attributes);
256 if (block_attributes.size() < nb_attributes) {
257 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
258 "In block %s expected %d attributes, but given %ld",
259 it->getName().c_str(), nb_attributes, block_attributes.size());
260 }
261 Range faces;
262 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
263 true);
264 bc_vec_ptr->emplace_back(it->getName(), block_attributes, faces);
265 }
267 }
IFACE getInterface() const
Get interface pointer to pointer of interface.

◆ getExternalStrain()

MoFEMErrorCode EshelbianCore::getExternalStrain ( )
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 6896 of file EshelbianPlasticity.cpp.

6896 {
6898
6899 auto getExternalStrain = [&](boost::shared_ptr<ExternalStrainVec> &ext_strain_vec_ptr,
6900 const std::string block_name,
6901 const int nb_attributes) {
6903 for (auto it : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
6904 std::regex((boost::format("%s(.*)") % block_name).str()))
6905 ) {
6906 std::vector<double> block_attributes;
6907 CHKERR it->getAttributes(block_attributes);
6908 if (block_attributes.size() < nb_attributes) {
6909 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
6910 "In block %s expected %d attributes, but given %ld",
6911 it->getName().c_str(), nb_attributes, block_attributes.size());
6912 }
6913
6914 auto get_block_ents = [&]() {
6915 Range ents;
6916 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, ents,
6917 true);
6918 return ents;
6919 };
6920 auto Ents = get_block_ents();
6921 ext_strain_vec_ptr->emplace_back(it->getName(), block_attributes,
6922 get_block_ents());
6923 }
6925 };
6926
6927 externalStrainVecPtr = boost::make_shared<ExternalStrainVec>();
6928 CHKERR getExternalStrain(externalStrainVecPtr, "EXTERNALSTRAIN", 2);
6929
6930 auto ts_pre_stretch =
6931 boost::make_shared<DynamicRelaxationTimeScale>("externalstrain_history.txt");
6932 for (auto &ext_strain_block: *externalStrainVecPtr) {
6933 MOFEM_LOG("EP", Sev::noisy)
6934 << "Add time scaling external strain: " << ext_strain_block.blockName;
6935 timeScaleMap[ext_strain_block.blockName] =
6937 ts_pre_stretch, "externalstrain_history", ".txt", ext_strain_block.blockName);
6938 }
6939
6941}
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
std::map< std::string, boost::shared_ptr< ScalingMethod > > timeScaleMap
MoFEMErrorCode getExternalStrain()
static boost::shared_ptr< ScalingMethod > get(boost::shared_ptr< ScalingMethod > ts, std::string file_prefix, std::string file_suffix, std::string block_name, Args &&...args)

◆ getOptions()

MoFEMErrorCode EshelbianCore::getOptions ( )
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 922 of file EshelbianPlasticity.cpp.

922 {
924 const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
925 const char *list_symm[] = {"symm", "not_symm"};
926 const char *list_release[] = {"griffith_force", "griffith_skeleton"};
927 const char *list_stretches[] = {"linear", "log", "log_quadratic"};
928 PetscInt choice_rot = EshelbianCore::rotSelector;
929 PetscInt choice_grad = EshelbianCore::gradApproximator;
930 PetscInt choice_symm = EshelbianCore::symmetrySelector;
931 PetscInt choice_release = EshelbianCore::energyReleaseSelector;
932 PetscInt choice_stretch = StretchSelector::LOG;
933 char analytical_expr_file_name[255] = "analytical_expr.py";
934
935 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
936 "none");
937 CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
938 spaceOrder, &spaceOrder, PETSC_NULLPTR);
939 CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
940 spaceH1Order, &spaceH1Order, PETSC_NULLPTR);
941 CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
942 "", materialOrder, &materialOrder, PETSC_NULLPTR);
943 CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
944 &alphaU, PETSC_NULLPTR);
945 CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
946 &alphaW, PETSC_NULLPTR);
947 CHKERR PetscOptionsScalar("-viscosity_alpha_omega", "rot viscosity", "",
948 alphaOmega, &alphaOmega, PETSC_NULLPTR);
949 CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
950 &alphaRho, PETSC_NULLPTR);
951 CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
952 LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
953 PETSC_NULLPTR);
954 CHKERR PetscOptionsEList("-grad", "gradient of defamation approximate", "",
955 list_rots, NO_H1_CONFIGURATION + 1,
956 list_rots[choice_grad], &choice_grad, PETSC_NULLPTR);
957 CHKERR PetscOptionsEList("-symm", "symmetric variant", "", list_symm, 2,
958 list_symm[choice_symm], &choice_symm, PETSC_NULLPTR);
959
960 CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
961 &EshelbianCore::exponentBase, PETSC_NULLPTR);
962 CHKERR PetscOptionsEList("-stretches", "stretches", "", list_stretches,
963 StretchSelector::STRETCH_SELECTOR_LAST,
964 list_stretches[choice_stretch], &choice_stretch,
965 PETSC_NULLPTR);
966
967 CHKERR PetscOptionsBool("-no_stretch", "do not solve for stretch", "",
968 noStretch, &noStretch, PETSC_NULLPTR);
969 CHKERR PetscOptionsBool("-set_singularity", "set singularity", "",
970 setSingularity, &setSingularity, PETSC_NULLPTR);
971 CHKERR PetscOptionsBool("-l2_user_base_scale", "streach scale", "",
972 l2UserBaseScale, &l2UserBaseScale, PETSC_NULLPTR);
973
974 // dynamic relaxation
975 CHKERR PetscOptionsBool("-dynamic_relaxation", "dynamic time relaxation", "",
976 dynamicRelaxation, &dynamicRelaxation, PETSC_NULLPTR);
977
978 // contact parameters
979 CHKERR PetscOptionsInt("-contact_max_post_proc_ref_level", "refinement level",
981 PETSC_NULLPTR);
982
983 // cracking parameters
984 CHKERR PetscOptionsBool("-cracking_on", "cracking ON", "", crackingOn,
985 &crackingOn, PETSC_NULLPTR);
986 CHKERR PetscOptionsScalar("-cracking_start_time", "cracking start time", "",
987 crackingStartTime, &crackingStartTime, PETSC_NULLPTR);
988 CHKERR PetscOptionsScalar("-griffith_energy", "Griffith energy", "",
989 griffithEnergy, &griffithEnergy, PETSC_NULLPTR);
990 CHKERR PetscOptionsEList("-energy_release_variant", "energy release variant",
991 "", list_release, 2, list_release[choice_release],
992 &choice_release, PETSC_NULLPTR);
993 CHKERR PetscOptionsInt("-nb_J_integral_levels", "Number of J integarl levels",
994 "", nbJIntegralLevels, &nbJIntegralLevels, PETSC_NULLPTR);
995
996 // internal stress
997 char tag_name[255] = "";
998 CHKERR PetscOptionsString("-internal_stress_tag_name",
999 "internal stress tag name", "", "", tag_name, 255,
1000 PETSC_NULLPTR);
1001 internalStressTagName = string(tag_name);
1002 CHKERR PetscOptionsInt(
1003 "-internal_stress_interp_order", "internal stress interpolation order",
1005 CHKERR PetscOptionsBool("-internal_stress_voigt", "Voigt index notation", "",
1007 PETSC_NULLPTR);
1008
1009 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-analytical_expr_file",
1010 analytical_expr_file_name, 255, PETSC_NULLPTR);
1011
1012 PetscOptionsEnd();
1013
1015 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
1016 "Unsupported internal stress interpolation order %d",
1018 }
1019
1020 if (setSingularity) {
1021 l2UserBaseScale = PETSC_TRUE;
1022 }
1023
1024 EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
1025 EshelbianCore::gradApproximator = static_cast<RotSelector>(choice_grad);
1026 EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
1027 EshelbianCore::symmetrySelector = static_cast<SymmetrySelector>(choice_symm);
1029 static_cast<EnergyReleaseSelector>(choice_release);
1030
1032 case StretchSelector::LINEAR:
1039 break;
1040 case StretchSelector::LOG:
1041 if (std::fabs(EshelbianCore::exponentBase - exp(1)) >
1042 std::numeric_limits<float>::epsilon()) {
1049 } else {
1056 }
1057 break;
1058 case StretchSelector::LOG_QUADRATIC:
1062 EshelbianCore::inv_f = [](const double x) {
1064 "No logarithmic quadratic stretch for this case");
1065 return 0;
1066 };
1069 break; // no stretch, do not use stretch functions
1070 default:
1071 SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
1072 break;
1073 };
1074
1075 MOFEM_LOG("EP", Sev::inform) << "spaceOrder: -space_order " << spaceOrder;
1076 MOFEM_LOG("EP", Sev::inform)
1077 << "spaceH1Order: -space_h1_order " << spaceH1Order;
1078 MOFEM_LOG("EP", Sev::inform)
1079 << "materialOrder: -material_order " << materialOrder;
1080 MOFEM_LOG("EP", Sev::inform) << "alphaU: -viscosity_alpha_u " << alphaU;
1081 MOFEM_LOG("EP", Sev::inform) << "alphaW: -viscosity_alpha_w " << alphaW;
1082 MOFEM_LOG("EP", Sev::inform)
1083 << "alphaOmega: -viscosity_alpha_omega " << alphaOmega;
1084 MOFEM_LOG("EP", Sev::inform) << "alphaRho: -density_alpha_rho " << alphaRho;
1085 MOFEM_LOG("EP", Sev::inform)
1086 << "Rotations: -rotations " << list_rots[EshelbianCore::rotSelector];
1087 MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
1088 << list_rots[EshelbianCore::gradApproximator];
1089 MOFEM_LOG("EP", Sev::inform)
1090 << "Symmetry: -symm " << list_symm[EshelbianCore::symmetrySelector];
1091 if (exponentBase != exp(1))
1092 MOFEM_LOG("EP", Sev::inform)
1093 << "Base exponent: -exponent_base " << EshelbianCore::exponentBase;
1094 else
1095 MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
1096 MOFEM_LOG("EP", Sev::inform)
1097 << "Stretch: -stretches " << list_stretches[choice_stretch];
1098 MOFEM_LOG("EP", Sev::inform) << "No stretch: -no_stretch " << (noStretch)
1099 ? "yes"
1100 : "no";
1101
1102 MOFEM_LOG("EP", Sev::inform)
1103 << "Dynamic relaxation: -dynamic_relaxation " << (dynamicRelaxation)
1104 ? "yes"
1105 : "no";
1106 MOFEM_LOG("EP", Sev::inform)
1107 << "Singularity: -set_singularity " << (setSingularity)
1108 ? "yes"
1109 : "no";
1110 MOFEM_LOG("EP", Sev::inform)
1111 << "L2 user base scale: -l2_user_base_scale " << (l2UserBaseScale)
1112 ? "yes"
1113 : "no";
1114
1115 MOFEM_LOG("EP", Sev::inform) << "Cracking on: -cracking_on " << (crackingOn)
1116 ? "yes"
1117 : "no";
1118 MOFEM_LOG("EP", Sev::inform)
1119 << "Cracking start time: -cracking_start_time " << crackingStartTime;
1120 MOFEM_LOG("EP", Sev::inform)
1121 << "Griffith energy: -griffith_energy " << griffithEnergy;
1122 MOFEM_LOG("EP", Sev::inform)
1123 << "Energy release variant: -energy_release_variant "
1124 << list_release[EshelbianCore::energyReleaseSelector];
1125 MOFEM_LOG("EP", Sev::inform)
1126 << "Number of J integral levels: -nb_J_integral_levels "
1128
1129#ifdef ENABLE_PYTHON_BINDING
1130 auto file_exists = [](std::string myfile) {
1131 std::ifstream file(myfile.c_str());
1132 if (file) {
1133 return true;
1134 }
1135 return false;
1136 };
1137
1138 if (file_exists(analytical_expr_file_name)) {
1139 MOFEM_LOG("EP", Sev::inform) << analytical_expr_file_name << " file found";
1140
1141 AnalyticalExprPythonPtr = boost::make_shared<AnalyticalExprPython>();
1142 CHKERR AnalyticalExprPythonPtr->analyticalExprInit(
1143 analytical_expr_file_name);
1144 AnalyticalExprPythonWeakPtr = AnalyticalExprPythonPtr;
1145 } else {
1146 MOFEM_LOG("EP", Sev::warning)
1147 << analytical_expr_file_name << " file NOT found";
1148 }
1149#endif
1150
1151 if (spaceH1Order == -1)
1153
1155}
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
static enum StretchSelector stretchSelector
static double inv_f_linear(const double v)
static double dd_f_log(const double v)
static boost::function< double(const double)> inv_dd_f
static double inv_d_f_log(const double v)
static PetscBool l2UserBaseScale
static int internalStressInterpOrder
static PetscBool crackingOn
static double dd_f_log_e(const double v)
static enum RotSelector rotSelector
static enum RotSelector gradApproximator
static PetscBool dynamicRelaxation
static double d_f_log(const double v)
static double crackingStartTime
static double inv_f_log(const double v)
static double dd_f_linear(const double v)
boost::shared_ptr< AnalyticalExprPython > AnalyticalExprPythonPtr
static double inv_d_f_linear(const double v)
static double inv_dd_f_log(const double v)
static std::string internalStressTagName
static enum SymmetrySelector symmetrySelector
static PetscBool internalStressVoigt
static double inv_dd_f_linear(const double v)
static PetscBool setSingularity
static double d_f_log_e(const double v)
static boost::function< double(const double)> dd_f
static double f_log_e(const double v)
static boost::function< double(const double)> d_f
static boost::function< double(const double)> inv_d_f
static double d_f_linear(const double v)
static double f_log(const double v)
static boost::function< double(const double)> inv_f
static double f_linear(const double v)

◆ getSpatialDispBc()

MoFEMErrorCode EshelbianCore::getSpatialDispBc ( )

[Getting norms]

Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 6718 of file EshelbianPlasticity.cpp.

6718 {
6720
6721 auto bc_mng = mField.getInterface<BcManager>();
6723 "", piolaStress, false, false);
6724
6725 bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
6726 for (auto bc : bc_mng->getBcMapByBlockName()) {
6727 if (auto disp_bc = bc.second->dispBcPtr) {
6728
6729 auto [field_name, block_name] =
6731 MOFEM_LOG("EP", Sev::inform)
6732 << "Field name: " << field_name << " Block name: " << block_name;
6733 MOFEM_LOG("EP", Sev::noisy) << "Displacement BC: " << *disp_bc;
6734
6735 std::vector<double> block_attributes(6, 0.);
6736 if (disp_bc->data.flag1 == 1) {
6737 block_attributes[0] = disp_bc->data.value1;
6738 block_attributes[3] = 1;
6739 }
6740 if (disp_bc->data.flag2 == 1) {
6741 block_attributes[1] = disp_bc->data.value2;
6742 block_attributes[4] = 1;
6743 }
6744 if (disp_bc->data.flag3 == 1) {
6745 block_attributes[2] = disp_bc->data.value3;
6746 block_attributes[5] = 1;
6747 }
6748 auto faces = bc.second->bcEnts.subset_by_dimension(2);
6749 bcSpatialDispVecPtr->emplace_back(block_name, block_attributes, faces);
6750 }
6751 }
6752 // old way of naming blocksets for displacement BCs
6753 CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
6754
6756 boost::make_shared<NormalDisplacementBcVec>();
6757 for (auto bc : bc_mng->getBcMapByBlockName()) {
6758 auto block_name = "(.*)NORMAL_DISPLACEMENT(.*)";
6759 std::regex reg_name(block_name);
6760 if (std::regex_match(bc.first, reg_name)) {
6761 auto [field_name, block_name] =
6763 MOFEM_LOG("EP", Sev::inform)
6764 << "Field name: " << field_name << " Block name: " << block_name;
6766 block_name, bc.second->bcAttributes,
6767 bc.second->bcEnts.subset_by_dimension(2));
6768 }
6769 }
6770
6772 boost::make_shared<AnalyticalDisplacementBcVec>();
6773
6774 for (auto bc : bc_mng->getBcMapByBlockName()) {
6775 auto block_name = "(.*)ANALYTICAL_DISPLACEMENT(.*)";
6776 std::regex reg_name(block_name);
6777 if (std::regex_match(bc.first, reg_name)) {
6778 auto [field_name, block_name] =
6780 MOFEM_LOG("EP", Sev::inform)
6781 << "Field name: " << field_name << " Block name: " << block_name;
6783 block_name, bc.second->bcAttributes,
6784 bc.second->bcEnts.subset_by_dimension(2));
6785 }
6786 }
6787
6788 auto ts_displacement =
6789 boost::make_shared<DynamicRelaxationTimeScale>("disp_history.txt");
6790 for (auto &bc : *bcSpatialDispVecPtr) {
6791 MOFEM_LOG("EP", Sev::noisy)
6792 << "Add time scaling displacement BC: " << bc.blockName;
6793 timeScaleMap[bc.blockName] =
6795 ts_displacement, "disp_history", ".txt", bc.blockName);
6796 }
6797
6798 auto ts_normal_displacement =
6799 boost::make_shared<DynamicRelaxationTimeScale>("normal_disp_history.txt");
6800 for (auto &bc : *bcSpatialNormalDisplacementVecPtr) {
6801 MOFEM_LOG("EP", Sev::noisy)
6802 << "Add time scaling normal displacement BC: " << bc.blockName;
6803 timeScaleMap[bc.blockName] =
6805 ts_normal_displacement, "normal_disp_history", ".txt",
6806 bc.blockName);
6807 }
6808
6810}
MoFEMErrorCode pushMarkDOFsOnEntities(const std::string problem_name, const std::string block_name, const std::string field_name, int lo, int hi, bool get_low_dim_ents=true)
Mark DOFs on block entities for boundary conditions.
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
Boundary condition manager for finite element problem setup.
static std::pair< std::string, std::string > extractStringFromBlockId(const std::string block_id, const std::string prb_name)
Extract block name and block name from block id.
Definition of the displacement bc data structure.
Definition BCData.hpp:72

◆ getSpatialRotationBc()

MoFEMErrorCode EshelbianCore::getSpatialRotationBc ( )
inline
Examples
ep.cpp, and mofem/users_modules/eshelbian_plasticity/ep.cpp.

Definition at line 271 of file EshelbianCore.hpp.

271 {
273 bcSpatialRotationVecPtr = boost::make_shared<BcRotVec>();
274 CHKERR getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_BC", 4);
275 CHKERR getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_AXIS_BC", 7);
276
277 auto ts_rotation =
278 boost::make_shared<DynamicRelaxationTimeScale>("rotation_history.txt");
279 for (auto &bc : *bcSpatialRotationVecPtr) {
280 timeScaleMap[bc.blockName] =
281 GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
282 ts_rotation, "rotation_history", ".txt", bc.blockName);
283 }
284
286 }

◆ getSpatialTractionBc()

MoFEMErrorCode EshelbianCore::getSpatialTractionBc ( )
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 6812 of file EshelbianPlasticity.cpp.

6812 {
6814
6815 auto bc_mng = mField.getInterface<BcManager>();
6817 false, false);
6818
6819 bcSpatialTractionVecPtr = boost::make_shared<TractionBcVec>();
6820
6821 for (auto bc : bc_mng->getBcMapByBlockName()) {
6822 if (auto force_bc = bc.second->forceBcPtr) {
6823
6824 auto [field_name, block_name] =
6826 MOFEM_LOG("EP", Sev::inform)
6827 << "Field name: " << field_name << " Block name: " << block_name;
6828 MOFEM_LOG("EP", Sev::noisy) << "Force BC: " << *force_bc;
6829
6830 std::vector<double> block_attributes(6, 0.);
6831 block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
6832 block_attributes[3] = 1;
6833 block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
6834 block_attributes[4] = 1;
6835 block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
6836 block_attributes[5] = 1;
6837 auto faces = bc.second->bcEnts.subset_by_dimension(2);
6838 bcSpatialTractionVecPtr->emplace_back(block_name, block_attributes,
6839 faces);
6840 }
6841 }
6842 CHKERR getBc(bcSpatialTractionVecPtr, "SPATIAL_TRACTION_BC", 6);
6843
6844 bcSpatialPressureVecPtr = boost::make_shared<PressureBcVec>();
6845 for (auto bc : bc_mng->getBcMapByBlockName()) {
6846 auto block_name = "(.*)PRESSURE(.*)";
6847 std::regex reg_name(block_name);
6848 if (std::regex_match(bc.first, reg_name)) {
6849
6850 auto [field_name, block_name] =
6852 MOFEM_LOG("EP", Sev::inform)
6853 << "Field name: " << field_name << " Block name: " << block_name;
6854 bcSpatialPressureVecPtr->emplace_back(
6855 block_name, bc.second->bcAttributes,
6856 bc.second->bcEnts.subset_by_dimension(2));
6857 }
6858 }
6859
6861 boost::make_shared<AnalyticalTractionBcVec>();
6862
6863 for (auto bc : bc_mng->getBcMapByBlockName()) {
6864 auto block_name = "(.*)ANALYTICAL_TRACTION(.*)";
6865 std::regex reg_name(block_name);
6866 if (std::regex_match(bc.first, reg_name)) {
6867 auto [field_name, block_name] =
6869 MOFEM_LOG("EP", Sev::inform)
6870 << "Field name: " << field_name << " Block name: " << block_name;
6872 block_name, bc.second->bcAttributes,
6873 bc.second->bcEnts.subset_by_dimension(2));
6874 }
6875 }
6876
6877 auto ts_traction =
6878 boost::make_shared<DynamicRelaxationTimeScale>("traction_history.txt");
6879 for (auto &bc : *bcSpatialTractionVecPtr) {
6880 timeScaleMap[bc.blockName] =
6882 ts_traction, "traction_history", ".txt", bc.blockName);
6883 }
6884
6885 auto ts_pressure =
6886 boost::make_shared<DynamicRelaxationTimeScale>("pressure_history.txt");
6887 for (auto &bc : *bcSpatialPressureVecPtr) {
6888 timeScaleMap[bc.blockName] =
6890 ts_pressure, "pressure_history", ".txt", bc.blockName);
6891 }
6892
6894}
Definition of the force bc data structure.
Definition BCData.hpp:135

◆ getSpatialTractionFreeBc()

MoFEMErrorCode EshelbianCore::getSpatialTractionFreeBc ( const EntityHandle  meshset = 0)
inline
Examples
ep.cpp, and mofem/users_modules/eshelbian_plasticity/ep.cpp.

Definition at line 305 of file EshelbianCore.hpp.

305 {
307 boost::shared_ptr<TractionFreeBc>(new TractionFreeBc());
308 return getTractionFreeBc(meshset, bcSpatialFreeTractionVecPtr, "CONTACT");
309 }
std::vector< Range > TractionFreeBc
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.

◆ gettingNorms()

MoFEMErrorCode EshelbianCore::gettingNorms ( )

[Getting norms]

Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 6650 of file EshelbianPlasticity.cpp.

6650 {
6652
6653 auto post_proc_norm_fe =
6654 boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
6655
6656 auto post_proc_norm_rule_hook = [](int, int, int p) -> int {
6657 return 2 * (p);
6658 };
6659 post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
6660
6661 post_proc_norm_fe->getUserPolynomialBase() =
6662 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
6663
6664 CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
6665 post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV}, materialH1Positions,
6667
6668 enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
6669 auto norms_vec =
6670 createVectorMPI(mField.get_comm(), LAST_NORM, PETSC_DETERMINE);
6671 CHKERR VecZeroEntries(norms_vec);
6672
6673 auto u_l2_ptr = boost::make_shared<MatrixDouble>();
6674 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
6675 post_proc_norm_fe->getOpPtrVector().push_back(
6677 post_proc_norm_fe->getOpPtrVector().push_back(
6679 post_proc_norm_fe->getOpPtrVector().push_back(
6680 new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_NORM_L2));
6681 post_proc_norm_fe->getOpPtrVector().push_back(
6682 new OpCalcNormL2Tensor1<SPACE_DIM>(u_h1_ptr, norms_vec, U_NORM_H1));
6683 post_proc_norm_fe->getOpPtrVector().push_back(
6684 new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_ERROR_L2,
6685 u_h1_ptr));
6686
6687 auto piola_ptr = boost::make_shared<MatrixDouble>();
6688 post_proc_norm_fe->getOpPtrVector().push_back(
6690 post_proc_norm_fe->getOpPtrVector().push_back(
6692 MBMAXTYPE));
6693
6694 post_proc_norm_fe->getOpPtrVector().push_back(
6695 new OpCalcNormL2Tensor2<3, 3>(piola_ptr, norms_vec, PIOLA_NORM));
6696
6697 TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
6699 *post_proc_norm_fe);
6700 TetPolynomialBase::switchCacheBaseOff<HDIV>({post_proc_norm_fe.get()});
6701
6702 CHKERR VecAssemblyBegin(norms_vec);
6703 CHKERR VecAssemblyEnd(norms_vec);
6704 const double *norms;
6705 CHKERR VecGetArrayRead(norms_vec, &norms);
6706 MOFEM_LOG("EP", Sev::inform) << "norm_u: " << std::sqrt(norms[U_NORM_L2]);
6707 MOFEM_LOG("EP", Sev::inform) << "norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
6708 MOFEM_LOG("EP", Sev::inform)
6709 << "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
6710 MOFEM_LOG("EP", Sev::inform)
6711 << "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
6712 CHKERR VecRestoreArrayRead(norms_vec, &norms);
6713
6715}
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
CGG User Polynomial Base.
Get norm of input MatrixDouble for Tensor1.
Get norm of input MatrixDouble for Tensor2.

◆ 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
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 2073 of file EshelbianPlasticity.cpp.

2075 {
2077
2078 // get skin from all tets
2079 Range tets;
2080 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
2081 Range tets_skin_part;
2082 Skinner skin(&mField.get_moab());
2083 CHKERR skin.find_skin(0, tets, false, tets_skin_part);
2084 ParallelComm *pcomm =
2085 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2086 Range tets_skin;
2087 CHKERR pcomm->filter_pstatus(tets_skin_part,
2088 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2089 PSTATUS_NOT, -1, &tets_skin);
2090
2091 bc_ptr->resize(3);
2092 for (int dd = 0; dd != 3; ++dd)
2093 (*bc_ptr)[dd] = tets_skin;
2094
2095 // Do not remove dofs on which traction is applied
2097 for (auto &v : *bcSpatialDispVecPtr) {
2098 if (v.flags[0])
2099 (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
2100 if (v.flags[1])
2101 (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
2102 if (v.flags[2])
2103 (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
2104 }
2105
2106 // Do not remove dofs on which rotation is applied
2108 for (auto &v : *bcSpatialRotationVecPtr) {
2109 (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
2110 (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
2111 (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
2112 }
2113
2115 for (auto &v : *bcSpatialNormalDisplacementVecPtr) {
2116 (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
2117 (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
2118 (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
2119 }
2120
2123 if (v.flags[0])
2124 (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
2125 if (v.flags[1])
2126 (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
2127 if (v.flags[2])
2128 (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
2129 }
2130
2132 for (auto &v : *bcSpatialTractionVecPtr) {
2133 (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
2134 (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
2135 (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
2136 }
2137
2139 for (auto &v : *bcSpatialAnalyticalTractionVecPtr) {
2140 (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
2141 (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
2142 (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
2143 }
2144
2146 for (auto &v : *bcSpatialPressureVecPtr) {
2147 (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
2148 (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
2149 (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
2150 }
2151
2152 // remove contact
2153 for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
2154 std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
2155 Range faces;
2156 CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
2157 true);
2158 (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
2159 (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
2160 (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
2161 }
2162
2164}
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33

◆ inv_d_f_linear()

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

◆ inv_d_f_log()

static double EshelbianCore::inv_d_f_log ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 151 of file EshelbianCore.hpp.

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

◆ inv_d_f_log_e()

static double EshelbianCore::inv_d_f_log_e ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 124 of file EshelbianCore.hpp.

124 {
125 if (v > std::exp(v_min))
126 return 1. / v;
127 else
128 return 1. / exp(v_min);
129 }
static constexpr double v_min

◆ inv_dd_f_linear()

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

◆ inv_dd_f_log()

static double EshelbianCore::inv_dd_f_log ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 154 of file EshelbianCore.hpp.

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

◆ inv_dd_f_log_e()

static double EshelbianCore::inv_dd_f_log_e ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 130 of file EshelbianCore.hpp.

130 {
131 if (v > std::exp(v_min))
132 return -1. / (v * v);
133 else
134 return 0.;
135 }

◆ inv_f_linear()

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

◆ inv_f_log()

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

◆ inv_f_log_e()

static double EshelbianCore::inv_f_log_e ( const double  v)
inlinestatic
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 115 of file EshelbianCore.hpp.

115 {
116 if (v > std::exp(v_min))
117 return std::log(v);
118 else
119 // y = exp(v_min) * v + exp(v_min) * (1 - v_min);
120 // y/exp(v_min) = v + (1 - v_min);
121 // y/exp(v_min) - (1 - v_min) = v;
122 return v / exp(v_min) - (1 - v_min);
123 }

◆ 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
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 3592 of file EshelbianPlasticity.cpp.

3594 {
3596
3597 // mark crack surface
3598 if (crackingOn) {
3599 Tag th;
3600 rval = mField.get_moab().tag_get_handle("CRACK", th);
3601 if (rval == MB_SUCCESS) {
3602 CHKERR mField.get_moab().tag_delete(th);
3603 }
3604 int def_val[] = {0};
3605 CHKERR mField.get_moab().tag_get_handle(
3606 "CRACK", 1, MB_TYPE_INTEGER, th, MB_TAG_SPARSE | MB_TAG_CREAT, def_val);
3607 int mark[] = {1};
3608 CHKERR mField.get_moab().tag_clear_data(th, *crackFaces, mark);
3609 tags_to_transfer.push_back(th);
3610
3611 auto get_tag = [&](auto name, auto dim) {
3612 auto &mob = mField.get_moab();
3613 Tag tag;
3614 double def_val[] = {0., 0., 0.};
3615 CHK_MOAB_THROW(mob.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3616 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3617 "create tag");
3618 return tag;
3619 };
3620
3621 tags_to_transfer.push_back(get_tag("MaterialForce", 3));
3622 // tags_to_transfer.push_back(get_tag("GriffithForce", 1));
3623 }
3624
3625 // add tags to transfer
3626 for (auto t : listTagsToTransfer) {
3627 tags_to_transfer.push_back(t);
3628 }
3629
3630 if (!dataAtPts) {
3631 dataAtPts =
3632 boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
3633 }
3634
3635 struct exclude_sdf {
3636 exclude_sdf(Range &&r) : map(r) {}
3637 bool operator()(FEMethod *fe_method_ptr) {
3638 auto ent = fe_method_ptr->getFEEntityHandle();
3639 if (map.find(ent) != map.end()) {
3640 return false;
3641 }
3642 return true;
3643 }
3644
3645 private:
3646 Range map;
3647 };
3648
3649 contactTreeRhs->exeTestHook =
3650 exclude_sdf(get_range_from_block(mField, "CONTACT_SDF", SPACE_DIM - 1));
3651
3653
3654 auto get_post_proc = [&](auto &post_proc_mesh, auto sense) {
3656 auto post_proc_ptr =
3657 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
3658 mField, post_proc_mesh);
3659 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
3660 post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
3662
3663 auto domain_ops = [&](auto &fe, int sense) {
3665 fe.getUserPolynomialBase() =
3666 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
3667 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3668 fe.getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions,
3670 auto piola_scale_ptr = boost::make_shared<double>(1.0);
3671 fe.getOpPtrVector().push_back(physicalEquations->returnOpSetScale(
3672 piola_scale_ptr, physicalEquations));
3673 fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
3674 piolaStress, dataAtPts->getApproxPAtPts(), piola_scale_ptr));
3675 fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
3676 bubbleField, dataAtPts->getApproxPAtPts(), piola_scale_ptr,
3677 SmartPetscObj<Vec>(), MBMAXTYPE));
3678 if (noStretch) {
3679 fe.getOpPtrVector().push_back(
3680 physicalEquations->returnOpCalculateStretchFromStress(
3682 } else {
3683 fe.getOpPtrVector().push_back(
3685 stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
3686 }
3687 if (var_vector) {
3688 auto vec = SmartPetscObj<Vec>(var_vector, true);
3689 fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
3690 piolaStress, dataAtPts->getVarPiolaPts(),
3691 boost::make_shared<double>(1), vec));
3692 fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
3693 bubbleField, dataAtPts->getVarPiolaPts(),
3694 boost::make_shared<double>(1), vec, MBMAXTYPE));
3695 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3696 rotAxis, dataAtPts->getVarRotAxisPts(), vec, MBTET));
3697 if (noStretch) {
3698 fe.getOpPtrVector().push_back(
3699 physicalEquations->returnOpCalculateVarStretchFromStress(
3701 } else {
3702 fe.getOpPtrVector().push_back(
3704 stretchTensor, dataAtPts->getVarLogStreachPts(), vec, MBTET));
3705 }
3706 }
3707
3708 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3709 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
3710 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3711 rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
3712
3713 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3714 spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
3715 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3716 spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
3717 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
3718 spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
3719 // evaluate derived quantities
3720 fe.getOpPtrVector().push_back(
3722
3723 // evaluate integration points
3724 fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
3725 tag, true, false, dataAtPts, physicalEquations));
3726 if (auto op =
3727 physicalEquations->returnOpCalculateEnergy(dataAtPts, nullptr)) {
3728 fe.getOpPtrVector().push_back(op);
3729 fe.getOpPtrVector().push_back(new OpCalculateEshelbyStress(dataAtPts));
3730 }
3731
3732 // // post-proc
3735 VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator>;
3736
3737 struct OpSidePPMap : public OpPPMap {
3738 OpSidePPMap(moab::Interface &post_proc_mesh,
3739 std::vector<EntityHandle> &map_gauss_pts,
3740 DataMapVec data_map_scalar, DataMapMat data_map_vec,
3741 DataMapMat data_map_mat, DataMapMat data_symm_map_mat,
3742 int sense)
3743 : OpPPMap(post_proc_mesh, map_gauss_pts, data_map_scalar,
3744 data_map_vec, data_map_mat, data_symm_map_mat),
3745 tagSense(sense) {}
3746
3747 MoFEMErrorCode doWork(int side, EntityType type,
3750
3751 if (tagSense != OpPPMap::getSkeletonSense())
3753
3754 CHKERR OpPPMap::doWork(side, type, data);
3756 }
3757
3758 private:
3759 int tagSense;
3760 };
3761
3762 OpPPMap::DataMapMat vec_fields;
3763 vec_fields["SpatialDisplacementL2"] = dataAtPts->getSmallWL2AtPts();
3764 vec_fields["SpatialDisplacementH1"] = dataAtPts->getSmallWH1AtPts();
3765 vec_fields["Omega"] = dataAtPts->getRotAxisAtPts();
3766 vec_fields["ContactDisplacement"] = dataAtPts->getContactL2AtPts();
3767 vec_fields["AngularMomentum"] = dataAtPts->getLeviKirchhoffAtPts();
3768 vec_fields["X"] = dataAtPts->getLargeXH1AtPts();
3769 if (!noStretch) {
3770 vec_fields["EiegnLogStreach"] = dataAtPts->getEigenValsAtPts();
3771 }
3772 if (var_vector) {
3773 auto vec = SmartPetscObj<Vec>(var_vector, true);
3774 vec_fields["VarOmega"] = dataAtPts->getVarRotAxisPts();
3775 vec_fields["VarSpatialDisplacementL2"] =
3776 boost::make_shared<MatrixDouble>();
3777 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3778 spatialL2Disp, vec_fields["VarSpatialDisplacementL2"], vec, MBTET));
3779 }
3780 if (f_residual) {
3781 auto vec = SmartPetscObj<Vec>(f_residual, true);
3782 vec_fields["ResSpatialDisplacementL2"] =
3783 boost::make_shared<MatrixDouble>();
3784 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3785 spatialL2Disp, vec_fields["ResSpatialDisplacementL2"], vec, MBTET));
3786 vec_fields["ResOmega"] = boost::make_shared<MatrixDouble>();
3787 fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3788 rotAxis, vec_fields["ResOmega"], vec, MBTET));
3789 }
3790
3791 OpPPMap::DataMapMat mat_fields;
3792 mat_fields["PiolaStress"] = dataAtPts->getApproxPAtPts();
3793 if (var_vector) {
3794 mat_fields["VarPiolaStress"] = dataAtPts->getVarPiolaPts();
3795 }
3796 if (f_residual) {
3797 auto vec = SmartPetscObj<Vec>(f_residual, true);
3798 mat_fields["ResPiolaStress"] = boost::make_shared<MatrixDouble>();
3799 fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
3800 piolaStress, mat_fields["ResPiolaStress"],
3801 boost::make_shared<double>(1), vec));
3802 fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
3803 bubbleField, mat_fields["ResPiolaStress"],
3804 boost::make_shared<double>(1), vec, MBMAXTYPE));
3805 }
3806 if (!internalStressTagName.empty()) {
3807 mat_fields[internalStressTagName] = dataAtPts->getInternalStressAtPts();
3808 switch (internalStressInterpOrder) {
3809 case 0:
3810 fe.getOpPtrVector().push_back(
3812 break;
3813 case 1:
3814 fe.getOpPtrVector().push_back(
3816 break;
3817 default:
3818 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
3819 "Unsupported internal stress interpolation order %d",
3821 }
3822 }
3823
3824 OpPPMap::DataMapMat mat_fields_symm;
3825 mat_fields_symm["LogSpatialStretch"] =
3826 dataAtPts->getLogStretchTensorAtPts();
3827 mat_fields_symm["SpatialStretch"] = dataAtPts->getStretchTensorAtPts();
3828 if (var_vector) {
3829 mat_fields_symm["VarLogSpatialStretch"] =
3830 dataAtPts->getVarLogStreachPts();
3831 }
3832 if (f_residual) {
3833 auto vec = SmartPetscObj<Vec>(f_residual, true);
3834 if (!noStretch) {
3835 mat_fields_symm["ResLogSpatialStretch"] =
3836 boost::make_shared<MatrixDouble>();
3837 fe.getOpPtrVector().push_back(
3839 stretchTensor, mat_fields_symm["ResLogSpatialStretch"], vec,
3840 MBTET));
3841 }
3842 }
3843
3844 fe.getOpPtrVector().push_back(
3845
3846 new OpSidePPMap(
3847
3848 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3849
3850 {},
3851
3852 vec_fields,
3853
3854 mat_fields,
3855
3856 mat_fields_symm,
3857
3858 sense
3859
3860 )
3861
3862 );
3863
3864 fe.getOpPtrVector().push_back(new OpPostProcDataStructure(
3865 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3866 dataAtPts, sense));
3867
3869 };
3870
3871 post_proc_ptr->getOpPtrVector().push_back(
3873 dataAtPts->getContactL2AtPts()));
3874
3875 auto X_h1_ptr = boost::make_shared<MatrixDouble>();
3876 // H1 material positions
3877 post_proc_ptr->getOpPtrVector().push_back(
3879 dataAtPts->getLargeXH1AtPts()));
3880
3881 // domain
3884 domain_ops(*(op_loop_side->getSideFEPtr()), sense);
3885 post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
3886
3887 return post_proc_ptr;
3888 };
3889
3890 // contact
3891 auto calcs_side_traction_and_displacements = [&](auto &post_proc_ptr,
3892 auto &pip) {
3894 auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
3895 // evaluate traction
3896 using EleOnSide =
3898 using SideEleOp = EleOnSide::UserDataOperator;
3899 auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
3900 mField, elementVolumeName, SPACE_DIM, Sev::noisy);
3901 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
3902 boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
3903 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3904 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
3906 op_loop_domain_side->getOpPtrVector().push_back(
3908 piolaStress, contact_common_data_ptr->contactTractionPtr(),
3909 boost::make_shared<double>(1.0)));
3910 pip.push_back(op_loop_domain_side);
3911 // evaluate contact displacement and contact conditions
3912 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
3913 pip.push_back(new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
3915 contactDisp, contact_common_data_ptr->contactDispPtr()));
3916 pip.push_back(new OpTreeSearch(
3917 contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
3918 get_range_from_block(mField, "CONTACT", SPACE_DIM - 1),
3919 &post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
3920
3921 if (f_residual) {
3922
3923 using BoundaryEle =
3925 auto op_this = new OpLoopThis<BoundaryEle>(mField, contactElement);
3926 pip.push_back(op_this);
3927 auto contact_residual = boost::make_shared<MatrixDouble>();
3928 op_this->getOpPtrVector().push_back(
3930 contactDisp, contact_residual,
3931 SmartPetscObj<Vec>(f_residual, true)));
3933 op_this->getOpPtrVector().push_back(
3934
3935 new OpPPMap(
3936
3937 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3938
3939 {},
3940
3941 {{"res_contact", contact_residual}},
3942
3943 {},
3944
3945 {}
3946
3947 )
3948
3949 );
3950 }
3951
3953 };
3954
3955 auto post_proc_mesh = boost::make_shared<moab::Core>();
3956 auto post_proc_ptr = get_post_proc(post_proc_mesh, 1);
3957 auto post_proc_negative_sense_ptr = get_post_proc(post_proc_mesh, -1);
3958 CHKERR calcs_side_traction_and_displacements(post_proc_ptr,
3959 post_proc_ptr->getOpPtrVector());
3960
3961 auto own_tets =
3963 .subset_by_dimension(SPACE_DIM);
3964 Range own_faces;
3965 CHKERR mField.get_moab().get_adjacencies(own_tets, SPACE_DIM - 1, true,
3966 own_faces, moab::Interface::UNION);
3967
3968 auto get_post_negative = [&](auto &&ents) {
3969 auto crack_faces_pos = ents;
3970 auto crack_faces_neg = crack_faces_pos;
3971 auto skin = get_skin(mField, own_tets);
3972 auto crack_on_proc_skin = intersect(crack_faces_pos, skin);
3973 for (auto f : crack_on_proc_skin) {
3974 Range tet;
3975 CHKERR mField.get_moab().get_adjacencies(&f, 1, SPACE_DIM, false, tet);
3976 tet = intersect(tet, own_tets);
3977 int side_number, sense, offset;
3978 CHKERR mField.get_moab().side_number(tet[0], f, side_number, sense,
3979 offset);
3980 if (sense == 1) {
3981 crack_faces_neg.erase(f);
3982 } else {
3983 crack_faces_pos.erase(f);
3984 }
3985 }
3986 return std::make_pair(crack_faces_pos, crack_faces_neg);
3987 };
3988
3989 auto get_crack_faces = [&](auto crack_faces) {
3990 auto get_adj = [&](auto e, auto dim) {
3991 Range adj;
3992 CHKERR mField.get_moab().get_adjacencies(e, dim, true, adj,
3993 moab::Interface::UNION);
3994 return adj;
3995 };
3996 auto tets = get_adj(crack_faces, 3);
3997 auto faces = subtract(get_adj(tets, 2), crack_faces);
3998 tets = subtract(tets, get_adj(faces, 3));
3999 return subtract(crack_faces, get_adj(tets, 2));
4000 };
4001
4002 auto [crack_faces_pos, crack_faces_neg] =
4003 get_post_negative(intersect(own_faces, get_crack_faces(*crackFaces)));
4004
4005 auto only_crack_faces = [&](FEMethod *fe_method_ptr) {
4006 auto ent = fe_method_ptr->getFEEntityHandle();
4007 if (crack_faces_pos.find(ent) == crack_faces_pos.end()) {
4008 return false;
4009 }
4010 return true;
4011 };
4012
4013 auto only_negative_crack_faces = [&](FEMethod *fe_method_ptr) {
4014 auto ent = fe_method_ptr->getFEEntityHandle();
4015 if (crack_faces_neg.find(ent) == crack_faces_neg.end()) {
4016 return false;
4017 }
4018 return true;
4019 };
4020
4021 auto get_adj_front = [&]() {
4022 auto skeleton_faces = *skeletonFaces;
4023 Range adj_front;
4024 CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2, true, adj_front,
4025 moab::Interface::UNION);
4026
4027 adj_front = intersect(adj_front, skeleton_faces);
4028 adj_front = subtract(adj_front, *crackFaces);
4029 adj_front = intersect(own_faces, adj_front);
4030 return adj_front;
4031 };
4032
4033 post_proc_ptr->setTagsToTransfer(tags_to_transfer);
4034 post_proc_negative_sense_ptr->setTagsToTransfer(tags_to_transfer);
4035
4036 auto post_proc_begin =
4038 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
4040 post_proc_ptr->exeTestHook = only_crack_faces;
4041 post_proc_negative_sense_ptr->exeTestHook = only_negative_crack_faces;
4043 dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
4045 post_proc_negative_sense_ptr, 0,
4047
4048 constexpr bool debug = false;
4049 if (debug) {
4050
4051 auto [adj_front_pos, adj_front_neg] =
4052 get_post_negative(filter_owners(mField, get_adj_front()));
4053
4054 auto only_front_faces_pos = [adj_front_pos](FEMethod *fe_method_ptr) {
4055 auto ent = fe_method_ptr->getFEEntityHandle();
4056 if (adj_front_pos.find(ent) == adj_front_pos.end()) {
4057 return false;
4058 }
4059 return true;
4060 };
4061
4062 auto only_front_faces_neg = [adj_front_neg](FEMethod *fe_method_ptr) {
4063 auto ent = fe_method_ptr->getFEEntityHandle();
4064 if (adj_front_neg.find(ent) == adj_front_neg.end()) {
4065 return false;
4066 }
4067 return true;
4068 };
4069
4070 post_proc_ptr->exeTestHook = only_front_faces_pos;
4072 dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
4073 post_proc_negative_sense_ptr->exeTestHook = only_front_faces_neg;
4075 post_proc_negative_sense_ptr, 0,
4077 }
4078 auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
4079 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
4080
4081 CHKERR post_proc_end.writeFile(file.c_str());
4083}
static auto filter_owners(MoFEM::Interface &m_field, Range skin)
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:546
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:536
PostProcBrokenMeshInMoabBaseEndImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseEnd
Enable to run stack of post-processing elements. Use this to end stack.
PostProcBrokenMeshInMoabBaseBeginImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseBegin
Enable to run stack of post-processing elements. Use this to begin stack.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
std::vector< Tag > listTagsToTransfer
list of tags to transfer to postprocessor
boost::shared_ptr< ContactTree > contactTreeRhs
Make a contact tree.
static Range getPartEntities(moab::Interface &moab, int part)
Data on single entity (This is passed as argument to DataOperator::doWork)
Structure for user loop methods on finite elements.
EntityHandle getFEEntityHandle() const
Get the entity handle of the current finite element.
Calculate trace of vector (Hdiv/Hcurl) space.
Execute "this" element in the operator.
Post post-proc data at points from hash maps.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Template struct for dimension-specific finite element types.
intrusive_ptr for managing petsc objects

◆ postProcessSkeletonResults()

MoFEMErrorCode EshelbianCore::postProcessSkeletonResults ( const int  tag,
const std::string  file,
Vec  f_residual = PETSC_NULLPTR,
std::vector< Tag tags_to_transfer = {} 
)
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 4086 of file EshelbianPlasticity.cpp.

4088 {
4090
4092
4093 auto post_proc_mesh = boost::make_shared<moab::Core>();
4094 auto post_proc_ptr =
4095 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
4096 mField, post_proc_mesh);
4097 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM - 1, SPACE_DIM>::add(
4098 post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
4100
4101 auto hybrid_disp = boost::make_shared<MatrixDouble>();
4102 post_proc_ptr->getOpPtrVector().push_back(
4104 post_proc_ptr->getOpPtrVector().push_back(
4106 hybridSpatialDisp, dataAtPts->getGradHybridDispAtPts()));
4107
4108 auto op_loop_domain_side =
4110 mField, elementVolumeName, SPACE_DIM, Sev::noisy);
4111 post_proc_ptr->getOpPtrVector().push_back(op_loop_domain_side);
4112
4113 // evaluated in side domain, that is op_loop_domain_side
4114 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
4115 boost::make_shared<CGGUserPolynomialBase>();
4116 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
4117 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
4119 op_loop_domain_side->getOpPtrVector().push_back(
4121 piolaStress, dataAtPts->getApproxPAtPts()));
4122 op_loop_domain_side->getOpPtrVector().push_back(
4124 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
4125 op_loop_domain_side->getOpPtrVector().push_back(
4127 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
4128 if (noStretch) {
4129 op_loop_domain_side->getOpPtrVector().push_back(
4130 physicalEquations->returnOpCalculateStretchFromStress(
4132 } else {
4133 op_loop_domain_side->getOpPtrVector().push_back(
4135 stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
4136 }
4137
4139
4140 OpPPMap::DataMapMat vec_fields;
4141 vec_fields["HybridDisplacement"] = hybrid_disp;
4142 vec_fields["Omega"] = dataAtPts->getRotAxisAtPts();
4143 OpPPMap::DataMapMat mat_fields;
4144 mat_fields["PiolaStress"] = dataAtPts->getApproxPAtPts();
4145 mat_fields["HybridDisplacementGradient"] =
4146 dataAtPts->getGradHybridDispAtPts();
4147 OpPPMap::DataMapMat mat_fields_symm;
4148 mat_fields_symm["LogSpatialStretch"] = dataAtPts->getLogStretchTensorAtPts();
4149
4150 post_proc_ptr->getOpPtrVector().push_back(
4151
4152 new OpPPMap(
4153
4154 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4155
4156 {},
4157
4158 vec_fields,
4159
4160 mat_fields,
4161
4162 mat_fields_symm
4163
4164 )
4165
4166 );
4167
4168 if (f_residual) {
4169 auto hybrid_res = boost::make_shared<MatrixDouble>();
4170 post_proc_ptr->getOpPtrVector().push_back(
4172 hybridSpatialDisp, hybrid_res,
4173 SmartPetscObj<Vec>(f_residual, true)));
4175 post_proc_ptr->getOpPtrVector().push_back(
4176
4177 new OpPPMap(
4178
4179 post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
4180
4181 {},
4182
4183 {{"res_hybrid", hybrid_res}},
4184
4185 {},
4186
4187 {}
4188
4189 )
4190
4191 );
4192 }
4193
4194 post_proc_ptr->setTagsToTransfer(tags_to_transfer);
4195
4196 auto post_proc_begin =
4198 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
4200 auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
4201 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
4202
4203 CHKERR post_proc_end.writeFile(file.c_str());
4204
4206}

◆ projectGeometry()

MoFEMErrorCode EshelbianCore::projectGeometry ( const EntityHandle  meshset = 0)
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 1406 of file EshelbianPlasticity.cpp.

1406 {
1408
1409 Range meshset_ents;
1410 CHKERR mField.get_moab().get_entities_by_handle(meshset, meshset_ents);
1411
1412 auto project_ho_geometry = [&](auto field) {
1414 return mField.loop_dofs(field, ent_method);
1415 };
1416 CHKERR project_ho_geometry(materialH1Positions);
1417
1418 auto get_adj_front_edges = [&](auto &front_edges) {
1419 Range front_crack_nodes;
1420 Range crack_front_edges_with_both_nodes_not_at_front;
1421
1422 if (mField.get_comm_rank() == 0) {
1423 auto &moab = mField.get_moab();
1424 CHKERR moab.get_connectivity(front_edges, front_crack_nodes, true);
1425 Range crack_front_edges;
1426 CHKERR moab.get_adjacencies(front_crack_nodes, SPACE_DIM - 2, false,
1427 crack_front_edges, moab::Interface::UNION);
1428 crack_front_edges_with_both_nodes_not_at_front =
1429 subtract(crack_front_edges, front_edges);
1430 }
1431
1432 front_crack_nodes = send_type(mField, front_crack_nodes, MBVERTEX);
1433 crack_front_edges_with_both_nodes_not_at_front = send_type(
1434 mField, crack_front_edges_with_both_nodes_not_at_front, MBEDGE);
1435
1436 return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
1437 boost::make_shared<Range>(
1438 crack_front_edges_with_both_nodes_not_at_front));
1439 };
1440
1441 crackFaces = boost::make_shared<Range>(
1442 get_range_from_block(mField, "CRACK", SPACE_DIM - 1));
1443 frontEdges =
1444 boost::make_shared<Range>(get_crack_front_edges(mField, *crackFaces));
1445 auto [front_vertices, front_adj_edges] = get_adj_front_edges(*frontEdges);
1446 frontVertices = front_vertices;
1447 frontAdjEdges = front_adj_edges;
1448
1449#ifndef NDEBUG
1450 if (crackingOn) {
1451 auto rank = mField.get_comm_rank();
1452 // CHKERR save_range(mField.get_moab(),
1453 // (boost::format("meshset_ents_%d.vtk") % rank).str(),
1454 // meshset_ents);
1456 (boost::format("crack_faces_%d.vtk") % rank).str(),
1457 *crackFaces);
1459 (boost::format("front_edges_%d.vtk") % rank).str(),
1460 *frontEdges);
1461 // CHKERR save_range(mField.get_moab(),
1462 // (boost::format("front_vertices_%d.vtk") % rank).str(),
1463 // *frontVertices);
1464 // CHKERR save_range(mField.get_moab(),
1465 // (boost::format("front_adj_edges_%d.vtk") % rank).str(),
1466 // *frontAdjEdges);
1467 }
1468#endif // NDEBUG
1469
1470 auto set_singular_dofs = [&](auto &front_adj_edges, auto &front_vertices) {
1472 auto &moab = mField.get_moab();
1473
1474 double eps = 1;
1475 double beta = 0;
1476 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "-singularity_eps", &beta,
1477 PETSC_NULLPTR);
1478 MOFEM_LOG("EP", Sev::inform) << "Singularity eps " << beta;
1479 eps -= beta;
1480
1481 for (auto edge : front_adj_edges) {
1482 int num_nodes;
1483 const EntityHandle *conn;
1484 CHKERR moab.get_connectivity(edge, conn, num_nodes, false);
1485 double coords[6];
1486 CHKERR moab.get_coords(conn, num_nodes, coords);
1487 const double dir[3] = {coords[3] - coords[0], coords[4] - coords[1],
1488 coords[5] - coords[2]};
1489 double dof[3] = {0, 0, 0};
1490 if (front_vertices.find(conn[0]) != front_vertices.end()) {
1491 for (int dd = 0; dd != 3; dd++) {
1492 dof[dd] = -dir[dd] * eps;
1493 }
1494 } else if (front_vertices.find(conn[1]) != front_vertices.end()) {
1495 for (int dd = 0; dd != 3; dd++) {
1496 dof[dd] = +dir[dd] * eps;
1497 }
1498 }
1500 mField, materialH1Positions, edge, dit)) {
1501 const int idx = dit->get()->getEntDofIdx();
1502 if (idx > 2) {
1503 dit->get()->getFieldData() = 0;
1504 } else {
1505 dit->get()->getFieldData() = dof[idx];
1506 }
1507 }
1508 }
1509
1511 };
1512
1513 if (setSingularity)
1514 CHKERR set_singular_dofs(*frontAdjEdges, *frontVertices);
1515
1517}
static auto get_crack_front_edges(MoFEM::Interface &m_field, Range crack_faces)
#define _IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(MFIELD, NAME, ENT, IT)
loop over all dofs from a moFEM field and particular field
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Projection of edge entities with one mid-node on hierarchical basis.

◆ query_interface()

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

Getting interface of core database.

Parameters
uuidunique ID of interface
ifacereturned pointer to interface
Returns
error code
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 898 of file EshelbianPlasticity.cpp.

899 {
900 *iface = const_cast<EshelbianCore *>(this);
901 return 0;
902}

◆ saveOrgCoords()

MoFEMErrorCode EshelbianCore::saveOrgCoords ( )
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 6618 of file EshelbianPlasticity.cpp.

6618 {
6620 auto crack_faces =
6621 get_range_from_block(mField, "CRACK_COMPUTED", SPACE_DIM - 1);
6622 Range conn;
6623 CHKERR mField.get_moab().get_connectivity(crack_faces, conn, true);
6624 Range verts;
6625 CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, verts);
6626 verts = subtract(verts, conn);
6627 std::vector<double> coords(3 * verts.size());
6628 CHKERR mField.get_moab().get_coords(verts, coords.data());
6629 double def_coords[] = {0., 0., 0.};
6630 Tag th_org_coords;
6631 CHKERR mField.get_moab().tag_get_handle(
6632 "ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
6633 MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
6634 CHKERR mField.get_moab().tag_set_data(th_org_coords, verts, coords.data());
6636}

◆ 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
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 2189 of file EshelbianPlasticity.cpp.

2192 {
2194
2195 auto bubble_cache =
2196 boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0, MatrixDouble());
2197 fe->getUserPolynomialBase() =
2198 boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
2199 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2200 fe->getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions, frontAdjEdges);
2201
2202 // set integration rule
2203 fe->getRuleHook = [](int, int, int) { return -1; };
2204 fe->setRuleHook = SetIntegrationAtFrontVolume(frontVertices, frontAdjEdges);
2205 // fe->getRuleHook = VolRule();
2206
2207 if (!dataAtPts) {
2208 dataAtPts =
2209 boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
2210 dataAtPts->physicsPtr = physicalEquations;
2211 }
2212
2213 // calculate fields values
2214 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
2215 piolaStress, dataAtPts->getApproxPAtPts()));
2216 fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
2217 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
2218 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
2219 piolaStress, dataAtPts->getDivPAtPts()));
2220
2221 if (noStretch) {
2222 fe->getOpPtrVector().push_back(
2223 physicalEquations->returnOpCalculateStretchFromStress(
2225 } else {
2226 fe->getOpPtrVector().push_back(
2228 stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
2229 }
2230
2231 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2232 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
2233 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2234 rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
2235 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2236 spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
2237
2238 // H1 displacements
2239 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2240 spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
2241 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
2242 spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
2243
2244 // velocities
2245 if (calc_rates) {
2246 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
2247 spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
2248 if (noStretch) {
2249 } else {
2250 fe->getOpPtrVector().push_back(
2252 stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
2253 fe->getOpPtrVector().push_back(
2255 stretchTensor, dataAtPts->getGradLogStretchDotTensorAtPts(),
2256 MBTET));
2257 }
2258 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
2259 rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
2260 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradientDot<3, 3>(
2261 rotAxis, dataAtPts->getRotAxisGradDotAtPts(), MBTET));
2262
2263 // acceleration
2264 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
2265 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
2266 spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
2267 }
2268 }
2269
2270 // variations
2271 if (var_vec.use_count()) {
2272 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
2273 piolaStress, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
2274 var_vec));
2275 fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
2276 bubbleField, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
2277 var_vec, MBMAXTYPE));
2278 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2279 rotAxis, dataAtPts->getVarRotAxisPts(), var_vec, MBTET));
2280 fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
2281 piolaStress, dataAtPts->getDivVarPiolaPts(), var_vec));
2282 fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2283 spatialL2Disp, dataAtPts->getVarWL2Pts(), var_vec, MBTET));
2284
2285 if (noStretch) {
2286 fe->getOpPtrVector().push_back(
2287 physicalEquations->returnOpCalculateVarStretchFromStress(
2289 } else {
2290 fe->getOpPtrVector().push_back(
2292 stretchTensor, dataAtPts->getVarLogStreachPts(), var_vec, MBTET));
2293 }
2294 }
2295
2296 // calculate other derived quantities
2297 fe->getOpPtrVector().push_back(new OpCalculateRotationAndSpatialGradient(
2298 dataAtPts, ((do_rhs || do_lhs) && calc_rates) ? alphaOmega : 0.));
2299
2300 // evaluate integration points
2301 if (noStretch) {
2302 } else {
2303 fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
2304 tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
2305 }
2306
2308}
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
Calculate divergence of tonsorial field using vectorial base.
Calculate symmetric tensor field rates ant integratio pts.
Get field gradients time derivative at integration pts for scalar field rank 0, i....
Approximate field values for given petsc vector.

◆ setBlockTagsOnSkin()

MoFEMErrorCode EshelbianCore::setBlockTagsOnSkin ( )
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 3531 of file EshelbianPlasticity.cpp.

3531 {
3533
3534 auto set_block = [&](auto name, int dim) {
3535 std::map<int, Range> map;
3536 auto set_tag_impl = [&](auto name) {
3538 auto mesh_mng = mField.getInterface<MeshsetsManager>();
3539 auto bcs = mesh_mng->getCubitMeshsetPtr(
3540
3541 std::regex((boost::format("%s(.*)") % name).str())
3542
3543 );
3544 for (auto bc : bcs) {
3545 Range r;
3546 CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), dim, r,
3547 true);
3548 map[bc->getMeshsetId()] = r;
3549 }
3551 };
3552
3553 CHKERR set_tag_impl(name);
3554
3555 return std::make_pair(name, map);
3556 };
3557
3558 auto set_skin = [&](auto &&map) {
3559 for (auto &m : map.second) {
3560 auto s = filter_true_skin(mField, get_skin(mField, m.second));
3561 m.second.swap(s);
3562 }
3563 return map;
3564 };
3565
3566 auto set_tag = [&](auto &&map) {
3567 Tag th;
3568 auto name = map.first;
3569 int def_val[] = {-1};
3571 mField.get_moab().tag_get_handle(name, 1, MB_TYPE_INTEGER, th,
3572 MB_TAG_SPARSE | MB_TAG_CREAT, def_val),
3573 "create tag");
3574 for (auto &m : map.second) {
3575 int id = m.first;
3576 CHK_MOAB_THROW(mField.get_moab().tag_clear_data(th, m.second, &id),
3577 "clear tag");
3578 }
3579 return th;
3580 };
3581
3582 listTagsToTransfer.push_back(set_tag(set_skin(set_block("BODY", 3))));
3583 listTagsToTransfer.push_back(set_tag(set_skin(set_block("MAT_ELASTIC", 3))));
3584 listTagsToTransfer.push_back(
3585 set_tag(set_skin(set_block("MAT_NEOHOOKEAN", 3))));
3586 listTagsToTransfer.push_back(set_tag(set_block("CONTACT", 2)));
3587
3589}
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset

◆ setContactElementRhsOps()

MoFEMErrorCode EshelbianCore::setContactElementRhsOps ( boost::shared_ptr< ContactTree > &  fe_contact_tree)

Contact requires that body is marked

Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 2957 of file EshelbianPlasticity.cpp.

2961 {
2963
2964 /** Contact requires that body is marked */
2965 auto get_body_range = [this](auto name, int dim, auto sev) {
2966 std::map<int, Range> map;
2967
2968 for (auto m_ptr :
2969 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
2970
2971 (boost::format("%s(.*)") % name).str()
2972
2973 ))
2974
2975 ) {
2976 Range ents;
2977 CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
2978 dim, ents, true),
2979 "by dim");
2980 map[m_ptr->getMeshsetId()] = ents;
2981 MOFEM_LOG("EPSYNC", sev) << "Meshset: " << m_ptr->getMeshsetId() << " "
2982 << ents.size() << " entities";
2983 }
2984
2986 return map;
2987 };
2988
2989 auto get_map_skin = [this](auto &&map) {
2990 ParallelComm *pcomm =
2991 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2992
2993 Skinner skin(&mField.get_moab());
2994 for (auto &m : map) {
2995 Range skin_faces;
2996 CHKERR skin.find_skin(0, m.second, false, skin_faces);
2997 CHK_MOAB_THROW(pcomm->filter_pstatus(skin_faces,
2998 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2999 PSTATUS_NOT, -1, nullptr),
3000 "filter");
3001 m.second.swap(skin_faces);
3002 }
3003 return map;
3004 };
3005
3006 /* The above code is written in C++ and it appears to be defining and using
3007 various operations on boundary elements and side elements. */
3009 using BoundaryEleOp = BoundaryEle::UserDataOperator;
3010
3011 auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
3012
3013 auto calcs_side_traction = [&](auto &pip) {
3015 using EleOnSide =
3017 using SideEleOp = EleOnSide::UserDataOperator;
3018 auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
3019 mField, elementVolumeName, SPACE_DIM, Sev::noisy);
3020 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
3021 boost::make_shared<CGGUserPolynomialBase>();
3022 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3023 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
3025 op_loop_domain_side->getOpPtrVector().push_back(
3027 piolaStress, contact_common_data_ptr->contactTractionPtr(),
3028 boost::make_shared<double>(1.0)));
3029 pip.push_back(op_loop_domain_side);
3031 };
3032
3033 auto add_contact_three = [&]() {
3035 auto tree_moab_ptr = boost::make_shared<moab::Core>();
3036 fe_contact_tree = boost::make_shared<ContactTree>(
3037 mField, tree_moab_ptr, spaceOrder,
3038 get_body_range("CONTACT", SPACE_DIM - 1, Sev::inform));
3039 fe_contact_tree->getOpPtrVector().push_back(
3041 contactDisp, contact_common_data_ptr->contactDispPtr()));
3042 CHKERR calcs_side_traction(fe_contact_tree->getOpPtrVector());
3043 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
3044 fe_contact_tree->getOpPtrVector().push_back(
3046 fe_contact_tree->getOpPtrVector().push_back(
3047 new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
3049 };
3050
3051 CHKERR add_contact_three();
3052
3054}

◆ setElasticElementOps()

MoFEMErrorCode EshelbianCore::setElasticElementOps ( const int  tag)
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 3056 of file EshelbianPlasticity.cpp.

3056 {
3058
3059 // Add contact operators. Note that only for rhs. THe lhs is assembled with
3060 // volume element, to enable schur complement evaluation.
3062
3065
3066 auto adj_cache =
3067 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
3068
3069 auto get_op_contact_bc = [&]() {
3071 auto op_loop_side = new OpLoopSide<SideEle>(
3072 mField, contactElement, SPACE_DIM - 1, Sev::noisy, adj_cache);
3073 return op_loop_side;
3074 };
3075
3077}
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition plastic.cpp:61
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
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)
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
MoFEMErrorCode setContactElementRhsOps(boost::shared_ptr< ContactTree > &fe_contact_tree)
MoFEMErrorCode setFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)

◆ setElasticElementToTs()

MoFEMErrorCode EshelbianCore::setElasticElementToTs ( DM  dm)
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 3079 of file EshelbianPlasticity.cpp.

3079 {
3081 boost::shared_ptr<FEMethod> null;
3082
3083 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
3084
3086 null);
3088 null);
3090 null);
3092 null);
3093
3094 } else {
3096 null);
3098 null);
3100 null);
3102 null);
3103 }
3104
3106}
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition DMMoFEM.cpp:790
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition DMMoFEM.cpp:843
PetscErrorCode DMMoFEMTSSetI2Jacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition DMMoFEM.cpp:1007
PetscErrorCode DMMoFEMTSSetI2Function(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS implicit function evaluation function
Definition DMMoFEM.cpp:965

◆ 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
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 2798 of file EshelbianPlasticity.cpp.

2801 {
2803
2804 fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
2805 fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
2806
2807 // set integration rule
2808 // fe_rhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
2809 // fe_lhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
2810 fe_rhs->getRuleHook = [](int, int, int) { return -1; };
2811 fe_lhs->getRuleHook = [](int, int, int) { return -1; };
2812 fe_rhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2813 fe_lhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2814
2815 CHKERR
2816 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2817 fe_rhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
2818 CHKERR
2819 EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2820 fe_lhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
2821
2822 if (add_elastic) {
2823
2824 auto get_broken_op_side = [this](auto &pip) {
2825 using EleOnSide =
2827 using SideEleOp = EleOnSide::UserDataOperator;
2828 // Iterate over domain FEs adjacent to boundary.
2829 auto broken_data_ptr =
2830 boost::make_shared<std::vector<BrokenBaseSideData>>();
2831 // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2832 auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
2833 mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2834 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2835 boost::make_shared<CGGUserPolynomialBase>();
2836 CHKERR
2837 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2838 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2840 op_loop_domain_side->getOpPtrVector().push_back(
2841 new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
2842 boost::shared_ptr<double> piola_scale_ptr(new double);
2843 op_loop_domain_side->getOpPtrVector().push_back(
2844 physicalEquations->returnOpSetScale(piola_scale_ptr,
2846
2847 auto piola_stress_mat_ptr = boost::make_shared<MatrixDouble>();
2848 op_loop_domain_side->getOpPtrVector().push_back(
2850 piola_stress_mat_ptr));
2851 pip.push_back(op_loop_domain_side);
2852 return std::make_tuple(broken_data_ptr, piola_scale_ptr,
2853 piola_stress_mat_ptr);
2854 };
2855
2856 auto set_rhs = [&]() {
2858
2859 auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
2860 get_broken_op_side(fe_rhs->getOpPtrVector());
2861
2862 fe_rhs->getOpPtrVector().push_back(
2863 new OpDispBc(broken_data_ptr, bcSpatialDispVecPtr, timeScaleMap));
2864 fe_rhs->getOpPtrVector().push_back(new OpAnalyticalDispBc(
2866 timeScaleMap));
2867 fe_rhs->getOpPtrVector().push_back(new OpRotationBc(
2868 broken_data_ptr, bcSpatialRotationVecPtr, timeScaleMap));
2869
2870 fe_rhs->getOpPtrVector().push_back(
2872 piola_scale_ptr, timeScaleMap));
2873 auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
2874 // if you push gradient of L2 base to physical element, it will not work.
2875 fe_rhs->getOpPtrVector().push_back(
2877 hybridSpatialDisp, hybrid_grad_ptr));
2878 fe_rhs->getOpPtrVector().push_back(new OpBrokenPressureBc(
2880 hybrid_grad_ptr, timeScaleMap));
2881 fe_rhs->getOpPtrVector().push_back(new OpBrokenAnalyticalTractionBc(
2883 timeScaleMap));
2884
2885 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
2886 fe_rhs->getOpPtrVector().push_back(
2888 hybrid_ptr));
2889 fe_rhs->getOpPtrVector().push_back(new OpNormalDispRhsBc(
2890 hybridSpatialDisp, hybrid_ptr, piola_stress_mat_ptr,
2892
2893 auto get_normal_disp_bc_faces = [&]() {
2894 auto faces =
2895 get_range_from_block(mField, "NORMAL_DISPLACEMENT", SPACE_DIM - 1);
2896 return boost::make_shared<Range>(faces);
2897 };
2898
2899 using BoundaryEle =
2901 using BdyEleOp = BoundaryEle::UserDataOperator;
2903 GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
2904 fe_rhs->getOpPtrVector().push_back(new OpC_dBroken(
2905 broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0),
2906 get_normal_disp_bc_faces()));
2907
2909 };
2910
2911 auto set_lhs = [&]() {
2913
2914 auto [broken_data_ptr, piola_scale_ptr, piola_stress_mat_ptr] =
2915 get_broken_op_side(fe_lhs->getOpPtrVector());
2916
2917 fe_lhs->getOpPtrVector().push_back(new OpNormalDispLhsBc_dU(
2919 fe_lhs->getOpPtrVector().push_back(new OpNormalDispLhsBc_dP(
2921 timeScaleMap));
2922
2923 auto hybrid_grad_ptr = boost::make_shared<MatrixDouble>();
2924 // if you push gradient of L2 base to physical element, it will not work.
2925 fe_rhs->getOpPtrVector().push_back(
2927 hybridSpatialDisp, hybrid_grad_ptr));
2928 fe_lhs->getOpPtrVector().push_back(new OpBrokenPressureBcLhs_dU(
2930 timeScaleMap));
2931
2932 auto get_normal_disp_bc_faces = [&]() {
2933 auto faces =
2934 get_range_from_block(mField, "NORMAL_DISPLACEMENT", SPACE_DIM - 1);
2935 return boost::make_shared<Range>(faces);
2936 };
2937
2938 using BoundaryEle =
2940 using BdyEleOp = BoundaryEle::UserDataOperator;
2942 GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
2943 fe_lhs->getOpPtrVector().push_back(new OpC(
2944 hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0),
2945 true, true, get_normal_disp_bc_faces()));
2946
2948 };
2949
2950 CHKERR set_rhs();
2951 CHKERR set_lhs();
2952 }
2953
2955}
@ GAUSS
Gaussian quadrature integration.
constexpr AssemblyType A
Apply rotation boundary condition.
BoundaryEle::UserDataOperator BdyEleOp

◆ setNewFrontCoordinates()

MoFEMErrorCode EshelbianCore::setNewFrontCoordinates ( )
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 5842 of file EshelbianPlasticity.cpp.

5842 {
5844
5845 if (!maxMovedFaces)
5847
5848 Tag th_front_position;
5849 auto rval =
5850 mField.get_moab().tag_get_handle("FrontPosition", th_front_position);
5851 if (rval == MB_SUCCESS && maxMovedFaces) {
5852 Range verts;
5853 CHKERR mField.get_moab().get_connectivity(*maxMovedFaces, verts, true);
5854 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
5855 std::vector<double> coords(3 * verts.size());
5856 CHKERR mField.get_moab().get_coords(verts, coords.data());
5857 std::vector<double> pos(3 * verts.size());
5858 CHKERR mField.get_moab().tag_get_data(th_front_position, verts, pos.data());
5859 for (int i = 0; i != 3 * verts.size(); ++i) {
5860 coords[i] += pos[i];
5861 }
5862 CHKERR mField.get_moab().set_coords(verts, coords.data());
5863 double zero[] = {0., 0., 0.};
5864 CHKERR mField.get_moab().tag_clear_data(th_front_position, verts, zero);
5865 }
5866
5867#ifndef NDEBUG
5868 constexpr bool debug = false;
5869 if (debug) {
5870
5872 mField.get_moab(),
5873 "set_coords_faces_" +
5874 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
5875 *maxMovedFaces);
5876 }
5877#endif
5879}

◆ 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
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 2310 of file EshelbianPlasticity.cpp.

2313 {
2315
2316 /** Contact requires that body is marked */
2317 auto get_body_range = [this](auto name, int dim) {
2318 std::map<int, Range> map;
2319
2320 for (auto m_ptr :
2321 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
2322
2323 (boost::format("%s(.*)") % name).str()
2324
2325 ))
2326
2327 ) {
2328 Range ents;
2329 CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
2330 dim, ents, true),
2331 "by dim");
2332 map[m_ptr->getMeshsetId()] = ents;
2333 }
2334
2335 return map;
2336 };
2337
2338 auto rule_contact = [](int, int, int o) { return -1; };
2340
2341 auto set_rule_contact = [refine](
2342
2343 ForcesAndSourcesCore *fe_raw_ptr, int order_row,
2344 int order_col, int order_data
2345
2346 ) {
2348 auto rule = 2 * order_data;
2349 fe_raw_ptr->gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
2351 };
2352
2353 // Right hand side
2354 fe_rhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
2355 CHKERR setBaseVolumeElementOps(tag, true, false, true, SmartPetscObj<Vec>(),
2356 fe_rhs);
2357
2358 // elastic
2359 if (add_elastic) {
2360
2361 fe_rhs->getOpPtrVector().push_back(
2363 fe_rhs->getOpPtrVector().push_back(
2365 if (noStretch) {
2366 // do nothing - no stretch approximation
2367 } else {
2368 if (!internalStressTagName.empty()) {
2369 switch (internalStressInterpOrder) {
2370 case 0:
2371 fe_rhs->getOpPtrVector().push_back(
2373 break;
2374 case 1:
2375 fe_rhs->getOpPtrVector().push_back(
2377 break;
2378 default:
2379 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
2380 "Unsupported internal stress interpolation order %d",
2382 }
2383 if (internalStressVoigt) {
2384 fe_rhs->getOpPtrVector().push_back(
2386 dataAtPts));
2387 } else {
2388 fe_rhs->getOpPtrVector().push_back(
2390 dataAtPts));
2391 }
2392 }
2393 if (auto op = physicalEquations->returnOpSpatialPhysicalExternalStrain(
2395 fe_rhs->getOpPtrVector().push_back(op);
2396 } else if (externalStrainVecPtr && !externalStrainVecPtr->empty()) {
2397 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
2398 "OpSpatialPhysicalExternalStrain not implemented for this "
2399 "material");
2400 }
2401
2402 fe_rhs->getOpPtrVector().push_back(
2403 physicalEquations->returnOpSpatialPhysical(stretchTensor, dataAtPts,
2404 alphaU));
2405 }
2406 fe_rhs->getOpPtrVector().push_back(
2408 fe_rhs->getOpPtrVector().push_back(
2410 fe_rhs->getOpPtrVector().push_back(
2412
2413 auto set_hybridisation = [&](auto &pip) {
2415
2416 using BoundaryEle =
2418 using EleOnSide =
2420 using SideEleOp = EleOnSide::UserDataOperator;
2421 using BdyEleOp = BoundaryEle::UserDataOperator;
2422
2423 // First: Iterate over skeleton FEs adjacent to Domain FEs
2424 // Note: BoundaryEle, i.e. uses skeleton interation rule
2425 auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
2426 mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
2427 // op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
2428 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2429 return -1;
2430 };
2431 op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2432 SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2433
2434 CHKERR EshelbianPlasticity::
2435 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2436 op_loop_skeleton_side->getOpPtrVector(), {L2},
2438
2439 // Second: Iterate over domain FEs adjacent to skelton, particularly one
2440 // domain element.
2441 auto broken_data_ptr =
2442 boost::make_shared<std::vector<BrokenBaseSideData>>();
2443 // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2444 auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
2445 mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2446 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2447 boost::make_shared<CGGUserPolynomialBase>();
2448 CHKERR
2449 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2450 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2452 op_loop_domain_side->getOpPtrVector().push_back(
2453 new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
2454 auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
2455 op_loop_domain_side->getOpPtrVector().push_back(
2457 flux_mat_ptr));
2458 op_loop_domain_side->getOpPtrVector().push_back(
2459 new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
2460
2461 // Assemble on skeleton
2462 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2464 GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
2466 GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
2467 op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dHybrid(
2468 hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0)));
2469 auto hybrid_ptr = boost::make_shared<MatrixDouble>();
2470 op_loop_skeleton_side->getOpPtrVector().push_back(
2472 hybrid_ptr));
2473 op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dBroken(
2474 broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0)));
2475
2476 // Add skeleton to domain pipeline
2477 pip.push_back(op_loop_skeleton_side);
2478
2480 };
2481
2482 auto set_contact = [&](auto &pip) {
2484
2485 using BoundaryEle =
2487 using EleOnSide =
2489 using SideEleOp = EleOnSide::UserDataOperator;
2490 using BdyEleOp = BoundaryEle::UserDataOperator;
2491
2492 // First: Iterate over skeleton FEs adjacent to Domain FEs
2493 // Note: BoundaryEle, i.e. uses skeleton interation rule
2494 auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
2495 mField, contactElement, SPACE_DIM - 1, Sev::noisy);
2496
2497 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
2498 op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
2499 CHKERR EshelbianPlasticity::
2500 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2501 op_loop_skeleton_side->getOpPtrVector(), {L2},
2503
2504 // Second: Iterate over domain FEs adjacent to skelton, particularly
2505 // one domain element.
2506 auto broken_data_ptr =
2507 boost::make_shared<std::vector<BrokenBaseSideData>>();
2508
2509 // Data storing contact fields
2510 auto contact_common_data_ptr =
2511 boost::make_shared<ContactOps::CommonData>();
2512
2513 auto add_ops_domain_side = [&](auto &pip) {
2515 // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2516 auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
2517 mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2518 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2519 boost::make_shared<CGGUserPolynomialBase>();
2520 CHKERR
2521 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2522 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2524 op_loop_domain_side->getOpPtrVector().push_back(
2526 broken_data_ptr));
2527 op_loop_domain_side->getOpPtrVector().push_back(
2529 piolaStress, contact_common_data_ptr->contactTractionPtr()));
2530 pip.push_back(op_loop_domain_side);
2532 };
2533
2534 auto add_ops_contact_rhs = [&](auto &pip) {
2536 // get body id and SDF range
2537 auto contact_sfd_map_range_ptr =
2538 boost::make_shared<std::map<int, Range>>(
2539 get_body_range("CONTACT_SDF", SPACE_DIM - 1));
2540
2541 pip.push_back(new OpCalculateVectorFieldValues<3>(
2542 contactDisp, contact_common_data_ptr->contactDispPtr()));
2543 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2544 pip.push_back(
2546 pip.push_back(new OpTreeSearch(
2547 contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
2548 get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
2549 nullptr));
2551 contactDisp, contact_common_data_ptr, contactTreeRhs,
2552 contact_sfd_map_range_ptr));
2553 pip.push_back(
2555 broken_data_ptr, contact_common_data_ptr, contactTreeRhs));
2556
2558 };
2559
2560 // push ops to face/side pipeline
2561 CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
2562 CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
2563
2564 // Add skeleton to domain pipeline
2565 pip.push_back(op_loop_skeleton_side);
2566
2568 };
2569
2570 CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
2571 CHKERR set_contact(fe_rhs->getOpPtrVector());
2572
2573 // Body forces
2574 using BodyNaturalBC =
2576 Assembly<PETSC>::LinearForm<GAUSS>;
2577 using OpBodyForce =
2578 BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
2579
2580 auto body_time_scale =
2581 boost::make_shared<DynamicRelaxationTimeScale>("body_force.txt");
2582 CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
2583 fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {body_time_scale},
2584 "BODY_FORCE", Sev::inform);
2585 }
2586
2587 // Left hand side
2588 fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
2589 CHKERR setBaseVolumeElementOps(tag, true, true, true, SmartPetscObj<Vec>(),
2590 fe_lhs);
2591
2592 // elastic
2593 if (add_elastic) {
2594
2595 if (noStretch) {
2596 fe_lhs->getOpPtrVector().push_back(
2598 fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dP(
2600 fe_lhs->getOpPtrVector().push_back(
2602 dataAtPts));
2603 } else {
2604 fe_lhs->getOpPtrVector().push_back(
2605 physicalEquations->returnOpSpatialPhysical_du_du(
2607 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
2609 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
2611 fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
2613 symmetrySelector == SYMMETRIC ? true : false));
2614 }
2615
2616 fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
2618 fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
2620
2621 fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
2623 symmetrySelector == SYMMETRIC ? true : false));
2624 fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
2626 symmetrySelector == SYMMETRIC ? true : false));
2627
2628 if (symmetrySelector > SYMMETRIC) {
2629 if (!noStretch) {
2630 fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_du(
2631 rotAxis, stretchTensor, dataAtPts, false));
2632 }
2633 fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
2634 rotAxis, piolaStress, dataAtPts, false));
2635 fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
2636 rotAxis, bubbleField, dataAtPts, false));
2637 }
2638 fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_domega(
2640
2641 auto set_hybridisation = [&](auto &pip) {
2643
2644 using BoundaryEle =
2646 using EleOnSide =
2648 using SideEleOp = EleOnSide::UserDataOperator;
2649 using BdyEleOp = BoundaryEle::UserDataOperator;
2650
2651 // First: Iterate over skeleton FEs adjacent to Domain FEs
2652 // Note: BoundaryEle, i.e. uses skeleton interation rule
2653 auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
2654 mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
2655 // op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
2656 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2657 return -1;
2658 };
2659 op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2660 SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2661 CHKERR EshelbianPlasticity::
2662 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2663 op_loop_skeleton_side->getOpPtrVector(), {L2},
2665
2666 // Second: Iterate over domain FEs adjacent to skelton, particularly one
2667 // domain element.
2668 auto broken_data_ptr =
2669 boost::make_shared<std::vector<BrokenBaseSideData>>();
2670 // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2671 auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
2672 mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2673 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2674 boost::make_shared<CGGUserPolynomialBase>();
2675 CHKERR
2676 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2677 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2679 op_loop_domain_side->getOpPtrVector().push_back(
2680 new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
2681
2682 op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2684 GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
2685 op_loop_skeleton_side->getOpPtrVector().push_back(
2686 new OpC(hybridSpatialDisp, broken_data_ptr,
2687 boost::make_shared<double>(1.0), true, false));
2688
2689 pip.push_back(op_loop_skeleton_side);
2690
2692 };
2693
2694 auto set_contact = [&](auto &pip) {
2696
2697 using BoundaryEle =
2699 using EleOnSide =
2701 using SideEleOp = EleOnSide::UserDataOperator;
2702 using BdyEleOp = BoundaryEle::UserDataOperator;
2703
2704 // First: Iterate over skeleton FEs adjacent to Domain FEs
2705 // Note: BoundaryEle, i.e. uses skeleton interation rule
2706 auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
2707 mField, contactElement, SPACE_DIM - 1, Sev::noisy);
2708
2709 op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
2710 op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
2711 CHKERR EshelbianPlasticity::
2712 AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2713 op_loop_skeleton_side->getOpPtrVector(), {L2},
2715
2716 // Second: Iterate over domain FEs adjacent to skelton, particularly
2717 // one domain element.
2718 auto broken_data_ptr =
2719 boost::make_shared<std::vector<BrokenBaseSideData>>();
2720
2721 // Data storing contact fields
2722 auto contact_common_data_ptr =
2723 boost::make_shared<ContactOps::CommonData>();
2724
2725 auto add_ops_domain_side = [&](auto &pip) {
2727 // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2728 auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
2729 mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2730 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2731 boost::make_shared<CGGUserPolynomialBase>();
2732 CHKERR
2733 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2734 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2736 op_loop_domain_side->getOpPtrVector().push_back(
2738 broken_data_ptr));
2739 op_loop_domain_side->getOpPtrVector().push_back(
2741 piolaStress, contact_common_data_ptr->contactTractionPtr()));
2742 pip.push_back(op_loop_domain_side);
2744 };
2745
2746 auto add_ops_contact_lhs = [&](auto &pip) {
2748 pip.push_back(new OpCalculateVectorFieldValues<3>(
2749 contactDisp, contact_common_data_ptr->contactDispPtr()));
2750 auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2751 pip.push_back(
2753 pip.push_back(new OpTreeSearch(
2754 contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
2755 get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
2756 nullptr));
2757
2758 // get body id and SDF range
2759 auto contact_sfd_map_range_ptr =
2760 boost::make_shared<std::map<int, Range>>(
2761 get_body_range("CONTACT_SDF", SPACE_DIM - 1));
2762
2764 contactDisp, contactDisp, contact_common_data_ptr, contactTreeRhs,
2765 contact_sfd_map_range_ptr));
2766 pip.push_back(
2768 contactDisp, broken_data_ptr, contact_common_data_ptr,
2769 contactTreeRhs, contact_sfd_map_range_ptr));
2770 pip.push_back(
2772 broken_data_ptr, contactDisp, contact_common_data_ptr,
2774
2776 };
2777
2778 // push ops to face/side pipeline
2779 CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
2780 CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
2781
2782 // Add skeleton to domain pipeline
2783 pip.push_back(op_loop_skeleton_side);
2784
2786 };
2787
2788 CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
2789 CHKERR set_contact(fe_lhs->getOpPtrVector());
2790 }
2791
2792 if (add_material) {
2793 }
2794
2796}
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)
structure to get information from mofem into EntitiesFieldData
MatrixDouble gaussPts
Matrix of integration points.
Natural boundary conditions.
Definition Natural.hpp:57
Operator for broken loop side.
static RefineTrianglesReturn refineTriangle(int nb_levels)
create uniform triangle mesh of refined elements
Definition Tools.cpp:724
static MatrixDouble refineTriangleIntegrationPts(MatrixDouble pts, RefineTrianglesReturn refined)
generate integration points for refined triangle mesh for last level
Definition Tools.cpp:791

◆ solveDynamicRelaxation()

MoFEMErrorCode EshelbianCore::solveDynamicRelaxation ( TS  ts,
Vec  x 
)
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 3440 of file EshelbianPlasticity.cpp.

3440 {
3442
3443 auto storage = solve_elastic_setup::setup(this, ts, x, false);
3444
3445 double final_time = 1;
3446 double delta_time = 0.1;
3447 int max_it = 10;
3448 PetscBool ts_h1_update = PETSC_FALSE;
3449
3450 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options",
3451 "none");
3452
3453 CHKERR PetscOptionsScalar("-dynamic_final_time",
3454 "dynamic relaxation final time", "", final_time,
3455 &final_time, PETSC_NULLPTR);
3456 CHKERR PetscOptionsScalar("-dynamic_delta_time",
3457 "dynamic relaxation final time", "", delta_time,
3458 &delta_time, PETSC_NULLPTR);
3459 CHKERR PetscOptionsInt("-dynamic_max_it", "dynamic relaxation iterations", "",
3460 max_it, &max_it, PETSC_NULLPTR);
3461 CHKERR PetscOptionsBool("-dynamic_h1_update", "update each ts step", "",
3462 ts_h1_update, &ts_h1_update, PETSC_NULLPTR);
3463
3464 PetscOptionsEnd();
3465
3466 auto setup_ts_monitor = [&]() {
3467 auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
3468 return monitor_ptr;
3469 };
3470 auto monitor_ptr = setup_ts_monitor();
3471
3472 TetPolynomialBase::switchCacheBaseOn<HDIV>(
3473 {elasticFeLhs.get(), elasticFeRhs.get()});
3474 CHKERR TSSetUp(ts);
3476
3477 double ts_delta_time;
3478 CHKERR TSGetTimeStep(ts, &ts_delta_time);
3479
3480 if (ts_h1_update) {
3481 CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
3482 CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
3483 }
3484
3487
3488 dynamicTime = 0;
3489 dynamicStep = 0;
3490 monitor_ptr->ts_u = PETSC_NULLPTR;
3491 monitor_ptr->ts_t = dynamicTime;
3492 monitor_ptr->ts_step = dynamicStep;
3494 MOFEM_LOG("EP", Sev::inform) << "IT " << dynamicStep << " Time "
3495 << dynamicTime << " delta time " << delta_time;
3496 dynamicTime += delta_time;
3497 ++dynamicStep;
3498
3499 for (; dynamicTime < final_time; dynamicTime += delta_time) {
3500 MOFEM_LOG("EP", Sev::inform) << "IT " << dynamicStep << " Time "
3501 << dynamicTime << " delta time " << delta_time;
3502
3503 CHKERR TSSetStepNumber(ts, 0);
3504 CHKERR TSSetTime(ts, 0);
3505 CHKERR TSSetTimeStep(ts, ts_delta_time);
3506 if (!ts_h1_update) {
3508 }
3509 CHKERR TSSolve(ts, PETSC_NULLPTR);
3510 if (!ts_h1_update) {
3512 }
3513
3514 monitor_ptr->ts_u = PETSC_NULLPTR;
3515 monitor_ptr->ts_t = dynamicTime;
3516 monitor_ptr->ts_step = dynamicStep;
3518
3519 ++dynamicStep;
3520 if (dynamicStep > max_it)
3521 break;
3522 }
3523
3525 TetPolynomialBase::switchCacheBaseOff<HDIV>(
3526 {elasticFeLhs.get(), elasticFeRhs.get()});
3527
3529}
static double dynamicTime
static int dynamicStep
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
static MoFEMErrorCode postStepDestroy()
static MoFEMErrorCode preStepFun(TS ts)
static MoFEMErrorCode postStepFun(TS ts)

◆ solveElastic()

MoFEMErrorCode EshelbianCore::solveElastic ( TS  ts,
Vec  x 
)
Examples
ep.cpp, mofem/users_modules/eshelbian_plasticity/ep.cpp, and mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 3310 of file EshelbianPlasticity.cpp.

3310 {
3312
3313 PetscBool debug_model = PETSC_FALSE;
3314 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-debug_model", &debug_model,
3315 PETSC_NULLPTR);
3316
3317 if (debug_model == PETSC_TRUE) {
3318 auto ts_ctx_ptr = getDMTsCtx(dmElastic);
3319 auto post_proc = [&](TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, Vec F,
3320 void *ctx) {
3322
3323 SNES snes;
3324 CHKERR TSGetSNES(ts, &snes);
3325 int it;
3326 CHKERR SNESGetIterationNumber(snes, &it);
3327 std::string file_name = "snes_iteration_" + std::to_string(it) + ".h5m";
3328 CHKERR postProcessResults(1, file_name, F, u_t);
3329 std::string file_skel_name =
3330 "snes_iteration_skel_" + std::to_string(it) + ".h5m";
3331
3332 auto get_material_force_tag = [&]() {
3333 auto &moab = mField.get_moab();
3334 Tag tag;
3335 CHK_MOAB_THROW(moab.tag_get_handle("MaterialForce", tag),
3336 "can't get tag");
3337 return tag;
3338 };
3339
3341 CHKERR postProcessSkeletonResults(1, file_skel_name, F,
3342 {get_material_force_tag()});
3343
3345 };
3346 ts_ctx_ptr->tsDebugHook = post_proc;
3347 }
3348
3349 auto storage = solve_elastic_setup::setup(this, ts, x, true);
3350
3351 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
3352 Vec xx;
3353 CHKERR VecDuplicate(x, &xx);
3354 CHKERR VecZeroEntries(xx);
3355 CHKERR TS2SetSolution(ts, x, xx);
3356 CHKERR VecDestroy(&xx);
3357 } else {
3358 CHKERR TSSetSolution(ts, x);
3359 }
3360
3361 TetPolynomialBase::switchCacheBaseOn<HDIV>(
3362 {elasticFeLhs.get(), elasticFeRhs.get()});
3363 CHKERR TSSetUp(ts);
3364 CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
3365 CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
3367 CHKERR TSSolve(ts, PETSC_NULLPTR);
3369 TetPolynomialBase::switchCacheBaseOff<HDIV>(
3370 {elasticFeLhs.get(), elasticFeRhs.get()});
3371
3372 SNES snes;
3373 CHKERR TSGetSNES(ts, &snes);
3374 int lin_solver_iterations;
3375 CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
3376 MOFEM_LOG("EP", Sev::inform)
3377 << "Number of linear solver iterations " << lin_solver_iterations;
3378
3379 PetscBool test_cook_flg = PETSC_FALSE;
3380 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_cook", &test_cook_flg,
3381 PETSC_NULLPTR);
3382 if (test_cook_flg) {
3383 constexpr int expected_lin_solver_iterations = 11;
3384 if (lin_solver_iterations > expected_lin_solver_iterations)
3385 SETERRQ(
3386 PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
3387 "Expected number of iterations is different than expected %d > %d",
3388 lin_solver_iterations, expected_lin_solver_iterations);
3389 }
3390
3391 PetscBool test_sslv116_flag = PETSC_FALSE;
3392 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_sslv116",
3393 &test_sslv116_flag, PETSC_NULLPTR);
3394
3395 if (test_sslv116_flag) {
3396 double max_val = 0.0;
3397 double min_val = 0.0;
3398 auto field_min_max = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
3400 auto ent_type = ent_ptr->getEntType();
3401 if (ent_type == MBVERTEX) {
3402 max_val = std::max(ent_ptr->getEntFieldData()[SPACE_DIM - 1], max_val);
3403 min_val = std::min(ent_ptr->getEntFieldData()[SPACE_DIM - 1], min_val);
3404 }
3406 };
3407 CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
3408 field_min_max, spatialH1Disp);
3409
3410 double global_max_val = 0.0;
3411 double global_min_val = 0.0;
3412 MPI_Allreduce(&max_val, &global_max_val, 1, MPI_DOUBLE, MPI_MAX,
3413 mField.get_comm());
3414 MPI_Allreduce(&min_val, &global_min_val, 1, MPI_DOUBLE, MPI_MIN,
3415 mField.get_comm());
3416 MOFEM_LOG("EP", Sev::inform)
3417 << "Max " << spatialH1Disp << " value: " << global_max_val;
3418 MOFEM_LOG("EP", Sev::inform)
3419 << "Min " << spatialH1Disp << " value: " << global_min_val;
3420
3421 double ref_max_val = 0.00767;
3422 double ref_min_val = -0.00329;
3423 if (std::abs(global_max_val - ref_max_val) > 1e-5) {
3424 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
3425 "Incorrect max value of the displacement field: %f != %f",
3426 global_max_val, ref_max_val);
3427 }
3428 if (std::abs(global_min_val - ref_min_val) > 4e-5) {
3429 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
3430 "Incorrect min value of the displacement field: %f != %f",
3431 global_min_val, ref_min_val);
3432 }
3433 }
3434
3436
3438}
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ F
const FTensor::Tensor2< T, Dim, Dim > Vec
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1276
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode calculateFaceMaterialForce(const int tag, TS ts)
MoFEMErrorCode gettingNorms()
[Getting norms]
MoFEMErrorCode postProcessSkeletonResults(const int tag, const std::string file, Vec f_residual=PETSC_NULLPTR, std::vector< Tag > tags_to_transfer={})
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={})
Basic algebra on fields.
Definition FieldBlas.hpp:21

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

◆ alphaOmega

double EshelbianCore::alphaOmega = 0

◆ alphaRho

double EshelbianCore::alphaRho = 0

◆ alphaU

double EshelbianCore::alphaU = 0

◆ alphaW

double EshelbianCore::alphaW = 0

◆ AnalyticalExprPythonPtr

boost::shared_ptr<AnalyticalExprPython> EshelbianCore::AnalyticalExprPythonPtr

◆ aoS

AO EshelbianCore::aoS = PETSC_NULLPTR

Definition at line 444 of file EshelbianCore.hpp.

◆ bcSpatialAnalyticalDisplacementVecPtr

boost::shared_ptr<AnalyticalDisplacementBcVec> EshelbianCore::bcSpatialAnalyticalDisplacementVecPtr

◆ bcSpatialAnalyticalTractionVecPtr

boost::shared_ptr<AnalyticalTractionBcVec> EshelbianCore::bcSpatialAnalyticalTractionVecPtr

◆ bcSpatialDispVecPtr

boost::shared_ptr<BcDispVec> EshelbianCore::bcSpatialDispVecPtr

◆ bcSpatialFreeTractionVecPtr

boost::shared_ptr<TractionFreeBc> EshelbianCore::bcSpatialFreeTractionVecPtr

◆ bcSpatialNormalDisplacementVecPtr

boost::shared_ptr<NormalDisplacementBcVec> EshelbianCore::bcSpatialNormalDisplacementVecPtr

◆ bcSpatialPressureVecPtr

boost::shared_ptr<PressureBcVec> EshelbianCore::bcSpatialPressureVecPtr

◆ bcSpatialRotationVecPtr

boost::shared_ptr<BcRotVec> EshelbianCore::bcSpatialRotationVecPtr

◆ bcSpatialTractionVecPtr

boost::shared_ptr<TractionBcVec> EshelbianCore::bcSpatialTractionVecPtr

◆ bitAdjEnt

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

bit ref level for parent

Examples
mofem/users_modules/eshelbian_plasticity/src/impl/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
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 430 of file EshelbianCore.hpp.

◆ bitAdjParent

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

bit ref level for parent

Examples
mofem/users_modules/eshelbian_plasticity/src/impl/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
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 427 of file EshelbianCore.hpp.

◆ bubbleField

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

◆ contactDisp

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

◆ contactElement

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

◆ contactFaces

boost::shared_ptr<Range> EshelbianCore::contactFaces

◆ contactRefinementLevels

int EshelbianCore::contactRefinementLevels = 1

◆ contactTreeRhs

boost::shared_ptr<ContactTree> EshelbianCore::contactTreeRhs

Make a contact tree.

Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 185 of file EshelbianCore.hpp.

◆ crackFaces

boost::shared_ptr<Range> EshelbianCore::crackFaces

◆ crackHybridIs

SmartPetscObj<IS> EshelbianCore::crackHybridIs

Definition at line 445 of file EshelbianCore.hpp.

◆ crackingOn

PetscBool EshelbianCore::crackingOn = PETSC_FALSE
inlinestatic

◆ crackingStartTime

double EshelbianCore::crackingStartTime = 0
inlinestatic

◆ d_f

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

◆ dataAtPts

boost::shared_ptr<DataAtIntegrationPts> EshelbianCore::dataAtPts

◆ dd_f

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

◆ dM

SmartPetscObj<DM> EshelbianCore::dM

Coupled problem all fields.

Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 187 of file EshelbianCore.hpp.

◆ dmElastic

SmartPetscObj<DM> EshelbianCore::dmElastic

◆ dmMaterial

SmartPetscObj<DM> EshelbianCore::dmMaterial

Material problem.

Definition at line 189 of file EshelbianCore.hpp.

◆ dmPrjSpatial

SmartPetscObj<DM> EshelbianCore::dmPrjSpatial

Projection spatial displacement.

Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 190 of file EshelbianCore.hpp.

◆ dynamicRelaxation

PetscBool EshelbianCore::dynamicRelaxation
inlinestatic

◆ dynamicStep

int EshelbianCore::dynamicStep
inlinestatic

◆ dynamicTime

double EshelbianCore::dynamicTime
inlinestatic

◆ edgeExchange

CommInterface::EntitiesPetscVector EshelbianCore::edgeExchange

◆ elasticBcLhs

boost::shared_ptr<FaceElementForcesAndSourcesCore> EshelbianCore::elasticBcLhs

◆ elasticBcRhs

boost::shared_ptr<FaceElementForcesAndSourcesCore> EshelbianCore::elasticBcRhs

◆ elasticFeLhs

boost::shared_ptr<VolumeElementForcesAndSourcesCore> EshelbianCore::elasticFeLhs

◆ elasticFeRhs

boost::shared_ptr<VolumeElementForcesAndSourcesCore> EshelbianCore::elasticFeRhs

◆ elementVolumeName

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

◆ energyReleaseSelector

enum EnergyReleaseSelector EshelbianCore::energyReleaseSelector
inlinestatic

◆ exponentBase

double EshelbianCore::exponentBase = exp(1)
static

◆ externalStrainVecPtr

boost::shared_ptr<ExternalStrainVec> EshelbianCore::externalStrainVecPtr

◆ f

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

◆ faceExchange

CommInterface::EntitiesPetscVector EshelbianCore::faceExchange

◆ frontAdjEdges

boost::shared_ptr<Range> EshelbianCore::frontAdjEdges

◆ frontEdges

boost::shared_ptr<Range> EshelbianCore::frontEdges

◆ frontLayers

int EshelbianCore::frontLayers = 3

Definition at line 225 of file EshelbianCore.hpp.

◆ frontVertices

boost::shared_ptr<Range> EshelbianCore::frontVertices

◆ gradApproximator

enum RotSelector EshelbianCore::gradApproximator = LARGE_ROT
inlinestatic

◆ griffithEnergy

double EshelbianCore::griffithEnergy = 1
inlinestatic

◆ hybridSpatialDisp

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

◆ internalStressInterpOrder

int EshelbianCore::internalStressInterpOrder
inlinestatic

◆ internalStressTagName

std::string EshelbianCore::internalStressTagName
inlinestatic

◆ internalStressVoigt

PetscBool EshelbianCore::internalStressVoigt
inlinestatic
Initial value:
=
PETSC_FALSE
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/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:
=
static double inv_d_f_log_e(const double v)
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/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:
=
static double inv_dd_f_log_e(const double v)
Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 48 of file EshelbianCore.hpp.

◆ inv_f

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

◆ l2UserBaseScale

PetscBool EshelbianCore::l2UserBaseScale = PETSC_TRUE
inlinestatic

◆ listTagsToTransfer

std::vector<Tag> EshelbianCore::listTagsToTransfer

list of tags to transfer to postprocessor

Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 441 of file EshelbianCore.hpp.

◆ materialH1Positions

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

◆ materialOrder

int EshelbianCore::materialOrder = 1

◆ materialSkeletonFaces

boost::shared_ptr<Range> EshelbianCore::materialSkeletonFaces

Definition at line 420 of file EshelbianCore.hpp.

◆ maxMovedFaces

boost::shared_ptr<Range> EshelbianCore::maxMovedFaces

◆ mField

MoFEM::Interface& EshelbianCore::mField

◆ naturalBcElement

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

◆ nbCrackFaces

int EshelbianCore::nbCrackFaces = 0

◆ nbJIntegralLevels

int EshelbianCore::nbJIntegralLevels
inlinestatic

◆ noStretch

PetscBool EshelbianCore::noStretch = PETSC_FALSE
inlinestatic

◆ parentAdjSkeletonFunctionDim2

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

◆ physicalEquations

boost::shared_ptr<PhysicalEquations> EshelbianCore::physicalEquations

◆ piolaStress

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

◆ rotAxis

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

◆ rotSelector

enum RotSelector EshelbianCore::rotSelector = LARGE_ROT
inlinestatic

◆ S

Mat EshelbianCore::S = PETSC_NULLPTR

Definition at line 443 of file EshelbianCore.hpp.

◆ setSingularity

PetscBool EshelbianCore::setSingularity = PETSC_TRUE
inlinestatic

◆ skeletonElement

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

◆ skeletonFaces

boost::shared_ptr<Range> EshelbianCore::skeletonFaces

◆ skinElement

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

◆ solTSStep

SmartPetscObj<Vec> EshelbianCore::solTSStep

◆ spaceH1Order

int EshelbianCore::spaceH1Order = -1

◆ spaceOrder

int EshelbianCore::spaceOrder = 2

◆ spatialH1Disp

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

◆ spatialL2Disp

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

◆ stretchSelector

enum StretchSelector EshelbianCore::stretchSelector = LOG
inlinestatic

◆ stretchTensor

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

◆ symmetrySelector

enum SymmetrySelector EshelbianCore::symmetrySelector = NOT_SYMMETRIC
inlinestatic

◆ timeScaleMap

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

◆ 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

◆ volumeExchange

CommInterface::EntitiesPetscVector EshelbianCore::volumeExchange

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