v0.13.0
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
*
*/
/* 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;
OpSetContravariantPiolaTransformOnEdge2D;
};
template <> struct ElementsAndOps<3> {
static constexpr FieldSpace CONTACT_SPACE = HDIV;
OpHOSetContravariantPiolaTransformOnFace3D;
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
constexpr EntityType boundary_ent = SPACE_DIM == 3 ? MBTRI : MBEDGE;
FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
FormsIntegrators<BoundaryEleOp>::Assembly<PETSC>::OpBase;
//! [Operators used for contact]
using OpMixDivULhs = FormsIntegrators<DomainEleOp>::Assembly<
using OpLambdaGraULhs = FormsIntegrators<DomainEleOp>::Assembly<
using OpMixDivURhs = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
using OpMixLambdaGradURhs = FormsIntegrators<DomainEleOp>::Assembly<
using OpMixUTimesDivLambdaRhs = FormsIntegrators<DomainEleOp>::Assembly<
using OpMixUTimesLambdaRhs = FormsIntegrators<DomainEleOp>::Assembly<
using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<
//! [Operators used for contact]
//! [Body force]
using OpBodyForce = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpSource<1, SPACE_DIM>;
//! [Body force]
//! [Only used with Hooke equation (linear material model)]
using OpKCauchy = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
using OpInternalForceCauchy = FormsIntegrators<DomainEleOp>::Assembly<
//! [Only used with Hooke equation (linear material model)]
//! [Only used for dynamics]
using OpMass = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
//! [Only used for dynamics]
//! [Essential boundary conditions]
using OpBoundaryMass = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpBoundaryVec = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpBoundaryInternal = FormsIntegrators<BoundaryEleOp>::Assembly<
//! [Essential boundary conditions]
// Only used with Hencky/nonlinear material
using OpKPiola = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
using OpInternalForcePiola = FormsIntegrators<DomainEleOp>::Assembly<
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) {}
MoFEMErrorCode runProblem();
private:
MoFEMErrorCode setupProblem();
MoFEMErrorCode createCommonData();
MoFEMErrorCode tsSolve();
MoFEMErrorCode postProcess();
MoFEMErrorCode checkResults();
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]
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR bC();
CHKERR OPs();
CHKERR tsSolve();
CHKERR postProcess();
CHKERR checkResults();
}
//! [Run problem]
//! [Set up problem]
Simple *simple = mField.getInterface<Simple>();
// 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:
base = LASTBASE;
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>();
auto simple = mField.getInterface<Simple>();
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]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
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>();
if (SPACE_DIM == 2) {
pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
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 OpCalculateHOJacVolume(jac_ptr));
pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(
new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
pipeline.push_back(new OpSetHOInvJacToScalarBases(H1, 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.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));
pipeline.push_back(new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
"U", commonDataPtr->mStrainPtr, commonDataPtr->mStressPtr,
commonDataPtr->mDPtr));
pipeline.push_back(
new OpInternalForceCauchy("U", commonDataPtr->mStressPtr));
}
pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", commonDataPtr->contactDispPtr));
pipeline.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
"SIGMA", commonDataPtr->contactStressPtr));
pipeline.push_back(
new OpCalculateHVecTensorDivergence<SPACE_DIM, SPACE_DIM>(
"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(
new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>("U",
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());
pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", commonDataPtr->contactDispPtr));
pipeline.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
"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(
new OpConstrainBoundaryLhs_dU("SIGMA", "U", commonDataPtr));
pipeline.push_back(
new OpConstrainBoundaryLhs_dTraction("SIGMA", "SIGMA", commonDataPtr));
pipeline.push_back(new OpSpringLhs(
"U", "U",
[this](double, double, double) { return spring_stiffness; }
));
if (boundaryMarker)
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) {
for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
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]
Simple *simple = mField.getInterface<Simple>();
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
ISManager *is_manager = mField.getInterface<ISManager>();
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) {
SmartPetscObj<IS> is;
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),
SmartPetscObj<VecScatter>(scatter));
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, commonDataPtr, uXScatter, uYScatter, uZScatter));
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::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]
Simple *simple = m_field.getInterface<Simple>();
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
int main(int argc, char *argv[])
static char help[]
EntitiesFieldData::EntData EntData
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
constexpr double spring_stiffness
Definition: contact.cpp:139
ElementsAndOps< SPACE_DIM >::OpSetPiolaTransformOnBoundary OpSetPiolaTransformOnBoundary
Definition: contact.cpp:75
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixTensorTimesGradU< SPACE_DIM > OpMixLambdaGradURhs
Definition: contact.cpp:86
constexpr FieldSpace CONTACT_SPACE
Definition: contact.cpp:76
constexpr double rho
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
constexpr bool is_quasi_static
constexpr double young_modulus
@ 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 auto t_kd
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
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:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#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< 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
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:994
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
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
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
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
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:24
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
PetscBool is_large_strains
Definition: plastic.cpp:101
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:93
double cn
Definition: plastic.cpp:111
constexpr EntityType boundary_ent
Definition: plastic.cpp:54
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
[Example]
Range getEntsOnMeshSkin()
[Check]
Definition: contact.cpp:696
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
Definition: plastic.cpp:787
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:387
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode postProcess()
[Solve]
Definition: contact.cpp:689
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:319
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)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
Postprocess on face.
Post processing.
/* 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/>. */
/**
* \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);
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);
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);
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)
: AssemblyBoundaryEleOp(field_name, field_name, DomainEleOp::OPROW),
commonDataPtr(common_data_ptr) {}
OpConstrainBoundaryRhs::iNtegrate(EntitiesFieldData::EntData &data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &nf = AssemblyBoundaryEleOp::locF;
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(
const size_t nb_gauss_pts = getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
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(
const size_t nb_gauss_pts = getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
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
constexpr double alpha
double normal_traction(FTensor::Tensor1< T, SPACE_DIM > &t_traction, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:98
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:77
FTensor::Index< 'k', SPACE_DIM > k
Definition: ContactOps.hpp:41
double gap0(FTensor::Tensor1< T, 3 > &t_coords, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:86
FTensor::Index< 'j', SPACE_DIM > j
Definition: ContactOps.hpp:40
double w(const double g, const double t)
Definition: ContactOps.hpp:112
double diff_constrains_dgap(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:122
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: ContactOps.hpp:39
double diff_constrains_dtraction(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:118
FTensor::Index< 'l', SPACE_DIM > l
Definition: ContactOps.hpp:42
double constrian(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:114
double gap(FTensor::Tensor1< T, SPACE_DIM > &t_disp, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:92
double sign(double x)
Definition: ContactOps.hpp:103
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:76
constexpr double g
boost::shared_ptr< MatrixDouble > curlContactStressPtr
Definition: ContactOps.hpp:35
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ContactOps.hpp:25
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: ContactOps.hpp:26
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: ContactOps.hpp:27
boost::shared_ptr< MatrixDouble > contactTractionPtr
Definition: ContactOps.hpp:32
boost::shared_ptr< MatrixDouble > contactDispPtr
Definition: ContactOps.hpp:33
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: ContactOps.hpp:28
boost::shared_ptr< MatrixDouble > contactStressPtr
Definition: ContactOps.hpp:30
boost::shared_ptr< MatrixDouble > contactStressDivergencePtr
Definition: ContactOps.hpp:31
OpConstrainBoundaryLhs_dTraction(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:277
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:72
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:286
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:61
OpConstrainBoundaryLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:199
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:208
OpConstrainBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:126
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:50
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Definition: ContactOps.hpp:132