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