v0.13.2
Loading...
Searching...
No Matches
ADV-1: Contact problem
Note
Prerequisites of this tutorial include VEC-0: Linear elasticity


Note
Intended learning outcome:
  • general structure of a program developed using MoFEM
  • idea of Simple Interface in MoFEM and how to use it
  • idea of Domain element in MoFEM and how to use it
  • Use default Forms Integrals
  • how to push the developed UDOs to the Pipeline
  • utilisation of tools to convert outputs (MOAB) and visualise them (Paraview)
/**
* \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_rhs = [&]() {
mField, pipeline_mng->getDomainRhsFE(),
{boost::make_shared<TimeScale>()});
return hook;
};
auto get_bc_hook_lhs = [&]() {
mField, pipeline_mng->getDomainLhsFE(),
{boost::make_shared<TimeScale>()});
return hook;
};
pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
}
//! [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
constexpr double cn
Definition: contact.cpp:124
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:60
constexpr double bulk_modulus_K
Definition: elastic.cpp:59
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: 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.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
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: 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 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< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
Definition: plastic.cpp:52
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:73
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:54
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:75
constexpr auto size_symm
Definition: plastic.cpp:33
PetscBool is_large_strains
Definition: plastic.cpp:94
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
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:141
Range getEntsOnMeshSkin()
[Check]
Definition: contact.cpp:669
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: contact.cpp:150
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
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:666
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:662
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:188
boost::shared_ptr< PostProcEle > postProcFe
Definition: contact.cpp:151
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:397
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:160
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:200
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 values 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:27
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.
/**
* \file ContactOps.hpp
* \example ContactOps.hpp
*/
namespace ContactOps {
//! [Common data]
struct CommonData {
boost::shared_ptr<MatrixDouble> mDPtr;
boost::shared_ptr<MatrixDouble> mGradPtr;
boost::shared_ptr<MatrixDouble> mStrainPtr;
boost::shared_ptr<MatrixDouble> mStressPtr;
boost::shared_ptr<MatrixDouble> contactStressPtr;
boost::shared_ptr<MatrixDouble> contactStressDivergencePtr;
boost::shared_ptr<MatrixDouble> contactTractionPtr;
boost::shared_ptr<MatrixDouble> contactDispPtr;
boost::shared_ptr<MatrixDouble> curlContactStressPtr;
};
//! [Common data]
struct OpConstrainBoundaryRhs : public AssemblyBoundaryEleOp {
OpConstrainBoundaryRhs(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpConstrainBoundaryLhs_dU : public AssemblyBoundaryEleOp {
OpConstrainBoundaryLhs_dU(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
struct OpConstrainBoundaryLhs_dTraction : public AssemblyBoundaryEleOp {
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <typename T1, typename T2>
t_normal(i) = 0;
t_normal(1) = 1.;
return t_normal;
}
template <typename T>
inline double gap0(FTensor::Tensor1<T, 3> &t_coords,
return (-0.5 - t_coords(1)) * t_normal(1);
}
template <typename T>
inline double gap(FTensor::Tensor1<T, SPACE_DIM> &t_disp,
return t_disp(i) * t_normal(i);
}
template <typename T>
return -t_traction(i) * t_normal(i);
}
inline double sign(double x) {
if (x == 0)
return 0;
else if (x > 0)
return 1;
else
return -1;
};
inline double w(const double g, const double t) { return g - cn * t; }
inline double constrian(double &&g0, double &&g, double &&t) {
return (w(g - g0, t) + std::abs(w(g - g0, t))) / 2 + g0;
};
inline double diff_constrains_dtraction(double &&g0, double &&g, double &&t) {
return -cn * (1 + sign(w(g - g0, t))) / 2;
}
inline double diff_constrains_dgap(double &&g0, double &&g, double &&t) {
return (1 + sign(w(g - g0, t))) / 2;
}
const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {}
OpConstrainBoundaryRhs::iNtegrate(EntitiesFieldData::EntData &data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto t_normal = getFTensor1Normal();
t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
auto t_w = getFTensor0IntegrationWeight();
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
size_t nb_base_functions = data.getN().size2() / 3;
auto t_base = data.getFTensor1N<3>();
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
auto t_nf = getFTensor1FromPtr<SPACE_DIM>(&nf[0]);
const double alpha = t_w * getMeasure();
auto t_contact_normal = normal(t_coords, t_disp);
t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
t_rhs_constrains(i) =
t_contact_normal(i) *
constrian(gap0(t_coords, t_contact_normal),
gap(t_disp, t_contact_normal),
normal_traction(t_traction, t_contact_normal));
t_rhs_tangent_traction;
t_rhs_tangent_disp(i) = t_Q(i, j) * t_disp(j);
t_rhs_tangent_traction(i) = cn * t_Q(i, j) * t_traction(j);
size_t bb = 0;
for (; bb != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++bb) {
const double beta = alpha * (t_base(i) * t_normal(i));
t_nf(i) -= beta * t_rhs_constrains(i);
t_nf(i) -= beta * t_rhs_tangent_disp(i);
t_nf(i) += beta * t_rhs_tangent_traction(i);
++t_nf;
++t_base;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
}
}
OpConstrainBoundaryLhs_dU::OpConstrainBoundaryLhs_dU(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
: AssemblyBoundaryEleOp(row_field_name, col_field_name,
DomainEleOp::OPROWCOL),
commonDataPtr(common_data_ptr) {
sYmm = false;
}
MoFEMErrorCode OpConstrainBoundaryLhs_dU::iNtegrate(
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const size_t nb_gauss_pts = getGaussPts().size2();
auto t_normal = getFTensor1Normal();
t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
const double alpha = t_w * getMeasure();
auto t_contact_normal = normal(t_coords, t_disp);
t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
auto diff_constrain = diff_constrains_dgap(
gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
normal_traction(t_traction, t_contact_normal));
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
locMat, SPACE_DIM * rr);
const double row_base = t_row_base(i) * t_normal(i);
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / SPACE_DIM;
++cc) {
const double beta = alpha * row_base * t_col_base;
t_mat(i, j) -= (beta * diff_constrain) * t_P(i, j);
t_mat(i, j) -= beta * t_Q(i, j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nb_face_functions; ++rr)
++t_row_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
}
}
OpConstrainBoundaryLhs_dTraction::OpConstrainBoundaryLhs_dTraction(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
: AssemblyBoundaryEleOp(row_field_name, col_field_name,
DomainEleOp::OPROWCOL),
commonDataPtr(common_data_ptr) {
sYmm = false;
}
MoFEMErrorCode OpConstrainBoundaryLhs_dTraction::iNtegrate(
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const size_t nb_gauss_pts = getGaussPts().size2();
auto t_normal = getFTensor1Normal();
t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
const double alpha = t_w * getMeasure();
auto t_contact_normal = normal(t_coords, t_disp);
t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
const double diff_traction = diff_constrains_dtraction(
gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
normal_traction(t_traction, t_contact_normal));
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
locMat, SPACE_DIM * rr);
const double row_base = t_row_base(i) * t_normal(i);
auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / SPACE_DIM;
++cc) {
const double col_base = t_col_base(i) * t_normal(i);
const double beta = alpha * row_base * col_base;
t_mat(i, j) += (beta * diff_traction) * t_P(i, j);
t_mat(i, j) += beta * cn * t_Q(i, j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nb_face_functions; ++rr)
++t_row_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
}
}
}; // namespace ContactOps
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
double normal_traction(FTensor::Tensor1< T, SPACE_DIM > &t_traction, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:86
FTensor::Index< 'k', SPACE_DIM > k
Definition: ContactOps.hpp:29
double gap0(FTensor::Tensor1< T, 3 > &t_coords, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:74
FTensor::Index< 'j', SPACE_DIM > j
Definition: ContactOps.hpp:28
double w(const double g, const double t)
Definition: ContactOps.hpp:100
double diff_constrains_dgap(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:110
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:65
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: ContactOps.hpp:27
double diff_constrains_dtraction(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:106
FTensor::Index< 'l', SPACE_DIM > l
Definition: ContactOps.hpp:30
double constrian(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:102
double gap(FTensor::Tensor1< T, SPACE_DIM > &t_disp, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:80
double sign(double x)
Definition: ContactOps.hpp:91
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
double gap0(double g, FTensor::Tensor1< T, 3 > &t_disp, FTensor::Tensor1< T, 3 > &t_normal)
double normal_traction(Tensor1< T1, 3 > &t_traction, Tensor1< T2, 3 > &t_normal)
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
constexpr double g
boost::shared_ptr< MatrixDouble > curlContactStressPtr
Definition: ContactOps.hpp:23
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ContactOps.hpp:13
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: ContactOps.hpp:14
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: ContactOps.hpp:15
boost::shared_ptr< MatrixDouble > contactTractionPtr
Definition: ContactOps.hpp:20
boost::shared_ptr< MatrixDouble > contactDispPtr
Definition: ContactOps.hpp:21
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: ContactOps.hpp:16
boost::shared_ptr< MatrixDouble > contactStressPtr
Definition: ContactOps.hpp:18
boost::shared_ptr< MatrixDouble > contactStressDivergencePtr
Definition: ContactOps.hpp:19
OpConstrainBoundaryLhs_dTraction(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:265
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:60
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:274
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:49
OpConstrainBoundaryLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:187
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:196
OpConstrainBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:114
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:38
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Definition: ContactOps.hpp:120
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column