42constexpr double rho = 1;
110 auto project_ho_geometry = [&]() {
114 CHKERR project_ho_geometry();
128 simple->getProblemName(),
"U");
149 auto rho_ptr = boost::make_shared<double>(
rho);
151 auto add_rho_block = [
this, rho_ptr](
auto &pip,
auto block_name,
Sev sev) {
155 OpMatRhoBlocks(boost::shared_ptr<double> rho_ptr,
157 std::vector<const CubitMeshSets *> meshset_vec_ptr)
160 "Can not get data from block");
168 for (
auto &b : blockData) {
169 if (
b.blockEnts.find(getFEEntityHandle()) !=
b.blockEnts.end()) {
186 std::vector<BlockData> blockData;
190 std::vector<const CubitMeshSets *> meshset_vec_ptr,
194 for (
auto m : meshset_vec_ptr) {
196 std::vector<double> block_data;
197 CHKERR m->getAttributes(block_data);
198 if (block_data.size() < 1) {
200 "Expected that block has two attributes");
202 auto get_block_ents = [&]() {
205 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
210 <<
m->getName() <<
": ro = " << block_data[0];
212 blockData.push_back({block_data[0], get_block_ents()});
219 boost::shared_ptr<double> rhoPtr;
222 pip.push_back(
new OpMatRhoBlocks(
228 (boost::format(
"%s(.*)") % block_name).str()
242 auto get_time_scale = [
this](
const double time) {
243 return sin(time *
omega * M_PI);
246 auto apply_lhs = [&](
auto &pip) {
250 CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
251 mField, pip,
"U",
"MAT_ELASTIC", Sev::verbose);
252 CHKERR add_rho_block(pip,
"MAT_RHO", Sev::verbose);
254 pip.push_back(
new OpMass(
"U",
"U", get_rho));
255 static_cast<OpMass &
>(pip.back()).feScalingFun =
260 auto apply_rhs = [&](
auto &pip) {
264 pipeline_mng->getOpDomainRhsPipeline(), {H1},
"GEOMETRY");
265 CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
266 mField, pip,
"U",
"MAT_ELASTIC", Sev::inform);
267 CHKERR add_rho_block(pip,
"MAT_RHO", Sev::inform);
269 auto mat_acceleration_ptr = boost::make_shared<MatrixDouble>();
272 "U", mat_acceleration_ptr));
273 pip.push_back(
new OpInertiaForce(
"U", mat_acceleration_ptr, get_rho));
277 {boost::make_shared<TimeScaleVector<SPACE_DIM>>(
"-time_vector_file",
279 "BODY_FORCE", Sev::inform);
284 CHKERR apply_lhs(pipeline_mng->getOpDomainLhsPipeline());
285 CHKERR apply_rhs(pipeline_mng->getOpDomainRhsPipeline());
293 auto get_bc_hook = [&]() {
295 mField, pipeline_mng->getDomainRhsFE(),
296 {boost::make_shared<TimeScale>()});
300 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook();
321 "out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
328 boost::shared_ptr<PostProcEle>
postProc;
337 auto dm =
simple->getDM();
339 ts = pipeline_mng->createTSIM2();
342 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
344 post_proc_fe->getOpPtrVector(), {H1});
345 auto common_ptr = commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
346 mField, post_proc_fe->getOpPtrVector(),
"U",
"MAT_ELASTIC", Sev::inform);
348 auto u_ptr = boost::make_shared<MatrixDouble>();
349 post_proc_fe->getOpPtrVector().push_back(
351 auto X_ptr = boost::make_shared<MatrixDouble>();
352 post_proc_fe->getOpPtrVector().push_back(
357 post_proc_fe->getOpPtrVector().push_back(
361 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
365 {{
"U", u_ptr}, {
"GEOMETRY", X_ptr}},
367 {{
"GRAD", common_ptr->matGradPtr},
368 {
"FIRST_PIOLA", common_ptr->getMatFirstPiolaStress()}},
377 boost::shared_ptr<FEMethod> null_fe;
378 auto monitor_ptr = boost::make_shared<Monitor>(dm, post_proc_fe);
380 null_fe, monitor_ptr);
384 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
390 CHKERR TS2SetSolution(ts, T, TT);
391 CHKERR TSSetFromOptions(ts);
394 CHKERR TSGetTime(ts, &ftime);
396 PetscInt steps, snesfails, rejects, nonlinits, linits;
397#if PETSC_VERSION_GE(3, 8, 0)
398 CHKERR TSGetStepNumber(ts, &steps);
400 CHKERR TSGetTimeStepNumber(ts, &steps);
402 CHKERR TSGetSNESFailures(ts, &snesfails);
403 CHKERR TSGetStepRejections(ts, &rejects);
404 CHKERR TSGetSNESIterations(ts, &nonlinits);
405 CHKERR TSGetKSPIterations(ts, &linits);
407 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
409 steps, rejects, snesfails, ftime, nonlinits, linits);
418 PetscBool test_flg = PETSC_FALSE;
426 CHKERR VecNorm(T, NORM_2, &nrm2);
427 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Regression norm " << nrm2;
428 constexpr double regression_value = 0.0194561;
429 if (fabs(nrm2 - regression_value) > 1e-2)
431 "Regression test failed; wrong norm value.");
446int main(
int argc,
char *argv[]) {
453 auto core_log = logging::core::get();
462 DMType dm_name =
"DMMOFEM";
467 moab::Core mb_instance;
468 moab::Interface &moab = mb_instance;
#define MOFEM_LOG_C(channel, severity, format,...)
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_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_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
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.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#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
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
static char help[]
[Check]
static double * ts_aa_ptr
constexpr double poisson_ratio
static double * ts_time_ptr
constexpr double young_modulus
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
static constexpr int approx_order
FTensor::Index< 'i', 3 > i
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
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.
structure for User Loop Methods on finite elements
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.
Interface for managing meshsets containing materials and boundary conditions.
Approximate field values for given petsc vector.
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.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
PetscReal ts_aa
shift for U_tt shift for U_tt
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
[Push operators to pipeline]
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcEle > postProc