#include <users_modules/eshelbian_plasticty/src/EshelbianPlasticity.hpp>

Inheritance diagram for EshelbianPlasticity::EshelbianCore:

Collaboration diagram for EshelbianPlasticity::EshelbianCore:

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	addFields (const EntityHandle meshset=0)

MoFEMErrorCode	addVolumeFiniteElement (const EntityHandle meshset=0)

MoFEMErrorCode	addBoundaryFiniteElement (const EntityHandle meshset=0)

MoFEMErrorCode	addDMs (const BitRefLevel bit=BitRefLevel().set(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_Hencky (double E, double nu)

MoFEMErrorCode	setBaseVolumeElementOps (const int tag, const bool do_rhs, const bool do_lhs, 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	setContactElementOps (boost::shared_ptr< ContactTree > &fe_contact_tree, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)

MoFEMErrorCode	setElasticElementOps (const int tag)

MoFEMErrorCode	setElasticElementToTs (DM dm)

MoFEMErrorCode	solveElastic (TS ts, Mat m, Vec f, Vec x)

MoFEMErrorCode	postProcessResults (const int tag, const std::string file)

MoFEMErrorCode	gettingNorms ()
	[Getting norms] More...

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 refernce 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 (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< FaceElementForcesAndSourcesCore >	contactRhs

boost::shared_ptr< ContactTree >	contactTreeRhs
	Make a contact treee. More...

SmartPetscObj< DM >	dM
	Coupled problem all fields. More...

SmartPetscObj< DM >	dmElastic
	Elastic problem. More...

SmartPetscObj< DM >	dmPrjSpatial
	Projection spatial displacement. More...

const std::string	piolaStress

const std::string	eshelbyStress

const std::string	spatialL2Disp

const std::string	spatialH1Disp

const std::string	materialL2Disp

const std::string	stretchTensor

const std::string	rotAxis

const std::string	materialGradient

const std::string	tauField

const std::string	lambdaField

const std::string	bubbleField

const std::string	elementVolumeName

const std::string	naturalBcElement

const std::string	essentialBcElement

const std::string	skinElement

const std::string	contactElement

int	spaceOrder

int	spaceH1Order

int	materialOrder

double	alphaU

double	alphaW

double	alphaRho

double	precEps

double	precEpsOmega

double	precEpsW

boost::shared_ptr< BcDispVec >	bcSpatialDispVecPtr

boost::shared_ptr< BcRotVec >	bcSpatialRotationVecPtr

boost::shared_ptr< TractionBcVec >	bcSpatialTraction

boost::shared_ptr< TractionFreeBc >	bcSpatialFreeTraction

Static Public Attributes
static enum RotSelector	rotSelector

static enum RotSelector	gradApperoximator

static enum StretchSelector	stretchSelector

static double	exponentBase = exp(1)

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

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

Detailed Description

Examples: ep.cpp.

Definition at line 823 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

◆ EshelbianCore()

EshelbianPlasticity::EshelbianCore::EshelbianCore ( MoFEM::Interface & m_field )

Examples: EshelbianPlasticity.cpp.

Definition at line 91 of file EshelbianPlasticity.cpp.

    : mField(m_field), piolaStress("P"), eshelbyStress("S"),
      spatialL2Disp("wL2"), spatialH1Disp("wH1"), materialL2Disp("W"),
      stretchTensor("u"), rotAxis("omega"), materialGradient("G"),
      tauField("TAU"), lambdaField("LAMBDA"), bubbleField("BUBBLE"),
      elementVolumeName("EP"), naturalBcElement("NATURAL_BC"),
      essentialBcElement("ESSENTIAL_BC"), skinElement("SKIN_ELEMENT"),
      contactElement("CONTACT_ELEMENT") {
 
  ierr = getOptions();
  CHKERRABORT(PETSC_COMM_WORLD, ierr);
}

◆ ~EshelbianCore()

EshelbianPlasticity::EshelbianCore::~EshelbianCore ( )

virtualdefault

Member Function Documentation

◆ addBoundaryFiniteElement()

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

Examples: ep.cpp.

Definition at line 396 of file EshelbianPlasticity.cpp.

                                                                  {
  MoFEMFunctionBegin;
 
  auto bc_elements_add_to_range = [&](const std::string block_name, Range &r) {
    MoFEMFunctionBegin;
 
    auto mesh_mng = mField.getInterface<MeshsetsManager>();
    auto bcs = mesh_mng->getCubitMeshsetPtr(
 
        std::regex((boost::format("%s(.*)") % block_name).str())
 
    );
 
    for (auto bc : bcs) {
      Range faces;
      CHKERR bc->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                                 true);
      r.merge(faces);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  // set finite element fields
  auto add_field_to_fe = [this](const std::string fe,
                                const std::string field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
    MoFEMFunctionReturn(0);
  };
 
  Range natural_bc_elements;
  if (bcSpatialDispVecPtr) {
    for (auto &v : *bcSpatialDispVecPtr) {
      natural_bc_elements.merge(v.faces);
    }
  }
  if (bcSpatialRotationVecPtr) {
    for (auto &v : *bcSpatialRotationVecPtr) {
      natural_bc_elements.merge(v.faces);
    }
  }
  Range essentail_bc_elements;
  if (bcSpatialTraction) {
    for (auto &v : *bcSpatialTraction) {
      essentail_bc_elements.merge(v.faces);
    }
  }
 
  CHKERR mField.add_finite_element(naturalBcElement, MF_ZERO);
  CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
                                                   naturalBcElement);
  CHKERR add_field_to_fe(naturalBcElement, piolaStress);
  // CHKERR add_field_to_fe(naturalBcElement, eshelbyStress);
  CHKERR mField.build_finite_elements(naturalBcElement);
 
  CHKERR mField.add_finite_element(essentialBcElement, MF_ZERO);
  CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
                                                   essentialBcElement);
  CHKERR add_field_to_fe(essentialBcElement, piolaStress);
  // CHKERR add_field_to_fe(essentialBcElement, eshelbyStress);
  CHKERR mField.build_finite_elements(essentialBcElement);
 
  auto get_skin = [&]() {
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
    Skinner skin(&mField.get_moab());
    Range skin_ents;
    CHKERR skin.find_skin(0, body_ents, false, skin_ents);
    return skin_ents;
  };
 
  auto filter_true_skin = [&](auto &&skin) {
    Range boundary_ents;
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &boundary_ents);
    return boundary_ents;
  };
 
  auto skin = filter_true_skin(get_skin());
 
  CHKERR mField.add_finite_element(skinElement, MF_ZERO);
  CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
  CHKERR add_field_to_fe(skinElement, piolaStress);
  CHKERR add_field_to_fe(skinElement, spatialH1Disp);
  // CHKERR add_field_to_fe(skinElement, eshelbyStress);
  CHKERR mField.build_finite_elements(skinElement);
 
  Range contact_range;
  CHKERR bc_elements_add_to_range("CONTACT", contact_range);
  MOFEM_LOG("EP", Sev::inform) << "Contact elements " << contact_range.size();
 
  CHKERR mField.add_finite_element(contactElement, MF_ZERO);
  CHKERR mField.add_ents_to_finite_element_by_type(contact_range, MBTRI,
                                                   contactElement);
  CHKERR add_field_to_fe(contactElement, piolaStress);
  CHKERR mField.build_finite_elements(contactElement);
 
  MoFEMFunctionReturn(0);
}

◆ addDMs()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addDMs ( const BitRefLevel bit = BitRefLevel().set(0) )

Examples: ep.cpp.

Definition at line 501 of file EshelbianPlasticity.cpp.

                                                          {
  MoFEMFunctionBegin;
 
  // find adjacencies between finite elements and dofs
  CHKERR mField.build_adjacencies(bit, QUIET);
 
  // Create coupled problem
  dM = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
                            BitRefLevel().set());
  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
  CHKERR DMMoFEMAddElement(dM, elementVolumeName);
  CHKERR DMMoFEMAddElement(dM, naturalBcElement);
  CHKERR DMMoFEMAddElement(dM, essentialBcElement);
  CHKERR DMMoFEMAddElement(dM, contactElement);
  CHKERR DMMoFEMAddElement(dM, skinElement);
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
  CHKERR DMSetUp(dM);
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
 
  auto remove_dofs_on_essential_spatial_stress_boundary =
      [&](const std::string prb_name) {
        MoFEMFunctionBegin;
        for (int d : {0, 1, 2})
          CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
              prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, 0,
              MAX_DOFS_ON_ENTITY, NOISY, true);
        MoFEMFunctionReturn(0);
      };
  CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
 
  // Create elastic sub-problem
  dmElastic = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
  CHKERR DMMoFEMSetDestroyProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, piolaStress.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp.c_str());
  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
  CHKERR DMMoFEMAddElement(dmElastic, essentialBcElement);
  CHKERR DMMoFEMAddElement(dmElastic, contactElement);
  CHKERR DMMoFEMAddElement(dmElastic, skinElement);
  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
  CHKERR DMSetUp(dmElastic);
 
  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
      "ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
      "ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
      "ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
      "ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
      "ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
      "ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
      "ELASTIC_PROBLEM", bubbleField, bubbleField);
  // CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
  //     "ELASTIC_PROBLEM", piolaStress, piolaStress);
  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
      "ELASTIC_PROBLEM", bubbleField, piolaStress);
  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
      "ELASTIC_PROBLEM", piolaStress, bubbleField);
  {
    PetscSection section;
    CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
                                                           &section);
    CHKERR DMSetSection(dmElastic, section);
    CHKERR DMSetGlobalSection(dmElastic, section);
    CHKERR PetscSectionDestroy(&section);
  }
 
  dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
  CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp.c_str());
  CHKERR DMMoFEMAddElement(dmPrjSpatial, elementVolumeName);
  CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
  CHKERR DMSetUp(dmPrjSpatial);
 
  CHKERR mField.getInterface<BcManager>()
      ->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
          "PROJECT_SPATIAL", spatialH1Disp, true, false);
 
  MoFEMFunctionReturn(0);
}

◆ addFields()

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

Examples: ep.cpp.

Definition at line 223 of file EshelbianPlasticity.cpp.

                                                                  {
  MoFEMFunctionBegin;
 
  Range tets;
  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
  Range tets_skin_part;
  Skinner skin(&mField.get_moab());
  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  Range tets_skin;
  CHKERR pcomm->filter_pstatus(tets_skin_part,
                               PSTATUS_SHARED | PSTATUS_MULTISHARED,
                               PSTATUS_NOT, -1, &tets_skin);
 
  if (bcSpatialDispVecPtr)
    for (auto &v : *bcSpatialDispVecPtr) {
      tets_skin = subtract(tets_skin, v.faces);
    }
  if (bcSpatialRotationVecPtr)
    for (auto &v : *bcSpatialRotationVecPtr) {
      tets_skin = subtract(tets_skin, v.faces);
    }
  if (bcSpatialTraction)
    for (auto &v : *bcSpatialTraction) {
      tets_skin = subtract(tets_skin, v.faces);
    }
 
  auto subtract_faces_where_displacements_are_applied =
      [&](const std::string block_name) {
        MoFEMFunctionBegin;
        auto msh_mng = mField.getInterface<MeshsetsManager>();
        auto reg_exp = std::regex((boost::format("%s(.*)") % block_name).str());
        for (auto m : msh_mng->getCubitMeshsetPtr(reg_exp)) {
          Range faces;
          CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(),
                                                    SPACE_DIM - 1, faces, true);
          tets_skin = subtract(tets_skin, faces);
          MOFEM_LOG("EP", Sev::inform)
              << "Subtracting " << m->getName() << " " << faces.size();
        }
        MoFEMFunctionReturn(0);
      };
  CHKERR subtract_faces_where_displacements_are_applied("CONTACT");
 
  Range faces;
  CHKERR mField.get_moab().get_adjacencies(tets, 2, true, faces,
                                           moab::Interface::UNION);
  Range faces_not_on_the_skin = subtract(faces, tets_skin);
 
  auto add_hdiv_field = [&](const std::string field_name, const int order,
                            const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, HDIV, DEMKOWICZ_JACOBI_BASE, dim,
                            MB_TAG_SPARSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    CHKERR mField.set_field_order(faces_not_on_the_skin, field_name, order);
    CHKERR mField.set_field_order(tets_skin, field_name, 0);
    MoFEMFunctionReturn(0);
  };
 
  auto add_l2_field = [this, meshset](const std::string field_name,
                                      const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_h1_field = [this, meshset](const std::string field_name,
                                      const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
    CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_bubble_field = [this, meshset](const std::string field_name,
                                          const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
                            MF_ZERO);
    // Modify field
    auto field_ptr = mField.get_field_structure(field_name);
    auto field_order_table =
        const_cast<Field *>(field_ptr)->getFieldOrderTable();
    auto get_cgg_bubble_order_zero = [](int p) { return 0; };
    auto get_cgg_bubble_order_tet = [](int p) {
      return NBVOLUMETET_CCG_BUBBLE(p);
    };
    field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
    field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
    field_order_table[MBTRI] = get_cgg_bubble_order_zero;
    field_order_table[MBTET] = get_cgg_bubble_order_tet;
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  // spatial fields
  CHKERR add_hdiv_field(piolaStress, spaceOrder, 3);
  CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
  CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
  CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
  CHKERR add_l2_field(stretchTensor, spaceOrder, 6);
 
  // material fields
  // CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
  // CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
  // CHKERR add_l2_field(materialL2Disp, materialOrder - 1, 3);
  // CHKERR add_l2_field(tauField, materialOrder - 1, 1);
  // CHKERR add_l2_field(lambdaField, materialOrder - 1, 1);
 
  // Add history filedes
  // CHKERR add_l2_field(materialGradient + "0", materialOrder - 1, 9);
 
  // spatial displacement
  CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
 
  CHKERR mField.build_fields();
 
  MoFEMFunctionReturn(0);
}

◆ addMaterial_Hencky()

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

Examples: ep.cpp.

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

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

◆ addMaterial_HMHHStVenantKirchhoff()

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

Examples: ep.cpp.

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

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

◆ addMaterial_HMHMooneyRivlin()

MoFEMErrorCode EshelbianPlasticity::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 455 of file EshelbianADOL-C.cpp.

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

◆ addVolumeFiniteElement()

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

Examples: ep.cpp.

Definition at line 357 of file EshelbianPlasticity.cpp.

                                                                {
  MoFEMFunctionBegin;
 
  // set finite element fields
  auto add_field_to_fe = [this](const std::string fe,
                                const std::string field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
    MoFEMFunctionReturn(0);
  };
 
  if (!mField.check_finite_element(elementVolumeName)) {
    CHKERR mField.add_finite_element(elementVolumeName, MF_ZERO);
    CHKERR mField.add_ents_to_finite_element_by_type(meshset, MBTET,
                                                     elementVolumeName);
 
    CHKERR add_field_to_fe(elementVolumeName, piolaStress);
    CHKERR add_field_to_fe(elementVolumeName, bubbleField);
    // CHKERR add_field_to_fe(elementVolumeName, eshelbyStress);
    CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
    CHKERR add_field_to_fe(elementVolumeName, rotAxis);
    CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
    CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
    CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
    // CHKERR add_field_to_fe(elementVolumeName, materialGradient);
    // CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
    //                                                    materialGradient +
    //                                                    "0");
  }
 
  // build finite elements data structures
  CHKERR mField.build_finite_elements(elementVolumeName);
 
  MoFEMFunctionReturn(0);
}

◆ d_f_linear()

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

inlinestatic

Definition at line 859 of file EshelbianPlasticity.hpp.

859{ return 1; }

◆ d_f_log()

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

inlinestatic

Examples: EshelbianPlasticity.cpp.

Definition at line 840 of file EshelbianPlasticity.hpp.

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

◆ dd_f_linear()

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

inlinestatic

Definition at line 860 of file EshelbianPlasticity.hpp.

860{ return 0; }

◆ dd_f_log()

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

inlinestatic

Examples: EshelbianPlasticity.cpp.

Definition at line 843 of file EshelbianPlasticity.hpp.

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

◆ f_linear()

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

inlinestatic

Definition at line 858 of file EshelbianPlasticity.hpp.

858{ return v + 1; }

◆ f_log()

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

inlinestatic

Examples: EshelbianPlasticity.cpp.

Definition at line 837 of file EshelbianPlasticity.hpp.

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

◆ getBc()

template<typename BC >

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

inline

Definition at line 930 of file EshelbianPlasticity.hpp.

                                                                            {
    MoFEMFunctionBegin;
    for (auto it :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % block_name).str()
 
                 ))
 
    ) {
      std::vector<double> block_attributes;
      CHKERR it->getAttributes(block_attributes);
      if (block_attributes.size() != nb_attributes) {
        SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "In block %s expected %d attributes, but given %d",
                 it->getName().c_str(), nb_attributes, block_attributes.size());
      }
      Range faces;
      CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                                 true);
      bc_vec_ptr->emplace_back(it->getName(), block_attributes, faces);
    }
    MoFEMFunctionReturn(0);
  }

◆ getOptions()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getOptions ( )

Definition at line 106 of file EshelbianPlasticity.cpp.

                                         {
  MoFEMFunctionBegin;
  const char *list_rots[] = {"small", "moderate", "large"};
  PetscInt choice_rot = EshelbianCore::rotSelector;
  PetscInt choice_grad = EshelbianCore::gradApperoximator;
 
  const char *list_stretches[] = {"linear", "log"};
  PetscInt choice_stretch = StretchSelector::LOG;
 
  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
                           "none");
 
  spaceOrder = 2;
  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
                         spaceOrder, &spaceOrder, PETSC_NULL);
  spaceH1Order = -1;
  CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
                         spaceH1Order, &spaceH1Order, PETSC_NULL);
 
  materialOrder = 1;
  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
                         "", materialOrder, &materialOrder, PETSC_NULL);
 
  alphaU = 0;
  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
                            &alphaU, PETSC_NULL);
 
  alphaW = 0;
  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
                            &alphaW, PETSC_NULL);
 
  alphaRho = 0;
  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
                            &alphaRho, PETSC_NULL);
 
  precEps = 0;
  CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
                            precEps, &precEps, PETSC_NULL);
 
  precEpsOmega = 0;
  CHKERR PetscOptionsScalar("-preconditioner_eps_omega", "preconditioner_eps",
                            "", precEpsOmega, &precEpsOmega, PETSC_NULL);
  precEpsW = 0;
  CHKERR PetscOptionsScalar("-preconditioner_eps_w", "preconditioner_eps", "",
                            precEpsW, &precEpsW, PETSC_NULL);
 
  CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
                           LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
                           PETSC_NULL);
  CHKERR PetscOptionsEList("-grad", "gradient of defomation approximator", "",
                           list_rots, LARGE_ROT + 1, list_rots[choice_grad],
                           &choice_grad, PETSC_NULL);
 
  CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
                            &EshelbianCore::exponentBase, PETSC_NULL);
  CHKERR PetscOptionsEList(
      "-stretches", "stretches", "", list_stretches, StretchSelector::LOG + 1,
      list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
 
  ierr = PetscOptionsEnd();
  CHKERRG(ierr);
 
  EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
  EshelbianCore::gradApperoximator = static_cast<RotSelector>(choice_grad);
  EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
 
  switch (EshelbianCore::stretchSelector) {
  case StretchSelector::LINEAR:
    EshelbianCore::f = f_linear;
    EshelbianCore::d_f = d_f_linear;
    EshelbianCore::dd_f = dd_f_linear;
    EshelbianCore::inv_f = inv_f_linear;
    EshelbianCore::inv_d_f = inv_d_f_linear;
    EshelbianCore::inv_dd_f = inv_dd_f_linear;
    break;
  case StretchSelector::LOG:
    EshelbianCore::f = f_log;
    EshelbianCore::d_f = d_f_log;
    EshelbianCore::dd_f = dd_f_log;
    EshelbianCore::inv_f = inv_f_log;
    EshelbianCore::inv_d_f = inv_d_f_log;
    EshelbianCore::inv_dd_f = inv_dd_f_log;
    break;
  default:
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
    break;
  };
 
  precEpsOmega += precEps;
  precEpsW += precEps;
 
  MOFEM_LOG("EP", Sev::inform) << "spaceOrder " << spaceOrder;
  MOFEM_LOG("EP", Sev::inform) << "spaceH1Order " << spaceH1Order;
  MOFEM_LOG("EP", Sev::inform) << "materialOrder " << materialOrder;
  MOFEM_LOG("EP", Sev::inform) << "alphaU " << alphaU;
  MOFEM_LOG("EP", Sev::inform) << "alphaW " << alphaW;
  MOFEM_LOG("EP", Sev::inform) << "alphaRho " << alphaRho;
  MOFEM_LOG("EP", Sev::inform) << "precEps " << precEps;
  MOFEM_LOG("EP", Sev::inform) << "precEpsOmega " << precEpsOmega;
  MOFEM_LOG("EP", Sev::inform) << "precEpsW " << precEpsW;
  MOFEM_LOG("EP", Sev::inform)
      << "Rotations " << list_rots[EshelbianCore::rotSelector];
  MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
                               << list_rots[EshelbianCore::gradApperoximator];
  if (exponentBase != exp(1))
    MOFEM_LOG("EP", Sev::inform)
        << "Base exponent " << EshelbianCore::exponentBase;
  else
    MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
  MOFEM_LOG("EP", Sev::inform) << "Stretch " << list_stretches[choice_stretch];
 
  if (spaceH1Order == -1)
    spaceH1Order = spaceOrder;
 
  MoFEMFunctionReturn(0);
}

◆ getSpatialDispBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialDispBc ( )

[Getting norms]

Examples: ep.cpp.

Definition at line 1530 of file EshelbianPlasticity.cpp.

                                               {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      "", piolaStress, false, false);
 
  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto disp_bc = bc.second->dispBcPtr) {
 
      MOFEM_LOG("EP", Sev::noisy) << *disp_bc;
 
      std::vector<double> block_attributes(6, 0.);
      if (disp_bc->data.flag1 == 1) {
        block_attributes[0] = disp_bc->data.value1;
        block_attributes[3] = 1;
      }
      if (disp_bc->data.flag2 == 1) {
        block_attributes[1] = disp_bc->data.value2;
        block_attributes[4] = 1;
      }
      if (disp_bc->data.flag3 == 1) {
        block_attributes[2] = disp_bc->data.value3;
        block_attributes[5] = 1;
      }
      auto faces = bc.second->bcEnts.subset_by_dimension(2);
      bcSpatialDispVecPtr->emplace_back(bc.first, block_attributes, faces);
    }
  }
 
  // old way of naming blocksets for displacement BCs
  CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
 
  MoFEMFunctionReturn(0);
}

◆ getSpatialRotationBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialRotationBc ( )

inline

Examples: ep.cpp.

Definition at line 958 of file EshelbianPlasticity.hpp.

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

◆ getSpatialTractionBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialTractionBc ( )

Examples: ep.cpp.

Definition at line 1568 of file EshelbianPlasticity.cpp.

                                                   {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
                                                          false, false);
 
  bcSpatialTraction = boost::make_shared<TractionBcVec>();
 
  for (auto bc : bc_mng->getBcMapByBlockName()) {
    if (auto force_bc = bc.second->forceBcPtr) {
      std::vector<double> block_attributes(6, 0.);
      block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
      block_attributes[3] = 1;
      block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
      block_attributes[4] = 1;
      block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
      block_attributes[5] = 1;
      auto faces = bc.second->bcEnts.subset_by_dimension(2);
      bcSpatialTraction->emplace_back(bc.first, block_attributes, faces);
    }
  }
 
  CHKERR getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
 
  MoFEMFunctionReturn(0);
}

◆ getSpatialTractionFreeBc()

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

inline

Examples: ep.cpp.

Definition at line 980 of file EshelbianPlasticity.hpp.

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

◆ gettingNorms()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::gettingNorms ( )

[Getting norms]

Definition at line 1469 of file EshelbianPlasticity.cpp.

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

◆ getTractionFreeBc()

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

Remove all, but entities where kinematic constrains are applied.

Parameters

meshset
bc_ptr
disp_block_set_name
rot_block_set_name
contact_set_name

Returns: MoFEMErrorCode

Definition at line 641 of file EshelbianPlasticity.cpp.

                                                                   {
  MoFEMFunctionBegin;
 
  // get skin from all tets
  Range tets;
  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
  Range tets_skin_part;
  Skinner skin(&mField.get_moab());
  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  Range tets_skin;
  CHKERR pcomm->filter_pstatus(tets_skin_part,
                               PSTATUS_SHARED | PSTATUS_MULTISHARED,
                               PSTATUS_NOT, -1, &tets_skin);
 
  bc_ptr->resize(3);
  for (int dd = 0; dd != 3; ++dd)
    (*bc_ptr)[dd] = tets_skin;
 
  if (bcSpatialDispVecPtr)
    for (auto &v : *bcSpatialDispVecPtr) {
      if (v.flags[0])
        (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      if (v.flags[1])
        (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      if (v.flags[2])
        (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  if (bcSpatialRotationVecPtr)
    for (auto &v : *bcSpatialRotationVecPtr) {
      (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
      (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
      (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
    }
 
  // remove contact
  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
           std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
    Range faces;
    CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                              true);
    (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
    (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
    (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
  }
 
  // for (int dd = 0; dd != 3; ++dd) {
  //   EntityHandle meshset;
  //   CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
  //   CHKERR mField.get_moab().add_entities(meshset, (*bc_ptr)[dd]);
  //   std::string file_name = disp_block_set_name +
  //       "_traction_free_bc_" + boost::lexical_cast<std::string>(dd) + ".vtk";
  //   CHKERR mField.get_moab().write_file(file_name.c_str(), " VTK ", "",
  //                                       &meshset, 1);
  //   CHKERR mField.get_moab().delete_entities(&meshset, 1);
  // }
 
  MoFEMFunctionReturn(0);
}

◆ inv_d_f_linear()

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

inlinestatic

Definition at line 863 of file EshelbianPlasticity.hpp.

863{ return 0; }

◆ inv_d_f_log()

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

inlinestatic

Examples: EshelbianPlasticity.cpp.

Definition at line 851 of file EshelbianPlasticity.hpp.

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

◆ inv_dd_f_linear()

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

inlinestatic

Definition at line 864 of file EshelbianPlasticity.hpp.

864{ return 0; }

◆ inv_dd_f_log()

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

inlinestatic

Examples: EshelbianPlasticity.cpp.

Definition at line 854 of file EshelbianPlasticity.hpp.

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

◆ inv_f_linear()

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

inlinestatic

Definition at line 862 of file EshelbianPlasticity.hpp.

862{ return v - 1; }

◆ inv_f_log()

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

inlinestatic

Examples: EshelbianPlasticity.cpp.

Definition at line 848 of file EshelbianPlasticity.hpp.

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

◆ postProcessResults()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::postProcessResults	(	const int	tag,
		const std::string	file
	)

Definition at line 1401 of file EshelbianPlasticity.cpp.

                                                                       {
  MoFEMFunctionBegin;
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
  }
 
  PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore> post_proc(mField);
  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
 
  auto domain_ops = [&](auto &fe) {
    MoFEMFunctionBegin;
    fe.getUserPolynomialBase() =
        boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(fe.getOpPtrVector(),
                                                   {HDIV, H1, L2});
    fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
        piolaStress, dataAtPts->getApproxPAtPts()));
    fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
        bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
    fe.getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
        stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
        spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
    // evaluate derived quantities
    fe.getOpPtrVector().push_back(
        new OpCalculateRotationAndSpatialGradient(dataAtPts));
 
    // evaluate integration points
    fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
        tag, true, false, dataAtPts, physicalEquations));
 
    // contact
    fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
        spatialL2Disp, contact_common_data_ptr->contactDispPtr()));
 
    MoFEMFunctionReturn(0);
  };
 
  auto op_loop_side = new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
      mField, elementVolumeName, SPACE_DIM);
  CHKERR domain_ops(*(op_loop_side->getSideFEPtr()));
  post_proc.getOpPtrVector().push_back(op_loop_side);
  post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
      post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), spatialL2Disp,
      dataAtPts));
 
  // contact
  post_proc.getOpPtrVector().push_back(new OpTreeSearch(
      contactTreeRhs, contact_common_data_ptr, &post_proc.getPostProcMesh(),
      &post_proc.getMapGaussPts()));
 
  CHKERR DMoFEMLoopFiniteElements(dM, skinElement.c_str(),
                                  contactTreeRhs.get());
  CHKERR DMoFEMLoopFiniteElements(dM, skinElement.c_str(), &post_proc);
  CHKERR post_proc.writeFile(file.c_str());
  MoFEMFunctionReturn(0);
}

◆ query_interface()

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

Getting interface of core database.

Parameters

uuid	unique ID of interface
iface	returned pointer to interface

Returns: error code

Examples: EshelbianPlasticity.cpp.

Definition at line 73 of file EshelbianPlasticity.cpp.

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

◆ setBaseVolumeElementOps()

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

Definition at line 827 of file EshelbianPlasticity.cpp.

                                                            {
  MoFEMFunctionBegin;
  fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
 
  fe->getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(fe->getOpPtrVector(),
                                                 {HDIV, H1, L2});
 
  // set integration rule
  fe->getRuleHook = VolRule();
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
    dataAtPts->physicsPtr = physicalEquations;
  }
 
  // calculate fields values
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      piolaStress, dataAtPts->getApproxPAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
      bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
      piolaStress, dataAtPts->getDivPAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
      stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
 
  // velocities
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
      spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
  fe->getOpPtrVector().push_back(
      new OpCalculateTensor2SymmetricFieldValuesDot<3>(
          stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
      rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
 
  // acceleration
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
    fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
        spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
  }
 
  // H1 displacements
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
      spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
 
  // calculate other derived quantities
  fe->getOpPtrVector().push_back(
      new OpCalculateRotationAndSpatialGradient(dataAtPts));
 
  // evaluate integration points
  fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
      tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
 
  MoFEMFunctionReturn(0);
}

◆ setContactElementOps()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setContactElementOps	(	boost::shared_ptr< ContactTree > &	fe_contact_tree,
		boost::shared_ptr< FaceElementForcesAndSourcesCore > &	fe_rhs,
		boost::shared_ptr< FaceElementForcesAndSourcesCore > &	fe_lhs
	)

Contact requires that body is marked

Definition at line 1026 of file EshelbianPlasticity.cpp.

                                                              {
  MoFEMFunctionBegin;
 
  /** Contact requires that body is marked */
  auto get_body_range = [this](auto name) {
    map<int, Range> map;
 
    for (auto m_ptr :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % name).str()
 
                 ))
 
    ) {
      Range ents;
      CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
                                                            3, ents, true),
                     "by dim");
      map[m_ptr->getMeshsetId()] = ents;
    }
 
    return map;
  };
 
  auto get_map_skin = [this](auto &&map) {
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
 
    Skinner skin(&mField.get_moab());
    for(auto &m : map) {
      Range skin_faces;
      CHKERR skin.find_skin(0, m.second, false, skin_faces);
      CHK_MOAB_THROW(pcomm->filter_pstatus(skin_faces,
                                           PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                           PSTATUS_NOT, -1, nullptr),
                     "filter");
      m.second.swap(skin_faces);
    }
    return map;
  };
 
  /* The above code is written in C++ and it appears to be defining and using
  various operations on boundary elements and side elements. */
  using BoundaryEle = FaceElementForcesAndSourcesCore;
  using BoundaryEleOp = BoundaryEle::UserDataOperator;
 
  fe_rhs = boost::make_shared<BoundaryEle>(mField);
  fe_lhs = boost::make_shared<BoundaryEle>(mField);
  fe_rhs->getRuleHook = FaceRule();
  fe_lhs->getRuleHook = FaceRule();
 
  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(fe_rhs->getOpPtrVector(),
                                                     {HDIV});
  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(fe_lhs->getOpPtrVector(),
                                                     {HDIV});
 
  auto adj_cache =
      boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
 
  auto get_op_side = [&](auto disp_ptr, auto col_ent_data_ptr) {
    using SideEle = VolumeElementForcesAndSourcesCoreOnSide;
    auto op_loop_side = new OpLoopSide<SideEle>(
        mField, elementVolumeName, SPACE_DIM, Sev::noisy, adj_cache);
    op_loop_side->getSideFEPtr()->getUserPolynomialBase() =
        boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        op_loop_side->getOpPtrVector(), {HDIV, H1, L2});
    op_loop_side->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, disp_ptr));
    if (col_ent_data_ptr != nullptr) {
      op_loop_side->getOpPtrVector().push_back(
          new EshelbianPlasticity::OpLoopSideGetDataForSideEle(
              spatialL2Disp, col_ent_data_ptr));
    }
    return op_loop_side;
  };
 
  auto add_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    auto col_ent_data_ptr = boost::make_shared<EntitiesFieldData::EntData>();
    auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
    pip.push_back(get_op_side(contact_common_data_ptr->contactDispPtr(),
                              col_ent_data_ptr));
    pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
        piolaStress, contact_common_data_ptr->contactTractionPtr()));
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryLhs_dU(
        piolaStress, contact_common_data_ptr, col_ent_data_ptr,
        fe_contact_tree));
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryLhs_dP(
        piolaStress, contact_common_data_ptr, fe_contact_tree));
    MoFEMFunctionReturn(0);
  };
 
  auto add_contact_three = [&]() {
    MoFEMFunctionBegin;
    auto tree_moab_ptr = boost::make_shared<moab::Core>();
    fe_contact_tree =
        boost::make_shared<ContactTree>(mField, tree_moab_ptr, spaceOrder,
                                        get_map_skin(get_body_range("BODY")));
    auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
 
    using SideEle = VolumeElementForcesAndSourcesCoreOnSide;
    auto op_loop_side = new OpLoopSide<SideEle>(
        mField, elementVolumeName, SPACE_DIM, Sev::noisy, adj_cache);
    op_loop_side->getSideFEPtr()->getUserPolynomialBase() =
        boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        op_loop_side->getOpPtrVector(), {HDIV, H1, L2});
    op_loop_side->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>(
            spatialL2Disp, contact_common_data_ptr->contactDispPtr()));
    auto col_ent_data_ptr = boost::make_shared<EntitiesFieldData::EntData>();
    op_loop_side->getOpPtrVector().push_back(
        new EshelbianPlasticity::OpLoopSideGetDataForSideEle(spatialL2Disp,
                                                             col_ent_data_ptr));
 
    fe_contact_tree->getOpPtrVector().push_back(op_loop_side);
 
    auto u_h1_ptr = boost::make_shared<MatrixDouble>();
    fe_contact_tree->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
    fe_contact_tree->getOpPtrVector().push_back(new OpMoveNode(
        fe_contact_tree, contact_common_data_ptr, u_h1_ptr, col_ent_data_ptr));
    MoFEMFunctionReturn(0);
  };
 
  auto add_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
    auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
    pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
        piolaStress, contact_common_data_ptr->contactTractionPtr()));
    pip.push_back(
        get_op_side(contact_common_data_ptr->contactDispPtr(), nullptr));
    pip.push_back(new OpTreeSearch(fe_contact_tree, contact_common_data_ptr,
                                   nullptr, nullptr));
    pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryRhs(
        piolaStress, contact_common_data_ptr, fe_contact_tree));
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_contact_three();
  CHKERR add_ops_lhs(fe_lhs->getOpPtrVector());
  CHKERR add_ops_rhs(fe_rhs->getOpPtrVector());
 
  MoFEMFunctionReturn(0);
}

◆ setElasticElementOps()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setElasticElementOps ( const int tag )

Examples: ep.cpp.

Definition at line 1177 of file EshelbianPlasticity.cpp.

                                                                {
  MoFEMFunctionBegin;
  CHKERR setVolumeElementOps(tag, true, false, elasticFeRhs, elasticFeLhs);
  CHKERR setFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
 
  auto adj_cache =
      boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
 
  auto get_op_contact_bc = [&]() {
    using SideEle = FaceElementForcesAndSourcesCore;
    auto op_loop_side = new OpLoopSide<SideEle>(
        mField, contactElement, SPACE_DIM - 1, Sev::noisy, adj_cache);
    return op_loop_side;
  };
 
  auto op_contact_bc = get_op_contact_bc();
  elasticFeLhs->getOpPtrVector().push_back(op_contact_bc);
  CHKERR setContactElementOps(contactTreeRhs, contactRhs,
                              op_contact_bc->getSideFEPtr());
  MoFEMFunctionReturn(0);
}

◆ setElasticElementToTs()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setElasticElementToTs ( DM dm )

Examples: ep.cpp.

Definition at line 1199 of file EshelbianPlasticity.cpp.

                                                         {
  MoFEMFunctionBegin;
  boost::shared_ptr<FEMethod> null;
  boost::shared_ptr<FeTractionBc> spatial_traction_bc(
      new FeTractionBc(mField, piolaStress, bcSpatialTraction));
 
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
    
    CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
                                  null);
    CHKERR DMMoFEMTSSetI2Function(dm, elementVolumeName, elasticFeRhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Function(dm, naturalBcElement, elasticBcRhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Function(dm, skinElement, contactTreeRhs, null, null);
    CHKERR DMMoFEMTSSetI2Function(dm, contactElement, contactRhs, null, null);
    CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, null,
                                  spatial_traction_bc);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
                                  null);
 
  } else {
    CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
                                 null);
    CHKERR DMMoFEMTSSetIFunction(dm, elementVolumeName, elasticFeRhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIFunction(dm, naturalBcElement, elasticBcRhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIFunction(dm, skinElement, contactTreeRhs, null, null);
    CHKERR DMMoFEMTSSetIFunction(dm, contactElement, contactRhs, null, null);
    CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, null,
                                 spatial_traction_bc);
    CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
                                 null);
  }
  MoFEMFunctionReturn(0);
}

◆ setFaceElementOps()

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

Definition at line 998 of file EshelbianPlasticity.cpp.

                                                              {
  MoFEMFunctionBegin;
 
  fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
  fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
 
  // set integration rule
  fe_rhs->getRuleHook = FaceRule();
  fe_lhs->getRuleHook = FaceRule();
 
  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(fe_rhs->getOpPtrVector(),
                                                     {HDIV});
  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(fe_lhs->getOpPtrVector(),
                                                     {HDIV});
 
  if (add_elastic) {
    fe_rhs->getOpPtrVector().push_back(
        new OpDispBc(piolaStress, dataAtPts, bcSpatialDispVecPtr));
    fe_rhs->getOpPtrVector().push_back(
        new OpRotationBc(piolaStress, dataAtPts, bcSpatialRotationVecPtr));
  }
 
  MoFEMFunctionReturn(0);
}

◆ setVolumeElementOps()

MoFEMErrorCode EshelbianPlasticity::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
	)

Definition at line 893 of file EshelbianPlasticity.cpp.

                                                                {
  MoFEMFunctionBegin;
 
  auto time_scale = boost::make_shared<TimeScale>();
 
  // Right hand side
  CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
 
  // elastic
  if (add_elastic) {
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialEquilibrium(spatialL2Disp, dataAtPts, alphaW, alphaRho));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialRotation(rotAxis, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialPhysical(stretchTensor, dataAtPts, alphaU));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyP(piolaStress, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
 
    // Body forces
    using BodyNaturalBC =
        NaturalBC<VolumeElementForcesAndSourcesCore::UserDataOperator>::
            Assembly<PETSC>::LinearForm<GAUSS>;
    using OpBodyForce =
        BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
    CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
        fe_rhs->getOpPtrVector(), mField, "w", {time_scale}, "BODY_FORCE",
        Sev::inform);
  }
 
  // Left hand side
  CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
 
  // elastic
  if (add_elastic) {
 
    const bool symmetric_system =
        (gradApperoximator <= MODERATE_ROT && rotSelector == SMALL_ROT);
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
        stretchTensor, stretchTensor, dataAtPts, alphaU));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
        stretchTensor, piolaStress, dataAtPts, true));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
        stretchTensor, bubbleField, dataAtPts, true));
    if (!symmetric_system) {
      fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
          stretchTensor, rotAxis, dataAtPts, false));
    }
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
        spatialL2Disp, piolaStress, dataAtPts, true));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
        spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
        piolaStress, rotAxis, dataAtPts, symmetric_system));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
        bubbleField, rotAxis, dataAtPts, symmetric_system));
 
    if (!symmetric_system) {
      fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
          rotAxis, piolaStress, dataAtPts, false));
      fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
          rotAxis, bubbleField, dataAtPts, false));
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialRotation_domega_domega(rotAxis, rotAxis, dataAtPts));
    }
 
    // Stabilisation
    if constexpr (A == AssemblyType::SCHUR) {
      // Note that we assemble to AMat, however Assembly<SCHUR> assemble by
      // default to PMat
      using OpMassStab =
          FormsIntegrators<VolUserDataOperatorStabAssembly>::Assembly<
              SCHUR>::BiLinearForm<GAUSS>::OpMass<1, SPACE_DIM>;
      if (precEpsOmega > std::numeric_limits<double>::epsilon()) {
        fe_lhs->getOpPtrVector().push_back(
            new OpMassStab(rotAxis, rotAxis, [this](double, double, double) {
              return precEpsOmega;
            }));
      }
      if (std::abs(alphaRho + alphaW) <
          std::numeric_limits<double>::epsilon()) {
 
        fe_lhs->getOpPtrVector().push_back(new OpMassStab(
            spatialL2Disp, spatialL2Disp,
            [this](double, double, double) { return precEpsW; }));
      }
    }
 
    if (add_material) {
    }
  }
 
  MoFEMFunctionReturn(0);
}

◆ solveElastic()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::solveElastic	(	TS	ts,
		Mat	m,
		Vec	f,
		Vec	x
	)

Examples: ep.cpp.

Definition at line 1245 of file EshelbianPlasticity.cpp.

                                                                     {
  MoFEMFunctionBegin;
 
#ifdef PYTHON_SFD
 
  boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
 
  auto file_exists = [](std::string myfile) {
    std::ifstream file(myfile.c_str());
    if (file) {
      return true;
    }
    return false;
  };
 
  if (file_exists("sdf.py")) {
    MOFEM_LOG("EP", Sev::inform) << "sdf.py file found";
    sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
    CHKERR sdf_python_ptr->sdfInit("sdf.py");
    ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
  } else {
    MOFEM_LOG("EP", Sev::warning) << "sdf.py file NOT found";
  }
#else
#endif
 
  boost::shared_ptr<TsCtx> ts_ctx;
  CHKERR DMMoFEMGetTsCtx(dmElastic, ts_ctx);
 
  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
 
  boost::shared_ptr<FEMethod> monitor_ptr(new EshelbianMonitor(*this));
  ts_ctx->getLoopsMonitor().push_back(
      TsCtx::PairNameFEMethodPtr(elementVolumeName, monitor_ptr));
 
  CHKERR TSAppendOptionsPrefix(ts, "elastic_");
  CHKERR TSSetFromOptions(ts);
 
  CHKERR TSSetDM(ts, dmElastic);
 
  SNES snes;
  CHKERR TSGetSNES(ts, &snes);
 
  PetscViewerAndFormat *vf;
  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
                                    PETSC_VIEWER_DEFAULT, &vf);
  CHKERR SNESMonitorSet(
      snes,
      (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
      vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
 
  PetscSection section;
  CHKERR DMGetSection(dmElastic, &section);
  int num_fields;
  CHKERR PetscSectionGetNumFields(section, &num_fields);
  for (int ff = 0; ff != num_fields; ff++) {
    const char *field_name;
    CHKERR PetscSectionGetFieldName(section, ff, &field_name);
    MOFEM_LOG_C("EP", Sev::inform, "Field %d name %s", ff, field_name);
  }
 
  CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
 
  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
 
  // Adding field split solver
  boost::shared_ptr<SetUpSchur> schur_ptr;
 
  if constexpr (A == AssemblyType::SCHUR) {
 
    auto p = matDuplicate(m, MAT_DO_NOT_COPY_VALUES);
 
    auto ts_ctx_ptr = getDMTsCtx(dmElastic);
    ts_ctx_ptr->zeroMatrix = false;
 
    // If density is larger than zero, use dynamic time solver
    if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
      CHKERR TSSetI2Function(ts, f, PETSC_NULL, PETSC_NULL);
      CHKERR TSSetI2Jacobian(ts, m, p, PETSC_NULL, PETSC_NULL);
    } else {
      CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
      CHKERR TSSetIJacobian(ts, m, p, PETSC_NULL, PETSC_NULL);
    }
 
    KSP ksp;
    CHKERR SNESGetKSP(snes, &ksp);
    PC pc;
    CHKERR KSPGetPC(ksp, &pc);
    schur_ptr = SetUpSchur::createSetUpSchur(
        mField, SmartPetscObj<Mat>(m, true), p, &*this);
    CHKERR schur_ptr->setUp(ksp);
 
    elasticFeLhs->preProcessHook = [&]() {
      MoFEMFunctionBegin;
      *(elasticFeLhs->matAssembleSwitch) = false;
      CHKERR schur_ptr->preProc();
      MoFEMFunctionReturn(0);
    };
    elasticFeLhs->postProcessHook = [&]() {
      MoFEMFunctionBegin;
      MoFEMFunctionReturn(0);
    };
    elasticBcLhs->preProcessHook = [&]() {
      MoFEMFunctionBegin;
      MoFEMFunctionReturn(0);
    };
    elasticBcLhs->preProcessHook = [&]() {
      MoFEMFunctionBegin;
      *(elasticBcLhs->matAssembleSwitch) = false;
      CHKERR schur_ptr->postProc();
      MoFEMFunctionReturn(0);
    };
  }
 
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
    Vec xx;
    CHKERR VecDuplicate(x, &xx);
    CHKERR VecZeroEntries(xx);
    CHKERR TS2SetSolution(ts, x, xx);
    CHKERR VecDestroy(&xx);
  } else {
    CHKERR TSSetSolution(ts, x);
  }
 
  CHKERR TSSetUp(ts);
  CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
  CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
  CHKERR TSElasticPostStep::postStepInitialise(this);
  CHKERR TSSolve(ts, PETSC_NULL);
  CHKERR TSElasticPostStep::postStepDestroy();
 
  // CHKERR TSGetSNES(ts, &snes);
  int lin_solver_iterations;
  CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
  MOFEM_LOG("EP", Sev::inform)
      << "Number of linear solver iterations " << lin_solver_iterations;
 
  PetscBool test_cook_flg = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_cook", &test_cook_flg,
                             PETSC_NULL);
  if (test_cook_flg) {
    constexpr int expected_lin_solver_iterations = 11;
    // if (lin_solver_iterations != expected_lin_solver_iterations)
    //   SETERRQ2(
    //       PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
    //       "Expected number of iterations is different than expected %d !=
    //       %d", lin_solver_iterations, expected_lin_solver_iterations);
  }
 
  CHKERR gettingNorms();
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

◆ alphaRho

double EshelbianPlasticity::EshelbianCore::alphaRho

Definition at line 917 of file EshelbianPlasticity.hpp.

◆ alphaU

double EshelbianPlasticity::EshelbianCore::alphaU

Definition at line 915 of file EshelbianPlasticity.hpp.

◆ alphaW

double EshelbianPlasticity::EshelbianCore::alphaW

Definition at line 916 of file EshelbianPlasticity.hpp.

◆ bcSpatialDispVecPtr

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

Definition at line 924 of file EshelbianPlasticity.hpp.

◆ bcSpatialFreeTraction

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

Definition at line 927 of file EshelbianPlasticity.hpp.

◆ bcSpatialRotationVecPtr

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

Definition at line 925 of file EshelbianPlasticity.hpp.

◆ bcSpatialTraction

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

Definition at line 926 of file EshelbianPlasticity.hpp.

◆ bubbleField

const std::string EshelbianPlasticity::EshelbianCore::bubbleField

Definition at line 901 of file EshelbianPlasticity.hpp.

◆ contactElement

const std::string EshelbianPlasticity::EshelbianCore::contactElement

Definition at line 907 of file EshelbianPlasticity.hpp.

◆ contactRhs

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

Definition at line 884 of file EshelbianPlasticity.hpp.

◆ contactTreeRhs

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

Make a contact treee.

Definition at line 885 of file EshelbianPlasticity.hpp.

◆ d_f

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

static

Initial value:

=

EshelbianCore::d_f_log

Examples: EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 831 of file EshelbianPlasticity.hpp.

◆ dataAtPts

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

Definition at line 877 of file EshelbianPlasticity.hpp.

◆ dd_f

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

static

Initial value:

=

EshelbianCore::dd_f_log

Examples: EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 832 of file EshelbianPlasticity.hpp.

◆ dM

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dM

Coupled problem all fields.

Definition at line 887 of file EshelbianPlasticity.hpp.

◆ dmElastic

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElastic

Elastic problem.

Examples: ep.cpp.

Definition at line 888 of file EshelbianPlasticity.hpp.

◆ dmPrjSpatial

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmPrjSpatial

Projection spatial displacement.

Definition at line 889 of file EshelbianPlasticity.hpp.

◆ elasticBcLhs

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

Definition at line 882 of file EshelbianPlasticity.hpp.

◆ elasticBcRhs

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

Definition at line 883 of file EshelbianPlasticity.hpp.

◆ elasticFeLhs

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

Definition at line 881 of file EshelbianPlasticity.hpp.

◆ elasticFeRhs

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

Definition at line 880 of file EshelbianPlasticity.hpp.

◆ elementVolumeName

const std::string EshelbianPlasticity::EshelbianCore::elementVolumeName

Definition at line 903 of file EshelbianPlasticity.hpp.

◆ eshelbyStress

const std::string EshelbianPlasticity::EshelbianCore::eshelbyStress

Definition at line 892 of file EshelbianPlasticity.hpp.

◆ essentialBcElement

const std::string EshelbianPlasticity::EshelbianCore::essentialBcElement

Definition at line 905 of file EshelbianPlasticity.hpp.

◆ exponentBase

double EshelbianPlasticity::EshelbianCore::exponentBase = exp(1)

static

Examples: EshelbianPlasticity.cpp.

Definition at line 829 of file EshelbianPlasticity.hpp.

◆ f

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

static

Examples: EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 830 of file EshelbianPlasticity.hpp.

◆ gradApperoximator

enum RotSelector EshelbianPlasticity::EshelbianCore::gradApperoximator

static

Examples: EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 826 of file EshelbianPlasticity.hpp.

◆ inv_d_f

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

static

Initial value:

=

EshelbianCore::inv_d_f_log

Examples: EshelbianPlasticity.cpp.

Definition at line 834 of file EshelbianPlasticity.hpp.

◆ inv_dd_f

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

static

Initial value:

=

EshelbianCore::inv_dd_f_log

Examples: EshelbianPlasticity.cpp.

Definition at line 835 of file EshelbianPlasticity.hpp.

◆ inv_f

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

static

Initial value:

=

EshelbianCore::inv_f_log

Examples: EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 833 of file EshelbianPlasticity.hpp.

◆ lambdaField

const std::string EshelbianPlasticity::EshelbianCore::lambdaField

Definition at line 900 of file EshelbianPlasticity.hpp.

◆ materialGradient

const std::string EshelbianPlasticity::EshelbianCore::materialGradient

Definition at line 898 of file EshelbianPlasticity.hpp.

◆ materialL2Disp

const std::string EshelbianPlasticity::EshelbianCore::materialL2Disp

Definition at line 895 of file EshelbianPlasticity.hpp.

◆ materialOrder

int EshelbianPlasticity::EshelbianCore::materialOrder

Definition at line 914 of file EshelbianPlasticity.hpp.

◆ mField

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

Definition at line 875 of file EshelbianPlasticity.hpp.

◆ naturalBcElement

const std::string EshelbianPlasticity::EshelbianCore::naturalBcElement

Definition at line 904 of file EshelbianPlasticity.hpp.

◆ physicalEquations

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

Definition at line 878 of file EshelbianPlasticity.hpp.

◆ piolaStress

const std::string EshelbianPlasticity::EshelbianCore::piolaStress

Definition at line 891 of file EshelbianPlasticity.hpp.

◆ precEps

double EshelbianPlasticity::EshelbianCore::precEps

Definition at line 918 of file EshelbianPlasticity.hpp.

◆ precEpsOmega

double EshelbianPlasticity::EshelbianCore::precEpsOmega

Definition at line 919 of file EshelbianPlasticity.hpp.

◆ precEpsW

double EshelbianPlasticity::EshelbianCore::precEpsW

Definition at line 920 of file EshelbianPlasticity.hpp.

◆ rotAxis

const std::string EshelbianPlasticity::EshelbianCore::rotAxis

Definition at line 897 of file EshelbianPlasticity.hpp.

◆ rotSelector

enum RotSelector EshelbianPlasticity::EshelbianCore::rotSelector

static

Examples: EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 825 of file EshelbianPlasticity.hpp.

◆ skinElement

const std::string EshelbianPlasticity::EshelbianCore::skinElement

Definition at line 906 of file EshelbianPlasticity.hpp.

◆ spaceH1Order

int EshelbianPlasticity::EshelbianCore::spaceH1Order

Definition at line 913 of file EshelbianPlasticity.hpp.

◆ spaceOrder

int EshelbianPlasticity::EshelbianCore::spaceOrder

Definition at line 912 of file EshelbianPlasticity.hpp.

◆ spatialH1Disp

const std::string EshelbianPlasticity::EshelbianCore::spatialH1Disp

Definition at line 894 of file EshelbianPlasticity.hpp.

◆ spatialL2Disp

const std::string EshelbianPlasticity::EshelbianCore::spatialL2Disp

Definition at line 893 of file EshelbianPlasticity.hpp.

◆ stretchSelector

enum StretchSelector EshelbianPlasticity::EshelbianCore::stretchSelector

static

Examples: EshelbianOperators.cpp, and EshelbianPlasticity.cpp.

Definition at line 827 of file EshelbianPlasticity.hpp.

◆ stretchTensor

const std::string EshelbianPlasticity::EshelbianCore::stretchTensor

Definition at line 896 of file EshelbianPlasticity.hpp.

◆ tauField

const std::string EshelbianPlasticity::EshelbianCore::tauField

Definition at line 899 of file EshelbianPlasticity.hpp.

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

users_modules/eshelbian_plasticty/src/EshelbianPlasticity.hpp
users_modules/eshelbian_plasticty/src/impl/EshelbianADOL-C.cpp
users_modules/eshelbian_plasticty/src/impl/EshelbianPlasticity.cpp

Classes

Public Member Functions

Static Public Member Functions

Public Attributes

Static Public Attributes

Detailed Description

Constructor & Destructor Documentation

◆ EshelbianCore()

◆ ~EshelbianCore()

Member Function Documentation

◆ addBoundaryFiniteElement()

◆ addDMs()

◆ addFields()

◆ addMaterial_Hencky()

◆ addMaterial_HMHHStVenantKirchhoff()

◆ addMaterial_HMHMooneyRivlin()

◆ addVolumeFiniteElement()

◆ d_f_linear()

◆ d_f_log()

◆ dd_f_linear()

◆ dd_f_log()

◆ f_linear()

◆ f_log()

◆ getBc()

◆ getOptions()

◆ getSpatialDispBc()

◆ getSpatialRotationBc()

◆ getSpatialTractionBc()

◆ getSpatialTractionFreeBc()

◆ gettingNorms()

◆ getTractionFreeBc()

◆ inv_d_f_linear()

◆ inv_d_f_log()

◆ inv_dd_f_linear()

◆ inv_dd_f_log()

◆ inv_f_linear()

◆ inv_f_log()

◆ postProcessResults()

◆ query_interface()

◆ setBaseVolumeElementOps()

◆ setContactElementOps()

◆ setElasticElementOps()

◆ setElasticElementToTs()

◆ setFaceElementOps()

◆ setVolumeElementOps()

◆ solveElastic()

Member Data Documentation

◆ alphaRho

◆ alphaU

◆ alphaW

◆ bcSpatialDispVecPtr

◆ bcSpatialFreeTraction

◆ bcSpatialRotationVecPtr

◆ bcSpatialTraction

◆ bubbleField

◆ contactElement

◆ contactRhs

◆ contactTreeRhs

◆ d_f

◆ dataAtPts

◆ dd_f

◆ dM

◆ dmElastic

◆ dmPrjSpatial

◆ elasticBcLhs

◆ elasticBcRhs

◆ elasticFeLhs

◆ elasticFeRhs

◆ elementVolumeName

◆ eshelbyStress

◆ essentialBcElement

◆ exponentBase

◆ f

◆ gradApperoximator

◆ inv_d_f

◆ inv_dd_f

◆ inv_f

◆ lambdaField

◆ materialGradient

◆ materialL2Disp