v0.14.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
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
50static double *ts_time_ptr;
51static double *ts_aa_ptr;
52
53struct Example {
54
55 Example(MoFEM::Interface &m_field) : mField(m_field) {}
56
58
59private:
61
69};
70
71//! [Run problem]
82}
83//! [Run problem]
84
85//! [Read mesh]
89 CHKERR simple->getOptions();
90 CHKERR simple->loadFile();
92}
93//! [Read mesh]
94
95//! [Set up problem]
99 // Add field
100 CHKERR simple->addDomainField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
101 SPACE_DIM);
102 CHKERR simple->addDataField("GEOMETRY", H1, AINSWORTH_LEGENDRE_BASE,
103 SPACE_DIM);
104 int order = 2;
105 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
106 CHKERR simple->setFieldOrder("U", order);
107 CHKERR simple->setFieldOrder("GEOMETRY", order);
108 CHKERR simple->setUp();
109
110 auto project_ho_geometry = [&]() {
111 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
112 return mField.loop_dofs("GEOMETRY", ent_method);
113 };
114 CHKERR project_ho_geometry();
115
117}
118//! [Set up problem]
119
120//! [Boundary condition]
123
125 auto bc_mng = mField.getInterface<BcManager>();
126
127 CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
128 simple->getProblemName(), "U");
129
131}
132//! [Boundary condition]
133
134//! [Push operators to pipeline]
137 auto *simple = mField.getInterface<Simple>();
138 auto *pipeline_mng = mField.getInterface<PipelineManager>();
139
140 auto get_body_force = [this](const double, const double, const double) {
143 t_source(i) = 0.;
144 t_source(0) = 0.1;
145 t_source(1) = 1.;
146 return t_source;
147 };
148
149 auto rho_ptr = boost::make_shared<double>(rho);
150
151 auto add_rho_block = [this, rho_ptr](auto &pip, auto block_name, Sev sev) {
153
154 struct OpMatRhoBlocks : public DomainEleOp {
155 OpMatRhoBlocks(boost::shared_ptr<double> rho_ptr,
156 MoFEM::Interface &m_field, Sev sev,
157 std::vector<const CubitMeshSets *> meshset_vec_ptr)
158 : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), rhoPtr(rho_ptr) {
159 CHK_THROW_MESSAGE(extractRhoData(m_field, meshset_vec_ptr, sev),
160 "Can not get data from block");
161 }
162
163 MoFEMErrorCode doWork(int side, EntityType type,
165
167
168 for (auto &b : blockData) {
169 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
170 *rhoPtr = b.rho;
172 }
173 }
174
175 *rhoPtr = rho;
176
178 }
179
180 private:
181 struct BlockData {
182 double rho;
183 Range blockEnts;
184 };
185
186 std::vector<BlockData> blockData;
187
189 extractRhoData(MoFEM::Interface &m_field,
190 std::vector<const CubitMeshSets *> meshset_vec_ptr,
191 Sev sev) {
193
194 for (auto m : meshset_vec_ptr) {
195 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Rho Block") << *m;
196 std::vector<double> block_data;
197 CHKERR m->getAttributes(block_data);
198 if (block_data.size() < 1) {
199 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
200 "Expected that block has two attributes");
201 }
202 auto get_block_ents = [&]() {
203 Range ents;
204 CHKERR
205 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
206 return ents;
207 };
208
209 MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Rho Block")
210 << m->getName() << ": ro = " << block_data[0];
211
212 blockData.push_back({block_data[0], get_block_ents()});
213 }
214 MOFEM_LOG_CHANNEL("WORLD");
215 MOFEM_LOG_CHANNEL("WORLD");
217 }
218
219 boost::shared_ptr<double> rhoPtr;
220 };
221
222 pip.push_back(new OpMatRhoBlocks(
223 rho_ptr, mField, sev,
224
225 // Get blockset using regular expression
227
228 (boost::format("%s(.*)") % block_name).str()
229
230 ))
231
232 ));
234 };
235
236 // Get pointer to U_tt shift in domain element
237 auto get_rho = [rho_ptr](const double, const double, const double) {
238 return *rho_ptr;
239 };
240
241 // specific time scaling
242 auto get_time_scale = [this](const double time) {
243 return sin(time * omega * M_PI);
244 };
245
246 auto apply_lhs = [&](auto &pip) {
249 "GEOMETRY");
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);
253
254 pip.push_back(new OpMass("U", "U", get_rho));
255 static_cast<OpMass &>(pip.back()).feScalingFun =
256 [](const FEMethod *fe_ptr) { return fe_ptr->ts_aa; };
258 };
259
260 auto apply_rhs = [&](auto &pip) {
262
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);
268
269 auto mat_acceleration_ptr = boost::make_shared<MatrixDouble>();
270 // Apply inertia
272 "U", mat_acceleration_ptr));
273 pip.push_back(new OpInertiaForce("U", mat_acceleration_ptr, get_rho));
274
276 pip, mField, "U", {},
277 {boost::make_shared<TimeScaleVector<SPACE_DIM>>("-time_vector_file",
278 true)},
279 "BODY_FORCE", Sev::inform);
280
282 };
283
284 CHKERR apply_lhs(pipeline_mng->getOpDomainLhsPipeline());
285 CHKERR apply_rhs(pipeline_mng->getOpDomainRhsPipeline());
286
287 auto integration_rule = [](int, int, int approx_order) {
288 return 2 * approx_order;
289 };
290 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
291 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
292
293 auto get_bc_hook = [&]() {
295 mField, pipeline_mng->getDomainRhsFE(),
296 {boost::make_shared<TimeScale>()});
297 return hook;
298 };
299
300 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook();
301
303}
304//! [Push operators to pipeline]
305
306/**
307 * @brief Monitor solution
308 *
309 * This functions is called by TS solver at the end of each step. It is used
310 * to output results to the hard drive.
311 */
312struct Monitor : public FEMethod {
313 Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
314 : dM(dm), postProc(post_proc){};
317 constexpr int save_every_nth_step = 1;
318 if (ts_step % save_every_nth_step == 0) {
320 CHKERR postProc->writeFile(
321 "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
322 }
324 }
325
326private:
328 boost::shared_ptr<PostProcEle> postProc;
329};
330
331//! [Solve]
334 auto *simple = mField.getInterface<Simple>();
335 auto *pipeline_mng = mField.getInterface<PipelineManager>();
336
337 auto dm = simple->getDM();
339 ts = pipeline_mng->createTSIM2();
340
341 // Setup postprocessing
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);
347
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(
353 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
354
356
357 post_proc_fe->getOpPtrVector().push_back(
358
359 new OpPPMap(
360
361 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
362
363 {},
364
365 {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
366
367 {{"GRAD", common_ptr->matGradPtr},
368 {"FIRST_PIOLA", common_ptr->getMatFirstPiolaStress()}},
369
370 {}
371
372 )
373
374 );
375
376 // Add monitor to time solver
377 boost::shared_ptr<FEMethod> null_fe;
378 auto monitor_ptr = boost::make_shared<Monitor>(dm, post_proc_fe);
379 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
380 null_fe, monitor_ptr);
381
382 double ftime = 1;
383 // CHKERR TSSetMaxTime(ts, ftime);
384 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
385
386 auto T = createDMVector(simple->getDM());
387 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
388 SCATTER_FORWARD);
389 auto TT = vectorDuplicate(T);
390 CHKERR TS2SetSolution(ts, T, TT);
391 CHKERR TSSetFromOptions(ts);
392
393 CHKERR TSSolve(ts, NULL);
394 CHKERR TSGetTime(ts, &ftime);
395
396 PetscInt steps, snesfails, rejects, nonlinits, linits;
397#if PETSC_VERSION_GE(3, 8, 0)
398 CHKERR TSGetStepNumber(ts, &steps);
399#else
400 CHKERR TSGetTimeStepNumber(ts, &steps);
401#endif
402 CHKERR TSGetSNESFailures(ts, &snesfails);
403 CHKERR TSGetStepRejections(ts, &rejects);
404 CHKERR TSGetSNESIterations(ts, &nonlinits);
405 CHKERR TSGetKSPIterations(ts, &linits);
406 MOFEM_LOG_C("EXAMPLE", Sev::inform,
407 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
408 "%d, linits %d\n",
409 steps, rejects, snesfails, ftime, nonlinits, linits);
410
412}
413//! [Solve]
414
415//! [Postprocess results]
418 PetscBool test_flg = PETSC_FALSE;
419 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
420 if (test_flg) {
421 auto *simple = mField.getInterface<Simple>();
422 auto T = createDMVector(simple->getDM());
423 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
424 SCATTER_FORWARD);
425 double nrm2;
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)
430 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
431 "Regression test failed; wrong norm value.");
432 }
434}
435//! [Postprocess results]
436
437//! [Check]
441}
442//! [Check]
443
444static char help[] = "...\n\n";
445
446int main(int argc, char *argv[]) {
447
448 // Initialisation of MoFEM/PETSc and MOAB data structures
449 const char param_file[] = "param_file.petsc";
451
452 // Add logging channel for example
453 auto core_log = logging::core::get();
454 core_log->add_sink(
456 LogManager::setLog("EXAMPLE");
457 MOFEM_LOG_TAG("EXAMPLE", "example");
458
459 try {
460
461 //! [Register MoFEM discrete manager in PETSc]
462 DMType dm_name = "DMMOFEM";
463 CHKERR DMRegister_MoFEM(dm_name);
464 //! [Register MoFEM discrete manager in PETSc
465
466 //! [Create MoAB]
467 moab::Core mb_instance; ///< mesh database
468 moab::Interface &moab = mb_instance; ///< mesh database interface
469 //! [Create MoAB]
470
471 //! [Create MoFEM]
472 MoFEM::Core core(moab); ///< finite element database
473 MoFEM::Interface &m_field = core; ///< finite element database interface
474 //! [Create MoFEM]
475
476 //! [Example]
477 Example ex(m_field);
478 CHKERR ex.runProblem();
479 //! [Example]
480 }
482
484}
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ 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 ...
Definition: definitions.h:346
@ 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()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
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:509
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
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:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
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
const double T
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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:1042
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]
static double * ts_aa_ptr
constexpr int SPACE_DIM
constexpr double poisson_ratio
constexpr double omega
static double * ts_time_ptr
constexpr double young_modulus
int save_every_nth_step
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double rho
Definition: plastic.cpp:191
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
static constexpr int approx_order
constexpr double omega
FTensor::Index< 'i', 3 > i
[Example]
Definition: plastic.cpp:228
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:235
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
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:112
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
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.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
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 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.
Definition: Simple.hpp:27
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.
[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< PostProcFace > postProc