v0.13.1
contact.cpp

Example of contact problem

/**
* \file contact.cpp
* \example contact.cpp
*
* Example of contact problem
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#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;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
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(simple->getProblemName(), "REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_X", "U",
0, 0);
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Y", "U",
1, 1);
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Z", "U",
2, 2);
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
"U", 0, 3);
boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
}
//! [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(
pipeline.push_back(new OpSetHOWeights(det_ptr));
if (SPACE_DIM == 2) {
pipeline.push_back(new OpMakeHdivFromHcurl());
pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
} else {
pipeline.push_back(
new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac_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.push_back(new OpSetBc("U", true, boundaryMarker));
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));
pipeline.push_back(new OpUnSetBc("U"));
};
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 OpSetBc("U", true, boundaryMarker));
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; }));
}
pipeline.push_back(new OpUnSetBc("U"));
};
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(new OpUnSetBc("U"));
}
}
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
pipeline.push_back(
pipeline.push_back(new OpSpringLhs(
"U", "U",
[this](double, double, double) { return spring_stiffness; }
));
pipeline.push_back(new OpUnSetBc("U"));
};
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 OpUnSetBc("U"));
}
}
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(new OpConstrainBoundaryRhs("SIGMA", commonDataPtr));
pipeline.push_back(new OpSpringRhs(
"U", commonDataPtr->contactDispPtr,
[this](double, double, double) { return spring_stiffness; }));
pipeline.push_back(new OpUnSetBc("U"));
};
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);
}
//! [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
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
constexpr double spring_stiffness
Definition: contact.cpp:139
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< SPACE_DIM > OpMixUTimesDivLambdaRhs
Definition: contact.cpp:88
int main(int argc, char *argv[])
Definition: contact.cpp:708
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: contact.cpp:129
EntitiesFieldData::EntData EntData
Definition: contact.cpp:63
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Essential boundary conditions]
Definition: contact.cpp:127
constexpr double rho
Definition: contact.cpp:137
static char help[]
Definition: contact.cpp:706
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: contact.cpp:111
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesVec< SPACE_DIM > OpMixDivULhs
[Operators used for contact]
Definition: contact.cpp:80
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: contact.cpp:106
#define EXECUTABLE_DIMENSION
Definition: contact.cpp:24
constexpr int SPACE_DIM
Definition: contact.cpp:59
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpMixUTimesLambdaRhs
Definition: contact.cpp:90
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:94
ElementsAndOps< SPACE_DIM >::OpSetPiolaTransformOnBoundary OpSetPiolaTransformOnBoundary
Definition: contact.cpp:75
constexpr double poisson_ratio
Definition: contact.cpp:136
constexpr bool is_large_strains
Definition: contact.cpp:132
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: contact.cpp:120
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< SPACE_DIM > OpLambdaGraULhs
Definition: contact.cpp:82
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, SPACE_DIM, SPACE_DIM > OpMixDivURhs
Definition: contact.cpp:84
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: contact.cpp:104
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: contact.cpp:122
constexpr double cn
Definition: contact.cpp:138
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixTensorTimesGradU< SPACE_DIM > OpMixLambdaGradURhs
Definition: contact.cpp:86
constexpr int order
Definition: contact.cpp:134
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used for dynamics]
Definition: contact.cpp:118
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpSpringLhs
Definition: contact.cpp:92
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
[Operators used for contact]
Definition: contact.cpp:99
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: contact.cpp:113
constexpr bool is_quasi_static
Definition: contact.cpp:131
constexpr double young_modulus
Definition: contact.cpp:135
constexpr FieldSpace CONTACT_SPACE
Definition: contact.cpp:76
constexpr EntityType boundary_ent
Definition: contact.cpp:62
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ LASTBASE
Definition: definitions.h:82
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
FieldSpace
approximation spaces
Definition: definitions.h:95
@ H1
continuous field
Definition: definitions.h:98
@ HCURL
field with continuous tangents
Definition: definitions.h:99
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
constexpr double shear_modulus_G
Definition: elastic.cpp:66
constexpr double bulk_modulus_K
Definition: elastic.cpp:65
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
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:482
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
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:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
FTensor::Index< 'i', SPACE_DIM > i
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 >::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:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
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:1015
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
double A
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:88
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:86
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:74
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:95
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: plastic.cpp:72
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:97
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:93
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
[Example]
Definition: plastic.cpp:130
Range getEntsOnMeshSkin()
[Check]
Definition: contact.cpp:696
boost::shared_ptr< PostProcEle > postProcFe
Definition: plastic.cpp:147
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: contact.cpp:165
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:788
FieldApproximationBase base
Definition: plot_base.cpp:83
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:153
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: contact.cpp:693
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:195
Example(MoFEM::Interface &m_field)
Definition: plastic.cpp:132
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:387
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:163
MoFEM::Interface & mField
Definition: plastic.cpp:137
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:154
MoFEMErrorCode postProcess()
[Solve]
Definition: contact.cpp:689
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:156
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:175
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:319
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:152
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:129
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:33
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 geometr.
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:35
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.