v0.14.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
* \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>
#include <MatrixFunction.hpp>
using namespace MoFEM;
#include <GenericElementInterface.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
constexpr AssemblyType AT =
(SCHUR_ASSEMBLE) ? AssemblyType::BLOCK_SCHUR
: AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType IT =
IntegrationType::GAUSS; //< selected integration type
template <int DIM> struct ElementsAndOps;
template <> struct ElementsAndOps<2> : PipelineManager::ElementsAndOpsByDim<2> {
static constexpr FieldSpace CONTACT_SPACE = HCURL;
};
template <> struct ElementsAndOps<3> : PipelineManager::ElementsAndOpsByDim<3> {
static constexpr FieldSpace CONTACT_SPACE = HDIV;
};
constexpr FieldSpace ElementsAndOps<2>::CONTACT_SPACE;
constexpr FieldSpace ElementsAndOps<3>::CONTACT_SPACE;
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. */
using EntData = EntitiesFieldData::EntData;
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
using SideEle = ElementsAndOps<SPACE_DIM>::FaceSideEle;
//! [Specialisation for assembly]
constexpr FieldSpace CONTACT_SPACE = ElementsAndOps<SPACE_DIM>::CONTACT_SPACE;
//! [Operators used for contact]
using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<AT>::BiLinearForm<
IT>::OpMass<1, SPACE_DIM>;
using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<AT>::LinearForm<
IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
//! [Operators used for contact]
PetscBool is_quasi_static = PETSC_TRUE;
int order = 2; //< Order of displacements in the domain
int contact_order = 2; //< Order of displacements in boundary and side elements
int sigma_order = 1; //< Order of Lagrange multiplier in side elements
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 is_axisymmetric = PETSC_FALSE; //< Axisymmetric model
// #define HENCKY_SMALL_STRAIN
int atom_test = 0;
namespace ContactOps {
double cn_contact = 0.1;
}; // namespace ContactOps
#include <HenckyOps.hpp>
using namespace HenckyOps;
#include <ContactOps.hpp>
#include <PostProcContact.hpp>
#include <ContactNaturalBC.hpp>
#ifdef WITH_MODULE_MFRONT_INTERFACE
#include <MFrontMoFEMInterface.hpp>
#endif
using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
using OpDomainRhsBCs =
DomainRhsBCs::OpFlux<ContactOps::DomainBCs, 1, SPACE_DIM>;
using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
using OpBoundaryRhsBCs =
BoundaryRhsBCs::OpFlux<ContactOps::BoundaryBCs, 1, SPACE_DIM>;
using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;
using OpBoundaryLhsBCs =
BoundaryLhsBCs::OpFlux<ContactOps::BoundaryBCs, 1, SPACE_DIM>;
using namespace ContactOps;
struct Contact {
Contact(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEM::Interface &mField;
MoFEMErrorCode setupProblem();
MoFEMErrorCode createCommonData();
MoFEMErrorCode bC();
MoFEMErrorCode OPs();
MoFEMErrorCode tsSolve();
MoFEMErrorCode checkResults();
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 {
using MoFEM::TimeScale::TimeScale;
double getScale(const double time) {
return scale * MoFEM::TimeScale::getScale(time);
};
};
};
//! [Run problem]
MoFEMErrorCode Contact::runProblem() {
MoFEMFunctionBegin;
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR bC();
CHKERR OPs();
CHKERR tsSolve();
CHKERR checkResults();
MoFEMFunctionReturn(0);
}
//! [Run problem]
//! [Set up problem]
MoFEMErrorCode Contact::setupProblem() {
MoFEMFunctionBegin;
Simple *simple = mField.getInterface<Simple>();
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-contact_order", &contact_order,
PETSC_NULL);
sigma_order = std::max(order, contact_order) - 1;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-sigma_order", &sigma_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) << "Contact order " << contact_order;
MOFEM_LOG("CONTACT", Sev::inform) << "Sigma order " << sigma_order;
MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
// 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);
FieldApproximationBase base;
switch (choice_base_value) {
case AINSWORTH:
base = AINSWORTH_LEGENDRE_BASE;
MOFEM_LOG("CONTACT", Sev::inform)
<< "Set AINSWORTH_LEGENDRE_BASE for displacements";
break;
case DEMKOWICZ:
base = DEMKOWICZ_JACOBI_BASE;
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,
SPACE_DIM);
CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
SPACE_DIM);
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 :
mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % "CONTACT").str()
))
) {
is_contact_block =
true; ///< blocs interation is collective, so that is set irrespective
///< if there are entities 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();
MOFEM_LOG_SYNCHRONISE(mField.get_comm());
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", sigma_order, &boundary_ents);
if (contact_order > order) {
Range ho_ents;
CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false, ho_ents,
moab::Interface::UNION);
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
CHKERR simple->setFieldOrder("U", contact_order, &ho_ents);
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
}
CHKERR simple->setUp();
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();
}
MoFEMFunctionReturn(0);
} //! [Set up problem]
//! [Create common data]
MoFEMErrorCode Contact::createCommonData() {
MoFEMFunctionBegin;
PetscBool use_mfront = PETSC_FALSE;
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);
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 (!mfrontInterface) {
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) {
scale /= young_modulus;
}
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");
#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
if (is_axisymmetric) {
if (SPACE_DIM == 3) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Use executable contact_2d with axisymmetric model");
} else {
if (!use_mfront) {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"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) {
#ifndef WITH_MODULE_MFRONT_INTERFACE
SETERRQ(
PETSC_COMM_SELF, MOFEM_NOT_FOUND,
"MFrontInterface module was not found while use_mfront was set to 1");
#else
if (SPACE_DIM == 3) {
mfrontInterface =
boost::make_shared<MFrontMoFEMInterface<TRIDIMENSIONAL>>(
mField, "U", "GEOMETRY", true, is_quasi_static);
} else if (SPACE_DIM == 2) {
if (is_axisymmetric) {
mfrontInterface =
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();
}
Simple *simple = mField.getInterface<Simple>();
auto dm = simple->getDM();
monitorPtr =
boost::make_shared<Monitor>(dm, scale, mfrontInterface, is_axisymmetric);
if (use_mfront) {
mfrontInterface->setMonitorPtr(monitorPtr);
}
MoFEMFunctionReturn(0);
}
//! [Create common data]
//! [Boundary condition]
MoFEMErrorCode Contact::bC() {
MoFEMFunctionBegin;
auto bc_mng = mField.getInterface<BcManager>();
auto simple = mField.getInterface<Simple>();
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);
// Note remove has to be always before push. Then node marking will be
// corrupted.
CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
MoFEMFunctionReturn(0);
}
//! [Boundary condition]
//! [Push operators to pip]
MoFEMErrorCode Contact::OPs() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
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) {
MoFEMFunctionBegin;
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
"GEOMETRY");
MoFEMFunctionReturn(0);
};
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) {
MoFEMFunctionBegin;
//! [Only used for dynamics]
using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
GAUSS>::OpMass<1, SPACE_DIM>;
//! [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 (!mfrontInterface) {
CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
} else {
CHKERR mfrontInterface->opFactoryDomainLhs(pip);
}
MoFEMFunctionReturn(0);
};
auto add_domain_ops_rhs = [&](auto &pip) {
MoFEMFunctionBegin;
CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
pip, mField, "U", {body_force_time_scale}, Sev::inform);
//! [Only used for dynamics]
using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
//! [Only used for dynamics]
// only in case of dynamics
if (is_quasi_static == PETSC_FALSE) {
auto mat_acceleration = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
"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(
new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
pip.push_back(
new OpInertiaForce("U", mat_velocity, [](double, double, double) {
return alpha_damping * scale;
}));
}
if (!mfrontInterface) {
CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
} else {
CHKERR mfrontInterface->opFactoryDomainRhs(pip);
}
CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
pip, "SIGMA", "U", is_axisymmetric);
MoFEMFunctionReturn(0);
};
auto add_boundary_base_ops = [&](auto &pip) {
MoFEMFunctionBegin;
CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
"GEOMETRY");
// We have to integrate on curved face geometry, thus integration weight
// have to adjusted.
pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
MoFEMFunctionReturn(0);
};
auto add_boundary_ops_lhs = [&](auto &pip) {
MoFEMFunctionBegin;
//! [Operators used for contact]
using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<
AT>::BiLinearForm<IT>::OpMass<1, SPACE_DIM>;
//! [Operators used for contact]
// Add Natural BCs to LHS
CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
pip, mField, "U", Sev::inform);
if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
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) {
return spring_stiffness * scale +
(vis_spring_stiffness * scale) * fe_boundary_lhs->ts_a;
}
));
}
CHKERR
ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
pip, "SIGMA", "U", is_axisymmetric);
CHKERR ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, GAUSS,
DomainEle>(
mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
integration_rule_vol, is_axisymmetric);
MoFEMFunctionReturn(0);
};
auto add_boundary_ops_rhs = [&](auto &pip) {
MoFEMFunctionBegin;
//! [Operators used for contact]
using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<
AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
//! [Operators used for contact]
// Add Natural BCs to RHS
CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
pip, mField, "U", {time_scale}, Sev::inform);
if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
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(
new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", dot_u_disp));
pip.push_back(
new OpSpringRhs("U", u_disp, [this](double, double, double) {
return spring_stiffness * scale;
}));
pip.push_back(
new OpSpringRhs("U", dot_u_disp, [this](double, double, double) {
return vis_spring_stiffness * scale;
}));
}
CHKERR
ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
pip, "SIGMA", "U", is_axisymmetric);
MoFEMFunctionReturn(0);
};
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 (mfrontInterface) {
CHKERR mfrontInterface->setUpdateElementVariablesOperators();
}
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);
MoFEMFunctionReturn(0);
}
//! [Push operators to pip]
//! [Solve]
struct SetUpSchur {
static boost::shared_ptr<SetUpSchur>
createSetUpSchur(MoFEM::Interface &m_field);
virtual MoFEMErrorCode setUp(SmartPetscObj<TS> solver) = 0;
protected:
SetUpSchur() = default;
};
MoFEMErrorCode Contact::tsSolve() {
MoFEMFunctionBegin;
Simple *simple = mField.getInterface<Simple>();
PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
ISManager *is_manager = mField.getInterface<ISManager>();
auto set_section_monitor = [&](auto solver) {
MoFEMFunctionBegin;
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);
MoFEMFunctionReturn(0);
};
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) {
MoFEMFunctionBegin;
monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
boost::shared_ptr<ForcesAndSourcesCore> null;
CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
monitorPtr, null, null);
MoFEMFunctionReturn(0);
};
auto set_essential_bc = [&]() {
MoFEMFunctionBegin;
// 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 = [&]() {
EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
{time_scale}, false);
return hook;
};
pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
auto get_post_proc_hook_rhs = [&]() {
return EssentialPostProcRhs<DisplacementCubitBcData>(
mField, post_proc_rhs_ptr, 1.);
};
auto get_post_proc_hook_lhs = [&]() {
return EssentialPostProcLhs<DisplacementCubitBcData>(
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);
post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
MoFEMFunctionReturn(0);
};
// Set up Schur preconditioner
auto set_schur_pc = [&](auto solver) {
boost::shared_ptr<SetUpSchur> schur_ptr;
if (AT == AssemblyType::BLOCK_SCHUR) {
// Set up Schur preconditioner
schur_ptr = SetUpSchur::createSetUpSchur(mField);
CHK_MOAB_THROW(schur_ptr->setUp(solver), "SetUpSchur::setUp");
}
return schur_ptr;
};
auto dm = simple->getDM();
auto D = createDMVector(dm);
ContactOps::CommonData::createTotalTraction(mField);
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 schur_pc_ptr = set_schur_pc(solver);
auto D = createDMVector(dm);
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 schur_pc_ptr = set_schur_pc(solver);
auto dm = simple->getDM();
auto D = createDMVector(dm);
auto DD = vectorDuplicate(D);
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TS2SetSolution(solver, D, DD);
CHKERR TSSetUp(solver);
CHKERR TSSolve(solver, NULL);
}
MoFEMFunctionReturn(0);
}
//! [Solve]
//! [Check]
MoFEMErrorCode Contact::checkResults() {
MoFEMFunctionBegin;
if (atom_test && !mField.get_comm_rank()) {
const double *t_ptr;
CHKERR VecGetArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
double hertz_force;
double fem_force;
double analytical_active_area = 1.0;
double norm = 1e-5;
double tol_force = 1e-3;
double tol_norm = 7.5; // change when analytical functions are updated
double tol_area = 3e-2;
double fem_active_area = t_ptr[3];
switch (atom_test) {
case 1: // plane stress
hertz_force = 3.927;
fem_force = t_ptr[1];
break;
case 2: // plane strain
hertz_force = 4.675;
fem_force = t_ptr[1];
norm = monitorPtr->getErrorNorm(1);
break;
case 3: // Hertz 3D
hertz_force = 3.968;
tol_force = 2e-3;
fem_force = t_ptr[2];
analytical_active_area = M_PI / 4;
tol_area = 0.2;
break;
case 4: // axisymmetric
tol_force = 5e-3;
tol_area = 0.2;
// analytical_active_area = M_PI;
case 5: // axisymmetric
hertz_force = 15.873;
tol_force = 5e-3;
fem_force = t_ptr[1];
norm = monitorPtr->getErrorNorm(1);
analytical_active_area = M_PI;
break;
case 6: // wavy 2d
hertz_force = 0.374;
fem_force = t_ptr[1];
break;
case 7: // wavy 3d
hertz_force = 0.5289;
fem_force = t_ptr[2];
break;
default:
SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
"atom test %d does not exist", atom_test);
}
if (fabs(fem_force - hertz_force) / hertz_force > tol_force) {
SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test %d failed: Wrong FORCE output: %3.4e != %3.4e",
atom_test, fem_force, hertz_force);
}
if (norm > tol_norm) {
SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test %d failed: Wrong NORM output: %3.4e > %3.4e",
atom_test, norm, tol_norm);
}
if (fabs(fem_active_area - analytical_active_area) > tol_area) {
SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test %d failed: AREA computed %3.4e but should be %3.4e",
atom_test, fem_active_area, analytical_active_area);
}
CHKERR VecRestoreArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
}
ContactOps::CommonData::totalTraction.reset();
MoFEMFunctionReturn(0);
}
//! [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";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// 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";
CHKERR DMRegister_MoFEM(dm_name);
DMType dm_name_mg = "DMMOFEM_MG";
CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
//! [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]
Simple *simple = m_field.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile("");
//! [Load mesh]
//! [CONTACT]
Contact ex(m_field);
CHKERR ex.runProblem();
//! [CONTACT]
}
CATCH_ERRORS;
CHKERR MoFEM::Core::Finalize();
#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() {}
MoFEMErrorCode setUp(SmartPetscObj<TS> solver);
private:
MoFEMErrorCode createSubDM();
MoFEMErrorCode setOperator();
MoFEMErrorCode setPC(PC pc);
MoFEMErrorCode setDiagonalPC(PC pc);
MoFEM::Interface &mField;
SmartPetscObj<Mat> S;
SmartPetscObj<DM> schurDM;
SmartPetscObj<DM> blockDM;
};
MoFEMErrorCode SetUpSchurImpl::setUp(SmartPetscObj<TS> solver) {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
auto pip = mField.getInterface<PipelineManager>();
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 (S) {
CHK_THROW_MESSAGE(
MOFEM_DATA_INCONSISTENCY,
"It is expected that Schur matrix is not allocated. This is "
"possible only if PC is set up twice");
}
CHKERR createSubDM();
// Add data to DM storage
S = createDMMatrix(schurDM);
CHKERR MatSetBlockSize(S, SPACE_DIM);
// CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
// Set DM to use shell block matrix
DM solver_dm;
CHKERR TSGetDM(solver, &solver_dm);
CHKERR DMSetMatType(solver_dm, MATSHELL);
auto ts_ctx_ptr = getDMTsCtx(solver_dm);
auto A = createDMBlockMat(simple->getDM());
auto P = createDMNestSchurMat(simple->getDM());
if (is_quasi_static == PETSC_TRUE) {
auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a,
Mat A, Mat B, void *ctx) {
return TsSetIJacobian(ts, t, u, u_t, a, B, A, ctx);
};
CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
} else {
auto swap_assemble = [](TS ts, PetscReal t, Vec u, Vec u_t, Vec utt,
PetscReal a, PetscReal aa, Mat A, Mat B,
void *ctx) {
return TsSetI2Jacobian(ts, t, u, u_t, utt, a, aa, B, A, ctx);
};
CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
}
CHKERR KSPSetOperators(ksp, A, P);
CHKERR setOperator();
CHKERR setPC(pc);
CHKERR TSSetUp(solver);
CHKERR KSPSetUp(ksp);
CHKERR setDiagonalPC(pc);
} else {
MOFEM_LOG("CONTACT", Sev::inform) << "No Schur PC";
pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
pip->getOpBoundaryLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
pip->getOpDomainLhsPipeline().push_back(createOpSchurAssembleEnd({}, {}));
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode SetUpSchurImpl::createSubDM() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
auto create_dm = [&](const char *name, const char *field_name, auto dm_type) {
auto dm = createDM(mField.get_comm(), dm_type);
auto create_dm_imp = [&]() {
MoFEMFunctionBegin;
CHKERR DMMoFEMCreateSubDM(dm, simple->getDM(), name);
CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
CHKERR DMMoFEMAddElement(dm, simple->getDomainFEName());
CHKERR DMMoFEMAddSubFieldRow(dm, field_name);
CHKERR DMMoFEMAddSubFieldCol(dm, field_name);
CHKERR DMSetUp(dm);
MoFEMFunctionReturn(0);
};
CHK_THROW_MESSAGE(
create_dm_imp(),
"Error in creating schurDM. It is possible that schurDM is "
"already created");
return dm;
};
// Note: here we can make block with bubbles of "U" and "SIGMA" fields. See
// vec-0 where bubbles are added.
schurDM = create_dm("SCHUR", "U", "DMMOFEM_MG");
blockDM = create_dm("BLOCK", "SIGMA", "DMMOFEM");
if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
auto block_mat_data = createBlockMatStructure(
simple->getDM(),
{{
simple->getDomainFEName(),
{
{"U", "U"}, {"SIGMA", "U"}, {"U", "SIGMA"}, {"SIGMA", "SIGMA"}
}}}
);
return createSchurNestedMatrixStruture(
{schur_dm, block_dm}, block_mat_data,
{"SIGMA"}, {nullptr}, true
);
};
auto nested_mat_data = get_nested_mat_data(schurDM, blockDM);
CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
} else {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Only BLOCK_SCHUR is implemented");
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode SetUpSchurImpl::setOperator() {
MoFEMFunctionBegin;
double eps_stab = 1e-4;
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
PETSC_NULL);
using B = FormsIntegrators<BoundaryEleOp>::Assembly<
BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<IT>;
using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
auto simple = mField.getInterface<Simple>();
auto pip = mField.getInterface<PipelineManager>();
auto dm_is = getDMSubData(schurDM)->getSmartRowIs();
auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
// Boundary
pip->getOpBoundaryLhsPipeline().push_front(createOpSchurAssembleBegin());
pip->getOpBoundaryLhsPipeline().push_back(
new OpMassStab("SIGMA", "SIGMA",
[eps_stab](double, double, double) { return eps_stab; }));
pip->getOpBoundaryLhsPipeline().push_back(
createOpSchurAssembleEnd({"SIGMA"}, {nullptr}, ao_up, S, false, false)
);
// Domain
pip->getOpDomainLhsPipeline().push_front(createOpSchurAssembleBegin());
pip->getOpDomainLhsPipeline().push_back(
createOpSchurAssembleEnd({"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]() {
MoFEMFunctionBegin;
CHKERR MatZeroEntries(S);
MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Begin";
MoFEMFunctionReturn(0);
};
post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
post_proc_schur_lhs_ptr]() {
MoFEMFunctionBegin;
MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble End";
auto print_mat_norm = [this](auto a, std::string prefix) {
MoFEMFunctionBegin;
double nrm;
CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
MOFEM_LOG("CONTACT", Sev::noisy) << prefix << " norm = " << nrm;
MoFEMFunctionReturn(0);
};
CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
#ifndef NDEBUG
CHKERR print_mat_norm(S, "S");
#endif // NDEBUG
MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Finish";
MoFEMFunctionReturn(0);
};
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);
MoFEMFunctionReturn(0);
}
MoFEMErrorCode SetUpSchurImpl::setPC(PC pc) {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
auto block_is = getDMSubData(blockDM)->getSmartRowIs();
CHKERR PCFieldSplitSetIS(pc, NULL, block_is);
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
MoFEMFunctionReturn(0);
}
MoFEMErrorCode SetUpSchurImpl::setDiagonalPC(PC pc) {
MoFEMFunctionBegin;
KSP *subksp;
CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
auto get_pc = [](auto ksp) {
PC pc_raw;
CHKERR KSPGetPC(ksp, &pc_raw);
return SmartPetscObj<PC>(pc_raw, true); // bump reference
};
CHKERR setSchurA00MatSolvePC(get_pc(subksp[0]));
auto set_pc_p_mg = [](auto dm, auto pc, auto S) {
MoFEMFunctionBegin;
CHKERR PCSetDM(pc, dm);
PetscBool same = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
if (same) {
CHKERR PCMGSetUpViaApproxOrders(
pc, createPCMGSetUpViaApproxOrdersCtx(dm, S, true));
CHKERR PCSetFromOptions(pc);
}
MoFEMFunctionReturn(0);
};
CHKERR set_pc_p_mg(schurDM, get_pc(subksp[1]), S);
CHKERR PetscFree(subksp);
MoFEMFunctionReturn(0);
}
boost::shared_ptr<SetUpSchur>
SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field) {
return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
}
/**
* \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;
MatrixDouble contactDispGrad;
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
VectorDouble constraintVals;
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, 3, 4};
totalTraction =
createGhostVector(m_field.get_comm(),
(m_field.get_comm_rank() == 0) ? 5 : 0, 5,
(m_field.get_comm_rank() == 0) ? 0 : 5, ghosts);
return totalTraction;
}
static auto getFTensor1TotalTraction() {
if (CommonData::totalTraction) {
const double *t_ptr;
CHK_THROW_MESSAGE(VecGetArrayRead(CommonData::totalTraction, &t_ptr),
"get array");
FTensor::Tensor1<double, 5> t{t_ptr[0], t_ptr[1], t_ptr[2], t_ptr[3],
t_ptr[4]};
CHK_THROW_MESSAGE(VecRestoreArrayRead(CommonData::totalTraction, &t_ptr),
"restore array");
return t;
} else {
return FTensor::Tensor1<double, 5>{0., 0., 0., 0., 0.};
}
}
inline auto contactTractionPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
&contactTraction);
}
inline auto contactDispPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactDisp);
}
inline auto contactDispGradPtr() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
&contactDispGrad);
}
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);
}
inline auto constraintPtr() {
return boost::shared_ptr<VectorDouble>(shared_from_this(), &constraintVals);
}
};
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) {
MoFEMFunctionBegin;
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();
CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
}
MoFEMFunctionReturn(0);
};
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 delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
np::ndarray &sdf
) {
MoFEMFunctionBegin;
try {
// call python function
sdf = bp::extract<np::ndarray>(
sdfFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode evalGradSdf(
double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
np::ndarray &grad_sdf
) {
MoFEMFunctionBegin;
try {
// call python function
grad_sdf = bp::extract<np::ndarray>(
sdfGradFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
}
MoFEMFunctionReturn(0);
}
MoFEMErrorCode evalHessSdf(
double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
np::ndarray &hess_sdf
) {
MoFEMFunctionBegin;
try {
// call python function
hess_sdf = bp::extract<np::ndarray>(
sdfHessFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
} catch (bp::error_already_set const &) {
// print all other errors to stderr
PyErr_Print();
CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
}
MoFEMFunctionReturn(0);
}
private:
bp::object mainNamespace;
bp::object sdfFun;
bp::object sdfGradFun;
bp::object sdfHessFun;
};
static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
inline np::ndarray convert_to_numpy(VectorDouble &data, int nb_gauss_pts,
int id) {
auto dtype = np::dtype::get_builtin<double>();
auto size = bp::make_tuple(nb_gauss_pts);
auto stride = bp::make_tuple(3 * sizeof(double));
return (np::from_data(&data[id], dtype, size, stride, bp::object()));
};
#endif
//! [Surface distance function from python]
using SurfaceDistanceFunction = boost::function<VectorDouble(
double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
MatrixDouble &normals_at_pts, int block_id)>;
using GradSurfaceDistanceFunction = boost::function<MatrixDouble(
double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
MatrixDouble &normals_at_pts, int block_id)>;
using HessSurfaceDistanceFunction = boost::function<MatrixDouble(
double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
MatrixDouble &normals_at_pts, int block_id)>;
inline 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) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
VectorDouble v_spatial_coords = m_spatial_coords.data();
VectorDouble v_normal_at_pts = m_normals_at_pts.data();
bp::list python_coords;
bp::list python_normals;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(
convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
python_normals.append(
convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
}
np::ndarray np_sdf = np::empty(bp::make_tuple(nb_gauss_pts),
np::dtype::get_builtin<double>());
CHK_MOAB_THROW(sdf_ptr->evalSdf(delta_t, t,
bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]),
block_id, np_sdf),
"Failed python call");
double *sdf_val_ptr = reinterpret_cast<double *>(np_sdf.get_data());
VectorDouble v_sdf;
v_sdf.resize(nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
v_sdf[gg] = *(sdf_val_ptr + gg);
return v_sdf;
}
#endif
VectorDouble v_sdf;
v_sdf.resize(nb_gauss_pts, false);
auto t_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
v_sdf[gg] = -t_coords(2) - 0.1;
++t_coords;
}
return v_sdf;
}
inline 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) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
VectorDouble v_spatial_coords = m_spatial_coords.data();
VectorDouble v_normal_at_pts = m_normals_at_pts.data();
bp::list python_coords;
bp::list python_normals;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(
convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
python_normals.append(
convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
}
np::ndarray np_grad_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 3),
np::dtype::get_builtin<double>());
CHK_MOAB_THROW(sdf_ptr->evalGradSdf(
delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]), block_id,
np_grad_sdf),
"Failed python call");
double *grad_ptr = reinterpret_cast<double *>(np_grad_sdf.get_data());
MatrixDouble m_grad_sdf;
m_grad_sdf.resize(3, nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
for (int idx = 0; idx < 3; ++idx)
m_grad_sdf(idx, gg) = *(grad_ptr + (3 * gg + idx));
}
return m_grad_sdf;
}
#endif
MatrixDouble m_grad_sdf;
m_grad_sdf.resize(3, nb_gauss_pts, false);
FTensor::Index<'i', 3> i;
FTensor::Tensor1<double, 3> t_grad_sdf_set{0.0, 0.0, -1.0};
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
t_grad_sdf(i) = t_grad_sdf_set(i);
++t_grad_sdf;
}
return m_grad_sdf;
}
inline 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) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
VectorDouble v_spatial_coords = m_spatial_coords.data();
VectorDouble v_normal_at_pts = m_normals_at_pts.data();
bp::list python_coords;
bp::list python_normals;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(
convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
python_normals.append(
convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
};
np::ndarray np_hess_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 6),
np::dtype::get_builtin<double>());
CHK_MOAB_THROW(sdf_ptr->evalHessSdf(
delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]), block_id,
np_hess_sdf),
"Failed python call");
double *hess_ptr = reinterpret_cast<double *>(np_hess_sdf.get_data());
MatrixDouble m_hess_sdf;
m_hess_sdf.resize(6, nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
for (int idx = 0; idx < 6; ++idx)
m_hess_sdf(idx, gg) =
*(hess_ptr + (6 * gg + idx));
}
return m_hess_sdf;
}
#endif
MatrixDouble m_hess_sdf;
m_hess_sdf.resize(6, nb_gauss_pts, false);
FTensor::Index<'i', 3> i;
FTensor::Index<'j', 3> j;
FTensor::Tensor2_symmetric<double, 3> t_hess_sdf_set{0., 0., 0., 0., 0., 0.};
auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
t_hess_sdf(i, j) = t_hess_sdf_set(i, j);
++t_hess_sdf;
}
return m_hess_sdf;
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl;
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpAssembleTotalContactAreaImpl;
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 <typename T1, typename T2, int DIM1, int DIM2>
inline auto get_spatial_coords(FTensor::Tensor1<T1, DIM1> &&t_coords,
FTensor::Tensor1<T2, DIM2> &&t_disp,
size_t nb_gauss_pts) {
MatrixDouble m_spatial_coords(nb_gauss_pts, 3);
m_spatial_coords.clear();
auto t_spatial_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
FTensor::Index<'i', DIM2> i;
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_spatial_coords(i) = t_coords(i) + t_disp(i);
++t_spatial_coords;
++t_coords;
++t_disp;
}
return m_spatial_coords;
}
template <typename T1, int DIM1>
inline auto get_normalize_normals(FTensor::Tensor1<T1, DIM1> &&t_normal_at_pts,
size_t nb_gauss_pts) {
MatrixDouble m_normals_at_pts(3, nb_gauss_pts);
m_normals_at_pts.clear();
FTensor::Index<'i', DIM1> i;
auto t_set_normal = getFTensor1FromMat<3>(m_normals_at_pts);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_set_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
++t_set_normal;
++t_normal_at_pts;
}
return m_normals_at_pts;
}
template <int DIM, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>
: public BoundaryEleOp {
OpAssembleTotalContactTractionImpl(
boost::shared_ptr<CommonData> common_data_ptr, double scale = 1,
bool is_axisymmetric = false);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
const double scaleTraction;
bool isAxisymmetric;
};
template <int DIM, typename BoundaryEleOp>
struct OpAssembleTotalContactAreaImpl<DIM, GAUSS, BoundaryEleOp>
: public BoundaryEleOp {
OpAssembleTotalContactAreaImpl(
boost::shared_ptr<CommonData> common_data_ptr,
bool is_axisymmetric = false,
boost::shared_ptr<Range> contact_range_ptr = nullptr);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
grad_surface_distance_function;
private:
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
boost::shared_ptr<Range> contactRange;
};
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;
SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
grad_surface_distance_function;
HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
hess_surface_distance_function;
};
template <int DIM, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
: public AssemblyBoundaryEleOp {
OpConstrainBoundaryRhsImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr,
bool is_axisymmetric = false);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
grad_surface_distance_function;
private:
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
};
template <int DIM, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
: public AssemblyBoundaryEleOp {
OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
bool is_axisymmetric = false);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
grad_surface_distance_function;
HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
hess_surface_distance_function;
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
};
template <int DIM, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
: public AssemblyBoundaryEleOp {
OpConstrainBoundaryLhs_dTractionImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
bool is_axisymmetric = false);
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data);
SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
grad_surface_distance_function;
private:
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
};
template <typename BoundaryEleOp> struct ContactIntegrators {
template <int DIM, IntegrationType I>
using OpAssembleTotalContactTraction =
OpAssembleTotalContactTractionImpl<DIM, I, BoundaryEleOp>;
template <int DIM, IntegrationType I>
using OpAssembleTotalContactArea =
OpAssembleTotalContactAreaImpl<DIM, I, BoundaryEleOp>;
template <int DIM, IntegrationType I>
using OpEvaluateSDF = OpEvaluateSDFImpl<DIM, I, BoundaryEleOp>;
template <AssemblyType A> struct Assembly {
using AssemblyBoundaryEleOp =
typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
template <int DIM, IntegrationType I>
using OpConstrainBoundaryRhs =
OpConstrainBoundaryRhsImpl<DIM, I, AssemblyBoundaryEleOp>;
template <int DIM, IntegrationType I>
using OpConstrainBoundaryLhs_dU =
OpConstrainBoundaryLhs_dUImpl<DIM, I, AssemblyBoundaryEleOp>;
template <int DIM, IntegrationType I>
using OpConstrainBoundaryLhs_dTraction =
OpConstrainBoundaryLhs_dTractionImpl<DIM, I, AssemblyBoundaryEleOp>;
};
};
inline double sign(double x) {
constexpr auto eps = std::numeric_limits<float>::epsilon();
if (std::abs(x) < eps)
return 0;
else if (x > eps)
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 sdf 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,
bool is_axisymmetric)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
commonDataPtr(common_data_ptr), scaleTraction(scale),
isAxisymmetric(is_axisymmetric) {}
template <int DIM, typename BoundaryEleOp>
MoFEMErrorCode
OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * BoundaryEleOp::getMeasure();
t_sum_t(i) += alpha * t_traction(i);
++t_w;
++t_traction;
++t_coords;
}
t_sum_t(i) *= scaleTraction;
constexpr int ind[] = {0, 1, 2};
CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
ADD_VALUES);
MoFEMFunctionReturn(0);
}
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactAreaImpl<DIM, GAUSS, BoundaryEleOp>::
OpAssembleTotalContactAreaImpl(
boost::shared_ptr<CommonData> common_data_ptr, bool is_axisymmetric,
boost::shared_ptr<Range> contact_range_ptr)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric),
contactRange(contact_range_ptr) {}
template <int DIM, typename BoundaryEleOp>
MoFEMErrorCode
OpAssembleTotalContactAreaImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
auto fe_type = BoundaryEleOp::getFEType();
const auto fe_ent = BoundaryEleOp::getFEEntityHandle();
if (contactRange->find(fe_ent) != contactRange->end()) {
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Tensor1<double, 2> t_sum_a{0., 0.};
auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_grad = getFTensor2FromMat<DIM, DIM>(commonDataPtr->contactDispGrad);
auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
auto m_spatial_coords = get_spatial_coords(
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
auto m_normals_at_pts = get_normalize_normals(
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = BoundaryEleOp::getTStime();
auto ts_time_step = BoundaryEleOp::getTStimeStep();
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf = gradSurfaceDistanceFunction(
ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0); // Axisymmetric Jacobian
}
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(t_sdf, tn);
double alpha = t_w * jacobian;
FTensor::Tensor2<double, DIM, DIM> F;
FTensor::Tensor2<double, DIM, DIM> invF;
FTensor::Tensor1<double, DIM> t_normal_current;
F(i, j) = t_grad(i, j) + kronecker_delta(i, j);
auto det = determinantTensor(F);
CHKERR invertTensor(F, det, invF);
t_normal_current(i) = det * (invF(j, i) * t_normal_at_pts(j));
alpha *= sqrt(t_normal_current(i) * t_normal_current(i));
if (fe_type == MBTRI) {
alpha /= 2;
}
if (c > 1e-12) {
t_sum_a(0) += alpha; // real area
}
t_sum_a(1) += alpha; // Potential area
++t_w;
++t_traction;
++t_coords;
++t_sdf;
++t_grad_sdf;
++t_grad;
++t_normal_at_pts;
}
constexpr int ind[] = {3, 4};
CHKERR VecSetValues(commonDataPtr->totalTraction, 2, ind, &t_sum_a(0),
ADD_VALUES);
}
MoFEMFunctionReturn(0);
}
template <int DIM, typename BoundaryEleOp>
OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::OpEvaluateSDFImpl(
boost::shared_ptr<CommonData> common_data_ptr)
: BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
commonDataPtr(common_data_ptr) {}
template <int DIM, typename BoundaryEleOp>
MoFEMErrorCode
OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::doWork(int side, EntityType type,
EntData &data) {
MoFEMFunctionBegin;
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
auto &sdf_vec = commonDataPtr->sdfVals;
auto &grad_mat = commonDataPtr->gradsSdf;
auto &hess_mat = commonDataPtr->hessSdf;
auto &constraint_vec = commonDataPtr->constraintVals;
auto &contactTraction_mat = commonDataPtr->contactTraction;
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);
constraint_vec.resize(nb_gauss_pts, false);
auto t_traction = getFTensor1FromMat<DIM>(contactTraction_mat);
auto t_sdf = getFTensor0FromVec(sdf_vec);
auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
auto t_constraint = getFTensor0FromVec(constraint_vec);
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
auto ts_time = BoundaryEleOp::getTStime();
auto ts_time_step = BoundaryEleOp::getTStimeStep();
auto m_spatial_coords = get_spatial_coords(
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
auto m_normals_at_pts = get_normalize_normals(
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
// placeholder to pass boundary block id to python
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_hess_sdf =
hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf_v = getFTensor0FromVec(v_sdf);
auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
auto next = [&]() {
++t_sdf;
++t_sdf_v;
++t_grad_sdf;
++t_grad_sdf_v;
++t_hess_sdf;
++t_hess_sdf_v;
++t_disp;
++t_traction;
++t_constraint;
};
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
auto tn = -t_traction(i) * t_grad_sdf_v(i);
auto c = constrain(t_sdf_v, tn);
t_sdf = t_sdf_v;
t_grad_sdf(i) = t_grad_sdf_v(i);
t_hess_sdf(i, j) = t_hess_sdf_v(i, j);
t_constraint = c;
next();
}
MoFEMFunctionReturn(0);
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
OpConstrainBoundaryRhsImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data_ptr,
bool is_axisymmetric)
: AssemblyBoundaryEleOp(field_name, field_name,
AssemblyBoundaryEleOp::OPROW),
commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {}
template <int DIM, typename AssemblyBoundaryEleOp>
MoFEMErrorCode
OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
EntitiesFieldData::EntData &data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
FTensor::Index<'l', DIM> l;
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &nf = AssemblyBoundaryEleOp::locF;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
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>();
auto m_spatial_coords = get_spatial_coords(
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
auto m_normals_at_pts = get_normalize_normals(
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
// placeholder to pass boundary block id to python
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(t_sdf, tn);
FTensor::Tensor2<double, DIM, DIM> t_cP;
t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
FTensor::Tensor2<double, DIM, DIM> t_cQ;
t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
FTensor::Tensor1<double, DIM> t_rhs;
t_rhs(i) =
t_cQ(i, j) * (t_disp(j) - cn_contact * t_traction(j))
+
t_cP(i, j) * t_disp(j) +
c * (t_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;
++t_sdf;
++t_grad_sdf;
}
MoFEMFunctionReturn(0);
}
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,
bool is_axisymmetric)
: AssemblyBoundaryEleOp(row_field_name, col_field_name,
AssemblyBoundaryEleOp::OPROWCOL),
commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
AssemblyBoundaryEleOp::sYmm = false;
}
template <int DIM, typename AssemblyBoundaryEleOp>
MoFEMErrorCode
OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
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>();
auto m_spatial_coords = get_spatial_coords(
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
auto m_normals_at_pts = get_normalize_normals(
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
// placeholder to pass boundary block id to python
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_hess_sdf =
hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(t_sdf, tn);
FTensor::Tensor2<double, DIM, DIM> t_cP;
t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
FTensor::Tensor2<double, DIM, DIM> t_cQ;
t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
FTensor::Tensor2<double, DIM, DIM> t_res_dU;
t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
if (c > 0) {
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 * t_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_traction;
++t_coords;
++t_w;
++t_normal;
++t_sdf;
++t_grad_sdf;
++t_hess_sdf;
}
MoFEMFunctionReturn(0);
}
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, bool is_axisymmetric)
: AssemblyBoundaryEleOp(row_field_name, col_field_name,
AssemblyBoundaryEleOp::OPROWCOL),
commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
AssemblyBoundaryEleOp::sYmm = false;
}
template <int DIM, typename AssemblyBoundaryEleOp>
MoFEMErrorCode
OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
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;
auto m_spatial_coords = get_spatial_coords(
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
auto m_normals_at_pts = get_normalize_normals(
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
// placeholder to pass boundary block id to python
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
auto tn = -t_traction(i) * t_grad_sdf(i);
auto c = constrain(t_sdf, tn);
FTensor::Tensor2<double, DIM, DIM> t_cP;
t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
FTensor::Tensor2<double, DIM, DIM> t_cQ;
t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
FTensor::Tensor2<double, DIM, DIM> t_res_dt;
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_traction;
++t_coords;
++t_w;
++t_normal;
++t_sdf;
++t_grad_sdf;
}
MoFEMFunctionReturn(0);
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode opFactoryDomainRhs(
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u, bool is_axisymmetric = false) {
MoFEMFunctionBegin;
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template LinearForm<I>;
using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
using OpMixDivUCylRhs =
typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
using OpMixUTimesDivLambdaRhs =
typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
using OpMixUTimesLambdaRhs =
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>();
auto jacobian = [is_axisymmetric](const double r, const double,
const double) {
if (is_axisymmetric)
return 2. * M_PI * r;
else
return 1.;
};
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(
new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
if (!is_axisymmetric) {
pip.push_back(
new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
} else {
pip.push_back(new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
sigma, div_stress_ptr));
}
pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
if (!is_axisymmetric) {
pip.push_back(
new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
} else {
pip.push_back(new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
jacobian));
}
pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr, jacobian));
pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr, jacobian));
MoFEMFunctionReturn(0);
}
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) {
MoFEMFunctionBegin;
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);
}
MoFEMFunctionReturn(0);
}
};
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
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) {
MoFEMFunctionBegin;
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 OpMixDivUCylLhs =
typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
using OpMixDivUCylLhsSide = OpMixLhsSide<OpMixDivUCylLhs>;
using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
auto unity = []() { return 1; };
auto jacobian = [is_axisymmetric](const double r, const double,
const double) {
if (is_axisymmetric)
return 2. * M_PI * r;
else
return 1.;
};
if (!is_axisymmetric) {
op_loop_side->getOpPtrVector().push_back(
new OpMixDivULhsSide(sigma, u, unity, jacobian, true));
} else {
op_loop_side->getOpPtrVector().push_back(
new OpMixDivUCylLhsSide(sigma, u, unity, jacobian, true));
}
op_loop_side->getOpPtrVector().push_back(
new OpLambdaGraULhsSide(sigma, u, unity, jacobian, true));
op_loop_side->getSideFEPtr()->getRuleHook = rule;
MoFEMFunctionReturn(0);
}
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
MoFEMErrorCode opFactoryBoundaryLhs(
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u, bool is_axisymmetric = false) {
MoFEMFunctionBegin;
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, is_axisymmetric));
pip.push_back(new typename C::template Assembly<A>::
template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
sigma, sigma, common_data_ptr, is_axisymmetric));
MoFEMFunctionReturn(0);
}
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
MoFEMErrorCode opFactoryBoundaryRhs(
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, std::string u, bool is_axisymmetric = false) {
MoFEMFunctionBegin;
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, is_axisymmetric));
MoFEMFunctionReturn(0);
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
MoFEMErrorCode opFactoryCalculateTraction(
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string sigma, bool is_axisymmetric = false) {
MoFEMFunctionBegin;
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, is_axisymmetric));
MoFEMFunctionReturn(0);
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
MoFEMErrorCode opFactoryCalculateArea(
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
OpLoopSide<SideEle> *op_loop_side, std::string sigma, std::string u,
bool is_axisymmetric = false,
boost::shared_ptr<Range> contact_range_ptr = nullptr) {
MoFEMFunctionBegin;
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
op_loop_side->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
"U", common_data_ptr->contactDispGradPtr()));
if (contact_range_ptr) {
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(op_loop_side);
pip.push_back(new typename C::template OpAssembleTotalContactArea<DIM, I>(
common_data_ptr, is_axisymmetric, contact_range_ptr));
}
MoFEMFunctionReturn(0);
}
}; // namespace ContactOps
#endif // __CONTACTOPS_HPP__
ContactOps::get_spatial_coords
auto get_spatial_coords(FTensor::Tensor1< T1, DIM1 > &&t_coords, FTensor::Tensor1< T2, DIM2 > &&t_disp, size_t nb_gauss_pts)
Definition: ContactOps.hpp:414
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
ContactOps::CommonData::contactTractionPtr
auto contactTractionPtr()
Definition: ContactOps.hpp:57
MoFEM::createAOMappingIS
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
Definition: PetscSmartObj.hpp:318
OpMixLhs
SetUpSchurImpl
Definition: test_broken_space.cpp:515
ContactOps::sign
double sign(double x)
Definition: ContactOps.hpp:580
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
ContactOps::CommonData::constraintVals
VectorDouble constraintVals
Definition: ContactOps.hpp:26
ContactOps::CommonData::contactDisp
MatrixDouble contactDisp
Definition: ContactOps.hpp:17
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
sdf.hess_sdf
def hess_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:19
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
spring_stiffness
double spring_stiffness
Definition: contact.cpp:87
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
Contact::bC
MoFEMErrorCode bC()
[Create common data]
Definition: contact.cpp:435
LASTBASE
@ LASTBASE
Definition: definitions.h:69
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
ContactOps::hess_surface_distance_function
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:338
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
ContactOps::w
double w(const double sdf, const double tn)
Definition: ContactOps.hpp:590
ContactOps
Definition: contact.cpp:99
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
Definition: test_broken_space.cpp:528
rho
double rho
Definition: plastic.cpp:140
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
ContactOps::constrain
double constrain(double sdf, double tn)
constrain function
Definition: ContactOps.hpp:603
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
is_axisymmetric
PetscBool is_axisymmetric
Definition: contact.cpp:93
PlasticOps::opFactoryDomainRhs
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
Definition: PlasticOps.hpp:500
MoFEM::BLOCK_PRECONDITIONER_SCHUR
@ BLOCK_PRECONDITIONER_SCHUR
Definition: FormsIntegrators.hpp:109
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:134
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::NaturalBC::Assembly::BiLinearForm
Definition: Natural.hpp:74
MoFEM::createPCMGSetUpViaApproxOrdersCtx
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
Definition: PCMGSetUpViaApproxOrders.cpp:630
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
ContactOps::CommonData::contactDispGrad
MatrixDouble contactDispGrad
Definition: ContactOps.hpp:18
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
ContactOps::ContactIntegrators::Assembly::OpConstrainBoundaryLhs_dU
OpConstrainBoundaryLhs_dUImpl< DIM, I, AssemblyBoundaryEleOp > OpConstrainBoundaryLhs_dU
Definition: ContactOps.hpp:572
sdf
Definition: sdf.py:1
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
ContactOps::CommonData::contactDispGradPtr
auto contactDispGradPtr()
Definition: ContactOps.hpp:66
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
SetUpSchurImpl::setPC
MoFEMErrorCode setPC(PC pc)
Definition: contact.cpp:1226
ContactOps::SurfaceDistanceFunction
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:213
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:113
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::TsSetI2Jacobian
PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, PetscReal a, PetscReal aa, Mat A, Mat B, void *ctx)
Calculation Jacobian for second order PDE in time.
Definition: TsCtx.cpp:511
ContactOps::HessSurfaceDistanceFunction
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:221
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
ContactOps::ContactIntegrators::OpAssembleTotalContactArea
OpAssembleTotalContactAreaImpl< DIM, I, BoundaryEleOp > OpAssembleTotalContactArea
Definition: ContactOps.hpp:556
HenckyOps
Definition: HenckyOps.hpp:11
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
scale
double scale
Definition: plastic.cpp:119
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
HenckyOps.hpp
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
ContactOps::opFactoryBoundaryToDomainLhs
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)
Definition: ContactOps.hpp:1265
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1157
geom_order
int geom_order
Order if fixed.
Definition: plastic.cpp:137
SetUpSchurImpl::createSubDM
MoFEMErrorCode createSubDM()
Definition: contact.cpp:1083
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
sdf.r
int r
Definition: sdf.py:8
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1056
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::createDMBlockMat
auto createDMBlockMat(DM dm)
Definition: DMMoFEM.hpp:1076
SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE
Definition: contact.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: dynamic_first_order_con_law.cpp:63
atom_test
int atom_test
Definition: contact.cpp:97
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
ContactOps::ContactIntegrators::OpAssembleTotalContactTraction
OpAssembleTotalContactTractionImpl< DIM, I, BoundaryEleOp > OpAssembleTotalContactTraction
Definition: ContactOps.hpp:552
MoFEM::createOpSchurAssembleEnd
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range >> field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
Definition: Schur.cpp:2186
AT
constexpr AssemblyType AT
Definition: test_broken_space.cpp:20
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:23
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::createGhostVector
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
Definition: PetscSmartObj.hpp:179
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
vis_spring_stiffness
double vis_spring_stiffness
Definition: contact.cpp:88
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
Definition: ScalingMethod.cpp:22
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
BoundaryEleOp
ContactOps::CommonData::totalTraction
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:29
OpSpringRhs
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::LinearForm< IT >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:75
is_quasi_static
PetscBool is_quasi_static
Definition: plastic.cpp:139
ContactOps::ContactIntegrators::Assembly::OpConstrainBoundaryLhs_dTraction
OpConstrainBoundaryLhs_dTractionImpl< DIM, I, AssemblyBoundaryEleOp > OpConstrainBoundaryLhs_dTraction
Definition: ContactOps.hpp:576
ContactOps::get_normalize_normals
auto get_normalize_normals(FTensor::Tensor1< T1, DIM1 > &&t_normal_at_pts, size_t nb_gauss_pts)
Definition: ContactOps.hpp:431
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::DMMoFEMCreateSubDM
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:215
double
ContactOps::opFactoryBoundaryRhs
MoFEMErrorCode opFactoryBoundaryRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1343
convert.type
type
Definition: convert.py:64
ContactOps::opFactoryCalculateArea
MoFEMErrorCode opFactoryCalculateArea(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, OpLoopSide< SideEle > *op_loop_side, std::string sigma, std::string u, bool is_axisymmetric=false, boost::shared_ptr< Range > contact_range_ptr=nullptr)
Definition: ContactOps.hpp:1381
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
MoFEM::DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:302
MoFEM::TimeScale
Force scale operator for reading two columns.
Definition: ScalingMethod.hpp:32
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::invertTensor
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
Definition: Templates.hpp:1771
Contact
Definition: contact.cpp:125
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:197
MoFEM::determinantTensor
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
Definition: Templates.hpp:1571
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
ContactOps::CommonData::hessSdf
MatrixDouble hessSdf
Definition: ContactOps.hpp:23
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2181
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
ContactOps::CommonData::gradSdfPtr
auto gradSdfPtr()
Definition: ContactOps.hpp:75
contact_order
int contact_order
Definition: contact.cpp:81
MatrixFunction.hpp
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
MoFEM::createDMNestSchurMat
auto createDMNestSchurMat(DM dm)
Definition: DMMoFEM.hpp:1083
MoFEM::TsSetIJacobian
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:165
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
SetUpSchurImpl::setOperator
MoFEMErrorCode setOperator()
Definition: contact.cpp:1150
ContactOps::CommonData::sdfVals
VectorDouble sdfVals
size is equal to number of gauss points on element
Definition: ContactOps.hpp:20
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
Contact::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: contact.cpp:314
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
ContactOps::CommonData::gradsSdf
MatrixDouble gradsSdf
Definition: ContactOps.hpp:21
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
ContactOps::CommonData::sdfPtr
auto sdfPtr()
Definition: ContactOps.hpp:71
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
MoFEM::OpSetHOWeightsOnSubDim
Definition: HODataOperators.hpp:145
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
ContactOps::CommonData::contactDispPtr
auto contactDispPtr()
Definition: ContactOps.hpp:62
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
MoFEM::TimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: ScalingMethod.cpp:133
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
ContactOps.hpp
FTensor::Index< 'i', 3 >
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
ContactOps::CommonData::createTotalTraction
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:31
ContactOps::cn_contact
double cn_contact
Definition: contact.cpp:100
MoFEM::setSchurA00MatSolvePC
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
Definition: Schur.cpp:2223
ContactOps::GradSurfaceDistanceFunction
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:217
Contact::OPs
MoFEMErrorCode OPs()
[Boundary condition]
Definition: contact.cpp:472
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ElementsAndOps
Definition: child_and_parent.cpp:18
Range
DomainEleOp
ContactOps::opFactoryBoundaryLhs
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1319
MoFEM::CoreTmp< 0 >::Initialize
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
Contact::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: contact.cpp:174
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
ContactOps::CommonData::constraintPtr
auto constraintPtr()
Definition: ContactOps.hpp:83
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
ContactOps::CommonData::contactTraction
MatrixDouble contactTraction
Definition: ContactOps.hpp:16
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
FTensor::kronecker_delta
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
Definition: Kronecker_Delta.hpp:81
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
CONTACT_SPACE
constexpr FieldSpace CONTACT_SPACE
Definition: plastic.cpp:52
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
approx_order
int approx_order
Definition: test_broken_space.cpp:50
MoFEM::createSchurNestedMatrixStruture
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM >> dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range >> field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition: Schur.cpp:1944
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
MoFEM::DMMoFEMTSSetMonitor
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:1056
IT
constexpr IntegrationType IT
Definition: test_broken_space.cpp:24
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
ContactOps::CommonData::getFTensor1TotalTraction
static auto getFTensor1TotalTraction()
Definition: ContactOps.hpp:42
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
Contact::tsSolve
MoFEMErrorCode tsSolve()
Definition: contact.cpp:702
ContactOps::ContactIntegrators::Assembly::OpConstrainBoundaryRhs
OpConstrainBoundaryRhsImpl< DIM, I, AssemblyBoundaryEleOp > OpConstrainBoundaryRhs
Definition: ContactOps.hpp:568
ContactOps::surface_distance_function
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:223
Contact::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: contact.cpp:161
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Definition: FormsIntegrators.hpp:108
ContactNaturalBC.hpp
SetUpSchur
[Push operators to pipeline]
Definition: test_broken_space.cpp:40
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
AssemblyBoundaryEleOp
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::DMMoFEMSetNestSchurData
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition: DMMoFEM.cpp:1562
MoFEM::createBlockMatStructure
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition: Schur.cpp:1009
MoFEM::PCMGSetUpViaApproxOrders
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
Definition: PCMGSetUpViaApproxOrders.cpp:634
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
sdf_wavy_2d.ind
float ind
Definition: sdf_wavy_2d.py:7
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:578
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< VecScatter >
SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
Definition: test_broken_space.cpp:768
OpSpringLhs
FormsIntegrators< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT >::OpMass< 1, SPACE_DIM > OpSpringLhs
[Operators used for contact]
Definition: contact.cpp:73
ContactOps::grad_surface_distance_function
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:281
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
ContactOps::opFactoryCalculateTraction
MoFEMErrorCode opFactoryCalculateTraction(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1364
sigma_order
int sigma_order
Definition: contact.cpp:82
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
SetUpSchurImpl::setDiagonalPC
MoFEMErrorCode setDiagonalPC(PC pc)
Definition: contact.cpp:1235
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
PostProcContact.hpp
ContactOps::ContactIntegrators::OpEvaluateSDF
OpEvaluateSDFImpl< DIM, I, BoundaryEleOp > OpEvaluateSDF
Definition: ContactOps.hpp:559
ContactOps::CommonData::hessSdfPtr
auto hessSdfPtr()
Definition: ContactOps.hpp:79
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
sdf.grad_sdf
def grad_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:15
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
Contact::checkResults
MoFEMErrorCode checkResults()
[Solve]
Definition: contact.cpp:839
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
alpha_damping
double alpha_damping
Definition: plastic.cpp:141
F
@ F
Definition: free_surface.cpp:394
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow