v0.13.1
Loading...
Searching...
No Matches
nonlinear_elastic.cpp
Go to the documentation of this file.
1/**
2 * \file nonlinear_elastic.cpp
3 * \example nonlinear_elastic.cpp
4 *
5 * Plane stress elastic dynamic problem
6 *
7 */
8
9#include <MoFEM.hpp>
10#include <MatrixFunction.hpp>
11
12using namespace MoFEM;
13
14template <int DIM> struct ElementsAndOps {};
15
16template <> struct ElementsAndOps<2> {
19};
20
21template <> struct ElementsAndOps<3> {
24};
25
26constexpr int SPACE_DIM = 3; //< Space dimension of problem, mesh
27
33
35
37 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
40
45
49
50constexpr double young_modulus = 100;
51constexpr double poisson_ratio = 0.3;
52constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
53constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
54
55#include <HenckyOps.hpp>
56using namespace HenckyOps;
57
58struct Example {
59
60 Example(MoFEM::Interface &m_field) : mField(m_field) {}
61
63
64private:
66
75
76 boost::shared_ptr<MatrixDouble> matGradPtr;
77 boost::shared_ptr<MatrixDouble> matStrainPtr;
78 boost::shared_ptr<MatrixDouble> matStressPtr;
79 boost::shared_ptr<MatrixDouble> matDPtr;
80 boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
81};
82
83//! [Create common data]
86
87 auto set_material_stiffness = [&]() {
94 constexpr double A =
95 (SPACE_DIM == 2) ? 2 * shear_modulus_G /
96 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
97 : 1;
98 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
99 t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
100 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
101 t_kd(i, j) * t_kd(k, l);
103 };
104
105 matGradPtr = boost::make_shared<MatrixDouble>();
106 matStrainPtr = boost::make_shared<MatrixDouble>();
107 matStressPtr = boost::make_shared<MatrixDouble>();
108 matDPtr = boost::make_shared<MatrixDouble>();
109
110 commonHenckyDataPtr = boost::make_shared<HenckyOps::CommonData>();
111 commonHenckyDataPtr->matGradPtr = matGradPtr;
112 commonHenckyDataPtr->matDPtr = matDPtr;
113
114 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
115 matDPtr->resize(size_symm * size_symm, 1);
116
117 CHKERR set_material_stiffness();
118
120}
121//! [Create common data]
122
123//! [Run problem]
135}
136//! [Run problem]
137
138//! [Read mesh]
142 CHKERR simple->getOptions();
143 CHKERR simple->loadFile();
145}
146//! [Read mesh]
147
148//! [Set up problem]
152 // Add field
153 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
154 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
155 int order = 2;
156 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
157 CHKERR simple->setFieldOrder("U", order);
158 CHKERR simple->setUp();
160}
161//! [Set up problem]
162
163//! [Boundary condition]
166 auto *pipeline_mng = mField.getInterface<PipelineManager>();
168 auto bc_mng = mField.getInterface<BcManager>();
169
171 pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U",
172 {boost::make_shared<TimeScale>()}, "FORCE", Sev::inform);
173
174 //! [Define gravity vector]
176 pipeline_mng->getOpDomainRhsPipeline(), mField, "U",
177 {boost::make_shared<TimeScale>()}, "BODY_FORCE", Sev::inform);
178
179 // Essential BC
180 CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
181 simple->getProblemName(), "U");
182
183 auto get_bc_hook = [&]() {
185 mField, pipeline_mng->getDomainRhsFE(),
186 {boost::make_shared<TimeScale>()});
187 return hook;
188 };
189
190 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook();
191
193}
194//! [Boundary condition]
195
196//! [Push operators to pipeline]
199 auto *simple = mField.getInterface<Simple>();
200 auto *pipeline_mng = mField.getInterface<PipelineManager>();
201
202 auto add_domain_base_ops = [&](auto &pipeline) {
204
205 auto det_ptr = boost::make_shared<VectorDouble>();
206 auto jac_ptr = boost::make_shared<MatrixDouble>();
207 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
208 pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
209 pipeline.push_back(
210 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
211 pipeline.push_back(
212 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
213 pipeline.push_back(new OpSetHOWeights(det_ptr));
214
216 "U", matGradPtr));
217 pipeline.push_back(
219 pipeline.push_back(
221 pipeline.push_back(
223 pipeline.push_back(
225 pipeline.push_back(
227
229 };
230
231 auto add_domain_ops_lhs = [&](auto &pipeline) {
233
234 pipeline.push_back(
236 pipeline.push_back(new OpK("U", "U", commonHenckyDataPtr->getMatTangent()));
237
239 };
240
241 auto add_domain_ops_rhs = [&](auto &pipeline) {
243 // Calculate internal forece
244 pipeline.push_back(new OpInternalForce(
245 "U", commonHenckyDataPtr->getMatFirstPiolaStress()));
246
248 };
249
250 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
251 CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
252 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
253 CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
254
255 auto integration_rule = [](int, int, int approx_order) {
256 return 2 * (approx_order - 1);
257 };
258
259 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
260 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
261 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
262
264}
265//! [Push operators to pipeline]
266
267/**
268 * @brief Monitor solution
269 *
270 * This functions is called by TS solver at the end of each step. It is used
271 * to output results to the hard drive.
272 */
273struct Monitor : public FEMethod {
274 Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
275 : dM(dm), postProc(post_proc){};
278 constexpr int save_every_nth_step = 1;
279 if (ts_step % save_every_nth_step == 0) {
281 CHKERR postProc->writeFile(
282 "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
283 }
285 }
286
287private:
289 boost::shared_ptr<PostProcEle> postProc;
290};
291
292//! [Solve]
295 auto *simple = mField.getInterface<Simple>();
296 auto *pipeline_mng = mField.getInterface<PipelineManager>();
297
298 auto dm = simple->getDM();
299 auto ts = pipeline_mng->createTSIM();
300
301 // Setup postprocessing
302 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
303
304 auto det_ptr = boost::make_shared<VectorDouble>();
305 auto jac_ptr = boost::make_shared<MatrixDouble>();
306 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
307 post_proc_fe->getOpPtrVector().push_back(
308 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
309 post_proc_fe->getOpPtrVector().push_back(
310 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
311 post_proc_fe->getOpPtrVector().push_back(
312 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
313
314 post_proc_fe->getOpPtrVector().push_back(
316 "U", commonHenckyDataPtr->matGradPtr));
317 post_proc_fe->getOpPtrVector().push_back(
319 post_proc_fe->getOpPtrVector().push_back(
321 post_proc_fe->getOpPtrVector().push_back(
323 post_proc_fe->getOpPtrVector().push_back(
325 post_proc_fe->getOpPtrVector().push_back(
327
328 auto u_ptr = boost::make_shared<MatrixDouble>();
329 post_proc_fe->getOpPtrVector().push_back(
331
333
334 post_proc_fe->getOpPtrVector().push_back(
335
336 new OpPPMap(
337
338 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
339
340 {},
341
342 {{"U", u_ptr}},
343
344 {{"GRAD", commonHenckyDataPtr->matGradPtr},
345 {"FIRST_PIOLA", commonHenckyDataPtr->getMatFirstPiolaStress()}},
346
347 {}
348
349 )
350
351 );
352
353 // Add monitor to time solver
354 boost::shared_ptr<FEMethod> null_fe;
355 auto monitor_ptr = boost::make_shared<Monitor>(dm, post_proc_fe);
356 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
357 null_fe, monitor_ptr);
358
359 double ftime = 1;
360 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
361 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
362
363 auto D = smartCreateDMVector(simple->getDM());
364
365 CHKERR TSSetSolution(ts, D);
366 CHKERR TSSetFromOptions(ts);
367
368 CHKERR TSSolve(ts, NULL);
369 CHKERR TSGetTime(ts, &ftime);
370
371 PetscInt steps, snesfails, rejects, nonlinits, linits;
372 CHKERR TSGetTimeStepNumber(ts, &steps);
373 CHKERR TSGetSNESFailures(ts, &snesfails);
374 CHKERR TSGetStepRejections(ts, &rejects);
375 CHKERR TSGetSNESIterations(ts, &nonlinits);
376 CHKERR TSGetKSPIterations(ts, &linits);
377 MOFEM_LOG_C("EXAMPLE", Sev::inform,
378 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
379 "%D, linits %D\n",
380 steps, rejects, snesfails, ftime, nonlinits, linits);
381
383}
384//! [Solve]
385
386//! [Postprocess results]
389 PetscBool test_flg = PETSC_FALSE;
390 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
391 if (test_flg) {
392 auto *simple = mField.getInterface<Simple>();
393 auto T = smartCreateDMVector(simple->getDM());
394 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
395 SCATTER_FORWARD);
396 double nrm2;
397 CHKERR VecNorm(T, NORM_2, &nrm2);
398 MOFEM_LOG("EXAMPLE", Sev::inform) << "Regression norm " << nrm2;
399 constexpr double regression_value = 1.09572;
400 if (fabs(nrm2 - regression_value) > 1e-2)
401 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
402 "Regression test field; wrong norm value.");
403 }
405}
406//! [Postprocess results]
407
408//! [Check]
412}
413//! [Check]
414
415static char help[] = "...\n\n";
416
417int main(int argc, char *argv[]) {
418
419 // Initialisation of MoFEM/PETSc and MOAB data structures
420 const char param_file[] = "param_file.petsc";
422
423 // Add logging channel for example
424 auto core_log = logging::core::get();
425 core_log->add_sink(
427 LogManager::setLog("EXAMPLE");
428 MOFEM_LOG_TAG("EXAMPLE", "example");
429
430 try {
431
432 //! [Register MoFEM discrete manager in PETSc]
433 DMType dm_name = "DMMOFEM";
434 CHKERR DMRegister_MoFEM(dm_name);
435 //! [Register MoFEM discrete manager in PETSc
436
437 //! [Create MoAB]
438 moab::Core mb_instance; ///< mesh database
439 moab::Interface &moab = mb_instance; ///< mesh database interface
440 //! [Create MoAB]
441
442 //! [Create MoFEM]
443 MoFEM::Core core(moab); ///< finite element database
444 MoFEM::Interface &m_field = core; ///< finite element database interface
445 //! [Create MoFEM]
446
447 //! [Example]
448 Example ex(m_field);
449 CHKERR ex.runProblem();
450 //! [Example]
451 }
453
455}
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
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
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#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
#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
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
Definition: elastic.cpp:44
constexpr double shear_modulus_G
Definition: elastic.cpp:63
constexpr double bulk_modulus_K
Definition: elastic.cpp:62
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
auto integration_rule
constexpr auto t_kd
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.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: DMMMoFEM.cpp:533
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
double D
FTensor::Index< 'k', 3 > k
const double T
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
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: DMMMoFEM.cpp:1003
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)
int save_every_nth_step
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto size_symm
Definition: plastic.cpp:33
constexpr AssemblyType A
Definition: plastic.cpp:35
static constexpr int approx_order
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
[Example]
Definition: plastic.cpp:137
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
boost::shared_ptr< MatrixDouble > matGradPtr
boost::shared_ptr< MatrixDouble > matStressPtr
boost::shared_ptr< MatrixDouble > matStrainPtr
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:144
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
boost::shared_ptr< MatrixDouble > matDPtr
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: plastic.cpp:153
Simple interface for fast problem set-up.
Definition: BcManager.hpp:23
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:72
Data on single entity (This is passed as argument to DataOperator::doWork)
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:39
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:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
Assembly methods.
Definition: Natural.hpp:57
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Set inverse jacobian to base functions.
PipelineManager interface.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.
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
static char help[]
[Check]
constexpr int SPACE_DIM
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
constexpr double young_modulus