#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
IntegrationType::GAUSS;
};
};
GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
boost::shared_ptr<MatrixDouble> mat_D_Ptr,
Sev sev);
};
constexpr double rho = 0;
constexpr double cn = 0.01;
private:
};
}
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
PetscInt choice_base_value = AINSWORTH;
LASBASETOPT, &choice_base_value, PETSC_NULL);
switch (choice_base_value) {
case AINSWORTH:
<< "Set AINSWORTH_LEGENDRE_BASE for displacements";
break;
case DEMKOWICZ:
<< "Set DEMKOWICZ_JACOBI_BASE for displacements";
break;
default:
break;
}
auto skin_edges = getEntsOnMeshSkin<SPACE_DIM>();
ParallelComm *pcomm =
if (pcomm == NULL) {
"Communicator not created");
}
CHKERR pcomm->filter_pstatus(skin_edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
}
commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
boost::make_shared<MatrixDouble>();
commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->contactDispPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->curlContactStressPtr = boost::make_shared<MatrixDouble>();
}
for (auto f : {"U", "SIGMA"}) {
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
}
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"NO_CONTACT", "SIGMA", 0, 3);
simple->getProblemName(),
"U");
}
auto time_scale = boost::make_shared<TimeScale>();
auto add_domain_base_ops = [&](auto &pip) {
};
auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
auto add_domain_ops_lhs = [&](auto &pip) {
Sev::verbose);
pip.push_back(
pip.push_back(
pip.push_back(
pip.push_back(
pip.push_back(
new OpKPiola(
"U",
"U", henky_common_data_ptr->getMatTangent()));
return rho * fe_domain_lhs->ts_aa;
};
pip.push_back(
new OpMass(
"U",
"U", get_rho));
}
auto unity = []() { return 1; };
};
auto add_domain_ops_rhs = [&](auto &pip) {
CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
pip,
mField,
"U", {time_scale}, Sev::inform);
Sev::inform);
pip.push_back(
pip.push_back(
pip.push_back(
pip.push_back(
"U", henky_common_data_ptr->getMatFirstPiolaStress()));
[](double, double, double) { return 1; }));
pip.push_back(
auto mat_acceleration = boost::make_shared<MatrixDouble>();
"U", mat_acceleration));
"U", mat_acceleration, [](
double,
double,
double) {
return rho; }));
}
};
auto add_boundary_base_ops = [&](auto &pip) {
};
auto add_boundary_ops_lhs = [&](auto &pip) {
CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
pip,
mField,
"U", Sev::inform);
pip.push_back(
"U", "U",
));
};
auto add_boundary_ops_rhs = [&](auto &pip) {
CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
pip,
mField,
"U", {time_scale}, Sev::inform);
[this](double, double, double) { return spring_stiffness; }));
};
};
};
auto get_bc_hook_rhs = [&]() {
return hook;
};
auto get_bc_hook_lhs = [&]() {
return hook;
};
}
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
(
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
};
auto scatter_create = [&](
auto D,
auto coeff) {
CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
ROW,
"U", coeff, coeff, is);
int loc_size;
CHKERR ISGetLocalSize(is, &loc_size);
VecScatter scatter;
CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
boost::shared_ptr<ForcesAndSourcesCore> null;
monitor_ptr, null, null);
};
auto solver = pip_mng->createTSIM();
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetFromOptions(solver);
} else {
auto solver = pip_mng->createTSIM2();
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TS2SetSolution(solver,
D, DD);
CHKERR TSSetFromOptions(solver);
}
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
}
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
boost::shared_ptr<MatrixDouble> mat_D_Ptr,
Sev sev) {
OpMatBlocks(std::string
field_name, boost::shared_ptr<MatrixDouble>
m,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
"Can not get data from block");
}
for (auto &b : blockData) {
if (
b.blockEnts.find(getFEEntityHandle()) !=
b.blockEnts.end()) {
CHKERR getMatDPtr(matDPtr,
b.bulkModulusK,
b.shearModulusG);
}
}
CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
}
private:
boost::shared_ptr<MatrixDouble> matDPtr;
double bulkModulusK;
double shearModulusG;
};
double bulkModulusKDefault;
double shearModulusGDefault;
std::vector<BlockData> blockData;
std::vector<const CubitMeshSets *> meshset_vec_ptr,
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() != 2) {
"Expected that block has two attribute");
}
auto get_block_ents = [&]() {
m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
return ents;
};
blockData.push_back(
}
}
auto set_material_stiffness = [&]() {
: 1;
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
};
set_material_stiffness();
}
};
pipeline.push_back(new OpMatBlocks(
(boost::format("%s(.*)") % block_name).str()
))
));
}
#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.
NaturalBC< BoundaryEleOp >::Assembly< A >::LinearForm< I > BoundaryRhsBCs
NaturalBC< DomainEleOp >::Assembly< A >::LinearForm< I > DomainRhsBCs
BoundaryLhsBCs::OpFlux< BoundaryBCs, 1, SPACE_DIM > OpBoundaryLhsBCs
NaturalBC< BoundaryEleOp >::Assembly< A >::BiLinearForm< I > BoundaryLhsBCs
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
DomainRhsBCs::OpFlux< DomainBCs, 1, SPACE_DIM > OpDomainRhsBCs
BoundaryRhsBCs::OpFlux< BoundaryBCs, 1, SPACE_DIM > OpBoundaryRhsBCs
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 smartCreateDMVector(DM dm)
Get smart vector from DM.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
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
const FTensor::Tensor2< T, Dim, Dim > Vec
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.
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
constexpr IntegrationType I
#define EXECUTABLE_DIMENSION
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Range getEntsOnMeshSkin()
[Check]
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
MoFEMErrorCode tsSolve()
[Push operators to pip]
FieldApproximationBase base
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
MoFEMErrorCode postProcess()
[Solve]
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode addMatBlockOps(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
MoFEMErrorCode bC()
[Create common data]
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
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)
Class (Function) to enforce essential constrains.
default operator for TRI element
Section manager is used to create indexes and sections.
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.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.