42constexpr double rho = 1;
107 auto project_ho_geometry = [&]() {
111 CHKERR project_ho_geometry();
125 simple->getProblemName(),
"U");
145 auto rho_ptr = boost::make_shared<double>(
rho);
147 auto add_rho_block = [
this, rho_ptr](
auto &pip,
auto block_name,
Sev sev) {
151 OpMatRhoBlocks(boost::shared_ptr<double> rho_ptr,
153 std::vector<const CubitMeshSets *> meshset_vec_ptr)
156 "Can not get data from block");
164 for (
auto &b : blockData) {
165 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
182 std::vector<BlockData> blockData;
186 std::vector<const CubitMeshSets *> meshset_vec_ptr,
190 for (
auto m : meshset_vec_ptr) {
192 std::vector<double> block_data;
193 CHKERR m->getAttributes(block_data);
194 if (block_data.size() < 1) {
196 "Expected that block has two attributes");
198 auto get_block_ents = [&]() {
201 m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
206 <<
m->getName() <<
": ro = " << block_data[0];
208 blockData.push_back({block_data[0], get_block_ents()});
215 boost::shared_ptr<double> rhoPtr;
218 pip.push_back(
new OpMatRhoBlocks(
224 (boost::format(
"%s(.*)") % block_name).str()
238 auto get_time_scale = [
this](
const double time) {
239 return sin(time *
omega * M_PI);
242 auto apply_lhs = [&](
auto &pip) {
247 mField, pip,
"U",
"MAT_ELASTIC", Sev::verbose);
248 CHKERR add_rho_block(pip,
"MAT_RHO", Sev::verbose);
250 pip.push_back(
new OpMass(
"U",
"U", get_rho));
251 static_cast<OpMass &
>(pip.back()).feScalingFun =
252 [](
const FEMethod *fe_ptr) {
return fe_ptr->ts_aa; };
256 auto apply_rhs = [&](
auto &pip) {
260 pipeline_mng->getOpDomainRhsPipeline(), {H1},
"GEOMETRY");
262 mField, pip,
"U",
"MAT_ELASTIC", Sev::inform);
263 CHKERR add_rho_block(pip,
"MAT_RHO", Sev::inform);
265 auto mat_acceleration_ptr = boost::make_shared<MatrixDouble>();
268 "U", mat_acceleration_ptr));
269 pip.push_back(
new OpInertiaForce(
"U", mat_acceleration_ptr, get_rho));
273 {boost::make_shared<TimeScaleVector<SPACE_DIM>>(
"-time_vector_file",
275 "BODY_FORCE", Sev::inform);
280 CHKERR apply_lhs(pipeline_mng->getOpDomainLhsPipeline());
281 CHKERR apply_rhs(pipeline_mng->getOpDomainRhsPipeline());
289 auto get_bc_hook = [&]() {
291 mField, pipeline_mng->getDomainRhsFE(),
292 {boost::make_shared<TimeScale>()});
296 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook();
317 "out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
324 boost::shared_ptr<PostProcEle>
postProc;
333 auto dm =
simple->getDM();
335 ts = pipeline_mng->createTSIM2();
338 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
340 post_proc_fe->getOpPtrVector(), {H1});
342 mField, post_proc_fe->getOpPtrVector(),
"U",
"MAT_ELASTIC", Sev::inform);
344 auto u_ptr = boost::make_shared<MatrixDouble>();
345 post_proc_fe->getOpPtrVector().push_back(
347 auto X_ptr = boost::make_shared<MatrixDouble>();
348 post_proc_fe->getOpPtrVector().push_back(
353 post_proc_fe->getOpPtrVector().push_back(
357 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
361 {{
"U", u_ptr}, {
"GEOMETRY", X_ptr}},
363 {{
"GRAD", common_ptr->matGradPtr},
364 {
"FIRST_PIOLA", common_ptr->getMatFirstPiolaStress()}},
373 boost::shared_ptr<FEMethod> null_fe;
374 auto monitor_ptr = boost::make_shared<Monitor>(dm, post_proc_fe);
376 null_fe, monitor_ptr);
380 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
386 CHKERR TS2SetSolution(ts, T, TT);
388 CHKERR TSSetI2Jacobian(ts,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
389 CHKERR TSSetFromOptions(ts);
392 CHKERR TSGetTime(ts, &ftime);
394 PetscInt steps, snesfails, rejects, nonlinits, linits;
395#if PETSC_VERSION_GE(3, 8, 0)
396 CHKERR TSGetStepNumber(ts, &steps);
398 CHKERR TSGetTimeStepNumber(ts, &steps);
400 CHKERR TSGetSNESFailures(ts, &snesfails);
401 CHKERR TSGetStepRejections(ts, &rejects);
402 CHKERR TSGetSNESIterations(ts, &nonlinits);
403 CHKERR TSGetKSPIterations(ts, &linits);
405 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
407 steps, rejects, snesfails, ftime, nonlinits, linits);
416 PetscBool test_flg = PETSC_FALSE;
424 CHKERR VecNorm(T, NORM_2, &nrm2);
425 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Regression norm " << nrm2;
426 constexpr double regression_value = 0.0194561;
427 if (fabs(nrm2 - regression_value) > 1e-2)
429 "Regression test failed; wrong norm value.");
444int main(
int argc,
char *argv[]) {
447 const char param_file[] =
"param_file.petsc";
451 auto core_log = logging::core::get();
460 DMType dm_name =
"DMMOFEM";
465 moab::Core mb_instance;
466 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)
constexpr int SPACE_DIM
[Define dimension]
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
@ 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 >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
constexpr double omega
Save field DOFS on vertices/tags.
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.
auto createDMMatrix(DM dm)
Get smart matrix 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
FTensor::Index< 'i', SPACE_DIM > i
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
auto commonDataFactory(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
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]
constexpr double poisson_ratio
constexpr double young_modulus
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
static constexpr int approx_order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
FTensor::Index< 'm', 3 > m
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 field 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.
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode getOptions()
get options
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference 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