v0.14.0
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
* \CONTACT contact.cpp
*
* CONTACT of contact problem
*
* @copyright Copyright (c) 2023
*/
/* The above code is a preprocessor directive in C++ that checks if the macro
"EXECUTABLE_DIMENSION" has been defined. If it has not been defined, it is set
to 3" */
#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
#ifndef SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE 0
#endif
#include <MoFEM.hpp>
#ifdef PYTHON_SDF
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
using namespace MoFEM;
constexpr AssemblyType AT =
(SCHUR_ASSEMBLE) ? AssemblyType::SCHUR
: AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType IT =
IntegrationType::GAUSS; //< selected integration type
template <int DIM> struct ElementsAndOps;
static constexpr FieldSpace CONTACT_SPACE = HCURL;
};
static constexpr FieldSpace CONTACT_SPACE = HDIV;
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
/* The above code is defining an alias `EntData` for the type
`EntitiesFieldData::EntData`. This is a C++ syntax for creating a new name for
an existing type. */
//! [Specialisation for assembly]
// Assemble to A matrix, by default, however, some terms are assembled only to
// preconditioning.
template <>
const EntitiesFieldData::EntData &row_data,
return MatSetValues<AssemblyTypeSelector<SCHUR>>(
op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
};
template <>
const EntitiesFieldData::EntData &row_data,
return MatSetValues<AssemblyTypeSelector<SCHUR>>(
op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
};
/**
* @brief Element used to specialise assembly
*
*/
using BoundaryEleOp::BoundaryEleOp;
};
/**
* @brief Specialise assembly for Stabilised matrix
*
* @tparam
*/
template <>
const EntitiesFieldData::EntData &row_data,
return MatSetValues<AssemblyTypeSelector<SCHUR>>(
op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
};
//! [Specialisation for assembly]
//! [Operators used for contact]
IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
//! [Operators used for contact]
PetscBool is_quasi_static = PETSC_TRUE;
int order = 2;
int geom_order = 1;
double young_modulus = 100;
double poisson_ratio = 0.25;
double rho = 0.0;
double spring_stiffness = 0.0;
double vis_spring_stiffness = 0.0;
double alpha_damping = 0;
double scale = 1.;
PetscBool use_mfront = PETSC_FALSE;
PetscBool is_axisymmetric = PETSC_FALSE;
int atom_test = 0;
namespace ContactOps {
double cn_contact = 0.1;
}; // namespace ContactOps
//#define HENCKY_SMALL_STRAIN
#include <HenckyOps.hpp>
using namespace HenckyOps;
#include <ContactOps.hpp>
#include <PostProcContact.hpp>
#ifdef WITH_MODULE_MFRONT_INTERFACE
#include <MFrontMoFEMInterface.hpp>
#endif
using namespace ContactOps;
struct Contact {
Contact(MoFEM::Interface &m_field) : mField(m_field) {
mfrontInterface = nullptr;
}
private:
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
boost::shared_ptr<GenericElementInterface> mfrontInterface;
boost::shared_ptr<Monitor> monitorPtr;
#ifdef PYTHON_SDF
boost::shared_ptr<SDFPython> sdfPythonPtr;
#endif
struct ScaledTimeScale : public MoFEM::TimeScale {
double getScale(const double time) {
};
};
};
//! [Run problem]
}
//! [Run problem]
//! [Set up problem]
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
PETSC_NULL);
MOFEM_LOG("CONTACT", Sev::inform) << "Order " << order;
MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-use_mfront", &use_mfront,
PETSC_NULL);
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_axisymmetric",
&is_axisymmetric, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-atom_test", &atom_test, PETSC_NULL);
// 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("CONTACT", Sev::inform)
<< "Set AINSWORTH_LEGENDRE_BASE for displacements";
break;
case DEMKOWICZ:
MOFEM_LOG("CONTACT", Sev::inform)
<< "Set DEMKOWICZ_JACOBI_BASE for displacements";
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->addDataField("GEOMETRY", H1, base, SPACE_DIM);
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
auto get_skin = [&]() {
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
Skinner skin(&mField.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, body_ents, false, skin_ents);
return skin_ents;
};
auto filter_blocks = [&](auto skin) {
bool is_contact_block = false;
Range contact_range;
for (auto m :
(boost::format("%s(.*)") % "CONTACT").str()
))
) {
is_contact_block =
true; ///< bloks interation is collectibe, so that is set irrespective
///< if there are enerities in given rank or not in the block
MOFEM_LOG("CONTACT", Sev::inform)
<< "Find contact block set: " << m->getName();
auto meshset = m->getMeshset();
Range contact_meshset_range;
CHKERR mField.get_moab().get_entities_by_dimension(
meshset, SPACE_DIM - 1, contact_meshset_range, true);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
contact_meshset_range);
contact_range.merge(contact_meshset_range);
}
if (is_contact_block) {
MOFEM_LOG("SYNC", Sev::inform)
<< "Nb entities in contact surface: " << contact_range.size();
skin = intersect(skin, contact_range);
}
return skin;
};
auto filter_true_skin = [&](auto skin) {
Range boundary_ents;
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
CHKERR simple->setFieldOrder("SIGMA", 0);
CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
if (SPACE_DIM == 3) {
SETERRQ(PETSC_COMM_SELF, 1,
"Use executable contact_2d with axisymmetric model");
} else {
if (!use_mfront) {
SETERRQ(PETSC_COMM_SELF, 1,
"Axisymmetric model is only available with MFront (set "
"use_mfront to 1)");
} else {
MOFEM_LOG("CONTACT", Sev::inform) << "Using axisymmetric model";
}
}
} else {
if (SPACE_DIM == 2) {
MOFEM_LOG("CONTACT", Sev::inform) << "Using plane strain model";
}
}
if (!use_mfront) {
CHKERR simple->setUp();
} else {
#ifndef WITH_MODULE_MFRONT_INTERFACE
SETERRQ(
PETSC_COMM_SELF, 1,
"MFrontInterface module was not found while use_mfront was set to 1");
#else
SETERRQ(PETSC_COMM_SELF, 1,
"MFrontInterface module is not compatible with Schur assembly");
}
if (SPACE_DIM == 3) {
boost::make_shared<MFrontMoFEMInterface<TRIDIMENSIONAL>>(
mField, "U", "GEOMETRY", true, is_quasi_static);
} else if (SPACE_DIM == 2) {
boost::make_shared<MFrontMoFEMInterface<AXISYMMETRICAL>>(
mField, "U", "GEOMETRY", true, is_quasi_static);
} else {
mfrontInterface = boost::make_shared<MFrontMoFEMInterface<PLANESTRAIN>>(
mField, "U", "GEOMETRY", true, is_quasi_static);
}
}
#endif
CHKERR mfrontInterface->getCommandLineParameters();
CHKERR mfrontInterface->addElementFields();
CHKERR mfrontInterface->createElements();
CHKERR simple->defineFiniteElements();
CHKERR simple->defineProblem(PETSC_TRUE);
CHKERR simple->buildFields();
CHKERR simple->buildFiniteElements();
CHKERR mfrontInterface->addElementsToDM(simple->getDM());
CHKERR simple->buildProblem();
}
auto dm = simple->getDM();
monitorPtr = boost::make_shared<Monitor>(dm, use_mfront, mfrontInterface,
if (use_mfront) {
mfrontInterface->setMonitorPtr(monitorPtr);
}
auto project_ho_geometry = [&]() {
Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
return mField.loop_dofs("GEOMETRY", ent_method);
};
PetscBool project_geometry = PETSC_TRUE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
&project_geometry, PETSC_NULL);
if (project_geometry){
CHKERR project_ho_geometry();
}
} //! [Set up problem]
//! [Create common data]
auto get_options = [&]() {
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
&young_modulus, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
&poisson_ratio, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn_contact,
PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-spring_stiffness",
&spring_stiffness, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-vis_spring_stiffness",
&vis_spring_stiffness, PETSC_NULL);
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping",
&alpha_damping, PETSC_NULL);
if (!use_mfront) {
MOFEM_LOG("CONTACT", Sev::inform) << "Young modulus " << young_modulus;
MOFEM_LOG("CONTACT", Sev::inform) << "Poisson_ratio " << poisson_ratio;
} else {
MOFEM_LOG("CONTACT", Sev::inform) << "Using MFront for material model";
}
MOFEM_LOG("CONTACT", Sev::inform) << "Density " << rho;
MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << cn_contact;
MOFEM_LOG("CONTACT", Sev::inform)
<< "Spring stiffness " << spring_stiffness;
MOFEM_LOG("CONTACT", Sev::inform)
<< "Vis spring_stiffness " << vis_spring_stiffness;
MOFEM_LOG("CONTACT", Sev::inform) << "alpha_damping " << alpha_damping;
PetscBool use_scale = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-use_scale", &use_scale,
PETSC_NULL);
if (use_scale) {
}
MOFEM_LOG("CONTACT", Sev::inform) << "Scale " << scale;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_quasi_static",
&is_quasi_static, PETSC_NULL);
MOFEM_LOG("CONTACT", Sev::inform)
<< "Is quasi-static: " << (is_quasi_static ? "true" : "false");
};
CHKERR get_options();
#ifdef PYTHON_SDF
char sdf_file_name[255];
CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-sdf_file",
sdf_file_name, 255, PETSC_NULL);
sdfPythonPtr = boost::make_shared<SDFPython>();
CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
sdfPythonWeakPtr = sdfPythonPtr;
#endif
}
//! [Create common data]
//! [Boundary condition]
auto bc_mng = mField.getInterface<BcManager>();
for (auto f : {"U", "SIGMA"}) {
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_X", f, 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_Y", f, 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_Z", f, 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_ALL", f, 0, 3);
}
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
"SIGMA", 0, 0, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
"SIGMA", 1, 1, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
"SIGMA", 2, 2, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
"SIGMA", 0, 3, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(), "REMOVE_X", "SIGMA", 0, 3, false, true);
}
// Note remove has to be always before push. Then node marking will be
// corrupted.
CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
}
//! [Boundary condition]
//! [Push operators to pip]
auto *pip_mng = mField.getInterface<PipelineManager>();
auto bc_mng = mField.getInterface<BcManager>();
auto time_scale = boost::make_shared<ScaledTimeScale>();
auto body_force_time_scale =
boost::make_shared<ScaledTimeScale>("body_force_hist.txt");
auto integration_rule_vol = [](int, int, int approx_order) {
return 2 * approx_order + geom_order - 1;
};
auto integration_rule_boundary = [](int, int, int approx_order) {
return 2 * approx_order + geom_order - 1;
};
auto add_domain_base_ops = [&](auto &pip) {
"GEOMETRY");
};
auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
henky_common_data_ptr->matDPtr = boost::make_shared<MatrixDouble>();
henky_common_data_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
auto add_domain_ops_lhs = [&](auto &pip) {
//! [Only used for dynamics]
//! [Only used for dynamics]
if (is_quasi_static == PETSC_FALSE) {
auto *pip_mng = mField.getInterface<PipelineManager>();
auto fe_domain_lhs = pip_mng->getDomainLhsFE();
auto get_inertia_and_mass_damping =
[this, fe_domain_lhs](const double, const double, const double) {
return (rho * scale) * fe_domain_lhs->ts_aa +
(alpha_damping * scale) * fe_domain_lhs->ts_a;
};
pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
} else {
auto *pip_mng = mField.getInterface<PipelineManager>();
auto fe_domain_lhs = pip_mng->getDomainLhsFE();
auto get_inertia_and_mass_damping =
[this, fe_domain_lhs](const double, const double, const double) {
return (alpha_damping * scale) * fe_domain_lhs->ts_a;
};
pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
}
if (!use_mfront) {
CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
}
};
auto add_domain_ops_rhs = [&](auto &pip) {
pip, mField, "U", {body_force_time_scale}, Sev::inform);
//! [Only used for dynamics]
//! [Only used for dynamics]
// only in case of dynamics
if (is_quasi_static == PETSC_FALSE) {
auto mat_acceleration = boost::make_shared<MatrixDouble>();
"U", mat_acceleration));
pip.push_back(
new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
return rho * scale;
}));
}
// only in case of viscosity
if (alpha_damping > 0) {
auto mat_velocity = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(
new OpInertiaForce("U", mat_velocity, [](double, double, double) {
return alpha_damping * scale;
}));
}
if (!use_mfront) {
CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
}
CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
pip, "SIGMA", "U", is_axisymmetric);
};
auto add_boundary_base_ops = [&](auto &pip) {
"GEOMETRY");
// We have to integrate on curved face geometry, thus integration weight
// have to adjusted.
pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
};
auto add_boundary_ops_lhs = [&](auto &pip) {
//! [Operators used for contact]
//! [Operators used for contact]
// Add Natural BCs to LHS
pip, mField, "U", Sev::inform);
auto *pip_mng = mField.getInterface<PipelineManager>();
auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
pip.push_back(new OpSpringLhs(
"U", "U",
[this, fe_boundary_lhs](double, double, double) {
(vis_spring_stiffness * scale) * fe_boundary_lhs->ts_a;
}
));
}
ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
pip, "SIGMA", "U", is_axisymmetric);
mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
integration_rule_vol, is_axisymmetric);
};
auto add_boundary_ops_rhs = [&](auto &pip) {
//! [Operators used for contact]
//! [Operators used for contact]
// Add Natural BCs to RHS
pip, mField, "U", {time_scale}, Sev::inform);
auto u_disp = boost::make_shared<MatrixDouble>();
auto dot_u_disp = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
pip.push_back(
pip.push_back(
new OpSpringRhs("U", u_disp, [this](double, double, double) {
}));
pip.push_back(
new OpSpringRhs("U", dot_u_disp, [this](double, double, double) {
}));
}
ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
pip, "SIGMA", "U", is_axisymmetric);
};
CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
if (use_mfront) {
if (is_quasi_static == PETSC_TRUE)
mfrontInterface->setOperators();
mfrontInterface->setupSolverFunctionTS(t_type);
mfrontInterface->setupSolverJacobianTS(t_type);
}
CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
}
//! [Push operators to pip]
//! [Solve]
struct SetUpSchur {
static boost::shared_ptr<SetUpSchur>
virtual MoFEMErrorCode setUp(SmartPetscObj<TS> solver) = 0;
protected:
SetUpSchur() = default;
};
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) {
monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
boost::shared_ptr<ForcesAndSourcesCore> null;
CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
monitorPtr, null, null);
};
auto set_essential_bc = [&]() {
// This is low level pushing finite elements (pipelines) to solver
auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
auto pre_proc_ptr = boost::make_shared<FEMethod>();
auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
// Add boundary condition scaling
auto time_scale = boost::make_shared<TimeScale>();
auto get_bc_hook_rhs = [&]() {
{time_scale}, false);
return hook;
};
pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
auto get_post_proc_hook_rhs = [&]() {
mField, post_proc_rhs_ptr, 1.);
};
auto get_post_proc_hook_lhs = [&]() {
mField, post_proc_lhs_ptr, 1.);
};
post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
if (AT != AssemblyType::SCHUR) {
post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
}
};
auto set_schur_pc = [&](auto solver) {
boost::shared_ptr<SetUpSchur> schur_ptr;
if (AT == AssemblyType::SCHUR) {
CHK_MOAB_THROW(schur_ptr->setUp(solver), "SetUpSchur::setUp");
}
return schur_ptr;
};
auto dm = simple->getDM();
auto D = createDMVector(dm);
uXScatter = scatter_create(D, 0);
uYScatter = scatter_create(D, 1);
if (SPACE_DIM == 3)
uZScatter = scatter_create(D, 2);
// Add extra finite elements to SNES solver pipelines to resolve essential
// boundary conditions
CHKERR set_essential_bc();
if (is_quasi_static == PETSC_TRUE) {
auto solver = pip_mng->createTSIM();
CHKERR TSSetFromOptions(solver);
auto D = createDMVector(dm);
auto schur_pc_ptr = set_schur_pc(solver);
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetSolution(solver, D);
CHKERR TSSetUp(solver);
CHKERR TSSolve(solver, NULL);
} else {
auto solver = pip_mng->createTSIM2();
CHKERR TSSetFromOptions(solver);
auto dm = simple->getDM();
auto D = createDMVector(dm);
auto DD = vectorDuplicate(D);
auto schur_pc_ptr = set_schur_pc(solver);
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TS2SetSolution(solver, D, DD);
CHKERR TSSetUp(solver);
CHKERR TSSolve(solver, NULL);
}
}
//! [Solve]
//! [Check]
//! [Check]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
#ifdef PYTHON_SDF
Py_Initialize();
np::initialize();
#endif
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Add logging channel for CONTACT
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
LogManager::setLog("CONTACT");
MOFEM_LOG_TAG("CONTACT", "Indent");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database interface
//! [Create MoFEM]
//! [Load mesh]
CHKERR simple->getOptions();
CHKERR simple->loadFile("");
//! [Load mesh]
//! [CONTACT]
Contact ex(m_field);
CHKERR ex.runProblem();
//! [CONTACT]
}
#ifdef PYTHON_SDF
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
return 0;
}
struct SetUpSchurImpl : public SetUpSchur {
SetUpSchurImpl(MoFEM::Interface &m_field) : SetUpSchur(), mField(m_field) {}
virtual ~SetUpSchurImpl() {
A.reset();
P.reset();
S.reset();
}
private:
};
auto dm = simple->getDM();
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR KSPSetFromOptions(ksp);
PC pc;
CHKERR KSPGetPC(ksp, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
MOFEM_LOG("CONTACT", Sev::inform) << "Setup Schur pc";
if (A || P || S) {
"It is expected that Schur matrix is not allocated. This is "
"possible only if PC is set up twice");
}
P = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
CHKERR MatSetBlockSize(S, SPACE_DIM);
auto ts_ctx_ptr = getDMTsCtx(dm);
CHKERR TSSetIJacobian(solver, A, P, TsSetIJacobian, ts_ctx_ptr.get());
} else {
MOFEM_LOG("CONTACT", Sev::inform) << "No Schur pc";
pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
pip->getOpBoundaryLhsPipeline().push_back(
new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
pip->getOpDomainLhsPipeline().push_back(
new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
post_proc_schur_lhs_ptr->postProcessHook = [this,
post_proc_schur_lhs_ptr]() {
mField, post_proc_schur_lhs_ptr, 1.)();
};
auto ts_ctx_ptr = getDMTsCtx(dm);
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
}
}
auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
auto set_up = [&]() {
CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "SUB");
CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
CHKERR DMSetUp(sub_dm);
};
CHK_THROW_MESSAGE(set_up(), "sey up dm");
return sub_dm;
}
double eps_stab = 1e-4;
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
PETSC_NULL);
using B =
using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
// Boundary
auto dm_is = getDMSubData(subDM)->getSmartRowIs();
auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
pip->getOpBoundaryLhsPipeline().push_back(
new OpMassStab("SIGMA", "SIGMA",
[eps_stab](double, double, double) { return eps_stab; }));
pip->getOpBoundaryLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DGESV>(
{"SIGMA"}, {nullptr}, {ao_up}, {S}, {false}, false));
// Domain
pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DGESV>(
{"SIGMA"}, {nullptr}, {ao_up}, {S}, {false}, false));
auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
CHKERR MatZeroEntries(A);
CHKERR MatZeroEntries(P);
CHKERR MatZeroEntries(S);
MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Begin";
};
post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
post_proc_schur_lhs_ptr]() {
MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble End";
*post_proc_schur_lhs_ptr->matAssembleSwitch = false;
auto print_mat_norm = [this](auto a, std::string prefix) {
double nrm;
CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
MOFEM_LOG("CONTACT", Sev::noisy) << prefix << " norm = " << nrm;
};
CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
mField, post_proc_schur_lhs_ptr, 1, A)();
CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
CHKERR MatAXPY(P, 1, A, SAME_NONZERO_PATTERN);
CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
#ifndef NDEBUG
CHKERR print_mat_norm(A, "A");
CHKERR print_mat_norm(P, "P");
CHKERR print_mat_norm(S, "S");
#endif // NDEBUG
MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Finish";
};
auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
}
mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
simple->getProblemName(), ROW, "SIGMA", 0, SPACE_DIM, is);
CHKERR PCFieldSplitSetIS(pc, NULL, is);
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
}
boost::shared_ptr<SetUpSchur>
return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
}
std::string param_file
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static char help[]
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
double young_modulus
Definition: contact.cpp:127
constexpr AssemblyType AT
Definition: contact.cpp:34
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
double spring_stiffness
Definition: contact.cpp:130
double rho
Definition: contact.cpp:129
int atom_test
Definition: contact.cpp:139
PetscBool is_quasi_static
[Operators used for contact]
Definition: contact.cpp:123
double alpha_damping
Definition: contact.cpp:132
PetscBool use_mfront
Definition: contact.cpp:136
constexpr int SPACE_DIM
Definition: contact.cpp:53
double poisson_ratio
Definition: contact.cpp:128
double vis_spring_stiffness
Definition: contact.cpp:131
double scale
Definition: contact.cpp:134
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT >::OpMass< 1, SPACE_DIM > OpSpringLhs
[Operators used for contact]
Definition: contact.cpp:118
int order
Definition: contact.cpp:125
PetscBool is_axisymmetric
Definition: contact.cpp:137
int geom_order
Definition: contact.cpp:126
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:120
constexpr FieldSpace CONTACT_SPACE
[Specialisation for assembly]
Definition: contact.cpp:114
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
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
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ 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 int order
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:219
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:442
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:988
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
double D
const double v
phase velocity of light in medium (cm/ns)
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
MoFEMErrorCode opFactoryBoundaryToDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string fe_domain_name, std::string sigma, std::string u, std::string geom, ForcesAndSourcesCore::RuleHookFun rule, bool is_axisymmetric=false)
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:1042
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1061
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
double young_modulus
Young modulus.
Definition: plastic.cpp:172
constexpr AssemblyType AT
Definition: plastic.cpp:44
constexpr IntegrationType IT
Definition: plastic.cpp:47
double rho
Definition: plastic.cpp:191
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
PetscBool is_quasi_static
Definition: plastic.cpp:190
double alpha_damping
Definition: plastic.cpp:192
double scale
Definition: plastic.cpp:170
int geom_order
Order if fixed.
Definition: plastic.cpp:188
constexpr FieldSpace CONTACT_SPACE
Definition: plastic.cpp:52
static constexpr int approx_order
Element used to specialise assembly.
Definition: contact.cpp:94
double getScale(const double time)
Get scaling at a given time.
Definition: contact.cpp:200
boost::shared_ptr< Monitor > monitorPtr
Definition: contact.cpp:192
MoFEMErrorCode tsSolve()
Definition: contact.cpp:760
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: contact.cpp:188
MoFEMErrorCode runProblem()
[Run problem]
Definition: contact.cpp:207
MoFEMErrorCode setupProblem()
[Run problem]
Definition: contact.cpp:220
MoFEMErrorCode OPs()
[Boundary condition]
Definition: contact.cpp:526
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: contact.cpp:414
MoFEMErrorCode checkResults()
[Solve]
Definition: contact.cpp:900
boost::shared_ptr< GenericElementInterface > mfrontInterface
Definition: contact.cpp:191
MoFEM::Interface & mField
Definition: contact.cpp:178
MoFEMErrorCode bC()
[Create common data]
Definition: contact.cpp:484
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: contact.cpp:187
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: contact.cpp:189
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:28
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:30
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Managing BitRefLevels.
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:76
Data on single entity (This is passed as argument to DataOperator::doWork)
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition: Essential.hpp:33
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition: Essential.hpp:41
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:25
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition: Natural.hpp:65
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Clear Schur complement internal data.
Definition: Schur.hpp:22
Assemble Schur complement.
Definition: Schur.hpp:104
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
intrusive_ptr for managing petsc objects
Force scale operator for reading two columns.
double getScale(const double time)
Get scaling at a given time.
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
Definition: plastic.cpp:762
SetUpSchur()=default
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
Create data structure for handling Schur complement.
Definition: plastic.cpp:1482
virtual MoFEMErrorCode setUp(TS solver)=0
SmartPetscObj< Mat > A
Definition: contact.cpp:983
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1304
MoFEMErrorCode setUp(KSP solver)
SmartPetscObj< Mat > S
SmartPetscObj< Mat > P
SmartPetscObj< DM > createSubDM()
Definition: contact.cpp:1056
MoFEMErrorCode setEntities()
Definition: elastic.cpp:942
virtual ~SetUpSchurImpl()
MoFEMErrorCode setPC(PC pc)
Definition: contact.cpp:1160
MoFEMErrorCode setOperator()
Definition: contact.cpp:1072
MoFEM::Interface & mField
/**
* \file ContactOps.hpp
* \example ContactOps.hpp
*/
#ifndef __CONTACTOPS_HPP__
#define __CONTACTOPS_HPP__
namespace ContactOps {
//! [Common data]
struct CommonData : public boost::enable_shared_from_this<CommonData> {
// MatrixDouble contactStress;
MatrixDouble contactTraction;
MatrixDouble contactDisp;
VectorDouble sdfVals; ///< size is equal to number of gauss points on element
MatrixDouble gradsSdf; ///< nb of rows is equals to dimension, and nb of cols
///< is equals to number of gauss points on element
MatrixDouble hessSdf; ///< nb of rows is equals to nb of element of symmetric
///< matrix, and nb of cols is equals to number of gauss
///< points on element
static SmartPetscObj<Vec>
totalTraction; // User have to release and create vector when appropiate.
static auto createTotalTraction(MoFEM::Interface &m_field) {
constexpr int ghosts[] = {0, 1, 2};
createGhostVector(m_field.get_comm(),
(m_field.get_comm_rank() == 0) ? 3 : 0, 3,
(m_field.get_comm_rank() == 0) ? 0 : 3, ghosts);
return totalTraction;
}
static auto getFTensor1TotalTraction() {
const double *t_ptr;
CHK_THROW_MESSAGE(VecGetArrayRead(CommonData::totalTraction, &t_ptr),
"get array");
FTensor::Tensor1<double, 3> t{t_ptr[0], t_ptr[1], t_ptr[2]};
CHK_THROW_MESSAGE(VecRestoreArrayRead(CommonData::totalTraction, &t_ptr),
"restore array");
return t;
} else {
return FTensor::Tensor1<double, 3>{0., 0., 0.};
}
}
inline auto contactTractionPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
inline auto contactDispPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactDisp);
}
inline auto sdfPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &sdfVals);
}
inline auto gradSdfPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradsSdf);
}
inline auto hessSdfPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &hessSdf);
}
};
SmartPetscObj<Vec> CommonData::totalTraction;
//! [Common data]
//! [Surface distance function from python]
#ifdef PYTHON_SDF
struct SDFPython {
SDFPython() = default;
virtual ~SDFPython() = default;
MoFEMErrorCode sdfInit(const std::string py_file) {
try {
// create main module
auto main_module = bp::import("__main__");
mainNamespace = main_module.attr("__dict__");
bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
// create a reference to python function
sdfFun = mainNamespace["sdf"];
sdfGradFun = mainNamespace["grad_sdf"];
sdfHessFun = mainNamespace["hess_sdf"];
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
};
template <typename T>
inline std::vector<T>
py_list_to_std_vector(const boost::python::object &iterable) {
return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
boost::python::stl_input_iterator<T>());
}
MoFEMErrorCode evalSdf(
double t, double x, double y, double z, double tx, double ty, double tz,
double &sdf
) {
try {
// call python function
sdf = bp::extract<double>(sdfFun(t, x, y, z, tx, ty, tz));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
}
MoFEMErrorCode evalGradSdf(
double t, double x, double y, double z, double tx, double ty, double tz,
std::vector<double> &grad_sdf
) {
try {
// call python function
py_list_to_std_vector<double>(sdfGradFun(t, x, y, z, tx, ty, tz));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
if (grad_sdf.size() != 3) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Expected size 3");
}
}
MoFEMErrorCode evalHessSdf(
double t, double x, double y, double z, double tx, double ty, double tz,
std::vector<double> &hess_sdf
) {
try {
// call python function
py_list_to_std_vector<double>(sdfHessFun(t, x, y, z, tx, ty, tz));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
}
if (hess_sdf.size() != 6) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Expected size 6");
}
}
private:
bp::object mainNamespace;
bp::object sdfFun;
bp::object sdfGradFun;
bp::object sdfHessFun;
};
static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
#endif
//! [Surface distance function from python]
using SurfaceDistanceFunction = boost::function<double(
double t, double x, double y, double z, double tx, double ty, double tz)>;
using GradSurfaceDistanceFunction = boost::function<FTensor::Tensor1<double, 3>(
double t, double x, double y, double z, double tx, double ty, double tz)>;
boost::function<FTensor::Tensor2_symmetric<double, 3>(
double t, double x, double y, double z, double tx, double ty,
double tz)>;
inline double surface_distance_function(double t, double x, double y, double z,
double tx, double ty, double tz) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
double sdf;
CHK_MOAB_THROW(sdf_ptr->evalSdf(t, x, y, z, tx, ty, tz, sdf),
"Failed python call");
return sdf;
}
#endif
return y + 0.5;
}
grad_surface_distance_function(double t, double x, double y, double z,
double tx, double ty, double tz) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
std::vector<double> grad_sdf;
CHK_MOAB_THROW(sdf_ptr->evalGradSdf(t, x, y, z, tx, ty, tz, grad_sdf),
"Failed python call");
}
#endif
return FTensor::Tensor1<double, 3>{0., 1., 0.};
}
hess_surface_distance_function(double t, double x, double y, double z,
double tx, double ty, double tz) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
std::vector<double> hess_sdf;
CHK_MOAB_THROW(sdf_ptr->evalHessSdf(t, x, y, z, tx, ty, tz, hess_sdf),
"Failed python call");
hess_sdf[4], hess_sdf[5]};
}
#endif
return FTensor::Tensor2_symmetric<double, 3>{0., 0., 0., 0., 0., 0.};
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl;
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpEvaluateSDFImpl;
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryRhsImpl;
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dTractionImpl;
template <int DIM, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>
: public BoundaryEleOp {
OpAssembleTotalContactTractionImpl(
boost::shared_ptr<CommonData> common_data_ptr, double scale = 1);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
const double scaleTraction;
};
template <int DIM, typename BoundaryEleOp>
struct OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp> : public BoundaryEleOp {
OpEvaluateSDFImpl(boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
};
template <int DIM, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dUImpl(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);
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dTractionImpl(
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);
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <typename BoundaryEleOp> struct ContactIntegrators {
template <int DIM, IntegrationType I>
OpAssembleTotalContactTractionImpl<DIM, I, BoundaryEleOp>;
template <int DIM, IntegrationType I>
using OpEvaluateSDF = OpEvaluateSDFImpl<DIM, I, BoundaryEleOp>;
template <AssemblyType A> struct Assembly {
typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
template <int DIM, IntegrationType I>
OpConstrainBoundaryRhsImpl<DIM, I, AssemblyBoundaryEleOp>;
template <int DIM, IntegrationType I>
OpConstrainBoundaryLhs_dUImpl<DIM, I, AssemblyBoundaryEleOp>;
template <int DIM, IntegrationType I>
OpConstrainBoundaryLhs_dTractionImpl<DIM, I, AssemblyBoundaryEleOp>;
};
};
inline double sign(double x) {
if (x == 0)
return 0;
else if (x > 0)
return 1;
else
return -1;
};
inline double w(const double sdf, const double tn) {
return sdf - cn_contact * tn;
}
/**
* @brief constrain function
*
* return 1 if negative sdn or positive tn
*
* @param sdf signed distance
* @param tn traction
* @return double
*/
inline double constrain(double sdf, double tn) {
const auto s = sign(w(sdf, tn));
return (1 - s) / 2;
}
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::
OpAssembleTotalContactTractionImpl(
boost::shared_ptr<CommonData> common_data_ptr, double scale)
commonDataPtr(common_data_ptr), scaleTraction(scale) {}
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
int side, EntityType type, EntData &data) {
FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
const double alpha = t_w * BoundaryEleOp::getMeasure();
t_sum_t(i) += alpha * t_traction(i);
++t_w;
++t_traction;
}
t_sum_t(i) *= scaleTraction;
constexpr int ind[] = {0, 1, 2};
CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
ADD_VALUES);
}
template <int DIM, typename BoundaryEleOp>
OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::OpEvaluateSDFImpl(
boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {}
template <int DIM, typename BoundaryEleOp>
OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
int side, EntityType type, EntData &data) {
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
auto &sdf_vec = commonDataPtr->sdfVals;
auto &grad_mat = commonDataPtr->gradsSdf;
auto &hess_mat = commonDataPtr->hessSdf;
sdf_vec.resize(nb_gauss_pts, false);
grad_mat.resize(DIM, nb_gauss_pts, false);
hess_mat.resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
auto t_total_traction = CommonData::getFTensor1TotalTraction();
auto t_sdf = getFTensor0FromVec(sdf_vec);
auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto next = [&]() {
++t_sdf;
++t_grad_sdf;
++t_hess_sdf;
++t_disp;
++t_coords;
};
auto ts_time = BoundaryEleOp::getTStime();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
t_spatial_coords(i) = t_coords(i) + t_disp(i);
auto sdf_v = surfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto t_grad_sdf_v = gradSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto t_hess_sdf_v = hessSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
t_sdf = sdf_v;
t_grad_sdf(i) = t_grad_sdf_v(i);
t_hess_sdf(i, j) = t_hess_sdf_v(i, j);
next();
}
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
OpConstrainBoundaryRhsImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
EntitiesFieldData::EntData &data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto t_total_traction = CommonData::getFTensor1TotalTraction();
auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = AssemblyBoundaryEleOp::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) {
t_normal(i) =
t_normal_at_pts(i) / std::sqrt(t_normal_at_pts(i) * t_normal_at_pts(i));
auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
t_spatial_coords(i) = t_coords(i) + t_disp(i);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto sdf = surfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto t_grad_sdf = gradSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(sdf, tn);
t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
t_rhs(i) =
t_cQ(i, j) * (t_disp(j) - cn_contact * t_traction(j))
+
t_cP(i, j) * t_disp(j) +
c * (sdf * t_grad_sdf(i)); // add gap0 displacements
size_t bb = 0;
for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
const double beta = alpha * (t_base(i) * t_normal(i));
t_nf(i) -= beta * t_rhs(i);
++t_nf;
++t_base;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
++t_normal_at_pts;
}
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
OpConstrainBoundaryLhs_dUImpl(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,
commonDataPtr(common_data_ptr) {
AssemblyBoundaryEleOp::sYmm = false;
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto t_total_traction = CommonData::getFTensor1TotalTraction();
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
t_normal(i) =
t_normal_at_pts(i) / std::sqrt(t_normal_at_pts(i) * t_normal_at_pts(i));
const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
t_spatial_coords(i) = t_coords(i) + t_disp(i);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto sdf = surfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto t_grad_sdf = gradSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(sdf, tn);
t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
if (c > 0) {
auto t_hess_sdf = hessSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1),
t_spatial_coords(2), t_total_traction(0), t_total_traction(1),
t_total_traction(2));
t_res_dU(i, j) +=
(c * cn_contact) *
(t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
c * sdf * t_hess_sdf(i, j);
}
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, 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 / DIM; ++cc) {
const double beta = alpha * row_base * t_col_base;
t_mat(i, j) -= beta * t_res_dU(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;
++t_normal_at_pts;
}
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
OpConstrainBoundaryLhs_dTractionImpl(
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,
commonDataPtr(common_data_ptr) {
AssemblyBoundaryEleOp::sYmm = false;
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto t_total_traction = CommonData::getFTensor1TotalTraction();
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_w = AssemblyBoundaryEleOp::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) {
t_normal(i) =
t_normal_at_pts(i) / std::sqrt(t_normal_at_pts(i) * t_normal_at_pts(i));
const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
t_spatial_coords(i) = t_coords(i) + t_disp(i);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto sdf = surfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto t_grad_sdf = gradSurfaceDistanceFunction(
ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
t_total_traction(0), t_total_traction(1), t_total_traction(2));
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(sdf, tn);
t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
t_res_dt(i, j) = -cn_contact * t_cQ(i, j);
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, 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 / 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 * t_res_dt(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;
++t_normal_at_pts;
}
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u) {
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template LinearForm<I>;
using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
typename B::template OpGradTimesTensor<1, DIM, DIM>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
auto div_stress_ptr = boost::make_shared<MatrixDouble>();
auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(
new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
pip.push_back(
new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
pip.push_back(
new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(),
[](double, double, double) constexpr { return 1; }));
pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr));
pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr));
pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr));
}
template <typename OpMixLhs> struct OpMixLhsSide : public OpMixLhs {
using OpMixLhs::OpMixLhs;
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
// Only assemble side which correspond to edge entity on boundary
if (side_fe_entity == side_fe_data) {
CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
col_data);
}
}
};
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
MoFEM::Interface &m_field,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string fe_domain_name, std::string sigma, std::string u,
std::string geom, ForcesAndSourcesCore::RuleHookFun rule) {
using DomainEleOp = typename DomainEle::UserDataOperator;
auto op_loop_side = new OpLoopSide<DomainEle>(
m_field, fe_domain_name, DIM, Sev::noisy,
boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
pip.push_back(op_loop_side);
CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
{H1, HDIV}, geom);
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template BiLinearForm<I>;
using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
auto unity = []() { return 1; };
op_loop_side->getOpPtrVector().push_back(
new OpMixDivULhsSide(sigma, u, unity, true));
op_loop_side->getOpPtrVector().push_back(
new OpLambdaGraULhsSide(sigma, u, unity, true));
op_loop_side->getSideFEPtr()->getRuleHook = rule;
}
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u) {
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
sigma, common_data_ptr->contactTractionPtr()));
pip.push_back(
new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
DIM, GAUSS>(sigma, u, common_data_ptr));
pip.push_back(new typename C::template Assembly<A>::
template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
sigma, sigma, common_data_ptr));
}
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u) {
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
sigma, common_data_ptr->contactTractionPtr()));
pip.push_back(
new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
DIM, GAUSS>(sigma, common_data_ptr));
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma) {
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
sigma, common_data_ptr->contactTractionPtr()));
pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
common_data_ptr, 1. / scale));
}
}; // namespace ContactOps
#endif // __CONTACTOPS_HPP__
Kronecker Delta class.
@ NOSPACE
Definition: definitions.h:83
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
constexpr auto t_kd
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, string field_name, Pip &pip, std::string block_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, Sev sev=Sev::inform)
Assemble the left hand side, i.e. tangent matrix.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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 >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
MatrixDouble grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:273
MatrixDouble hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:330
MoFEMErrorCode opFactoryBoundaryRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
double w(const double sdf, const double tn)
Definition: ContactOps.hpp:559
boost::function< MatrixDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> HessSurfaceDistanceFunction
Definition: ContactOps.hpp:214
double constrain(double sdf, double tn)
constrain function
Definition: ContactOps.hpp:572
VectorDouble surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:216
boost::function< VectorDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> SurfaceDistanceFunction
[Common data]
Definition: ContactOps.hpp:206
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
boost::function< MatrixDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> GradSurfaceDistanceFunction
Definition: ContactOps.hpp:210
MoFEMErrorCode opFactoryCalculateTraction(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, bool is_axisymmetric=false)
double sign(double x)
Definition: ContactOps.hpp:549
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: sdf.py:1
def grad_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:15
def hess_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:19
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr AssemblyType A
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
VectorDouble sdfVals
size is equal to number of gauss points on element
Definition: ContactOps.hpp:19
MatrixDouble contactTraction
Definition: ContactOps.hpp:16
static auto getFTensor1TotalTraction()
Definition: ContactOps.hpp:41
MatrixDouble contactDisp
Definition: ContactOps.hpp:17
MatrixDouble hessSdf
Definition: ContactOps.hpp:22
MatrixDouble gradsSdf
Definition: ContactOps.hpp:20
OpConstrainBoundaryLhs_dUImpl< DIM, I, AssemblyBoundaryEleOp > OpConstrainBoundaryLhs_dU
Definition: ContactOps.hpp:541
OpConstrainBoundaryLhs_dTractionImpl< DIM, I, AssemblyBoundaryEleOp > OpConstrainBoundaryLhs_dTraction
Definition: ContactOps.hpp:545
OpConstrainBoundaryRhsImpl< DIM, I, AssemblyBoundaryEleOp > OpConstrainBoundaryRhs
Definition: ContactOps.hpp:537
OpEvaluateSDFImpl< DIM, I, BoundaryEleOp > OpEvaluateSDF
Definition: ContactOps.hpp:528
OpAssembleTotalContactTractionImpl< DIM, I, BoundaryEleOp > OpAssembleTotalContactTraction
Definition: ContactOps.hpp:525
virtual int get_comm_rank() const =0
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
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