v0.15.0
Loading...
Searching...
No Matches
nonlinear_dynamic_elastic.cpp
Go to the documentation of this file.
1/**
2 * \file dynamic_elastic.cpp
3 * \example dynamic_elastic.cpp
4 *
5 * Plane stress elastic dynamic problem
6 *
7 */
8
9#include <MoFEM.hpp>
10#include <MatrixFunction.hpp>
11using namespace MoFEM;
12
13template <int DIM> struct ElementsAndOps {};
14
15template <> struct ElementsAndOps<2> {
17};
18
19template <> struct ElementsAndOps<3> {
21};
22
23constexpr int SPACE_DIM =
24 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
25
28using DomainEleOp = DomainEle::UserDataOperator;
30
35
40 SPACE_DIM>;
41
42constexpr double rho = 1;
43constexpr double omega = 2.4;
44constexpr double young_modulus = 1;
45constexpr double poisson_ratio = 0.25;
46
47#include <HenckyOps.hpp>
48using namespace HenckyOps;
49
50struct Example {
51
52 Example(MoFEM::Interface &m_field) : mField(m_field) {}
53
55
56private:
58
66};
67
68//! [Run problem]
79}
80//! [Run problem]
81
82//! [Read mesh]
87 CHKERR simple->loadFile();
89}
90//! [Read mesh]
91
92//! [Set up problem]
96 // Add field
98 SPACE_DIM);
99 CHKERR simple->addDataField("GEOMETRY", H1, AINSWORTH_LEGENDRE_BASE,
100 SPACE_DIM);
101 int order = 2;
102 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
103 CHKERR simple->setFieldOrder("U", order);
104 CHKERR simple->setFieldOrder("GEOMETRY", order);
105 CHKERR simple->setUp();
106
107 auto project_ho_geometry = [&]() {
108 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
109 return mField.loop_dofs("GEOMETRY", ent_method);
110 };
111 CHKERR project_ho_geometry();
112
114}
115//! [Set up problem]
116
117//! [Boundary condition]
120
122 auto bc_mng = mField.getInterface<BcManager>();
123
124 CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
125 simple->getProblemName(), "U");
126
128}
129//! [Boundary condition]
130
131//! [Push operators to pipeline]
134 auto *pipeline_mng = mField.getInterface<PipelineManager>();
135
136 auto get_body_force = [this](const double, const double, const double) {
139 t_source(i) = 0.;
140 t_source(0) = 0.1;
141 t_source(1) = 1.;
142 return t_source;
143 };
144
145 auto rho_ptr = boost::make_shared<double>(rho);
146
147 auto add_rho_block = [this, rho_ptr](auto &pip, auto block_name, Sev sev) {
149
150 struct OpMatRhoBlocks : public DomainEleOp {
151 OpMatRhoBlocks(boost::shared_ptr<double> rho_ptr,
152 MoFEM::Interface &m_field, Sev sev,
153 std::vector<const CubitMeshSets *> meshset_vec_ptr)
154 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), rhoPtr(rho_ptr) {
155 CHK_THROW_MESSAGE(extractRhoData(m_field, meshset_vec_ptr, sev),
156 "Can not get data from block");
157 }
158
159 MoFEMErrorCode doWork(int side, EntityType type,
161
163
164 for (auto &b : blockData) {
165 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
166 *rhoPtr = b.rho;
168 }
169 }
170
171 *rhoPtr = rho;
172
174 }
175
176 private:
177 struct BlockData {
178 double rho;
179 Range blockEnts;
180 };
181
182 std::vector<BlockData> blockData;
183
185 extractRhoData(MoFEM::Interface &m_field,
186 std::vector<const CubitMeshSets *> meshset_vec_ptr,
187 Sev sev) {
189
190 for (auto m : meshset_vec_ptr) {
191 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Rho Block") << *m;
192 std::vector<double> block_data;
193 CHKERR m->getAttributes(block_data);
194 if (block_data.size() < 1) {
195 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
196 "Expected that block has two attributes");
197 }
198 auto get_block_ents = [&]() {
199 Range ents;
200 CHKERR
201 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
202 return ents;
203 };
204
205 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Rho Block")
206 << m->getName() << ": ro = " << block_data[0];
207
208 blockData.push_back({block_data[0], get_block_ents()});
209 }
210 MOFEM_LOG_CHANNEL("WORLD");
211 MOFEM_LOG_CHANNEL("WORLD");
213 }
214
215 boost::shared_ptr<double> rhoPtr;
216 };
217
218 pip.push_back(new OpMatRhoBlocks(
219 rho_ptr, mField, sev,
220
221 // Get blockset using regular expression
223
224 (boost::format("%s(.*)") % block_name).str()
225
226 ))
227
228 ));
230 };
231
232 // Get pointer to U_tt shift in domain element
233 auto get_rho = [rho_ptr](const double, const double, const double) {
234 return *rho_ptr;
235 };
236
237 // specific time scaling
238 auto get_time_scale = [this](const double time) {
239 return sin(time * omega * M_PI);
240 };
241
242 auto apply_lhs = [&](auto &pip) {
245 "GEOMETRY");
247 mField, pip, "U", "MAT_ELASTIC", Sev::verbose);
248 CHKERR add_rho_block(pip, "MAT_RHO", Sev::verbose);
249
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; };
254 };
255
256 auto apply_rhs = [&](auto &pip) {
258
260 pipeline_mng->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
262 mField, pip, "U", "MAT_ELASTIC", Sev::inform);
263 CHKERR add_rho_block(pip, "MAT_RHO", Sev::inform);
264
265 auto mat_acceleration_ptr = boost::make_shared<MatrixDouble>();
266 // Apply inertia
268 "U", mat_acceleration_ptr));
269 pip.push_back(new OpInertiaForce("U", mat_acceleration_ptr, get_rho));
270
272 pip, mField, "U", {},
273 {boost::make_shared<TimeScaleVector<SPACE_DIM>>("-time_vector_file",
274 true)},
275 "BODY_FORCE", Sev::inform);
276
278 };
279
280 CHKERR apply_lhs(pipeline_mng->getOpDomainLhsPipeline());
281 CHKERR apply_rhs(pipeline_mng->getOpDomainRhsPipeline());
282
283 auto integration_rule = [](int, int, int approx_order) {
284 return 2 * approx_order;
285 };
286 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
287 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
288
289 auto get_bc_hook = [&]() {
291 mField, pipeline_mng->getDomainRhsFE(),
292 {boost::make_shared<TimeScale>()});
293 return hook;
294 };
295
296 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook();
297
299}
300//! [Push operators to pipeline]
301
302/**
303 * @brief Monitor solution
304 *
305 * This functions is called by TS solver at the end of each step. It is used
306 * to output results to the hard drive.
307 */
308struct Monitor : public FEMethod {
309 Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
310 : dM(dm), postProc(post_proc){};
313 constexpr int save_every_nth_step = 1;
314 if (ts_step % save_every_nth_step == 0) {
316 CHKERR postProc->writeFile(
317 "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
318 }
320 }
321
322private:
324 boost::shared_ptr<PostProcEle> postProc;
325};
326
327//! [Solve]
330 auto *simple = mField.getInterface<Simple>();
331 auto *pipeline_mng = mField.getInterface<PipelineManager>();
332
333 auto dm = simple->getDM();
335 ts = pipeline_mng->createTSIM2();
336
337 // Setup postprocessing
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);
343
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(
349 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
350
352
353 post_proc_fe->getOpPtrVector().push_back(
354
355 new OpPPMap(
356
357 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
358
359 {},
360
361 {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
362
363 {{"GRAD", common_ptr->matGradPtr},
364 {"FIRST_PIOLA", common_ptr->getMatFirstPiolaStress()}},
365
366 {}
367
368 )
369
370 );
371
372 // Add monitor to time solver
373 boost::shared_ptr<FEMethod> null_fe;
374 auto monitor_ptr = boost::make_shared<Monitor>(dm, post_proc_fe);
375 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
376 null_fe, monitor_ptr);
377
378 double ftime = 1;
379 // CHKERR TSSetMaxTime(ts, ftime);
380 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
381
382 auto T = createDMVector(simple->getDM());
383 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
384 SCATTER_FORWARD);
385 auto TT = vectorDuplicate(T);
386 CHKERR TS2SetSolution(ts, T, TT);
387 auto B = createDMMatrix(dm);
388 CHKERR TSSetI2Jacobian(ts, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
389 CHKERR TSSetFromOptions(ts);
390
391 CHKERR TSSolve(ts, NULL);
392 CHKERR TSGetTime(ts, &ftime);
393
394 PetscInt steps, snesfails, rejects, nonlinits, linits;
395#if PETSC_VERSION_GE(3, 8, 0)
396 CHKERR TSGetStepNumber(ts, &steps);
397#else
398 CHKERR TSGetTimeStepNumber(ts, &steps);
399#endif
400 CHKERR TSGetSNESFailures(ts, &snesfails);
401 CHKERR TSGetStepRejections(ts, &rejects);
402 CHKERR TSGetSNESIterations(ts, &nonlinits);
403 CHKERR TSGetKSPIterations(ts, &linits);
404 MOFEM_LOG_C("EXAMPLE", Sev::inform,
405 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
406 "%d, linits %d\n",
407 steps, rejects, snesfails, ftime, nonlinits, linits);
408
410}
411//! [Solve]
412
413//! [Postprocess results]
416 PetscBool test_flg = PETSC_FALSE;
417 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test", &test_flg, PETSC_NULLPTR);
418 if (test_flg) {
419 auto *simple = mField.getInterface<Simple>();
420 auto T = createDMVector(simple->getDM());
421 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
422 SCATTER_FORWARD);
423 double nrm2;
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)
428 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
429 "Regression test failed; wrong norm value.");
430 }
432}
433//! [Postprocess results]
434
435//! [Check]
439}
440//! [Check]
441
442static char help[] = "...\n\n";
443
444int main(int argc, char *argv[]) {
445
446 // Initialisation of MoFEM/PETSc and MOAB data structures
447 const char param_file[] = "param_file.petsc";
448 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
449
450 // Add logging channel for example
451 auto core_log = logging::core::get();
452 core_log->add_sink(
454 LogManager::setLog("EXAMPLE");
455 MOFEM_LOG_TAG("EXAMPLE", "example");
456
457 try {
458
459 //! [Register MoFEM discrete manager in PETSc]
460 DMType dm_name = "DMMOFEM";
461 CHKERR DMRegister_MoFEM(dm_name);
462 //! [Register MoFEM discrete manager in PETSc
463
464 //! [Create MoAB]
465 moab::Core mb_instance; ///< mesh database
466 moab::Interface &moab = mb_instance; ///< mesh database interface
467 //! [Create MoAB]
468
469 //! [Create MoFEM]
470 MoFEM::Core core(moab); ///< finite element database
471 MoFEM::Interface &m_field = core; ///< finite element database interface
472 //! [Create MoFEM]
473
474 //! [Example]
475 Example ex(m_field);
476 CHKERR ex.runProblem();
477 //! [Example]
478 }
480
482}
#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)
Definition acoustic.cpp:69
constexpr int SPACE_DIM
[Define dimension]
Definition adjoint.cpp:22
int main()
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
#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()
@ H1
continuous field
Definition definitions.h:85
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
constexpr double omega
Save field DOFS on vertices/tags.
auto integration_rule
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1059
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
Definition Common.hpp:10
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.
Definition DMMoFEM.cpp:1046
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.
constexpr double rho
static char help[]
[Check]
constexpr int SPACE_DIM
constexpr double poisson_ratio
constexpr double omega
constexpr double young_modulus
int save_every_nth_step
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double rho
Definition plastic.cpp:144
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13
static constexpr int approx_order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition seepage.cpp:55
FTensor::Index< 'm', 3 > m
[Example]
Definition plastic.cpp:216
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition plastic.cpp:226
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition BCData.hpp:72
Data on single entity (This is passed as argument to DataOperator::doWork)
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
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.
Assembly methods.
Definition Natural.hpp:65
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.
Definition Simple.hpp:27
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.
Definition Simple.cpp:261
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
[Push operators to pipeline]
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
SmartPetscObj< DM > dM
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcEle > postProc