EshelbianPlasticity::EshelbianCore Struct Reference

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

Inheritance diagram for EshelbianPlasticity::EshelbianCore:

Collaboration diagram for EshelbianPlasticity::EshelbianCore:

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_set_name, const int nb_attributes)
MoFEMErrorCode	getSpatialDispBc ()
MoFEMErrorCode	getSpatialRotationBc ()
MoFEMErrorCode	getSpatialTractionBc ()
MoFEMErrorCode	getTractionFreeBc (const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string disp_block_set_name, const std::string rot_block_set_name)
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	setBaseVolumeElementOps (const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe)
MoFEMErrorCode	setGenericVolumeElementOps (const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe_rhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > &fe_lhs)
MoFEMErrorCode	setGenericFaceElementOps (const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &fe_rhs, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > &fe_lhs)
MoFEMErrorCode	setElasticElementOps (const int tag)
MoFEMErrorCode	setElasticElementToTs (DM dm)
MoFEMErrorCode	setUpTSElastic (TS ts, Mat m, Vec f)
MoFEMErrorCode	solveElastic (TS ts, Vec x)
MoFEMErrorCode	postProcessResults (const int tag, const std::string file)
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

Public Attributes
MoFEM::Interface &	mField
boost::shared_ptr< DataAtIntegrationPts >	dataAtPts
boost::shared_ptr< PhysicalEquations >	physicalEquations
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > >	elasticFeRhs
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > >	elasticFeLhs
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > >	elasticBcLhs
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > >	elasticBcRhs
boost::shared_ptr< EpFEMethod >	schurAssembly
SmartPetscObj< DM >	dM
	Coupled problem all fields. More...
SmartPetscObj< DM >	dmElastic
	Elastic problem. More...
SmartPetscObj< DM >	dmElasticSchurStreach
	Sub problem of dmElastic Schur. More...
SmartPetscObj< DM >	dmElasticSchurBubble
	Sub problem of dmElastic Schur. More...
SmartPetscObj< DM >	dmElasticSchurSpatialDisp
	Sub problem of dmElastic Schur. More...
SmartPetscObj< DM >	dmElasticSchurOmega
	Sub problem of dmElastic Schur. More...
SmartPetscObj< DM >	dmMaterial
	Material problem. More...
const std::string	piolaStress
const std::string	eshelbyStress
const std::string	spatialDisp
const std::string	materialDisp
const std::string	streachTensor
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
int	spaceOrder
int	materialOrder
double	alphaU
double	alphaW
double	alphaRho
double	precEps
boost::shared_ptr< BcDispVec >	bcSpatialDispVecPtr
boost::shared_ptr< BcRotVec >	bcSpatialRotationVecPtr
boost::shared_ptr< TractionBcVec >	bcSpatialTraction
boost::shared_ptr< TractionFreeBc >	bcSpatialFreeTraction

Additional Inherited Members
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...

Detailed Description

Examples

ep.cpp.

Definition at line 1157 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

EshelbianCore()

EshelbianPlasticity::EshelbianCore::EshelbianCore

(

MoFEM::Interface &

m_field

)

Examples

EshelbianPlasticity.cpp.

Definition at line 40 of file EshelbianPlasticity.cpp.

    : mField(m_field), piolaStress("P"), eshelbyStress("S"), spatialDisp("w"),
      materialDisp("W"), streachTensor("u"), rotAxis("omega"),
      materialGradient("G"), tauField("TAU"), lambdaField("LAMBDA"),
      bubbleField("BUBBLE"), elementVolumeName("EP"),
      naturalBcElement("NATURAL_BC"), essentialBcElement("ESSENTIAL_BC") {
 
  ierr = getOptions();
  CHKERRABORT(PETSC_COMM_WORLD, ierr);
}

~EshelbianCore()

EshelbianPlasticity::EshelbianCore::~EshelbianCore ( )

virtual

Examples

EshelbianPlasticity.cpp.

Definition at line 51 of file EshelbianPlasticity.cpp.

{}

Member Function Documentation

addBoundaryFiniteElement()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addBoundaryFiniteElement

(

const EntityHandle

meshset = 0

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 254 of file EshelbianPlasticity.cpp.

                                                                  {
  MoFEMFunctionBegin;
 
  auto bc_elements_add_to_range = [&](const std::string disp_block_set_name,
                                      Range &r) {
    MoFEMFunctionBegin;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, disp_block_set_name.length(),
                                disp_block_set_name) == 0) {
        Range faces;
        CHKERR it->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;
  CHKERR bc_elements_add_to_range("SPATIAL_DISP_BC", natural_bc_elements);
  CHKERR bc_elements_add_to_range("SPATIAL_ROTATION_BC", natural_bc_elements);
  Range essentail_bc_elements;
  CHKERR bc_elements_add_to_range("SPATIAL_TRACTION_BC", essentail_bc_elements);
 
  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);
 
  MoFEMFunctionReturn(0);
}

addDMs()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addDMs

(

const BitRefLevel

bit = BitRefLevel().set(0)

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 305 of file EshelbianPlasticity.cpp.

                                                          {
  MoFEMFunctionBegin;
 
  // find adjacencies between finite elements and dofs
  CHKERR mField.build_adjacencies(bit, QUIET);
 
  // Create coupled problem
  dM = createSmartDM(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.c_str());
  CHKERR DMMoFEMAddElement(dM, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dM, essentialBcElement.c_str());
  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, 100,
              NOISY, true);
        MoFEMFunctionReturn(0);
      };
  CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
 
  // Create elastic sub-problem
  dmElastic = createSmartDM(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, streachTensor.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialDisp.c_str());
  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName.c_str());
  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dmElastic, essentialBcElement.c_str());
  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
  CHKERR DMSetUp(dmElastic);
 
  // Create elastic streach-problem
  dmElasticSchurStreach = createSmartDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElasticSchurStreach, dmElastic,
                            "ELASTIC_PROBLEM_STREACH_SCHUR");
  CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurStreach, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, piolaStress.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, bubbleField.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, rotAxis.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurStreach, spatialDisp.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurStreach, elementVolumeName.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurStreach, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurStreach, essentialBcElement.c_str());
  CHKERR DMMoFEMSetSquareProblem(dmElasticSchurStreach, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurStreach, PETSC_TRUE);
  CHKERR DMSetUp(dmElasticSchurStreach);
 
  // Create elastic bubble-problem
  dmElasticSchurBubble = createSmartDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElasticSchurBubble, dmElasticSchurStreach,
                            "ELASTIC_PROBLEM_BUBBLE_SCHUR");
  CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurBubble, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, piolaStress.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, rotAxis.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurBubble, spatialDisp.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurBubble, elementVolumeName.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurBubble, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurBubble, essentialBcElement.c_str());
  CHKERR DMMoFEMSetSquareProblem(dmElasticSchurBubble, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurBubble, PETSC_TRUE);
  CHKERR DMSetUp(dmElasticSchurBubble);
 
  // Create elastic omega-problem
  dmElasticSchurOmega = createSmartDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElasticSchurOmega, dmElasticSchurBubble,
                            "ELASTIC_PROBLEM_OMEGA_SCHUR");
  CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurOmega, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurOmega, piolaStress.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurOmega, spatialDisp.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurOmega, elementVolumeName.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurOmega, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurOmega, essentialBcElement.c_str());
  CHKERR DMMoFEMSetSquareProblem(dmElasticSchurOmega, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurOmega, PETSC_TRUE);
  CHKERR DMSetUp(dmElasticSchurOmega);
 
  // Create elastic tet_stress-problem
  dmElasticSchurSpatialDisp = createSmartDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(dmElasticSchurSpatialDisp, dmElasticSchurOmega,
                            "ELASTIC_PROBLEM_SPATIAL_DISP_SCHUR");
  CHKERR DMMoFEMSetDestroyProblem(dmElasticSchurSpatialDisp, PETSC_TRUE);
  CHKERR DMMoFEMAddSubFieldRow(dmElasticSchurSpatialDisp, piolaStress.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp,
                           elementVolumeName.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp, naturalBcElement.c_str());
  CHKERR DMMoFEMAddElement(dmElasticSchurSpatialDisp,
                           essentialBcElement.c_str());
  CHKERR DMMoFEMSetSquareProblem(dmElasticSchurSpatialDisp, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dmElasticSchurSpatialDisp, PETSC_TRUE);
  CHKERR DMSetUp(dmElasticSchurSpatialDisp);
 
  {
    PetscSection section;
    CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
                                                           &section);
    CHKERR DMSetDefaultSection(dmElastic, section);
    CHKERR DMSetDefaultGlobalSection(dmElastic, section);
    CHKERR PetscSectionDestroy(&section);
  }
 
  MoFEMFunctionReturn(0);
}

addFields()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addFields

(

const EntityHandle

meshset = 0

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 86 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);
 
  auto subtract_faces_where_displacements_are_applied =
      [&](const std::string disp_block_set_name) {
        MoFEMFunctionBegin;
        for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
          if (it->getName().compare(0, disp_block_set_name.length(),
                                    disp_block_set_name) == 0) {
            Range faces;
            CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
                                                       faces, true);
            tets_skin = subtract(tets_skin, faces);
          }
        }
        MoFEMFunctionReturn(0);
      };
  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_DISP_BC");
  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_ROTATION_BC");
  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_TRACTION_BC");
 
  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_hdiv_rt_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_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, 0);
    CHKERR mField.set_field_order(tets_skin, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_l2_field = [this, meshset](const std::string field_name,
                                      const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_h1_field = [this, meshset](const std::string field_name,
                                      const int order, const int dim) {
    MoFEMFunctionBegin;
    CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
                            MB_TAG_DENSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
    CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
    CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
    CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
    MoFEMFunctionReturn(0);
  };
 
  auto add_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(spatialDisp, spaceOrder - 1, 3);
  CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
  CHKERR add_l2_field(streachTensor, spaceOrder, 6);
 
  // material fields
  CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
  CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
  // CHKERR add_l2_field(materialDisp, 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);
 
  CHKERR mField.build_fields();
 
  MoFEMFunctionReturn(0);
}

addMaterial_HMHHStVenantKirchhoff()

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

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

                          {
  MoFEMFunctionBegin;
  physicalEquations = boost::shared_ptr<HMHStVenantKirchhoff>(
      new HMHStVenantKirchhoff(lambda, mu, sigma_y));
  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 643 of file EshelbianADOL-C.cpp.

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

addVolumeFiniteElement()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::addVolumeFiniteElement

(

const EntityHandle

meshset = 0

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 216 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, streachTensor);
    CHKERR add_field_to_fe(elementVolumeName, rotAxis);
    CHKERR add_field_to_fe(elementVolumeName, spatialDisp);
    CHKERR add_field_to_fe(elementVolumeName, streachTensor);
    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);
}

getBc()

template<typename BC >

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

inline

Definition at line 1221 of file EshelbianPlasticity.hpp.

                                                {
    MoFEMFunctionBegin;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      auto block_name = it->getName();
      if (block_name.compare(0, block_set_name.length(), block_set_name) == 0) {
        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(block_name, block_attributes, faces);
      }
    }
    MoFEMFunctionReturn(0);
  }

getOptions()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getOptions

(

)

Examples

EshelbianPlasticity.cpp.

Definition at line 53 of file EshelbianPlasticity.cpp.

                                         {
  MoFEMFunctionBegin;
  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
                           "none");
 
  spaceOrder = 1;
  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
                         spaceOrder, &spaceOrder, 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 = 1e-3;
  CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
                            precEps, &precEps, PETSC_NULL);
 
  ierr = PetscOptionsEnd();
  CHKERRG(ierr);
  MoFEMFunctionReturn(0);
}

getSpatialDispBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialDispBc ( )

inline

Examples

ep.cpp.

Definition at line 1245 of file EshelbianPlasticity.hpp.

                                           {
    bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
    return getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
  }

getSpatialRotationBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialRotationBc ( )

inline

Examples

ep.cpp.

Definition at line 1250 of file EshelbianPlasticity.hpp.

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

getSpatialTractionBc()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::getSpatialTractionBc ( )

inline

Examples

ep.cpp.

Definition at line 1255 of file EshelbianPlasticity.hpp.

                                               {
    bcSpatialTraction = boost::make_shared<TractionBcVec>();
    return getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
  }

getSpatialTractionFreeBc()

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

inline

Examples

ep.cpp.

Definition at line 1265 of file EshelbianPlasticity.hpp.

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

getTractionFreeBc()

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

Examples

EshelbianPlasticity.cpp.

Definition at line 454 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);
 
  bc_ptr->resize(3);
  for (int dd = 0; dd != 3; ++dd)
    (*bc_ptr)[dd] = tets_skin;
 
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    if (it->getName().compare(0, disp_block_set_name.length(),
                              disp_block_set_name) == 0) {
      std::vector<double> block_attributes;
      CHKERR it->getAttributes(block_attributes);
      if (block_attributes.size() != 6) {
        SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "In block %s six attributes are required for given BC "
                 "blockset (3 values + "
                 "3 flags) != %d",
                 it->getName().c_str(), block_attributes.size());
      }
      Range faces;
      CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                                 true);
      if (block_attributes[3] != 0)
        (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
      if (block_attributes[4] != 0)
        (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
      if (block_attributes[5] != 0)
        (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
    }
    if (it->getName().compare(0, rot_block_set_name.length(),
                              rot_block_set_name) == 0) {
      Range faces;
      CHKERR it->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);
}

postProcessResults()

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1471 of file EshelbianPlasticity.cpp.

                                                                       {
  MoFEMFunctionBegin;
 
  if (!dataAtPts) {
    dataAtPts =
        boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
  }
 
  PostProcVolumeOnRefinedMesh post_proc(mField);
  post_proc.getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
  post_proc.getOpPtrVector().push_back(new OpL2Transform());
 
  CHKERR post_proc.generateReferenceElementMesh();
  post_proc.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      piolaStress, dataAtPts->getApproxPAtPts()));
  post_proc.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
      bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  post_proc.getOpPtrVector().push_back(
      new OpCalculateTensor2SymmetricFieldValues<3>(
          streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  post_proc.getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
      materialGradient, dataAtPts->getBigGAtPts(), MBTET));
  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
 
  // evaluate derived quantities
  post_proc.getOpPtrVector().push_back(
      new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
 
  // evaluate integration points
  post_proc.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
      spatialDisp, tag, true, false, dataAtPts, physicalEquations));
  post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
      post_proc.postProcMesh, post_proc.mapGaussPts, spatialDisp, dataAtPts));
 
  CHKERR DMoFEMLoopFiniteElements(dM, elementVolumeName.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 20 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< EpElement< VolumeElementForcesAndSourcesCore > > &	fe
	)

Examples

EshelbianPlasticity.cpp.

Definition at line 630 of file EshelbianPlasticity.cpp.

                                                                       {
  MoFEMFunctionBegin;
  fe = boost::make_shared<EpElement<VolumeElementForcesAndSourcesCore>>(mField);
 
  fe->getUserPolynomialBase() =
      boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
  fe->getOpPtrVector().push_back(new OpL2Transform());
 
  // 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 OpCalculateHVecTensorField<3, 3>(
      eshelbyStress, dataAtPts->getApproxSigmaAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
      eshelbyStress, dataAtPts->getDivSigmaAtPts()));
  fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
      streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
      materialGradient, dataAtPts->getBigGAtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
      materialGradient + "0", dataAtPts->getBigG0AtPts(), MBTET));
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
      spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
 
  // velocities
  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
      spatialDisp, dataAtPts->getSmallWDotAtPts(), MBTET));
  fe->getOpPtrVector().push_back(
      new OpCalculateTensor2SymmetricFieldValuesDot<3>(
          streachTensor, dataAtPts->getLogStreachDotTensorAtPts(), 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>(
        spatialDisp, dataAtPts->getSmallWDotDotAtPts(), MBTET));
  }
 
  // calculate other derived quantities
  fe->getOpPtrVector().push_back(
      new OpCalculateRotationAndSpatialGradient(rotAxis, dataAtPts));
 
  // evaluate integration points
  fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
      spatialDisp, tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
 
  MoFEMFunctionReturn(0);
}

setElasticElementOps()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setElasticElementOps

(

const int

tag

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 811 of file EshelbianPlasticity.cpp.

                                                                {
  MoFEMFunctionBegin;
  CHKERR setGenericVolumeElementOps(tag, true, false, elasticFeRhs,
                                    elasticFeLhs);
  CHKERR setGenericFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
  MoFEMFunctionReturn(0);
}

setElasticElementToTs()

MoFEMErrorCode EshelbianPlasticity::EshelbianCore::setElasticElementToTs

(

DM

dm

)

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 819 of file EshelbianPlasticity.cpp.

                                                         {
  MoFEMFunctionBegin;
  boost::shared_ptr<FEMethod> null;
  boost::shared_ptr<EpElement<FeTractionBc>> spatial_traction_bc(
      new EpElement<FeTractionBc>(mField, piolaStress, bcSpatialTraction));
  schurAssembly = boost::make_shared<EpFEMethod>();
 
  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, DM_NO_ELEMENT, null, null,
                                  spatial_traction_bc);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, DM_NO_ELEMENT, null, schurAssembly, null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
                                  null);
    CHKERR DMMoFEMTSSetI2Jacobian(dm, DM_NO_ELEMENT, null, null, schurAssembly);
  } 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, DM_NO_ELEMENT, null, null,
                                 spatial_traction_bc);
    CHKERR DMMoFEMTSSetIJacobian(dm, DM_NO_ELEMENT, null, schurAssembly, null);
    CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
                                 null);
    CHKERR DMMoFEMTSSetIJacobian(dm, DM_NO_ELEMENT, null, null, schurAssembly);
  }
  MoFEMFunctionReturn(0);
}

setGenericFaceElementOps()

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

Examples

EshelbianPlasticity.cpp.

Definition at line 781 of file EshelbianPlasticity.cpp.

                                                                         {
  MoFEMFunctionBegin;
 
  fe_rhs =
      boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
  fe_lhs =
      boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
 
  // set integration rule
  fe_rhs->getRuleHook = FaceRule();
  fe_lhs->getRuleHook = FaceRule();
 
  fe_rhs->getOpPtrVector().push_back(
      new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
  fe_lhs->getOpPtrVector().push_back(
      new OpHOSetContravariantPiolaTransformOnFace3D(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);
}

setGenericVolumeElementOps()

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

Examples

EshelbianPlasticity.cpp.

Definition at line 697 of file EshelbianPlasticity.cpp.

                                                                           {
  MoFEMFunctionBegin;
 
  // Right hand side
  CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
 
  // elastic
  if (add_elastic) {
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialEquilibrium(spatialDisp, dataAtPts, alphaW, alphaRho));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialRotation(rotAxis, dataAtPts));
    fe_rhs->getOpPtrVector().push_back(
        new OpSpatialPhysical(streachTensor, 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));
  }
 
  // Left hand side
  CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
 
  // elastic
  if (add_elastic) {
 
    // Schur
    fe_lhs->getOpPtrVector().push_back(
        new OpSpatialSchurBegin(spatialDisp, dataAtPts));
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
        streachTensor, streachTensor, dataAtPts, alphaU));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
        streachTensor, bubbleField, dataAtPts, true));
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
        spatialDisp, piolaStress, dataAtPts, true));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
        spatialDisp, spatialDisp, dataAtPts, alphaW, alphaRho));
 
    fe_lhs->getOpPtrVector().push_back(
        new OpSpatialConsistency_dP_domega(piolaStress, rotAxis, dataAtPts));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
        bubbleField, rotAxis, dataAtPts));
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
        streachTensor, piolaStress, dataAtPts, true));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
        streachTensor, rotAxis, dataAtPts, false));
 
    fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
        rotAxis, piolaStress, dataAtPts, false));
    fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
        rotAxis, bubbleField, dataAtPts, false));
    fe_lhs->getOpPtrVector().push_back(
        new OpSpatialRotation_domega_domega(rotAxis, rotAxis, dataAtPts));
 
    // Schur
    dataAtPts->ooMatPtr = boost::make_shared<MatrixDouble>();
    fe_lhs->getOpPtrVector().push_back(
        new OpSpatialPreconditionMass(rotAxis, dataAtPts->ooMatPtr));
    if (alphaW < std::numeric_limits<double>::epsilon() &&
        alphaRho < std::numeric_limits<double>::epsilon()) {
      dataAtPts->wwMatPtr = boost::make_shared<MatrixDouble>();
      fe_lhs->getOpPtrVector().push_back(
          new OpSpatialPreconditionMass(spatialDisp, dataAtPts->wwMatPtr));
    } else {
      dataAtPts->wwMatPtr.reset();
    }
    fe_lhs->getOpPtrVector().push_back(
        new OpSpatialSchurEnd(spatialDisp, dataAtPts, precEps));
 
    if (add_material) {
    }
  }
 
  MoFEMFunctionReturn(0);
}

setUpTSElastic()

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

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 860 of file EshelbianPlasticity.cpp.

                                                                {
  MoFEMFunctionBegin;
  boost::shared_ptr<TsCtx> ts_ctx;
  CHKERR DMMoFEMGetTsCtx(dmElastic, ts_ctx);
 
  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
    CHKERR TSSetI2Function(ts, f, PETSC_NULL, PETSC_NULL);
    CHKERR TSSetI2Jacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
  } else {
    CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
    CHKERR TSSetIJacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
  }
 
  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
 
  auto add_schur_streach_op = [this](auto &list, SmartPetscObj<Mat> S,
                                     SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
        fe_cast->addStreachSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_streach_pre = [this](auto &list, SmartPetscObj<Mat> S,
                                      SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
        fe_cast->addStreachSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_bubble_op = [this](auto &list, SmartPetscObj<Mat> S,
                                    SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
        fe_cast->addBubbleSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_bubble_pre = [this](auto &list, SmartPetscObj<Mat> S,
                                     SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
        fe_cast->addBubbleSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_spatial_disp_op = [this](auto &list, SmartPetscObj<Mat> S,
                                          SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
        fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_spatial_disp_pre = [this](auto &list, SmartPetscObj<Mat> S,
                                           SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
        fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_omega_op = [this](auto &list, SmartPetscObj<Mat> S,
                                   SmartPetscObj<AO> aoS) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
        fe_cast->addOmegaSchurMatrix(S, aoS);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  auto add_schur_omega_pre = [this](auto &list, SmartPetscObj<Mat> S,
                                    SmartPetscObj<AO> ao) {
    MoFEMFunctionBegin;
    for (auto &fe : list)
      if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
        fe_cast->addOmegaSchurMatrix(S, ao);
      else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
    MoFEMFunctionReturn(0);
  };
 
  const MoFEM::Problem *schur_streach_prb_ptr;
  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_streach_prb_ptr);
  if (auto sub_data = schur_streach_prb_ptr->subProblemData) {
    SmartPetscObj<Mat> Suu;
    CHKERR DMCreateMatrix_MoFEM(dmElasticSchurStreach, Suu);
    SmartPetscObj<AO> aoSuu = sub_data->getSmartRowMap();
 
    CHKERR add_schur_streach_op(ts_ctx->loopsIJacobian, Suu, aoSuu);
    CHKERR add_schur_streach_pre(ts_ctx->preProcessIJacobian, Suu, aoSuu);
    CHKERR add_schur_streach_pre(ts_ctx->postProcessIJacobian, Suu, aoSuu);
 
    const MoFEM::Problem *schur_bubble_prb_ptr;
    CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_prb_ptr);
    if (auto bubble_data = schur_bubble_prb_ptr->subProblemData) {
      SmartPetscObj<Mat> SBubble;
      CHKERR DMCreateMatrix_MoFEM(dmElasticSchurBubble, SBubble);
      SmartPetscObj<AO> aoSBubble = bubble_data->getSmartRowMap();
      CHKERR add_schur_bubble_op(ts_ctx->loopsIJacobian, SBubble,
                                 aoSBubble);
      CHKERR add_schur_bubble_pre(ts_ctx->preProcessIJacobian, SBubble,
                                  aoSBubble);
      CHKERR add_schur_bubble_pre(ts_ctx->postProcessIJacobian, SBubble,
                                  aoSBubble);
 
      const MoFEM::Problem *schur_omega_prb_ptr;
      CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_prb_ptr);
      if (auto tet_stress_data = schur_omega_prb_ptr->subProblemData) {
        SmartPetscObj<Mat> SOmega;
        CHKERR DMCreateMatrix_MoFEM(dmElasticSchurOmega, SOmega);
        SmartPetscObj<AO> aoSOmega = tet_stress_data->getSmartRowMap();
        CHKERR add_schur_omega_op(ts_ctx->loopsIJacobian, SOmega,
                                  aoSOmega);
        CHKERR add_schur_omega_pre(ts_ctx->preProcessIJacobian, SOmega,
                                   aoSOmega);
        CHKERR add_schur_omega_pre(ts_ctx->postProcessIJacobian, SOmega,
                                   aoSOmega);
 
        const MoFEM::Problem *schur_spatial_disp_prb_ptr;
        CHKERR DMMoFEMGetProblemPtr(dmElasticSchurSpatialDisp,
                                    &schur_spatial_disp_prb_ptr);
        if (auto spatial_disp_data =
                schur_spatial_disp_prb_ptr->subProblemData) {
          SmartPetscObj<Mat> Sw;
          CHKERR DMCreateMatrix_MoFEM(dmElasticSchurSpatialDisp, Sw);
          SmartPetscObj<AO> aoSw = spatial_disp_data->getSmartRowMap();
          CHKERR add_schur_spatial_disp_op(ts_ctx->loopsIJacobian, Sw,
                                           aoSw);
          CHKERR add_schur_spatial_disp_pre(ts_ctx->preProcessIJacobian, Sw,
                                            aoSw);
          CHKERR add_schur_spatial_disp_pre(ts_ctx->postProcessIJacobian, Sw,
                                            aoSw);
        } else
          SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                  "Problem does not have sub-problem data");
 
      } else
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                "Problem does not have sub-problem data");
 
    } else
      SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
              "Problem does not have sub-problem data");
 
  } else
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "Problem does not have sub-problem data");
 
  struct Monitor : public FEMethod {
 
    using Ele = ForcesAndSourcesCore;
    using VolEle = VolumeElementForcesAndSourcesCore;
    using VolOp = VolumeElementForcesAndSourcesCore::UserDataOperator;
    using SetPtsData = FieldEvaluatorInterface::SetPtsData;
 
    EshelbianCore &eP;
    boost::shared_ptr<SetPtsData> dataFieldEval;
    boost::shared_ptr<VolEle> volPostProcEnergy;
    boost::shared_ptr<double> gEnergy;
 
    Monitor(EshelbianCore &ep)
        : eP(ep),
          dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
                            ->getData<VolEle>()),
          volPostProcEnergy(new VolEle(ep.mField)), gEnergy(new double) {
      ierr = ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
          dataFieldEval, "EP");
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
 
      auto no_rule = [](int, int, int) { return -1; };
 
      auto set_element_for_field_eval = [&]() {
        boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
        vol_ele->getRuleHook = no_rule;
        vol_ele->getUserPolynomialBase() =
            boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
        vol_ele->getOpPtrVector().push_back(new OpL2Transform());
 
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateHVecTensorField<3, 3>(
                eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateHTensorTensorField<3, 3>(
                eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<3>(
                eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
                MBTET));
        vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateTensor2FieldValues<3, 3>(
                eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
        vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
            eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
                                                      eP.dataAtPts));
      };
 
      auto set_element_for_post_process = [&]() {
        volPostProcEnergy->getRuleHook = VolRule();
        volPostProcEnergy->getUserPolynomialBase() =
            boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
        volPostProcEnergy->getOpPtrVector().push_back(new OpL2Transform());
 
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateHVecTensorField<3, 3>(
                eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateHTensorTensorField<3, 3>(
                eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<3>(
                eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(),
                MBTET));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValues<3>(
                eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateTensor2FieldValues<3, 3>(
                eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValues<3>(
                eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateRotationAndSpatialGradient(eP.rotAxis,
                                                      eP.dataAtPts));
        volPostProcEnergy->getOpPtrVector().push_back(
            new OpCalculateStrainEnergy(eP.spatialDisp, eP.dataAtPts, gEnergy));
      };
 
      set_element_for_field_eval();
      set_element_for_post_process();
    }
 
    MoFEMErrorCode preProcess() { return 0; }
 
    MoFEMErrorCode operator()() { return 0; }
 
    MoFEMErrorCode postProcess() {
      MoFEMFunctionBegin;
 
      // auto get_str_time = [](auto ts_t) {
      //   std::ostringstream ss;
      //   ss << boost::str(boost::format("%d") %
      //                    static_cast<int>(std::ceil(ts_t * 1e6)));
      //   std::string s = ss.str();
      //   return s;
      // };
 
      auto get_step = [](auto ts_step) {
        std::ostringstream ss;
        ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
        std::string s = ss.str();
        return s;
      };
 
      PetscViewer viewer;
      CHKERR PetscViewerBinaryOpen(
          PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
          FILE_MODE_WRITE, &viewer);
      CHKERR VecView(ts_u, viewer);
      CHKERR PetscViewerDestroy(&viewer);
 
      CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_step(ts_step) +
                                          ".h5m");
 
      // Loop boundary elements with traction boundary conditions
      *gEnergy = 0;
      CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
                                            *volPostProcEnergy);
 
      double body_energy;
      MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
                    eP.mField.get_comm());
      MOFEM_LOG_C("EP", Sev::inform, "Step %d time %3.4g strain energy %3.6e",
                  ts_step, ts_t, body_energy);
 
      auto post_proc_at_points = [&](std::array<double, 3> point,
                                     std::string str) {
        MoFEMFunctionBegin;
 
        dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
 
        struct OpPrint : public VolOp {
 
          EshelbianCore &eP;
          std::array<double, 3> point;
          std::string str;
 
          OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
                  std::string &str)
              : VolOp(ep.spatialDisp, VolOp::OPROW), eP(ep), point(point),
                str(str) {}
 
          MoFEMErrorCode doWork(int side, EntityType type,
                                DataForcesAndSourcesCore::EntData &data) {
            MoFEMFunctionBegin;
            if (type == MBTET) {
              if (getGaussPts().size2()) {
 
                auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
                auto t_approx_P =
                    getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
 
                FTensor::Index<'i', 3> i;
                FTensor::Index<'j', 3> j;
                FTensor::Index<'k', 3> k;
                const double jac = determinantTensor3by3(t_h);
                FTensor::Tensor2<double, 3, 3> t_cauchy;
                t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
 
                auto add = [&]() {
                  std::ostringstream s;
                  s << str << " elem " << getFEEntityHandle() << " ";
                  return s.str();
                };
 
                auto print_tensor = [](auto &t) {
                  std::ostringstream s;
                  s << t;
                  return s.str();
                };
 
                std::ostringstream print;
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "comm rank " << eP.mField.get_comm_rank();
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "point " << getVectorAdaptor(point.data(), 3);
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "coords at gauss pts " << getCoordsAtGaussPts();
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "w " << eP.dataAtPts->wAtPts;
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "Piola " << eP.dataAtPts->approxPAtPts;
                MOFEM_LOG("EPSYNC", Sev::inform)
                    << add() << "Cauchy " << print_tensor(t_cauchy);
              }
            }
            MoFEMFunctionReturn(0);
          }
        };
 
        if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
 
          fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
          CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
              ->evalFEAtThePoint3D(
                  point.data(), 1e-12, problemPtr->getName(), "EP",
                  dataFieldEval, eP.mField.get_comm_rank(),
                  eP.mField.get_comm_rank(), nullptr, MF_EXIST, QUIET);
          fe_ptr->getOpPtrVector().pop_back();
        }
 
        MoFEMFunctionReturn(0);
      };
 
      // Points for Cook beam
      std::array<double, 3> pointA = {48., 60., 5.};
      CHKERR post_proc_at_points(pointA, "Point A");
      MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
 
      std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
      CHKERR post_proc_at_points(pointB, "Point B");
      MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
 
      std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
      CHKERR post_proc_at_points(pointC, "Point C");
      MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
 
      MoFEMFunctionReturn(0);
    }
  };
 
  boost::shared_ptr<FEMethod> monitor_ptr(new Monitor(*this));
  ts_ctx->getLoopsMonitor().push_back(
      TsCtx::PairNameFEMethodPtr(elementVolumeName, monitor_ptr));
 
  CHKERR TSAppendOptionsPrefix(ts, "elastic_");
  CHKERR TSSetFromOptions(ts);
  MoFEMFunctionReturn(0);
}

solveElastic()

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

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 1266 of file EshelbianPlasticity.cpp.

                                                       {
  MoFEMFunctionBegin;
 
  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 DMGetDefaultSection(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
 
  KSP ksp;
  CHKERR SNESGetKSP(snes, &ksp);
  PC pc;
  CHKERR KSPGetPC(ksp, &pc);
  PetscBool is_uu_field_split;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_uu_field_split);
  if (is_uu_field_split) {
 
    const MoFEM::Problem *schur_uu_ptr;
    CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_uu_ptr);
    if (auto uu_data = schur_uu_ptr->subProblemData) {
 
      const MoFEM::Problem *prb_ptr;
      CHKERR DMMoFEMGetProblemPtr(dmElastic, &prb_ptr);
      map<std::string, IS> is_map;
      for (int ff = 0; ff != num_fields; ff++) {
        const char *field_name;
        CHKERR PetscSectionGetFieldName(section, ff, &field_name);
        CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
            prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
            &is_map[field_name]);
      }
      // CHKERR mField.getInterface<ISManager>()
      //     ->isCreateProblemFieldAndEntityType(
      //         prb_ptr->getName(), ROW, piolaStress, MBTET, MBTET, 0,
      //         MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TETS"]);
      // CHKERR mField.getInterface<ISManager>()
      //     ->isCreateProblemFieldAndEntityType(
      //         prb_ptr->getName(), ROW, piolaStress, MBTRI, MBTRI, 0,
      //         MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TRIS"]);
 
      CHKERR uu_data->getRowIs(&is_map["E_IS_SUU"]);
 
      CHKERR PCFieldSplitSetIS(pc, NULL, is_map[streachTensor]);
      CHKERR PCFieldSplitSetIS(pc, NULL, is_map["E_IS_SUU"]);
 
      CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
                                     schurAssembly->Suu);
 
      CHKERR PCSetUp(pc);
      PetscInt n;
      KSP *uu_ksp;
      CHKERR PCFieldSplitGetSubKSP(pc, &n, &uu_ksp);
      PC bubble_pc;
      CHKERR KSPGetPC(uu_ksp[1], &bubble_pc);
      PetscBool is_bubble_field_split;
      PetscObjectTypeCompare((PetscObject)bubble_pc, PCFIELDSPLIT,
                             &is_bubble_field_split);
      if (is_bubble_field_split) {
 
        const MoFEM::Problem *schur_bubble_ptr;
        CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_ptr);
        if (auto bubble_data = schur_bubble_ptr->subProblemData) {
 
          CHKERR bubble_data->getRowIs(&is_map["E_IS_BUBBLE"]);
 
          AO uu_ao;
          CHKERR uu_data->getRowMap(&uu_ao);
 
          CHKERR AOApplicationToPetscIS(uu_ao, is_map[bubbleField]);
          CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map[bubbleField]);
          CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map["E_IS_BUBBLE"]);
          CHKERR PCFieldSplitSetSchurPre(
              bubble_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->SBubble);
 
          CHKERR PCSetUp(bubble_pc);
          PetscInt bubble_n;
          KSP *bubble_ksp;
          CHKERR PCFieldSplitGetSubKSP(bubble_pc, &bubble_n, &bubble_ksp);
          PC omega_pc;
          CHKERR KSPGetPC(bubble_ksp[1], &omega_pc);
          PetscBool is_omega_field_split;
          PetscObjectTypeCompare((PetscObject)omega_pc, PCFIELDSPLIT,
                                 &is_omega_field_split);
 
          if (is_omega_field_split) {
 
            const MoFEM::Problem *schur_omega_ptr;
            CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_ptr);
            if (auto omega_data = schur_omega_ptr->subProblemData) {
 
              AO bubble_ao;
              CHKERR bubble_data->getRowMap(&bubble_ao);
 
              CHKERR AOApplicationToPetscIS(uu_ao, is_map[rotAxis]);
              CHKERR AOApplicationToPetscIS(bubble_ao, is_map[rotAxis]);
              CHKERR omega_data->getRowIs(&is_map["E_IS_OMEGA"]);
 
              CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map[rotAxis]);
              CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map["E_IS_OMEGA"]);
              CHKERR PCFieldSplitSetSchurPre(omega_pc,
                                             PC_FIELDSPLIT_SCHUR_PRE_USER,
                                             schurAssembly->SOmega);
 
              CHKERR PCSetUp(omega_pc);
              PetscInt omega_n;
              KSP *omega_ksp;
              CHKERR PCFieldSplitGetSubKSP(omega_pc, &omega_n, &omega_ksp);
              PC w_pc;
              CHKERR KSPGetPC(omega_ksp[1], &w_pc);
              PetscBool is_w_field_split;
              PetscObjectTypeCompare((PetscObject)w_pc, PCFIELDSPLIT,
                                     &is_w_field_split);
              if (is_w_field_split) {
 
                const MoFEM::Problem *schur_w_ptr;
                CHKERR DMMoFEMGetProblemPtr(dmElasticSchurSpatialDisp,
                                            &schur_w_ptr);
                if (auto w_data = schur_w_ptr->subProblemData) {
 
                  AO omega_ao;
                  CHKERR omega_data->getRowMap(&omega_ao);
 
                  CHKERR AOApplicationToPetscIS(uu_ao, is_map[spatialDisp]);
                  CHKERR AOApplicationToPetscIS(bubble_ao, is_map[spatialDisp]);
                  CHKERR AOApplicationToPetscIS(omega_ao, is_map[spatialDisp]);
                  CHKERR w_data->getRowIs(&is_map["E_IS_W"]);
 
                  CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map[spatialDisp]);
                  CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map["E_IS_W"]);
                  CHKERR PCFieldSplitSetSchurPre(
                      w_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->Sw);
 
                  CHKERR AODestroy(&omega_ao);
                }
              }
 
              CHKERR PetscFree(omega_ksp);
            }
          }
          CHKERR PetscFree(bubble_ksp);
          CHKERR AODestroy(&uu_ao);
        }
      }
      CHKERR PetscFree(uu_ksp);
 
      for (auto &m : is_map)
        CHKERR ISDestroy(&m.second);
    }
  }
 
  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 TSSolve(ts, PETSC_NULL);
 
  // 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)
    if (lin_solver_iterations != 2)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Expected number of iterations is different than expected");
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

alphaRho

double EshelbianPlasticity::EshelbianCore::alphaRho

Examples

EshelbianPlasticity.cpp.

Definition at line 1210 of file EshelbianPlasticity.hpp.

alphaU

double EshelbianPlasticity::EshelbianCore::alphaU

Examples

EshelbianPlasticity.cpp.

Definition at line 1208 of file EshelbianPlasticity.hpp.

alphaW

double EshelbianPlasticity::EshelbianCore::alphaW

Examples

EshelbianPlasticity.cpp.

Definition at line 1209 of file EshelbianPlasticity.hpp.

bcSpatialDispVecPtr

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1215 of file EshelbianPlasticity.hpp.

bcSpatialFreeTraction

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1218 of file EshelbianPlasticity.hpp.

bcSpatialRotationVecPtr

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1216 of file EshelbianPlasticity.hpp.

bcSpatialTraction

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1217 of file EshelbianPlasticity.hpp.

bubbleField

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1197 of file EshelbianPlasticity.hpp.

dataAtPts

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1170 of file EshelbianPlasticity.hpp.

dM

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dM

Coupled problem all fields.

Examples

EshelbianPlasticity.cpp.

Definition at line 1179 of file EshelbianPlasticity.hpp.

dmElastic

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElastic

Elastic problem.

Examples

EshelbianPlasticity.cpp, and ep.cpp.

Definition at line 1180 of file EshelbianPlasticity.hpp.

dmElasticSchurBubble

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElasticSchurBubble

Sub problem of dmElastic Schur.

Examples

EshelbianPlasticity.cpp.

Definition at line 1182 of file EshelbianPlasticity.hpp.

dmElasticSchurOmega

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElasticSchurOmega

Sub problem of dmElastic Schur.

Examples

EshelbianPlasticity.cpp.

Definition at line 1185 of file EshelbianPlasticity.hpp.

dmElasticSchurSpatialDisp

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElasticSchurSpatialDisp

Sub problem of dmElastic Schur.

Examples

EshelbianPlasticity.cpp.

Definition at line 1184 of file EshelbianPlasticity.hpp.

dmElasticSchurStreach

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmElasticSchurStreach

Sub problem of dmElastic Schur.

Examples

EshelbianPlasticity.cpp.

Definition at line 1181 of file EshelbianPlasticity.hpp.

dmMaterial

SmartPetscObj<DM> EshelbianPlasticity::EshelbianCore::dmMaterial

Material problem.

Definition at line 1186 of file EshelbianPlasticity.hpp.

elasticBcLhs

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1175 of file EshelbianPlasticity.hpp.

elasticBcRhs

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1176 of file EshelbianPlasticity.hpp.

elasticFeLhs

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1174 of file EshelbianPlasticity.hpp.

elasticFeRhs

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1173 of file EshelbianPlasticity.hpp.

elementVolumeName

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1199 of file EshelbianPlasticity.hpp.

eshelbyStress

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1189 of file EshelbianPlasticity.hpp.

essentialBcElement

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1201 of file EshelbianPlasticity.hpp.

lambdaField

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

Definition at line 1196 of file EshelbianPlasticity.hpp.

materialDisp

const std::string EshelbianPlasticity::EshelbianCore::materialDisp

Definition at line 1191 of file EshelbianPlasticity.hpp.

materialGradient

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1194 of file EshelbianPlasticity.hpp.

materialOrder

int EshelbianPlasticity::EshelbianCore::materialOrder

Examples

EshelbianPlasticity.cpp.

Definition at line 1207 of file EshelbianPlasticity.hpp.

mField

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1168 of file EshelbianPlasticity.hpp.

naturalBcElement

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1200 of file EshelbianPlasticity.hpp.

physicalEquations

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1171 of file EshelbianPlasticity.hpp.

piolaStress

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1188 of file EshelbianPlasticity.hpp.

precEps

double EshelbianPlasticity::EshelbianCore::precEps

Examples

EshelbianPlasticity.cpp.

Definition at line 1211 of file EshelbianPlasticity.hpp.

rotAxis

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

Examples

EshelbianPlasticity.cpp.

Definition at line 1193 of file EshelbianPlasticity.hpp.

schurAssembly

boost::shared_ptr<EpFEMethod> EshelbianPlasticity::EshelbianCore::schurAssembly

Examples

EshelbianPlasticity.cpp.

Definition at line 1177 of file EshelbianPlasticity.hpp.

spaceOrder

int EshelbianPlasticity::EshelbianCore::spaceOrder

Examples

EshelbianPlasticity.cpp.

Definition at line 1206 of file EshelbianPlasticity.hpp.

spatialDisp

const std::string EshelbianPlasticity::EshelbianCore::spatialDisp

Examples

EshelbianPlasticity.cpp.

Definition at line 1190 of file EshelbianPlasticity.hpp.

streachTensor

const std::string EshelbianPlasticity::EshelbianCore::streachTensor

Examples

EshelbianPlasticity.cpp.

Definition at line 1192 of file EshelbianPlasticity.hpp.

tauField

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

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