#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
};
};
GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
GAUSS>::OpSource<1, SPACE_DIM>;
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
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 displacents";
break;
default:
break;
}
auto skin_edges = getEntsOnMeshSkin<SPACE_DIM>();
Range boundary_ents;
ParallelComm *pcomm =
if (pcomm == NULL) {
"Communicator not created");
}
CHKERR pcomm->filter_pstatus(skin_edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
}
auto set_matrial_stiffness = [&]() {
: 1;
auto t_D =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
commonDataPtr->mDPtr);
};
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>();
CHKERR set_matrial_stiffness();
}
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"NO_CONTACT", "SIGMA", 0, 3);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_X",
"U",
0, 0);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_Y",
"U",
1, 1);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_Z",
"U",
2, 2);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_ALL",
"U", 0, 3);
boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{
"FIX_"});
}
auto add_domain_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
pipeline.push_back(
} else {
pipeline.push_back(
}
};
auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
auto add_domain_ops_lhs = [&](auto &pipeline) {
new OpKPiola(
"U",
"U", henky_common_data_ptr->getMatTangent()));
} else {
}
return rho * fe_domain_lhs->ts_aa;
};
new OpMass(
"U",
"U", get_rho));
}
auto unity = []() { return 1; };
pipeline.push_back(
new OpMixDivULhs(
"SIGMA",
"U", unity,
true));
};
auto add_domain_ops_rhs = [&](auto &pipeline) {
const auto time = fe_domain_rhs->ts_t;
t_source(1) = 1.0 * time;
return t_source;
};
pipeline.push_back(
new OpBodyForce(
"U", get_body_force));
"U", henky_common_data_ptr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(
}
pipeline.push_back(
pipeline.push_back(
[](double, double, double) { return 1; }));
pipeline.push_back(
pipeline.push_back(
auto mat_acceleration = boost::make_shared<MatrixDouble>();
mat_acceleration));
"U", mat_acceleration, [](
double,
double,
double) {
return rho; }));
}
};
auto add_boundary_base_ops = [&](auto &pipeline) {
};
auto add_boundary_ops_lhs = [&](auto &pipeline) {
auto &bc_map = bc_mng->getBcMapByBlockName();
for (auto bc : bc_map) {
if (bc_mng->checkBlock(bc, "FIX_")) {
<< "Set boundary matrix for " << bc.first;
pipeline.push_back(
new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
"U", "U", [](double, double, double) { return 1.; },
bc.second->getBcEntsPtr()));
}
}
pipeline.push_back(
pipeline.push_back(
"U", "U",
));
};
return -fe_domain_rhs->ts_t;
};
auto add_boundary_ops_rhs = [&](auto &pipeline) {
if (bc_mng->checkBlock(bc, "FIX_")) {
<< "Set boundary residual for " << bc.first;
pipeline.push_back(
new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
auto attr_vec = boost::make_shared<MatrixDouble>(
SPACE_DIM, 1);
attr_vec->clear();
if (bc.second->bcAttributes.size() !=
SPACE_DIM)
"Wrong size of boundary attributes vector. Set right block "
"size attributes. Size of attributes %d",
bc.second->bcAttributes.size());
std::copy(&bc.second->bcAttributes[0],
attr_vec->data().begin());
bc.second->getBcEntsPtr()));
[](double, double, double) { return 1.; },
bc.second->getBcEntsPtr()));
}
}
[this](double, double, double) { return spring_stiffness; }));
};
};
};
}
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 = pipeline_mng->createTSIM();
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetFromOptions(solver);
} else {
auto solver = pipeline_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);
}
Range body_ents;
Range skin_ents;
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";
}
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#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.
#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< 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)
DeprecatedCoreInterface Interface
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
Range getEntsOnMeshSkin()
[Check]
boost::shared_ptr< PostProcEle > postProcFe
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
FieldApproximationBase base
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode createCommonData()
[Set up problem]
Example(MoFEM::Interface &m_field)
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
MoFEMErrorCode postProcess()
[Solve]
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode bC()
[Create common data]
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.
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.
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 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()
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.
Volume finite element base.