v0.14.0
dynamic_first_order_con_law.cpp
Go to the documentation of this file.
1 /**
2  * \file dynamic_first_order_con_law.cpp
3  * \example dynamic_first_order_con_law.cpp
4  *
5  * Explicit first order conservation laws for solid dynamics
6  *
7  */
8 
9 #include <MoFEM.hpp>
10 #include <MatrixFunction.hpp>
11 
12 using namespace MoFEM;
13 
14 template <typename T> inline double trace(FTensor::Tensor2<T, 2, 2> &t_stress) {
15  constexpr double third = boost::math::constants::third<double>();
16  return (t_stress(0, 0) + t_stress(1, 1));
17 };
18 
19 template <typename T> inline double trace(FTensor::Tensor2<T, 3, 3> &t_stress) {
20  return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2));
21 };
22 
23 constexpr int SPACE_DIM =
24  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
25 
26 using EntData = EntitiesFieldData::EntData;
27 using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
28 using DomainEleOp = DomainEle::UserDataOperator;
29 using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
30 using PostProcFaceEle =
31  PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore>;
32 
33 using BoundaryEle =
34  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
35 using BoundaryEleOp = BoundaryEle::UserDataOperator;
36 using SetPtsData = FieldEvaluatorInterface::SetPtsData;
37 
38 template <int DIM> struct PostProcEleByDim;
39 
40 template <> struct PostProcEleByDim<2> {
41  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
42  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
43  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
44 };
45 
46 template <> struct PostProcEleByDim<3> {
47  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
48  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
49  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
50 };
51 
52 using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;
53 using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;
54 using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;
55 
56 using OpMassV = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
57  GAUSS>::OpMass<1, SPACE_DIM>;
58 using OpMassF = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
59  GAUSS>::OpMass<1, SPACE_DIM * SPACE_DIM>;
60 
61 using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
62  PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 1>;
63 
64 using OpBodyForce = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
65  GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 0>;
66 
67 using DomainNaturalBC =
68  NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
69 using OpBodyForceVector =
70  DomainNaturalBC::OpFlux<NaturalMeshsetTypeVectorScaling<BLOCKSET>, 1,
71  SPACE_DIM>;
72 
73 using OpGradTimesTensor2 =
74  FormsIntegrators<DomainEleOp>::Assembly<AssemblyType::PETSC>::LinearForm<
75  IntegrationType::GAUSS>::OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;
76 
77 using OpRhsTestPiola =
78  FormsIntegrators<DomainEleOp>::Assembly<AssemblyType::PETSC>::LinearForm<
79  IntegrationType::GAUSS>::OpBaseTimesVector<1, SPACE_DIM * SPACE_DIM, 1>;
80 
81 using OpGradTimesPiola =
82  FormsIntegrators<DomainEleOp>::Assembly<AssemblyType::PETSC>::LinearForm<
83  IntegrationType::GAUSS>::OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;
84 
85 using BoundaryNaturalBC =
86  NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
87 using OpForce = BoundaryNaturalBC::OpFlux<NaturalForceMeshsets, 1, SPACE_DIM>;
88 
89 /** \brief Save field DOFS on vertices/tags
90  */
91 
92 constexpr double omega = 1.;
93 constexpr double young_modulus = 1.;
94 constexpr double poisson_ratio = 0.;
95 double bulk_modulus_K = young_modulus / (3. * (1. - 2. * poisson_ratio));
96 double shear_modulus_G = young_modulus / (2. * (1. + poisson_ratio));
97 double mu = young_modulus / (2. * (1. + poisson_ratio));
98 double lamme_lambda = young_modulus * poisson_ratio /
99  ((1. + poisson_ratio) * (1. - 2. * poisson_ratio));
100 
101 // Operator to Calculate F
102 template <int DIM_0, int DIM_1>
103 struct OpCalculateFStab : public ForcesAndSourcesCore::UserDataOperator {
104  OpCalculateFStab(boost::shared_ptr<MatrixDouble> def_grad_ptr,
105  boost::shared_ptr<MatrixDouble> def_grad_stab_ptr,
106  boost::shared_ptr<MatrixDouble> def_grad_dot_ptr,
107  double tau_F_ptr, double xi_F_ptr,
108  boost::shared_ptr<MatrixDouble> grad_x_ptr,
109  boost::shared_ptr<MatrixDouble> grad_vel_ptr)
110  : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),
111  defGradPtr(def_grad_ptr), defGradStabPtr(def_grad_stab_ptr),
112  defGradDotPtr(def_grad_dot_ptr), tauFPtr(tau_F_ptr), xiF(xi_F_ptr),
113  gradxPtr(grad_x_ptr), gradVelPtr(grad_vel_ptr) {}
114 
115  MoFEMErrorCode doWork(int side, EntityType type,
116  DataForcesAndSourcesCore::EntData &data) {
117  MoFEMFunctionBegin;
118  // Define Indicies
119  FTensor::Index<'i', SPACE_DIM> i;
120  FTensor::Index<'j', SPACE_DIM> j;
121 
122  // Number of Gauss points
123  const size_t nb_gauss_pts = getGaussPts().size2();
124 
125  defGradStabPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false);
126  defGradStabPtr->clear();
127 
128  // Extract matrix from data matrix
129  auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
130  auto t_Fstab = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradStabPtr);
131  auto t_F_dot = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradDotPtr);
132 
133  // tau_F = alpha deltaT
134  auto tau_F = tauFPtr;
135  double xi_F = xiF;
136  auto t_gradx = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradxPtr);
137  auto t_gradVel = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradVelPtr);
138 
139  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
140  // Stabilised Deformation Gradient
141  t_Fstab(i, j) = t_F(i, j) + tau_F * (t_gradVel(i, j) - t_F_dot(i, j)) +
142  xi_F * (t_gradx(i, j) - t_F(i, j));
143 
144  ++t_F;
145  ++t_Fstab;
146  ++t_gradVel;
147  ++t_F_dot;
148 
149  ++t_gradx;
150  }
151 
152  MoFEMFunctionReturn(0);
153  }
154 
155 private:
156  double tauFPtr;
157  double xiF;
158  boost::shared_ptr<MatrixDouble> defGradPtr;
159  boost::shared_ptr<MatrixDouble> defGradStabPtr;
160  boost::shared_ptr<MatrixDouble> defGradDotPtr;
161  boost::shared_ptr<MatrixDouble> gradxPtr;
162  boost::shared_ptr<MatrixDouble> gradVelPtr;
163 };
164 
165 // Operator to Calculate P
166 template <int DIM_0, int DIM_1>
167 struct OpCalculatePiola : public ForcesAndSourcesCore::UserDataOperator {
168  OpCalculatePiola(double shear_modulus, double bulk_modulus, double m_u,
169  double lambda_lamme,
170  boost::shared_ptr<MatrixDouble> first_piola_ptr,
171  boost::shared_ptr<MatrixDouble> def_grad_ptr)
172  : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),
173  shearModulus(shear_modulus), bulkModulus(bulk_modulus), mU(m_u),
174  lammeLambda(lambda_lamme), firstPiolaPtr(first_piola_ptr),
175  defGradPtr(def_grad_ptr) {}
176 
177  MoFEMErrorCode doWork(int side, EntityType type,
178  DataForcesAndSourcesCore::EntData &data) {
179  MoFEMFunctionBegin;
180  // Define Indicies
181  FTensor::Index<'i', SPACE_DIM> i;
182  FTensor::Index<'j', SPACE_DIM> j;
183  FTensor::Index<'k', SPACE_DIM> k;
184 
185  // Define Kronecker Delta
186  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
187 
188  // Number of Gauss points
189  const size_t nb_gauss_pts = getGaussPts().size2();
190 
191  // Resize Piola
192  firstPiolaPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false); // ignatios check
193  firstPiolaPtr->clear();
194 
195  // Extract matrix from data matrix
196  auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*firstPiolaPtr);
197  auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
198  const double two_o_three = 2. / 3.;
199  const double trace_t_dk = DIM_0;
200  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
201 
202  t_P(i, j) = shearModulus * (t_F(i, j) + t_F(j, i) - 2. * t_kd(i, j) -
203  two_o_three * trace(t_F) * t_kd(i, j) +
204  two_o_three * trace_t_dk * t_kd(i, j)) +
205  bulkModulus * trace(t_F) * t_kd(i, j) -
206  bulkModulus * trace_t_dk * t_kd(i, j);
207 
208  ++t_F;
209  ++t_P;
210  }
211 
212  MoFEMFunctionReturn(0);
213  }
214 
215 private:
216  double shearModulus;
217  double bulkModulus;
218  double mU;
219  double lammeLambda;
220  boost::shared_ptr<MatrixDouble> firstPiolaPtr;
221  boost::shared_ptr<MatrixDouble> defGradPtr;
222 };
223 
224 template <int DIM>
225 struct OpCalculateDisplacement : public ForcesAndSourcesCore::UserDataOperator {
226  OpCalculateDisplacement(boost::shared_ptr<MatrixDouble> spatial_pos_ptr,
227  boost::shared_ptr<MatrixDouble> reference_pos_ptr,
228  boost::shared_ptr<MatrixDouble> u_ptr)
229  : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),
230  xPtr(spatial_pos_ptr), XPtr(reference_pos_ptr), uPtr(u_ptr) {}
231 
232  MoFEMErrorCode doWork(int side, EntityType type,
233  DataForcesAndSourcesCore::EntData &data) {
234  MoFEMFunctionBegin;
235  // Define Indicies
236  FTensor::Index<'i', DIM> i;
237 
238  // Number of Gauss points
239  const size_t nb_gauss_pts = getGaussPts().size2();
240 
241  uPtr->resize(DIM, nb_gauss_pts, false); // ignatios check
242  uPtr->clear();
243 
244  // Extract matrix from data matrix
245  auto t_x = getFTensor1FromMat<DIM>(*xPtr);
246  auto t_X = getFTensor1FromMat<DIM>(*XPtr);
247  auto t_u = getFTensor1FromMat<DIM>(*uPtr);
248  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
249 
250  t_u(i) = t_x(i) - t_X(i);
251  ++t_u;
252  ++t_x;
253  ++t_X;
254  }
255 
256  MoFEMFunctionReturn(0);
257  }
258 
259 private:
260  boost::shared_ptr<MatrixDouble> xPtr;
261  boost::shared_ptr<MatrixDouble> XPtr;
262  boost::shared_ptr<MatrixDouble> uPtr;
263 };
264 
265 template <int DIM_0, int DIM_1>
266 struct OpCalculatePiolaIncompressibleNH
267  : public ForcesAndSourcesCore::UserDataOperator {
268  OpCalculatePiolaIncompressibleNH(
269  double shear_modulus, double bulk_modulus, double m_u,
270  double lambda_lamme, boost::shared_ptr<MatrixDouble> first_piola_ptr,
271  boost::shared_ptr<MatrixDouble> def_grad_ptr,
272  boost::shared_ptr<MatrixDouble> inv_def_grad_ptr,
273  boost::shared_ptr<VectorDouble> det)
274  : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),
275  shearModulus(shear_modulus), bulkModulus(bulk_modulus), mU(m_u),
276  lammeLambda(lambda_lamme), firstPiolaPtr(first_piola_ptr),
277  defGradPtr(def_grad_ptr), invDefGradPtr(inv_def_grad_ptr), dEt(det) {}
278 
279  MoFEMErrorCode doWork(int side, EntityType type,
280  DataForcesAndSourcesCore::EntData &data) {
281  MoFEMFunctionBegin;
282  // Define Indicies
283  FTensor::Index<'i', SPACE_DIM> i;
284  FTensor::Index<'j', SPACE_DIM> j;
285  FTensor::Index<'k', SPACE_DIM> k;
286  FTensor::Index<'l', SPACE_DIM> l;
287 
288  // Define Kronecker Delta
289  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
290 
291  // Number of Gauss points
292  const size_t nb_gauss_pts = getGaussPts().size2();
293 
294  // Resize Piola
295  firstPiolaPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false); // ignatios check
296  firstPiolaPtr->clear();
297 
298  // Extract matrix from data matrix
299  auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*firstPiolaPtr);
300  auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
301  auto t_inv_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*invDefGradPtr);
302  auto t_det = getFTensor0FromVec<1>(*dEt);
303  const double two_o_three = 2. / 3.;
304  const double one_o_three = 1. / 3.;
305  const double bulk_mod = bulkModulus;
306  const double shear_mod = shearModulus;
307  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
308 
309  // Nearly incompressible NH
310  // volumetric part
311  t_P(i, j) = bulk_mod * (t_det - 1.) * t_det * t_inv_F(j, i);
312  // deviatoric part
313  t_P(i, j) +=
314  shear_mod * pow(t_det, two_o_three) *
315  (t_F(i, j) - one_o_three * (t_F(l, k) * t_F(l, k)) * t_inv_F(j, i));
316 
317  ++t_F;
318  ++t_P;
319  ++t_inv_F;
320  ++t_det;
321  }
322 
323  MoFEMFunctionReturn(0);
324  }
325 
326 private:
327  double shearModulus;
328  double bulkModulus;
329  double mU;
330  double lammeLambda;
331  boost::shared_ptr<MatrixDouble> firstPiolaPtr;
332  boost::shared_ptr<MatrixDouble> defGradPtr;
333  boost::shared_ptr<MatrixDouble> invDefGradPtr;
334  boost::shared_ptr<VectorDouble> dEt;
335 };
336 
337 template <int DIM_0, int DIM_1>
338 struct OpCalculateDeformationGradient
339  : public ForcesAndSourcesCore::UserDataOperator {
340  OpCalculateDeformationGradient(
341  boost::shared_ptr<MatrixDouble> def_grad_ptr,
342  boost::shared_ptr<MatrixDouble> grad_tensor_ptr)
343  : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),
344  defGradPtr(def_grad_ptr), gradTensorPtr(grad_tensor_ptr) {}
345 
346  MoFEMErrorCode doWork(int side, EntityType type,
347  DataForcesAndSourcesCore::EntData &data) {
348  MoFEMFunctionBegin;
349  // Define Indicies
350  FTensor::Index<'i', SPACE_DIM> i;
351  FTensor::Index<'j', SPACE_DIM> j;
352 
353  // Define Kronecker Delta
354  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
355 
356  // Number of Gauss points
357  const size_t nb_gauss_pts = getGaussPts().size2();
358 
359  // Resize Piola
360  defGradPtr->resize(DIM_0 * DIM_1, nb_gauss_pts, false);
361  defGradPtr->clear();
362 
363  // Extract matrix from data matrix
364  auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*defGradPtr);
365  auto t_H = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*gradTensorPtr);
366  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
367 
368  t_F(i, j) = t_H(i, j) + t_kd(i, j);
369 
370  ++t_F;
371  ++t_H;
372  }
373 
374  MoFEMFunctionReturn(0);
375  }
376 
377 private:
378  boost::shared_ptr<MatrixDouble> gradTensorPtr;
379  boost::shared_ptr<MatrixDouble> defGradPtr;
380 };
381 
382 struct Example;
383 struct TSPrePostProc {
384 
385  TSPrePostProc() = default;
386  virtual ~TSPrePostProc() = default;
387 
388  /**
389  * @brief Used to setup TS solver
390  *
391  * @param ts
392  * @return MoFEMErrorCode
393  */
394  MoFEMErrorCode tsSetUp(TS ts);
395 
396  // SmartPetscObj<VecScatter> getScatter(Vec x, Vec y, enum FR fr);
397  Example *fsRawPtr;
398  static MoFEMErrorCode tsPostStage(TS ts, PetscReal stagetime,
399  PetscInt stageindex, Vec *Y);
400  static MoFEMErrorCode tsPostStep(TS ts);
401  static MoFEMErrorCode tsPreStep(TS ts);
402 };
403 
404 static boost::weak_ptr<TSPrePostProc> tsPrePostProc;
405 
406 struct LinMomTimeScale : public MoFEM::TimeScale {
407  using MoFEM::TimeScale::TimeScale;
408  double getScale(const double time) {
409  return sin(2. * M_PI * MoFEM::TimeScale::getScale(time));
410  };
411 };
412 
413 struct CommonData {
414  SmartPetscObj<Mat> M; ///< Mass matrix
415  SmartPetscObj<KSP> ksp; ///< Linear solver
416 };
417 
418 struct Example {
419 
420  Example(MoFEM::Interface &m_field) : mField(m_field) {}
421 
422  MoFEMErrorCode runProblem();
423 
424 private:
425  MoFEM::Interface &mField;
426 
427  MoFEMErrorCode readMesh();
428  MoFEMErrorCode setupProblem();
429  MoFEMErrorCode boundaryCondition();
430  MoFEMErrorCode assembleSystem();
431  MoFEMErrorCode solveSystem();
432  MoFEMErrorCode outputResults();
433  MoFEMErrorCode checkResults();
434  friend struct TSPrePostProc;
435 
436  struct DynamicFirstOrderConsSinusTimeScale : public MoFEM::TimeScale {
437  using MoFEM::TimeScale::TimeScale;
438  double getScale(const double time) { return 0.001 * sin(0.1 * time); };
439  };
440 
441  struct DynamicFirstOrderConsConstantTimeScale : public MoFEM::TimeScale {
442  using MoFEM::TimeScale::TimeScale;
443  double getScale(const double time) { return 0.001; };
444  };
445 };
446 
447 //! [Run problem]
448 MoFEMErrorCode Example::runProblem() {
449  MoFEMFunctionBegin;
450  CHKERR readMesh();
451  CHKERR setupProblem();
452  CHKERR boundaryCondition();
453  CHKERR assembleSystem();
454  CHKERR solveSystem();
455  CHKERR outputResults();
456  CHKERR checkResults();
457  MoFEMFunctionReturn(0);
458 }
459 //! [Run problem]
460 
461 //! [Read mesh]
462 MoFEMErrorCode Example::readMesh() {
463  MoFEMFunctionBegin;
464  auto simple = mField.getInterface<Simple>();
465 
466  CHKERR simple->getOptions();
467  CHKERR simple->loadFile();
468  MoFEMFunctionReturn(0);
469 }
470 //! [Read mesh]
471 
472 //! [Set up problem]
473 MoFEMErrorCode Example::setupProblem() {
474  MoFEMFunctionBegin;
475  Simple *simple = mField.getInterface<Simple>();
476  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
477  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
478  PetscInt choice_base_value = AINSWORTH;
479  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
480  LASBASETOPT, &choice_base_value, PETSC_NULL);
481 
482  FieldApproximationBase base;
483  switch (choice_base_value) {
484  case AINSWORTH:
485  base = AINSWORTH_LEGENDRE_BASE;
486  MOFEM_LOG("WORLD", Sev::inform)
487  << "Set AINSWORTH_LEGENDRE_BASE for displacements";
488  break;
489  case DEMKOWICZ:
490  base = DEMKOWICZ_JACOBI_BASE;
491  MOFEM_LOG("WORLD", Sev::inform)
492  << "Set DEMKOWICZ_JACOBI_BASE for displacements";
493  break;
494  default:
495  base = LASTBASE;
496  break;
497  }
498  // Add field
499  CHKERR simple->addDomainField("V", H1, base, SPACE_DIM);
500  CHKERR simple->addBoundaryField("V", H1, base, SPACE_DIM);
501  CHKERR simple->addDomainField("F", H1, base, SPACE_DIM * SPACE_DIM);
502  CHKERR simple->addDataField("x_1", H1, base, SPACE_DIM);
503  CHKERR simple->addDataField("x_2", H1, base, SPACE_DIM);
504  CHKERR simple->addDataField("F_0", H1, base, SPACE_DIM * SPACE_DIM);
505  CHKERR simple->addDataField("F_dot", H1, base, SPACE_DIM * SPACE_DIM);
506 
507  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
508  int order = 2;
509  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
510  CHKERR simple->setFieldOrder("V", order);
511  CHKERR simple->setFieldOrder("F", order);
512  CHKERR simple->setFieldOrder("F_0", order);
513  CHKERR simple->setFieldOrder("F_dot", order);
514  CHKERR simple->setFieldOrder("x_1", order);
515  CHKERR simple->setFieldOrder("x_2", order);
516  CHKERR simple->setFieldOrder("GEOMETRY", order);
517  CHKERR simple->setUp();
518 
519  auto project_ho_geometry = [&]() {
520  Projection10NodeCoordsOnField ent_method_x(mField, "x_1");
521  CHKERR mField.loop_dofs("x_1", ent_method_x);
522  Projection10NodeCoordsOnField ent_method_x_2(mField, "x_2");
523  CHKERR mField.loop_dofs("x_2", ent_method_x_2);
524 
525  Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
526  return mField.loop_dofs("GEOMETRY", ent_method);
527  };
528  CHKERR project_ho_geometry();
529 
530  MoFEMFunctionReturn(0);
531 }
532 //! [Set up problem]
533 
534 //! [Boundary condition]
535 MoFEMErrorCode Example::boundaryCondition() {
536  MoFEMFunctionBegin;
537 
538  auto simple = mField.getInterface<Simple>();
539  auto bc_mng = mField.getInterface<BcManager>();
540  auto *pipeline_mng = mField.getInterface<PipelineManager>();
541  auto time_scale = boost::make_shared<TimeScale>();
542 
543  PetscBool sin_time_function = PETSC_FALSE;
544  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-sin_time_function",
545  &sin_time_function, PETSC_NULL);
546 
547  if (sin_time_function)
548  time_scale = boost::make_shared<DynamicFirstOrderConsSinusTimeScale>();
549  else
550  time_scale = boost::make_shared<DynamicFirstOrderConsConstantTimeScale>();
551 
552  pipeline_mng->getBoundaryExplicitRhsFE().reset();
553  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
554  pipeline_mng->getOpBoundaryExplicitRhsPipeline(), {NOSPACE}, "GEOMETRY");
555 
556  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
557  pipeline_mng->getOpBoundaryExplicitRhsPipeline(), mField, "V",
558  {time_scale}, "FORCE", Sev::inform);
559 
560  auto integration_rule = [](int, int, int approx_order) {
561  return 2 * approx_order;
562  };
563 
564  CHKERR pipeline_mng->setBoundaryExplicitRhsIntegrationRule(integration_rule);
565  CHKERR pipeline_mng->setDomainExplicitRhsIntegrationRule(integration_rule);
566 
567  CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
568  simple->getProblemName(), "V");
569 
570  auto get_pre_proc_hook = [&]() {
571  return EssentialPreProc<DisplacementCubitBcData>(
572  mField, pipeline_mng->getDomainExplicitRhsFE(), {time_scale});
573  };
574  pipeline_mng->getDomainExplicitRhsFE()->preProcessHook = get_pre_proc_hook();
575 
576  MoFEMFunctionReturn(0);
577 }
578 //! [Boundary condition]
579 
580 MoFEMErrorCode TSPrePostProc::tsPostStage(TS ts, PetscReal stagetime,
581  PetscInt stageindex, Vec *Y) {
582  MoFEMFunctionBegin;
583  // cerr << "tsPostStage " <<"\n";
584  if (auto ptr = tsPrePostProc.lock()) {
585  auto &m_field = ptr->fsRawPtr->mField;
586 
587  auto fb = m_field.getInterface<FieldBlas>();
588  double dt;
589  CHKERR TSGetTimeStep(ts, &dt);
590  double time;
591  CHKERR TSGetTime(ts, &time);
592  PetscInt num_stages;
593  Vec *stage_solutions;
594 
595  CHKERR TSGetStages(ts, &num_stages, &stage_solutions);
596  PetscPrintf(PETSC_COMM_WORLD, "Check timestep %d time %e dt %e\n",
597  num_stages, time, dt);
598 
599  const double inv_num_step = (double)num_stages;
600  CHKERR fb->fieldCopy(1., "x_1", "x_2");
601  CHKERR fb->fieldAxpy(dt, "V", "x_2");
602  CHKERR fb->fieldCopy(1., "x_2", "x_1");
603 
604  CHKERR fb->fieldCopy(-inv_num_step / dt, "F_0", "F_dot");
605  CHKERR fb->fieldAxpy(inv_num_step / dt, "F", "F_dot");
606  CHKERR fb->fieldCopy(1., "F", "F_0");
607  }
608  MoFEMFunctionReturn(0);
609 }
610 
611 MoFEMErrorCode TSPrePostProc::tsPostStep(TS ts) {
612  MoFEMFunctionBegin;
613 
614  if (auto ptr = tsPrePostProc.lock()) {
615  auto &m_field = ptr->fsRawPtr->mField;
616  double dt;
617  CHKERR TSGetTimeStep(ts, &dt);
618  double time;
619  CHKERR TSGetTime(ts, &time);
620  }
621  MoFEMFunctionReturn(0);
622 }
623 
624 MoFEMErrorCode TSPrePostProc::tsPreStep(TS ts) {
625  MoFEMFunctionBegin;
626 
627  if (auto ptr = tsPrePostProc.lock()) {
628  auto &m_field = ptr->fsRawPtr->mField;
629  double dt;
630  CHKERR TSGetTimeStep(ts, &dt);
631  double time;
632  CHKERR TSGetTime(ts, &time);
633  int step_num;
634  CHKERR TSGetStepNumber(ts, &step_num);
635  }
636  MoFEMFunctionReturn(0);
637 }
638 
639 //! [Push operators to pipeline]
640 MoFEMErrorCode Example::assembleSystem() {
641  MoFEMFunctionBegin;
642  auto get_body_force = [this](const double, const double, const double) {
643  FTensor::Index<'i', SPACE_DIM> i;
644  FTensor::Tensor1<double, SPACE_DIM> t_source;
645  t_source(i) = 0.;
646  t_source(0) = 0.1;
647  t_source(1) = 1.;
648  return t_source;
649  };
650 
651  // specific time scaling
652  auto get_time_scale = [this](const double time) {
653  return sin(time * omega * M_PI);
654  };
655 
656  auto apply_rhs = [&](auto &pip) {
657  MoFEMFunctionBegin;
658 
659  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
660  "GEOMETRY");
661 
662  // Calculate Gradient of velocity
663  auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
664  pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
665  "V", mat_v_grad_ptr));
666 
667  auto gravity_vector_ptr = boost::make_shared<MatrixDouble>();
668  gravity_vector_ptr->resize(SPACE_DIM, 1);
669  auto set_body_force = [&]() {
670  FTensor::Index<'i', SPACE_DIM> i;
671  MoFEMFunctionBegin;
672  auto t_force = getFTensor1FromMat<SPACE_DIM, 0>(*gravity_vector_ptr);
673  double unit_weight = 0.;
674  CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-unit_weight", &unit_weight,
675  PETSC_NULL);
676  t_force(i) = 0;
677  if (SPACE_DIM == 2) {
678  t_force(1) = -unit_weight;
679  } else if (SPACE_DIM == 3) {
680  t_force(2) = unit_weight;
681  }
682  MoFEMFunctionReturn(0);
683  };
684 
685  CHKERR set_body_force();
686  pip.push_back(new OpBodyForce("V", gravity_vector_ptr,
687  [](double, double, double) { return 1.; }));
688 
689  // Calculate unknown F
690  auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
691  pip.push_back(new OpCalculateTensor2FieldValues<SPACE_DIM, SPACE_DIM>(
692  "F", mat_H_tensor_ptr));
693 
694  // // Calculate F
695  double tau = 0.2;
696  CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-tau", &tau, PETSC_NULL);
697 
698  double xi = 0.;
699  CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-xi", &xi, PETSC_NULL);
700 
701  // Calculate P stab
702  auto one = [&](const double, const double, const double) {
703  return 3. * bulk_modulus_K;
704  };
705  auto minus_one = [](const double, const double, const double) {
706  return -1.;
707  };
708 
709  auto mat_dot_F_tensor_ptr = boost::make_shared<MatrixDouble>();
710  pip.push_back(new OpCalculateTensor2FieldValues<SPACE_DIM, SPACE_DIM>(
711  "F_dot", mat_dot_F_tensor_ptr));
712 
713  // Calculate Gradient of Spatial Positions
714  auto mat_x_grad_ptr = boost::make_shared<MatrixDouble>();
715  pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
716  "x_2", mat_x_grad_ptr));
717 
718  auto mat_F_tensor_ptr = boost::make_shared<MatrixDouble>();
719  pip.push_back(new OpCalculateDeformationGradient<SPACE_DIM, SPACE_DIM>(
720  mat_F_tensor_ptr, mat_H_tensor_ptr));
721 
722  auto mat_F_stab_ptr = boost::make_shared<MatrixDouble>();
723  pip.push_back(new OpCalculateFStab<SPACE_DIM, SPACE_DIM>(
724  mat_F_tensor_ptr, mat_F_stab_ptr, mat_dot_F_tensor_ptr, tau, xi,
725  mat_x_grad_ptr, mat_v_grad_ptr));
726 
727  PetscBool is_linear_elasticity = PETSC_TRUE;
728  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_linear_elasticity",
729  &is_linear_elasticity, PETSC_NULL);
730 
731  auto mat_P_stab_ptr = boost::make_shared<MatrixDouble>();
732  if (is_linear_elasticity) {
733  pip.push_back(new OpCalculatePiola<SPACE_DIM, SPACE_DIM>(
734  shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_stab_ptr,
735  mat_F_stab_ptr));
736  } else {
737  auto inv_F = boost::make_shared<MatrixDouble>();
738  auto det_ptr = boost::make_shared<VectorDouble>();
739 
740  pip.push_back(
741  new OpInvertMatrix<SPACE_DIM>(mat_F_stab_ptr, det_ptr, inv_F));
742 
743  // OpCalculatePiolaIncompressibleNH
744  pip.push_back(new OpCalculatePiolaIncompressibleNH<SPACE_DIM, SPACE_DIM>(
745  shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_stab_ptr,
746  mat_F_stab_ptr, inv_F, det_ptr));
747  }
748 
749  pip.push_back(new OpGradTimesTensor2("V", mat_P_stab_ptr, minus_one));
750  pip.push_back(new OpRhsTestPiola("F", mat_v_grad_ptr, one));
751 
752  MoFEMFunctionReturn(0);
753  };
754 
755  auto *pipeline_mng = mField.getInterface<PipelineManager>();
756  CHKERR apply_rhs(pipeline_mng->getOpDomainExplicitRhsPipeline());
757 
758  auto integration_rule = [](int, int, int approx_order) {
759  return 2 * approx_order;
760  };
761  CHKERR pipeline_mng->setDomainExplicitRhsIntegrationRule(integration_rule);
762 
763  MoFEMFunctionReturn(0);
764 }
765 //! [Push operators to pipeline]
766 
767 /**
768  * @brief Monitor solution
769  *
770  * This functions is called by TS solver at the end of each step. It is used
771  * to output results to the hard drive.
772  */
773 
774 struct Monitor : public FEMethod {
775  MoFEM::Interface &mField;
776  Monitor(SmartPetscObj<DM> dm, MoFEM::Interface &m_field,
777  boost::shared_ptr<PostProcEle> post_proc,
778  boost::shared_ptr<PostProcFaceEle> post_proc_bdry,
779  boost::shared_ptr<MatrixDouble> velocity_field_ptr,
780  boost::shared_ptr<MatrixDouble> x2_field_ptr,
781  boost::shared_ptr<MatrixDouble> geometry_field_ptr,
782  std::array<double, 3> pass_field_eval_coords,
783  boost::shared_ptr<SetPtsData> pass_field_eval_data)
784  : dM(dm), mField(m_field), postProc(post_proc),
785  postProcBdy(post_proc_bdry), velocityFieldPtr(velocity_field_ptr),
786  x2FieldPtr(x2_field_ptr), geometryFieldPtr(geometry_field_ptr),
787  fieldEvalCoords(pass_field_eval_coords),
788  fieldEvalData(pass_field_eval_data){};
789  MoFEMErrorCode postProcess() {
790  MoFEMFunctionBegin;
791 
792  // cerr << "wagawaga\n";
793  auto *simple = mField.getInterface<Simple>();
794 
795  if (SPACE_DIM == 3) {
796  CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint3D(
797  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
798  simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
799  mField.get_comm_rank(), getCacheWeakPtr().lock(), MF_EXIST, QUIET);
800  } else {
801  CHKERR mField.getInterface<FieldEvaluatorInterface>()->evalFEAtThePoint2D(
802  fieldEvalCoords.data(), 1e-12, simple->getProblemName(),
803  simple->getDomainFEName(), fieldEvalData, mField.get_comm_rank(),
804  mField.get_comm_rank(), getCacheWeakPtr().lock(), MF_EXIST, QUIET);
805  }
806 
807  if (velocityFieldPtr->size1()) {
808  auto t_vel = getFTensor1FromMat<SPACE_DIM>(*velocityFieldPtr);
809  auto t_x2_field = getFTensor1FromMat<SPACE_DIM>(*x2FieldPtr);
810  auto t_geom = getFTensor1FromMat<SPACE_DIM>(*geometryFieldPtr);
811 
812  double u_x = t_x2_field(0) - t_geom(0);
813  double u_y = t_x2_field(1) - t_geom(1);
814  double u_z = t_x2_field(2) - t_geom(2);
815 
816  MOFEM_LOG("SYNC", Sev::inform)
817  << "Velocities x: " << t_vel(0) << " y: " << t_vel(1)
818  << " z: " << t_vel(2) << "\n";
819  MOFEM_LOG("SYNC", Sev::inform) << "Displacement x: " << u_x
820  << " y: " << u_y << " z: " << u_z << "\n";
821  }
822 
823  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
824  std::regex((boost::format("%s(.*)") % "Data_Vertex").str()))) {
825  Range ents;
826  mField.get_moab().get_entities_by_dimension(m->getMeshset(), 0, ents,
827  true);
828  auto print_vets = [](boost::shared_ptr<FieldEntity> ent_ptr) {
829  MoFEMFunctionBegin;
830  if (!(ent_ptr->getPStatus() & PSTATUS_NOT_OWNED)) {
831  MOFEM_LOG("SYNC", Sev::inform)
832  << "Velocities: " << ent_ptr->getEntFieldData()[0] << " "
833  << ent_ptr->getEntFieldData()[1] << " "
834  << ent_ptr->getEntFieldData()[2] << "\n";
835  }
836  MoFEMFunctionReturn(0);
837  };
838  CHKERR mField.getInterface<FieldBlas>()->fieldLambdaOnEntities(
839  print_vets, "V", &ents);
840  }
841  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::inform);
842 
843  PetscBool print_volume = PETSC_FALSE;
844  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-print_volume", &print_volume,
845  PETSC_NULL);
846 
847  PetscBool print_skin = PETSC_FALSE;
848  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-print_skin", &print_skin,
849  PETSC_NULL);
850 
851  int save_every_nth_step = 1;
852  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_step",
853  &save_every_nth_step, PETSC_NULL);
854  if (ts_step % save_every_nth_step == 0) {
855 
856  if (print_volume) {
857  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
858  CHKERR postProc->writeFile(
859  "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
860  }
861 
862  if (print_skin) {
863  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcBdy);
864  CHKERR postProcBdy->writeFile(
865  "out_boundary_" + boost::lexical_cast<std::string>(ts_step) +
866  ".h5m");
867  }
868  }
869  MoFEMFunctionReturn(0);
870  }
871 
872 private:
873  SmartPetscObj<DM> dM;
874  boost::shared_ptr<PostProcEle> postProc;
875  boost::shared_ptr<PostProcFaceEle> postProcBdy;
876  boost::shared_ptr<MatrixDouble> velocityFieldPtr;
877  boost::shared_ptr<MatrixDouble> x2FieldPtr;
878  boost::shared_ptr<MatrixDouble> geometryFieldPtr;
879  std::array<double, 3> fieldEvalCoords;
880  boost::shared_ptr<SetPtsData> fieldEvalData;
881 };
882 
883 //! [Solve]
884 MoFEMErrorCode Example::solveSystem() {
885  MoFEMFunctionBegin;
886  auto *simple = mField.getInterface<Simple>();
887  auto *pipeline_mng = mField.getInterface<PipelineManager>();
888 
889  auto dm = simple->getDM();
890 
891  auto calculate_stress_ops = [&](auto &pip) {
892  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
893 
894  auto v_ptr = boost::make_shared<MatrixDouble>();
895  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("V", v_ptr));
896  auto X_ptr = boost::make_shared<MatrixDouble>();
897  pip.push_back(
898  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
899 
900  auto x_ptr = boost::make_shared<MatrixDouble>();
901  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("x_1", x_ptr));
902 
903  // Calculate unknown F
904  auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
905  pip.push_back(new OpCalculateTensor2FieldValues<SPACE_DIM, SPACE_DIM>(
906  "F", mat_H_tensor_ptr));
907 
908  auto u_ptr = boost::make_shared<MatrixDouble>();
909  pip.push_back(new OpCalculateDisplacement<SPACE_DIM>(x_ptr, X_ptr, u_ptr));
910  // Calculate P
911 
912  auto mat_F_ptr = boost::make_shared<MatrixDouble>();
913  pip.push_back(new OpCalculateDeformationGradient<SPACE_DIM, SPACE_DIM>(
914  mat_F_ptr, mat_H_tensor_ptr));
915 
916  PetscBool is_linear_elasticity = PETSC_TRUE;
917  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_linear_elasticity",
918  &is_linear_elasticity, PETSC_NULL);
919 
920  auto mat_P_ptr = boost::make_shared<MatrixDouble>();
921  if (is_linear_elasticity) {
922  pip.push_back(new OpCalculatePiola<SPACE_DIM, SPACE_DIM>(
923  shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_ptr,
924  mat_F_ptr));
925  } else {
926  auto inv_F = boost::make_shared<MatrixDouble>();
927  auto det_ptr = boost::make_shared<VectorDouble>();
928 
929  pip.push_back(new OpInvertMatrix<SPACE_DIM>(mat_F_ptr, det_ptr, inv_F));
930 
931  pip.push_back(new OpCalculatePiolaIncompressibleNH<SPACE_DIM, SPACE_DIM>(
932  shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_ptr,
933  mat_F_ptr, inv_F, det_ptr));
934  }
935 
936  auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
937  pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
938  "V", mat_v_grad_ptr));
939 
940  return boost::make_tuple(v_ptr, X_ptr, x_ptr, mat_P_ptr, mat_F_ptr, u_ptr);
941  };
942 
943  auto post_proc_boundary = [&]() {
944  auto boundary_post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
945 
946  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
947  boundary_post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
948  auto op_loop_side =
949  new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
950  // push ops to side element, through op_loop_side operator
951  auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,
952  boundary_mat_F_ptr, boundary_u_ptr] =
953  calculate_stress_ops(op_loop_side->getOpPtrVector());
954  boundary_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
955 
956  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
957 
958  boundary_post_proc_fe->getOpPtrVector().push_back(
959 
960  new OpPPMap(
961 
962  boundary_post_proc_fe->getPostProcMesh(),
963  boundary_post_proc_fe->getMapGaussPts(),
964 
965  OpPPMap::DataMapVec{},
966 
967  OpPPMap::DataMapMat{{"V", boundary_v_ptr},
968  {"GEOMETRY", boundary_X_ptr},
969  {"x", boundary_x_ptr},
970  {"U", boundary_u_ptr}},
971 
972  OpPPMap::DataMapMat{{"FIRST_PIOLA", boundary_mat_P_ptr},
973  {"F", boundary_mat_F_ptr}},
974 
975  OpPPMap::DataMapMat{}
976 
977  )
978 
979  );
980  return boundary_post_proc_fe;
981  };
982 
983  // Add monitor to time solver
984 
985  double rho = 1.;
986  CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-density", &rho, PETSC_NULL);
987  auto get_rho = [rho](const double, const double, const double) {
988  return rho;
989  };
990 
991  SmartPetscObj<Mat> M; ///< Mass matrix
992  SmartPetscObj<KSP> ksp; ///< Linear solver
993 
994  auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
995  tsPrePostProc = ts_pre_post_proc;
996 
997  CHKERR DMCreateMatrix_MoFEM(dm, M);
998  CHKERR MatZeroEntries(M);
999 
1000  boost::shared_ptr<DomainEle> vol_mass_ele(new DomainEle(mField));
1001 
1002  vol_mass_ele->B = M;
1003 
1004  auto integration_rule = [](int, int, int approx_order) {
1005  return 2 * approx_order;
1006  };
1007 
1008  vol_mass_ele->getRuleHook = integration_rule;
1009 
1010  vol_mass_ele->getOpPtrVector().push_back(new OpMassV("V", "V", get_rho));
1011  vol_mass_ele->getOpPtrVector().push_back(new OpMassF("F", "F"));
1012 
1013  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), vol_mass_ele);
1014  CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
1015  CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
1016 
1017  auto lumpVec = createDMVector(simple->getDM());
1018  CHKERR MatGetRowSum(M, lumpVec);
1019 
1020  CHKERR MatZeroEntries(M);
1021  CHKERR MatDiagonalSet(M, lumpVec, INSERT_VALUES);
1022 
1023  // Create and septup KSP (linear solver), we need this to calculate g(t,u) =
1024  // M^-1G(t,u)
1025  ksp = createKSP(mField.get_comm());
1026  CHKERR KSPSetOperators(ksp, M, M);
1027  CHKERR KSPSetFromOptions(ksp);
1028  CHKERR KSPSetUp(ksp);
1029 
1030  auto solve_boundary_for_g = [&]() {
1031  MoFEMFunctionBegin;
1032  if (*(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch)) {
1033 
1034  CHKERR VecGhostUpdateBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
1035  ADD_VALUES, SCATTER_REVERSE);
1036  CHKERR VecGhostUpdateEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
1037  ADD_VALUES, SCATTER_REVERSE);
1038  CHKERR VecAssemblyBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1039  CHKERR VecAssemblyEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1040  *(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch) = false;
1041 
1042  auto D =
1043  smartVectorDuplicate(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1044  CHKERR KSPSolve(ksp, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F, D);
1045  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1046  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1047  CHKERR VecCopy(D, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1048  }
1049 
1050  MoFEMFunctionReturn(0);
1051  };
1052 
1053  pipeline_mng->getBoundaryExplicitRhsFE()->postProcessHook =
1054  solve_boundary_for_g;
1055 
1056  MoFEM::SmartPetscObj<TS> ts;
1057  ts = pipeline_mng->createTSEX(dm);
1058 
1059  // Field eval
1060  PetscBool field_eval_flag = PETSC_TRUE;
1061  boost::shared_ptr<MatrixDouble> velocity_field_ptr;
1062  boost::shared_ptr<MatrixDouble> geometry_field_ptr;
1063  boost::shared_ptr<MatrixDouble> spatial_position_field_ptr;
1064  boost::shared_ptr<SetPtsData> field_eval_data;
1065 
1066  std::array<double, 3> field_eval_coords = {0.5, 0.5, 5.};
1067  int dim = 3;
1068  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
1069  field_eval_coords.data(), &dim,
1070  &field_eval_flag);
1071 
1072  if (field_eval_flag) {
1073  field_eval_data =
1074  mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
1075  if (SPACE_DIM == 3) {
1076  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
1077  field_eval_data, simple->getDomainFEName());
1078  } else {
1079  CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(
1080  field_eval_data, simple->getDomainFEName());
1081  }
1082 
1083  field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1084 
1085  auto no_rule = [](int, int, int) { return -1; };
1086 
1087  auto fe_ptr = field_eval_data->feMethodPtr.lock();
1088  fe_ptr->getRuleHook = no_rule;
1089  velocity_field_ptr = boost::make_shared<MatrixDouble>();
1090  geometry_field_ptr = boost::make_shared<MatrixDouble>();
1091  spatial_position_field_ptr = boost::make_shared<MatrixDouble>();
1092  fe_ptr->getOpPtrVector().push_back(
1093  new OpCalculateVectorFieldValues<SPACE_DIM>("V", velocity_field_ptr));
1094  fe_ptr->getOpPtrVector().push_back(
1095  new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY",
1096  geometry_field_ptr));
1097  fe_ptr->getOpPtrVector().push_back(
1098  new OpCalculateVectorFieldValues<SPACE_DIM>(
1099  "x_2", spatial_position_field_ptr));
1100  }
1101 
1102  auto post_proc_domain = [&]() {
1103  auto post_proc_fe_vol = boost::make_shared<PostProcEle>(mField);
1104 
1105  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
1106 
1107  auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,
1108  boundary_mat_F_ptr, boundary_u_ptr] =
1109  calculate_stress_ops(post_proc_fe_vol->getOpPtrVector());
1110 
1111  post_proc_fe_vol->getOpPtrVector().push_back(
1112 
1113  new OpPPMap(
1114 
1115  post_proc_fe_vol->getPostProcMesh(),
1116  post_proc_fe_vol->getMapGaussPts(),
1117 
1118  {},
1119 
1120  {{"V", boundary_v_ptr},
1121  {"GEOMETRY", boundary_X_ptr},
1122  {"x", boundary_x_ptr},
1123  {"U", boundary_u_ptr}},
1124 
1125  {{"FIRST_PIOLA", boundary_mat_P_ptr}, {"F", boundary_mat_F_ptr}},
1126 
1127  {}
1128 
1129  )
1130 
1131  );
1132  return post_proc_fe_vol;
1133  };
1134 
1135  boost::shared_ptr<FEMethod> null_fe;
1136  auto monitor_ptr = boost::make_shared<Monitor>(
1137  SmartPetscObj<DM>(dm, true), mField, post_proc_domain(),
1138  post_proc_boundary(), velocity_field_ptr, spatial_position_field_ptr,
1139  geometry_field_ptr, field_eval_coords, field_eval_data);
1140 
1141  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
1142  null_fe, monitor_ptr);
1143 
1144  double ftime = 1;
1145  // CHKERR TSSetMaxTime(ts, ftime);
1146  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
1147 
1148  auto T = createDMVector(simple->getDM());
1149  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
1150  SCATTER_FORWARD);
1151  CHKERR TSSetSolution(ts, T);
1152  CHKERR TSSetFromOptions(ts);
1153 
1154  CHKERR TSSetPostStage(ts, TSPrePostProc::tsPostStage);
1155  CHKERR TSSetPostStep(ts, TSPrePostProc::tsPostStep);
1156  CHKERR TSSetPreStep(ts, TSPrePostProc::tsPreStep);
1157 
1158  boost::shared_ptr<ForcesAndSourcesCore> null;
1159 
1160  if (auto ptr = tsPrePostProc.lock()) {
1161  ptr->fsRawPtr = this;
1162  CHKERR TSSetUp(ts);
1163  CHKERR TSSolve(ts, NULL);
1164  CHKERR TSGetTime(ts, &ftime);
1165  }
1166 
1167  MoFEMFunctionReturn(0);
1168 }
1169 //! [Solve]
1170 
1171 //! [Postprocess results]
1172 MoFEMErrorCode Example::outputResults() {
1173  MoFEMFunctionBegin;
1174  PetscBool test_flg = PETSC_FALSE;
1175  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
1176  if (test_flg) {
1177  auto *simple = mField.getInterface<Simple>();
1178  auto T = createDMVector(simple->getDM());
1179  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
1180  SCATTER_FORWARD);
1181  double nrm2;
1182  CHKERR VecNorm(T, NORM_2, &nrm2);
1183  MOFEM_LOG("EXAMPLE", Sev::inform) << "Regression norm " << nrm2;
1184  constexpr double regression_value = 0.0194561;
1185  if (fabs(nrm2 - regression_value) > 1e-2)
1186  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1187  "Regression test failed; wrong norm value.");
1188  }
1189  MoFEMFunctionReturn(0);
1190 }
1191 //! [Postprocess results]
1192 
1193 //! [Check]
1194 MoFEMErrorCode Example::checkResults() {
1195  MoFEMFunctionBegin;
1196  MoFEMFunctionReturn(0);
1197 }
1198 //! [Check]
1199 
1200 static char help[] = "...\n\n";
1201 
1202 int main(int argc, char *argv[]) {
1203 
1204  // Initialisation of MoFEM/PETSc and MOAB data structures
1205  const char param_file[] = "param_file.petsc";
1206  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
1207 
1208  // Add logging channel for example
1209  auto core_log = logging::core::get();
1210  core_log->add_sink(
1211  LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
1212  LogManager::setLog("EXAMPLE");
1213  MOFEM_LOG_TAG("EXAMPLE", "example");
1214 
1215  try {
1216 
1217  //! [Register MoFEM discrete manager in PETSc]
1218  DMType dm_name = "DMMOFEM";
1219  CHKERR DMRegister_MoFEM(dm_name);
1220  //! [Register MoFEM discrete manager in PETSc
1221 
1222  //! [Create MoAB]
1223  moab::Core mb_instance; ///< mesh database
1224  moab::Interface &moab = mb_instance; ///< mesh database interface
1225  //! [Create MoAB]
1226 
1227  //! [Create MoFEM]
1228  MoFEM::Core core(moab); ///< finite element database
1229  MoFEM::Interface &m_field = core; ///< finite element database interface
1230  //! [Create MoFEM]
1231 
1232  //! [Example]
1233  Example ex(m_field);
1234  CHKERR ex.runProblem();
1235  //! [Example]
1236  }
1237  CATCH_ERRORS;
1238 
1239  CHKERR MoFEM::Core::Finalize();
1240 }
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
OpCalculatePiolaIncompressibleNH::shearModulus
double shearModulus
Definition: dynamic_first_order_con_law.cpp:327
OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBodyForce
Definition: dynamic_first_order_con_law.cpp:65
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:712
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: dynamic_first_order_con_law.cpp:1194
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
young_modulus
constexpr double young_modulus
Definition: dynamic_first_order_con_law.cpp:93
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
omega
constexpr double omega
Save field DOFS on vertices/tags.
Definition: dynamic_first_order_con_law.cpp:92
OpCalculateDeformationGradient::defGradPtr
boost::shared_ptr< MatrixDouble > defGradPtr
Definition: dynamic_first_order_con_law.cpp:379
OpMassV
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMassV
Definition: dynamic_first_order_con_law.cpp:57
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:640
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
OpGradTimesTensor2
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpGradTimesTensor2
Definition: dynamic_first_order_con_law.cpp:75
OpCalculatePiolaIncompressibleNH::dEt
boost::shared_ptr< VectorDouble > dEt
Definition: dynamic_first_order_con_law.cpp:334
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
OpCalculateFStab::gradxPtr
boost::shared_ptr< MatrixDouble > gradxPtr
Definition: dynamic_first_order_con_law.cpp:161
OpCalculateFStab::OpCalculateFStab
OpCalculateFStab(boost::shared_ptr< MatrixDouble > def_grad_ptr, boost::shared_ptr< MatrixDouble > def_grad_stab_ptr, boost::shared_ptr< MatrixDouble > def_grad_dot_ptr, double tau_F_ptr, double xi_F_ptr, boost::shared_ptr< MatrixDouble > grad_x_ptr, boost::shared_ptr< MatrixDouble > grad_vel_ptr)
Definition: dynamic_first_order_con_law.cpp:104
LASTBASE
@ LASTBASE
Definition: definitions.h:69
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
Example::Example
Example(MoFEM::Interface &m_field)
Definition: dynamic_first_order_con_law.cpp:420
MoFEM::TimeScale::TimeScale
TimeScale(std::string file_name="", bool error_if_file_not_given=false, ScalingFun def_scaling_fun=[](double time) { return time;})
TimeScale constructor.
Definition: ScalingMethod.cpp:22
main
int main(int argc, char *argv[])
Definition: dynamic_first_order_con_law.cpp:1202
Monitor::fieldEvalCoords
std::array< double, 3 > fieldEvalCoords
Definition: dynamic_first_order_con_law.cpp:879
MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:913
OpMassF
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM *SPACE_DIM > OpMassF
Definition: dynamic_first_order_con_law.cpp:59
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:373
rho
double rho
Definition: plastic.cpp:144
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
OpRhsTestPiola
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpBaseTimesVector< 1, SPACE_DIM *SPACE_DIM, 1 > OpRhsTestPiola
Definition: dynamic_first_order_con_law.cpp:79
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
OpCalculateFStab::tauFPtr
double tauFPtr
Definition: dynamic_first_order_con_law.cpp:156
CommonData::ksp
SmartPetscObj< KSP > ksp
Linear solver.
Definition: dynamic_first_order_con_law.cpp:415
OpGradTimesPiola
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpGradTimesPiola
Definition: dynamic_first_order_con_law.cpp:83
OpCalculatePiolaIncompressibleNH::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: dynamic_first_order_con_law.cpp:279
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
SPACE_DIM
constexpr int SPACE_DIM
Definition: dynamic_first_order_con_law.cpp:23
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468
PostProcEleByDim< 2 >::SideEle
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
Definition: dynamic_first_order_con_law.cpp:43
TSPrePostProc::tsPostStage
static MoFEMErrorCode tsPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
[Boundary condition]
Definition: dynamic_first_order_con_law.cpp:580
OpCalculateFStab::defGradStabPtr
boost::shared_ptr< MatrixDouble > defGradStabPtr
Definition: dynamic_first_order_con_law.cpp:159
OpCalculateDisplacement::uPtr
boost::shared_ptr< MatrixDouble > uPtr
Definition: dynamic_first_order_con_law.cpp:262
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: dynamic_first_order_con_law.cpp:52
OpCalculatePiola::OpCalculatePiola
OpCalculatePiola(double shear_modulus, double bulk_modulus, double m_u, double lambda_lamme, boost::shared_ptr< MatrixDouble > first_piola_ptr, boost::shared_ptr< MatrixDouble > def_grad_ptr)
Definition: dynamic_first_order_con_law.cpp:168
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
OpCalculateDeformationGradient
Definition: dynamic_first_order_con_law.cpp:338
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
bulk_modulus_K
double bulk_modulus_K
Definition: dynamic_first_order_con_law.cpp:95
MoFEM::DMoFEMMeshToLocalVector
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:523
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
OpCalculatePiola::defGradPtr
boost::shared_ptr< MatrixDouble > defGradPtr
Definition: dynamic_first_order_con_law.cpp:221
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:261
tsPrePostProc
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
Definition: dynamic_first_order_con_law.cpp:404
Monitor::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: dynamic_first_order_con_law.cpp:789
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
OpCalculatePiolaIncompressibleNH
Definition: dynamic_first_order_con_law.cpp:266
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
EntData
EntitiesFieldData::EntData EntData
Definition: dynamic_first_order_con_law.cpp:26
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
OpCalculateDisplacement::XPtr
boost::shared_ptr< MatrixDouble > XPtr
Definition: dynamic_first_order_con_law.cpp:261
poisson_ratio
constexpr double poisson_ratio
Definition: dynamic_first_order_con_law.cpp:94
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: dynamic_first_order_con_law.cpp:462
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: dynamic_first_order_con_law.cpp:62
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:677
PostProcEleByDim
Definition: dynamic_first_order_con_law.cpp:38
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
trace
double trace(FTensor::Tensor2< T, 2, 2 > &t_stress)
Definition: dynamic_first_order_con_law.cpp:14
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1197
TSPrePostProc
Set of functions called by PETSc solver used to refine and update mesh.
Definition: dynamic_first_order_con_law.cpp:383
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
OpCalculateDeformationGradient::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: dynamic_first_order_con_law.cpp:346
OpCalculateDisplacement::OpCalculateDisplacement
OpCalculateDisplacement(boost::shared_ptr< MatrixDouble > spatial_pos_ptr, boost::shared_ptr< MatrixDouble > reference_pos_ptr, boost::shared_ptr< MatrixDouble > u_ptr)
Definition: dynamic_first_order_con_law.cpp:226
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
TSPrePostProc::tsPostStep
static MoFEMErrorCode tsPostStep(TS ts)
Definition: dynamic_first_order_con_law.cpp:611
Monitor::velocityFieldPtr
boost::shared_ptr< MatrixDouble > velocityFieldPtr
Definition: dynamic_first_order_con_law.cpp:876
Monitor::mField
MoFEM::Interface & mField
Definition: dynamic_first_order_con_law.cpp:775
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpCalculatePiola::bulkModulus
double bulkModulus
Definition: dynamic_first_order_con_law.cpp:217
OpCalculateDisplacement
Definition: dynamic_first_order_con_law.cpp:225
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
BoundaryEleOp
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:152
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::FieldEvaluatorInterface::SetPtsData
Definition: FieldEvaluator.hpp:29
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Example
[Example]
Definition: plastic.cpp:216
double
convert.type
type
Definition: convert.py:64
DomainEle
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomainEle DomainEle
Definition: dynamic_first_order_con_law.cpp:27
CommonData::M
SmartPetscObj< Mat > M
Mass matrix.
Definition: dynamic_first_order_con_law.cpp:414
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::TimeScale
Force scale operator for reading two columns.
Definition: ScalingMethod.hpp:32
OpCalculateFStab
Definition: dynamic_first_order_con_law.cpp:103
Monitor::fieldEvalData
boost::shared_ptr< SetPtsData > fieldEvalData
Definition: dynamic_first_order_con_law.cpp:880
help
static char help[]
[Check]
Definition: dynamic_first_order_con_law.cpp:1200
OpCalculateDisplacement::xPtr
boost::shared_ptr< MatrixDouble > xPtr
Definition: dynamic_first_order_con_law.cpp:260
OpCalculatePiola::mU
double mU
Definition: dynamic_first_order_con_law.cpp:218
MoFEM::PostProcBrokenMeshInMoabBaseContImpl
Definition: PostProcBrokenMeshInMoabBase.hpp:890
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: dynamic_first_order_con_law.cpp:884
OpCalculatePiola
Definition: dynamic_first_order_con_law.cpp:167
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
OpCalculatePiolaIncompressibleNH::invDefGradPtr
boost::shared_ptr< MatrixDouble > invDefGradPtr
Definition: dynamic_first_order_con_law.cpp:333
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
SetPtsData
FieldEvaluatorInterface::SetPtsData SetPtsData
Definition: dynamic_first_order_con_law.cpp:36
OpCalculatePiola::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: dynamic_first_order_con_law.cpp:177
MatrixFunction.hpp
MoFEM::smartVectorDuplicate
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
Definition: PetscSmartObj.hpp:230
Example::DynamicFirstOrderConsConstantTimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: dynamic_first_order_con_law.cpp:443
Monitor::dM
SmartPetscObj< DM > dM
Definition: dynamic_first_order_con_law.cpp:873
OpCalculatePiolaIncompressibleNH::defGradPtr
boost::shared_ptr< MatrixDouble > defGradPtr
Definition: dynamic_first_order_con_law.cpp:332
lamme_lambda
double lamme_lambda
Definition: dynamic_first_order_con_law.cpp:98
Monitor::postProcBdy
boost::shared_ptr< PostProcFaceEle > postProcBdy
Definition: dynamic_first_order_con_law.cpp:875
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
LinMomTimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: dynamic_first_order_con_law.cpp:408
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
Example::DynamicFirstOrderConsSinusTimeScale
Definition: dynamic_first_order_con_law.cpp:436
OpCalculatePiolaIncompressibleNH::OpCalculatePiolaIncompressibleNH
OpCalculatePiolaIncompressibleNH(double shear_modulus, double bulk_modulus, double m_u, double lambda_lamme, boost::shared_ptr< MatrixDouble > first_piola_ptr, boost::shared_ptr< MatrixDouble > def_grad_ptr, boost::shared_ptr< MatrixDouble > inv_def_grad_ptr, boost::shared_ptr< VectorDouble > det)
Definition: dynamic_first_order_con_law.cpp:268
OpCalculatePiolaIncompressibleNH::bulkModulus
double bulkModulus
Definition: dynamic_first_order_con_law.cpp:328
PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
Definition: dynamic_first_order_con_law.cpp:54
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
Example::DynamicFirstOrderConsSinusTimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: dynamic_first_order_con_law.cpp:438
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:139
BiLinearForm
Example::DynamicFirstOrderConsConstantTimeScale
Definition: dynamic_first_order_con_law.cpp:441
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
MoFEM::TimeScale::getScale
double getScale(const double time)
Get scaling at a given time.
Definition: ScalingMethod.cpp:137
FTensor::Index< 'i', SPACE_DIM >
OpCalculateFStab::xiF
double xiF
Definition: dynamic_first_order_con_law.cpp:157
OpCalculateFStab::defGradPtr
boost::shared_ptr< MatrixDouble > defGradPtr
Definition: dynamic_first_order_con_law.cpp:158
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:417
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1561
OpCalculateDeformationGradient::gradTensorPtr
boost::shared_ptr< MatrixDouble > gradTensorPtr
Definition: dynamic_first_order_con_law.cpp:378
OpCalculatePiolaIncompressibleNH::lammeLambda
double lammeLambda
Definition: dynamic_first_order_con_law.cpp:330
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
integration_rule
auto integration_rule
Definition: free_surface.cpp:187
OpCalculateDisplacement::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: dynamic_first_order_con_law.cpp:232
Range
DomainEleOp
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
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
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:96
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
Monitor::x2FieldPtr
boost::shared_ptr< MatrixDouble > x2FieldPtr
Definition: dynamic_first_order_con_law.cpp:877
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:275
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: dynamic_first_order_con_law.cpp:535
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
Monitor
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:774
approx_order
int approx_order
Definition: test_broken_space.cpp:54
OpCalculateDeformationGradient::OpCalculateDeformationGradient
OpCalculateDeformationGradient(boost::shared_ptr< MatrixDouble > def_grad_ptr, boost::shared_ptr< MatrixDouble > grad_tensor_ptr)
Definition: dynamic_first_order_con_law.cpp:340
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
MoFEM::DMMoFEMTSSetMonitor
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:1056
mu
double mu
Definition: dynamic_first_order_con_law.cpp:97
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:64
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
Monitor::postProc
boost::shared_ptr< PostProcEle > postProc
Definition: dynamic_first_order_con_law.cpp:874
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3684
OpCalculateFStab::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: dynamic_first_order_con_law.cpp:115
OpCalculatePiolaIncompressibleNH::mU
double mU
Definition: dynamic_first_order_con_law.cpp:329
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:256
OpCalculatePiola::firstPiolaPtr
boost::shared_ptr< MatrixDouble > firstPiolaPtr
Definition: dynamic_first_order_con_law.cpp:220
TSPrePostProc::tsPreStep
static MoFEMErrorCode tsPreStep(TS ts)
Definition: dynamic_first_order_con_law.cpp:624
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
Example::mField
MoFEM::Interface & mField
Definition: plastic.cpp:226
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
Monitor::Monitor
Monitor(SmartPetscObj< DM > dm, MoFEM::Interface &m_field, boost::shared_ptr< PostProcEle > post_proc, boost::shared_ptr< PostProcFaceEle > post_proc_bdry, boost::shared_ptr< MatrixDouble > velocity_field_ptr, boost::shared_ptr< MatrixDouble > x2_field_ptr, boost::shared_ptr< MatrixDouble > geometry_field_ptr, std::array< double, 3 > pass_field_eval_coords, boost::shared_ptr< SetPtsData > pass_field_eval_data)
Definition: dynamic_first_order_con_law.cpp:776
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
OpCalculatePiolaIncompressibleNH::firstPiolaPtr
boost::shared_ptr< MatrixDouble > firstPiolaPtr
Definition: dynamic_first_order_con_law.cpp:331
OpCalculatePiola::shearModulus
double shearModulus
Definition: dynamic_first_order_con_law.cpp:216
LinMomTimeScale
Definition: dynamic_first_order_con_law.cpp:406
QUIET
@ QUIET
Definition: definitions.h:221
OpCalculatePiola::lammeLambda
double lammeLambda
Definition: dynamic_first_order_con_law.cpp:219
third
constexpr double third
Definition: EshelbianADOL-C.cpp:17
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:589
dt
double dt
Definition: heat_method.cpp:26
MoFEM::SmartPetscObj< Mat >
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition: dynamic_first_order_con_law.cpp:1172
MoFEM::DMoFEMLoopFiniteElements
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:586
OpCalculateFStab::defGradDotPtr
boost::shared_ptr< MatrixDouble > defGradDotPtr
Definition: dynamic_first_order_con_law.cpp:160
MF_EXIST
@ MF_EXIST
Definition: definitions.h:113
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
PostProcEleByDim< 3 >::SideEle
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
Definition: dynamic_first_order_con_law.cpp:49
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1291
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
Monitor::geometryFieldPtr
boost::shared_ptr< MatrixDouble > geometryFieldPtr
Definition: dynamic_first_order_con_law.cpp:878
TSPrePostProc::fsRawPtr
Example * fsRawPtr
Definition: dynamic_first_order_con_law.cpp:397
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:711
OpCalculateFStab::gradVelPtr
boost::shared_ptr< MatrixDouble > gradVelPtr
Definition: dynamic_first_order_con_law.cpp:162
MoFEM::PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:687
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:709
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow