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

Example of contact problem

/**
* \file contact.cpp
* \example contact.cpp
*
* Example of contact problem
*
* @copyright Copyright (c) 2023
*/
#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
#include <MoFEM.hpp>
using namespace MoFEM;
constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType I =
IntegrationType::GAUSS; //< selected integration type
template <int DIM> struct ElementsAndOps {};
static constexpr FieldSpace CONTACT_SPACE = HCURL;
};
static constexpr FieldSpace CONTACT_SPACE = HDIV;
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
//! [Operators used for contact]
GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
//! [Operators used for contact]
//! [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>;
namespace ContactOps {
struct DomainBCs {};
struct BoundaryBCs {};
using OpDomainRhsBCs = DomainRhsBCs::OpFlux<DomainBCs, 1, SPACE_DIM>;
using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
using OpBoundaryLhsBCs = BoundaryLhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
MoFEM::Interface &m_field,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string field_name, std::string block_name,
boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev);
}; // namespace ContactOps
constexpr bool is_quasi_static = 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 HenckyOps;
using namespace ContactOps;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
private:
boost::shared_ptr<ContactOps::CommonData> commonDataPtr;
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 displacements";
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]
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>();
commonDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
}
//! [Create common data]
//! [Boundary condition]
auto bc_mng = mField.getInterface<BcManager>();
for (auto f : {"U", "SIGMA"}) {
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_X", f, 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_Y", f, 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_Z", f, 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_ALL", f, 0, 3);
}
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"NO_CONTACT", "SIGMA", 0, 3);
CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
}
//! [Boundary condition]
//! [Push operators to pip]
auto bc_mng = mField.getInterface<BcManager>();
auto time_scale = boost::make_shared<TimeScale>();
auto add_domain_base_ops = [&](auto &pip) {
};
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 &pip) {
CHKERR addMatBlockOps(mField, pip, "U", "MAT_ELASTIC", commonDataPtr->mDPtr,
Sev::verbose);
"U", commonDataPtr->mGradPtr));
pip.push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(
new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(new OpHenckyTangent<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(
new OpKPiola("U", "U", henky_common_data_ptr->getMatTangent()));
// Get pointer to U_tt shift in domain element
auto get_rho = [this](const double, const double, const double) {
auto *pip_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_lhs = pip_mng->getDomainLhsFE();
return rho * fe_domain_lhs->ts_aa;
};
pip.push_back(new OpMass("U", "U", get_rho));
}
auto unity = []() { return 1; };
pip.push_back(new OpMixDivULhs("SIGMA", "U", unity, true));
pip.push_back(new OpLambdaGraULhs("SIGMA", "U", unity, true));
};
auto add_domain_ops_rhs = [&](auto &pip) {
CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
pip, mField, "U", {time_scale}, Sev::inform);
CHKERR addMatBlockOps(mField, pip, "U", "MAT_ELASTIC", commonDataPtr->mDPtr,
Sev::inform);
"U", commonDataPtr->mGradPtr));
pip.push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(
new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
pip.push_back(new OpInternalForcePiola(
"U", henky_common_data_ptr->getMatFirstPiolaStress()));
"U", commonDataPtr->contactDispPtr));
"SIGMA", commonDataPtr->contactStressPtr));
"SIGMA", commonDataPtr->contactStressDivergencePtr));
pip.push_back(new OpMixDivURhs("SIGMA", commonDataPtr->contactDispPtr,
[](double, double, double) { return 1; }));
pip.push_back(new OpMixLambdaGradURhs("SIGMA", commonDataPtr->mGradPtr));
pip.push_back(new OpMixUTimesDivLambdaRhs(
"U", commonDataPtr->contactStressDivergencePtr));
pip.push_back(
new OpMixUTimesLambdaRhs("U", commonDataPtr->contactStressPtr));
// only in case of dynamics
auto mat_acceleration = boost::make_shared<MatrixDouble>();
"U", mat_acceleration));
pip.push_back(new OpInertiaForce(
"U", mat_acceleration, [](double, double, double) { return rho; }));
}
};
auto add_boundary_base_ops = [&](auto &pip) {
"U", commonDataPtr->contactDispPtr));
"SIGMA", commonDataPtr->contactTractionPtr));
};
auto add_boundary_ops_lhs = [&](auto &pip) {
CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
pip, mField, "U", Sev::inform);
pip.push_back(new OpConstrainBoundaryLhs_dU("SIGMA", "U", commonDataPtr));
pip.push_back(
pip.push_back(new OpSpringLhs(
"U", "U",
[this](double, double, double) { return spring_stiffness; }
));
};
auto add_boundary_ops_rhs = [&](auto &pip) {
CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
pip, mField, "U", {time_scale}, Sev::inform);
pip.push_back(new OpConstrainBoundaryRhs("SIGMA", commonDataPtr));
pip.push_back(new OpSpringRhs(
"U", commonDataPtr->contactDispPtr,
[this](double, double, double) { return spring_stiffness; }));
};
add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
auto integration_rule_vol = [](int, int, int approx_order) {
return 3 * order;
};
CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
auto integration_rule_boundary = [](int, int, int approx_order) {
return 3 * order;
};
CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
// Enforce non-homogenous boundary conditions
auto get_bc_hook_rhs = [&]() {
mField, pip_mng->getDomainRhsFE(), {boost::make_shared<TimeScale>()});
return hook;
};
auto get_bc_hook_lhs = [&]() {
mField, pip_mng->getDomainLhsFE(), {boost::make_shared<TimeScale>()});
return hook;
};
pip_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
pip_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
}
//! [Push operators to pip]
//! [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 = pip_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 = pip_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 interface
//! [Create MoFEM]
//! [Load mesh]
CHKERR simple->getOptions();
CHKERR simple->loadFile("");
//! [Load mesh]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
MoFEM::Interface &m_field,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string field_name, std::string block_name,
boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev) {
struct OpMatBlocks : public DomainEleOp {
OpMatBlocks(std::string field_name, boost::shared_ptr<MatrixDouble> m,
MoFEM::Interface &m_field, Sev sev,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
: DomainEleOp(field_name, DomainEleOp::OPROW), matDPtr(m),
bulkModulusKDefault(bulk_modulus_K),
shearModulusGDefault(shear_modulus_G) {
std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
"Can not get data from block");
}
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
}
}
CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
}
private:
boost::shared_ptr<MatrixDouble> matDPtr;
struct BlockData {
double bulkModulusK;
double shearModulusG;
Range blockEnts;
};
double bulkModulusKDefault;
double shearModulusGDefault;
std::vector<BlockData> blockData;
extractBlockData(MoFEM::Interface &m_field,
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (auto m : meshset_vec_ptr) {
MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() != 2) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Expected that block has two attribute");
}
auto get_block_ents = [&]() {
Range ents;
m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
return ents;
};
double young_modulus = block_data[0];
double poisson_ratio = block_data[1];
double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
<< "E = " << young_modulus << " nu = " << poisson_ratio;
blockData.push_back(
{bulk_modulus_K, shear_modulus_G, get_block_ents()});
}
}
MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
double bulk_modulus_K, double shear_modulus_G) {
//! [Calculate elasticity tensor]
auto set_material_stiffness = [&]() {
double A = (SPACE_DIM == 2)
: 1;
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
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);
};
//! [Calculate elasticity tensor]
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
mat_D_ptr->resize(size_symm * size_symm, 1);
set_material_stiffness();
}
};
double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
pipeline.push_back(new OpMatBlocks(
field_name, mat_D_Ptr, bulk_modulus_K, shear_modulus_G, m_field, sev,
// Get blockset using regular expression
(boost::format("%s(.*)") % block_name).str()
))
));
}
std::string param_file
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:345
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:119
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< I >::OpMixTensorTimesGradU< SPACE_DIM > OpMixLambdaGradURhs
Definition: contact.cpp:59
constexpr double cn
Definition: contact.cpp:118
constexpr bool is_quasi_static
Definition: contact.cpp:112
constexpr FieldSpace CONTACT_SPACE
Definition: contact.cpp:49
@ 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
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
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
NaturalBC< BoundaryEleOp >::Assembly< A >::LinearForm< I > BoundaryRhsBCs
Definition: elastic.cpp:43
NaturalBC< DomainEleOp >::Assembly< A >::LinearForm< I > DomainRhsBCs
Definition: elastic.cpp:41
BoundaryLhsBCs::OpFlux< BoundaryBCs, 1, SPACE_DIM > OpBoundaryLhsBCs
Definition: elastic.cpp:46
NaturalBC< BoundaryEleOp >::Assembly< A >::BiLinearForm< I > BoundaryLhsBCs
Definition: elastic.cpp:45
constexpr double shear_modulus_G
Definition: elastic.cpp:60
constexpr double bulk_modulus_K
Definition: elastic.cpp:59
DomainRhsBCs::OpFlux< DomainBCs, 1, SPACE_DIM > OpDomainRhsBCs
Definition: elastic.cpp:42
BoundaryRhsBCs::OpFlux< BoundaryBCs, 1, SPACE_DIM > OpBoundaryRhsBCs
Definition: elastic.cpp:44
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
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: DMMoFEM.cpp:511
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:995
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:277
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
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 >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
Definition: contact.cpp:609
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: DMMoFEM.cpp:1044
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)
constexpr IntegrationType I
constexpr AssemblyType A
double young_modulus
Definition: plastic.cpp:99
double rho
Definition: plastic.cpp:101
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:10
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:68
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:73
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:75
constexpr auto size_symm
Definition: plastic.cpp:33
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:66
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:77
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
[Example]
Definition: plastic.cpp:141
Range getEntsOnMeshSkin()
[Check]
Definition: contact.cpp:551
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: contact.cpp:147
MoFEMErrorCode tsSolve()
[Push operators to pip]
Definition: plastic.cpp:752
FieldApproximationBase base
Definition: plot_base.cpp:68
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:161
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: contact.cpp:548
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:261
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:459
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:176
MoFEM::Interface & mField
Definition: plastic.cpp:148
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:162
MoFEMErrorCode postProcess()
[Solve]
Definition: contact.cpp:544
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:188
MoFEMErrorCode addMatBlockOps(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
Definition: elastic.cpp:86
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:397
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:160
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
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
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
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
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition: Natural.hpp:57
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 values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
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:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.