11#ifndef EXECUTABLE_DIMENSION
12#define EXECUTABLE_DIMENSION 3
70 GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
85 GAUSS>::OpSource<1, SPACE_DIM>;
90 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
113 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
124constexpr double cn = 0.01;
179 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
180 const char *list_bases[LASBASETOPT] = {
"ainsworth",
"demkowicz"};
181 PetscInt choice_base_value = AINSWORTH;
183 LASBASETOPT, &choice_base_value, PETSC_NULL);
186 switch (choice_base_value) {
190 <<
"Set AINSWORTH_LEGENDRE_BASE for displacements";
195 <<
"Set DEMKOWICZ_JACOBI_BASE for displacents";
215 auto skin_edges = getEntsOnMeshSkin<SPACE_DIM>();
219 ParallelComm *pcomm =
223 "Communicator not created");
226 CHKERR pcomm->filter_pstatus(skin_edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
227 PSTATUS_NOT, -1, &boundary_ents);
246 auto set_matrial_stiffness = [&]() {
262 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
commonDataPtr->mDPtr);
269 commonDataPtr = boost::make_shared<ContactOps::CommonData>();
271 commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
272 commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
273 commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
274 commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
276 boost::make_shared<MatrixDouble>();
277 commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
278 commonDataPtr->contactDispPtr = boost::make_shared<MatrixDouble>();
279 commonDataPtr->curlContactStressPtr = boost::make_shared<MatrixDouble>();
286 CHKERR set_matrial_stiffness();
297 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
298 "NO_CONTACT",
"SIGMA", 0, 3);
300 simple->getProblemName(),
"U");
312 auto add_domain_base_ops = [&](
auto &pipeline) {
313 auto det_ptr = boost::make_shared<VectorDouble>();
314 auto jac_ptr = boost::make_shared<MatrixDouble>();
315 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
336 auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
340 auto add_domain_ops_lhs = [&](
auto &pipeline) {
359 new OpKPiola(
"U",
"U", henky_common_data_ptr->getMatTangent()));
369 return rho * fe_domain_lhs->ts_aa;
372 new OpMass(
"U",
"U", get_rho));
375 auto unity = []() {
return 1; };
376 pipeline.push_back(
new OpMixDivULhs(
"SIGMA",
"U", unity,
true));
381 auto add_domain_ops_rhs = [&](
auto &pipeline) {
387 const auto time = fe_domain_rhs->ts_t;
391 t_source(1) = 1.0 * time;
395 pipeline.push_back(
new OpBodyForce(
"U", get_body_force));
411 "U", henky_common_data_ptr->getMatFirstPiolaStress()));
433 [](
double,
double,
double) { return 1; }));
444 auto mat_acceleration = boost::make_shared<MatrixDouble>();
449 "U", mat_acceleration, [](
double,
double,
double) {
return rho; }));
454 auto add_boundary_base_ops = [&](
auto &pipeline) {
464 auto add_boundary_ops_lhs = [&](
auto &pipeline) {
466 auto &bc_map = bc_mng->getBcMapByBlockName();
467 for (
auto bc : bc_map) {
468 if (bc_mng->checkBlock(bc,
"FIX_")) {
470 <<
"Set boundary matrix for " << bc.first;
472 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
474 "U",
"U", [](
double,
double,
double) {
return 1.; },
475 bc.second->getBcEntsPtr()));
496 return -fe_domain_rhs->ts_t;
499 auto add_boundary_ops_rhs = [&](
auto &pipeline) {
502 if (bc_mng->checkBlock(bc,
"FIX_")) {
504 <<
"Set boundary residual for " << bc.first;
506 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
507 auto attr_vec = boost::make_shared<MatrixDouble>(
SPACE_DIM, 1);
509 if (bc.second->bcAttributes.size() !=
SPACE_DIM)
511 "Wrong size of boundary attributes vector. Set right block "
512 "size attributes. Size of attributes %d",
513 bc.second->bcAttributes.size());
514 std::copy(&bc.second->bcAttributes[0],
516 attr_vec->data().begin());
518 pipeline.push_back(
new OpBoundaryVec(
"U", attr_vec, time_scaled,
519 bc.second->getBcEntsPtr()));
522 [](
double,
double,
double) { return 1.; },
523 bc.second->getBcEntsPtr()));
531 [
this](
double,
double,
double) { return spring_stiffness; }));
557 auto get_bc_hook_rhs = [&]() {
560 {boost::make_shared<TimeScale>()});
564 auto get_bc_hook_lhs = [&]() {
567 {boost::make_shared<TimeScale>()});
571 pipeline_mng->
getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
572 pipeline_mng->
getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
586 auto set_section_monitor = [&](
auto solver) {
589 CHKERR TSGetSNES(solver, &snes);
590 PetscViewerAndFormat *vf;
591 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
592 PETSC_VIEWER_DEFAULT, &vf);
595 (
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
600 auto scatter_create = [&](
auto D,
auto coeff) {
602 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
603 ROW,
"U", coeff, coeff, is);
605 CHKERR ISGetLocalSize(is, &loc_size);
609 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
614 auto set_time_monitor = [&](
auto dm,
auto solver) {
616 boost::shared_ptr<Monitor> monitor_ptr(
618 boost::shared_ptr<ForcesAndSourcesCore> null;
620 monitor_ptr, null, null);
624 auto dm =
simple->getDM();
632 auto solver = pipeline_mng->createTSIM();
634 CHKERR set_section_monitor(solver);
635 CHKERR set_time_monitor(dm, solver);
636 CHKERR TSSetSolution(solver,
D);
637 CHKERR TSSetFromOptions(solver);
639 CHKERR TSSolve(solver, NULL);
641 auto solver = pipeline_mng->createTSIM2();
642 auto dm =
simple->getDM();
645 CHKERR set_section_monitor(solver);
646 CHKERR set_time_monitor(dm, solver);
647 CHKERR TS2SetSolution(solver,
D, DD);
648 CHKERR TSSetFromOptions(solver);
650 CHKERR TSSolve(solver, NULL);
653 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
654 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
674 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
681int main(
int argc,
char *argv[]) {
688 auto core_log = logging::core::get();
697 DMType dm_name =
"DMMOFEM";
702 moab::Core mb_instance;
703 moab::Interface &moab = mb_instance;
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 .
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
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_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.
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.
const double v
phase velocity of light in medium (cm/ns)
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 >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
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.
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
PetscBool is_large_strains
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
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
Range getEntsOnMeshSkin()
[Check]
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
FieldApproximationBase base
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
MoFEMErrorCode postProcess()
[Solve]
MoFEMErrorCode setupProblem()
boost::shared_ptr< PostProcEle > postProcFe
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Simple interface for fast problem set-up.
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
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.
Deprecated interface functions.
Definition of the displacement bc data structure.
default operator for EDGE element
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.
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 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.
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 base.