v0.13.1
Loading...
Searching...
No Matches
contact.cpp

Example of contact problem

/**
* \file contact.cpp
* \example contact.cpp
*
* Example of contact problem
*
*/
#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
#include <MoFEM.hpp>
using namespace MoFEM;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
static constexpr FieldSpace CONTACT_SPACE = HCURL;
};
template <> struct ElementsAndOps<3> {
static constexpr FieldSpace CONTACT_SPACE = HDIV;
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
constexpr EntityType boundary_ent = SPACE_DIM == 3 ? MBTRI : MBEDGE;
//! [Operators used for contact]
GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
//! [Operators used for contact]
//! [Body force]
GAUSS>::OpSource<1, SPACE_DIM>;
//! [Body force]
//! [Only used with Hooke equation (linear material model)]
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
//! [Only used with Hooke equation (linear material model)]
//! [Only used for dynamics]
//! [Only used for dynamics]
//! [Essential boundary conditions]
//! [Essential boundary conditions]
// Only used with Hencky/nonlinear material
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
constexpr bool is_quasi_static = true;
constexpr bool is_large_strains = true;
constexpr int order = 2;
constexpr double young_modulus = 100;
constexpr double poisson_ratio = 0.25;
constexpr double rho = 0;
constexpr double cn = 0.01;
constexpr double spring_stiffness = 0.1;
#include <ContactOps.hpp>
#include <HenckyOps.hpp>
using namespace ContactOps;
using namespace HenckyOps;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
private:
boost::shared_ptr<ContactOps::CommonData> commonDataPtr;
boost::shared_ptr<PostProcEle> postProcFe;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
template <int DIM> Range getEntsOnMeshSkin();
};
//! [Run problem]
}
//! [Run problem]
//! [Set up problem]
// Select base
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
PetscInt choice_base_value = AINSWORTH;
CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
LASBASETOPT, &choice_base_value, PETSC_NULL);
switch (choice_base_value) {
case AINSWORTH:
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Set AINSWORTH_LEGENDRE_BASE for displacements";
break;
case DEMKOWICZ:
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Set DEMKOWICZ_JACOBI_BASE for displacents";
break;
default:
break;
}
// Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
// Demkowicz base. We need to implement Demkowicz H1 base on tet.
CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setFieldOrder("SIGMA", 0);
auto skin_edges = getEntsOnMeshSkin<SPACE_DIM>();
// filter not owned entities, those are not on boundary
Range boundary_ents;
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
if (pcomm == NULL) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Communicator not created");
}
CHKERR pcomm->filter_pstatus(skin_edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
// Range adj_edges;
// CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false,
// adj_edges,
// moab::Interface::UNION);
// adj_edges.merge(boundary_ents);
// CHKERR simple->setFieldOrder("U", order, &adj_edges);
CHKERR simple->setUp();
}
//! [Set up problem]
//! [Create common data]
auto set_matrial_stiffness = [&]() {
constexpr double bulk_modulus_K =
young_modulus / (3 * (1 - 2 * poisson_ratio));
constexpr double shear_modulus_G =
constexpr double A =
(SPACE_DIM == 2) ? 2 * shear_modulus_G /
: 1;
auto t_D =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*commonDataPtr->mDPtr);
t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
t_kd(i, j) * t_kd(k, l);
};
commonDataPtr = boost::make_shared<ContactOps::CommonData>();
commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->contactStressDivergencePtr =
boost::make_shared<MatrixDouble>();
commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->contactDispPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->curlContactStressPtr = boost::make_shared<MatrixDouble>();
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
commonDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->mDPtr->resize(size_symm * size_symm, 1);
CHKERR set_matrial_stiffness();
}
//! [Create common data]
//! [Boundary condition]
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"NO_CONTACT", "SIGMA", 0, 3);
CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
}
//! [Boundary condition]
//! [Push operators to pipeline]
auto bc_mng = mField.getInterface<BcManager>();
auto add_domain_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
pipeline.push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(
if (SPACE_DIM == 2) {
pipeline.push_back(new OpMakeHdivFromHcurl());
pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
pipeline.push_back(new OpSetHOWeightsOnFace());
} else {
pipeline.push_back(
new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
pipeline.push_back(new OpSetHOWeights(det_ptr));
}
};
auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
henky_common_data_ptr->matGradPtr = commonDataPtr->mGradPtr;
henky_common_data_ptr->matDPtr = commonDataPtr->mDPtr;
auto add_domain_ops_lhs = [&](auto &pipeline) {
pipeline_mng->getOpDomainLhsPipeline().push_back(
"U", commonDataPtr->mGradPtr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpHenckyTangent<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpKPiola("U", "U", henky_common_data_ptr->getMatTangent()));
} else {
pipeline.push_back(new OpKCauchy("U", "U", commonDataPtr->mDPtr));
}
// Get pointer to U_tt shift in domain element
auto get_rho = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
return rho * fe_domain_lhs->ts_aa;
};
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMass("U", "U", get_rho));
}
auto unity = []() { return 1; };
pipeline.push_back(new OpMixDivULhs("SIGMA", "U", unity, true));
pipeline.push_back(new OpLambdaGraULhs("SIGMA", "U", unity, true));
};
auto add_domain_ops_rhs = [&](auto &pipeline) {
auto get_body_force = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
const auto time = fe_domain_rhs->ts_t;
// hardcoded gravity load
t_source(i) = 0;
t_source(1) = 1.0 * time;
return t_source;
};
pipeline.push_back(new OpBodyForce("U", get_body_force));
"U", commonDataPtr->mGradPtr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpInternalForcePiola(
"U", henky_common_data_ptr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
"U", commonDataPtr->mGradPtr, commonDataPtr->mStrainPtr));
"U", commonDataPtr->mStrainPtr, commonDataPtr->mStressPtr,
commonDataPtr->mDPtr));
pipeline.push_back(
new OpInternalForceCauchy("U", commonDataPtr->mStressPtr));
}
"U", commonDataPtr->contactDispPtr));
"SIGMA", commonDataPtr->contactStressPtr));
pipeline.push_back(
"SIGMA", commonDataPtr->contactStressDivergencePtr));
pipeline.push_back(
new OpMixDivURhs("SIGMA", commonDataPtr->contactDispPtr,
[](double, double, double) { return 1; }));
pipeline.push_back(
new OpMixLambdaGradURhs("SIGMA", commonDataPtr->mGradPtr));
pipeline.push_back(new OpMixUTimesDivLambdaRhs(
"U", commonDataPtr->contactStressDivergencePtr));
pipeline.push_back(
new OpMixUTimesLambdaRhs("U", commonDataPtr->contactStressPtr));
// only in case of dynamics
auto mat_acceleration = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
mat_acceleration));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpInertiaForce(
"U", mat_acceleration, [](double, double, double) { return rho; }));
}
};
auto add_boundary_base_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetPiolaTransformOnBoundary(CONTACT_SPACE));
if (SPACE_DIM == 3)
pipeline.push_back(new OpSetHOWeightsOnFace());
"U", commonDataPtr->contactDispPtr));
"SIGMA", commonDataPtr->contactTractionPtr));
};
auto add_boundary_ops_lhs = [&](auto &pipeline) {
auto &bc_map = bc_mng->getBcMapByBlockName();
for (auto bc : bc_map) {
if (bc_mng->checkBlock(bc, "FIX_")) {
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Set boundary matrix for " << bc.first;
pipeline.push_back(
new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
pipeline.push_back(new OpBoundaryMass(
"U", "U", [](double, double, double) { return 1.; },
bc.second->getBcEntsPtr()));
}
}
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(new OpSpringLhs(
"U", "U",
[this](double, double, double) { return spring_stiffness; }
));
};
auto time_scaled = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getDomainRhsFE();
return -fe_domain_rhs->ts_t;
};
auto add_boundary_ops_rhs = [&](auto &pipeline) {
if (bc_mng->checkBlock(bc, "FIX_")) {
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Set boundary residual for " << bc.first;
pipeline.push_back(
new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
auto attr_vec = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
attr_vec->clear();
if (bc.second->bcAttributes.size() != SPACE_DIM)
SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Wrong size of boundary attributes vector. Set right block "
"size attributes. Size of attributes %d",
bc.second->bcAttributes.size());
std::copy(&bc.second->bcAttributes[0],
&bc.second->bcAttributes[SPACE_DIM],
attr_vec->data().begin());
pipeline.push_back(new OpBoundaryVec("U", attr_vec, time_scaled,
bc.second->getBcEntsPtr()));
pipeline.push_back(new OpBoundaryInternal(
"U", commonDataPtr->contactDispPtr,
[](double, double, double) { return 1.; },
bc.second->getBcEntsPtr()));
}
}
pipeline.push_back(new OpConstrainBoundaryRhs("SIGMA", commonDataPtr));
pipeline.push_back(new OpSpringRhs(
"U", commonDataPtr->contactDispPtr,
[this](double, double, double) { return spring_stiffness; }));
};
add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
CHKERR add_boundary_ops_lhs(pipeline_mng->getOpBoundaryLhsPipeline());
CHKERR add_boundary_ops_rhs(pipeline_mng->getOpBoundaryRhsPipeline());
auto integration_rule_vol = [](int, int, int approx_order) {
return 3 * order;
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_vol);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_vol);
auto integration_rule_boundary = [](int, int, int approx_order) {
return 3 * order;
};
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
// Enforce non-homegonus boundary conditions
auto get_bc_hook = [&]() {
mField, pipeline_mng->getDomainRhsFE(),
{boost::make_shared<TimeScale>()});
return hook;
};
pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook();
pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook();
}
//! [Push operators to pipeline]
//! [Solve]
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &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);
};
auto scatter_create = [&](auto D, auto coeff) {
CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
ROW, "U", coeff, coeff, is);
int loc_size;
CHKERR ISGetLocalSize(is, &loc_size);
Vec v;
CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
VecScatter scatter;
CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
return std::make_tuple(SmartPetscObj<Vec>(v),
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
boost::shared_ptr<ForcesAndSourcesCore> null;
CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
monitor_ptr, null, null);
};
auto dm = simple->getDM();
auto D = smartCreateDMVector(dm);
uXScatter = scatter_create(D, 0);
uYScatter = scatter_create(D, 1);
if (SPACE_DIM == 3)
uZScatter = scatter_create(D, 2);
auto solver = pipeline_mng->createTSIM();
auto D = smartCreateDMVector(dm);
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetSolution(solver, D);
CHKERR TSSetFromOptions(solver);
CHKERR TSSetUp(solver);
CHKERR TSSolve(solver, NULL);
} else {
auto solver = pipeline_mng->createTSIM2();
auto dm = simple->getDM();
auto D = smartCreateDMVector(dm);
auto DD = smartVectorDuplicate(D);
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TS2SetSolution(solver, D, DD);
CHKERR TSSetFromOptions(solver);
CHKERR TSSetUp(solver);
CHKERR TSSolve(solver, NULL);
}
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Solve]
//! [Postprocess results]
//! [Postprocess results]
//! [Check]
//! [Check]
template <int DIM> Range Example::getEntsOnMeshSkin() {
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, DIM, body_ents);
Skinner skin(&mField.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, body_ents, false, skin_ents);
return skin_ents;
};
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Load mesh]
CHKERR simple->getOptions();
CHKERR simple->loadFile("");
//! [Load mesh]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
std::string param_file
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static char help[]
int main()
Definition: adol-c_atom.cpp:46
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
constexpr double spring_stiffness
Definition: contact.cpp:125
ElementsAndOps< SPACE_DIM >::OpSetPiolaTransformOnBoundary OpSetPiolaTransformOnBoundary
Definition: contact.cpp:61
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixTensorTimesGradU< SPACE_DIM > OpMixLambdaGradURhs
Definition: contact.cpp:72
constexpr bool is_quasi_static
Definition: contact.cpp:117
constexpr FieldSpace CONTACT_SPACE
Definition: contact.cpp:62
constexpr EntityType boundary_ent
Definition: contact.cpp:48
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
FieldSpace
approximation spaces
Definition: definitions.h:82
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr double shear_modulus_G
Definition: elastic.cpp:63
constexpr double bulk_modulus_K
Definition: elastic.cpp:62
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
constexpr auto t_kd
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'k', 3 > k
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMMoFEM.cpp:1003
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
double young_modulus
Definition: plastic.cpp:107
double rho
Definition: plastic.cpp:109
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:10
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:77
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
Definition: plastic.cpp:61
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:82
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:63
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:84
constexpr auto size_symm
Definition: plastic.cpp:44
constexpr AssemblyType A
Definition: plastic.cpp:46
PetscBool is_large_strains
Definition: plastic.cpp:103
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:75
double cn
Definition: plastic.cpp:113
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:86
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
[Example]
Definition: plastic.cpp:139
Range getEntsOnMeshSkin()
[Check]
Definition: contact.cpp:662
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: contact.cpp:150
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:747
FieldApproximationBase base
Definition: plot_base.cpp:68
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:159
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: contact.cpp:659
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:259
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:455
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:174
MoFEM::Interface & mField
Definition: plastic.cpp:146
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:160
MoFEMErrorCode postProcess()
[Solve]
Definition: contact.cpp:655
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:186
boost::shared_ptr< PostProcEle > postProcFe
Definition: contact.cpp:151
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:393
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:158
Simple interface for fast problem set-up.
Definition: BcManager.hpp:23
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:165
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:72
Data on single entity (This is passed as argument to DataOperator::doWork)
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:39
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Approximate field valuse for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
transform Hdiv base fluxes from reference element to physical triangle
Make Hdiv space from Hcurl space in 2d.
Set indices on entities on finite element.
Apply contravariant (Piola) transfer to Hdiv space for HO geometry.
Set inverse jacobian to base functions.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.
Postprocess on face.
Post processing.