v0.10.0
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 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #include <MoFEM.hpp>
24 
25 using namespace MoFEM;
26 
27 template <int DIM> struct ElementsAndOps {};
28 
29 template <> struct ElementsAndOps<2> {
33 };
34 
35 template <> struct ElementsAndOps<3> {
39 };
40 
41 constexpr int SPACE_DIM = 2; //< Space dimension of problem, mesh
42 
47 
49  GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
55  GAUSS>::OpSource<1, SPACE_DIM>;
58 
59 constexpr double rho = 1;
60 constexpr double omega = 2.4;
61 constexpr double young_modulus = 1;
62 constexpr double poisson_ratio = 0.25;
63 constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
64 constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
65 
66 #include <OpPostProcElastic.hpp>
67 using namespace Tutorial;
68 
69 static double *ts_time_ptr;
70 static double *ts_aa_ptr;
71 
72 struct Example {
73 
74  Example(MoFEM::Interface &m_field) : mField(m_field) {}
75 
77 
78 private:
79  MoFEM::Interface &mField;
80 
89 
91  boost::shared_ptr<MatrixDouble> matGradPtr;
92  boost::shared_ptr<MatrixDouble> matStrainPtr;
93  boost::shared_ptr<MatrixDouble> matStressPtr;
94  boost::shared_ptr<MatrixDouble> matAccelerationPtr;
95  boost::shared_ptr<MatrixDouble> matInertiaPtr;
96  boost::shared_ptr<MatrixDouble> matDPtr;
97 };
98 
99 //! [Create common data]
102 
103  auto set_matrial_stiffens = [&]() {
108  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
110  constexpr double A =
111  (SPACE_DIM == 2) ? 2 * shear_modulus_G /
112  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
113  : 1;
114  auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
115  t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
116  A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
117  t_kd(i, j) * t_kd(k, l);
119  };
120 
121  matGradPtr = boost::make_shared<MatrixDouble>();
122  matStrainPtr = boost::make_shared<MatrixDouble>();
123  matStressPtr = boost::make_shared<MatrixDouble>();
124  matAccelerationPtr = boost::make_shared<MatrixDouble>();
125  matInertiaPtr = boost::make_shared<MatrixDouble>();
126  matDPtr = boost::make_shared<MatrixDouble>();
127 
128  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
129  matDPtr->resize(size_symm * size_symm, 1);
130 
131  CHKERR set_matrial_stiffens();
132 
134 }
135 //! [Create common data]
136 
137 //! [Run problem]
140  CHKERR readMesh();
141  CHKERR setupProblem();
142  CHKERR createCommonData();
143  CHKERR boundaryCondition();
144  CHKERR assembleSystem();
145  CHKERR solveSystem();
146  CHKERR outputResults();
147  CHKERR checkResults();
149 }
150 //! [Run problem]
151 
152 //! [Read mesh]
155  auto simple = mField.getInterface<Simple>();
156  CHKERR simple->getOptions();
157  CHKERR simple->loadFile();
159 }
160 //! [Read mesh]
161 
162 //! [Set up problem]
165  Simple *simple = mField.getInterface<Simple>();
166  // Add field
167  CHKERR simple->addDomainField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
168  SPACE_DIM);
169  int order = 3;
170  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
171  CHKERR simple->setFieldOrder("U", order);
172  CHKERR simple->setUp();
174 }
175 //! [Set up problem]
176 
177 //! [Boundary condition]
180 
181  auto fix_disp = [&](const std::string blockset_name) {
182  Range fix_ents;
184  if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
185  0) {
186  CHKERR mField.get_moab().get_entities_by_handle(it->meshset, fix_ents,
187  true);
188  }
189  }
190  return fix_ents;
191  };
192 
193  auto remove_ents = [&](const Range &&ents, const int lo, const int hi) {
194  auto prb_mng = mField.getInterface<ProblemsManager>();
195  auto simple = mField.getInterface<Simple>();
197  Range verts;
198  CHKERR mField.get_moab().get_connectivity(ents, verts, true);
199  verts.merge(ents);
200  if (SPACE_DIM == 3) {
201  Range adj;
202  CHKERR mField.get_moab().get_adjacencies(ents, 1, false, adj,
203  moab::Interface::UNION);
204  verts.merge(adj);
205  };
206  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
207  CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(), "U", verts,
208  lo, hi);
210  };
211 
212  CHKERR remove_ents(fix_disp("FIX_X"), 0, 0);
213  CHKERR remove_ents(fix_disp("FIX_Y"), 1, 1);
214  CHKERR remove_ents(fix_disp("FIX_Z"), 2, 2);
215  CHKERR remove_ents(fix_disp("FIX_ALL"), 0, 3);
216 
218 }
219 //! [Boundary condition]
220 
221 //! [Push operators to pipeline]
224  auto *simple = mField.getInterface<Simple>();
225  auto *pipeline_mng = mField.getInterface<PipelineManager>();
226 
227  if (SPACE_DIM == 2) {
228  pipeline_mng->getOpDomainLhsPipeline().push_back(
230  pipeline_mng->getOpDomainLhsPipeline().push_back(
232  pipeline_mng->getOpDomainRhsPipeline().push_back(
234  pipeline_mng->getOpDomainRhsPipeline().push_back(
236  }
237 
238  // Get pointer to U_tt shift in domain element
239  auto get_rho = [this](const double, const double, const double) {
240  auto *pipeline_mng = mField.getInterface<PipelineManager>();
241  auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
242  return rho * fe_domain_lhs->ts_aa;
243  };
244 
245  auto get_body_force = [this](const double, const double, const double) {
246  auto *pipeline_mng = mField.getInterface<PipelineManager>();
247  auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
248  auto t_source = FTensor::Tensor1<double, SPACE_DIM>{0.1, 1.};
250  const auto time = fe_domain_rhs->ts_t;
251  t_source(i) *= sin(time * omega * M_PI);
252  return t_source;
253  };
254 
255  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpK("U", "U", matDPtr));
256  pipeline_mng->getOpDomainLhsPipeline().push_back(
257  new OpMass("U", "U", get_rho));
258 
259  pipeline_mng->getOpDomainRhsPipeline().push_back(
260  new OpBodyForce("U", get_body_force));
261  pipeline_mng->getOpDomainRhsPipeline().push_back(
263  matGradPtr));
264  pipeline_mng->getOpDomainRhsPipeline().push_back(
265  new OpSymmetrizeTensor<SPACE_DIM>("U", matGradPtr, matStrainPtr));
266  pipeline_mng->getOpDomainRhsPipeline().push_back(
268  "U", matStrainPtr, matStressPtr, matDPtr));
269  pipeline_mng->getOpDomainRhsPipeline().push_back(
270  new OpInternalForce("U", matStressPtr));
271  pipeline_mng->getOpDomainRhsPipeline().push_back(
273  matAccelerationPtr));
274  pipeline_mng->getOpDomainRhsPipeline().push_back(
275  new OpScaleMatrix("U", rho, matAccelerationPtr, matInertiaPtr));
276  pipeline_mng->getOpDomainRhsPipeline().push_back(
277  new OpInertiaForce("U", matInertiaPtr));
278 
279  auto integration_rule = [](int, int, int approx_order) {
280  return 2 * approx_order;
281  };
282  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
283  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
284 
286 }
287 //! [Push operators to pipeline]
288 
289 /**
290  * @brief Monitor solution
291  *
292  * This functions is called by TS solver at the end of each step. It is used
293  * to output results to the hard drive.
294  */
295 struct Monitor : public FEMethod {
296  Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
297  : dM(dm), postProc(post_proc){};
300  constexpr int save_every_nth_step = 1;
301  if (ts_step % save_every_nth_step == 0) {
302  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
303  CHKERR postProc->writeFile(
304  "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
305  }
307  }
308 
309 private:
311  boost::shared_ptr<PostProcEle> postProc;
312 };
313 
314 //! [Solve]
317  auto *simple = mField.getInterface<Simple>();
318  auto *pipeline_mng = mField.getInterface<PipelineManager>();
319 
320  auto dm = simple->getDM();
321  auto ts = pipeline_mng->createTS2();
322 
323  // Setup postprocessing
324  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
325  post_proc_fe->generateReferenceElementMesh();
326  if (SPACE_DIM) {
327  post_proc_fe->getOpPtrVector().push_back(
329  post_proc_fe->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(invJac));
330  }
331  post_proc_fe->getOpPtrVector().push_back(
333  matGradPtr));
334  post_proc_fe->getOpPtrVector().push_back(
335  new OpSymmetrizeTensor<SPACE_DIM>("U", matGradPtr, matStrainPtr));
336  post_proc_fe->getOpPtrVector().push_back(
338  "U", matStrainPtr, matStressPtr, matDPtr));
339  post_proc_fe->getOpPtrVector().push_back(new OpPostProcElastic<SPACE_DIM>(
340  "U", post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts, matStrainPtr,
341  matStressPtr));
342  post_proc_fe->addFieldValuesPostProc("U");
343 
344  // Add monitor to time solver
345  boost::shared_ptr<FEMethod> null_fe;
346  auto monitor_ptr = boost::make_shared<Monitor>(dm, post_proc_fe);
347  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
348  null_fe, monitor_ptr);
349 
350  double ftime = 1;
351  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
352  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
353 
354  auto T = smartCreateDMVector(simple->getDM());
355  auto TT = smartVectorDuplicate(T);
356  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
357  SCATTER_FORWARD);
358 
359  CHKERR TS2SetSolution(ts, T, TT);
360  CHKERR TSSetFromOptions(ts);
361 
362  CHKERR TSSolve(ts, NULL);
363  CHKERR TSGetTime(ts, &ftime);
364 
365  PetscInt steps, snesfails, rejects, nonlinits, linits;
366  CHKERR TSGetTimeStepNumber(ts, &steps);
367  CHKERR TSGetSNESFailures(ts, &snesfails);
368  CHKERR TSGetStepRejections(ts, &rejects);
369  CHKERR TSGetSNESIterations(ts, &nonlinits);
370  CHKERR TSGetKSPIterations(ts, &linits);
371  MOFEM_LOG_C("EXAMPLE", Sev::inform,
372  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
373  "%D, linits %D\n",
374  steps, rejects, snesfails, ftime, nonlinits, linits);
375 
377 }
378 //! [Solve]
379 
380 //! [Postprocess results]
383  PetscBool test_flg = PETSC_FALSE;
384  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
385  if (test_flg) {
386  auto *simple = mField.getInterface<Simple>();
387  auto T = smartCreateDMVector(simple->getDM());
388  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
389  SCATTER_FORWARD);
390  double nrm2;
391  CHKERR VecNorm(T, NORM_2, &nrm2);
392  MOFEM_LOG("EXAMPLE", Sev::inform) << "Regression norm " << nrm2;
393  constexpr double regression_value = 1.09572;
394  if (fabs(nrm2 - regression_value) > 1e-2)
395  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
396  "Regression test faileed; wrong norm value.");
397  }
399 }
400 //! [Postprocess results]
401 
402 //! [Check]
406 }
407 //! [Check]
408 
409 static char help[] = "...\n\n";
410 
411 int main(int argc, char *argv[]) {
412 
413  // Initialisation of MoFEM/PETSc and MOAB data structures
414  const char param_file[] = "param_file.petsc";
415  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
416 
417  // Add logging channel for example
418  auto core_log = logging::core::get();
419  core_log->add_sink(
421  LogManager::setLog("EXAMPLE");
422  MOFEM_LOG_TAG("EXAMPLE", "example");
423 
424  try {
425 
426  //! [Register MoFEM discrete manager in PETSc]
427  DMType dm_name = "DMMOFEM";
428  CHKERR DMRegister_MoFEM(dm_name);
429  //! [Register MoFEM discrete manager in PETSc
430 
431  //! [Create MoAB]
432  moab::Core mb_instance; ///< mesh database
433  moab::Interface &moab = mb_instance; ///< mesh database interface
434  //! [Create MoAB]
435 
436  //! [Create MoFEM]
437  MoFEM::Core core(moab); ///< finite element database
438  MoFEM::Interface &m_field = core; ///< finite element database insterface
439  //! [Create MoFEM]
440 
441  //! [Example]
442  Example ex(m_field);
443  CHKERR ex.runProblem();
444  //! [Example]
445  }
446  CATCH_ERRORS;
447 
449 }
std::string param_file
MatrixDouble invJac
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:315
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:156
@ H1
continuous field
Definition: definitions.h:177
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ BLOCKSET
Definition: definitions.h:217
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:132
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
int main(int argc, char *argv[])
constexpr double rho
static char help[]
[Check]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
static double * ts_aa_ptr
constexpr int SPACE_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
constexpr double omega
DataForcesAndSourcesCore::EntData EntData
static double * ts_time_ptr
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
constexpr double young_modulus
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
Definition: elastic.cpp:39
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:939
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
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:445
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:368
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:312
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:343
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
const double T
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMass
VolumeElementForcesAndSourcesCore DomainEle
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:915
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)
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1129
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
Definition: plastic.cpp:59
static constexpr int approx_order
[Example]
Definition: plastic.cpp:117
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
boost::shared_ptr< MatrixDouble > matAccelerationPtr
boost::shared_ptr< MatrixDouble > matInertiaPtr
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
Managing BitRefLevels.
Core (interface) class.
Definition: Core.hpp:77
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Data on single entity (This is passed as argument to DataOperator::doWork)
Deprecated interface functions.
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:283
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:327
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field valuse for given petsc vector.
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return 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
Postprocess on face.
Post processing.
[Class definition]