32 IntegrationType::GAUSS;
42 I>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
44 I>::OpGradTimesSymTensor<1, SPACE_DIM, SPACE_DIM>;
83 boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator> &pipeline,
84 std::string
field_name, std::string block_name,
85 boost::shared_ptr<MatrixDouble> mat_D_Ptr,
Sev sev);
89 boost::ptr_vector<ForcesAndSourcesCore::UserDataOperator> &pipeline,
90 std::string
field_name, std::string block_name,
91 boost::shared_ptr<MatrixDouble> mat_D_Ptr,
Sev sev) {
95 OpMatBlocks(std::string
field_name, boost::shared_ptr<MatrixDouble>
m,
98 std::vector<const CubitMeshSets *> meshset_vec_ptr)
102 std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]),
false);
104 "Can not get data from block");
111 for (
auto &b : blockData) {
113 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
114 CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
119 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
124 boost::shared_ptr<MatrixDouble> matDPtr;
128 double shearModulusG;
132 double bulkModulusKDefault;
133 double shearModulusGDefault;
134 std::vector<BlockData> blockData;
138 std::vector<const CubitMeshSets *> meshset_vec_ptr,
142 for (
auto m : meshset_vec_ptr) {
144 std::vector<double> block_data;
145 CHKERR m->getAttributes(block_data);
146 if (block_data.size() != 2) {
148 "Expected that block has two attribute");
150 auto get_block_ents = [&]() {
153 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
172 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
176 auto set_material_stiffness = [&]() {
186 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
195 set_material_stiffness();
200 pipeline.push_back(
new OpMatBlocks(
206 (boost::format(
"%s(.*)") % block_name).str()
256 auto project_ho_geometry = [&]() {
260 CHKERR project_ho_geometry();
273 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
275 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
277 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
279 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
280 "REMOVE_ALL",
"U", 0, 3);
283 pipeline_mng->getOpDomainRhsPipeline(), {H1},
"GEOMETRY");
285 pipeline_mng->getOpBoundaryRhsPipeline(), {NOSPACE},
"GEOMETRY");
288 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
289 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
290 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
291 pipeline_mng->getOpDomainRhsPipeline().push_back(
294 pipeline_mng->getOpDomainRhsPipeline().push_back(
297 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
299 "MAT_ELASTIC", mat_D_ptr, Sev::inform);
300 pipeline_mng->getOpDomainRhsPipeline().push_back(
302 "U", mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
303 pipeline_mng->getOpDomainRhsPipeline().push_back(
305 [](
double,
double,
double)
constexpr {
return -1; }));
309 pipeline_mng->getOpDomainRhsPipeline(),
mField,
"U", Sev::inform);
312 pipeline_mng->getOpBoundaryRhsPipeline(),
mField,
"U", -1, Sev::inform);
314 pipeline_mng->getOpBoundaryLhsPipeline(),
mField,
"U", Sev::verbose);
318 simple->getProblemName(),
"U");
320 auto get_pre_proc_hook = [&]() {
322 mField, pipeline_mng->getDomainRhsFE(), {});
325 pipeline_mng->getDomainRhsFE()->preProcessHook = get_pre_proc_hook();
326 pipeline_mng->getBoundaryRhsFE()->preProcessHook = get_pre_proc_hook();
340 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
342 "MAT_ELASTIC", mat_D_ptr, Sev::verbose);
344 new OpK(
"U",
"U", mat_D_ptr));
363 auto solver = pipeline_mng->createKSP();
364 CHKERR KSPSetFromOptions(solver);
367 auto dm =
simple->getDM();
372 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
373 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
383 auto det_ptr = boost::make_shared<VectorDouble>();
384 auto jac_ptr = boost::make_shared<MatrixDouble>();
385 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
391 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
394 post_proc_fe->getOpPtrVector(), {H1},
"GEOMETRY");
396 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
397 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
398 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
400 post_proc_fe->getOpPtrVector().push_back(
403 post_proc_fe->getOpPtrVector().push_back(
406 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
408 mat_D_ptr, Sev::verbose);
409 post_proc_fe->getOpPtrVector().push_back(
411 "U", mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
413 auto u_ptr = boost::make_shared<MatrixDouble>();
414 post_proc_fe->getOpPtrVector().push_back(
416 auto x_ptr = boost::make_shared<MatrixDouble>();
417 post_proc_fe->getOpPtrVector().push_back(
422 post_proc_fe->getOpPtrVector().push_back(
426 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
430 {{
"U", u_ptr}, {
"GEOMETRY", x_ptr}},
434 {{
"STRAIN", mat_strain_ptr}, {
"STRESS", mat_stress_ptr}}
443 CHKERR post_proc_fe->writeFile(
"out_elastic.h5m");
454 pipeline_mng->getDomainRhsFE().reset();
455 pipeline_mng->getDomainLhsFE().reset();
456 pipeline_mng->getBoundaryRhsFE().reset();
457 pipeline_mng->getBoundaryLhsFE().reset();
464 pipeline_mng->getOpDomainRhsPipeline(), {H1},
"GEOMETRY");
466 pipeline_mng->getOpBoundaryRhsPipeline(), {NOSPACE},
"GEOMETRY");
468 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
469 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
470 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
472 pipeline_mng->getOpDomainRhsPipeline().push_back(
475 pipeline_mng->getOpDomainRhsPipeline().push_back(
478 auto mat_D_ptr = boost::make_shared<MatrixDouble>();
480 "MAT_ELASTIC", mat_D_ptr, Sev::verbose);
481 pipeline_mng->getOpDomainRhsPipeline().push_back(
483 "U", mat_strain_ptr, mat_stress_ptr, mat_D_ptr));
484 pipeline_mng->getOpDomainRhsPipeline().push_back(
488 pipeline_mng->getOpDomainRhsPipeline(),
mField,
"U", Sev::verbose);
490 pipeline_mng->getOpBoundaryRhsPipeline(),
mField,
"U", 1, Sev::verbose);
492 auto dm =
simple->getDM();
494 pipeline_mng->getDomainRhsFE()->ksp_f = res;
495 pipeline_mng->getBoundaryRhsFE()->ksp_f = res;
497 CHKERR VecZeroEntries(res);
500 CHKERR pipeline_mng->loopFiniteElements();
503 CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
504 CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
505 CHKERR VecAssemblyBegin(res);
506 CHKERR VecAssemblyEnd(res);
509 CHKERR VecNorm(res, NORM_2, &nrm2);
511 MOFEM_LOG_C(
"WORLD", Sev::inform,
"residual = %3.4e\n", nrm2);
512 constexpr double eps = 1e-8;
522int main(
int argc,
char *argv[]) {
531 DMType dm_name =
"DMMOFEM";
536 moab::Core mb_instance;
537 moab::Interface &moab = mb_instance;
Implementation of elastic spring bc.
#define MOFEM_LOG_C(channel, severity, format,...)
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Implementation of natural boundary conditions.
Boundary conditions in domain, i.e. body forces.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ 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()
#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.
static char help[]
[Check]
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
constexpr int SPACE_DIM
[Define dimension]
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr IntegrationType I
constexpr double bulk_modulus_K
constexpr double young_modulus
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.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SeverityLevel
Severity levels.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
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
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
constexpr auto field_name
static constexpr int approx_order
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode addMatBlockOps(boost::ptr_vector< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
virtual moab::Interface & get_moab()=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.
Data on single entity (This is passed as argument to DataOperator::doWork)
Class (Function) to enforce essential constrains.
default operator for TRI element
@ OPROW
operator doWork function is executed on FE rows
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.