10#ifndef EXECUTABLE_DIMENSION
11#define EXECUTABLE_DIMENSION 3
18#include <boost/python.hpp>
19#include <boost/python/def.hpp>
20namespace bp = boost::python;
27 IntegrationType::GAUSS;
62 GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
91 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
109 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
110 std::string
field_name, std::string block_name,
111 boost::shared_ptr<MatrixDouble> mat_D_Ptr,
Sev sev);
161 boost::shared_ptr<SDFPython> sdfPythonPtr;
191 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
192 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz"};
193 PetscInt choice_base_value = AINSWORTH;
195 LASBASETOPT, &choice_base_value, PETSC_NULL);
198 switch (choice_base_value) {
202 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
207 <<
"Set DEMKOWICZ_JACOBI_BASE for displacements";
229 auto get_skin = [&]() {
234 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
238 auto filter_blocks = [&](
auto skin) {
243 (boost::format(
"%s(.*)") %
"CONTACT").str()
249 <<
"Find contact block set: " <<
m->getName();
250 auto meshset =
m->getMeshset();
252 contact_range,
true);
255 <<
"Nb entities in contact surface: " << contact_range.size();
259 skin = intersect(skin, contact_range);
264 auto filter_true_skin = [&](
auto skin) {
266 ParallelComm *pcomm =
268 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
269 PSTATUS_NOT, -1, &boundary_ents);
270 return boundary_ents;
273 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
293 auto get_options = [&]() {
320 sdfPythonPtr = boost::make_shared<SDFPython>();
321 CHKERR sdfPythonPtr->sdfInit(
"sdf.py");
322 sdfPythonWeakPtr = sdfPythonPtr;
335 for (
auto f : {
"U",
"SIGMA"}) {
336 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
337 "REMOVE_X",
f, 0, 0);
338 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
339 "REMOVE_Y",
f, 1, 1);
340 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
341 "REMOVE_Z",
f, 2, 2);
342 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
343 "REMOVE_ALL",
f, 0, 3);
346 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_X",
347 "SIGMA", 0, 0,
false,
true);
348 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_Y",
349 "SIGMA", 1, 1,
false,
true);
350 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_Z",
351 "SIGMA", 2, 2,
false,
true);
352 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_ALL",
353 "SIGMA", 0, 3,
false,
true);
354 CHKERR bc_mng->removeBlockDOFsOnEntities(
355 simple->getProblemName(),
"NO_CONTACT",
"SIGMA", 0, 3,
false,
true);
360 simple->getProblemName(),
"U");
361 boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{
"FIX_",
"ROTATE_"});
373 auto time_scale = boost::make_shared<TimeScale>();
375 auto add_domain_base_ops = [&](
auto &pip) {
380 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
381 auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
382 henky_common_data_ptr->matGradPtr = common_data_ptr->mGradPtr();
383 henky_common_data_ptr->matDPtr = common_data_ptr->mDPtr();
385 auto add_domain_ops_lhs = [&](
auto &pip) {
389 common_data_ptr->mDPtr(), Sev::verbose);
392 "U", common_data_ptr->mGradPtr()));
404 new OpKPiola(
"U",
"U", henky_common_data_ptr->getMatTangent()));
407 auto get_inertia_and_mass_dumping = [
this](
const double,
const double,
410 auto &fe_domain_lhs = pip_mng->getDomainLhsFE();
413 pip.push_back(
new OpMass(
"U",
"U", get_inertia_and_mass_dumping));
415 auto get_mass_dumping = [
this](
const double,
const double,
418 auto &fe_domain_lhs = pip_mng->getDomainLhsFE();
421 pip.push_back(
new OpMass(
"U",
"U", get_mass_dumping));
424 auto unity = []() {
return 1; };
425 pip.push_back(
new OpMixDivULhs(
"SIGMA",
"U", unity,
true));
431 auto add_domain_ops_rhs = [&](
auto &pip) {
435 pip,
mField,
"U", {time_scale}, Sev::inform);
438 common_data_ptr->mDPtr(), Sev::inform);
440 "U", common_data_ptr->mGradPtr()));
452 "U", henky_common_data_ptr->getMatFirstPiolaStress()));
455 "U", common_data_ptr->contactDispPtr()));
458 "SIGMA", common_data_ptr->contactStressPtr()));
460 "SIGMA", common_data_ptr->contactStressDivergencePtr()));
462 pip.push_back(
new OpMixDivURhs(
"SIGMA", common_data_ptr->contactDispPtr(),
463 [](
double,
double,
double) { return 1; }));
468 "U", common_data_ptr->contactStressDivergencePtr()));
474 auto mat_acceleration = boost::make_shared<MatrixDouble>();
476 "U", mat_acceleration));
478 "U", mat_acceleration, [](
double,
double,
double) {
return rho; }));
481 auto mat_velocity = boost::make_shared<MatrixDouble>();
485 new OpInertiaForce(
"U", mat_velocity, [](
double,
double,
double) {
492 auto add_boundary_base_ops = [&](
auto &pip) {
496 "U", common_data_ptr->contactDispPtr()));
498 "SIGMA", common_data_ptr->contactTractionPtr()));
501 auto add_boundary_ops_lhs = [&](
auto &pip) {
508 pip,
mField,
"U", Sev::inform);
525 auto add_boundary_ops_rhs = [&](
auto &pip) {
528 auto u_mat_ptr = boost::make_shared<MatrixDouble>();
533 {boost::make_shared<TimeScale>()});
538 pip,
mField,
"U", {time_scale}, Sev::inform);
542 "U", common_data_ptr->contactDispPtr(),
543 [
this](
double,
double,
double) { return spring_stiffness; }));
548 add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
549 add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
550 add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
551 add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
553 add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
554 add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
555 CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
556 CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
558 auto integration_rule_vol = [](int, int,
int approx_order) {
561 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
562 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
563 auto integration_rule_boundary = [](int, int,
int approx_order) {
566 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
567 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
581 auto set_section_monitor = [&](
auto solver) {
584 CHKERR TSGetSNES(solver, &snes);
585 PetscViewerAndFormat *vf;
586 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
587 PETSC_VIEWER_DEFAULT, &vf);
590 (
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
595 auto scatter_create = [&](
auto D,
auto coeff) {
597 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
598 ROW,
"U", coeff, coeff, is);
600 CHKERR ISGetLocalSize(is, &loc_size);
604 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
609 auto set_time_monitor = [&](
auto dm,
auto solver) {
611 boost::shared_ptr<Monitor> monitor_ptr(
613 boost::shared_ptr<ForcesAndSourcesCore> null;
615 monitor_ptr, null, null);
619 auto set_fieldsplit_preconditioner = [&](
auto solver) {
623 CHKERR TSGetSNES(solver, &snes);
625 CHKERR SNESGetKSP(snes, &ksp);
627 CHKERR KSPGetPC(ksp, &pc);
628 PetscBool is_pcfs = PETSC_FALSE;
629 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
632 if (is_pcfs == PETSC_TRUE) {
634 auto name_prb =
simple->getProblemName();
635 auto is_all_bc = bc_mng->getBlockIS(name_prb,
"FIX_X",
"U", 0, 0);
636 is_all_bc = bc_mng->getBlockIS(name_prb,
"FIX_Y",
"U", 1, 1, is_all_bc);
637 is_all_bc = bc_mng->getBlockIS(name_prb,
"FIX_Z",
"U", 2, 2, is_all_bc);
638 is_all_bc = bc_mng->getBlockIS(name_prb,
"FIX_ALL",
"U", 0, 2, is_all_bc);
640 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
642 <<
"Field split block size " << is_all_bc_size;
643 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
650 auto dm =
simple->getDM();
661 auto solver = pip_mng->createTSIM();
663 CHKERR set_section_monitor(solver);
664 CHKERR set_time_monitor(dm, solver);
665 CHKERR TSSetSolution(solver,
D);
666 CHKERR TSSetFromOptions(solver);
667 CHKERR set_fieldsplit_preconditioner(solver);
669 CHKERR TSSolve(solver, NULL);
671 auto solver = pip_mng->createTSIM2();
672 auto dm =
simple->getDM();
675 CHKERR set_section_monitor(solver);
676 CHKERR set_time_monitor(dm, solver);
677 CHKERR TS2SetSolution(solver,
D, DD);
678 CHKERR TSSetFromOptions(solver);
679 CHKERR set_fieldsplit_preconditioner(solver);
681 CHKERR TSSolve(solver, NULL);
686 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
687 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
704int main(
int argc,
char *argv[]) {
715 auto core_log = logging::core::get();
724 DMType dm_name =
"DMMOFEM";
729 moab::Core mb_instance;
730 moab::Interface &moab = mb_instance;
754 if (Py_FinalizeEx() < 0) {
762 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
763 std::string
field_name, std::string block_name,
764 boost::shared_ptr<MatrixDouble> mat_D_Ptr,
Sev sev) {
768 OpMatBlocks(std::string
field_name, boost::shared_ptr<MatrixDouble>
m,
771 std::vector<const CubitMeshSets *> meshset_vec_ptr)
775 std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]),
false);
777 "Can not get data from block");
784 for (
auto &b : blockData) {
786 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
787 CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
792 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
797 boost::shared_ptr<MatrixDouble> matDPtr;
801 double shearModulusG;
805 double bulkModulusKDefault;
806 double shearModulusGDefault;
807 std::vector<BlockData> blockData;
811 std::vector<const CubitMeshSets *> meshset_vec_ptr,
815 for (
auto m : meshset_vec_ptr) {
817 std::vector<double> block_data;
818 CHKERR m->getAttributes(block_data);
819 if (block_data.size() != 2) {
821 "Expected that block has two attribute");
823 auto get_block_ents = [&]() {
826 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
845 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
849 auto set_material_stiffness = [&]() {
859 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
868 set_material_stiffness();
875 pipeline.push_back(
new OpMatBlocks(
881 (boost::format(
"%s(.*)") % block_name).str()
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto createDMVector(DM dm)
Get smart vector from DM.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
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.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DomainEle::UserDataOperator DomainEleOp
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Deprecated interface functions.
Definition of the displacement bc data structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
Essential boundary conditions.
default operator for TRI element
@ OPROW
operator doWork function is executed on FE rows
Section manager is used to create indexes and sections.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Interface for managing meshsets containing materials and boundary conditions.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Enforce essential constrains on lhs.
Enforce essential constrains on rhs.
Set indices on entities on finite element.
PipelineManager interface.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.