using DomainEle = FaceElementForcesAndSourcesCoreBase;
};
using DomainEle = VolumeElementForcesAndSourcesCoreBase;
};
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
GAUSS>::OpSource<1, SPACE_DIM>;
constexpr
double rho = 1;
constexpr
double omega = 2.4;
private:
boost::shared_ptr<MatrixDouble> matGradPtr;
boost::shared_ptr<MatrixDouble> matStrainPtr;
boost::shared_ptr<MatrixDouble> matStressPtr;
boost::shared_ptr<MatrixDouble> matAccelerationPtr;
boost::shared_ptr<MatrixDouble> matInertiaPtr;
boost::shared_ptr<MatrixDouble> matDPtr;
};
auto set_matrial_stiffens = [&]() {
constexpr double A =
: 1;
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
};
matGradPtr = boost::make_shared<MatrixDouble>();
matStrainPtr = boost::make_shared<MatrixDouble>();
matStressPtr = boost::make_shared<MatrixDouble>();
matAccelerationPtr = boost::make_shared<MatrixDouble>();
matInertiaPtr = boost::make_shared<MatrixDouble>();
matDPtr = boost::make_shared<MatrixDouble>();
matDPtr->resize(size_symm * size_symm, 1);
CHKERR set_matrial_stiffens();
}
}
auto simple = mField.getInterface<Simple>();
}
Simple *
simple = mField.getInterface<Simple>();
}
auto fix_disp = [&](const std::string blockset_name) {
Range fix_ents;
if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
0) {
CHKERR mField.get_moab().get_entities_by_handle(it->meshset, fix_ents,
true);
}
}
return fix_ents;
};
auto remove_ents = [&](const Range &&ents, const int lo, const int hi) {
auto prb_mng = mField.getInterface<ProblemsManager>();
auto simple = mField.getInterface<Simple>();
Range verts;
CHKERR mField.get_moab().get_connectivity(ents, verts,
true);
verts.merge(ents);
Range adj;
CHKERR mField.get_moab().get_adjacencies(ents, 1,
false, adj,
moab::Interface::UNION);
verts.merge(adj);
};
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
CHKERR prb_mng->removeDofsOnEntities(
simple->getProblemName(),
"U", verts,
lo, hi);
};
CHKERR remove_ents(fix_disp(
"FIX_X"), 0, 0);
CHKERR remove_ents(fix_disp(
"FIX_Y"), 1, 1);
CHKERR remove_ents(fix_disp(
"FIX_Z"), 2, 2);
CHKERR remove_ents(fix_disp(
"FIX_ALL"), 0, 3);
}
auto *
simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetInvJacH1ForFace(
invJac));
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacH1ForFace(
invJac));
}
auto get_rho = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
return rho * fe_domain_lhs->ts_aa;
};
auto get_body_force = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
const auto time = fe_domain_rhs->ts_t;
t_source(
i) *= sin(time *
omega * M_PI);
return t_source;
};
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpK(
"U",
"U", matDPtr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMass(
"U",
"U", get_rho));
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
matGradPtr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSymmetrizeTensor<SPACE_DIM>("U", matGradPtr, matStrainPtr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
"U", matStrainPtr, matStressPtr, matDPtr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>("U",
matAccelerationPtr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpScaleMatrix(
"U",
rho, matAccelerationPtr, matInertiaPtr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
}
Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
: dM(dm), postProc(post_proc){};
"out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
}
}
private:
SmartPetscObj<DM> dM;
boost::shared_ptr<PostProcEle> postProc;
};
auto *
simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto ts = pipeline_mng->createTS2();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
post_proc_fe->generateReferenceElementMesh();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(
invJac));
}
post_proc_fe->getOpPtrVector().push_back(
matGradPtr));
post_proc_fe->getOpPtrVector().push_back(
new OpSymmetrizeTensor<SPACE_DIM>("U", matGradPtr, matStrainPtr));
post_proc_fe->getOpPtrVector().push_back(
new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
"U", matStrainPtr, matStressPtr, matDPtr));
post_proc_fe->getOpPtrVector().push_back(new OpPostProcElastic<SPACE_DIM>(
"U", post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts, matStrainPtr,
matStressPtr));
post_proc_fe->addFieldValuesPostProc("U");
boost::shared_ptr<FEMethod> null_fe;
auto monitor_ptr = boost::make_shared<Monitor>(dm, post_proc_fe);
null_fe, monitor_ptr);
double ftime = 1;
CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
SCATTER_FORWARD);
PetscInt steps, snesfails, rejects, nonlinits, linits;
CHKERR TSGetTimeStepNumber(ts, &steps);
CHKERR TSGetSNESFailures(ts, &snesfails);
CHKERR TSGetStepRejections(ts, &rejects);
CHKERR TSGetSNESIterations(ts, &nonlinits);
CHKERR TSGetKSPIterations(ts, &linits);
"steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
"%D, linits %D\n",
steps, rejects, snesfails, ftime, nonlinits, linits);
}
PetscBool test_flg = PETSC_FALSE;
if (test_flg) {
auto *
simple = mField.getInterface<Simple>();
SCATTER_FORWARD);
double nrm2;
MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Regression norm " << nrm2;
constexpr double regression_value = 1.09572;
if (fabs(nrm2 - regression_value) > 1e-2)
"Regression test faileed; wrong norm value.");
}
}
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
try {
DMType dm_name = "DMMOFEM";
}
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
int main(int argc, char *argv[])
static char help[]
[Check]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
static double * ts_aa_ptr
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
DataForcesAndSourcesCore::EntData EntData
static double * ts_time_ptr
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
constexpr double young_modulus
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
auto smartCreateDMVector
Get smart vector from DM.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in 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.
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.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMass
VolumeElementForcesAndSourcesCore DomainEle
const int save_every_nth_step
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
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.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
static constexpr int approx_order
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Applying essential BC]
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode runProblem()
[Run problem]
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode outputResults()
[Solve]
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
Deprecated interface functions.
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.
default operator for TET element
[Push operators to pipeline]