v0.14.0
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_plasticit/src/EshelbianCore.hpp>

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

Classes

struct  SetUpSchur
 

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 Getting interface of core database. More...
 
 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] More...
 
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. More...
 
MoFEMErrorCode getSpatialTractionFreeBc (const EntityHandle meshset=0)
 
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_NULL, Vec var_vec=PETSC_NULL, std::vector< Tag > tags_to_transfer={})
 
MoFEMErrorCode postProcessSkeletonResults (const int tag, const std::string file, Vec f_residual=PETSC_NULL)
 
MoFEMErrorCode gettingNorms ()
 [Getting norms] More...
 
MoFEMErrorCode calculateEnergyRelease (const int tag, TS ts)
 
MoFEMErrorCode calculateOrientation (const int tag, bool set_orientation)
 
MoFEMErrorCode setNewFrontCoordinates ()
 
MoFEMErrorCode addCrackSurfaces ()
 
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. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface reference to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

Static Public Member Functions

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. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

Public Attributes

MoFEM::Interface & mField
 
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
 
boost::shared_ptr< PhysicalEquations > physicalEquations
 
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
 
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
 
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
 
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
 
boost::shared_ptr< ContactTree > contactTreeRhs
 Make a contact tree. More...
 
SmartPetscObj< DM > dM
 Coupled problem all fields. More...
 
SmartPetscObj< DM > dmElastic
 Elastic problem. More...
 
SmartPetscObj< DM > dmMaterial
 Material problem. More...
 
SmartPetscObj< DM > dmPrjSpatial
 Projection spatial displacement. More...
 
const std::string piolaStress = "P"
 
const std::string eshelbyStress = "S"
 
const std::string spatialL2Disp = "wL2"
 
const std::string materialL2Disp = "WL2"
 
const std::string spatialH1Disp = "wH1"
 
const std::string materialH1Positions = "XH1"
 
const std::string hybridSpatialDisp = "hybridSpatialDisp"
 
const std::string hybridMaterialDisp = "hybridMaterialDisp"
 
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"
 
const std::string materialVolumeElement = "MEP"
 
const std::string materialSkeletonElement = "MATERIAL_SKELETON"
 
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 > bcSpatialTraction
 
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
 
boost::shared_ptr< Range > contactFaces
 
boost::shared_ptr< Range > crackFaces
 
boost::shared_ptr< Range > frontEdges
 
boost::shared_ptr< Range > frontAdjEdges
 
boost::shared_ptr< Range > frontVertices
 
boost::shared_ptr< Range > skeletonFaces
 
boost::shared_ptr< Range > materialSkeletonFaces
 
boost::shared_ptr< Range > maxMovedFaces
 
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > > parentAdjSkeletonFunctionDim2
 
BitRefLevel bitAdjParent = BitRefLevel().set()
 bit ref level for parent More...
 
BitRefLevel bitAdjParentMask
 bit ref level for parent parent More...
 
BitRefLevel bitAdjEnt = BitRefLevel().set()
 bit ref level for parent More...
 
BitRefLevel bitAdjEntMask
 bit ref level for parent parent More...
 
SmartPetscObj< Vec > solTSStep
 
CommInterface::EntitiesPetscVector volumeExchange
 
CommInterface::EntitiesPetscVector faceExchange
 
CommInterface::EntitiesPetscVector edgeExchange
 
CommInterface::EntitiesPetscVector vertexExchange
 
std::vector< Tag > listTagsToTransfer
 list of tags to transfer to postprocessor More...
 
Mat S = PETSC_NULL
 
AO aoS = PETSC_NULL
 
SmartPetscObj< IS > crackHybridIs
 
std::vector< std::string > a00FieldList
 
std::vector< boost::shared_ptr< Range > > a00RangeList
 

Static Public Attributes

static enum SymmetrySelector symmetrySelector = 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 dynamicTime
 
static int dynamicStep
 
static PetscBool l2UserBaseScale = PETSC_FALSE
 
static int addCrackMeshsetId = 1000
 
static double griffithEnergy = 1
 Griffith energy. More...
 
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
 

Friends

struct solve_elastic_set_up
 

Detailed Description

Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 12 of file EshelbianCore.hpp.

Constructor & Destructor Documentation

◆ EshelbianCore()

EshelbianCore::EshelbianCore ( MoFEM::Interface &  m_field)
Examples
EshelbianPlasticity.cpp.

Definition at line 877 of file EshelbianPlasticity.cpp.

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

◆ ~EshelbianCore()

EshelbianCore::~EshelbianCore ( )
virtualdefault
Examples
EshelbianPlasticity.cpp.

Member Function Documentation

◆ addBoundaryFiniteElement()

MoFEMErrorCode EshelbianCore::addBoundaryFiniteElement ( const EntityHandle  meshset = 0)
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 1442 of file EshelbianPlasticity.cpp.

1442  {
1443  MoFEMFunctionBegin;
1444 
1445  Range meshset_ents;
1446  CHKERR mField.get_moab().get_entities_by_handle(meshset, meshset_ents);
1447 
1448  auto set_fe_adjacency = [&](auto fe_name) {
1449  MoFEMFunctionBegin;
1450  parentAdjSkeletonFunctionDim2 =
1451  boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
1452  bitAdjParent, bitAdjParentMask, bitAdjEnt, bitAdjEntMask);
1453  CHKERR mField.modify_finite_element_adjacency_table(
1454  fe_name, MBTRI, *parentAdjSkeletonFunctionDim2);
1455  MoFEMFunctionReturn(0);
1456  };
1457 
1458  // set finite element fields
1459  auto add_field_to_fe = [this](const std::string fe,
1460  const std::string field_name) {
1461  MoFEMFunctionBegin;
1462  CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
1463  CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
1464  CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
1465  CHKERR mField.modify_finite_element_add_field_data(fe, materialH1Positions);
1466  MoFEMFunctionReturn(0);
1467  };
1468 
1469  if (!mField.check_finite_element(naturalBcElement)) {
1470 
1471  Range natural_bc_elements;
1472  if (bcSpatialDispVecPtr) {
1473  for (auto &v : *bcSpatialDispVecPtr) {
1474  natural_bc_elements.merge(v.faces);
1475  }
1476  }
1477  if (bcSpatialRotationVecPtr) {
1478  for (auto &v : *bcSpatialRotationVecPtr) {
1479  natural_bc_elements.merge(v.faces);
1480  }
1481  }
1482  if (bcSpatialTraction) {
1483  for (auto &v : *bcSpatialTraction) {
1484  natural_bc_elements.merge(v.faces);
1485  }
1486  }
1487  natural_bc_elements = intersect(natural_bc_elements, meshset_ents);
1488 
1489  CHKERR mField.add_finite_element(naturalBcElement, MF_ZERO);
1490  CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
1491  naturalBcElement);
1492  CHKERR add_field_to_fe(naturalBcElement, hybridSpatialDisp);
1493  CHKERR mField.build_finite_elements(naturalBcElement);
1494  }
1495 
1496  auto get_skin = [&](auto &body_ents) {
1497  Skinner skin(&mField.get_moab());
1498  Range skin_ents;
1499  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
1500  return skin_ents;
1501  };
1502 
1503  auto filter_true_skin = [&](auto &&skin) {
1504  Range boundary_ents;
1505  ParallelComm *pcomm =
1506  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1507  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
1508  PSTATUS_NOT, -1, &boundary_ents);
1509  return boundary_ents;
1510  };
1511 
1512  if (!mField.check_finite_element(skinElement)) {
1513 
1514  Range body_ents;
1515  CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM,
1516  body_ents);
1517  auto skin = filter_true_skin(get_skin(body_ents));
1518 
1519  CHKERR mField.add_finite_element(skinElement, MF_ZERO);
1520  CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
1521  CHKERR mField.modify_finite_element_add_field_data(skinElement,
1522  spatialH1Disp);
1523  CHKERR mField.modify_finite_element_add_field_data(skinElement,
1524  materialH1Positions);
1525  CHKERR mField.modify_finite_element_add_field_data(skinElement,
1526  hybridSpatialDisp);
1527  CHKERR mField.modify_finite_element_add_field_data(skinElement,
1528  hybridMaterialDisp);
1529  CHKERR mField.modify_finite_element_add_field_data(skinElement,
1530  contactDisp);
1531  CHKERR mField.build_finite_elements(skinElement);
1532  }
1533 
1534  if (!mField.check_finite_element(contactElement)) {
1535  if (contactFaces) {
1536  MOFEM_LOG("EP", Sev::inform)
1537  << "Contact elements " << contactFaces->size();
1538  CHKERR mField.add_finite_element(contactElement, MF_ZERO);
1539  CHKERR mField.add_ents_to_finite_element_by_type(*contactFaces, MBTRI,
1540  contactElement);
1541  CHKERR add_field_to_fe(contactElement, piolaStress);
1542  CHKERR add_field_to_fe(contactElement, hybridSpatialDisp);
1543  CHKERR add_field_to_fe(contactElement, contactDisp);
1544  CHKERR add_field_to_fe(contactElement, spatialH1Disp);
1545  CHKERR set_fe_adjacency(contactElement);
1546  CHKERR mField.build_finite_elements(contactElement);
1547  }
1548  }
1549 
1550  if (!mField.check_finite_element(skeletonElement)) {
1551  if (!skeletonFaces)
1552  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
1553  MOFEM_LOG("EP", Sev::inform)
1554  << "Skeleton elements " << skeletonFaces->size();
1555  CHKERR mField.add_finite_element(skeletonElement, MF_ZERO);
1556  CHKERR mField.add_ents_to_finite_element_by_type(*skeletonFaces, MBTRI,
1557  skeletonElement);
1558  CHKERR add_field_to_fe(skeletonElement, piolaStress);
1559  CHKERR add_field_to_fe(skeletonElement, hybridSpatialDisp);
1560  CHKERR add_field_to_fe(skeletonElement, contactDisp);
1561  CHKERR add_field_to_fe(skeletonElement, spatialH1Disp);
1562  CHKERR set_fe_adjacency(skeletonElement);
1563  CHKERR mField.build_finite_elements(skeletonElement);
1564  }
1565 
1566  if (!mField.check_finite_element(materialSkeletonElement)) {
1567 
1568  Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM - 2);
1569 
1570  Range front_elements;
1571  for (auto l = 0; l != frontLayers; ++l) {
1572  Range front_elements_layer;
1573  CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM, true,
1574  front_elements_layer,
1575  moab::Interface::UNION);
1576  front_elements.merge(front_elements_layer);
1577  front_edges.clear();
1578  CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
1579  SPACE_DIM - 2, true, front_edges,
1580  moab::Interface::UNION);
1581  }
1582  Range body_ents;
1583  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
1584  Range front_elements_faces;
1585  CHKERR mField.get_moab().get_adjacencies(front_elements, SPACE_DIM - 1,
1586  true, front_elements_faces,
1587  moab::Interface::UNION);
1588  auto body_skin = filter_true_skin(get_skin(body_ents));
1589  auto front_elements_skin = filter_true_skin(get_skin(front_elements));
1590  Range material_skeleton_faces =
1591  subtract(front_elements_faces, front_elements_skin);
1592  material_skeleton_faces.merge(intersect(front_elements_skin, body_skin));
1593 
1594 #ifndef NDEBUG
1595  CHKERR save_range(mField.get_moab(), "front_elements.vtk", front_elements);
1596  CHKERR save_range(mField.get_moab(), "front_elements_skin.vtk",
1597  front_elements_skin);
1598  CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
1599  material_skeleton_faces);
1600 #endif
1601 
1602  CHKERR mField.add_finite_element(materialSkeletonElement, MF_ZERO);
1603  CHKERR mField.add_ents_to_finite_element_by_type(
1604  material_skeleton_faces, MBTRI, materialSkeletonElement);
1605  CHKERR add_field_to_fe(materialSkeletonElement, eshelbyStress);
1606  CHKERR add_field_to_fe(materialSkeletonElement, hybridMaterialDisp);
1607  CHKERR set_fe_adjacency(materialSkeletonElement);
1608  CHKERR mField.build_finite_elements(materialSkeletonElement);
1609  }
1610 
1611  MoFEMFunctionReturn(0);
1612 }

◆ addCrackSurfaces()

MoFEMErrorCode EshelbianCore::addCrackSurfaces ( )
Examples
EshelbianPlasticity.cpp.

Definition at line 5028 of file EshelbianPlasticity.cpp.

5028  {
5029  MoFEMFunctionBegin;
5030 
5031  constexpr bool potential_crack_debug = false;
5032  if constexpr (potential_crack_debug) {
5033 
5034  auto add_ents = get_range_from_block(mField, "POTENTIAL", SPACE_DIM - 1);
5035  Range crack_front_verts;
5036  CHKERR mField.get_moab().get_connectivity(*frontEdges, crack_front_verts,
5037  true);
5038  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
5039  crack_front_verts);
5040  Range crack_front_faces;
5041  CHKERR mField.get_moab().get_adjacencies(crack_front_verts, SPACE_DIM - 1,
5042  true, crack_front_faces,
5043  moab::Interface::UNION);
5044  crack_front_faces = intersect(crack_front_faces, add_ents);
5045  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
5046  crack_front_faces);
5047  CHKERR mField.getInterface<MeshsetsManager>()->addEntitiesToMeshset(
5048  BLOCKSET, addCrackMeshsetId, crack_front_faces);
5049  }
5050 
5051  auto get_crack_faces = [&]() {
5052  if (maxMovedFaces) {
5053  return unite(*crackFaces, *maxMovedFaces);
5054  } else {
5055  return *crackFaces;
5056  }
5057  };
5058 
5059  auto get_extended_crack_faces = [&]() {
5060  auto get_faces_of_crack_front_verts = [&](auto crack_faces_org) {
5061  ParallelComm *pcomm =
5062  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
5063 
5064  Range crack_faces;
5065 
5066  if (!pcomm->rank()) {
5067 
5068  auto get_nodes = [&](auto &&e) {
5069  Range nodes;
5070  CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
5071  "get connectivity");
5072  return nodes;
5073  };
5074 
5075  auto get_adj = [&](auto &&e, auto dim,
5076  auto t = moab::Interface::UNION) {
5077  Range adj;
5078  CHK_MOAB_THROW(
5079  mField.get_moab().get_adjacencies(e, dim, true, adj, t),
5080  "get adj");
5081  return adj;
5082  };
5083 
5084  Range body_ents;
5085  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
5086  body_ents);
5087  auto body_skin = get_skin(mField, body_ents);
5088  auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
5089 
5090  size_t s;
5091  do {
5092  s = crack_faces.size();
5093 
5094  auto crack_face_nodes = get_nodes(crack_faces_org);
5095  auto crack_faces_edges =
5096  get_adj(crack_faces_org, 1, moab::Interface::UNION);
5097 
5098  auto crack_skin =
5099  subtract(get_skin(mField, crack_faces_org), body_skin_edges);
5100  auto crack_skin_nodes = get_nodes(crack_skin);
5101 
5102  auto crack_skin_faces =
5103  get_adj(crack_skin, 2, moab::Interface::UNION);
5104  crack_skin_faces = subtract(crack_skin_faces, crack_faces_org);
5105 
5106  crack_faces = crack_faces_org;
5107  for (auto f : crack_skin_faces) {
5108  auto edges = intersect(
5109  get_adj(Range(f, f), 1, moab::Interface::UNION), crack_skin);
5110  auto edge_conn = get_nodes(Range(edges));
5111  if (edges.size() == 2) {
5112  auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
5113  crack_faces_org);
5114  if (faces.size() == 2) {
5115  auto edge0_conn = get_nodes(Range(edges[0], edges[0]));
5116  auto edge1_conn = get_nodes(Range(edges[1], edges[1]));
5117  auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
5118  crack_skin_nodes);
5119  if (edges_conn.size() == 1) {
5120 
5121  auto node_edges = subtract(intersect(
5122  get_adj(edges_conn, 1, moab::Interface::INTERSECT),
5123  crack_faces_edges), crack_skin);
5124 
5125  if (node_edges.size()) {
5126  FTENSOR_INDEX(SPACE_DIM, i);
5127  FTensor::Tensor1<double, SPACE_DIM> t_v0;
5128  CHKERR mField.get_moab().get_coords(edges_conn, &t_v0(0));
5129 
5130  auto get_t_dir = [&](auto e_conn) {
5131  auto other_node = subtract(e_conn, edges_conn);
5132  FTensor::Tensor1<double, SPACE_DIM> t_dir;
5133  CHKERR mField.get_moab().get_coords(other_node,
5134  &t_dir(0));
5135  t_dir(i) -= t_v0(i);
5136  return t_dir;
5137  };
5138 
5139  FTensor::Tensor1<double, SPACE_DIM> t_ave_dir;
5140  t_ave_dir(i) =
5141  get_t_dir(edge0_conn)(i) + get_t_dir(edge1_conn)(i);
5142 
5143  FTensor::Tensor1<double, SPACE_DIM> t_crack_surface_ave_dir;
5144  t_crack_surface_ave_dir(i) = 0;
5145  for (auto e : node_edges) {
5146  auto e_conn = get_nodes(Range(e, e));
5147  auto t_dir = get_t_dir(e_conn);
5148  t_crack_surface_ave_dir(i) += t_dir(i);
5149  }
5150 
5151  auto dot = t_ave_dir(i) * t_crack_surface_ave_dir(i);
5152  // ave edges is in opposite direction to crack surface, so
5153  // thus crack is not turning back
5154  if (dot < -std::numeric_limits<double>::epsilon()) {
5155  crack_faces.insert(f);
5156  }
5157  }
5158  }
5159  }
5160  } else if (edges.size() == 3) {
5161  crack_faces.insert(f);
5162  }
5163  }
5164 
5165  crack_faces_org = crack_faces;
5166 
5167  } while (s != crack_faces.size());
5168  };
5169 
5170  return send_type(mField, crack_faces, MBTRI);
5171  };
5172 
5173  return get_faces_of_crack_front_verts(get_crack_faces());
5174  };
5175 
5176 #ifndef NDEBUG
5177  constexpr bool debug = false;
5178  if (debug) {
5179  CHKERR save_range(
5180  mField.get_moab(),
5181  "new_crack_surface_" +
5182  boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
5183  get_extended_crack_faces());
5184  }
5185 #endif // NDEBUG
5186 
5187  auto new_crack_faces = subtract(get_extended_crack_faces(), *crackFaces);
5188  CHKERR mField.getInterface<MeshsetsManager>()->addEntitiesToMeshset(
5189  BLOCKSET, addCrackMeshsetId, new_crack_faces);
5190 
5191  MoFEMFunctionReturn(0);
5192 };

◆ addDMs()

MoFEMErrorCode EshelbianCore::addDMs ( const BitRefLevel  bit = BitRefLevel().set(0),
const EntityHandle  meshset = 0 
)
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 1614 of file EshelbianPlasticity.cpp.

1615  {
1616  MoFEMFunctionBegin;
1617 
1618  // find adjacencies between finite elements and dofs
1619  CHKERR mField.build_adjacencies(bit, QUIET);
1620 
1621  // Create coupled problem
1622  dM = createDM(mField.get_comm(), "DMMOFEM");
1623  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
1624  BitRefLevel().set());
1625  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
1626  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1627  CHKERR DMMoFEMAddElement(dM, elementVolumeName);
1628  CHKERR DMMoFEMAddElement(dM, naturalBcElement);
1629  CHKERR DMMoFEMAddElement(dM, skeletonElement);
1630  CHKERR DMMoFEMAddElement(dM, contactElement);
1631  CHKERR DMMoFEMAddElement(dM, skinElement);
1632 
1633  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
1634  CHKERR DMSetUp(dM);
1635  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
1636 
1637  auto remove_dofs_on_broken_skin = [&](const std::string prb_name) {
1638  MoFEMFunctionBegin;
1639  for (int d : {0, 1, 2}) {
1640  std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
1641  CHKERR mField.getInterface<ProblemsManager>()
1642  ->getSideDofsOnBrokenSpaceEntities(
1643  dofs_to_remove, prb_name, ROW, piolaStress,
1644  (*bcSpatialFreeTraction)[d], SPACE_DIM, d, d);
1645  // remove piola dofs, i.e. traction free boundary
1646  CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, ROW,
1647  dofs_to_remove);
1648  CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, COL,
1649  dofs_to_remove);
1650  }
1651  MoFEMFunctionReturn(0);
1652  };
1653  CHKERR remove_dofs_on_broken_skin("ESHELBY_PLASTICITY");
1654 
1655  // Create elastic sub-problem
1656  dmElastic = createDM(mField.get_comm(), "DMMOFEM");
1657  CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
1658  CHKERR DMMoFEMSetDestroyProblem(dmElastic, PETSC_TRUE);
1659  CHKERR DMMoFEMAddSubFieldRow(dmElastic, piolaStress);
1660  CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField);
1661  CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp);
1662  CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis);
1663  if (!noStretch) {
1664  CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor);
1665  }
1666  CHKERR DMMoFEMAddSubFieldRow(dmElastic, hybridSpatialDisp);
1667  CHKERR DMMoFEMAddSubFieldRow(dmElastic, contactDisp);
1668  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
1669  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
1670  CHKERR DMMoFEMAddElement(dmElastic, skeletonElement);
1671  CHKERR DMMoFEMAddElement(dmElastic, contactElement);
1672  CHKERR DMMoFEMAddElement(dmElastic, skinElement);
1673  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
1674  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
1675  CHKERR DMSetUp(dmElastic);
1676 
1677  // dmMaterial = createDM(mField.get_comm(), "DMMOFEM");
1678  // CHKERR DMMoFEMCreateSubDM(dmMaterial, dM, "MATERIAL_PROBLEM");
1679  // CHKERR DMMoFEMSetDestroyProblem(dmMaterial, PETSC_TRUE);
1680  // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, eshelbyStress);
1681  // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, materialL2Disp);
1682  // CHKERR DMMoFEMAddElement(dmMaterh elementVolumeName);
1683  // CHKERR DMMoFEMAddElement(dmMaterial, naturalBcElement);
1684  // CHKERR DMMoFEMAddElement(dmMaterial, skinElement);
1685  // CHKERR DMMoFEMSetSquareProblem(dmMaterial, PETSC_TRUE);
1686  // CHKERR DMMoFEMSetIsPartitioned(dmMaterial, PETSC_TRUE);
1687  // CHKERR DMSetUp(dmMaterial);
1688 
1689  auto set_zero_block = [&]() {
1690  MoFEMFunctionBegin;
1691  if (!noStretch) {
1692  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1693  "ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
1694  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1695  "ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
1696  }
1697  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1698  "ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
1699  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1700  "ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
1701  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1702  "ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
1703  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1704  "ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
1705  if (!noStretch) {
1706  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1707  "ELASTIC_PROBLEM", bubbleField, bubbleField);
1708  CHKERR
1709  mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1710  "ELASTIC_PROBLEM", piolaStress, piolaStress);
1711  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1712  "ELASTIC_PROBLEM", bubbleField, piolaStress);
1713  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1714  "ELASTIC_PROBLEM", piolaStress, bubbleField);
1715  }
1716 
1717  solTSStep = createDMVector(dmElastic);
1718  MoFEMFunctionReturn(0);
1719  };
1720 
1721  auto set_section = [&]() {
1722  MoFEMFunctionBegin;
1723  PetscSection section;
1724  CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
1725  &section);
1726  CHKERR DMSetSection(dmElastic, section);
1727  CHKERR DMSetGlobalSection(dmElastic, section);
1728  CHKERR PetscSectionDestroy(&section);
1729  MoFEMFunctionReturn(0);
1730  };
1731 
1732  CHKERR set_zero_block();
1733  CHKERR set_section();
1734 
1735  dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
1736  CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
1737  CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
1738  CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp);
1739  CHKERR DMMoFEMAddElement(dmPrjSpatial, elementVolumeName);
1740  CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
1741  CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
1742  CHKERR DMSetUp(dmPrjSpatial);
1743 
1744  CHKERR mField.getInterface<BcManager>()
1745  ->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
1746  "PROJECT_SPATIAL", spatialH1Disp, true, false);
1747 
1748  MoFEMFunctionReturn(0);
1749 }

◆ addFields()

MoFEMErrorCode EshelbianCore::addFields ( const EntityHandle  meshset = 0)
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 1027 of file EshelbianPlasticity.cpp.

1027  {
1028  MoFEMFunctionBegin;
1029 
1030  auto get_tets = [&]() {
1031  Range tets;
1032  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
1033  return tets;
1034  };
1035 
1036  auto get_tets_skin = [&]() {
1037  Range tets_skin_part;
1038  Skinner skin(&mField.get_moab());
1039  CHKERR skin.find_skin(0, get_tets(), false, tets_skin_part);
1040  ParallelComm *pcomm =
1041  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1042  Range tets_skin;
1043  CHKERR pcomm->filter_pstatus(tets_skin_part,
1044  PSTATUS_SHARED | PSTATUS_MULTISHARED,
1045  PSTATUS_NOT, -1, &tets_skin);
1046  return tets_skin;
1047  };
1048 
1049  auto subtract_boundary_conditions = [&](auto &&tets_skin) {
1050  // That mean, that hybrid field on all faces on which traction is applied,
1051  // on other faces, or enforcing displacements as
1052  // natural boundary condition.
1053  if (bcSpatialTraction)
1054  for (auto &v : *bcSpatialTraction) {
1055  tets_skin = subtract(tets_skin, v.faces);
1056  }
1057  return tets_skin;
1058  };
1059 
1060  auto add_blockset = [&](auto block_name, auto &&tets_skin) {
1061  auto crack_faces =
1062  get_range_from_block(mField, "block_name", SPACE_DIM - 1);
1063  tets_skin.merge(crack_faces);
1064  return tets_skin;
1065  };
1066 
1067  auto subtract_blockset = [&](auto block_name, auto &&tets_skin) {
1068  auto contact_range =
1069  get_range_from_block(mField, block_name, SPACE_DIM - 1);
1070  tets_skin = subtract(tets_skin, contact_range);
1071  return tets_skin;
1072  };
1073 
1074  auto get_stress_trace_faces = [&](auto &&tets_skin) {
1075  Range faces;
1076  CHKERR mField.get_moab().get_adjacencies(get_tets(), SPACE_DIM - 1, true,
1077  faces, moab::Interface::UNION);
1078  Range trace_faces = subtract(faces, tets_skin);
1079  return trace_faces;
1080  };
1081 
1082  auto tets = get_tets();
1083 
1084  // remove also contact faces, i.e. that is also kind of hybrid field but
1085  // named but used to enforce contact conditions
1086  auto trace_faces = get_stress_trace_faces(
1087 
1088  subtract_blockset("CONTACT",
1089  subtract_boundary_conditions(get_tets_skin()))
1090 
1091  );
1092 
1093  contactFaces = boost::make_shared<Range>(intersect(
1094  trace_faces, get_range_from_block(mField, "CONTACT", SPACE_DIM - 1)));
1095  skeletonFaces =
1096  boost::make_shared<Range>(subtract(trace_faces, *contactFaces));
1097  materialSkeletonFaces =
1098  boost::make_shared<Range>(get_stress_trace_faces(Range()));
1099 
1100 #ifndef NDEBUG
1101  if (contactFaces->size())
1102  CHKERR save_range(mField.get_moab(), "contact_faces.vtk", *contactFaces);
1103  if (skeletonFaces->size())
1104  CHKERR save_range(mField.get_moab(), "skeleton_faces.vtk", *skeletonFaces);
1105  if (skeletonFaces->size())
1106  CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
1107  *materialSkeletonFaces);
1108 #endif
1109 
1110  auto add_broken_hdiv_field = [this, meshset](const std::string field_name,
1111  const int order) {
1112  MoFEMFunctionBegin;
1113 
1114  constexpr FieldApproximationBase base = DEMKOWICZ_JACOBI_BASE;
1115 
1116  auto get_side_map_hdiv = [&]() {
1117  return std::vector<
1118 
1119  std::pair<EntityType,
1120  std::function<MoFEMErrorCode(BaseFunction::DofsSideMap &)>
1121 
1122  >>{
1123 
1124  {MBTET,
1125  [&](BaseFunction::DofsSideMap &dofs_side_map) -> MoFEMErrorCode {
1126  return TetPolynomialBase::setDofsSideMap(HDIV, DISCONTINUOUS, base,
1127  dofs_side_map);
1128  }}
1129 
1130  };
1131  };
1132 
1133  CHKERR mField.add_broken_field(field_name, HDIV, base, SPACE_DIM,
1134  get_side_map_hdiv(), MB_TAG_DENSE, MF_ZERO);
1135  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
1136  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
1137  MoFEMFunctionReturn(0);
1138  };
1139 
1140  auto add_l2_field = [this, meshset](const std::string field_name,
1141  const int order, const int dim) {
1142  MoFEMFunctionBegin;
1143  CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
1144  MB_TAG_DENSE, MF_ZERO);
1145  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
1146  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
1147  MoFEMFunctionReturn(0);
1148  };
1149 
1150  auto add_h1_field = [this, meshset](const std::string field_name,
1151  const int order, const int dim) {
1152  MoFEMFunctionBegin;
1153  CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
1154  MB_TAG_DENSE, MF_ZERO);
1155  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
1156  CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
1157  CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
1158  CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
1159  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
1160  MoFEMFunctionReturn(0);
1161  };
1162 
1163  auto add_l2_field_by_range = [this](const std::string field_name,
1164  const int order, const int dim,
1165  const int field_dim, Range &&r) {
1166  MoFEMFunctionBegin;
1167  CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, field_dim,
1168  MB_TAG_DENSE, MF_ZERO);
1169  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(r);
1170  CHKERR mField.add_ents_to_field_by_dim(r, dim, field_name);
1171  CHKERR mField.set_field_order(r, field_name, order);
1172  MoFEMFunctionReturn(0);
1173  };
1174 
1175  auto add_bubble_field = [this, meshset](const std::string field_name,
1176  const int order, const int dim) {
1177  MoFEMFunctionBegin;
1178  CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
1179  MF_ZERO);
1180  // Modify field
1181  auto field_ptr = mField.get_field_structure(field_name);
1182  auto field_order_table =
1183  const_cast<Field *>(field_ptr)->getFieldOrderTable();
1184  auto get_cgg_bubble_order_zero = [](int p) { return 0; };
1185  auto get_cgg_bubble_order_tet = [](int p) {
1186  return NBVOLUMETET_CCG_BUBBLE(p);
1187  };
1188  field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
1189  field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
1190  field_order_table[MBTRI] = get_cgg_bubble_order_zero;
1191  field_order_table[MBTET] = get_cgg_bubble_order_tet;
1192  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
1193  CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
1194  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
1195  MoFEMFunctionReturn(0);
1196  };
1197 
1198  auto add_user_l2_field = [this, meshset](const std::string field_name,
1199  const int order, const int dim) {
1200  MoFEMFunctionBegin;
1201  CHKERR mField.add_field(field_name, L2, USER_BASE, dim, MB_TAG_DENSE,
1202  MF_ZERO);
1203  // Modify field
1204  auto field_ptr = mField.get_field_structure(field_name);
1205  auto field_order_table =
1206  const_cast<Field *>(field_ptr)->getFieldOrderTable();
1207  auto zero_dofs = [](int p) { return 0; };
1208  auto dof_l2_tet = [](int p) { return NBVOLUMETET_L2(p); };
1209  field_order_table[MBVERTEX] = zero_dofs;
1210  field_order_table[MBEDGE] = zero_dofs;
1211  field_order_table[MBTRI] = zero_dofs;
1212  field_order_table[MBTET] = dof_l2_tet;
1213  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
1214  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
1215  MoFEMFunctionReturn(0);
1216  };
1217 
1218  // spatial fields
1219  CHKERR add_broken_hdiv_field(piolaStress, spaceOrder);
1220  CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
1221  CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
1222  CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
1223  CHKERR add_user_l2_field(stretchTensor, noStretch ? -1 : spaceOrder, 6);
1224 
1225  if (!skeletonFaces)
1226  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
1227  if (!contactFaces)
1228  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No contact faces");
1229 
1230  CHKERR add_l2_field_by_range(hybridSpatialDisp, spaceOrder - 1, 2, 3,
1231  Range(*skeletonFaces));
1232  CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
1233  Range(*contactFaces));
1234 
1235  // spatial displacement
1236  CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
1237  // material positions
1238  CHKERR add_h1_field(materialH1Positions, 2, 3);
1239 
1240  // Eshelby stress
1241  CHKERR add_broken_hdiv_field(eshelbyStress, spaceOrder);
1242  CHKERR add_l2_field(materialL2Disp, spaceOrder - 1, 3);
1243  CHKERR add_l2_field_by_range(hybridMaterialDisp, spaceOrder - 1, 2, 3,
1244  Range(*materialSkeletonFaces));
1245 
1246  CHKERR mField.build_fields();
1247 
1248  MoFEMFunctionReturn(0);
1249 }

◆ addMaterial_Hencky()

MoFEMErrorCode EshelbianCore::addMaterial_Hencky ( double  E,
double  nu 
)
Examples
ep.cpp.

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

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

◆ addMaterial_HMHHStVenantKirchhoff()

MoFEMErrorCode EshelbianCore::addMaterial_HMHHStVenantKirchhoff ( const int  tape,
const double  lambda,
const double  mu,
const double  sigma_y 
)
Examples
ep.cpp.

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

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

◆ addMaterial_HMHMooneyRivlin()

MoFEMErrorCode EshelbianCore::addMaterial_HMHMooneyRivlin ( const int  tape,
const double  alpha,
const double  beta,
const double  lambda,
const double  sigma_y 
)
Examples
ep.cpp.

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

509  {
510  MoFEMFunctionBegin;
511  physicalEquations =
512  boost::make_shared<HMHPMooneyRivlinWriggersEq63>(alpha, beta, lambda);
513  CHKERR physicalEquations->recordTape(tape, nullptr);
514  MoFEMFunctionReturn(0);
515 }

◆ addMaterial_HMHNeohookean()

MoFEMErrorCode EshelbianCore::addMaterial_HMHNeohookean ( const int  tape,
const double  c10,
const double  K 
)
Examples
ep.cpp.

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

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

◆ addVolumeFiniteElement()

MoFEMErrorCode EshelbianCore::addVolumeFiniteElement ( const EntityHandle  meshset = 0)
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 1380 of file EshelbianPlasticity.cpp.

1380  {
1381  MoFEMFunctionBegin;
1382 
1383  // set finite element fields
1384  auto add_field_to_fe = [this](const std::string fe,
1385  const std::string field_name) {
1386  MoFEMFunctionBegin;
1387  CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
1388  CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
1389  CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
1390  MoFEMFunctionReturn(0);
1391  };
1392 
1393  if (!mField.check_finite_element(elementVolumeName)) {
1394  CHKERR mField.add_finite_element(elementVolumeName, MF_ZERO);
1395  CHKERR mField.add_ents_to_finite_element_by_type(meshset, MBTET,
1396  elementVolumeName);
1397 
1398  CHKERR add_field_to_fe(elementVolumeName, piolaStress);
1399  CHKERR add_field_to_fe(elementVolumeName, bubbleField);
1400  if (!noStretch)
1401  CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
1402  CHKERR add_field_to_fe(elementVolumeName, rotAxis);
1403  CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
1404  CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
1405  CHKERR add_field_to_fe(elementVolumeName, contactDisp);
1406  CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
1407  materialH1Positions);
1408 
1409  // build finite elements data structures
1410  CHKERR mField.build_finite_elements(elementVolumeName);
1411  }
1412 
1413  if (!mField.check_finite_element(materialVolumeElement)) {
1414 
1415  Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM - 2);
1416 
1417  Range front_elements;
1418  for (auto l = 0; l != frontLayers; ++l) {
1419  Range front_elements_layer;
1420  CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM, true,
1421  front_elements_layer,
1422  moab::Interface::UNION);
1423  front_elements.merge(front_elements_layer);
1424  front_edges.clear();
1425  CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
1426  SPACE_DIM - 2, true, front_edges,
1427  moab::Interface::UNION);
1428  }
1429 
1430  CHKERR mField.add_finite_element(materialVolumeElement, MF_ZERO);
1431  CHKERR mField.add_ents_to_finite_element_by_type(front_elements, MBTET,
1432  materialVolumeElement);
1433  CHKERR add_field_to_fe(materialVolumeElement, eshelbyStress);
1434  CHKERR add_field_to_fe(materialVolumeElement, materialL2Disp);
1435  CHKERR mField.build_finite_elements(materialVolumeElement);
1436  }
1437 
1438  MoFEMFunctionReturn(0);
1439 }

◆ calculateEnergyRelease()

MoFEMErrorCode EshelbianCore::calculateEnergyRelease ( const int  tag,
TS  ts 
)
Examples
EshelbianPlasticity.cpp.

Definition at line 3638 of file EshelbianPlasticity.cpp.

3638  {
3639  MoFEMFunctionBegin;
3640  MOFEM_LOG_CHANNEL("SELF");
3641 
3642  constexpr bool debug = true;
3643 
3644  SmartPetscObj<IS> streach_is;
3645  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
3646  "ELASTIC_PROBLEM", ROW, stretchTensor, 0, 6, streach_is);
3647  SmartPetscObj<IS> piola_is;
3648  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
3649  "ELASTIC_PROBLEM", ROW, piolaStress, 0, SPACE_DIM, piola_is);
3650  SmartPetscObj<IS> bubble_is;
3651  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
3652  "ELASTIC_PROBLEM", ROW, bubbleField, 0, SPACE_DIM, bubble_is);
3653  SmartPetscObj<IS> rot_axis_is;
3654  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
3655  "ELASTIC_PROBLEM", ROW, rotAxis, 0, SPACE_DIM, rot_axis_is);
3656  SmartPetscObj<IS> spatial_l2_disp_is;
3657  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
3658  "ELASTIC_PROBLEM", ROW, spatialL2Disp, 0, SPACE_DIM, spatial_l2_disp_is);
3659  SmartPetscObj<IS> hybrid_is;
3660  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
3661  "ELASTIC_PROBLEM", ROW, hybridSpatialDisp, 0, SPACE_DIM, hybrid_is);
3662 
3663  auto get_tags_vec = [&]() {
3664  std::vector<Tag> tags(1);
3665 
3666  auto create_and_clean = [&]() {
3667  MoFEMFunctionBegin;
3668  auto &moab = mField.get_moab();
3669  auto rval = moab.tag_get_handle("ReleaseEnergy", tags[0]);
3670  if (rval == MB_SUCCESS) {
3671  moab.tag_delete(tags[0]);
3672  }
3673  double def_val[] = {0};
3674  CHKERR moab.tag_get_handle("ReleaseEnergy", 1, MB_TYPE_DOUBLE, tags[0],
3675  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
3676  MoFEMFunctionReturn(0);
3677  };
3678 
3679  CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
3680 
3681  return tags;
3682  };
3683 
3684  auto tags = get_tags_vec();
3685 
3686  auto get_adj_front = [&]() {
3687  auto skeleton_faces = *skeletonFaces;
3688  Range adj_front;
3689  CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2, true, adj_front,
3690  moab::Interface::UNION);
3691 
3692  adj_front = intersect(adj_front, skeleton_faces);
3693  adj_front = subtract(adj_front, *crackFaces);
3694  return adj_front;
3695  };
3696 
3697  auto get_nb_front_faces_on_proc = [&](auto &&adj_front) {
3698  int send_size = adj_front.size();
3699  std::vector<int> recv_data(mField.get_comm_size());
3700  MPI_Allgather(&send_size, 1, MPI_INT, recv_data.data(), 1, MPI_INT,
3701  mField.get_comm());
3702  recv_data[mField.get_comm_rank()] = send_size;
3703  return std::pair(recv_data, adj_front);
3704  };
3705 
3706  auto get_snes = [&](TS ts) {
3707  SNES snes;
3708  CHKERR TSGetSNES(ts, &snes);
3709  return snes;
3710  };
3711 
3712  auto get_ksp = [&](SNES snes) {
3713  KSP ksp;
3714  CHKERR SNESGetKSP(snes, &ksp);
3715  CHKERR KSPSetFromOptions(ksp);
3716  return ksp;
3717  };
3718 
3719  auto integrate_energy = [&](auto x, auto release_energy_ptr) {
3720  MoFEMFunctionBegin;
3721  auto set_volume = [&]() {
3722  auto fe_ptr =
3723  boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
3724  fe_ptr->ts_ctx = TSMethod::CTX_TSSETIJACOBIAN;
3725  CHKERR setBaseVolumeElementOps(tag, false, false, false, x, fe_ptr);
3726  fe_ptr->getOpPtrVector().push_back(
3727  new OpReleaseEnergy(dataAtPts, release_energy_ptr));
3728  return fe_ptr;
3729  };
3730  CHKERR DMoFEMLoopFiniteElements(dM, elementVolumeName, set_volume());
3731  MoFEMFunctionReturn(0);
3732  };
3733 
3734  auto calc_release_energy = [&](auto t, auto &&tuple_of_vectors,
3735  auto adj_front, double &energy, double &area) {
3736  MoFEMFunctionBegin;
3737  MOFEM_LOG("EP", Sev::inform) << "Calculate face release energy";
3738 
3739 
3740  auto copy_is = [&](auto is, auto a, auto b) {
3741  MoFEMFunctionBegin;
3742  const PetscInt *is_array;
3743  PetscInt is_size;
3744  CHKERR ISGetLocalSize(is, &is_size);
3745  CHKERR ISGetIndices(is, &is_array);
3746  double *pa;
3747  CHKERR VecGetArray(a, &pa);
3748  const double *pb;
3749  CHKERR VecGetArrayRead(b, &pb);
3750  for (int i = 0; i != is_size; ++i) {
3751  pa[is_array[i]] = pb[is_array[i]];
3752  }
3753  CHKERR VecRestoreArrayRead(b, &pb);
3754  CHKERR VecRestoreArray(a, &pa);
3755  CHKERR ISRestoreIndices(is, &is_array);
3756  MoFEMFunctionReturn(0);
3757  };
3758 
3759  auto axpy_is = [&](auto is, double alpha, auto a, auto b) {
3760  MoFEMFunctionBegin;
3761  const PetscInt *is_array;
3762  PetscInt is_size;
3763  CHKERR ISGetLocalSize(is, &is_size);
3764  CHKERR ISGetIndices(is, &is_array);
3765  double *pa;
3766  CHKERR VecGetArray(a, &pa);
3767  const double *pb;
3768  CHKERR VecGetArrayRead(b, &pb);
3769  for (int i = 0; i != is_size; ++i) {
3770  pa[is_array[i]] += alpha * pb[is_array[i]];
3771  }
3772  CHKERR VecRestoreArrayRead(b, &pb);
3773  CHKERR VecRestoreArray(a, &pa);
3774  CHKERR ISRestoreIndices(is, &is_array);
3775  MoFEMFunctionReturn(0);
3776  };
3777 
3778  auto zero_is = [&](auto is, auto a) {
3779  MoFEMFunctionBegin;
3780  const PetscInt *is_array;
3781  PetscInt is_size;
3782  CHKERR ISGetLocalSize(is, &is_size);
3783  CHKERR ISGetIndices(is, &is_array);
3784  double *pa;
3785  CHKERR VecGetArray(a, &pa);
3786  for (int i = 0; i != is_size; ++i) {
3787  pa[is_array[i]] = 0;
3788  }
3789  CHKERR VecRestoreArray(a, &pa);
3790  CHKERR ISRestoreIndices(is, &is_array);
3791  MoFEMFunctionReturn(0);
3792  };
3793 
3794  auto estimate_energy = [&](auto x) {
3795  auto release_energy_ptr = boost::make_shared<double>(0);
3796  CHK_THROW_MESSAGE(integrate_energy(x, release_energy_ptr),
3797  "integrate_energy");
3798 
3799  double area = 0;
3800  for (auto f : filter_owners(mField, adj_front)) {
3801  FTensor::Tensor1<double, SPACE_DIM> t_normal;
3802  CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
3803  area += t_normal.l2() / 2;
3804  }
3805 
3806  std::array<double, 2> send_data{area, *release_energy_ptr};
3807  std::array<double, 2> recv_data{0, 0};
3808 
3809  MPI_Allreduce(send_data.data(), recv_data.data(), 2, MPI_DOUBLE, MPI_SUM,
3810  mField.get_comm());
3811  return std::pair(recv_data[1], recv_data[0]);
3812  };
3813 
3814  auto set_tag = [&](auto &&release_energy, auto adj_front) {
3815  MoFEMFunctionBegin;
3816  auto own = filter_owners(mField, adj_front);
3817  MOFEM_LOG("EP", Sev::inform) << "Energy release " << release_energy;
3818  auto &moab = mField.get_moab();
3819  CHKERR moab.tag_clear_data(tags[0], own, &release_energy);
3820  MoFEMFunctionReturn(0);
3821  };
3822 
3823  auto solve = [&]() {
3824  MoFEMFunctionBegin;
3825 
3826  auto [x, f] = tuple_of_vectors;
3827 
3828  Mat mA;
3829  auto ksp = get_ksp(get_snes(ts));
3830  CHKERR KSPGetOperators(ksp, &mA, PETSC_NULL);
3831 
3832  Range faces = unite(adj_front, *crackFaces);
3833  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(faces);
3834 
3835  std::vector<boost::weak_ptr<NumeredDofEntity>> piola_skelton_dofs_vec;
3836  CHKERR mField.getInterface<ProblemsManager>()
3837  ->getSideDofsOnBrokenSpaceEntities(
3838  piola_skelton_dofs_vec, "ELASTIC_PROBLEM", ROW, piolaStress,
3839  faces, SPACE_DIM, 0, SPACE_DIM);
3840  SmartPetscObj<IS> skeleton_piola_is_local;
3841  CHKERR mField.getInterface<ISManager>()
3842  ->isCreateProblemBrokenFieldAndRankLocal(piola_skelton_dofs_vec,
3843  skeleton_piola_is_local);
3844  SmartPetscObj<IS> skeleton_piola_is_global;
3845  CHKERR mField.getInterface<ISManager>()
3846  ->isCreateProblemBrokenFieldAndRank(piola_skelton_dofs_vec,
3847  skeleton_piola_is_global);
3848 
3849  SmartPetscObj<IS> skeleton_hybrid_is_local;
3850  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRankLocal(
3851  "ELASTIC_PROBLEM", ROW, hybridSpatialDisp, 0, SPACE_DIM,
3852  skeleton_hybrid_is_local, &faces);
3853  SmartPetscObj<IS> skeleton_hybrid_is_global;
3854  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
3855  "ELASTIC_PROBLEM", ROW, hybridSpatialDisp, 0, SPACE_DIM,
3856  skeleton_hybrid_is_global, &faces);
3857 
3858  CHKERR VecZeroEntries(x);
3859  CHKERR copy_is(skeleton_piola_is_local, x, t);
3860  CHKERR MatMult(mA, x, f); // f = A*x
3861  CHKERR zero_is(skeleton_piola_is_local, f);
3862  CHKERR zero_is(skeleton_hybrid_is_local, f);
3863 
3864  auto mat_storage = getBlockMatStorageMat(mA);
3865  auto mat_precon_storage = getBlockMatPrconditionerStorageMat(mA);
3866  std::vector<double> save_mat_storage(*mat_storage);
3867  std::vector<double> save_mat_precon_storage(*mat_precon_storage);
3868 
3869  CHKERR MatZeroRowsColumnsIS(mA, skeleton_piola_is_global, 1, PETSC_NULL,
3870  PETSC_NULL);
3871  CHKERR MatZeroRowsColumnsIS(mA, skeleton_hybrid_is_global, 1, PETSC_NULL,
3872  PETSC_NULL);
3873 
3874  PetscBool factor_schur = PETSC_FALSE;
3875  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-crack_factor_schur",
3876  &factor_schur, PETSC_NULL);
3877  if (factor_schur == PETSC_TRUE) {
3878  SmartPetscObj<IS> schur_skeleton_hybrid_is;
3879  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
3880  "SUB_SCHUR", ROW, hybridSpatialDisp, 0, SPACE_DIM,
3881  schur_skeleton_hybrid_is, &faces);
3882  CHKERR MatZeroEntries(S);
3883  CHKERR assembleBlockMatSchur(mField, mA, S, a00FieldList, a00RangeList,
3884  SmartPetscObj<AO>(aoS, true));
3885  // Apply essential constrains to Schur complement
3886  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
3887  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
3888  CHKERR MatZeroRowsColumnsIS(S, schur_skeleton_hybrid_is, 1, PETSC_NULL,
3889  PETSC_NULL);
3890  }
3891 
3892  CHKERR copy_is(skeleton_piola_is_local, f, t);
3893  CHKERR VecScale(f, 1. / griffithEnergy);
3894  CHKERR KSPSolve(ksp, f, x);
3895 
3896  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
3897  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
3898 
3899  std::copy(save_mat_storage.begin(), save_mat_storage.end(),
3900  mat_storage->begin());
3901  std::copy(save_mat_precon_storage.begin(), save_mat_precon_storage.end(),
3902  mat_precon_storage->begin());
3903 
3904  MoFEMFunctionReturn(0);
3905  };
3906 
3907  CHKERR solve();
3908  auto [x, f] = tuple_of_vectors;
3909  std::tie(energy, area) = estimate_energy(x);
3910  CHKERR set_tag(energy, adj_front);
3911 
3912  MoFEMFunctionReturn(0);
3913  };
3914 
3915  auto run_faces = [&](auto &&p) {
3916  auto [recv_data, adj_front] = p;
3917  Vec t;
3918  CHKERR TSGetSolution(ts, &t);
3919  auto face_x = vectorDuplicate(t);
3920  auto face_f = vectorDuplicate(t);
3921 
3922  std::map<int, std::tuple<double, double, Range, SmartPetscObj<Vec>>>
3923  release_by_energy;
3924 
3925  auto solve_and_add = [&](auto &&face, auto id) {
3926  MoFEMFunctionBegin;
3927  MOFEM_LOG("SYNC", Sev::noisy) << face;
3928  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
3929  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(face);
3930  double energy, area;
3931  CHKERR calc_release_energy(t,
3932  std::make_tuple(face_x, face_f),
3933  face, energy, area);
3934 
3935  if (release_by_energy.find(id) == release_by_energy.end()) {
3936  auto vec = vectorDuplicate(face_x);
3937  CHKERR VecCopy(face_x, vec);
3938  release_by_energy[id] = std::make_tuple(energy, area, face, vec);
3939  } else {
3940  auto [e, a, faces, vec] = release_by_energy.at(id);
3941  if (e < energy) {
3942  CHKERR VecCopy(face_x, vec);
3943  faces.merge(face);
3944  release_by_energy[id] = std::make_tuple(energy, area, faces, vec);
3945  }
3946  }
3947 
3948  MoFEMFunctionReturn(0);
3949  };
3950 
3951  ParallelComm *pcomm =
3952  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
3953 
3954  for (auto p = 0; p < mField.get_comm_rank(); ++p) {
3955  for (auto f = 0; f < recv_data[p]; ++f) {
3956  int id;
3957  MPI_Bcast(&id, 1, MPI_INT, p, mField.get_comm());
3958  MOFEM_LOG_CHANNEL("SYNC");
3959  MOFEM_LOG("SYNC", Sev::noisy) << "edge id " << id;
3960  solve_and_add(Range(), id);
3961  }
3962  }
3963  for (auto f = 0; f < recv_data[mField.get_comm_rank()]; ++f) {
3964  auto face = Range(adj_front[f], adj_front[f]);
3965  Range edges;
3966  CHK_MOAB_THROW(mField.get_moab().get_adjacencies(face, 1, true, edges,
3967  moab::Interface::UNION),
3968  "get adj");
3969  edges = intersect(edges, *frontEdges);
3970  int owner;
3971  EntityHandle owner_handle;
3972  CHK_MOAB_THROW(pcomm->get_owner_handle(edges[0], owner, owner_handle),
3973  "get handle");
3974  int id = id_from_handle(owner_handle);
3975  MPI_Bcast(&id, 1, MPI_INT, mField.get_comm_rank(), mField.get_comm());
3976  MOFEM_LOG_CHANNEL("SYNC");
3977  MOFEM_LOG("SYNC", Sev::noisy) << "owner " << owner << "edge id " << id;
3978  solve_and_add(face, id);
3979  }
3980  for (auto p = mField.get_comm_rank() + 1; p < mField.get_comm_size(); ++p) {
3981  for (auto f = 0; f < recv_data[p]; ++f) {
3982  int id;
3983  MPI_Bcast(&id, 1, MPI_INT, p, mField.get_comm());
3984  MOFEM_LOG_CHANNEL("SYNC");
3985  MOFEM_LOG("SYNC", Sev::noisy) << "edge id " << id;
3986  solve_and_add(Range(), id);
3987  }
3988  }
3989  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
3990 
3991  return release_by_energy;
3992  };
3993 
3994  auto release_by_energy = run_faces(
3995 
3996  get_nb_front_faces_on_proc(filter_owners(mField, get_adj_front()))
3997 
3998  );
3999 
4000  std::map<double, std::tuple<double, EntityHandle, Range, SmartPetscObj<Vec>>>
4001  sorted_release_by_energy;
4002 
4003  for (auto &m : release_by_energy) {
4004  auto [energy, area, faces, vec] = m.second;
4005  sorted_release_by_energy[energy] =
4006  std::make_tuple(area, m.first, faces, vec);
4007  }
4008 
4009  double total_area = 0;
4010  double total_energy_sum = 0;
4011  Vec t;
4012  CHKERR TSGetSolution(ts, &t);
4013  auto x = vectorDuplicate(t);
4014  CHKERR VecZeroEntries(x);
4015  for (auto m = sorted_release_by_energy.rbegin();
4016  m != sorted_release_by_energy.rend(); ++m) {
4017  auto [area, id, faces, vec] = m->second;
4018  auto edge = ent_form_type_and_id(MBEDGE, id);
4019  CHKERR VecAXPY(x, 1., vec);
4020  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
4021  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
4022  auto release_energy_ptr = boost::make_shared<double>(0);
4023  CHKERR integrate_energy(x, release_energy_ptr);
4024  std::array<double, 2> send_data{area, *release_energy_ptr};
4025  std::array<double, 2> recv_data{0, 0};
4026  MPI_Allreduce(send_data.data(), recv_data.data(), 2, MPI_DOUBLE, MPI_SUM,
4027  mField.get_comm());
4028  if (recv_data[1] < total_energy_sum) {
4029  CHKERR VecAXPY(x, -1., vec);
4030  double release_energy = 0;
4031  CHKERR mField.get_moab().tag_clear_data(tags[0], faces, &release_energy);
4032  } else {
4033  total_energy_sum = recv_data[1];
4034  total_area += recv_data[0];
4035  }
4036  MOFEM_LOG("EP", Sev::inform)
4037  << "Aggregated energy release " << total_area << " " << recv_data[1];
4038  }
4039 
4040  CHKERR CommInterface::updateEntitiesPetscVector(mField.get_moab(),
4041  volumeExchange, tags[0]);
4042  CHKERR CommInterface::updateEntitiesPetscVector(mField.get_moab(),
4043  faceExchange, tags[0]);
4044 
4045  if (debug) {
4046  Range body_ents;
4047  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
4048  auto skin = filter_true_skin(mField, get_skin(mField, body_ents));
4049  for (auto f : skin) {
4050  Range tet;
4051  CHKERR mField.get_moab().get_adjacencies(&f, 1, SPACE_DIM, false, tet);
4052  if (tet.size() != 1) {
4053  MOFEM_LOG("EP", Sev::error) << "tet.size()!=1";
4054  continue;
4055  }
4056  double release_energy;
4057  CHKERR mField.get_moab().tag_get_data(tags[0], tet, &release_energy);
4058  release_energy /= total_energy_sum;
4059  CHKERR mField.get_moab().tag_set_data(tags[0], &f, 1, &release_energy);
4060  }
4061 
4062  CHKERR postProcessResults(1, "test_perturbation.h5m", PETSC_NULL, x, tags);
4063  }
4064 
4065  MoFEMFunctionReturn(0);
4066 }

◆ 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. The is then maximises across all processors.

For each front edge, find maximal face energy and orientation. This is done by solving quadratic equation.

Examples
EshelbianPlasticity.cpp.

Definition at line 4068 of file EshelbianPlasticity.cpp.

4069  {
4070  MoFEMFunctionBegin;
4071 
4072  constexpr bool debug = true;
4073  constexpr auto sev = Sev::inform;
4074 
4075  auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
4076  Range geometry_edges_verts;
4077  CHKERR mField.get_moab().get_connectivity(geometry_edges,
4078  geometry_edges_verts, true);
4079  Range crack_faces_verts;
4080  CHKERR mField.get_moab().get_connectivity(*crackFaces, crack_faces_verts,
4081  true);
4082  Range crack_faces_edges;
4083  CHKERR mField.get_moab().get_adjacencies(
4084  *crackFaces, 1, true, crack_faces_edges, moab::Interface::UNION);
4085 
4086  auto get_tags_vec = [&](auto tag_name, int dim) {
4087  std::vector<Tag> tags(1);
4088 
4089  if (dim > 3)
4090  CHK_THROW_MESSAGE(MOFEM_ATOM_TEST_INVALID, "dim>3");
4091 
4092  auto create_and_clean = [&]() {
4093  MoFEMFunctionBegin;
4094  auto &moab = mField.get_moab();
4095  auto rval = moab.tag_get_handle(tag_name, tags[0]);
4096  if (rval == MB_SUCCESS) {
4097  moab.tag_delete(tags[0]);
4098  }
4099  double def_val[] = {0., 0., 0.};
4100  CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
4101  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
4102  MoFEMFunctionReturn(0);
4103  };
4104 
4105  CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
4106 
4107  return tags;
4108  };
4109 
4110  auto get_adj_front = [&](bool subtract_crack) {
4111  Range adj_front;
4112  CHKERR mField.get_moab().get_adjacencies(*frontEdges, SPACE_DIM - 1, true,
4113  adj_front, moab::Interface::UNION);
4114  if (subtract_crack)
4115  adj_front = subtract(adj_front, *crackFaces);
4116  return adj_front;
4117  };
4118 
4119  MOFEM_LOG_CHANNEL("SELF");
4120 
4121  auto th_front_position = get_tags_vec("FrontPosition", 3);
4122  auto th_max_face_energy = get_tags_vec("MaxFaceEnergy", 1);
4123 
4124  if (mField.get_comm_rank() == 0) {
4125 
4126  auto get_crack_adj_tets = [&](auto r) {
4127  Range crack_faces_conn;
4128  CHKERR mField.get_moab().get_connectivity(r, crack_faces_conn);
4129  Range crack_faces_conn_tets;
4130  CHKERR mField.get_moab().get_adjacencies(crack_faces_conn, SPACE_DIM,
4131  true, crack_faces_conn_tets,
4132  moab::Interface::UNION);
4133  return crack_faces_conn_tets;
4134  };
4135 
4136  auto get_layers_for_sides = [&](auto &side) {
4137  std::vector<Range> layers;
4138  auto get = [&]() {
4139  MoFEMFunctionBegin;
4140 
4141  auto get_adj = [&](auto &r, int dim) {
4142  Range adj;
4143  CHKERR mField.get_moab().get_adjacencies(r, dim, true, adj,
4144  moab::Interface::UNION);
4145  return adj;
4146  };
4147 
4148  auto get_tets = [&](auto r) { return get_adj(r, SPACE_DIM); };
4149 
4150  Range front_nodes;
4151  CHKERR mField.get_moab().get_connectivity(*frontEdges, front_nodes,
4152  true);
4153  Range front_faces = get_adj(front_nodes, 2);
4154  front_faces = subtract(front_faces, *crackFaces);
4155  auto front_tets = get_tets(front_nodes);
4156  auto front_side = intersect(side, front_tets);
4157  layers.push_back(front_side);
4158  for (;;) {
4159  auto adj_faces = get_skin(mField, layers.back());
4160  adj_faces = intersect(adj_faces, front_faces);
4161  auto adj_faces_tets = get_tets(adj_faces);
4162  adj_faces_tets = intersect(adj_faces_tets, front_tets);
4163  layers.push_back(unite(layers.back(), adj_faces_tets));
4164  if (layers.back().size() == layers[layers.size() - 2].size()) {
4165  break;
4166  }
4167  }
4168  MoFEMFunctionReturn(0);
4169  };
4170  CHK_THROW_MESSAGE(get(), "get_layers_for_sides");
4171  return layers;
4172  };
4173 
4174  auto sides_pair = get_two_sides_of_crack_surface(mField, *crackFaces);
4175  auto layers_top = get_layers_for_sides(sides_pair.first);
4176  auto layers_bottom = get_layers_for_sides(sides_pair.second);
4177 
4178 #ifndef NDEBUG
4179  if (debug) {
4180  CHKERR save_range(
4181  mField.get_moab(),
4182  "crack_tets_" +
4183  boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
4184  get_crack_adj_tets(*crackFaces));
4185  CHKERR save_range(mField.get_moab(), "sides_first.vtk", sides_pair.first);
4186  CHKERR save_range(mField.get_moab(), "sides_second.vtk",
4187  sides_pair.second);
4188  MOFEM_LOG("SELF", sev) << "Nb. layers " << layers_top.size();
4189  int l = 0;
4190  for (auto &r : layers_top) {
4191  MOFEM_LOG("SELF", sev) << "Layer " << l << " size " << r.size();
4192  CHKERR save_range(
4193  mField.get_moab(),
4194  "layers_top_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
4195  ++l;
4196  }
4197 
4198  l = 0;
4199  for (auto &r : layers_bottom) {
4200  MOFEM_LOG("SELF", sev) << "Layer " << l << " size " << r.size();
4201  CHKERR save_range(
4202  mField.get_moab(),
4203  "layers_bottom_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
4204  ++l;
4205  }
4206  }
4207 #endif
4208 
4209  auto get_cross = [&](auto t_dir, auto f) {
4210  FTensor::Tensor1<double, SPACE_DIM> t_normal;
4211  CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
4212  t_normal.normalize();
4213  FTensor::Tensor1<double, SPACE_DIM> t_cross;
4214  FTENSOR_INDEX(SPACE_DIM, i);
4215  FTENSOR_INDEX(SPACE_DIM, j);
4216  FTENSOR_INDEX(SPACE_DIM, k);
4217  t_cross(i) = FTensor::levi_civita(i, j, k) * t_normal(j) * t_dir(k);
4218  return t_cross;
4219  };
4220 
4221  auto get_sense = [&](auto f, auto e) {
4222  int side, sense, offset;
4223  CHK_MOAB_THROW(mField.get_moab().side_number(f, e, side, sense, offset),
4224  "get sense");
4225  return std::make_tuple(side, sense, offset);
4226  };
4227 
4228  auto calculate_edge_direction = [&](auto e) {
4229  const EntityHandle *conn;
4230  int num_nodes;
4231  CHKERR mField.get_moab().get_connectivity(e, conn, num_nodes, true);
4232  std::array<double, 6> coords;
4233  CHKERR mField.get_moab().get_coords(conn, num_nodes, coords.data());
4234  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p0{
4235  &coords[0], &coords[1], &coords[2]};
4236  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p1{
4237  &coords[3], &coords[4], &coords[5]};
4238  FTensor::Tensor1<double, SPACE_DIM> t_dir;
4239  FTENSOR_INDEX(SPACE_DIM, i);
4240  t_dir(i) = t_p1(i) - t_p0(i);
4241  auto l = t_dir.l2();
4242  t_dir(i) /= l;
4243  return t_dir;
4244  };
4245 
4246  auto evaluate_face_energy_and_set_orientation = [&](auto front_edges,
4247  auto front_faces,
4248  auto &sides_pair,
4249  auto th_position) {
4250  MoFEMFunctionBegin;
4251 
4252  Tag th_face_energy;
4253  CHKERR mField.get_moab().tag_get_handle("ReleaseEnergy", th_face_energy);
4254 
4255  /**
4256  * Iterate over front edges, get adjacent faces, find maximal face energy.
4257  * Maximal face energy is stored in the edge. The is then maximises across
4258  * all processors.
4259  *
4260  */
4261  auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
4262  auto &edge_face_max_energy_map) {
4263  MoFEMFunctionBegin;
4264  for (auto e : front_edges) {
4265  Range faces;
4266  CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
4267  faces = intersect(faces, front_faces);
4268 
4269  std::vector<double> face_energy(faces.size());
4270  CHKERR mField.get_moab().tag_get_data(th_face_energy, faces,
4271  face_energy.data());
4272  auto max_energy_it =
4273  std::max_element(face_energy.begin(), face_energy.end());
4274  double max_energy =
4275  max_energy_it != face_energy.end() ? *max_energy_it : 0;
4276  CHKERR mField.get_moab().tag_set_data(th_face_energy, &e, 1,
4277  &max_energy);
4278  edge_face_max_energy_map[e] =
4279  std::make_tuple(faces[max_energy_it - face_energy.begin()],
4280  max_energy, static_cast<double>(0));
4281  MOFEM_LOG("SELF", sev)
4282  << "Edge " << e << " max energy " << max_energy;
4283  };
4284  MoFEMFunctionReturn(0);
4285  };
4286 
4287  /**
4288  * For each front edge, find maximal face energy and orientation. This is
4289  * done by solving quadratic equation.
4290  *
4291  */
4292  auto calculate_face_orientation = [&](auto &edge_face_max_energy_map) {
4293  MoFEMFunctionBegin;
4294 
4295  auto up_down_face = [&](
4296 
4297  auto &face_angle_map_up,
4298  auto &face_angle_map_down
4299 
4300  ) {
4301  MoFEMFunctionBegin;
4302 
4303  for (auto &m : edge_face_max_energy_map) {
4304  auto e = m.first;
4305  auto [max_face, energy, opt_angle] = m.second;
4306 
4307  Range faces;
4308  CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
4309  faces = intersect(faces, front_faces);
4310  Range adj_tets; // tetrahedrons adjacent to the face
4311  CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
4312  false, adj_tets,
4313  moab::Interface::UNION);
4314  if (adj_tets.size()) {
4315 
4316  Range adj_tets; // tetrahedrons adjacent to the face
4317  CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
4318  false, adj_tets,
4319  moab::Interface::UNION);
4320  if (adj_tets.size()) {
4321 
4322  Range adj_tets_faces;
4323  // get faces
4324  CHKERR mField.get_moab().get_adjacencies(
4325  adj_tets, SPACE_DIM - 1, false, adj_tets_faces,
4326  moab::Interface::UNION);
4327  adj_tets_faces = intersect(adj_tets_faces, faces);
4328  if (adj_tets_faces.size() != 3 && adj_tets_faces.size() != 2) {
4329  MOFEM_LOG("SELF", Sev::error)
4330  << "Wrong number of crack faces " << adj_tets_faces.size()
4331  << " " << adj_tets.size();
4332  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
4333  "Wrong number of crack faces");
4334  }
4335 
4336  FTENSOR_INDEX(SPACE_DIM, i);
4337 
4338  // cross product of face normal and edge direction
4339  auto t_cross_max =
4340  get_cross(calculate_edge_direction(e), max_face);
4341  auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
4342  t_cross_max(i) *= sense_max;
4343 
4344  for (auto t : adj_tets) {
4345  Range adj_tets_faces;
4346  CHKERR mField.get_moab().get_adjacencies(
4347  &t, 1, SPACE_DIM - 1, false, adj_tets_faces);
4348  adj_tets_faces = intersect(adj_tets_faces, faces);
4349  adj_tets_faces =
4350  subtract(adj_tets_faces, Range(max_face, max_face));
4351 
4352  if (adj_tets_faces.size() == 1) {
4353 
4354  // cross product of adjacent face normal and edge
4355  // direction
4356  auto t_cross = get_cross(calculate_edge_direction(e),
4357  adj_tets_faces[0]);
4358  auto [side, sense, offset] =
4359  get_sense(adj_tets_faces[0], e);
4360  t_cross(i) *= sense;
4361  double dot = t_cross(i) * t_cross_max(i);
4362  auto angle = std::acos(dot);
4363 
4364  double energy;
4365  CHKERR mField.get_moab().tag_get_data(
4366  th_face_energy, adj_tets_faces, &energy);
4367 
4368  auto [side_face, sense_face, offset_face] =
4369  get_sense(t, max_face);
4370 
4371  if (sense_face > 0) {
4372  face_angle_map_up[e] =
4373  std::make_tuple(energy, angle, adj_tets_faces[0]);
4374 
4375  } else {
4376  face_angle_map_down[e] =
4377  std::make_tuple(energy, -angle, adj_tets_faces[0]);
4378  }
4379  }
4380  }
4381  }
4382  }
4383  }
4384 
4385  MoFEMFunctionReturn(0);
4386  };
4387 
4388  MatrixDouble A(3, 3);
4389  VectorDouble b(3);
4390 
4391  auto calc_optimal_angle_impl = [&](double a0, double e0, double a1,
4392  double e1, double a2, double e2) {
4393  A(0, 0) = a1 * a1;
4394  A(1, 0) = a1;
4395  A(2, 0) = 1;
4396  A(0, 1) = a2 * a2;
4397  A(1, 1) = a2;
4398  A(2, 1) = 1;
4399  A(0, 2) = a0 * a0;
4400  A(1, 2) = a0;
4401  A(2, 2) = 1;
4402 
4403  b[0] = e1;
4404  b[1] = e2;
4405  b[2] = e0;
4406 
4407  CHKERR solveLinearSystem(A, b);
4408 
4409  if (b[0] > -std::numeric_limits<float>::epsilon()) {
4410  return std::make_pair(e0, 0.);
4411 
4412  } else {
4413 
4414  auto optimal_angle = -b[1] / (2 * b[0]);
4415  auto optimal_energy = b[0] * optimal_angle * optimal_angle +
4416  b[1] * optimal_angle + b[2];
4417 
4418 #ifndef NDEBUG
4419  if (debug) {
4420  MOFEM_LOG("SELF", sev)
4421  << "Optimal_energy " << optimal_energy << " ( " << e0 << " "
4422  << (optimal_energy - e0) / e0 << "% ) " << optimal_angle
4423  << " e1 : a1 " << e1 << " : " << a1 << " e2 : a2 " << e2
4424  << " : " << a2;
4425 
4426  double angle = -M_PI;
4427  for (double a = -M_PI; a < M_PI; a += M_PI / 20.) {
4428  double energy = b[0] * a * a + b[1] * a + b[2];
4429  MOFEM_LOG("SELF", sev) << "PlotAngles " << a << " " << energy;
4430  }
4431  MOFEM_LOG("SELF", sev) << "PlotAngles " << endl;
4432 
4433  MOFEM_LOG("SELF", sev) << "AnglePts " << a1 << " " << e1;
4434  MOFEM_LOG("SELF", sev) << "AnglePts " << a0 << " " << e0;
4435  MOFEM_LOG("SELF", sev) << "AnglePts " << a2 << " " << e2;
4436  MOFEM_LOG("SELF", sev) << "AnglePts " << endl;
4437  }
4438 #endif // NDEBUG
4439 
4440  return std::make_pair(optimal_energy, optimal_angle);
4441  }
4442  };
4443 
4444  if (debug) {
4445 #ifndef NDEBUG
4446  auto [opt_e0, opt_a0] = calc_optimal_angle_impl(1, 1, 0, 0, 3, -3);
4447  if (std::abs(opt_e0 - 1) > std::numeric_limits<double>::epsilon()) {
4448  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong optimal energy");
4449  }
4450  if (std::abs(opt_a0 - 1) > std::numeric_limits<double>::epsilon()) {
4451  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong optimal angle");
4452  }
4453 
4454 #endif // NDEBUG
4455  }
4456 
4457  auto calc_optimal_angle = [&](
4458 
4459  auto &face_angle_map_up,
4460  auto &face_angle_map_down
4461 
4462  ) {
4463  MoFEMFunctionBegin;
4464 
4465  for (auto &m : edge_face_max_energy_map) {
4466  auto e = m.first;
4467  auto &[max_face, e0, a0] = m.second;
4468 
4469  if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
4470 
4471  if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
4472  face_angle_map_down.find(e) == face_angle_map_down.end()) {
4473 
4474  } else {
4475  auto [e1, a1, face_up] = face_angle_map_up.at(e);
4476  auto [e2, a2, face_down] = face_angle_map_down.at(e);
4477 
4478  MOFEM_LOG("SELF", sev) << "Edge " << e;
4479  auto [optimal_energy, optimal_angle] =
4480  calc_optimal_angle_impl(a0, e0, a1, e1, a2, e2);
4481 
4482  e0 = optimal_energy;
4483  a0 = optimal_angle;
4484  }
4485  }
4486  }
4487 
4488  MoFEMFunctionReturn(0);
4489  };
4490 
4491  std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
4492  face_angle_map_up;
4493  std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
4494  face_angle_map_down;
4495  CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
4496  CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
4497 
4498  if (debug) {
4499  auto th_angle = get_tags_vec("Angle", 1);
4500  Range up;
4501  for (auto &m : face_angle_map_up) {
4502  auto [e, a, face] = m.second;
4503  up.insert(face);
4504  CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
4505  }
4506  Range down;
4507  for (auto &m : face_angle_map_down) {
4508  auto [e, a, face] = m.second;
4509  down.insert(face);
4510  CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
4511  }
4512 
4513  Range max_energy_faces;
4514  for (auto &m : edge_face_max_energy_map) {
4515  auto [face, e, angle] = m.second;
4516  max_energy_faces.insert(face);
4517  CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1,
4518  &angle);
4519  }
4520  if (mField.get_comm_rank() == 0) {
4521  CHKERR save_range(mField.get_moab(), "up_faces.vtk", up);
4522  CHKERR save_range(mField.get_moab(), "down_faces.vtk", down);
4523  CHKERR save_range(mField.get_moab(), "max_energy_faces.vtk",
4524  max_energy_faces);
4525  }
4526  }
4527 
4528  MoFEMFunctionReturn(0);
4529  };
4530 
4531  auto sort_by_energy = [&](auto &edge_face_max_energy_map) {
4532  std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
4533  sort_by_energy;
4534  for (auto m : edge_face_max_energy_map) {
4535  auto e = m.first;
4536  auto [max_face, energy, opt_angle] = m.second;
4537  sort_by_energy[energy] = std::make_tuple(e, max_face, opt_angle);
4538  }
4539  return sort_by_energy;
4540  };
4541 
4542  auto set_face_orientation = [&](auto &sort_by_energy, auto &layers,
4543  auto &set_vertices, double &all_quality) {
4544  MoFEMFunctionBegin;
4545 
4546  Range body_ents;
4547  CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
4548  auto body_skin = get_skin(mField, body_ents);
4549  Range body_skin_edges;
4550  CHKERR mField.get_moab().get_adjacencies(
4551  body_skin, 1, false, body_skin_edges, moab::Interface::UNION);
4552  Range boundary_skin_verts;
4553  CHKERR mField.get_moab().get_connectivity(body_skin_edges,
4554  boundary_skin_verts, true);
4555 
4556  auto get_rotated_normal = [&](auto e, auto f, auto angle) {
4557  FTENSOR_INDEX(SPACE_DIM, i);
4558  FTENSOR_INDEX(SPACE_DIM, j);
4559  auto t_edge_dir = calculate_edge_direction(e);
4560  auto [side, sense, offset] = get_sense(f, e);
4561  t_edge_dir(i) *= sense;
4562  t_edge_dir.normalize();
4563  t_edge_dir(i) *= angle;
4564  auto t_R = LieGroups::SO3::exp(t_edge_dir, angle);
4565  FTensor::Tensor1<double, SPACE_DIM> t_normal;
4566  mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
4567  FTensor::Tensor1<double, SPACE_DIM> t_rotated_normal;
4568  t_rotated_normal(i) = t_R(i, j) * t_normal(j);
4569  return std::make_tuple(t_normal, t_rotated_normal);
4570  };
4571 
4572  Range rotated_normals;
4573 
4574  // iterate over energies
4575  for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
4576  ++it) {
4577  auto energy = it->first;
4578  auto [e, max_face, opt_angle] = it->second;
4579 
4580  // calculate rotation of max energy face
4581  auto [t_normal, t_rotated_normal] =
4582  get_rotated_normal(e, max_face, opt_angle);
4583  rotated_normals.insert(max_face);
4584 
4585  Range front_vertex; // vertices at crack front
4586  CHKERR mField.get_moab().get_connectivity(&e, 1, front_vertex, true);
4587  Range adj_tets; // adjacent tets to max face
4588  CHKERR mField.get_moab().get_adjacencies(&max_face, 1, 3, false,
4589  adj_tets);
4590 
4591 #ifndef NDEBUG
4592  if (adj_tets.size() != 2) {
4593  MOFEM_LOG("SELF", sev)
4594  << "Wrong number of tets adjacent to max face "
4595  << adj_tets.size();
4596  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
4597  "Wrong number of crack faces");
4598  }
4599 #endif // NDEBUG
4600  Range adj_front_faces; // adjacent faces to front edge
4601  CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false,
4602  adj_front_faces);
4603  adj_front_faces = subtract(adj_front_faces, *crackFaces);
4604  Range adj_tet_faces; // adjacent faces to adjacent tets
4605  CHKERR mField.get_moab().get_adjacencies(
4606  adj_tets, 2, false, adj_tet_faces, moab::Interface::UNION);
4607  adj_tet_faces =
4608  intersect(adj_tet_faces, adj_front_faces); // only three faces
4609 #ifndef NDEBUG
4610  if (adj_tet_faces.size() != 3) {
4611  MOFEM_LOG("SELF", sev) << "Wrong number of faces adj. to test "
4612  << adj_tet_faces.size();
4613  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
4614  "Wrong number of crack faces");
4615  }
4616 #endif // NDEBUG
4617 
4618  Range adj_front_faces_edges; // edges adjacent to front faces
4619  CHKERR mField.get_moab().get_adjacencies(adj_front_faces, 1, false,
4620  adj_front_faces_edges,
4621  moab::Interface::UNION);
4622 
4623  // only process to faces
4624  std::array<EntityHandle, 3> processed_faces{0, max_face, 0};
4625  int set_index = 0;
4626  for (auto &l : layers) {
4627  auto adj_tets_layer = intersect(adj_tets, l);
4628  if (adj_tets_layer.empty()) {
4629  continue;
4630  }
4631  Range faces_layer; // adjacent faces to adjacent tets
4632  CHKERR mField.get_moab().get_adjacencies(
4633  adj_tets_layer, 2, false, faces_layer, moab::Interface::UNION);
4634  faces_layer = intersect(faces_layer, adj_front_faces);
4635  faces_layer = subtract(faces_layer, Range(max_face, max_face));
4636  if (set_index > 0) {
4637  faces_layer = subtract(
4638  faces_layer, Range(processed_faces[0], processed_faces[0]));
4639  }
4640 #ifndef NDEBUG
4641  if (faces_layer.size() != 1) {
4642  MOFEM_LOG("SELF", sev)
4643  << "Wrong number of faces to move " << faces_layer.size();
4644  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
4645  "Wrong number of crack faces");
4646  }
4647 #endif // NDEBUG
4648  processed_faces[set_index] = faces_layer[0];
4649  set_index = 2;
4650  }
4651 
4652  // choose face, if one of the vertices is already set
4653  for (auto f : {0, 1, 2}) {
4654 
4655  if (processed_faces[f] == 0)
4656  break;
4657 
4658  Range vertex; // vertex to move
4659  auto rval = mField.get_moab().get_connectivity(&processed_faces[f],
4660  1, vertex, true);
4661  if (rval != MB_SUCCESS) {
4662  MOFEM_LOG("SELF", Sev::error)
4663  << "get_connectivity failed " << f << " "
4664  << moab::CN::EntityTypeName(type_from_handle(f));
4665  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4666  "can noy get connectivity");
4667  }
4668  vertex = subtract(vertex, front_vertex);
4669  if (set_vertices.find(vertex[0]) != set_vertices.end()) {
4670  if (f == 0) {
4671  processed_faces[1] = 0;
4672  processed_faces[2] = 0;
4673  } else if (f == 1) {
4674  processed_faces[0] = 0;
4675  processed_faces[2] = 0;
4676  } else {
4677  processed_faces[0] = 0;
4678  processed_faces[1] = 0;
4679  }
4680  }
4681  }
4682 
4683  // [quality] -> [vertex, face, position]
4684  std::map<double, std::tuple<EntityHandle, EntityHandle,
4685  FTensor::Tensor1<double, SPACE_DIM>>>
4686  sort_by_quality;
4687 
4688  int face_idx = 0;
4689  for (auto f : processed_faces) {
4690 
4691  if (f == 0)
4692  continue;
4693 
4694  Range vertex; // vertex to move
4695  auto rval = mField.get_moab().get_connectivity(&f, 1, vertex, true);
4696  if (rval != MB_SUCCESS) {
4697  MOFEM_LOG("SELF", Sev::error)
4698  << "get_connectivity failed " << f << " "
4699  << moab::CN::EntityTypeName(type_from_handle(f));
4700  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4701  "can noy get connectivity");
4702  }
4703  vertex = subtract(vertex, front_vertex);
4704 #ifndef NDEBUG
4705  if (vertex.size() != 1) {
4706  MOFEM_LOG("SELF", sev)
4707  << "Wrong number of vertex to move " << vertex.size();
4708  }
4709 #endif // NDEBUG
4710  FTensor::Tensor1<double, 3>
4711  t_vertex_coords; // coords of vertex to move
4712  CHKERR mField.get_moab().get_coords(vertex, &t_vertex_coords(0));
4713 
4714  Range vertex_edges; // edges adjacent to vertex
4715  CHKERR mField.get_moab().get_adjacencies(vertex, 1, false,
4716  vertex_edges);
4717  vertex_edges = subtract(vertex_edges, adj_front_faces_edges);
4718  if (boundary_skin_verts.size() &&
4719  boundary_skin_verts.find(vertex[0]) !=
4720  boundary_skin_verts.end()) {
4721  MOFEM_LOG("SELF", sev) << "Boundary vertex";
4722  vertex_edges = intersect(vertex_edges, body_skin_edges);
4723  }
4724  if (geometry_edges_verts.size() &&
4725  geometry_edges_verts.find(vertex[0]) !=
4726  geometry_edges_verts.end()) {
4727  MOFEM_LOG("SELF", sev) << "Geometry edge vertex";
4728  vertex_edges = intersect(vertex_edges, geometry_edges);
4729  }
4730  if (crack_faces_verts.size() &&
4731  crack_faces_verts.find(vertex_edges[0]) !=
4732  crack_faces_verts.end()) {
4733  MOFEM_LOG("SELF", sev) << "Crack face vertex";
4734  vertex_edges = intersect(vertex_edges, crack_faces_edges);
4735  }
4736 
4737  EntityHandle f0 = front_vertex[0];
4738  EntityHandle f1 = front_vertex[1];
4739  FTensor::Tensor1<double, 3> t_v_e0, t_v_e1;
4740  CHKERR mField.get_moab().get_coords(&f0, 1, &t_v_e0(0));
4741  CHKERR mField.get_moab().get_coords(&f1, 1, &t_v_e1(0));
4742 
4743  FTENSOR_INDEX(SPACE_DIM, i);
4744  std::map<EntityHandle, FTensor::Tensor1<double, SPACE_DIM>>
4745  node_positions_map;
4746  for (auto e : vertex_edges) {
4747  Range edge_conn;
4748  CHKERR mField.get_moab().get_connectivity(&e, 1, edge_conn, true);
4749  edge_conn = subtract(edge_conn, vertex);
4750 #ifndef NDEBUG
4751  if (edge_conn.size() != 1) {
4752  MOFEM_LOG("SELF", sev)
4753  << "Wrong number of edge conn " << edge_conn.size();
4754  }
4755 #endif // NDEBUG
4756  auto it = set_vertices.find(vertex[0]);
4757  if (it != set_vertices.end()) {
4758  node_positions_map[e](i) = std::get<2>(it->second)(i);
4759  } else {
4760  FTensor::Tensor1<double, 3> t_v0;
4761  t_v0(i) = (t_v_e0(i) + t_v_e1(i)) / 2;
4762  FTensor::Tensor1<double, 3> t_v3;
4763  CHKERR mField.get_moab().get_coords(edge_conn, &t_v3(0));
4764  auto a = (t_v0(i) - t_vertex_coords(i)) * t_rotated_normal(i);
4765  auto b = (t_v3(i) - t_vertex_coords(i)) * t_rotated_normal(i);
4766  auto gamma = a / b;
4767  MOFEM_LOG("SELF", sev)
4768  << face_idx << " edge: " << e << " gamma: " << gamma;
4769  if (std::isnormal(gamma) &&
4770  gamma > -std::numeric_limits<double>::epsilon() &&
4771  gamma < 1 - std::numeric_limits<double>::epsilon()) {
4772  for (auto i : {0, 1, 2}) {
4773  node_positions_map[e](i) =
4774  gamma * (t_v3(i) - t_vertex_coords(i));
4775  }
4776  }
4777  }
4778  }
4779  if (vertex_edges.size() == 0) {
4780  node_positions_map[0](i) = 0;
4781  }
4782 
4783  Range adj_vertex_tets;
4784  CHKERR mField.get_moab().get_adjacencies(vertex, 3, false,
4785  adj_vertex_tets);
4786  for (auto &p : node_positions_map) {
4787  double min_quality = 1;
4788  for (auto t : adj_vertex_tets) {
4789  const EntityHandle *conn;
4790  int num_nodes;
4791  CHKERR mField.get_moab().get_connectivity(t, conn, num_nodes,
4792  true);
4793  std::array<double, 12> coords;
4794  CHKERR mField.get_moab().get_coords(conn, num_nodes,
4795  coords.data());
4796  int idx = 0;
4797  for (; idx < num_nodes; ++idx) {
4798  if (conn[idx] == vertex[0]) {
4799  break;
4800  }
4801  }
4802 
4803  if (set_vertices.find(vertex[0]) == set_vertices.end()) {
4804  for (auto ii : {0, 1, 2}) {
4805  coords[3 * idx + ii] += p.second(ii);
4806  }
4807  }
4808  for (auto vv : {0, 1, 2, 3}) {
4809  auto it = set_vertices.find(conn[vv]);
4810  if (it != set_vertices.end()) {
4811  auto &[f, energy, t_position] = it->second;
4812  for (auto ii : {0, 1, 2})
4813  coords[3 * vv + ii] += t_position(ii);
4814  }
4815  }
4816 
4817  double q = Tools::volumeLengthQuality(coords.data());
4818  if (!std::isnormal(q))
4819  q = -2;
4820  min_quality = std::min(min_quality, q);
4821  }
4822  sort_by_quality[min_quality] =
4823  std::make_tuple(vertex[0], f,
4824  FTensor::Tensor1<double, SPACE_DIM>{
4825  p.second(0), p.second(1), p.second(2)});
4826  }
4827 
4828  ++face_idx;
4829  }
4830 
4831  if (debug) {
4832  for (auto s : sort_by_quality) {
4833  auto &[vertex, face, t_position] = s.second;
4834  MOFEM_LOG("SELF", sev)
4835  << "Quality: " << s.first << " vertex " << vertex << " face "
4836  << face << " point: " << t_position;
4837  }
4838  }
4839  for (auto it = sort_by_quality.rbegin(); it != sort_by_quality.rend();
4840  ++it) {
4841  auto &[vertex, face, t_position] = it->second;
4842  all_quality = std::min(all_quality, it->first);
4843  MOFEM_LOG("SELF", sev)
4844  << "Set quality: " << it->first << " vertex " << vertex
4845  << " face " << face << " point: " << t_position;
4846  set_vertices[vertex] = std::make_tuple(face, energy, t_position);
4847  break;
4848  }
4849  }
4850 
4851  MoFEMFunctionReturn(0);
4852  };
4853 
4854  // map: {edge, {face, energy, optimal_angle}}
4855  std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
4856  edge_face_max_energy_map;
4857  CHKERR find_maximal_face_energy(front_edges, front_faces,
4858  edge_face_max_energy_map);
4859  CHKERR calculate_face_orientation(edge_face_max_energy_map);
4860 
4861  auto edge_map_sort_by_energy = sort_by_energy(edge_face_max_energy_map);
4862 
4863  MOFEM_LOG("SELF", sev) << "Top";
4864  double top_quality = 1;
4865  std::map<EntityHandle, std::tuple<EntityHandle, double,
4866  FTensor::Tensor1<double, SPACE_DIM>>>
4867  set_vertices_top;
4868  CHKERR set_face_orientation(edge_map_sort_by_energy, layers_top,
4869  set_vertices_top, top_quality);
4870  MOFEM_LOG("SELF", sev) << "Bottom";
4871  double bottom_quality = 1;
4872  std::map<EntityHandle, std::tuple<EntityHandle, double,
4873  FTensor::Tensor1<double, SPACE_DIM>>>
4874  set_vertices_bottom;
4875  CHKERR set_face_orientation(edge_map_sort_by_energy, layers_bottom,
4876  set_vertices_bottom, bottom_quality);
4877 
4878  MOFEM_LOG("SELF", sev) << "Top quality " << top_quality;
4879  MOFEM_LOG("SELF", sev) << "Bottom quality " << bottom_quality;
4880 
4881  auto set_tag = [&](auto &set_vertices, auto th_position) {
4882  MoFEMFunctionBegin;
4883  for (auto m : set_vertices) {
4884  auto v = m.first;
4885  auto &[f, energy, t_p] = m.second;
4886  CHKERR mField.get_moab().tag_set_data(th_position[0], &v, 1, &t_p(0));
4887  CHKERR mField.get_moab().tag_set_data(th_max_face_energy[0], &f, 1,
4888  &energy);
4889  }
4890  MoFEMFunctionReturn(0);
4891  };
4892 
4893  if (top_quality > bottom_quality) {
4894  CHKERR set_tag(set_vertices_top, th_position);
4895  } else {
4896  CHKERR set_tag(set_vertices_bottom, th_position);
4897  }
4898 
4899  if (debug) {
4900  auto th_front_top = get_tags_vec("FrontPositionTop", 3);
4901  auto th_front_bottom = get_tags_vec("FrontPositionBottom", 3);
4902 
4903  Range top_moved_faces;
4904  for (auto m : set_vertices_top) {
4905  auto v = m.first;
4906  auto &[f, energy, t_p] = m.second;
4907  top_moved_faces.insert(f);
4908  CHKERR mField.get_moab().tag_set_data(th_front_top[0], &v, 1,
4909  &t_p(0));
4910  }
4911  Range bottom_moved_faces;
4912  for (auto m : set_vertices_bottom) {
4913  auto v = m.first;
4914  auto &[f, energy, t_p] = m.second;
4915  bottom_moved_faces.insert(f);
4916  CHKERR mField.get_moab().tag_set_data(th_front_bottom[0], &v, 1,
4917  &t_p(0));
4918  }
4919 
4920  for (auto m : set_vertices_bottom) {
4921  auto &[f, energy, t_p] = m.second;
4922  bottom_moved_faces.insert(f);
4923  }
4924 
4925  CHKERR save_range(mField.get_moab(), "top_moved_faces.vtk",
4926  top_moved_faces);
4927  CHKERR save_range(mField.get_moab(), "bottom_moved_faces.vtk",
4928  bottom_moved_faces);
4929  auto front_faces = get_adj_front(false);
4930  CHKERR save_range(mField.get_moab(), "front_move.vtk", front_faces);
4931  }
4932 
4933  MoFEMFunctionReturn(0);
4934  };
4935 
4936  MOFEM_LOG("SELF", sev) << "Front edges " << frontEdges->size();
4937  CHKERR evaluate_face_energy_and_set_orientation(
4938  *frontEdges, get_adj_front(false), sides_pair, th_front_position);
4939  }
4940 
4941  CHKERR VecZeroEntries(vertexExchange.second);
4942  CHKERR VecGhostUpdateBegin(vertexExchange.second, INSERT_VALUES,
4943  SCATTER_FORWARD);
4944  CHKERR VecGhostUpdateEnd(vertexExchange.second, INSERT_VALUES,
4945  SCATTER_FORWARD);
4946  CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
4947  mField.get_moab(), vertexExchange, th_front_position[0]);
4948  CHKERR VecZeroEntries(faceExchange.second);
4949  CHKERR VecGhostUpdateBegin(faceExchange.second, INSERT_VALUES,
4950  SCATTER_FORWARD);
4951  CHKERR VecGhostUpdateEnd(faceExchange.second, INSERT_VALUES, SCATTER_FORWARD);
4952  CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
4953  mField.get_moab(), faceExchange, th_max_face_energy[0]);
4954 
4955  auto get_max_moved_faces = [&]() {
4956  Range max_moved_faces;
4957  auto adj_front = get_adj_front(false);
4958  std::vector<double> face_energy(adj_front.size());
4959  CHKERR mField.get_moab().tag_get_data(th_max_face_energy[0], adj_front,
4960  face_energy.data());
4961  for (int i = 0; i != adj_front.size(); ++i) {
4962  if (face_energy[i] > std::numeric_limits<double>::epsilon()) {
4963  max_moved_faces.insert(adj_front[i]);
4964  }
4965  }
4966 
4967  return boost::make_shared<Range>(max_moved_faces);
4968  };
4969 
4970  maxMovedFaces = get_max_moved_faces();
4971 
4972  if (debug) {
4973  Range tets;
4974  CHKERR mField.get_moab().get_entities_by_dimension(0, 3, tets);
4975  CHKERR save_range(
4976  mField.get_moab(),
4977  "projected_tets_" +
4978  boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
4979  tets);
4980  CHKERR save_range(
4981  mField.get_moab(),
4982  "max_moved_faces_" +
4983  boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
4984  *maxMovedFaces);
4985  }
4986 
4987  MoFEMFunctionReturn(0);
4988 }

◆ createCrackSurfaceMeshset()

MoFEMErrorCode EshelbianCore::createCrackSurfaceMeshset ( )
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 5214 of file EshelbianPlasticity.cpp.

5214  {
5215  MoFEMFunctionBegin;
5216  auto meshset_mng = mField.getInterface<MeshsetsManager>();
5217  while (meshset_mng->checkMeshset(BLOCKSET, addCrackMeshsetId))
5218  ++addCrackMeshsetId;
5219  MOFEM_LOG("EP", Sev::inform)
5220  << "Crack added surface meshset " << addCrackMeshsetId;
5221  CHKERR meshset_mng->addMeshset(BLOCKSET, addCrackMeshsetId, "CRACK_COMPUTED");
5222  MoFEMFunctionReturn(0);
5223 };

◆ createExchangeVectors()

MoFEMErrorCode EshelbianCore::createExchangeVectors ( Sev  sev)
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 5356 of file EshelbianPlasticity.cpp.

5356  {
5357  MoFEMFunctionBegin;
5358 
5359  volumeExchange = CommInterface::createEntitiesPetscVector(
5360  mField.get_comm(), mField.get_moab(), 3, 1, sev);
5361  faceExchange = CommInterface::createEntitiesPetscVector(
5362  mField.get_comm(), mField.get_moab(), 2, 1, Sev::inform);
5363  edgeExchange = CommInterface::createEntitiesPetscVector(
5364  mField.get_comm(), mField.get_moab(), 1, 1, Sev::inform);
5365  vertexExchange = CommInterface::createEntitiesPetscVector(
5366  mField.get_comm(), mField.get_moab(), 0, 3, Sev::inform);
5367 
5368  auto print_loc_size = [this](auto v, auto str, auto sev) {
5369  MoFEMFunctionBegin;
5370  int size;
5371  CHKERR VecGetLocalSize(v.second, &size);
5372  int low, high;
5373  CHKERR VecGetOwnershipRange(v.second, &low, &high);
5374  MOFEM_LOG("EPSYNC", sev) << str << " local size " << size << " ( " << low
5375  << " " << high << " ) ";
5376  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), sev);
5377  MoFEMFunctionReturn(0);
5378  };
5379 
5380  CHKERR print_loc_size(volumeExchange, "volumeExchange", sev);
5381  CHKERR print_loc_size(faceExchange, "faceExchange", sev);
5382  CHKERR print_loc_size(edgeExchange, "edgeExchange", sev);
5383  CHKERR print_loc_size(vertexExchange, "vertexExchange", sev);
5384 
5385  MoFEMFunctionReturn(0);
5386 }

◆ d_f_linear()

static double EshelbianCore::d_f_linear ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 68 of file EshelbianCore.hpp.

68 { return 1; }

◆ d_f_log()

static double EshelbianCore::d_f_log ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 49 of file EshelbianCore.hpp.

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

◆ d_f_log_e()

static double EshelbianCore::d_f_log_e ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 40 of file EshelbianCore.hpp.

40 { return std::exp(v); }

◆ dd_f_linear()

static double EshelbianCore::dd_f_linear ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 69 of file EshelbianCore.hpp.

69 { return 0; }

◆ dd_f_log()

static double EshelbianCore::dd_f_log ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 52 of file EshelbianCore.hpp.

52  {
53  return pow(EshelbianCore::exponentBase, v) *
54  log(EshelbianCore::exponentBase) * log(EshelbianCore::exponentBase);
55  }

◆ dd_f_log_e()

static double EshelbianCore::dd_f_log_e ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 41 of file EshelbianCore.hpp.

41 { return std::exp(v); }

◆ f_linear()

static double EshelbianCore::f_linear ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 67 of file EshelbianCore.hpp.

67 { return v + 1; }

◆ f_log()

static double EshelbianCore::f_log ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 46 of file EshelbianCore.hpp.

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

◆ f_log_e()

static double EshelbianCore::f_log_e ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 39 of file EshelbianCore.hpp.

39 { return std::exp(v); }

◆ getBc()

template<typename BC >
MoFEMErrorCode EshelbianCore::getBc ( boost::shared_ptr< BC > &  bc_vec_ptr,
const std::string  block_name,
const int  nb_attributes 
)
inline
Examples
EshelbianPlasticity.cpp.

Definition at line 143 of file EshelbianCore.hpp.

144  {
145  MoFEMFunctionBegin;
146  for (auto it :
147  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
148 
149  (boost::format("%s(.*)") % block_name).str()
150 
151  ))
152 
153  ) {
154  std::vector<double> block_attributes;
155  CHKERR it->getAttributes(block_attributes);
156  if (block_attributes.size() < nb_attributes) {
157  SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
158  "In block %s expected %d attributes, but given %d",
159  it->getName().c_str(), nb_attributes, block_attributes.size());
160  }
161  Range faces;
162  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
163  true);
164  bc_vec_ptr->emplace_back(it->getName(), block_attributes, faces);
165  }
166  MoFEMFunctionReturn(0);
167  }

◆ getOptions()

MoFEMErrorCode EshelbianCore::getOptions ( )
Examples
EshelbianPlasticity.cpp.

Definition at line 883 of file EshelbianPlasticity.cpp.

883  {
884  MoFEMFunctionBegin;
885  const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
886  const char *list_symm[] = {"symm", "not_symm"};
887  PetscInt choice_rot = EshelbianCore::rotSelector;
888  PetscInt choice_grad = EshelbianCore::gradApproximator;
889  PetscInt choice_symm = EshelbianCore::symmetrySelector;
890 
891  const char *list_stretches[] = {"linear", "log"};
892  PetscInt choice_stretch = StretchSelector::LOG;
893 
894  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
895  "none");
896  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
897  spaceOrder, &spaceOrder, PETSC_NULL);
898  CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
899  spaceH1Order, &spaceH1Order, PETSC_NULL);
900  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
901  "", materialOrder, &materialOrder, PETSC_NULL);
902  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
903  &alphaU, PETSC_NULL);
904  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
905  &alphaW, PETSC_NULL);
906  CHKERR PetscOptionsScalar("-viscosity_alpha_omega", "rot viscosity", "",
907  alphaOmega, &alphaOmega, PETSC_NULL);
908  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
909  &alphaRho, PETSC_NULL);
910  CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
911  LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
912  PETSC_NULL);
913  CHKERR PetscOptionsEList("-grad", "gradient of defamation approximate", "",
914  list_rots, NO_H1_CONFIGURATION + 1,
915  list_rots[choice_grad], &choice_grad, PETSC_NULL);
916  CHKERR PetscOptionsEList("-symm", "symmetric variant", "", list_symm, 2,
917  list_symm[choice_symm], &choice_symm, PETSC_NULL);
918 
919  CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
920  &EshelbianCore::exponentBase, PETSC_NULL);
921  CHKERR PetscOptionsEList(
922  "-stretches", "stretches", "", list_stretches, StretchSelector::LOG + 1,
923  list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
924 
925  CHKERR PetscOptionsBool("-no_stretch", "do not solve for stretch", "",
926  noStretch, &noStretch, PETSC_NULL);
927  CHKERR PetscOptionsBool("-set_singularity", "set singularity", "",
928  setSingularity, &setSingularity, PETSC_NULL);
929  CHKERR PetscOptionsBool("-l2_user_base_scale", "streach scale", "",
930  l2UserBaseScale, &l2UserBaseScale, PETSC_NULL);
931 
932  // dynamic relaxation
933  CHKERR PetscOptionsBool("-dynamic_relaxation", "dynamic time relaxation", "",
934  dynamicRelaxation, &dynamicRelaxation, PETSC_NULL);
935 
936  // contact parameters
937  CHKERR PetscOptionsInt("-contact_max_post_proc_ref_level", "refinement level",
938  "", contactRefinementLevels, &contactRefinementLevels,
939  PETSC_NULL);
940 
941  // cracking parameters
942  CHKERR PetscOptionsBool("-cracking_on", "cracking ON", "", crackingOn,
943  &crackingOn, PETSC_NULL);
944  CHKERR PetscOptionsScalar("-griffith_energy", "Griffith energy", "",
945  griffithEnergy, &griffithEnergy, PETSC_NULL);
946 
947  ierr = PetscOptionsEnd();
948  CHKERRG(ierr);
949 
950  if (setSingularity) {
951  l2UserBaseScale = PETSC_TRUE;
952  }
953 
954  EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
955  EshelbianCore::gradApproximator = static_cast<RotSelector>(choice_grad);
956  EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
957  EshelbianCore::symmetrySelector = static_cast<SymmetrySelector>(choice_symm);
958 
959  switch (EshelbianCore::stretchSelector) {
960  case StretchSelector::LINEAR:
961  EshelbianCore::f = f_linear;
962  EshelbianCore::d_f = d_f_linear;
963  EshelbianCore::dd_f = dd_f_linear;
964  EshelbianCore::inv_f = inv_f_linear;
965  EshelbianCore::inv_d_f = inv_d_f_linear;
966  EshelbianCore::inv_dd_f = inv_dd_f_linear;
967  break;
968  case StretchSelector::LOG:
969  if (EshelbianCore::exponentBase != exp(1)) {
970  EshelbianCore::f = f_log;
971  EshelbianCore::d_f = d_f_log;
972  EshelbianCore::dd_f = dd_f_log;
973  EshelbianCore::inv_f = inv_f_log;
974  EshelbianCore::inv_d_f = inv_d_f_log;
975  EshelbianCore::inv_dd_f = inv_dd_f_log;
976  } else {
977  EshelbianCore::f = f_log_e;
978  EshelbianCore::d_f = d_f_log_e;
979  EshelbianCore::dd_f = dd_f_log_e;
980  EshelbianCore::inv_f = inv_f_log;
981  EshelbianCore::inv_d_f = inv_d_f_log;
982  EshelbianCore::inv_dd_f = inv_dd_f_log;
983  }
984  break;
985  default:
986  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
987  break;
988  };
989 
990  MOFEM_LOG("EP", Sev::inform) << "spaceOrder " << spaceOrder;
991  MOFEM_LOG("EP", Sev::inform) << "spaceH1Order " << spaceH1Order;
992  MOFEM_LOG("EP", Sev::inform) << "materialOrder " << materialOrder;
993  MOFEM_LOG("EP", Sev::inform) << "alphaU " << alphaU;
994  MOFEM_LOG("EP", Sev::inform) << "alphaW " << alphaW;
995  MOFEM_LOG("EP", Sev::inform) << "alphaOmega " << alphaOmega;
996  MOFEM_LOG("EP", Sev::inform) << "alphaRho " << alphaRho;
997  MOFEM_LOG("EP", Sev::inform)
998  << "Rotations " << list_rots[EshelbianCore::rotSelector];
999  MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
1000  << list_rots[EshelbianCore::gradApproximator];
1001  MOFEM_LOG("EP", Sev::inform)
1002  << "Symmetry " << list_symm[EshelbianCore::symmetrySelector];
1003  if (exponentBase != exp(1))
1004  MOFEM_LOG("EP", Sev::inform)
1005  << "Base exponent " << EshelbianCore::exponentBase;
1006  else
1007  MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
1008  MOFEM_LOG("EP", Sev::inform) << "Stretch " << list_stretches[choice_stretch];
1009 
1010  MOFEM_LOG("EP", Sev::inform) << "Dynamic relaxation " << (dynamicRelaxation)
1011  ? "yes"
1012  : "no";
1013  MOFEM_LOG("EP", Sev::inform) << "Singularity " << (setSingularity) ? "yes"
1014  : "no";
1015  MOFEM_LOG("EP", Sev::inform) << "L2 user base scale " << (l2UserBaseScale)
1016  ? "yes"
1017  : "no";
1018 
1019  MOFEM_LOG("EP", Sev::inform) << "Griffith energy " << griffithEnergy;
1020 
1021  if (spaceH1Order == -1)
1022  spaceH1Order = spaceOrder;
1023 
1024  MoFEMFunctionReturn(0);
1025 }

◆ getSpatialDispBc()

MoFEMErrorCode EshelbianCore::getSpatialDispBc ( )

[Getting norms]

Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 5290 of file EshelbianPlasticity.cpp.

5290  {
5291  MoFEMFunctionBegin;
5292 
5293  auto bc_mng = mField.getInterface<BcManager>();
5294  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
5295  "", piolaStress, false, false);
5296 
5297  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
5298 
5299  for (auto bc : bc_mng->getBcMapByBlockName()) {
5300  if (auto disp_bc = bc.second->dispBcPtr) {
5301 
5302  MOFEM_LOG("EP", Sev::noisy) << *disp_bc;
5303 
5304  std::vector<double> block_attributes(6, 0.);
5305  if (disp_bc->data.flag1 == 1) {
5306  block_attributes[0] = disp_bc->data.value1;
5307  block_attributes[3] = 1;
5308  }
5309  if (disp_bc->data.flag2 == 1) {
5310  block_attributes[1] = disp_bc->data.value2;
5311  block_attributes[4] = 1;
5312  }
5313  if (disp_bc->data.flag3 == 1) {
5314  block_attributes[2] = disp_bc->data.value3;
5315  block_attributes[5] = 1;
5316  }
5317  auto faces = bc.second->bcEnts.subset_by_dimension(2);
5318  bcSpatialDispVecPtr->emplace_back(bc.first, block_attributes, faces);
5319  }
5320  }
5321 
5322  // old way of naming blocksets for displacement BCs
5323  CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
5324 
5325  MoFEMFunctionReturn(0);
5326 }

◆ getSpatialRotationBc()

MoFEMErrorCode EshelbianCore::getSpatialRotationBc ( )
inline
Examples
ep.cpp.

Definition at line 171 of file EshelbianCore.hpp.

171  {
172  bcSpatialRotationVecPtr = boost::make_shared<BcRotVec>();
173  return getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_BC", 4);
174  }

◆ getSpatialTractionBc()

MoFEMErrorCode EshelbianCore::getSpatialTractionBc ( )
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 5328 of file EshelbianPlasticity.cpp.

5328  {
5329  MoFEMFunctionBegin;
5330 
5331  auto bc_mng = mField.getInterface<BcManager>();
5332  CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
5333  false, false);
5334 
5335  bcSpatialTraction = boost::make_shared<TractionBcVec>();
5336 
5337  for (auto bc : bc_mng->getBcMapByBlockName()) {
5338  if (auto force_bc = bc.second->forceBcPtr) {
5339  std::vector<double> block_attributes(6, 0.);
5340  block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
5341  block_attributes[3] = 1;
5342  block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
5343  block_attributes[4] = 1;
5344  block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
5345  block_attributes[5] = 1;
5346  auto faces = bc.second->bcEnts.subset_by_dimension(2);
5347  bcSpatialTraction->emplace_back(bc.first, block_attributes, faces);
5348  }
5349  }
5350 
5351  CHKERR getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
5352 
5353  MoFEMFunctionReturn(0);
5354 }

◆ getSpatialTractionFreeBc()

MoFEMErrorCode EshelbianCore::getSpatialTractionFreeBc ( const EntityHandle  meshset = 0)
inline
Examples
ep.cpp.

Definition at line 193 of file EshelbianCore.hpp.

193  {
194  bcSpatialFreeTraction =
195  boost::shared_ptr<TractionFreeBc>(new TractionFreeBc());
196  return getTractionFreeBc(meshset, bcSpatialFreeTraction, "CONTACT");
197  }

◆ gettingNorms()

MoFEMErrorCode EshelbianCore::gettingNorms ( )

[Getting norms]

Examples
EshelbianPlasticity.cpp.

Definition at line 5226 of file EshelbianPlasticity.cpp.

5226  {
5227  MoFEMFunctionBegin;
5228 
5229  auto post_proc_norm_fe =
5230  boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
5231 
5232  auto post_proc_norm_rule_hook = [](int, int, int p) -> int {
5233  return 2 * (p);
5234  };
5235  post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
5236 
5237  post_proc_norm_fe->getUserPolynomialBase() =
5238  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
5239 
5240  CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
5241  post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV}, materialH1Positions,
5242  frontAdjEdges);
5243 
5244  enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
5245  auto norms_vec =
5246  createVectorMPI(mField.get_comm(), LAST_NORM, PETSC_DETERMINE);
5247  CHKERR VecZeroEntries(norms_vec);
5248 
5249  auto u_l2_ptr = boost::make_shared<MatrixDouble>();
5250  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
5251  post_proc_norm_fe->getOpPtrVector().push_back(
5252  new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, u_l2_ptr));
5253  post_proc_norm_fe->getOpPtrVector().push_back(
5254  new OpCalculateVectorFieldValues<SPACE_DIM>(spatialH1Disp, u_h1_ptr));
5255  post_proc_norm_fe->getOpPtrVector().push_back(
5256  new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_NORM_L2));
5257  post_proc_norm_fe->getOpPtrVector().push_back(
5258  new OpCalcNormL2Tensor1<SPACE_DIM>(u_h1_ptr, norms_vec, U_NORM_H1));
5259  post_proc_norm_fe->getOpPtrVector().push_back(
5260  new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_ERROR_L2,
5261  u_h1_ptr));
5262 
5263  auto piola_ptr = boost::make_shared<MatrixDouble>();
5264  post_proc_norm_fe->getOpPtrVector().push_back(
5265  new OpCalculateHVecTensorField<3, 3>(piolaStress, piola_ptr));
5266  post_proc_norm_fe->getOpPtrVector().push_back(
5267  new OpCalcNormL2Tensor2<3, 3>(piola_ptr, norms_vec, PIOLA_NORM));
5268 
5269  TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
5270  CHKERR mField.loop_finite_elements("ELASTIC_PROBLEM", elementVolumeName,
5271  *post_proc_norm_fe);
5272  TetPolynomialBase::switchCacheBaseOff<HDIV>({post_proc_norm_fe.get()});
5273 
5274  CHKERR VecAssemblyBegin(norms_vec);
5275  CHKERR VecAssemblyEnd(norms_vec);
5276  const double *norms;
5277  CHKERR VecGetArrayRead(norms_vec, &norms);
5278  MOFEM_LOG("EP", Sev::inform) << "norm_u: " << std::sqrt(norms[U_NORM_L2]);
5279  MOFEM_LOG("EP", Sev::inform) << "norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
5280  MOFEM_LOG("EP", Sev::inform)
5281  << "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
5282  MOFEM_LOG("EP", Sev::inform)
5283  << "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
5284  CHKERR VecRestoreArrayRead(norms_vec, &norms);
5285 
5286  MoFEMFunctionReturn(0);
5287 }

◆ getTractionFreeBc()

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

Remove all, but entities where kinematic constrains are applied.

Parameters
meshset
bc_ptr
disp_block_set_name
rot_block_set_name
contact_set_name
Returns
MoFEMErrorCode
Examples
EshelbianPlasticity.cpp.

Definition at line 1796 of file EshelbianPlasticity.cpp.

1798  {
1799  MoFEMFunctionBegin;
1800 
1801  // get skin from all tets
1802  Range tets;
1803  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
1804  Range tets_skin_part;
1805  Skinner skin(&mField.get_moab());
1806  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
1807  ParallelComm *pcomm =
1808  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1809  Range tets_skin;
1810  CHKERR pcomm->filter_pstatus(tets_skin_part,
1811  PSTATUS_SHARED | PSTATUS_MULTISHARED,
1812  PSTATUS_NOT, -1, &tets_skin);
1813 
1814  bc_ptr->resize(3);
1815  for (int dd = 0; dd != 3; ++dd)
1816  (*bc_ptr)[dd] = tets_skin;
1817 
1818  // Do not remove dofs on which traction is applied
1819  if (bcSpatialDispVecPtr)
1820  for (auto &v : *bcSpatialDispVecPtr) {
1821  if (v.flags[0])
1822  (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
1823  if (v.flags[1])
1824  (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
1825  if (v.flags[2])
1826  (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
1827  }
1828 
1829  // Do not remove dofs on which rotation is applied
1830  if (bcSpatialRotationVecPtr)
1831  for (auto &v : *bcSpatialRotationVecPtr) {
1832  (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
1833  (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
1834  (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
1835  }
1836 
1837  if (bcSpatialTraction)
1838  for (auto &v : *bcSpatialTraction) {
1839  (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
1840  (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
1841  (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
1842  }
1843 
1844  // remove contact
1845  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
1846  std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
1847  Range faces;
1848  CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
1849  true);
1850  (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
1851  (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
1852  (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
1853  }
1854 
1855  MoFEMFunctionReturn(0);
1856 }

◆ inv_d_f_linear()

static double EshelbianCore::inv_d_f_linear ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 72 of file EshelbianCore.hpp.

72 { return 0; }

◆ inv_d_f_log()

static double EshelbianCore::inv_d_f_log ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 60 of file EshelbianCore.hpp.

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

◆ inv_d_f_log_e()

static double EshelbianCore::inv_d_f_log_e ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 43 of file EshelbianCore.hpp.

43 { return (1. / v); }

◆ inv_dd_f_linear()

static double EshelbianCore::inv_dd_f_linear ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 73 of file EshelbianCore.hpp.

73 { return 0; }

◆ inv_dd_f_log()

static double EshelbianCore::inv_dd_f_log ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 63 of file EshelbianCore.hpp.

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

◆ inv_dd_f_log_e()

static double EshelbianCore::inv_dd_f_log_e ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 44 of file EshelbianCore.hpp.

44 { return -(1. / v / v); }

◆ inv_f_linear()

static double EshelbianCore::inv_f_linear ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 71 of file EshelbianCore.hpp.

71 { return v - 1; }

◆ inv_f_log()

static double EshelbianCore::inv_f_log ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 57 of file EshelbianCore.hpp.

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

◆ inv_f_log_e()

static double EshelbianCore::inv_f_log_e ( const double  v)
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 42 of file EshelbianCore.hpp.

42 { return std::log(v); }

◆ postProcessResults()

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

Definition at line 3148 of file EshelbianPlasticity.cpp.

3150  {
3151  MoFEMFunctionBegin;
3152 
3153  // mark crack surface
3154  if (crackingOn) {
3155  Tag th;
3156  rval = mField.get_moab().tag_get_handle("CRACK", th);
3157  if (rval == MB_SUCCESS) {
3158  CHKERR mField.get_moab().tag_delete(th);
3159  }
3160  int def_val[] = {0};
3161  CHKERR mField.get_moab().tag_get_handle(
3162  "CRACK", 1, MB_TYPE_INTEGER, th, MB_TAG_SPARSE | MB_TAG_CREAT, def_val);
3163  int mark[] = {1};
3164  CHKERR mField.get_moab().tag_clear_data(th, *crackFaces, mark);
3165  tags_to_transfer.push_back(th);
3166  }
3167 
3168  // add tags to transfer
3169  for (auto t : listTagsToTransfer) {
3170  tags_to_transfer.push_back(t);
3171  }
3172 
3173  if (!dataAtPts) {
3174  dataAtPts =
3175  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
3176  }
3177 
3178  CHKERR DMoFEMLoopFiniteElements(dM, contactElement, contactTreeRhs);
3179 
3180  auto get_post_proc = [&](auto &post_proc_mesh, auto sense) {
3181  using FaceEle = FaceElementForcesAndSourcesCore;
3182  auto post_proc_ptr =
3183  boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
3184  mField, post_proc_mesh);
3185  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
3186  post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
3187  frontAdjEdges);
3188 
3189  auto domain_ops = [&](auto &fe, int sense) {
3190  MoFEMFunctionBegin;
3191  fe.getUserPolynomialBase() =
3192  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
3193  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3194  fe.getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions,
3195  frontAdjEdges);
3196  auto piola_scale_ptr = boost::make_shared<double>(1.0);
3197  fe.getOpPtrVector().push_back(physicalEquations->returnOpSetScale(
3198  piola_scale_ptr, physicalEquations));
3199  fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
3200  piolaStress, dataAtPts->getApproxPAtPts(), piola_scale_ptr));
3201  fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
3202  bubbleField, dataAtPts->getApproxPAtPts(), piola_scale_ptr,
3203  SmartPetscObj<Vec>(), MBMAXTYPE));
3204  if (noStretch) {
3205  fe.getOpPtrVector().push_back(
3206  physicalEquations->returnOpCalculateStretchFromStress(
3207  dataAtPts, physicalEquations));
3208  } else {
3209  fe.getOpPtrVector().push_back(
3210  new OpCalculateTensor2SymmetricFieldValues<3>(
3211  stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
3212  }
3213  if (var_vector) {
3214  auto vec = SmartPetscObj<Vec>(var_vector, true);
3215  fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
3216  piolaStress, dataAtPts->getVarPiolaPts(),
3217  boost::make_shared<double>(1), vec));
3218  fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
3219  bubbleField, dataAtPts->getVarPiolaPts(),
3220  boost::make_shared<double>(1), vec, MBMAXTYPE));
3221  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3222  rotAxis, dataAtPts->getVarRotAxisPts(), vec, MBTET));
3223  if (noStretch) {
3224  fe.getOpPtrVector().push_back(
3225  physicalEquations->returnOpCalculateVarStretchFromStress(
3226  dataAtPts, physicalEquations));
3227  } else {
3228  fe.getOpPtrVector().push_back(
3229  new OpCalculateTensor2SymmetricFieldValues<3>(
3230  stretchTensor, dataAtPts->getVarLogStreachPts(), vec, MBTET));
3231  }
3232  }
3233 
3234  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3235  rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
3236  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3237  rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
3238 
3239  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3240  spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
3241  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3242  spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
3243  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
3244  spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
3245  // evaluate derived quantities
3246  fe.getOpPtrVector().push_back(
3247  new OpCalculateRotationAndSpatialGradient(dataAtPts));
3248 
3249  // evaluate integration points
3250  fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
3251  tag, true, false, dataAtPts, physicalEquations));
3252  if (auto op =
3253  physicalEquations->returnOpCalculateEnergy(dataAtPts, nullptr)) {
3254  fe.getOpPtrVector().push_back(op);
3255  fe.getOpPtrVector().push_back(new OpCalculateEshelbyStress(dataAtPts));
3256  }
3257 
3258  // // post-proc
3259  using OpPPMap = OpPostProcMapInMoab<
3260  SPACE_DIM, SPACE_DIM,
3261  VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator>;
3262 
3263  struct OpSidePPMap : public OpPPMap {
3264  OpSidePPMap(moab::Interface &post_proc_mesh,
3265  std::vector<EntityHandle> &map_gauss_pts,
3266  DataMapVec data_map_scalar, DataMapMat data_map_vec,
3267  DataMapMat data_map_mat, DataMapMat data_symm_map_mat,
3268  int sense)
3269  : OpPPMap(post_proc_mesh, map_gauss_pts, data_map_scalar,
3270  data_map_vec, data_map_mat, data_symm_map_mat),
3271  tagSense(sense) {}
3272 
3273  MoFEMErrorCode doWork(int side, EntityType type,
3274  EntitiesFieldData::EntData &data) {
3275  MoFEMFunctionBegin;
3276 
3277  if (tagSense != OpPPMap::getSkeletonSense())
3278  MoFEMFunctionReturnHot(0);
3279 
3280  CHKERR OpPPMap::doWork(side, type, data);
3281  MoFEMFunctionReturn(0);
3282  }
3283 
3284  private:
3285  int tagSense;
3286  };
3287 
3288  OpPPMap::DataMapMat vec_fields;
3289  vec_fields["SpatialDisplacementL2"] = dataAtPts->getSmallWL2AtPts();
3290  vec_fields["SpatialDisplacementH1"] = dataAtPts->getSmallWH1AtPts();
3291  vec_fields["Omega"] = dataAtPts->getRotAxisAtPts();
3292  vec_fields["ContactDisplacement"] = dataAtPts->getContactL2AtPts();
3293  vec_fields["AngularMomentum"] = dataAtPts->getLeviKirchhoffAtPts();
3294  vec_fields["X"] = dataAtPts->getLargeXH1AtPts();
3295  if (var_vector) {
3296  auto vec = SmartPetscObj<Vec>(var_vector, true);
3297  vec_fields["VarOmega"] = dataAtPts->getVarRotAxisPts();
3298  vec_fields["VarSpatialDisplacementL2"] =
3299  boost::make_shared<MatrixDouble>();
3300  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3301  spatialL2Disp, vec_fields["VarSpatialDisplacementL2"], vec, MBTET));
3302  }
3303  if (f_residual) {
3304  auto vec = SmartPetscObj<Vec>(f_residual, true);
3305  vec_fields["ResSpatialDisplacementL2"] =
3306  boost::make_shared<MatrixDouble>();
3307  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3308  spatialL2Disp, vec_fields["ResSpatialDisplacementL2"], vec, MBTET));
3309  vec_fields["ResOmega"] = boost::make_shared<MatrixDouble>();
3310  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3311  rotAxis, vec_fields["ResOmega"], vec, MBTET));
3312  }
3313 
3314  OpPPMap::DataMapMat mat_fields;
3315  mat_fields["PiolaStress"] = dataAtPts->getApproxPAtPts();
3316  if (var_vector) {
3317  mat_fields["VarPiolaStress"] = dataAtPts->getVarPiolaPts();
3318  }
3319  if (f_residual) {
3320  auto vec = SmartPetscObj<Vec>(f_residual, true);
3321  mat_fields["ResPiolaStress"] = boost::make_shared<MatrixDouble>();
3322  fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
3323  piolaStress, mat_fields["ResPiolaStress"],
3324  boost::make_shared<double>(1), vec));
3325  fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
3326  bubbleField, mat_fields["ResPiolaStress"],
3327  boost::make_shared<double>(1), vec, MBMAXTYPE));
3328  }
3329 
3330  OpPPMap::DataMapMat mat_fields_symm;
3331  mat_fields_symm["LogSpatialStretch"] =
3332  dataAtPts->getLogStretchTensorAtPts();
3333  mat_fields_symm["SpatialStretch"] = dataAtPts->getStretchTensorAtPts();
3334  if (var_vector) {
3335  mat_fields_symm["VarLogSpatialStretch"] =
3336  dataAtPts->getVarLogStreachPts();
3337  }
3338  if (f_residual) {
3339  auto vec = SmartPetscObj<Vec>(f_residual, true);
3340  if (!noStretch) {
3341  mat_fields_symm["ResLogSpatialStretch"] =
3342  boost::make_shared<MatrixDouble>();
3343  fe.getOpPtrVector().push_back(
3344  new OpCalculateTensor2SymmetricFieldValues<3>(
3345  stretchTensor, mat_fields_symm["ResLogSpatialStretch"], vec,
3346  MBTET));
3347  }
3348  }
3349 
3350  fe.getOpPtrVector().push_back(
3351 
3352  new OpSidePPMap(
3353 
3354  post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3355 
3356  {},
3357 
3358  vec_fields,
3359 
3360  mat_fields,
3361 
3362  mat_fields_symm,
3363 
3364  sense
3365 
3366  )
3367 
3368  );
3369 
3370  fe.getOpPtrVector().push_back(new OpPostProcDataStructure(
3371  post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3372  dataAtPts, sense));
3373 
3374  MoFEMFunctionReturn(0);
3375  };
3376 
3377  post_proc_ptr->getOpPtrVector().push_back(
3378  new OpCalculateVectorFieldValues<3>(contactDisp,
3379  dataAtPts->getContactL2AtPts()));
3380  auto X_h1_ptr = boost::make_shared<MatrixDouble>();
3381  // H1 material positions
3382  post_proc_ptr->getOpPtrVector().push_back(
3383  new OpCalculateVectorFieldValues<3>(materialH1Positions,
3384  dataAtPts->getLargeXH1AtPts()));
3385 
3386  // domain
3387  auto op_loop_side = new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
3388  mField, elementVolumeName, SPACE_DIM);
3389  domain_ops(*(op_loop_side->getSideFEPtr()), sense);
3390  post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
3391 
3392  return post_proc_ptr;
3393  };
3394 
3395  // contact
3396  auto calcs_side_traction_and_displacements = [&](auto &post_proc_ptr,
3397  auto &pip) {
3398  MoFEMFunctionBegin;
3399  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
3400  // evaluate traction
3401  using EleOnSide =
3402  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
3403  using SideEleOp = EleOnSide::UserDataOperator;
3404  auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
3405  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
3406  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
3407  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
3408  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
3409  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
3410  materialH1Positions, frontAdjEdges);
3411  op_loop_domain_side->getOpPtrVector().push_back(
3412  new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
3413  piolaStress, contact_common_data_ptr->contactTractionPtr(),
3414  boost::make_shared<double>(1.0)));
3415  pip.push_back(op_loop_domain_side);
3416  // evaluate contact displacement and contact conditions
3417  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
3418  pip.push_back(new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
3419  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
3420  contactDisp, contact_common_data_ptr->contactDispPtr()));
3421  pip.push_back(new OpTreeSearch(
3422  contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
3423  get_range_from_block(mField, "CONTACT", SPACE_DIM - 1),
3424  &post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
3425  MoFEMFunctionReturn(0);
3426  };
3427 
3428  auto post_proc_mesh = boost::make_shared<moab::Core>();
3429  auto post_proc_ptr = get_post_proc(post_proc_mesh, 1);
3430  auto post_proc_negative_sense_ptr = get_post_proc(post_proc_mesh, -1);
3431  CHKERR calcs_side_traction_and_displacements(post_proc_ptr,
3432  post_proc_ptr->getOpPtrVector());
3433 
3434  auto own_tets =
3435  CommInterface::getPartEntities(mField.get_moab(), mField.get_comm_rank())
3436  .subset_by_dimension(SPACE_DIM);
3437  Range own_faces;
3438  CHKERR mField.get_moab().get_adjacencies(own_tets, SPACE_DIM - 1, true,
3439  own_faces, moab::Interface::UNION);
3440 
3441  auto get_post_negative = [&](auto &&ents) {
3442  auto crack_faces_pos = ents;
3443  auto crack_faces_neg = crack_faces_pos;
3444  auto skin = get_skin(mField, own_tets);
3445  auto crack_on_proc_skin = intersect(crack_faces_pos, skin);
3446  for (auto f : crack_on_proc_skin) {
3447  Range tet;
3448  CHKERR mField.get_moab().get_adjacencies(&f, 1, SPACE_DIM, false, tet);
3449  tet = intersect(tet, own_tets);
3450  int side_number, sense, offset;
3451  CHKERR mField.get_moab().side_number(tet[0], f, side_number, sense,
3452  offset);
3453  if (sense == 1) {
3454  crack_faces_neg.erase(f);
3455  } else {
3456  crack_faces_pos.erase(f);
3457  }
3458  }
3459  return std::make_pair(crack_faces_pos, crack_faces_neg);
3460  };
3461 
3462  auto get_crack_faces = [&](auto crack_faces) {
3463  auto get_adj = [&](auto e, auto dim) {
3464  Range adj;
3465  CHKERR mField.get_moab().get_adjacencies(e, dim, true, adj,
3466  moab::Interface::UNION);
3467  return adj;
3468  };
3469  auto tets = get_adj(crack_faces, 3);
3470  auto faces = subtract(get_adj(tets, 2), crack_faces);
3471  tets = subtract(tets, get_adj(faces, 3));
3472  return subtract(crack_faces, get_adj(tets, 2));
3473  };
3474 
3475  auto [crack_faces_pos, crack_faces_neg] =
3476  get_post_negative(intersect(own_faces, get_crack_faces(*crackFaces)));
3477 
3478  auto only_crack_faces = [&](FEMethod *fe_method_ptr) {
3479  auto ent = fe_method_ptr->getFEEntityHandle();
3480  if (crack_faces_pos.find(ent) == crack_faces_pos.end()) {
3481  return false;
3482  }
3483  return true;
3484  };
3485 
3486  auto only_negative_crack_faces = [&](FEMethod *fe_method_ptr) {
3487  auto ent = fe_method_ptr->getFEEntityHandle();
3488  if (crack_faces_neg.find(ent) == crack_faces_neg.end()) {
3489  return false;
3490  }
3491  return true;
3492  };
3493 
3494  auto get_adj_front = [&]() {
3495  auto skeleton_faces = *skeletonFaces;
3496  Range adj_front;
3497  CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2, true, adj_front,
3498  moab::Interface::UNION);
3499 
3500  adj_front = intersect(adj_front, skeleton_faces);
3501  adj_front = subtract(adj_front, *crackFaces);
3502  adj_front = intersect(own_faces, adj_front);
3503  return adj_front;
3504  };
3505 
3506  post_proc_ptr->setTagsToTransfer(tags_to_transfer);
3507  post_proc_negative_sense_ptr->setTagsToTransfer(tags_to_transfer);
3508 
3509  auto post_proc_begin =
3510  PostProcBrokenMeshInMoabBaseBegin(mField, post_proc_mesh);
3511  CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
3512  CHKERR DMoFEMLoopFiniteElements(dM, skinElement, post_proc_ptr);
3513  post_proc_ptr->exeTestHook = only_crack_faces;
3514  post_proc_negative_sense_ptr->exeTestHook = only_negative_crack_faces;
3515  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
3516  dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
3517  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, skeletonElement,
3518  post_proc_negative_sense_ptr, 0,
3519  mField.get_comm_size());
3520 
3521  constexpr bool debug = false;
3522  if (debug) {
3523 
3524  auto [adj_front_pos, adj_front_neg] =
3525  get_post_negative(filter_owners(mField, get_adj_front()));
3526 
3527  auto only_front_faces_pos = [adj_front_pos](FEMethod *fe_method_ptr) {
3528  auto ent = fe_method_ptr->getFEEntityHandle();
3529  if (adj_front_pos.find(ent) == adj_front_pos.end()) {
3530  return false;
3531  }
3532  return true;
3533  };
3534 
3535  auto only_front_faces_neg = [adj_front_neg](FEMethod *fe_method_ptr) {
3536  auto ent = fe_method_ptr->getFEEntityHandle();
3537  if (adj_front_neg.find(ent) == adj_front_neg.end()) {
3538  return false;
3539  }
3540  return true;
3541  };
3542 
3543  post_proc_ptr->exeTestHook = only_front_faces_pos;
3544  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
3545  dM, skeletonElement, post_proc_ptr, 0, mField.get_comm_size());
3546  post_proc_negative_sense_ptr->exeTestHook = only_front_faces_neg;
3547  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, skeletonElement,
3548  post_proc_negative_sense_ptr, 0,
3549  mField.get_comm_size());
3550  }
3551  auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
3552  CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
3553 
3554  CHKERR post_proc_end.writeFile(file.c_str());
3555  MoFEMFunctionReturn(0);
3556 }

◆ postProcessSkeletonResults()

MoFEMErrorCode EshelbianCore::postProcessSkeletonResults ( const int  tag,
const std::string  file,
Vec  f_residual = PETSC_NULL 
)
Examples
EshelbianPlasticity.cpp.

Definition at line 3558 of file EshelbianPlasticity.cpp.

3560  {
3561  MoFEMFunctionBegin;
3562 
3563  using FaceEle = FaceElementForcesAndSourcesCore;
3564 
3565  auto post_proc_mesh = boost::make_shared<moab::Core>();
3566  auto post_proc_ptr =
3567  boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
3568  mField, post_proc_mesh);
3569  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
3570  post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
3571  frontAdjEdges);
3572 
3573  auto hybrid_disp = boost::make_shared<MatrixDouble>();
3574  post_proc_ptr->getOpPtrVector().push_back(
3575  new OpCalculateVectorFieldValues<3>(hybridSpatialDisp, hybrid_disp));
3576 
3577  auto hybrid_res = boost::make_shared<MatrixDouble>();
3578  post_proc_ptr->getOpPtrVector().push_back(
3579  new OpCalculateVectorFieldValues<3>(hybridSpatialDisp, hybrid_disp));
3580 
3581  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
3582  post_proc_ptr->getOpPtrVector().push_back(
3583 
3584  new OpPPMap(
3585 
3586  post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3587 
3588  {},
3589 
3590  {{"hybrid_disp", hybrid_disp}},
3591 
3592  {},
3593 
3594  {}
3595 
3596  )
3597 
3598  );
3599 
3600  if (f_residual) {
3601  auto hybrid_res = boost::make_shared<MatrixDouble>();
3602  post_proc_ptr->getOpPtrVector().push_back(
3603  new OpCalculateVectorFieldValues<3>(
3604  hybridSpatialDisp, hybrid_res,
3605  SmartPetscObj<Vec>(f_residual, true)));
3606  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
3607  post_proc_ptr->getOpPtrVector().push_back(
3608 
3609  new OpPPMap(
3610 
3611  post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
3612 
3613  {},
3614 
3615  {{"res_hybrid", hybrid_res}},
3616 
3617  {},
3618 
3619  {}
3620 
3621  )
3622 
3623  );
3624  }
3625 
3626  auto post_proc_begin =
3627  PostProcBrokenMeshInMoabBaseBegin(mField, post_proc_mesh);
3628  CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
3629  CHKERR DMoFEMLoopFiniteElements(dM, skeletonElement, post_proc_ptr);
3630  auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
3631  CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
3632 
3633  CHKERR post_proc_end.writeFile(file.c_str());
3634 
3635  MoFEMFunctionReturn(0);
3636 }

◆ projectGeometry()

MoFEMErrorCode EshelbianCore::projectGeometry ( const EntityHandle  meshset = 0)
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 1251 of file EshelbianPlasticity.cpp.

1251  {
1252  MoFEMFunctionBegin;
1253 
1254  Range meshset_ents;
1255  CHKERR mField.get_moab().get_entities_by_handle(meshset, meshset_ents);
1256 
1257  auto project_ho_geometry = [&](auto field) {
1258  Projection10NodeCoordsOnField ent_method(mField, materialH1Positions);
1259  return mField.loop_dofs(field, ent_method);
1260  };
1261  CHKERR project_ho_geometry(materialH1Positions);
1262 
1263  auto get_adj_front_edges = [&](auto &front_edges) {
1264  auto &moab = mField.get_moab();
1265  Range front_crack_nodes;
1266  CHKERR moab.get_connectivity(front_edges, front_crack_nodes, true);
1267  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1268  front_crack_nodes);
1269  Range crack_front_edges;
1270  CHKERR moab.get_adjacencies(front_crack_nodes, SPACE_DIM - 2, false,
1271  crack_front_edges, moab::Interface::UNION);
1272  crack_front_edges =
1273  intersect(subtract(crack_front_edges, front_edges), meshset_ents);
1274  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1275  crack_front_edges);
1276 
1277  Range crack_front_edges_nodes;
1278  CHKERR moab.get_connectivity(crack_front_edges, crack_front_edges_nodes,
1279  true);
1280  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1281  crack_front_edges_nodes);
1282  crack_front_edges_nodes =
1283  subtract(crack_front_edges_nodes, front_crack_nodes);
1284  Range crack_front_edges_with_both_nodes_not_at_front;
1285  CHKERR moab.get_adjacencies(crack_front_edges_nodes, SPACE_DIM - 2, false,
1286  crack_front_edges_with_both_nodes_not_at_front,
1287  moab::Interface::UNION);
1288  crack_front_edges_with_both_nodes_not_at_front =
1289  intersect(crack_front_edges_with_both_nodes_not_at_front, meshset_ents);
1290  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1291  crack_front_edges_with_both_nodes_not_at_front);
1292  crack_front_edges_with_both_nodes_not_at_front = intersect(
1293  crack_front_edges_with_both_nodes_not_at_front, crack_front_edges);
1294 
1295  return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
1296  boost::make_shared<Range>(
1297  crack_front_edges_with_both_nodes_not_at_front));
1298  };
1299 
1300 
1301  crackFaces = boost::make_shared<Range>(
1302  get_range_from_block(mField, "CRACK", SPACE_DIM - 1));
1303  frontEdges =
1304  boost::make_shared<Range>(get_crack_front_edges(mField, *crackFaces));
1305  auto [front_vertices, front_adj_edges] = get_adj_front_edges(*frontEdges);
1306  frontVertices = front_vertices;
1307  frontAdjEdges = front_adj_edges;
1308 
1309  // #ifndef NDEBUG
1310  if (crackingOn) {
1311  auto rank = mField.get_comm_rank();
1312  // CHKERR save_range(mField.get_moab(),
1313  // (boost::format("meshset_ents_%d.vtk") % rank).str(),
1314  // meshset_ents);
1315  CHKERR save_range(mField.get_moab(),
1316  (boost::format("crack_faces_%d.vtk") % rank).str(),
1317  *crackFaces);
1318  // CHKERR save_range(mField.get_moab(),
1319  // (boost::format("front_edges_%d.vtk") % rank).str(),
1320  // *frontEdges);
1321  // CHKERR save_range(mField.get_moab(),
1322  // (boost::format("front_vertices_%d.vtk") % rank).str(),
1323  // *frontVertices);
1324  // CHKERR save_range(mField.get_moab(),
1325  // (boost::format("front_adj_edges_%d.vtk") % rank).str(),
1326  // *frontAdjEdges);
1327  }
1328  // #endif // NDEBUG
1329 
1330  auto set_singular_dofs = [&](auto &front_adj_edges, auto &front_vertices) {
1331  MoFEMFunctionBegin;
1332  auto &moab = mField.get_moab();
1333 
1334  double eps = 1;
1335  double beta = 0;
1336  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-singularity_eps", &beta,
1337  PETSC_NULL);
1338  MOFEM_LOG("EP", Sev::inform) << "Singularity eps " << beta;
1339  eps -= beta;
1340 
1341  for (auto edge : front_adj_edges) {
1342  int num_nodes;
1343  const EntityHandle *conn;
1344  CHKERR moab.get_connectivity(edge, conn, num_nodes, false);
1345  double coords[6];
1346  CHKERR moab.get_coords(conn, num_nodes, coords);
1347  const double dir[3] = {coords[3] - coords[0], coords[4] - coords[1],
1348  coords[5] - coords[2]};
1349  double dof[3] = {0, 0, 0};
1350  if (front_vertices.find(conn[0]) != front_vertices.end()) {
1351  for (int dd = 0; dd != 3; dd++) {
1352  dof[dd] = -dir[dd] * eps;
1353  }
1354  } else if (front_vertices.find(conn[1]) != front_vertices.end()) {
1355  for (int dd = 0; dd != 3; dd++) {
1356  dof[dd] = +dir[dd] * eps;
1357  }
1358  }
1359  for (_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(
1360  mField, materialH1Positions, edge, dit)) {
1361  const int idx = dit->get()->getEntDofIdx();
1362  if (idx > 2) {
1363  dit->get()->getFieldData() = 0;
1364  } else {
1365  dit->get()->getFieldData() = dof[idx];
1366  }
1367  }
1368  }
1369 
1370  MoFEMFunctionReturn(0);
1371  };
1372 
1373  if (setSingularity)
1374  CHKERR set_singular_dofs(*frontAdjEdges, *frontVertices);
1375 
1376  MoFEMFunctionReturn(0);
1377 }

◆ 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
EshelbianPlasticity.cpp.

Definition at line 859 of file EshelbianPlasticity.cpp.

860  {
861  *iface = const_cast<EshelbianCore *>(this);
862  return 0;
863 }

◆ saveOrgCoords()

MoFEMErrorCode EshelbianCore::saveOrgCoords ( )
Examples
EshelbianPlasticity.cpp.

Definition at line 5194 of file EshelbianPlasticity.cpp.

5194  {
5195  MoFEMFunctionBegin;
5196  auto crack_faces =
5197  get_range_from_block(mField, "CRACK_COMPUTED", SPACE_DIM - 1);
5198  Range conn;
5199  CHKERR mField.get_moab().get_connectivity(crack_faces, conn, true);
5200  Range verts;
5201  CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, verts);
5202  verts = subtract(verts, conn);
5203  std::vector<double> coords(3 * verts.size());
5204  CHKERR mField.get_moab().get_coords(verts, coords.data());
5205  double def_coords[] = {0., 0., 0.};
5206  Tag th_org_coords;
5207  CHKERR mField.get_moab().tag_get_handle(
5208  "ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
5209  MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
5210  CHKERR mField.get_moab().tag_set_data(th_org_coords, verts, coords.data());
5211  MoFEMFunctionReturn(0);
5212 }

◆ setBaseVolumeElementOps()

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

Definition at line 1999 of file EshelbianPlasticity.cpp.

2002  {
2003  MoFEMFunctionBegin;
2004 
2005  auto bubble_cache =
2006  boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0, MatrixDouble());
2007  fe->getUserPolynomialBase() =
2008  boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
2009  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2010  fe->getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions, frontAdjEdges);
2011 
2012  // set integration rule
2013  fe->getRuleHook = [](int, int, int) { return -1; };
2014  fe->setRuleHook = SetIntegrationAtFrontVolume(frontVertices, frontAdjEdges);
2015  // fe->getRuleHook = VolRule();
2016 
2017  if (!dataAtPts) {
2018  dataAtPts =
2019  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
2020  dataAtPts->physicsPtr = physicalEquations;
2021  }
2022 
2023  // calculate fields values
2024  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
2025  piolaStress, dataAtPts->getApproxPAtPts()));
2026  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
2027  bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
2028  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
2029  piolaStress, dataAtPts->getDivPAtPts()));
2030 
2031  if (noStretch) {
2032  fe->getOpPtrVector().push_back(
2033  physicalEquations->returnOpCalculateStretchFromStress(
2034  dataAtPts, physicalEquations));
2035  } else {
2036  fe->getOpPtrVector().push_back(
2037  new OpCalculateTensor2SymmetricFieldValues<3>(
2038  stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
2039  }
2040 
2041  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2042  rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
2043  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2044  rotAxis, dataAtPts->getRotAxis0AtPts(), solTSStep, MBTET));
2045  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2046  spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
2047 
2048  // H1 displacements
2049  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2050  spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
2051  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
2052  spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
2053 
2054  // velocities
2055  if (calc_rates) {
2056  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
2057  spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
2058  if (noStretch) {
2059  } else {
2060  fe->getOpPtrVector().push_back(
2061  new OpCalculateTensor2SymmetricFieldValuesDot<3>(
2062  stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
2063  }
2064  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
2065  rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
2066  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradientDot<3, 3>(
2067  rotAxis, dataAtPts->getRotAxisGradDotAtPts(), MBTET));
2068 
2069  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradientDot<6, 3>(
2070  stretchTensor, dataAtPts->getGradLogStretchDotTensorAtPts(), MBTET));
2071 
2072  // acceleration
2073  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
2074  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
2075  spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
2076  }
2077  }
2078 
2079  // variations
2080  if (var_vec.use_count()) {
2081  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
2082  piolaStress, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
2083  var_vec));
2084  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
2085  bubbleField, dataAtPts->getVarPiolaPts(), boost::make_shared<double>(1),
2086  var_vec, MBMAXTYPE));
2087  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2088  rotAxis, dataAtPts->getVarRotAxisPts(), var_vec, MBTET));
2089  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
2090  piolaStress, dataAtPts->getDivVarPiolaPts(), var_vec));
2091  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2092  spatialL2Disp, dataAtPts->getVarWL2Pts(), var_vec, MBTET));
2093 
2094  if (noStretch) {
2095  fe->getOpPtrVector().push_back(
2096  physicalEquations->returnOpCalculateVarStretchFromStress(
2097  dataAtPts, physicalEquations));
2098  } else {
2099  fe->getOpPtrVector().push_back(
2100  new OpCalculateTensor2SymmetricFieldValues<3>(
2101  stretchTensor, dataAtPts->getVarLogStreachPts(), var_vec, MBTET));
2102  }
2103  }
2104 
2105  // calculate other derived quantities
2106  fe->getOpPtrVector().push_back(new OpCalculateRotationAndSpatialGradient(
2107  dataAtPts, ((do_rhs || do_lhs) && calc_rates) ? alphaOmega : 0.));
2108 
2109  // evaluate integration points
2110  if (noStretch) {
2111  } else {
2112  fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
2113  tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
2114  }
2115 
2116  MoFEMFunctionReturn(0);
2117 }

◆ setBlockTagsOnSkin()

MoFEMErrorCode EshelbianCore::setBlockTagsOnSkin ( )
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 3087 of file EshelbianPlasticity.cpp.

3087  {
3088  MoFEMFunctionBegin;
3089 
3090  auto set_block = [&](auto name, int dim) {
3091  std::map<int, Range> map;
3092  auto set_tag_impl = [&](auto name) {
3093  MoFEMFunctionBegin;
3094  auto mesh_mng = mField.getInterface<MeshsetsManager>();
3095  auto bcs = mesh_mng->getCubitMeshsetPtr(
3096 
3097  std::regex((boost::format("%s(.*)") % name).str())
3098 
3099  );
3100  for (auto bc : bcs) {
3101  Range r;
3102  CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), dim, r,
3103  true);
3104  map[bc->getMeshsetId()] = r;
3105  }
3106  MoFEMFunctionReturn(0);
3107  };
3108 
3109  CHKERR set_tag_impl(name);
3110 
3111  return std::make_pair(name, map);
3112  };
3113 
3114  auto set_skin = [&](auto &&map) {
3115  for (auto &m : map.second) {
3116  auto s = filter_true_skin(mField, get_skin(mField, m.second));
3117  m.second.swap(s);
3118  }
3119  return map;
3120  };
3121 
3122  auto set_tag = [&](auto &&map) {
3123  Tag th;
3124  auto name = map.first;
3125  int def_val[] = {-1};
3126  CHK_MOAB_THROW(
3127  mField.get_moab().tag_get_handle(name, 1, MB_TYPE_INTEGER, th,
3128  MB_TAG_SPARSE | MB_TAG_CREAT, def_val),
3129  "create tag");
3130  for (auto &m : map.second) {
3131  int id = m.first;
3132  CHK_MOAB_THROW(mField.get_moab().tag_clear_data(th, m.second, &id),
3133  "clear tag");
3134  }
3135  return th;
3136  };
3137 
3138  listTagsToTransfer.push_back(set_tag(set_skin(set_block("BODY", 3))));
3139  listTagsToTransfer.push_back(set_tag(set_skin(set_block("MAT_ELASTIC", 3))));
3140  listTagsToTransfer.push_back(
3141  set_tag(set_skin(set_block("MAT_NEOHOOKEAN", 3))));
3142  listTagsToTransfer.push_back(set_tag(set_block("CONTACT", 2)));
3143 
3144  MoFEMFunctionReturn(0);
3145 }

◆ setContactElementRhsOps()

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

Contact requires that body is marked

Examples
EshelbianPlasticity.cpp.

Definition at line 2653 of file EshelbianPlasticity.cpp.

2657  {
2658  MoFEMFunctionBegin;
2659 
2660  /** Contact requires that body is marked */
2661  auto get_body_range = [this](auto name, int dim, auto sev) {
2662  std::map<int, Range> map;
2663 
2664  for (auto m_ptr :
2665  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
2666 
2667  (boost::format("%s(.*)") % name).str()
2668 
2669  ))
2670 
2671  ) {
2672  Range ents;
2673  CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
2674  dim, ents, true),
2675  "by dim");
2676  map[m_ptr->getMeshsetId()] = ents;
2677  MOFEM_LOG("EPSYNC", sev) << "Meshset: " << m_ptr->getMeshsetId() << " "
2678  << ents.size() << " entities";
2679  }
2680 
2681  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), sev);
2682  return map;
2683  };
2684 
2685  auto get_map_skin = [this](auto &&map) {
2686  ParallelComm *pcomm =
2687  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2688 
2689  Skinner skin(&mField.get_moab());
2690  for (auto &m : map) {
2691  Range skin_faces;
2692  CHKERR skin.find_skin(0, m.second, false, skin_faces);
2693  CHK_MOAB_THROW(pcomm->filter_pstatus(skin_faces,
2694  PSTATUS_SHARED | PSTATUS_MULTISHARED,
2695  PSTATUS_NOT, -1, nullptr),
2696  "filter");
2697  m.second.swap(skin_faces);
2698  }
2699  return map;
2700  };
2701 
2702  /* The above code is written in C++ and it appears to be defining and using
2703  various operations on boundary elements and side elements. */
2704  using BoundaryEle = FaceElementForcesAndSourcesCore;
2705  using BoundaryEleOp = BoundaryEle::UserDataOperator;
2706 
2707  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
2708 
2709  auto calcs_side_traction = [&](auto &pip) {
2710  MoFEMFunctionBegin;
2711  using EleOnSide =
2712  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
2713  using SideEleOp = EleOnSide::UserDataOperator;
2714  auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
2715  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2716  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2717  boost::make_shared<CGGUserPolynomialBase>();
2718  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2719  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2720  materialH1Positions, frontAdjEdges);
2721  op_loop_domain_side->getOpPtrVector().push_back(
2722  new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
2723  piolaStress, contact_common_data_ptr->contactTractionPtr(),
2724  boost::make_shared<double>(1.0)));
2725  pip.push_back(op_loop_domain_side);
2726  MoFEMFunctionReturn(0);
2727  };
2728 
2729  auto add_contact_three = [&]() {
2730  MoFEMFunctionBegin;
2731  auto tree_moab_ptr = boost::make_shared<moab::Core>();
2732  fe_contact_tree = boost::make_shared<ContactTree>(
2733  mField, tree_moab_ptr, spaceOrder,
2734  get_body_range("CONTACT", SPACE_DIM - 1, Sev::inform));
2735  fe_contact_tree->getOpPtrVector().push_back(
2736  new OpCalculateVectorFieldValues<3>(
2737  contactDisp, contact_common_data_ptr->contactDispPtr()));
2738  CHKERR calcs_side_traction(fe_contact_tree->getOpPtrVector());
2739  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2740  fe_contact_tree->getOpPtrVector().push_back(
2741  new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
2742  fe_contact_tree->getOpPtrVector().push_back(
2743  new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
2744  MoFEMFunctionReturn(0);
2745  };
2746 
2747  CHKERR add_contact_three();
2748 
2749  MoFEMFunctionReturn(0);
2750 }

◆ setElasticElementOps()

MoFEMErrorCode EshelbianCore::setElasticElementOps ( const int  tag)
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 2752 of file EshelbianPlasticity.cpp.

2752  {
2753  MoFEMFunctionBegin;
2754 
2755  // Add contact operators. Note that only for rhs. THe lhs is assembled with
2756  // volume element, to enable schur complement evaluation.
2757  CHKERR setContactElementRhsOps(contactTreeRhs);
2758 
2759  CHKERR setVolumeElementOps(tag, true, false, elasticFeRhs, elasticFeLhs);
2760  CHKERR setFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
2761 
2762  auto adj_cache =
2763  boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
2764 
2765  auto get_op_contact_bc = [&]() {
2766  using SideEle = FaceElementForcesAndSourcesCore;
2767  auto op_loop_side = new OpLoopSide<SideEle>(
2768  mField, contactElement, SPACE_DIM - 1, Sev::noisy, adj_cache);
2769  return op_loop_side;
2770  };
2771 
2772  MoFEMFunctionReturn(0);
2773 }

◆ setElasticElementToTs()

MoFEMErrorCode EshelbianCore::setElasticElementToTs ( DM  dm)
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 2775 of file EshelbianPlasticity.cpp.

2775  {
2776  MoFEMFunctionBegin;
2777  boost::shared_ptr<FEMethod> null;
2778 
2779  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
2780 
2781  CHKERR DMMoFEMTSSetI2Function(dm, elementVolumeName, elasticFeRhs, null,
2782  null);
2783  CHKERR DMMoFEMTSSetI2Function(dm, naturalBcElement, elasticBcRhs, null,
2784  null);
2785  CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
2786  null);
2787  // CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
2788  // null);
2789 
2790  } else {
2791  CHKERR DMMoFEMTSSetIFunction(dm, elementVolumeName, elasticFeRhs, null,
2792  null);
2793  CHKERR DMMoFEMTSSetIFunction(dm, naturalBcElement, elasticBcRhs, null,
2794  null);
2795  CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
2796  null);
2797  // CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
2798  // null);
2799  }
2800 
2801  MoFEMFunctionReturn(0);
2802 }

◆ setFaceElementOps()

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

Definition at line 2570 of file EshelbianPlasticity.cpp.

2573  {
2574  MoFEMFunctionBegin;
2575 
2576  fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
2577  fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
2578 
2579  // set integration rule
2580  // fe_rhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
2581  // fe_lhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
2582  fe_rhs->getRuleHook = [](int, int, int) { return -1; };
2583  fe_lhs->getRuleHook = [](int, int, int) { return -1; };
2584  fe_rhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2585  fe_lhs->setRuleHook = SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2586 
2587  CHKERR
2588  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2589  fe_rhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
2590  CHKERR
2591  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2592  fe_lhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
2593 
2594  if (add_elastic) {
2595 
2596  auto get_broken_op_side = [this](auto &pip) {
2597  using EleOnSide =
2598  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
2599  using SideEleOp = EleOnSide::UserDataOperator;
2600  // Iterate over domain FEs adjacent to boundary.
2601  auto broken_data_ptr =
2602  boost::make_shared<std::vector<BrokenBaseSideData>>();
2603  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2604  auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
2605  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2606  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2607  boost::make_shared<CGGUserPolynomialBase>();
2608  CHKERR
2609  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2610  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2611  materialH1Positions, frontAdjEdges);
2612  op_loop_domain_side->getOpPtrVector().push_back(
2613  new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
2614  boost::shared_ptr<double> piola_scale_ptr(new double);
2615  op_loop_domain_side->getOpPtrVector().push_back(
2616  physicalEquations->returnOpSetScale(piola_scale_ptr,
2617  physicalEquations));
2618  pip.push_back(op_loop_domain_side);
2619  return std::make_pair(broken_data_ptr, piola_scale_ptr);
2620  };
2621 
2622  auto [broken_data_ptr, piola_scale_ptr] =
2623  get_broken_op_side(fe_rhs->getOpPtrVector());
2624 
2625  fe_rhs->getOpPtrVector().push_back(new OpDispBc(
2626  broken_data_ptr, bcSpatialDispVecPtr,
2627  {
2628 
2629  boost::make_shared<DynamicRelaxationTimeScale>("disp_history.txt")
2630 
2631  }));
2632  fe_rhs->getOpPtrVector().push_back(
2633  new OpRotationBc(broken_data_ptr, bcSpatialRotationVecPtr,
2634 
2635  {
2636 
2637  boost::make_shared<DynamicRelaxationTimeScale>(
2638  "rotation_history.txt")
2639 
2640  }));
2641 
2642  fe_rhs->getOpPtrVector().push_back(new OpBrokenTractionBc(
2643  hybridSpatialDisp, bcSpatialTraction, piola_scale_ptr,
2644 
2645  {boost::make_shared<DynamicRelaxationTimeScale>("traction_history.txt")}
2646 
2647  ));
2648  }
2649 
2650  MoFEMFunctionReturn(0);
2651 }

◆ setNewFrontCoordinates()

MoFEMErrorCode EshelbianCore::setNewFrontCoordinates ( )
Examples
EshelbianPlasticity.cpp.

Definition at line 4990 of file EshelbianPlasticity.cpp.

4990  {
4991  MoFEMFunctionBegin;
4992 
4993  if (!maxMovedFaces)
4994  MoFEMFunctionReturnHot(0);
4995 
4996  Tag th_front_position;
4997  auto rval =
4998  mField.get_moab().tag_get_handle("FrontPosition", th_front_position);
4999  if (rval == MB_SUCCESS && maxMovedFaces) {
5000  Range verts;
5001  CHKERR mField.get_moab().get_connectivity(*maxMovedFaces, verts, true);
5002  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
5003  std::vector<double> coords(3 * verts.size());
5004  CHKERR mField.get_moab().get_coords(verts, coords.data());
5005  std::vector<double> pos(3 * verts.size());
5006  CHKERR mField.get_moab().tag_get_data(th_front_position, verts, pos.data());
5007  for (int i = 0; i != 3 * verts.size(); ++i) {
5008  coords[i] += pos[i];
5009  }
5010  CHKERR mField.get_moab().set_coords(verts, coords.data());
5011  double zero[] = {0., 0., 0.};
5012  CHKERR mField.get_moab().tag_clear_data(th_front_position, verts, zero);
5013  }
5014 
5015  constexpr bool debug = true;
5016  if (debug) {
5017 
5018  CHKERR save_range(
5019  mField.get_moab(),
5020  "set_coords_faces_" +
5021  boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
5022  *maxMovedFaces);
5023  }
5024 
5025  MoFEMFunctionReturn(0);
5026 }

◆ setVolumeElementOps()

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

Contact requires that body is marked

Examples
EshelbianPlasticity.cpp.

Definition at line 2119 of file EshelbianPlasticity.cpp.

2122  {
2123  MoFEMFunctionBegin;
2124 
2125  /** Contact requires that body is marked */
2126  auto get_body_range = [this](auto name, int dim) {
2127  std::map<int, Range> map;
2128 
2129  for (auto m_ptr :
2130  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
2131 
2132  (boost::format("%s(.*)") % name).str()
2133 
2134  ))
2135 
2136  ) {
2137  Range ents;
2138  CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
2139  dim, ents, true),
2140  "by dim");
2141  map[m_ptr->getMeshsetId()] = ents;
2142  }
2143 
2144  return map;
2145  };
2146 
2147  auto rule_contact = [](int, int, int o) { return -1; };
2148  auto refine = Tools::refineTriangle(contactRefinementLevels);
2149 
2150  auto set_rule_contact = [refine](
2151 
2152  ForcesAndSourcesCore *fe_raw_ptr, int order_row,
2153  int order_col, int order_data
2154 
2155  ) {
2156  MoFEMFunctionBegin;
2157  auto rule = 2 * order_data;
2158  fe_raw_ptr->gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
2159  MoFEMFunctionReturn(0);
2160  };
2161 
2162  auto time_scale = boost::make_shared<DynamicRelaxationTimeScale>();
2163 
2164  // Right hand side
2165  fe_rhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
2166  CHKERR setBaseVolumeElementOps(tag, true, false, true, SmartPetscObj<Vec>(),
2167  fe_rhs);
2168 
2169  // elastic
2170  if (add_elastic) {
2171 
2172  fe_rhs->getOpPtrVector().push_back(
2173  new OpSpatialEquilibrium(spatialL2Disp, dataAtPts, alphaW, alphaRho));
2174  fe_rhs->getOpPtrVector().push_back(
2175  new OpSpatialRotation(rotAxis, dataAtPts, alphaOmega));
2176  if (noStretch) {
2177  // do nothing - no stretch approximation
2178  } else {
2179  fe_rhs->getOpPtrVector().push_back(
2180  physicalEquations->returnOpSpatialPhysical(stretchTensor, dataAtPts,
2181  alphaU));
2182  }
2183  fe_rhs->getOpPtrVector().push_back(
2184  new OpSpatialConsistencyP(piolaStress, dataAtPts));
2185  fe_rhs->getOpPtrVector().push_back(
2186  new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
2187  fe_rhs->getOpPtrVector().push_back(
2188  new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
2189 
2190  auto set_hybridisation = [&](auto &pip) {
2191  MoFEMFunctionBegin;
2192 
2193  using BoundaryEle =
2194  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
2195  using EleOnSide =
2196  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
2197  using SideEleOp = EleOnSide::UserDataOperator;
2198  using BdyEleOp = BoundaryEle::UserDataOperator;
2199 
2200  // First: Iterate over skeleton FEs adjacent to Domain FEs
2201  // Note: BoundaryEle, i.e. uses skeleton interation rule
2202  auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
2203  mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
2204  // op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
2205  op_loop_skeleton_side->getSideFEPtr()->getRuleHook = [](int, int, int) {
2206  return -1;
2207  };
2208  op_loop_skeleton_side->getSideFEPtr()->setRuleHook =
2209  SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
2210 
2211  CHKERR EshelbianPlasticity::
2212  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2213  op_loop_skeleton_side->getOpPtrVector(), {L2},
2214  materialH1Positions, frontAdjEdges);
2215 
2216  // Second: Iterate over domain FEs adjacent to skelton, particularly one
2217  // domain element.
2218  auto broken_data_ptr =
2219  boost::make_shared<std::vector<BrokenBaseSideData>>();
2220  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2221  auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
2222  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2223  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2224  boost::make_shared<CGGUserPolynomialBase>();
2225  CHKERR
2226  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2227  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2228  materialH1Positions, frontAdjEdges);
2229  op_loop_domain_side->getOpPtrVector().push_back(
2230  new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
2231  auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
2232  op_loop_domain_side->getOpPtrVector().push_back(
2233  new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(piolaStress,
2234  flux_mat_ptr));
2235  op_loop_domain_side->getOpPtrVector().push_back(
2236  new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
2237 
2238  // Assemble on skeleton
2239  op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2240  using OpC_dHybrid = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
2241  GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
2242  using OpC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
2243  GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
2244  op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dHybrid(
2245  hybridSpatialDisp, broken_data_ptr, boost::make_shared<double>(1.0)));
2246  auto hybrid_ptr = boost::make_shared<MatrixDouble>();
2247  op_loop_skeleton_side->getOpPtrVector().push_back(
2248  new OpCalculateVectorFieldValues<SPACE_DIM>(hybridSpatialDisp,
2249  hybrid_ptr));
2250  op_loop_skeleton_side->getOpPtrVector().push_back(new OpC_dBroken(
2251  broken_data_ptr, hybrid_ptr, boost::make_shared<double>(1.0)));
2252 
2253  // Add skeleton to domain pipeline
2254  pip.push_back(op_loop_skeleton_side);
2255 
2256  MoFEMFunctionReturn(0);
2257  };
2258 
2259  auto set_contact = [&](auto &pip) {
2260  MoFEMFunctionBegin;
2261 
2262  using BoundaryEle =
2263  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
2264  using EleOnSide =
2265  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
2266  using SideEleOp = EleOnSide::UserDataOperator;
2267  using BdyEleOp = BoundaryEle::UserDataOperator;
2268 
2269  // First: Iterate over skeleton FEs adjacent to Domain FEs
2270  // Note: BoundaryEle, i.e. uses skeleton interation rule
2271  auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
2272  mField, contactElement, SPACE_DIM - 1, Sev::noisy);
2273 
2274  op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
2275  op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
2276  CHKERR EshelbianPlasticity::
2277  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2278  op_loop_skeleton_side->getOpPtrVector(), {L2},
2279  materialH1Positions, frontAdjEdges);
2280 
2281  // Second: Iterate over domain FEs adjacent to skelton, particularly
2282  // one domain element.
2283  auto broken_data_ptr =
2284  boost::make_shared<std::vector<BrokenBaseSideData>>();
2285 
2286  // Data storing contact fields
2287  auto contact_common_data_ptr =
2288  boost::make_shared<ContactOps::CommonData>();
2289 
2290  auto add_ops_domain_side = [&](auto &pip) {
2291  MoFEMFunctionBegin;
2292  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2293  auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
2294  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2295  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2296  boost::make_shared<CGGUserPolynomialBase>();
2297  CHKERR
2298  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2299  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2300  materialH1Positions, frontAdjEdges);
2301  op_loop_domain_side->getOpPtrVector().push_back(
2302  new OpGetBrokenBaseSideData<SideEleOp>(piolaStress,
2303  broken_data_ptr));
2304  op_loop_domain_side->getOpPtrVector().push_back(
2305  new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
2306  piolaStress, contact_common_data_ptr->contactTractionPtr()));
2307  pip.push_back(op_loop_domain_side);
2308  MoFEMFunctionReturn(0);
2309  };
2310 
2311  auto add_ops_contact_rhs = [&](auto &pip) {
2312  MoFEMFunctionBegin;
2313  // get body id and SDF range
2314  auto contact_sfd_map_range_ptr =
2315  boost::make_shared<std::map<int, Range>>(
2316  get_body_range("CONTACT_SDF", SPACE_DIM - 1));
2317 
2318  pip.push_back(new OpCalculateVectorFieldValues<3>(
2319  contactDisp, contact_common_data_ptr->contactDispPtr()));
2320  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2321  pip.push_back(
2322  new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
2323  pip.push_back(new OpTreeSearch(
2324  contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
2325  get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
2326  nullptr));
2327  pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Rhs(
2328  contactDisp, contact_common_data_ptr, contactTreeRhs,
2329  contact_sfd_map_range_ptr));
2330  pip.push_back(
2331  new EshelbianPlasticity::OpConstrainBoundaryHDivRhs<A, GAUSS>(
2332  broken_data_ptr, contact_common_data_ptr, contactTreeRhs));
2333 
2334  MoFEMFunctionReturn(0);
2335  };
2336 
2337  // push ops to face/side pipeline
2338  CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
2339  CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
2340 
2341  // Add skeleton to domain pipeline
2342  pip.push_back(op_loop_skeleton_side);
2343 
2344  MoFEMFunctionReturn(0);
2345  };
2346 
2347  CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
2348  CHKERR set_contact(fe_rhs->getOpPtrVector());
2349 
2350  // Body forces
2351  using BodyNaturalBC =
2352  NaturalBC<VolumeElementForcesAndSourcesCore::UserDataOperator>::
2353  Assembly<PETSC>::LinearForm<GAUSS>;
2354  using OpBodyForce =
2355  BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
2356  CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
2357  fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {time_scale},
2358  "BODY_FORCE", Sev::inform);
2359  }
2360 
2361  // Left hand side
2362  fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
2363  CHKERR setBaseVolumeElementOps(tag, true, true, true, SmartPetscObj<Vec>(),
2364  fe_lhs);
2365 
2366  // elastic
2367  if (add_elastic) {
2368 
2369  if (noStretch) {
2370  fe_lhs->getOpPtrVector().push_back(
2371  new OpSpatialConsistency_dP_dP(piolaStress, piolaStress, dataAtPts));
2372  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dP(
2373  bubbleField, piolaStress, dataAtPts));
2374  fe_lhs->getOpPtrVector().push_back(
2375  new OpSpatialConsistency_dBubble_dBubble(bubbleField, bubbleField,
2376  dataAtPts));
2377  } else {
2378  fe_lhs->getOpPtrVector().push_back(
2379  physicalEquations->returnOpSpatialPhysical_du_du(
2380  stretchTensor, stretchTensor, dataAtPts, alphaU));
2381  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
2382  stretchTensor, piolaStress, dataAtPts, true));
2383  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
2384  stretchTensor, bubbleField, dataAtPts, true));
2385  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
2386  stretchTensor, rotAxis, dataAtPts,
2387  symmetrySelector == SYMMETRIC ? true : false));
2388  }
2389 
2390  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
2391  spatialL2Disp, piolaStress, dataAtPts, true));
2392  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
2393  spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
2394 
2395  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
2396  piolaStress, rotAxis, dataAtPts,
2397  symmetrySelector == SYMMETRIC ? true : false));
2398  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
2399  bubbleField, rotAxis, dataAtPts,
2400  symmetrySelector == SYMMETRIC ? true : false));
2401 
2402  if (symmetrySelector == SYMMETRIC ? false : true) {
2403  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_du(
2404  rotAxis, stretchTensor, dataAtPts, false));
2405  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
2406  rotAxis, piolaStress, dataAtPts, false));
2407  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
2408  rotAxis, bubbleField, dataAtPts, false));
2409  }
2410  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_domega(
2411  rotAxis, rotAxis, dataAtPts, alphaOmega));
2412 
2413  auto set_hybridisation = [&](auto &pip) {
2414  MoFEMFunctionBegin;
2415 
2416  using BoundaryEle =
2417  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
2418  using EleOnSide =
2419  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
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  CHKERR EshelbianPlasticity::
2434  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2435  op_loop_skeleton_side->getOpPtrVector(), {L2},
2436  materialH1Positions, frontAdjEdges);
2437 
2438  // Second: Iterate over domain FEs adjacent to skelton, particularly one
2439  // domain element.
2440  auto broken_data_ptr =
2441  boost::make_shared<std::vector<BrokenBaseSideData>>();
2442  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2443  auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
2444  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2445  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2446  boost::make_shared<CGGUserPolynomialBase>();
2447  CHKERR
2448  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2449  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2450  materialH1Positions, frontAdjEdges);
2451  op_loop_domain_side->getOpPtrVector().push_back(
2452  new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
2453 
2454  op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
2455  using OpC = FormsIntegrators<BdyEleOp>::Assembly<A>::BiLinearForm<
2456  GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
2457  op_loop_skeleton_side->getOpPtrVector().push_back(
2458  new OpC(hybridSpatialDisp, broken_data_ptr,
2459  boost::make_shared<double>(1.0), true, false));
2460 
2461  pip.push_back(op_loop_skeleton_side);
2462 
2463  MoFEMFunctionReturn(0);
2464  };
2465 
2466  auto set_contact = [&](auto &pip) {
2467  MoFEMFunctionBegin;
2468 
2469  using BoundaryEle =
2470  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
2471  using EleOnSide =
2472  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
2473  using SideEleOp = EleOnSide::UserDataOperator;
2474  using BdyEleOp = BoundaryEle::UserDataOperator;
2475 
2476  // First: Iterate over skeleton FEs adjacent to Domain FEs
2477  // Note: BoundaryEle, i.e. uses skeleton interation rule
2478  auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
2479  mField, contactElement, SPACE_DIM - 1, Sev::noisy);
2480 
2481  op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
2482  op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
2483  CHKERR EshelbianPlasticity::
2484  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2485  op_loop_skeleton_side->getOpPtrVector(), {L2},
2486  materialH1Positions, frontAdjEdges);
2487 
2488  // Second: Iterate over domain FEs adjacent to skelton, particularly
2489  // one domain element.
2490  auto broken_data_ptr =
2491  boost::make_shared<std::vector<BrokenBaseSideData>>();
2492 
2493  // Data storing contact fields
2494  auto contact_common_data_ptr =
2495  boost::make_shared<ContactOps::CommonData>();
2496 
2497  auto add_ops_domain_side = [&](auto &pip) {
2498  MoFEMFunctionBegin;
2499  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2500  auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
2501  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2502  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2503  boost::make_shared<CGGUserPolynomialBase>();
2504  CHKERR
2505  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2506  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2507  materialH1Positions, frontAdjEdges);
2508  op_loop_domain_side->getOpPtrVector().push_back(
2509  new OpGetBrokenBaseSideData<SideEleOp>(piolaStress,
2510  broken_data_ptr));
2511  op_loop_domain_side->getOpPtrVector().push_back(
2512  new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
2513  piolaStress, contact_common_data_ptr->contactTractionPtr()));
2514  pip.push_back(op_loop_domain_side);
2515  MoFEMFunctionReturn(0);
2516  };
2517 
2518  auto add_ops_contact_lhs = [&](auto &pip) {
2519  MoFEMFunctionBegin;
2520  pip.push_back(new OpCalculateVectorFieldValues<3>(
2521  contactDisp, contact_common_data_ptr->contactDispPtr()));
2522  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2523  pip.push_back(
2524  new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
2525  pip.push_back(new OpTreeSearch(
2526  contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
2527  get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
2528  nullptr));
2529 
2530  // get body id and SDF range
2531  auto contact_sfd_map_range_ptr =
2532  boost::make_shared<std::map<int, Range>>(
2533  get_body_range("CONTACT_SDF", SPACE_DIM - 1));
2534 
2535  pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU(
2536  contactDisp, contactDisp, contact_common_data_ptr, contactTreeRhs,
2537  contact_sfd_map_range_ptr));
2538  pip.push_back(
2539  new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP<A, GAUSS>(
2540  contactDisp, broken_data_ptr, contact_common_data_ptr,
2541  contactTreeRhs, contact_sfd_map_range_ptr));
2542  pip.push_back(
2543  new EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU<A, GAUSS>(
2544  broken_data_ptr, contactDisp, contact_common_data_ptr,
2545  contactTreeRhs));
2546 
2547  MoFEMFunctionReturn(0);
2548  };
2549 
2550  // push ops to face/side pipeline
2551  CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
2552  CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
2553 
2554  // Add skeleton to domain pipeline
2555  pip.push_back(op_loop_skeleton_side);
2556 
2557  MoFEMFunctionReturn(0);
2558  };
2559 
2560  CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
2561  CHKERR set_contact(fe_lhs->getOpPtrVector());
2562  }
2563 
2564  if (add_material) {
2565  }
2566 
2567  MoFEMFunctionReturn(0);
2568 }

◆ solveDynamicRelaxation()

MoFEMErrorCode EshelbianCore::solveDynamicRelaxation ( TS  ts,
Vec  x 
)
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 2995 of file EshelbianPlasticity.cpp.

2995  {
2996  MoFEMFunctionBegin;
2997 
2998  auto storage = solve_elastic_setup::setup(this, ts, x, false);
2999 
3000  double final_time = 1;
3001  double delta_time = 0.1;
3002  int max_it = 10;
3003  PetscBool ts_h1_update = PETSC_FALSE;
3004 
3005  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Dynamic Relaxation Options",
3006  "none");
3007 
3008  CHKERR PetscOptionsScalar("-dynamic_final_time",
3009  "dynamic relaxation final time", "", final_time,
3010  &final_time, PETSC_NULL);
3011  CHKERR PetscOptionsScalar("-dynamic_delta_time",
3012  "dynamic relaxation final time", "", delta_time,
3013  &delta_time, PETSC_NULL);
3014  CHKERR PetscOptionsInt("-dynamic_max_it", "dynamic relaxation iterations", "",
3015  max_it, &max_it, PETSC_NULL);
3016  CHKERR PetscOptionsBool("-dynamic_h1_update", "update each ts step", "",
3017  ts_h1_update, &ts_h1_update, PETSC_NULL);
3018 
3019  ierr = PetscOptionsEnd();
3020  CHKERRG(ierr);
3021 
3022  auto setup_ts_monitor = [&]() {
3023  auto monitor_ptr = boost::make_shared<EshelbianMonitor>(*this);
3024  return monitor_ptr;
3025  };
3026  auto monitor_ptr = setup_ts_monitor();
3027 
3028  TetPolynomialBase::switchCacheBaseOn<HDIV>(
3029  {elasticFeLhs.get(), elasticFeRhs.get()});
3030  CHKERR TSSetUp(ts);
3031  CHKERR TSElasticPostStep::postStepInitialise(this);
3032 
3033  double ts_delta_time;
3034  CHKERR TSGetTimeStep(ts, &ts_delta_time);
3035 
3036  if (ts_h1_update) {
3037  CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
3038  CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
3039  }
3040 
3041  CHKERR TSElasticPostStep::preStepFun(ts);
3042  CHKERR TSElasticPostStep::postStepFun(ts);
3043 
3044  dynamicTime = 0;
3045  dynamicStep = 0;
3046  monitor_ptr->ts_u = PETSC_NULL;
3047  monitor_ptr->ts_t = dynamicTime;
3048  monitor_ptr->ts_step = dynamicStep;
3049  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
3050  MOFEM_LOG("EP", Sev::inform) << "IT " << dynamicStep << " Time "
3051  << dynamicTime << " delta time " << delta_time;
3052  dynamicTime += delta_time;
3053  ++dynamicStep;
3054 
3055  for (; dynamicTime < final_time; dynamicTime += delta_time) {
3056  MOFEM_LOG("EP", Sev::inform) << "IT " << dynamicStep << " Time "
3057  << dynamicTime << " delta time " << delta_time;
3058 
3059  CHKERR TSSetStepNumber(ts, 0);
3060  CHKERR TSSetTime(ts, 0);
3061  CHKERR TSSetTimeStep(ts, ts_delta_time);
3062  if (!ts_h1_update) {
3063  CHKERR TSElasticPostStep::preStepFun(ts);
3064  }
3065  CHKERR TSSolve(ts, PETSC_NULL);
3066  if (!ts_h1_update) {
3067  CHKERR TSElasticPostStep::postStepFun(ts);
3068  }
3069 
3070  monitor_ptr->ts_u = PETSC_NULL;
3071  monitor_ptr->ts_t = dynamicTime;
3072  monitor_ptr->ts_step = dynamicStep;
3073  CHKERR DMoFEMLoopFiniteElements(dmElastic, elementVolumeName, monitor_ptr);
3074 
3075  ++dynamicStep;
3076  if (dynamicStep > max_it)
3077  break;
3078  }
3079 
3080  CHKERR TSElasticPostStep::postStepDestroy();
3081  TetPolynomialBase::switchCacheBaseOff<HDIV>(
3082  {elasticFeLhs.get(), elasticFeRhs.get()});
3083 
3084  MoFEMFunctionReturn(0);
3085 }

◆ solveElastic()

MoFEMErrorCode EshelbianCore::solveElastic ( TS  ts,
Vec  x 
)
Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 2920 of file EshelbianPlasticity.cpp.

2920  {
2921  MoFEMFunctionBegin;
2922 
2923  PetscBool debug_model = PETSC_FALSE;
2924  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-debug_model", &debug_model,
2925  PETSC_NULL);
2926 
2927  if (debug_model == PETSC_TRUE) {
2928  auto ts_ctx_ptr = getDMTsCtx(dmElastic);
2929  auto post_proc = [&](TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, Vec F,
2930  void *ctx) {
2931  MoFEMFunctionBegin;
2932 
2933  SNES snes;
2934  CHKERR TSGetSNES(ts, &snes);
2935  int it;
2936  CHKERR SNESGetIterationNumber(snes, &it);
2937  std::string file_name = "snes_iteration_" + std::to_string(it) + ".h5m";
2938  CHKERR postProcessResults(1, file_name, F);
2939  std::string file_skel_name =
2940  "snes_iteration_skel_" + std::to_string(it) + ".h5m";
2941  CHKERR postProcessSkeletonResults(1, file_skel_name, F);
2942 
2943  MoFEMFunctionReturn(0);
2944  };
2945  ts_ctx_ptr->tsDebugHook = post_proc;
2946  }
2947 
2948  auto storage = solve_elastic_setup::setup(this, ts, x, true);
2949 
2950  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
2951  Vec xx;
2952  CHKERR VecDuplicate(x, &xx);
2953  CHKERR VecZeroEntries(xx);
2954  CHKERR TS2SetSolution(ts, x, xx);
2955  CHKERR VecDestroy(&xx);
2956  } else {
2957  CHKERR TSSetSolution(ts, x);
2958  }
2959 
2960  TetPolynomialBase::switchCacheBaseOn<HDIV>(
2961  {elasticFeLhs.get(), elasticFeRhs.get()});
2962  CHKERR TSSetUp(ts);
2963  CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
2964  CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
2965  CHKERR TSElasticPostStep::postStepInitialise(this);
2966  CHKERR TSSolve(ts, PETSC_NULL);
2967  CHKERR TSElasticPostStep::postStepDestroy();
2968  TetPolynomialBase::switchCacheBaseOff<HDIV>(
2969  {elasticFeLhs.get(), elasticFeRhs.get()});
2970 
2971  SNES snes;
2972  CHKERR TSGetSNES(ts, &snes);
2973  int lin_solver_iterations;
2974  CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
2975  MOFEM_LOG("EP", Sev::inform)
2976  << "Number of linear solver iterations " << lin_solver_iterations;
2977 
2978  PetscBool test_cook_flg = PETSC_FALSE;
2979  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_cook", &test_cook_flg,
2980  PETSC_NULL);
2981  if (test_cook_flg) {
2982  constexpr int expected_lin_solver_iterations = 11;
2983  if (lin_solver_iterations > expected_lin_solver_iterations)
2984  SETERRQ2(
2985  PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
2986  "Expected number of iterations is different than expected %d > %d",
2987  lin_solver_iterations, expected_lin_solver_iterations);
2988  }
2989 
2990  CHKERR gettingNorms();
2991 
2992  MoFEMFunctionReturn(0);
2993 }

Friends And Related Function Documentation

◆ solve_elastic_set_up

friend struct solve_elastic_set_up
friend

Definition at line 248 of file EshelbianCore.hpp.

Member Data Documentation

◆ a00FieldList

std::vector<std::string> EshelbianCore::a00FieldList
Examples
EshelbianPlasticity.cpp.

Definition at line 315 of file EshelbianCore.hpp.

◆ a00RangeList

std::vector<boost::shared_ptr<Range> > EshelbianCore::a00RangeList
Examples
EshelbianPlasticity.cpp.

Definition at line 317 of file EshelbianCore.hpp.

◆ addCrackMeshsetId

int EshelbianCore::addCrackMeshsetId = 1000
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 28 of file EshelbianCore.hpp.

◆ alphaOmega

double EshelbianCore::alphaOmega = 0
Examples
EshelbianPlasticity.cpp.

Definition at line 129 of file EshelbianCore.hpp.

◆ alphaRho

double EshelbianCore::alphaRho = 0
Examples
EshelbianPlasticity.cpp.

Definition at line 130 of file EshelbianCore.hpp.

◆ alphaU

double EshelbianCore::alphaU = 0
Examples
EshelbianPlasticity.cpp.

Definition at line 127 of file EshelbianCore.hpp.

◆ alphaW

double EshelbianCore::alphaW = 0
Examples
EshelbianPlasticity.cpp.

Definition at line 128 of file EshelbianCore.hpp.

◆ aoS

AO EshelbianCore::aoS = PETSC_NULL
Examples
EshelbianPlasticity.cpp.

Definition at line 313 of file EshelbianCore.hpp.

◆ bcSpatialDispVecPtr

boost::shared_ptr<BcDispVec> EshelbianCore::bcSpatialDispVecPtr
Examples
EshelbianPlasticity.cpp.

Definition at line 137 of file EshelbianCore.hpp.

◆ bcSpatialFreeTraction

boost::shared_ptr<TractionFreeBc> EshelbianCore::bcSpatialFreeTraction
Examples
EshelbianPlasticity.cpp.

Definition at line 140 of file EshelbianCore.hpp.

◆ bcSpatialRotationVecPtr

boost::shared_ptr<BcRotVec> EshelbianCore::bcSpatialRotationVecPtr
Examples
EshelbianPlasticity.cpp.

Definition at line 138 of file EshelbianCore.hpp.

◆ bcSpatialTraction

boost::shared_ptr<TractionBcVec> EshelbianCore::bcSpatialTraction
Examples
EshelbianPlasticity.cpp.

Definition at line 139 of file EshelbianCore.hpp.

◆ bitAdjEnt

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

bit ref level for parent

Examples
EshelbianPlasticity.cpp.

Definition at line 298 of file EshelbianCore.hpp.

◆ bitAdjEntMask

BitRefLevel EshelbianCore::bitAdjEntMask
Initial value:
=
BitRefLevel().set()

bit ref level for parent parent

Examples
EshelbianPlasticity.cpp.

Definition at line 299 of file EshelbianCore.hpp.

◆ bitAdjParent

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

bit ref level for parent

Examples
EshelbianPlasticity.cpp.

Definition at line 295 of file EshelbianCore.hpp.

◆ bitAdjParentMask

BitRefLevel EshelbianCore::bitAdjParentMask
Initial value:
=
BitRefLevel().set()

bit ref level for parent parent

Examples
EshelbianPlasticity.cpp.

Definition at line 296 of file EshelbianCore.hpp.

◆ bubbleField

const std::string EshelbianCore::bubbleField = "bubble"
Examples
EshelbianPlasticity.cpp.

Definition at line 111 of file EshelbianCore.hpp.

◆ contactDisp

const std::string EshelbianCore::contactDisp = "contactDisp"
Examples
EshelbianPlasticity.cpp.

Definition at line 108 of file EshelbianCore.hpp.

◆ contactElement

const std::string EshelbianCore::contactElement = "CONTACT"
Examples
EshelbianPlasticity.cpp.

Definition at line 117 of file EshelbianCore.hpp.

◆ contactFaces

boost::shared_ptr<Range> EshelbianCore::contactFaces
Examples
EshelbianPlasticity.cpp.

Definition at line 283 of file EshelbianCore.hpp.

◆ contactRefinementLevels

int EshelbianCore::contactRefinementLevels = 1
Examples
EshelbianPlasticity.cpp.

Definition at line 132 of file EshelbianCore.hpp.

◆ contactTreeRhs

boost::shared_ptr<ContactTree> EshelbianCore::contactTreeRhs

Make a contact tree.

Examples
EshelbianPlasticity.cpp.

Definition at line 93 of file EshelbianCore.hpp.

◆ crackFaces

boost::shared_ptr<Range> EshelbianCore::crackFaces
Examples
EshelbianPlasticity.cpp.

Definition at line 284 of file EshelbianCore.hpp.

◆ crackHybridIs

SmartPetscObj<IS> EshelbianCore::crackHybridIs

Definition at line 314 of file EshelbianCore.hpp.

◆ crackingOn

PetscBool EshelbianCore::crackingOn = PETSC_FALSE
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 22 of file EshelbianCore.hpp.

◆ d_f

boost::function< double(const double)> EshelbianCore::d_f
static
Initial value:
=
EshelbianCore::d_f_log_e
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 33 of file EshelbianCore.hpp.

◆ dataAtPts

boost::shared_ptr<DataAtIntegrationPts> EshelbianCore::dataAtPts
Examples
EshelbianPlasticity.cpp.

Definition at line 86 of file EshelbianCore.hpp.

◆ dd_f

boost::function< double(const double)> EshelbianCore::dd_f
static
Initial value:
=
EshelbianCore::dd_f_log_e
Examples
EshelbianPlasticity.cpp.

Definition at line 34 of file EshelbianCore.hpp.

◆ dM

SmartPetscObj<DM> EshelbianCore::dM

Coupled problem all fields.

Examples
EshelbianPlasticity.cpp.

Definition at line 95 of file EshelbianCore.hpp.

◆ dmElastic

SmartPetscObj<DM> EshelbianCore::dmElastic

Elastic problem.

Examples
ep.cpp, and EshelbianPlasticity.cpp.

Definition at line 96 of file EshelbianCore.hpp.

◆ dmMaterial

SmartPetscObj<DM> EshelbianCore::dmMaterial

Material problem.

Definition at line 97 of file EshelbianCore.hpp.

◆ dmPrjSpatial

SmartPetscObj<DM> EshelbianCore::dmPrjSpatial

Projection spatial displacement.

Examples
EshelbianPlasticity.cpp.

Definition at line 98 of file EshelbianCore.hpp.

◆ dynamicRelaxation

PetscBool EshelbianCore::dynamicRelaxation
inlinestatic
Initial value:
=
PETSC_FALSE
Examples
ep.cpp, EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 20 of file EshelbianCore.hpp.

◆ dynamicStep

int EshelbianCore::dynamicStep
inlinestatic
Initial value:
=
0
Examples
EshelbianPlasticity.cpp.

Definition at line 25 of file EshelbianCore.hpp.

◆ dynamicTime

double EshelbianCore::dynamicTime
inlinestatic
Initial value:
=
0
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 23 of file EshelbianCore.hpp.

◆ edgeExchange

CommInterface::EntitiesPetscVector EshelbianCore::edgeExchange
Examples
EshelbianPlasticity.cpp.

Definition at line 306 of file EshelbianCore.hpp.

◆ elasticBcLhs

boost::shared_ptr<FaceElementForcesAndSourcesCore> EshelbianCore::elasticBcLhs
Examples
EshelbianPlasticity.cpp.

Definition at line 91 of file EshelbianCore.hpp.

◆ elasticBcRhs

boost::shared_ptr<FaceElementForcesAndSourcesCore> EshelbianCore::elasticBcRhs
Examples
EshelbianPlasticity.cpp.

Definition at line 92 of file EshelbianCore.hpp.

◆ elasticFeLhs

boost::shared_ptr<VolumeElementForcesAndSourcesCore> EshelbianCore::elasticFeLhs
Examples
EshelbianPlasticity.cpp.

Definition at line 90 of file EshelbianCore.hpp.

◆ elasticFeRhs

boost::shared_ptr<VolumeElementForcesAndSourcesCore> EshelbianCore::elasticFeRhs
Examples
EshelbianPlasticity.cpp.

Definition at line 89 of file EshelbianCore.hpp.

◆ elementVolumeName

const std::string EshelbianCore::elementVolumeName = "EP"
Examples
EshelbianPlasticity.cpp.

Definition at line 113 of file EshelbianCore.hpp.

◆ eshelbyStress

const std::string EshelbianCore::eshelbyStress = "S"
Examples
EshelbianPlasticity.cpp.

Definition at line 101 of file EshelbianCore.hpp.

◆ exponentBase

double EshelbianCore::exponentBase = exp(1)
static
Examples
EshelbianPlasticity.cpp.

Definition at line 31 of file EshelbianCore.hpp.

◆ f

boost::function< double(const double)> EshelbianCore::f = EshelbianCore::f_log_e
static
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 32 of file EshelbianCore.hpp.

◆ faceExchange

CommInterface::EntitiesPetscVector EshelbianCore::faceExchange
Examples
EshelbianPlasticity.cpp.

Definition at line 305 of file EshelbianCore.hpp.

◆ frontAdjEdges

boost::shared_ptr<Range> EshelbianCore::frontAdjEdges
Examples
EshelbianPlasticity.cpp.

Definition at line 286 of file EshelbianCore.hpp.

◆ frontEdges

boost::shared_ptr<Range> EshelbianCore::frontEdges
Examples
EshelbianPlasticity.cpp.

Definition at line 285 of file EshelbianCore.hpp.

◆ frontLayers

int EshelbianCore::frontLayers = 3
Examples
EshelbianPlasticity.cpp.

Definition at line 133 of file EshelbianCore.hpp.

◆ frontVertices

boost::shared_ptr<Range> EshelbianCore::frontVertices
Examples
EshelbianPlasticity.cpp.

Definition at line 287 of file EshelbianCore.hpp.

◆ gradApproximator

enum RotSelector EshelbianCore::gradApproximator = LARGE_ROT
inlinestatic
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 16 of file EshelbianCore.hpp.

◆ griffithEnergy

double EshelbianCore::griffithEnergy = 1
inlinestatic

Griffith energy.

Examples
EshelbianPlasticity.cpp.

Definition at line 29 of file EshelbianCore.hpp.

◆ hybridMaterialDisp

const std::string EshelbianCore::hybridMaterialDisp = "hybridMaterialDisp"
Examples
EshelbianPlasticity.cpp.

Definition at line 107 of file EshelbianCore.hpp.

◆ hybridSpatialDisp

const std::string EshelbianCore::hybridSpatialDisp = "hybridSpatialDisp"
Examples
EshelbianPlasticity.cpp.

Definition at line 106 of file EshelbianCore.hpp.

◆ inv_d_f

boost::function< double(const double)> EshelbianCore::inv_d_f
static
Initial value:
=
EshelbianCore::inv_d_f_log_e
Examples
EshelbianPlasticity.cpp.

Definition at line 36 of file EshelbianCore.hpp.

◆ inv_dd_f

boost::function< double(const double)> EshelbianCore::inv_dd_f
static
Initial value:
=
EshelbianCore::inv_dd_f_log_e
Examples
EshelbianPlasticity.cpp.

Definition at line 37 of file EshelbianCore.hpp.

◆ inv_f

boost::function< double(const double)> EshelbianCore::inv_f
static
Initial value:
=
EshelbianCore::inv_f_log_e
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 35 of file EshelbianCore.hpp.

◆ l2UserBaseScale

PetscBool EshelbianCore::l2UserBaseScale = PETSC_FALSE
inlinestatic
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 27 of file EshelbianCore.hpp.

◆ listTagsToTransfer

std::vector<Tag> EshelbianCore::listTagsToTransfer

list of tags to transfer to postprocessor

Examples
EshelbianPlasticity.cpp.

Definition at line 310 of file EshelbianCore.hpp.

◆ materialH1Positions

const std::string EshelbianCore::materialH1Positions = "XH1"
Examples
EshelbianPlasticity.cpp.

Definition at line 105 of file EshelbianCore.hpp.

◆ materialL2Disp

const std::string EshelbianCore::materialL2Disp = "WL2"
Examples
EshelbianPlasticity.cpp.

Definition at line 103 of file EshelbianCore.hpp.

◆ materialOrder

int EshelbianCore::materialOrder = 1
Examples
EshelbianPlasticity.cpp.

Definition at line 126 of file EshelbianCore.hpp.

◆ materialSkeletonElement

const std::string EshelbianCore::materialSkeletonElement = "MATERIAL_SKELETON"
Examples
EshelbianPlasticity.cpp.

Definition at line 119 of file EshelbianCore.hpp.

◆ materialSkeletonFaces

boost::shared_ptr<Range> EshelbianCore::materialSkeletonFaces
Examples
EshelbianPlasticity.cpp.

Definition at line 289 of file EshelbianCore.hpp.

◆ materialVolumeElement

const std::string EshelbianCore::materialVolumeElement = "MEP"
Examples
EshelbianPlasticity.cpp.

Definition at line 118 of file EshelbianCore.hpp.

◆ maxMovedFaces

boost::shared_ptr<Range> EshelbianCore::maxMovedFaces
Examples
EshelbianPlasticity.cpp.

Definition at line 290 of file EshelbianCore.hpp.

◆ mField

MoFEM::Interface& EshelbianCore::mField
Examples
EshelbianPlasticity.cpp.

Definition at line 84 of file EshelbianCore.hpp.

◆ naturalBcElement

const std::string EshelbianCore::naturalBcElement = "NATURAL_BC"
Examples
EshelbianPlasticity.cpp.

Definition at line 114 of file EshelbianCore.hpp.

◆ noStretch

PetscBool EshelbianCore::noStretch = PETSC_FALSE
inlinestatic
Examples
EshelbianPlasticity.cpp.

Definition at line 18 of file EshelbianCore.hpp.

◆ parentAdjSkeletonFunctionDim2

boost::shared_ptr<ParentFiniteElementAdjacencyFunctionSkeleton<2> > EshelbianCore::parentAdjSkeletonFunctionDim2
Examples
EshelbianPlasticity.cpp.

Definition at line 293 of file EshelbianCore.hpp.

◆ physicalEquations

boost::shared_ptr<PhysicalEquations> EshelbianCore::physicalEquations
Examples
EshelbianPlasticity.cpp.

Definition at line 87 of file EshelbianCore.hpp.

◆ piolaStress

const std::string EshelbianCore::piolaStress = "P"
Examples
EshelbianPlasticity.cpp.

Definition at line 100 of file EshelbianCore.hpp.

◆ rotAxis

const std::string EshelbianCore::rotAxis = "omega"
Examples
EshelbianPlasticity.cpp.

Definition at line 110 of file EshelbianCore.hpp.

◆ rotSelector

enum RotSelector EshelbianCore::rotSelector = LARGE_ROT
inlinestatic
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 15 of file EshelbianCore.hpp.

◆ S

Mat EshelbianCore::S = PETSC_NULL
Examples
EshelbianPlasticity.cpp.

Definition at line 312 of file EshelbianCore.hpp.

◆ setSingularity

PetscBool EshelbianCore::setSingularity = PETSC_TRUE
inlinestatic
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 19 of file EshelbianCore.hpp.

◆ skeletonElement

const std::string EshelbianCore::skeletonElement = "SKELETON"
Examples
EshelbianPlasticity.cpp.

Definition at line 116 of file EshelbianCore.hpp.

◆ skeletonFaces

boost::shared_ptr<Range> EshelbianCore::skeletonFaces
Examples
EshelbianPlasticity.cpp.

Definition at line 288 of file EshelbianCore.hpp.

◆ skinElement

const std::string EshelbianCore::skinElement = "SKIN"
Examples
EshelbianPlasticity.cpp.

Definition at line 115 of file EshelbianCore.hpp.

◆ solTSStep

SmartPetscObj<Vec> EshelbianCore::solTSStep
Examples
EshelbianPlasticity.cpp.

Definition at line 302 of file EshelbianCore.hpp.

◆ spaceH1Order

int EshelbianCore::spaceH1Order = -1
Examples
EshelbianPlasticity.cpp.

Definition at line 125 of file EshelbianCore.hpp.

◆ spaceOrder

int EshelbianCore::spaceOrder = 2
Examples
EshelbianPlasticity.cpp.

Definition at line 124 of file EshelbianCore.hpp.

◆ spatialH1Disp

const std::string EshelbianCore::spatialH1Disp = "wH1"
Examples
EshelbianPlasticity.cpp.

Definition at line 104 of file EshelbianCore.hpp.

◆ spatialL2Disp

const std::string EshelbianCore::spatialL2Disp = "wL2"
Examples
EshelbianPlasticity.cpp.

Definition at line 102 of file EshelbianCore.hpp.

◆ stretchSelector

enum StretchSelector EshelbianCore::stretchSelector = LOG
inlinestatic
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 17 of file EshelbianCore.hpp.

◆ stretchTensor

const std::string EshelbianCore::stretchTensor = "u"
Examples
EshelbianPlasticity.cpp.

Definition at line 109 of file EshelbianCore.hpp.

◆ symmetrySelector

enum SymmetrySelector EshelbianCore::symmetrySelector = SYMMETRIC
inlinestatic
Examples
EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 14 of file EshelbianCore.hpp.

◆ vertexExchange

CommInterface::EntitiesPetscVector EshelbianCore::vertexExchange
Examples
EshelbianPlasticity.cpp.

Definition at line 307 of file EshelbianCore.hpp.

◆ volumeExchange

CommInterface::EntitiesPetscVector EshelbianCore::volumeExchange
Examples
EshelbianPlasticity.cpp.

Definition at line 304 of file EshelbianCore.hpp.

The documentation for this struct was generated from the following files:
MoFEM::VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator
default operator for TET element
Definition: VolumeElementForcesAndSourcesCoreOnSide.hpp:107
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::CoreInterface::modify_finite_element_adjacency_table
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
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
a2
constexpr double a2
Definition: hcurl_check_approx_in_2d.cpp:39
OpSpatialConsistencyDivTerm
Definition: EshelbianOperators.hpp:265
EshelbianCore::materialL2Disp
const std::string materialL2Disp
Definition: EshelbianCore.hpp:103
MoFEM::id_from_handle
auto id_from_handle(const EntityHandle h)
Definition: Templates.hpp:1890
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:712
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
EshelbianCore::inv_d_f_log
static double inv_d_f_log(const double v)
Definition: EshelbianCore.hpp:60
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MoFEM::CoreInterface::loop_finite_elements
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.
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::OpPostProcMapInMoab::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: PostProcBrokenMeshInMoabBase.hpp:750
EshelbianPlasticity::LINEAR
@ LINEAR
Definition: EshelbianPlasticity.hpp:46
MoFEM::CoreInterface::loop_dofs
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.
EshelbianCore::inv_d_f_linear
static double inv_d_f_linear(const double v)
Definition: EshelbianCore.hpp:72
H1
@ H1
continuous field
Definition: definitions.h:85
EshelbianCore::dd_f
static boost::function< double(const double)> dd_f
Definition: EshelbianCore.hpp:34
EshelbianCore::setVolumeElementOps
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)
Definition: EshelbianPlasticity.cpp:2119
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
save_range
static auto save_range(moab::Interface &moab, const std::string name, const Range r, std::vector< Tag > tags={})
Definition: EshelbianPlasticity.cpp:128
EshelbianCore::bitAdjParent
BitRefLevel bitAdjParent
bit ref level for parent
Definition: EshelbianCore.hpp:295
EshelbianCore::a00FieldList
std::vector< std::string > a00FieldList
Definition: EshelbianCore.hpp:315
EshelbianCore::S
Mat S
Definition: EshelbianCore.hpp:312
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
EshelbianCore::spatialL2Disp
const std::string spatialL2Disp
Definition: EshelbianCore.hpp:102
EshelbianCore::contactFaces
boost::shared_ptr< Range > contactFaces
Definition: EshelbianCore.hpp:283
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
OpSpatialConsistencyBubble
Definition: EshelbianOperators.hpp:256
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:373
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
EshelbianCore::f
static boost::function< double(const double)> f
Definition: EshelbianCore.hpp:32
EleOnSide
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
Definition: scalar_check_approximation.cpp:27
EshelbianCore::elasticFeLhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
Definition: EshelbianCore.hpp:90
EshelbianCore::frontLayers
int frontLayers
Definition: EshelbianCore.hpp:133
MoFEM::BaseFunction::DofsSideMap
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.
Definition: BaseFunction.hpp:73
EshelbianCore::materialSkeletonFaces
boost::shared_ptr< Range > materialSkeletonFaces
Definition: EshelbianCore.hpp:289
EshelbianCore::bcSpatialFreeTraction
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
Definition: EshelbianCore.hpp:140
MoFEM::Tools
Auxiliary tools.
Definition: Tools.hpp:19
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
EshelbianCore::noStretch
static PetscBool noStretch
Definition: EshelbianCore.hpp:18
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot
Calculate symmetric tensor field rates ant integratio pts.
Definition: UserDataOperators.hpp:1097
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
TSElasticPostStep::postStepFun
static MoFEMErrorCode postStepFun(TS ts)
Definition: TSElasticPostStep.cpp:156
MoFEM::OpCalculateHTensorTensorField
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2855
EshelbianCore::naturalBcElement
const std::string naturalBcElement
Definition: EshelbianCore.hpp:114
MoFEM::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
EshelbianCore::stretchSelector
static enum StretchSelector stretchSelector
Definition: EshelbianCore.hpp:17
MoFEM::NaturalBC
Natural boundary conditions.
Definition: Natural.hpp:57
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
E
FTensor::Tensor1::l2
T l2() const
Definition: Tensor1_value.hpp:84
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468
_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_
#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
Definition: Interface.hpp:1917
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
EshelbianCore::alphaOmega
double alphaOmega
Definition: EshelbianCore.hpp:129
EshelbianCore::contactTreeRhs
boost::shared_ptr< ContactTree > contactTreeRhs
Make a contact tree.
Definition: EshelbianCore.hpp:93
EshelbianCore::materialVolumeElement
const std::string materialVolumeElement
Definition: EshelbianCore.hpp:118
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
EshelbianCore::postProcessSkeletonResults
MoFEMErrorCode postProcessSkeletonResults(const int tag, const std::string file, Vec f_residual=PETSC_NULL)
Definition: EshelbianPlasticity.cpp:3558
EshelbianCore::getTractionFreeBc
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.
Definition: EshelbianPlasticity.cpp:1796
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
EshelbianPlasticity::SymmetrySelector
SymmetrySelector
Definition: EshelbianPlasticity.hpp:44
OpPostProcDataStructure
Definition: EshelbianOperators.hpp:527
filter_true_skin
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
Definition: EshelbianPlasticity.cpp:142
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::ForceCubitBcData
Definition of the force bc data structure.
Definition: BCData.hpp:139
EshelbianCore::frontVertices
boost::shared_ptr< Range > frontVertices
Definition: EshelbianCore.hpp:287
OpSpatialConsistency_dBubble_dP
Definition: EshelbianOperators.hpp:487
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::CoreInterface::add_broken_field
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.
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::solveLinearSystem
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
Definition: Templates.hpp:1401
MoFEM::OpSetFlux
Definition: FormsBrokenSpaceConstraintImpl.hpp:139
NBVOLUMETET_L2
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:27
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition: Templates.hpp:2013
EshelbianCore::f_linear
static double f_linear(const double v)
Definition: EshelbianCore.hpp:67
EshelbianPlasticity::NO_H1_CONFIGURATION
@ NO_H1_CONFIGURATION
Definition: EshelbianPlasticity.hpp:45
FTensor::levi_civita
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
Definition: Levi_Civita.hpp:617
EshelbianCore::dd_f_log_e
static double dd_f_log_e(const double v)
Definition: EshelbianCore.hpp:41
EshelbianCore::physicalEquations
boost::shared_ptr< PhysicalEquations > physicalEquations
Definition: EshelbianCore.hpp:87
EshelbianCore::griffithEnergy
static double griffithEnergy
Griffith energy.
Definition: EshelbianCore.hpp:29
LieGroups::SO3::exp
static auto exp(A &&t_w_vee, B &&theta)
Definition: Lie.hpp:49
EshelbianCore::spaceOrder
int spaceOrder
Definition: EshelbianCore.hpp:124
sdf.r
int r
Definition: sdf.py:8
MoFEM::CoreInterface::add_ents_to_field_by_dim
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.
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:51
USER_BASE
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
EshelbianCore::frontEdges
boost::shared_ptr< Range > frontEdges
Definition: EshelbianCore.hpp:285
get_crack_front_edges
static auto get_crack_front_edges(MoFEM::Interface &m_field, Range crack_faces)
Definition: EshelbianPlasticity.cpp:170
OpRotationBc
Definition: EshelbianOperators.hpp:326
EshelbianCore::setSingularity
static PetscBool setSingularity
Definition: EshelbianCore.hpp:19
MoFEM::getBlockMatStorageMat
boost::shared_ptr< std::vector< double > > getBlockMatStorageMat(Mat B)
Get the Block Storage object.
Definition: Schur.cpp:2132
EshelbianCore::contactDisp
const std::string contactDisp
Definition: EshelbianCore.hpp:108
EshelbianCore::setBaseVolumeElementOps
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)
Definition: EshelbianPlasticity.cpp:1999
ROW
@ ROW
Definition: definitions.h:136
OpSpatialEquilibrium_dw_dP
Definition: EshelbianOperators.hpp:352
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
EshelbianCore::skeletonFaces
boost::shared_ptr< Range > skeletonFaces
Definition: EshelbianCore.hpp:288
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:597
EshelbianCore::maxMovedFaces
boost::shared_ptr< Range > maxMovedFaces
Definition: EshelbianCore.hpp:290
EshelbianCore::gradApproximator
static enum RotSelector gradApproximator
Definition: EshelbianCore.hpp:16
SideEleOp
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:438
EshelbianCore::dmElastic
SmartPetscObj< DM > dmElastic
Elastic problem.
Definition: EshelbianCore.hpp:96
EshelbianCore::parentAdjSkeletonFunctionDim2
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > > parentAdjSkeletonFunctionDim2
Definition: EshelbianCore.hpp:293
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
MoFEM::OpCalcNormL2Tensor2
Get norm of input MatrixDouble for Tensor2.
Definition: NormsOperators.hpp:72
MoFEM::OpGetBrokenBaseSideData
Definition: FormsBrokenSpaceConstraintImpl.hpp:68
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
EshelbianCore::alphaRho
double alphaRho
Definition: EshelbianCore.hpp:130
EshelbianPlasticity::SYMMETRIC
@ SYMMETRIC
Definition: EshelbianPlasticity.hpp:44
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::OpBrokenLoopSide
Definition: FormsBrokenSpaceConstraintImpl.hpp:15
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
send_type
static auto send_type(MoFEM::Interface &m_field, Range r, const EntityType type)
Definition: EshelbianPlasticity.cpp:52
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DMMoFEMTSSetIJacobian
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:853
EshelbianCore::inv_f_linear
static double inv_f_linear(const double v)
Definition: EshelbianCore.hpp:71
MoFEM::DMoFEMLoopFiniteElementsUpAndLowRank
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:567
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
EshelbianCore::bitAdjEnt
BitRefLevel bitAdjEnt
bit ref level for parent
Definition: EshelbianCore.hpp:298
OpSpatialConsistency_dP_dP
Definition: EshelbianOperators.hpp:460
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
TSElasticPostStep::preStepFun
static MoFEMErrorCode preStepFun(TS ts)
Definition: TSElasticPostStep.cpp:95
a1
constexpr double a1
Definition: hcurl_check_approx_in_2d.cpp:38
EshelbianCore::contactElement
const std::string contactElement
Definition: EshelbianCore.hpp:117
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
a
constexpr double a
Definition: approx_sphere.cpp:30
EshelbianCore::getOptions
MoFEMErrorCode getOptions()
Definition: EshelbianPlasticity.cpp:883
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
EshelbianCore::dM
SmartPetscObj< DM > dM
Coupled problem all fields.
Definition: EshelbianCore.hpp:95
EshelbianCore::rotAxis
const std::string rotAxis
Definition: EshelbianCore.hpp:110
BoundaryEleOp
MoFEM::DMMoFEMTSSetI2Jacobian
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:1017
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
EshelbianCore::alphaU
double alphaU
Definition: EshelbianCore.hpp:127
EshelbianPlasticity::TractionFreeBc
std::vector< Range > TractionFreeBc
Definition: EshelbianPlasticity.hpp:435
MoFEM::DMMoFEMTSSetI2Function
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:975
EshelbianCore::eshelbyStress
const std::string eshelbyStress
Definition: EshelbianCore.hpp:101
MoFEM::DMMoFEMCreateSubDM
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
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
EshelbianCore::materialH1Positions
const std::string materialH1Positions
Definition: EshelbianCore.hpp:105
EshelbianCore::gettingNorms
MoFEMErrorCode gettingNorms()
[Getting norms]
Definition: EshelbianPlasticity.cpp:5226
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
EshelbianCore::inv_dd_f_log
static double inv_dd_f_log(const double v)
Definition: EshelbianCore.hpp:63
EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU
Definition: EshelbianContact.hpp:198
EshelbianCore::elementVolumeName
const std::string elementVolumeName
Definition: EshelbianCore.hpp:113
EshelbianCore::d_f_log
static double d_f_log(const double v)
Definition: EshelbianCore.hpp:49
OpSpatialPhysical_du_dBubble
Definition: EshelbianOperators.hpp:388
EshelbianCore::spatialH1Disp
const std::string spatialH1Disp
Definition: EshelbianCore.hpp:104
OpCalculateEshelbyStress
Definition: HookeElement.hpp:177
COL
@ COL
Definition: definitions.h:136
OpReleaseEnergy
Definition: EshelbianOperators.hpp:577
EshelbianPlasticity::OpConstrainBoundaryHDivRhs
Definition: EshelbianContact.hpp:134
EshelbianCore::aoS
AO aoS
Definition: EshelbianCore.hpp:313
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2628
EshelbianPlasticity::StretchSelector
StretchSelector
Definition: EshelbianPlasticity.hpp:46
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
EshelbianCore::materialSkeletonElement
const std::string materialSkeletonElement
Definition: EshelbianCore.hpp:119
EshelbianCore::mField
MoFEM::Interface & mField
Definition: EshelbianCore.hpp:84
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::solve
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
Definition: Schur.cpp:1296
EshelbianCore::bitAdjEntMask
BitRefLevel bitAdjEntMask
bit ref level for parent parent
Definition: EshelbianCore.hpp:299
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU
Definition: EshelbianContact.hpp:155
get_skin
static auto get_skin(MoFEM::Interface &m_field, Range body_ents)
Definition: EshelbianPlasticity.cpp:163
EshelbianCore::addCrackMeshsetId
static int addCrackMeshsetId
Definition: EshelbianCore.hpp:28
OpSpatialPhysical_du_domega
Definition: EshelbianOperators.hpp:400
MoFEM::PostProcBrokenMeshInMoabBaseEnd
PostProcBrokenMeshInMoabBaseEndImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseEnd
Enable to run stack of post-processing elements. Use this to end stack.
Definition: PostProcBrokenMeshInMoabBase.hpp:962
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
MoFEM::OpCalcNormL2Tensor1
Get norm of input MatrixDouble for Tensor1.
Definition: NormsOperators.hpp:44
EshelbianCore::exponentBase
static double exponentBase
Definition: EshelbianCore.hpp:31
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1898
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
OpSpatialRotation_domega_domega
Definition: EshelbianOperators.hpp:446
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
OpBrokenTractionBc
Definition: EshelbianOperators.hpp:334
MoFEM::DMMoFEMCreateMoFEM
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
EshelbianCore::materialOrder
int materialOrder
Definition: EshelbianCore.hpp:126
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
EshelbianCore::spaceH1Order
int spaceH1Order
Definition: EshelbianCore.hpp:125
EshelbianCore::hybridSpatialDisp
const std::string hybridSpatialDisp
Definition: EshelbianCore.hpp:106
EshelbianCore::stretchTensor
const std::string stretchTensor
Definition: EshelbianCore.hpp:109
filter_owners
static auto filter_owners(MoFEM::Interface &m_field, Range skin)
Definition: EshelbianPlasticity.cpp:153
BiLinearForm
EshelbianCore::bcSpatialRotationVecPtr
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
Definition: EshelbianCore.hpp:138
EshelbianCore::inv_f_log_e
static double inv_f_log_e(const double v)
Definition: EshelbianCore.hpp:42
MoFEM::assembleBlockMatSchur
MoFEMErrorCode assembleBlockMatSchur(MoFEM::Interface &m_field, Mat B, Mat S, std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range >> field_ents, SmartPetscObj< AO > ao)
Assemble Schur matrix.
Definition: Schur.cpp:1817
EshelbianCore::faceExchange
CommInterface::EntitiesPetscVector faceExchange
Definition: EshelbianCore.hpp:305
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::OpCalculateHVecTensorTrace
Calculate trace of vector (Hdiv/Hcurl) space.
Definition: UserDataOperators.hpp:3184
FTensor::Tensor1::normalize
Tensor1< T, Tensor_Dim > normalize()
Definition: Tensor1_value.hpp:77
EshelbianCore::symmetrySelector
static enum SymmetrySelector symmetrySelector
Definition: EshelbianCore.hpp:14
EshelbianCore::elasticBcRhs
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
Definition: EshelbianCore.hpp:92
TSElasticPostStep::postStepDestroy
static MoFEMErrorCode postStepDestroy()
Definition: TSElasticPostStep.cpp:85
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1561
EshelbianCore::a00RangeList
std::vector< boost::shared_ptr< Range > > a00RangeList
Definition: EshelbianCore.hpp:317
OpCalculateRotationAndSpatialGradient
Definition: EshelbianOperators.hpp:173
OpSpatialConsistency_dP_domega
Definition: EshelbianOperators.hpp:501
EshelbianCore::bitAdjParentMask
BitRefLevel bitAdjParentMask
bit ref level for parent parent
Definition: EshelbianCore.hpp:296
EshelbianCore::skeletonElement
const std::string skeletonElement
Definition: EshelbianCore.hpp:116
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
OpSpatialEquilibrium
Definition: EshelbianOperators.hpp:221
OpSpatialConsistency_dBubble_domega
Definition: EshelbianOperators.hpp:514
EshelbianCore::setContactElementRhsOps
MoFEMErrorCode setContactElementRhsOps(boost::shared_ptr< ContactTree > &fe_contact_tree)
Definition: EshelbianPlasticity.cpp:2653
EshelbianCore::inv_dd_f_log_e
static double inv_dd_f_log_e(const double v)
Definition: EshelbianCore.hpp:44
EshelbianCore::listTagsToTransfer
std::vector< Tag > listTagsToTransfer
list of tags to transfer to postprocessor
Definition: EshelbianCore.hpp:310
MoFEM::OpCalculateHVecTensorDivergence
Calculate divergence of tonsorial field using vectorial base.
Definition: UserDataOperators.hpp:2950
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
EshelbianCore::elasticFeRhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
Definition: EshelbianCore.hpp:89
EshelbianCore::inv_d_f_log_e
static double inv_d_f_log_e(const double v)
Definition: EshelbianCore.hpp:43
EshelbianCore::d_f_log_e
static double d_f_log_e(const double v)
Definition: EshelbianCore.hpp:40
FTensor::dd
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
EshelbianCore::f_log
static double f_log(const double v)
Definition: EshelbianCore.hpp:46
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
EshelbianCore::alphaW
double alphaW
Definition: EshelbianCore.hpp:128
MoFEM::getBlockMatPrconditionerStorageMat
boost::shared_ptr< std::vector< double > > getBlockMatPrconditionerStorageMat(Mat B)
Get the Block Storage object.
Definition: Schur.cpp:2139
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
EshelbianPlasticity::OpConstrainBoundaryL2Rhs
Definition: EshelbianContact.hpp:118
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:45
EshelbianPlasticity::LOG
@ LOG
Definition: EshelbianPlasticity.hpp:46
EshelbianCore::dynamicRelaxation
static PetscBool dynamicRelaxation
Definition: EshelbianCore.hpp:20
OpSpatialRotation_domega_du
Definition: EshelbianOperators.hpp:412
EshelbianCore::postProcessResults
MoFEMErrorCode postProcessResults(const int tag, const std::string file, Vec f_residual=PETSC_NULL, Vec var_vec=PETSC_NULL, std::vector< Tag > tags_to_transfer={})
Definition: EshelbianPlasticity.cpp:3148
NBVOLUMETET_CCG_BUBBLE
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Definition: CGGTonsorialBubbleBase.hpp:19
EshelbianCore::bubbleField
const std::string bubbleField
Definition: EshelbianCore.hpp:111
EshelbianCore::frontAdjEdges
boost::shared_ptr< Range > frontAdjEdges
Definition: EshelbianCore.hpp:286
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
EshelbianCore::dd_f_log
static double dd_f_log(const double v)
Definition: EshelbianCore.hpp:52
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
EshelbianPlasticity::RotSelector
RotSelector
Definition: EshelbianPlasticity.hpp:45
DISCONTINUOUS
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
Definition: definitions.h:101
EshelbianCore::skinElement
const std::string skinElement
Definition: EshelbianCore.hpp:115
EshelbianCore::bcSpatialTraction
boost::shared_ptr< TractionBcVec > bcSpatialTraction
Definition: EshelbianCore.hpp:139
EshelbianCore::inv_f
static boost::function< double(const double)> inv_f
Definition: EshelbianCore.hpp:35
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
mu
double mu
Definition: dynamic_first_order_con_law.cpp:97
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::DMMoFEMTSSetIFunction
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:800
EshelbianCore::hybridMaterialDisp
const std::string hybridMaterialDisp
Definition: EshelbianCore.hpp:107
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
EshelbianCore::vertexExchange
CommInterface::EntitiesPetscVector vertexExchange
Definition: EshelbianCore.hpp:307
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition: UserDataOperators.hpp:1003
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:64
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::ForcesAndSourcesCore::gaussPts
MatrixDouble gaussPts
Matrix of integration points.
Definition: ForcesAndSourcesCore.hpp:109
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
EshelbianCore::bcSpatialDispVecPtr
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
Definition: EshelbianCore.hpp:137
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
EshelbianCore::elasticBcLhs
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
Definition: EshelbianCore.hpp:91
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::edges_conn
static constexpr int edges_conn[]
Definition: EntityRefine.cpp:18
TSElasticPostStep::postStepInitialise
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
Definition: TSElasticPostStep.cpp:11
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
get_two_sides_of_crack_surface
static auto get_two_sides_of_crack_surface(MoFEM::Interface &m_field, Range crack_faces)
Definition: EshelbianPlasticity.cpp:192
EshelbianCore::volumeExchange
CommInterface::EntitiesPetscVector volumeExchange
Definition: EshelbianCore.hpp:304
EshelbianCore::setFaceElementOps
MoFEMErrorCode setFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
Definition: EshelbianPlasticity.cpp:2570
EshelbianCore::piolaStress
const std::string piolaStress
Definition: EshelbianCore.hpp:100
BdyEleOp
BoundaryEle::UserDataOperator BdyEleOp
Definition: test_broken_space.cpp:41
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
a0
constexpr double a0
Definition: hcurl_check_approx_in_2d.cpp:37
EshelbianCore::inv_d_f
static boost::function< double(const double)> inv_d_f
Definition: EshelbianCore.hpp:36
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
EshelbianCore::inv_dd_f_linear
static double inv_dd_f_linear(const double v)
Definition: EshelbianCore.hpp:73
EshelbianCore::rotSelector
static enum RotSelector rotSelector
Definition: EshelbianCore.hpp:15
OpSpatialRotation_domega_dBubble
Definition: EshelbianOperators.hpp:434
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
EshelbianCore::getBc
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
Definition: EshelbianCore.hpp:143
QUIET
@ QUIET
Definition: definitions.h:221
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::CoreInterface::check_finite_element
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
EshelbianCore::crackFaces
boost::shared_ptr< Range > crackFaces
Definition: EshelbianCore.hpp:284
EshelbianCore::contactRefinementLevels
int contactRefinementLevels
Definition: EshelbianCore.hpp:132
MoFEM::CoreInterface::set_field_order
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.
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:589
EshelbianCore::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: EshelbianCore.hpp:86
EshelbianCore::dynamicStep
static int dynamicStep
Definition: EshelbianCore.hpp:25
OpSpatialRotation
Definition: EshelbianOperators.hpp:235
MoFEM::SmartPetscObj< IS >
OpSpatialRotation_domega_dP
Definition: EshelbianOperators.hpp:423
MoFEM::ent_form_type_and_id
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
Definition: Templates.hpp:1906
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP
Definition: EshelbianContact.hpp:172
EshelbianCore::crackingOn
static PetscBool crackingOn
Definition: EshelbianCore.hpp:22
get_range_from_block
static auto get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
Definition: EshelbianPlasticity.cpp:98
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EshelbianCore::d_f_linear
static double d_f_linear(const double v)
Definition: EshelbianCore.hpp:68
EshelbianCore::f_log_e
static double f_log_e(const double v)
Definition: EshelbianCore.hpp:39
EshelbianCore::solTSStep
SmartPetscObj< Vec > solTSStep
Definition: EshelbianCore.hpp:302
EshelbianCore::dynamicTime
static double dynamicTime
Definition: EshelbianCore.hpp:23
MoFEM::DMoFEMLoopFiniteElements
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:586
OpSpatialEquilibrium_dw_dw
Definition: EshelbianOperators.hpp:363
convert.int
int
Definition: convert.py:64
MoFEM::OpCalculateVectorFieldGradientDot
Get field gradients time derivative at integration pts for scalar filed rank 0, i....
Definition: UserDataOperators.hpp:1576
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
OpSpatialConsistency_dBubble_dBubble
Definition: EshelbianOperators.hpp:473
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
OpSpatialConsistencyP
Definition: EshelbianOperators.hpp:247
EshelbianCore::dmPrjSpatial
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
Definition: EshelbianCore.hpp:98
EshelbianCore::l2UserBaseScale
static PetscBool l2UserBaseScale
Definition: EshelbianCore.hpp:27
MoFEM::CoreInterface::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.
EshelbianCore::dd_f_linear
static double dd_f_linear(const double v)
Definition: EshelbianCore.hpp:69
OpSpatialPhysical_du_dP
Definition: EshelbianOperators.hpp:376
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1291
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::PostProcBrokenMeshInMoabBaseBegin
PostProcBrokenMeshInMoabBaseBeginImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseBegin
Enable to run stack of post-processing elements. Use this to begin stack.
Definition: PostProcBrokenMeshInMoabBase.hpp:944
F
@ F
Definition: free_surface.cpp:396
EshelbianCore::d_f
static boost::function< double(const double)> d_f
Definition: EshelbianCore.hpp:33
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
EshelbianCore::inv_dd_f
static boost::function< double(const double)> inv_dd_f
Definition: EshelbianCore.hpp:37
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:709
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
OpDispBc
Definition: EshelbianOperators.hpp:295
EshelbianCore
Definition: EshelbianCore.hpp:12
EshelbianCore::edgeExchange
CommInterface::EntitiesPetscVector edgeExchange
Definition: EshelbianCore.hpp:306
EshelbianCore::inv_f_log
static double inv_f_log(const double v)
Definition: EshelbianCore.hpp:57
Generated by Doxygen 1.8.17 and hosted at University of Glasgow