23 #ifndef EXECUTABLE_DIMENSION
24 #define EXECUTABLE_DIMENSION 3
30 using namespace MoFEM;
84 GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
99 GAUSS>::OpSource<1, SPACE_DIM>;
104 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
127 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
138 constexpr
double cn = 0.01;
166 boost::shared_ptr<PostProcEle> postProcFe;
171 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
173 template <
int DIM> Range getEntsOnMeshSkin();
180 CHKERR createCommonData();
196 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
197 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz"};
198 PetscInt choice_base_value = AINSWORTH;
200 LASBASETOPT, &choice_base_value, PETSC_NULL);
203 switch (choice_base_value) {
207 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
212 <<
"Set DEMKOWICZ_JACOBI_BASE for displacents";
232 auto skin_edges = getEntsOnMeshSkin<SPACE_DIM>();
236 ParallelComm *pcomm =
240 "Communicator not created");
243 CHKERR pcomm->filter_pstatus(skin_edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
244 PSTATUS_NOT, -1, &boundary_ents);
263 auto set_matrial_stiffness = [&]() {
279 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*commonDataPtr->mDPtr);
286 commonDataPtr = boost::make_shared<ContactOps::CommonData>();
288 commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
289 commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
290 commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
291 commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
292 commonDataPtr->contactStressDivergencePtr =
293 boost::make_shared<MatrixDouble>();
294 commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
295 commonDataPtr->contactDispPtr = boost::make_shared<MatrixDouble>();
296 commonDataPtr->curlContactStressPtr = boost::make_shared<MatrixDouble>();
300 commonDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
301 commonDataPtr->mDPtr->resize(size_symm * size_symm, 1);
303 CHKERR set_matrial_stiffness();
311 auto bc_mng = mField.getInterface<
BcManager>();
314 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
315 "NO_CONTACT",
"SIGMA", 0, 3);
317 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
319 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
321 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
323 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
324 "REMOVE_ALL",
"U", 0, 3);
326 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_X",
"U",
328 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_Y",
"U",
330 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_Z",
"U",
332 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_ALL",
335 boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{
"FIX_"});
347 auto add_domain_base_ops = [&](
auto &pipeline) {
348 auto det_ptr = boost::make_shared<VectorDouble>();
349 auto jac_ptr = boost::make_shared<MatrixDouble>();
350 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
370 auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
371 henky_common_data_ptr->matGradPtr = commonDataPtr->mGradPtr;
372 henky_common_data_ptr->matDPtr = commonDataPtr->mDPtr;
374 auto add_domain_ops_lhs = [&](
auto &pipeline) {
375 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
380 "U", commonDataPtr->mGradPtr));
394 new OpKPiola(
"U",
"U", henky_common_data_ptr->getMatTangent()));
396 pipeline.push_back(
new OpKCauchy(
"U",
"U", commonDataPtr->mDPtr));
401 auto get_rho = [
this](
const double,
const double,
const double) {
404 return rho * fe_domain_lhs->ts_aa;
407 new OpMass(
"U",
"U", get_rho));
410 auto unity = []() {
return 1; };
411 pipeline.push_back(
new OpMixDivULhs(
"SIGMA",
"U", unity,
true));
417 auto add_domain_ops_rhs = [&](
auto &pipeline) {
418 auto get_body_force = [
this](
const double,
const double,
const double) {
423 const auto time = fe_domain_rhs->ts_t;
427 t_source(1) = 1.0 * time;
431 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
433 pipeline.push_back(
new OpBodyForce(
"U", get_body_force));
435 "U", commonDataPtr->mGradPtr));
449 "U", henky_common_data_ptr->getMatFirstPiolaStress()));
452 "U", commonDataPtr->mGradPtr, commonDataPtr->mStrainPtr));
454 "U", commonDataPtr->mStrainPtr, commonDataPtr->mStressPtr,
455 commonDataPtr->mDPtr));
461 "U", commonDataPtr->contactDispPtr));
464 "SIGMA", commonDataPtr->contactStressPtr));
467 "SIGMA", commonDataPtr->contactStressDivergencePtr));
470 new OpMixDivURhs(
"SIGMA", commonDataPtr->contactDispPtr,
471 [](
double,
double,
double) { return 1; }));
476 "U", commonDataPtr->contactStressDivergencePtr));
482 auto mat_acceleration = boost::make_shared<MatrixDouble>();
487 "U", mat_acceleration, [](
double,
double,
double) {
return rho; }));
493 auto add_boundary_base_ops = [&](
auto &pipeline) {
498 "U", commonDataPtr->contactDispPtr));
500 "SIGMA", commonDataPtr->contactTractionPtr));
503 auto add_boundary_ops_lhs = [&](
auto &pipeline) {
505 auto &bc_map = bc_mng->getBcMapByBlockName();
506 for (
auto bc : bc_map) {
507 if (bc_mng->checkBlock(bc,
"FIX_")) {
509 <<
"Set boundary matrix for " << bc.first;
511 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
513 "U",
"U", [](
double,
double,
double) {
return 1.; },
514 bc.second->getBcEntsPtr()));
519 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
535 auto time_scaled = [&](double, double, double) {
538 return -fe_domain_rhs->ts_t;
541 auto add_boundary_ops_rhs = [&](
auto &pipeline) {
544 if (bc_mng->checkBlock(bc,
"FIX_")) {
546 <<
"Set boundary residual for " << bc.first;
548 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
549 auto attr_vec = boost::make_shared<MatrixDouble>(
SPACE_DIM, 1);
551 if (bc.second->bcAttributes.size() !=
SPACE_DIM)
553 "Wrong size of boundary attributes vector. Set right block "
554 "size attributes. Size of attributes %d",
555 bc.second->bcAttributes.size());
556 std::copy(&bc.second->bcAttributes[0],
558 attr_vec->data().begin());
560 pipeline.push_back(
new OpBoundaryVec(
"U", attr_vec, time_scaled,
561 bc.second->getBcEntsPtr()));
563 "U", commonDataPtr->contactDispPtr,
564 [](
double,
double,
double) { return 1.; },
565 bc.second->getBcEntsPtr()));
571 pipeline.push_back(
new OpSetBc(
"U",
true, boundaryMarker));
574 "U", commonDataPtr->contactDispPtr,
575 [
this](
double,
double,
double) { return spring_stiffness; }));
613 auto set_section_monitor = [&](
auto solver) {
616 CHKERR TSGetSNES(solver, &snes);
617 PetscViewerAndFormat *vf;
618 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
619 PETSC_VIEWER_DEFAULT, &vf);
622 (
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
627 auto scatter_create = [&](
auto D,
auto coeff) {
629 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
630 ROW,
"U", coeff, coeff, is);
632 CHKERR ISGetLocalSize(is, &loc_size);
634 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &
v);
636 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
641 auto set_time_monitor = [&](
auto dm,
auto solver) {
643 boost::shared_ptr<Monitor> monitor_ptr(
644 new Monitor(dm, commonDataPtr, uXScatter, uYScatter, uZScatter));
645 boost::shared_ptr<ForcesAndSourcesCore>
null;
647 monitor_ptr,
null,
null);
651 auto dm =
simple->getDM();
653 uXScatter = scatter_create(
D, 0);
654 uYScatter = scatter_create(
D, 1);
656 uZScatter = scatter_create(
D, 2);
659 auto solver = pipeline_mng->createTSIM();
661 CHKERR set_section_monitor(solver);
662 CHKERR set_time_monitor(dm, solver);
663 CHKERR TSSetSolution(solver,
D);
664 CHKERR TSSetFromOptions(solver);
666 CHKERR TSSolve(solver, NULL);
668 auto solver = pipeline_mng->createTSIM2();
669 auto dm =
simple->getDM();
672 CHKERR set_section_monitor(solver);
673 CHKERR set_time_monitor(dm, solver);
674 CHKERR TS2SetSolution(solver,
D, DD);
675 CHKERR TSSetFromOptions(solver);
677 CHKERR TSSolve(solver, NULL);
680 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
681 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
698 CHKERR mField.get_moab().get_entities_by_dimension(0, DIM, body_ents);
699 Skinner skin(&mField.get_moab());
701 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
706 static char help[] =
"...\n\n";
708 int main(
int argc,
char *argv[]) {
715 auto core_log = logging::core::get();
724 DMType dm_name =
"DMMOFEM";
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
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.
auto smartCreateDMVector
Get smart vector from DM.
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.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
PipelineManager::FaceEle DomainEle
Range getEntsOnMeshSkin()
[Check]
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
MoFEMErrorCode checkResults()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEMErrorCode postProcess()
[Solve]
MoFEMErrorCode setupProblem()
Simple interface for fast problem set-up.
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
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.
Deprecated interface functions.
default operator for EDGE element
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
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.
Calculate jacobian on Hex or other volume which is not simplex.
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 valuse for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Make Hdiv space from Hcurl space in 2d.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
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.
Volume finite element with switches.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator