v0.14.0
shallow_wave.cpp
Go to the documentation of this file.
1 /**
2  * \file shallow_wave.cpp
3  * \example shallow_wave.cpp
4  *
5  * Solving shallow wave equation on manifold
6  *
7  * The inital conditions are set following this paper \cite scott2016test.
8  *
9  */
10 
11 #include <MoFEM.hpp>
12 using namespace MoFEM;
13 
14 #include <boost/math/quadrature/gauss_kronrod.hpp>
15 using namespace boost::math::quadrature;
16 
17 template <int DIM> struct ElementsAndOps {};
18 
19 template <> struct ElementsAndOps<2> {
20  using DomainEle = FaceElementForcesAndSourcesCore;
21  using DomainEleOp = DomainEle::UserDataOperator;
22 };
23 
24 constexpr int FE_DIM = 2;
25 
26 using EntData = EntitiesFieldData::EntData;
27 using DomainEle = ElementsAndOps<FE_DIM>::DomainEle;
28 using DomainEleOp = ElementsAndOps<FE_DIM>::DomainEleOp;
29 using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
30 
31 using AssemblyDomainEleOp =
32  FormsIntegrators<DomainEleOp>::Assembly<PETSC>::OpBase;
33 
34 // Use forms iterators for Grad-Grad term
35 using OpBaseDivU = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
36  GAUSS>::OpMixScalarTimesDiv<3>;
37 using OpBaseGradH = FormsIntegrators<DomainEleOp>::Assembly<
38  PETSC>::BiLinearForm<GAUSS>::OpMixVectorTimesGrad<1, 3, 3>;
39 
40 // Use forms for Mass term
41 using OpMassUU = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
42  GAUSS>::OpMass<1, 3>;
43 using OpMassHH = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
44  GAUSS>::OpMass<1, 1>;
45 
46 using OpBaseTimesDotH = FormsIntegrators<DomainEleOp>::Assembly<
47  PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
48 using OpBaseTimesDivU = OpBaseTimesDotH;
49 
50 using OpSourceU = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
51  GAUSS>::OpSource<1, 3>;
52 using OpSourceH = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
53  GAUSS>::OpSource<1, 1>;
54 
55 using OpConvectiveH = FormsIntegrators<DomainEleOp>::Assembly<
56  PETSC>::LinearForm<GAUSS>::OpConvectiveTermRhs<1, 1, 3>;
57 using OpConvectiveH_dU = FormsIntegrators<DomainEleOp>::Assembly<
58  PETSC>::TriLinearForm<GAUSS>::OpConvectiveTermLhsDu<1, 1, 3>;
59 using OpConvectiveH_dGradH = FormsIntegrators<DomainEleOp>::Assembly<
60  PETSC>::TriLinearForm<GAUSS>::OpConvectiveTermLhsDy<1, 1, 3>;
61 
62 constexpr double omega = 7.292 * 1e-5;
63 constexpr double g = 9.80616;
64 constexpr double mu = 1e4;
65 constexpr double h0 = 1e4;
66 
67 constexpr double h_hat = 120;
68 constexpr double u_max = 80;
69 constexpr double phi_0 = M_PI / 7;
70 constexpr double phi_1 = M_PI / 2 - phi_0;
71 constexpr double phi_2 = M_PI / 4;
72 constexpr double alpha_montain = 1. / 3.;
73 constexpr double beta_montain = 1. / 15.;
74 
75 constexpr double penalty = 1;
76 
77 FTensor::Index<'i', 3> i;
78 FTensor::Index<'j', 3> j;
79 FTensor::Index<'l', 3> l;
80 FTensor::Index<'m', 3> m;
81 
82 struct OpURhs : public AssemblyDomainEleOp {
83 
84  OpURhs(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
85  boost::shared_ptr<MatrixDouble> u_dot_ptr,
86  boost::shared_ptr<MatrixDouble> grad_u_ptr,
87  boost::shared_ptr<MatrixDouble> grad_h_ptr)
88  : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
89  uPtr(u_ptr), uDotPtr(u_dot_ptr), uGradPtr(grad_u_ptr),
90  hGradPtr(grad_h_ptr) {}
91 
92  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
93  MoFEMFunctionBegin;
94 
95  const double vol = getMeasure();
96  auto t_w = getFTensor0IntegrationWeight();
97  auto t_row_base = row_data.getFTensor0N();
98  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
99  auto t_dot_u = getFTensor1FromMat<3>(*uDotPtr);
100  auto t_u = getFTensor1FromMat<3>(*uPtr);
101  auto t_grad_u = getFTensor2FromMat<3, 3>(*uGradPtr);
102  auto t_grad_h = getFTensor1FromMat<3>(*hGradPtr);
103  auto t_coords = getFTensor1CoordsAtGaussPts();
104  auto t_normal = getFTensor1NormalsAtGaussPts();
105 
106  for (int gg = 0; gg != nbIntegrationPts; gg++) {
107 
108  const double alpha = t_w * vol;
109  auto t_nf = getFTensor1FromArray<3, 3>(locF);
110 
111  FTensor::Tensor1<double, 3> t_rhs;
112  FTensor::Tensor2<double, 3, 3> t_rhs_grad;
113 
114  const auto a = std::sqrt(t_coords(i) * t_coords(i));
115  const auto sin_fi = t_coords(2) / a;
116  const auto f = 2 * omega * sin_fi;
117 
118  FTensor::Tensor1<double, 3> t_r;
119  t_r(i) = t_normal(i);
120  t_r.normalize();
121 
122  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
123  FTensor::Tensor2<double, 3, 3> t_P, t_Q;
124  t_P(i, j) = t_r(i) * t_r(j);
125  t_Q(i, j) = t_kd(i, j) - t_P(i, j);
126 
127  FTensor::Tensor2<double, 3, 3> t_A;
128  t_A(i, m) = levi_civita(i, j, m) * t_r(j);
129 
130  t_rhs(m) = t_Q(m, i) * (t_dot_u(i) + t_grad_u(i, j) * t_u(j) +
131  f * t_A(i, j) * t_u(j) + g * t_grad_h(i));
132  t_rhs_grad(m, j) = t_Q(m, i) * (mu * t_grad_u(i, j));
133 
134  t_rhs(m) += t_P(m, j) * t_u(j);
135 
136  int rr = 0;
137  for (; rr != nbRows / 3; ++rr) {
138  t_nf(i) += alpha * t_row_base * t_rhs(i);
139  t_nf(i) += alpha * t_row_diff_base(j) * t_rhs_grad(i, j);
140  ++t_row_base;
141  ++t_row_diff_base;
142  ++t_nf;
143  }
144  for (; rr < nbRowBaseFunctions; ++rr) {
145  ++t_row_base;
146  ++t_row_diff_base;
147  }
148 
149  ++t_w;
150  ++t_u;
151  ++t_dot_u;
152  ++t_grad_u;
153  ++t_grad_h;
154  ++t_coords;
155  ++t_normal;
156  }
157 
158  MoFEMFunctionReturn(0);
159  }
160 
161 private:
162  boost::shared_ptr<MatrixDouble> uPtr;
163  boost::shared_ptr<MatrixDouble> uDotPtr;
164  boost::shared_ptr<MatrixDouble> uGradPtr;
165  boost::shared_ptr<MatrixDouble> hGradPtr;
166 };
167 
168 struct OpULhs_dU : public AssemblyDomainEleOp {
169 
170  OpULhs_dU(const std::string field_name_row, const std::string field_name_col,
171  boost::shared_ptr<MatrixDouble> u_ptr,
172  boost::shared_ptr<MatrixDouble> grad_u_ptr)
173  : AssemblyDomainEleOp(field_name_row, field_name_col,
174  AssemblyDomainEleOp::OPROWCOL),
175  uPtr(u_ptr), uGradPtr(grad_u_ptr) {
176  this->sYmm = false;
177  }
178 
179  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
180  EntitiesFieldData::EntData &col_data) {
181  MoFEMFunctionBegin;
182 
183  const double vol = getMeasure();
184  auto t_w = getFTensor0IntegrationWeight();
185  auto t_row_base = row_data.getFTensor0N();
186  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
187  auto t_coords = getFTensor1CoordsAtGaussPts();
188  auto t_normal = getFTensor1NormalsAtGaussPts();
189 
190  auto t_u = getFTensor1FromMat<3>(*uPtr);
191  auto t_grad_u = getFTensor2FromMat<3, 3>(*uGradPtr);
192 
193  auto get_t_mat = [&](const int rr) {
194  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>{
195  &locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
196 
197  &locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
198 
199  &locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
200  };
201 
202  const auto ts_a = getFEMethod()->ts_a;
203 
204  for (int gg = 0; gg != nbIntegrationPts; gg++) {
205 
206  const auto a = std::sqrt(t_coords(i) * t_coords(i));
207  const auto sin_fi = t_coords(2) / a;
208  const auto f = 2 * omega * sin_fi;
209 
210  FTensor::Tensor1<double, 3> t_r;
211  t_r(i) = t_normal(i);
212  t_r.normalize();
213 
214  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
215  FTensor::Tensor2<double, 3, 3> t_P, t_Q;
216  t_P(i, j) = t_r(i) * t_r(j);
217  t_Q(i, j) = t_kd(i, j) - t_P(i, j);
218 
219  FTensor::Tensor2<double, 3, 3> t_A;
220  t_A(i, m) = levi_civita(i, j, m) * t_r(j);
221 
222  FTensor::Tensor2<double, 3, 3> t_rhs_du;
223  t_rhs_du(m, j) =
224  t_Q(m, i) * (ts_a * t_kd(i, j) + t_grad_u(i, j) + f * t_A(i, j)) +
225  t_P(m, j);
226 
227  const double alpha = t_w * vol;
228 
229  int rr = 0;
230  for (; rr != nbRows / 3; rr++) {
231 
232  auto t_col_base = col_data.getFTensor0N(gg, 0);
233  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
234  auto t_mat = get_t_mat(3 * rr);
235 
236  for (int cc = 0; cc != nbCols / 3; cc++) {
237  t_mat(i, j) += (alpha * t_row_base * t_col_base) * t_rhs_du(i, j);
238  t_mat(i, j) += (alpha * mu) * t_Q(i, j) *
239  (t_row_diff_base(m) * t_col_diff_base(m));
240  ++t_col_diff_base;
241  ++t_col_base;
242  ++t_mat;
243  }
244  ++t_row_base;
245  ++t_row_diff_base;
246  }
247  for (; rr < nbRowBaseFunctions; ++rr) {
248  ++t_row_base;
249  ++t_row_diff_base;
250  }
251 
252  ++t_w;
253  ++t_coords;
254  ++t_normal;
255  ++t_u;
256  ++t_grad_u;
257  }
258 
259  MoFEMFunctionReturn(0);
260  }
261 
262 private:
263  boost::shared_ptr<MatrixDouble> uPtr;
264  boost::shared_ptr<MatrixDouble> uGradPtr;
265 };
266 
267 struct OpULhs_dH : public AssemblyDomainEleOp {
268 
269  OpULhs_dH(const std::string field_name_row, const std::string field_name_col)
270  : AssemblyDomainEleOp(field_name_row, field_name_col,
271  AssemblyDomainEleOp::OPROWCOL) {
272  this->sYmm = false;
273  }
274 
275  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
276  EntitiesFieldData::EntData &col_data) {
277  MoFEMFunctionBegin;
278 
279  // get element volume
280  const double vol = getMeasure();
281  // get integration weights
282  auto t_w = getFTensor0IntegrationWeight();
283  // get base function gradient on rows
284  auto t_row_base = row_data.getFTensor0N();
285  // normal
286  auto t_normal = getFTensor1NormalsAtGaussPts();
287 
288  auto get_t_vec = [&](const int rr) {
289  return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>{
290  &locMat(rr + 0, 0),
291 
292  &locMat(rr + 1, 0),
293 
294  &locMat(rr + 2, 0)};
295  };
296 
297  // loop over integration points
298  for (int gg = 0; gg != nbIntegrationPts; gg++) {
299 
300  FTensor::Tensor1<double, 3> t_r;
301  t_r(i) = t_normal(i);
302  t_r.normalize();
303 
304  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
305  FTensor::Tensor2<double, 3, 3> t_P, t_Q;
306  t_P(i, j) = t_r(i) * t_r(j);
307  t_Q(i, j) = t_kd(i, j) - t_P(i, j);
308 
309  const double alpha = t_w * vol;
310 
311  int rr = 0;
312  for (; rr != nbRows / 3; rr++) {
313  auto t_vec = get_t_vec(3 * rr);
314  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
315  const double a = alpha * g * t_row_base;
316 
317  for (int cc = 0; cc != nbCols; cc++) {
318  t_vec(i) += a * (t_Q(i, m) * t_col_diff_base(m));
319  ++t_vec;
320  ++t_col_diff_base;
321  }
322 
323  ++t_row_base;
324  }
325 
326  for (; rr < nbRowBaseFunctions; ++rr)
327  ++t_row_base;
328 
329  ++t_w;
330  ++t_normal;
331  }
332 
333  MoFEMFunctionReturn(0);
334  }
335 };
336 
337 struct Example {
338 
339  Example(MoFEM::Interface &m_field) : mField(m_field) {}
340 
341  MoFEMErrorCode runProblem();
342 
343 private:
344  MoFEM::Interface &mField;
345 
346  MoFEMErrorCode readMesh();
347  MoFEMErrorCode setupProblem();
348  MoFEMErrorCode createCommonData();
349  MoFEMErrorCode boundaryCondition();
350  MoFEMErrorCode assembleSystem();
351  MoFEMErrorCode solveSystem();
352  MoFEMErrorCode outputResults();
353  MoFEMErrorCode checkResults();
354 
355  boost::shared_ptr<FEMethod> domianLhsFEPtr;
356  boost::shared_ptr<FEMethod> domianRhsFEPtr;
357 };
358 
359 //! [Create common data]
360 MoFEMErrorCode Example::createCommonData() {
361  MoFEMFunctionBegin;
362 
363  MoFEMFunctionReturn(0);
364 }
365 //! [Create common data]
366 
367 //! [Run problem]
368 MoFEMErrorCode Example::runProblem() {
369  MoFEMFunctionBegin;
370  CHKERR readMesh();
371  CHKERR setupProblem();
372  CHKERR createCommonData();
373  CHKERR boundaryCondition();
374  CHKERR assembleSystem();
375  CHKERR solveSystem();
376  CHKERR outputResults();
377  CHKERR checkResults();
378  MoFEMFunctionReturn(0);
379 }
380 //! [Run problem]
381 
382 //! [Read mesh]
383 MoFEMErrorCode Example::readMesh() {
384  MoFEMFunctionBegin;
385  auto simple = mField.getInterface<Simple>();
386  CHKERR simple->getOptions();
387  CHKERR simple->loadFile();
388  MoFEMFunctionReturn(0);
389 }
390 //! [Read mesh]
391 
392 //! [Set up problem]
393 MoFEMErrorCode Example::setupProblem() {
394  MoFEMFunctionBegin;
395  Simple *simple = mField.getInterface<Simple>();
396  // Add field
397  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 3);
398  CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
399  CHKERR simple->addDataField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
400 
401  int order = 2;
402  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
403  CHKERR simple->setFieldOrder("U", order);
404  CHKERR simple->setFieldOrder("H", order);
405  CHKERR simple->setFieldOrder("HO_POSITIONS", 3);
406 
407  CHKERR simple->setUp();
408  MoFEMFunctionReturn(0);
409 }
410 //! [Set up problem]
411 
412 //! [Boundary condition]
413 MoFEMErrorCode Example::boundaryCondition() {
414  MoFEMFunctionBegin;
415 
416  PetscBool is_restart = PETSC_FALSE;
417  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_restart", &is_restart,
418  PETSC_NULL);
419 
420  auto restart_vector = [&]() {
421  MoFEMFunctionBegin;
422  auto simple = mField.getInterface<Simple>();
423  auto dm = simple->getDM();
424  MOFEM_LOG("SW", Sev::inform)
425  << "reading vector in binary from vector.dat ...";
426  PetscViewer viewer;
427  PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_READ,
428  &viewer);
429  auto T = createDMVector(simple->getDM());
430  VecLoad(T, viewer);
431  CHKERR DMoFEMMeshToLocalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
432  MoFEMFunctionReturn(0);
433  };
434 
435  if (is_restart) {
436 
437  CHKERR restart_vector();
438 
439  } else {
440 
441  const double e_n = exp(-4 / pow(phi_1 - phi_0, 2));
442  const double u0 = u_max / e_n;
443 
444  FTensor::Tensor1<double, 3> t_k{0., 0., 1.};
445  FTensor::Tensor2<double, 3, 3> t_A;
446  t_A(i, m) = levi_civita(i, j, m) * t_k(j);
447 
448  auto get_phi = [&](const double x, const double y, const double z) {
449  FTensor::Tensor1<double, 3> t_r{x, y, 0.};
450  const double r = std::sqrt(t_r(i) * t_r(i));
451  return atan2(z, r);
452  };
453 
454  auto init_u_phi = [&](const double phi) {
455  if (phi > phi_0 && phi < phi_1) {
456  return u0 * exp(1. / ((phi - phi_0) * (phi - phi_1)));
457  } else {
458  return 0.;
459  }
460  };
461 
462  auto init_u = [&](const double x, const double y, const double z) {
463  FTensor::Tensor1<double, 3> t_u{0., 0., 0.};
464  const double u_phi = init_u_phi(get_phi(x, y, z));
465  if (u_phi > 0) {
466  FTensor::Tensor1<double, 3> t_r{x, y, 0.};
467  t_r.normalize();
468  t_u(i) = ((t_A(i, j) * t_r(j)) * u_phi);
469  }
470  return t_u;
471  };
472 
473  auto init_h = [&](const double x, const double y, const double z) {
474  const double a = std::sqrt(x * x + y * y + z * z);
475 
476  auto integral = [&](const double fi) {
477  const double u_phi = init_u_phi(fi);
478  const auto f = 2 * omega * sin(fi);
479  return a * u_phi * (f + (tan(fi) / a) * u_phi);
480  };
481 
482  auto montain = [&](const double lambda, const double fi) {
483  if (lambda > -M_PI && lambda < M_PI)
484  return h_hat * cos(fi) * exp(-pow(lambda / alpha_montain, 2)) *
485  exp(-pow((phi_2 - fi) / beta_montain, 2));
486  else
487  return 0.;
488  };
489 
490  const double fi = get_phi(x, y, z);
491  const double lambda = atan2(x, y);
492 
493  double h1 = 0;
494  if (fi > phi_0)
495  h1 = gauss_kronrod<double, 32>::integrate(
496  integral, phi_0, fi, 0, std::numeric_limits<float>::epsilon());
497 
498  return h0 + (montain(lambda, fi) - (h1 / g));
499  };
500 
501  auto set_domain_general = [&](auto &pipeline) {
502  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
503  pipeline.push_back(new OpGetHONormalsOnFace("HO_POSITIONS"));
504  pipeline.push_back(new OpCalculateHOCoords("HO_POSITIONS"));
505  pipeline.push_back(new OpSetHOWeightsOnFace());
506  };
507 
508  auto set_domain_rhs = [&](auto &pipeline) {
509  pipeline.push_back(new OpSourceU("U", init_u));
510  pipeline.push_back(new OpSourceH("H", init_h));
511  };
512 
513  auto set_domain_lhs = [&](auto &pipeline) {
514  pipeline.push_back(
515  new OpMassUU("U", "U", [](double, double, double) { return 1; }));
516  pipeline.push_back(
517  new OpMassHH("H", "H", [](double, double, double) { return 1; }));
518  };
519 
520  auto post_proc = [&]() {
521  MoFEMFunctionBegin;
522  auto simple = mField.getInterface<Simple>();
523  auto dm = simple->getDM();
524 
525  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
526 
527  auto det_ptr = boost::make_shared<VectorDouble>();
528  auto jac_ptr = boost::make_shared<MatrixDouble>();
529  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
530 
531  post_proc_fe->getOpPtrVector().push_back(
532  new OpGetHONormalsOnFace("HO_POSITIONS"));
533  post_proc_fe->getOpPtrVector().push_back(
534  new OpCalculateHOCoords("HO_POSITIONS"));
535  post_proc_fe->getOpPtrVector().push_back(
536  new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
537  post_proc_fe->getOpPtrVector().push_back(
538  new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
539 
540  auto u_ptr = boost::make_shared<MatrixDouble>();
541  auto h_ptr = boost::make_shared<VectorDouble>();
542  auto pos_ptr = boost::make_shared<MatrixDouble>();
543 
544  post_proc_fe->getOpPtrVector().push_back(
545  new OpCalculateVectorFieldValues<3>("U", u_ptr));
546  post_proc_fe->getOpPtrVector().push_back(
547  new OpCalculateScalarFieldValues("H", h_ptr));
548  post_proc_fe->getOpPtrVector().push_back(
549  new OpCalculateVectorFieldValues<3>("HO_POSITIONS", pos_ptr));
550 
551  using OpPPMap = OpPostProcMapInMoab<3, 3>;
552 
553  post_proc_fe->getOpPtrVector().push_back(
554 
555  new OpPPMap(post_proc_fe->getPostProcMesh(),
556  post_proc_fe->getMapGaussPts(),
557 
558  {{"H", h_ptr}},
559 
560  {{"U", u_ptr}, {"HO_POSITIONS", pos_ptr}},
561 
562  {}, {}
563 
564  )
565 
566  );
567 
568  CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
569  CHKERR post_proc_fe->writeFile("out_init.h5m");
570 
571  MoFEMFunctionReturn(0);
572  };
573 
574  auto solve_init = [&]() {
575  MoFEMFunctionBegin;
576  auto simple = mField.getInterface<Simple>();
577  auto pipeline_mng = mField.getInterface<PipelineManager>();
578 
579  auto solver = pipeline_mng->createKSP();
580  CHKERR KSPSetFromOptions(solver);
581  PC pc;
582  CHKERR KSPGetPC(solver, &pc);
583  PetscBool is_pcfs = PETSC_FALSE;
584  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
585  if (is_pcfs == PETSC_TRUE) {
586  auto bc_mng = mField.getInterface<BcManager>();
587  auto name_prb = simple->getProblemName();
588  SmartPetscObj<IS> is_u;
589  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
590  name_prb, ROW, "U", 0, 3, is_u);
591  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
592  CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_ADDITIVE);
593  }
594 
595  CHKERR KSPSetUp(solver);
596 
597  auto dm = simple->getDM();
598  auto D = createDMVector(dm);
599  auto F = vectorDuplicate(D);
600 
601  CHKERR KSPSolve(solver, F, D);
602  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
603  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
604  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
605 
606  MoFEMFunctionReturn(0);
607  };
608 
609  auto pipeline_mng = mField.getInterface<PipelineManager>();
610 
611  auto integration_rule = [](int, int, int approx_order) {
612  return 2 * approx_order + 4;
613  };
614  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
615  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
616 
617  set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
618  set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
619  set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
620  set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
621 
622  CHKERR solve_init();
623  CHKERR post_proc();
624 
625  // Clear pipelines
626  pipeline_mng->getOpDomainRhsPipeline().clear();
627  pipeline_mng->getOpDomainLhsPipeline().clear();
628  }
629 
630  MoFEMFunctionReturn(0);
631 }
632 //! [Boundary condition]
633 
634 //! [Push operators to pipeline]
635 MoFEMErrorCode Example::assembleSystem() {
636  MoFEMFunctionBegin;
637 
638  // Push element from reference configuration to current configuration in 3d
639  // space
640  auto set_domain_general = [&](auto &pipeline) {
641  auto det_ptr = boost::make_shared<VectorDouble>();
642  auto jac_ptr = boost::make_shared<MatrixDouble>();
643  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
644  pipeline.push_back(new OpGetHONormalsOnFace("HO_POSITIONS"));
645  pipeline.push_back(new OpCalculateHOCoords("HO_POSITIONS"));
646  pipeline.push_back(new OpSetHOWeightsOnFace());
647  pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
648  pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
649  pipeline.push_back(new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
650  };
651 
652  auto set_domain_rhs = [&](auto &pipeline) {
653  auto dot_u_ptr = boost::make_shared<MatrixDouble>();
654  auto u_ptr = boost::make_shared<MatrixDouble>();
655  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
656  auto div_u_ptr = boost::make_shared<VectorDouble>();
657  auto dot_h_ptr = boost::make_shared<VectorDouble>();
658  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
659 
660  pipeline.push_back(new OpCalculateVectorFieldValuesDot<3>("U", dot_u_ptr));
661  pipeline.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
662 
663  pipeline.push_back(new OpCalculateVectorFieldValues<3>("U", u_ptr));
664  pipeline.push_back(
665  new OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
666  pipeline.push_back(
667  new OpCalculateDivergenceVectorFieldValues<3>("U", div_u_ptr));
668  pipeline.push_back(new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
669 
670  pipeline.push_back(new OpBaseTimesDotH(
671  "H", dot_h_ptr, [](double, double, double) { return 1.; }));
672  pipeline.push_back(new OpBaseTimesDivU(
673  "H", div_u_ptr, [](double, double, double) { return h0; }));
674  pipeline.push_back(new OpConvectiveH("H", u_ptr, grad_h_ptr));
675  pipeline.push_back(
676  new OpURhs("U", u_ptr, dot_u_ptr, grad_u_ptr, grad_h_ptr));
677  };
678 
679  auto set_domain_lhs = [&](auto &pipeline) {
680  auto u_ptr = boost::make_shared<MatrixDouble>();
681  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
682  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
683 
684  pipeline.push_back(new OpCalculateVectorFieldValues<3>("U", u_ptr));
685  pipeline.push_back(
686  new OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
687  pipeline.push_back(new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
688 
689  pipeline.push_back(new OpMassHH("H", "H", [&](double, double, double) {
690  return domianLhsFEPtr->ts_a;
691  }));
692  pipeline.push_back(new OpBaseDivU(
693  "H", "U", [](const double, const double, const double) { return h0; },
694  false, false));
695  pipeline.push_back(
696  new OpConvectiveH_dU("H", "U", grad_h_ptr, []() { return 1; }));
697  pipeline.push_back(
698  new OpConvectiveH_dGradH("H", "H", u_ptr, []() { return 1; }));
699  pipeline.push_back(new OpULhs_dU("U", "U", u_ptr, grad_u_ptr));
700  pipeline.push_back(new OpULhs_dH("U", "H"));
701  };
702 
703  auto *pipeline_mng = mField.getInterface<PipelineManager>();
704 
705  auto integration_rule = [](int, int, int approx_order) {
706  return 2 * approx_order + 4;
707  };
708  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
709  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
710 
711  set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
712  set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
713 
714  set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
715  set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
716 
717  domianLhsFEPtr = pipeline_mng->getDomainLhsFE();
718  domianRhsFEPtr = pipeline_mng->getDomainRhsFE();
719 
720  MoFEMFunctionReturn(0);
721 }
722 //! [Push operators to pipeline]
723 
724 /**
725  * @brief Monitor solution
726  *
727  * This functions is called by TS solver at the end of each step. It is used
728  * to output results to the hard drive.
729  */
730 struct Monitor : public FEMethod {
731  Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
732  : dM(dm), postProc(post_proc){};
733  MoFEMErrorCode postProcess() {
734  MoFEMFunctionBegin;
735  constexpr int save_every_nth_step = 50;
736  if (ts_step % save_every_nth_step == 0) {
737  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
738  CHKERR postProc->writeFile(
739  "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
740  MOFEM_LOG("SW", Sev::verbose)
741  << "writing vector in binary to vector.dat ...";
742  PetscViewer viewer;
743  PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_WRITE,
744  &viewer);
745  VecView(ts_u, viewer);
746  PetscViewerDestroy(&viewer);
747  }
748  MoFEMFunctionReturn(0);
749  }
750 
751 private:
752  SmartPetscObj<DM> dM;
753  boost::shared_ptr<PostProcEle> postProc;
754 };
755 
756 //! [Solve]
757 MoFEMErrorCode Example::solveSystem() {
758  MoFEMFunctionBegin;
759  auto *simple = mField.getInterface<Simple>();
760  auto *pipeline_mng = mField.getInterface<PipelineManager>();
761 
762  auto dm = simple->getDM();
763 
764  auto set_initial_step = [&](auto ts) {
765  MoFEMFunctionBegin;
766  int step = 0;
767  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-step", &step,
768  PETSC_NULL);
769  CHKERR TSSetStepNumber(ts, step);
770  MoFEMFunctionReturn(0);
771  };
772 
773  auto set_fieldsplit_preconditioner_ksp = [&](auto ksp) {
774  MoFEMFunctionBeginHot;
775  PC pc;
776  CHKERR KSPGetPC(ksp, &pc);
777  PetscBool is_pcfs = PETSC_FALSE;
778  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
779  if (is_pcfs == PETSC_TRUE) {
780  auto bc_mng = mField.getInterface<BcManager>();
781  auto name_prb = simple->getProblemName();
782  SmartPetscObj<IS> is_u;
783  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
784  name_prb, ROW, "U", 0, 3, is_u);
785  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
786  }
787  MoFEMFunctionReturnHot(0);
788  };
789 
790  // Setup postprocessing
791  auto get_fe_post_proc = [&]() {
792  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
793 
794  auto det_ptr = boost::make_shared<VectorDouble>();
795  auto jac_ptr = boost::make_shared<MatrixDouble>();
796  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
797 
798  post_proc_fe->getOpPtrVector().push_back(
799  new OpGetHONormalsOnFace("HO_POSITIONS"));
800  post_proc_fe->getOpPtrVector().push_back(
801  new OpCalculateHOCoords("HO_POSITIONS"));
802  post_proc_fe->getOpPtrVector().push_back(
803  new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
804  post_proc_fe->getOpPtrVector().push_back(
805  new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
806  post_proc_fe->getOpPtrVector().push_back(
807  new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
808 
809  auto u_ptr = boost::make_shared<MatrixDouble>();
810  auto h_ptr = boost::make_shared<VectorDouble>();
811  auto pos_ptr = boost::make_shared<MatrixDouble>();
812 
813  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
814  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
815 
816  post_proc_fe->getOpPtrVector().push_back(
817  new OpCalculateVectorFieldValues<3>("U", u_ptr));
818  post_proc_fe->getOpPtrVector().push_back(
819  new OpCalculateScalarFieldValues("H", h_ptr));
820  post_proc_fe->getOpPtrVector().push_back(
821  new OpCalculateVectorFieldValues<3>("HO_POSITIONS", pos_ptr));
822 
823  post_proc_fe->getOpPtrVector().push_back(
824  new OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
825  post_proc_fe->getOpPtrVector().push_back(
826  new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
827 
828  using OpPPMap = OpPostProcMapInMoab<3, 3>;
829 
830  post_proc_fe->getOpPtrVector().push_back(
831 
832  new OpPPMap(
833  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
834 
835  {{"H", h_ptr}},
836 
837  {{"U", u_ptr}, {"HO_POSITIONS", pos_ptr}, {"GRAD_H", grad_h_ptr}},
838 
839  {{"GRAD_U", grad_u_ptr}}, {}
840 
841  )
842 
843  );
844 
845  return post_proc_fe;
846  };
847 
848  auto set_fieldsplit_preconditioner_ts = [&](auto solver) {
849  MoFEMFunctionBeginHot;
850  SNES snes;
851  CHKERR TSGetSNES(solver, &snes);
852  KSP ksp;
853  CHKERR SNESGetKSP(snes, &ksp);
854  CHKERR set_fieldsplit_preconditioner_ksp(ksp);
855  MoFEMFunctionReturnHot(0);
856  };
857 
858  MoFEM::SmartPetscObj<TS> ts;
859  ts = pipeline_mng->createTSIM();
860 
861  boost::shared_ptr<FEMethod> null_fe;
862  auto monitor_ptr = boost::make_shared<Monitor>(dm, get_fe_post_proc());
863  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
864  null_fe, monitor_ptr);
865 
866  // Add monitor to time solver
867  double ftime = 1;
868  // CHKERR TSSetMaxTime(ts, ftime);
869  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
870 
871  auto T = createDMVector(simple->getDM());
872  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
873  SCATTER_FORWARD);
874  CHKERR TSSetSolution(ts, T);
875  CHKERR TSSetFromOptions(ts);
876  CHKERR set_fieldsplit_preconditioner_ts(ts);
877  CHKERR TSSetUp(ts);
878  CHKERR set_initial_step(ts);
879  CHKERR TSSolve(ts, NULL);
880  CHKERR TSGetTime(ts, &ftime);
881 
882  MoFEMFunctionReturn(0);
883 }
884 //! [Solve]
885 
886 //! [Postprocess results]
887 MoFEMErrorCode Example::outputResults() {
888  MoFEMFunctionBegin;
889  MoFEMFunctionReturn(0);
890 }
891 //! [Postprocess results]
892 
893 //! [Check]
894 MoFEMErrorCode Example::checkResults() {
895  MoFEMFunctionBegin;
896  MoFEMFunctionReturn(0);
897 }
898 //! [Check]
899 
900 static char help[] = "...\n\n";
901 
902 int main(int argc, char *argv[]) {
903 
904  // Initialisation of MoFEM/PETSc and MOAB data structures
905  const char param_file[] = "param_file.petsc";
906  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
907 
908  // Add logging channel for example
909  auto core_log = logging::core::get();
910  core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "SW"));
911  LogManager::setLog("SW");
912  MOFEM_LOG_TAG("SW", "example");
913 
914  try {
915 
916  //! [Register MoFEM discrete manager in PETSc]
917  DMType dm_name = "DMMOFEM";
918  CHKERR DMRegister_MoFEM(dm_name);
919  //! [Register MoFEM discrete manager in PETSc
920 
921  //! [Create MoAB]
922  moab::Core mb_instance; ///< mesh database
923  moab::Interface &moab = mb_instance; ///< mesh database interface
924  //! [Create MoAB]
925 
926  //! [Create MoFEM]
927  MoFEM::Core core(moab); ///< finite element database
928  MoFEM::Interface &m_field = core; ///< finite element database insterface
929  //! [Create MoFEM]
930 
931  //! [Example]
932  Example ex(m_field);
933  CHKERR ex.runProblem();
934  //! [Example]
935  }
936  CATCH_ERRORS;
937 
938  CHKERR MoFEM::Core::Finalize();
939 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
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
g
constexpr double g
Definition: shallow_wave.cpp:63
Monitor::Monitor
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
Definition: shallow_wave.cpp:731
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
FTensor::Tensor1< double, 3 >
Example::Example
Example(MoFEM::Interface &m_field)
Definition: shallow_wave.cpp:339
MoFEM::OpCalculateHOJacForFaceEmbeddedIn3DSpace
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
Definition: HODataOperators.hpp:265
OpURhs
Definition: shallow_wave.cpp:82
OpBaseDivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< 3 > OpBaseDivU
Definition: shallow_wave.cpp:36
MoFEM::OpCalculateDivergenceVectorFieldValues
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:486
EntData
EntitiesFieldData::EntData EntData
Definition: shallow_wave.cpp:26
OpURhs::hGradPtr
boost::shared_ptr< MatrixDouble > hGradPtr
Definition: shallow_wave.cpp:165
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
init_h
auto init_h
Initialisation function.
Definition: free_surface.cpp:295
help
static char help[]
[Check]
Definition: shallow_wave.cpp:900
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:30
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
OpConvectiveTermRhs
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermRhs
Definition: operators_tests.cpp:52
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
OpConvectiveTermLhsDu
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDu
Definition: operators_tests.cpp:56
phi_1
constexpr double phi_1
Definition: shallow_wave.cpp:70
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
i
FTensor::Index< 'i', 3 > i
Definition: shallow_wave.cpp:77
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
OpSourceU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpSourceU
Definition: shallow_wave.cpp:51
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
OpConvectiveH_dGradH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::TriLinearForm< GAUSS >::OpConvectiveTermLhsDy< 1, 1, 3 > OpConvectiveH_dGradH
Definition: shallow_wave.cpp:60
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
Monitor::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: shallow_wave.cpp:733
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::OpCalculateHOCoords
Calculate HO coordinates at gauss points.
Definition: HODataOperators.hpp:41
Example::domianLhsFEPtr
boost::shared_ptr< FEMethod > domianLhsFEPtr
Definition: shallow_wave.cpp:355
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double, 3, 3 >
order
constexpr int order
Definition: dg_projection.cpp:18
OpMixScalarTimesDiv
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
Definition: free_surface.cpp:125
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
h0
constexpr double h0
Definition: shallow_wave.cpp:65
OpULhs_dH
Definition: shallow_wave.cpp:267
OpURhs::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
Class dedicated to integrate operator.
Definition: shallow_wave.cpp:92
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
OpBaseGradH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixVectorTimesGrad< 1, 3, 3 > OpBaseGradH
Definition: shallow_wave.cpp:38
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
OpURhs::uPtr
boost::shared_ptr< MatrixDouble > uPtr
Definition: shallow_wave.cpp:162
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: dynamic_first_order_con_law.cpp:463
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
u_max
constexpr double u_max
Definition: shallow_wave.cpp:68
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
main
int main(int argc, char *argv[])
Definition: shallow_wave.cpp:902
OpURhs::OpURhs
OpURhs(const std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > u_dot_ptr, boost::shared_ptr< MatrixDouble > grad_u_ptr, boost::shared_ptr< MatrixDouble > grad_h_ptr)
Definition: shallow_wave.cpp:84
ReactionDiffusionEquation::u0
const double u0
inital vale on blocksets
Definition: reaction_diffusion.cpp:24
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
phi_2
constexpr double phi_2
Definition: shallow_wave.cpp:71
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
a
constexpr double a
Definition: approx_sphere.cpp:30
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
OpURhs::uDotPtr
boost::shared_ptr< MatrixDouble > uDotPtr
Definition: shallow_wave.cpp:163
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Example
[Example]
Definition: plastic.cpp:177
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
OpURhs::uGradPtr
boost::shared_ptr< MatrixDouble > uGradPtr
Definition: shallow_wave.cpp:164
OpConvectiveTermLhsDy
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDy
Definition: operators_tests.cpp:60
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: dynamic_first_order_con_law.cpp:893
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::OpSetHOWeightsOnFace
Modify integration weights on face to take in account higher-order geometry.
Definition: HODataOperators.hpp:122
j
FTensor::Index< 'j', 3 > j
Definition: shallow_wave.cpp:78
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::OpGetHONormalsOnFace
Calculate normals at Gauss points of triangle element.
Definition: HODataOperators.hpp:281
OpSourceH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpSourceH
Definition: shallow_wave.cpp:53
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
beta_montain
constexpr double beta_montain
Definition: shallow_wave.cpp:73
penalty
constexpr double penalty
Definition: shallow_wave.cpp:75
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
OpBaseTimesDivU
OpBaseTimesDotH OpBaseTimesDivU
Definition: shallow_wave.cpp:48
BiLinearForm
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
FTensor::Tensor1::normalize
Tensor1< T, Tensor_Dim > normalize()
Definition: Tensor1_value.hpp:77
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
ElementsAndOps
Definition: child_and_parent.cpp:18
DomainEleOp
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::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:225
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: dynamic_first_order_con_law.cpp:536
Monitor
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:783
approx_order
int approx_order
Definition: test_broken_space.cpp:50
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
OpConvectiveH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpConvectiveTermRhs< 1, 1, 3 > OpConvectiveH
Definition: shallow_wave.cpp:56
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
OpConvectiveH_dU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::TriLinearForm< GAUSS >::OpConvectiveTermLhsDu< 1, 1, 3 > OpConvectiveH_dU
Definition: shallow_wave.cpp:58
OpULhs_dU::uGradPtr
boost::shared_ptr< MatrixDouble > uGradPtr
Definition: shallow_wave.cpp:264
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3249
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:213
init_u
constexpr double init_u
Definition: heat_equation.cpp:58
Example::domianRhsFEPtr
boost::shared_ptr< FEMethod > domianRhsFEPtr
Definition: shallow_wave.cpp:356
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
OpMassHH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpMassHH
Definition: shallow_wave.cpp:44
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
h_hat
constexpr double h_hat
Definition: shallow_wave.cpp:67
omega
constexpr double omega
Definition: shallow_wave.cpp:62
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
alpha_montain
constexpr double alpha_montain
Definition: shallow_wave.cpp:72
Example::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:396
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
OpULhs_dU::OpULhs_dU
OpULhs_dU(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > grad_u_ptr)
Definition: shallow_wave.cpp:170
MoFEM::SmartPetscObj< IS >
l
FTensor::Index< 'l', 3 > l
Definition: shallow_wave.cpp:79
FE_DIM
constexpr int FE_DIM
Definition: shallow_wave.cpp:24
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
MoFEM::OpSetInvJacH1ForFaceEmbeddedIn3DSpace
Definition: UserDataOperators.hpp:2932
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
OpULhs_dU::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
Definition: shallow_wave.cpp:179
OpMassUU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMassUU
Definition: shallow_wave.cpp:42
OpULhs_dU::uPtr
boost::shared_ptr< MatrixDouble > uPtr
Definition: shallow_wave.cpp:263
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
OpBaseDivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM > OpBaseDivU
[Essential boundary conditions]
Definition: seepage.cpp:79
phi_0
constexpr double phi_0
Definition: shallow_wave.cpp:69
OpULhs_dU
Definition: shallow_wave.cpp:168
F
@ F
Definition: free_surface.cpp:394
OpULhs_dH::OpULhs_dH
OpULhs_dH(const std::string field_name_row, const std::string field_name_col)
Definition: shallow_wave.cpp:269
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
OpULhs_dH::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
Definition: shallow_wave.cpp:275
mu
constexpr double mu
Definition: shallow_wave.cpp:64
OpBaseTimesDotH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseTimesDotH
Definition: shallow_wave.cpp:47
Generated by Doxygen 1.8.17 and hosted at University of Glasgow