v0.12.1
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 /* This file is part of MoFEM.
12  * MoFEM is free software: you can redistribute it and/or modify it under
13  * the terms of the GNU Lesser General Public License as published by the
14  * Free Software Foundation, either version 3 of the License, or (at your
15  * option) any later version.
16  *
17  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
18  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
19  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
20  * License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public
23  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
24 
25 #include <MoFEM.hpp>
26 using namespace MoFEM;
27 
28 #include <boost/math/quadrature/gauss_kronrod.hpp>
29 using namespace boost::math::quadrature;
30 
31 template <int DIM> struct ElementsAndOps {};
32 
33 template <> struct ElementsAndOps<2> {
37 };
38 
39 constexpr int FE_DIM = 2;
40 
45 
48 
49 // Use forms iterators for Grad-Grad term
54 
55 // Use forms for Mass term
60 
64 
66  GAUSS>::OpSource<1, 3>;
68  GAUSS>::OpSource<1, 1>;
69 
76 
77 constexpr double omega = 7.292 * 1e-5;
78 constexpr double g = 9.80616;
79 constexpr double mu = 1e4;
80 constexpr double h0 = 1e4;
81 
82 constexpr double h_hat = 120;
83 constexpr double u_max = 80;
84 constexpr double phi_0 = M_PI / 7;
85 constexpr double phi_1 = M_PI / 2 - phi_0;
86 constexpr double phi_2 = M_PI / 4;
87 constexpr double alpha_montain = 1. / 3.;
88 constexpr double beta_montain = 1. / 15.;
89 
90 constexpr double penalty = 1;
91 
96 
97 struct OpURhs : public AssemblyDomainEleOp {
98 
99  OpURhs(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
100  boost::shared_ptr<MatrixDouble> u_dot_ptr,
101  boost::shared_ptr<MatrixDouble> grad_u_ptr,
102  boost::shared_ptr<MatrixDouble> grad_h_ptr)
103  : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
104  uPtr(u_ptr), uDotPtr(u_dot_ptr), uGradPtr(grad_u_ptr),
105  hGradPtr(grad_h_ptr) {}
106 
109 
110  const double vol = getMeasure();
111  auto t_w = getFTensor0IntegrationWeight();
112  auto t_row_base = row_data.getFTensor0N();
113  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
114  auto t_dot_u = getFTensor1FromMat<3>(*uDotPtr);
115  auto t_u = getFTensor1FromMat<3>(*uPtr);
116  auto t_grad_u = getFTensor2FromMat<3, 3>(*uGradPtr);
117  auto t_grad_h = getFTensor1FromMat<3>(*hGradPtr);
118  auto t_coords = getFTensor1CoordsAtGaussPts();
119  auto t_normal = getFTensor1NormalsAtGaussPts();
120 
121  for (int gg = 0; gg != nbIntegrationPts; gg++) {
122 
123  const double alpha = t_w * vol;
124  auto t_nf = getFTensor1FromArray<3, 3>(locF);
125 
128 
129  const auto a = std::sqrt(t_coords(i) * t_coords(i));
130  const auto sin_fi = t_coords(2) / a;
131  const auto f = 2 * omega * sin_fi;
132 
134  t_r(i) = t_normal(i);
135  t_r.normalize();
136 
137  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
139  t_P(i, j) = t_r(i) * t_r(j);
140  t_Q(i, j) = t_kd(i, j) - t_P(i, j);
141 
143  t_A(i, m) = levi_civita(i, j, m) * t_r(j);
144 
145  t_rhs(m) = t_Q(m, i) * (t_dot_u(i) + t_grad_u(i, j) * t_u(j) +
146  f * t_A(i, j) * t_u(j) + g * t_grad_h(i));
147  t_rhs_grad(m, j) = t_Q(m, i) * (mu * t_grad_u(i, j));
148 
149  t_rhs(m) += t_P(m, j) * t_u(j);
150 
151  int rr = 0;
152  for (; rr != nbRows / 3; ++rr) {
153  t_nf(i) += alpha * t_row_base * t_rhs(i);
154  t_nf(i) += alpha * t_row_diff_base(j) * t_rhs_grad(i, j);
155  ++t_row_base;
156  ++t_row_diff_base;
157  ++t_nf;
158  }
159  for (; rr < nbRowBaseFunctions; ++rr) {
160  ++t_row_base;
161  ++t_row_diff_base;
162  }
163 
164  ++t_w;
165  ++t_u;
166  ++t_dot_u;
167  ++t_grad_u;
168  ++t_grad_h;
169  ++t_coords;
170  ++t_normal;
171  }
172 
174  }
175 
176 private:
177  boost::shared_ptr<MatrixDouble> uPtr;
178  boost::shared_ptr<MatrixDouble> uDotPtr;
179  boost::shared_ptr<MatrixDouble> uGradPtr;
180  boost::shared_ptr<MatrixDouble> hGradPtr;
181 };
182 
184 
185  OpULhs_dU(const std::string field_name_row, const std::string field_name_col,
186  boost::shared_ptr<MatrixDouble> u_ptr,
187  boost::shared_ptr<MatrixDouble> grad_u_ptr)
188  : AssemblyDomainEleOp(field_name_row, field_name_col,
189  AssemblyDomainEleOp::OPROWCOL),
190  uPtr(u_ptr), uGradPtr(grad_u_ptr) {
191  this->sYmm = false;
192  }
193 
197 
198  const double vol = getMeasure();
199  auto t_w = getFTensor0IntegrationWeight();
200  auto t_row_base = row_data.getFTensor0N();
201  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
202  auto t_coords = getFTensor1CoordsAtGaussPts();
203  auto t_normal = getFTensor1NormalsAtGaussPts();
204 
205  auto t_u = getFTensor1FromMat<3>(*uPtr);
206  auto t_grad_u = getFTensor2FromMat<3, 3>(*uGradPtr);
207 
208  auto get_t_mat = [&](const int rr) {
210  &locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
211 
212  &locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
213 
214  &locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
215  };
216 
217  const auto ts_a = getFEMethod()->ts_a;
218 
219  for (int gg = 0; gg != nbIntegrationPts; gg++) {
220 
221  const auto a = std::sqrt(t_coords(i) * t_coords(i));
222  const auto sin_fi = t_coords(2) / a;
223  const auto f = 2 * omega * sin_fi;
224 
226  t_r(i) = t_normal(i);
227  t_r.normalize();
228 
229  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
231  t_P(i, j) = t_r(i) * t_r(j);
232  t_Q(i, j) = t_kd(i, j) - t_P(i, j);
233 
235  t_A(i, m) = levi_civita(i, j, m) * t_r(j);
236 
238  t_rhs_du(m, j) =
239  t_Q(m, i) * (ts_a * t_kd(i, j) + t_grad_u(i, j) + f * t_A(i, j)) +
240  t_P(m, j);
241 
242  const double alpha = t_w * vol;
243 
244  int rr = 0;
245  for (; rr != nbRows / 3; rr++) {
246 
247  auto t_col_base = col_data.getFTensor0N(gg, 0);
248  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
249  auto t_mat = get_t_mat(3 * rr);
250 
251  for (int cc = 0; cc != nbCols / 3; cc++) {
252  t_mat(i, j) += (alpha * t_row_base * t_col_base) * t_rhs_du(i, j);
253  t_mat(i, j) += (alpha * mu) * t_Q(i, j) *
254  (t_row_diff_base(m) * t_col_diff_base(m));
255  ++t_col_diff_base;
256  ++t_col_base;
257  ++t_mat;
258  }
259  ++t_row_base;
260  ++t_row_diff_base;
261  }
262  for (; rr < nbRowBaseFunctions; ++rr) {
263  ++t_row_base;
264  ++t_row_diff_base;
265  }
266 
267  ++t_w;
268  ++t_coords;
269  ++t_normal;
270  ++t_u;
271  ++t_grad_u;
272  }
273 
275  }
276 
277 private:
278  boost::shared_ptr<MatrixDouble> uPtr;
279  boost::shared_ptr<MatrixDouble> uGradPtr;
280 };
281 
283 
284  OpULhs_dH(const std::string field_name_row, const std::string field_name_col)
285  : AssemblyDomainEleOp(field_name_row, field_name_col,
286  AssemblyDomainEleOp::OPROWCOL) {
287  this->sYmm = false;
288  }
289 
293 
294  // get element volume
295  const double vol = getMeasure();
296  // get integration weights
297  auto t_w = getFTensor0IntegrationWeight();
298  // get base function gradient on rows
299  auto t_row_base = row_data.getFTensor0N();
300  // normal
301  auto t_normal = getFTensor1NormalsAtGaussPts();
302 
303  auto get_t_vec = [&](const int rr) {
305  &locMat(rr + 0, 0),
306 
307  &locMat(rr + 1, 0),
308 
309  &locMat(rr + 2, 0)};
310  };
311 
312  // loop over integration points
313  for (int gg = 0; gg != nbIntegrationPts; gg++) {
314 
316  t_r(i) = t_normal(i);
317  t_r.normalize();
318 
319  constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
321  t_P(i, j) = t_r(i) * t_r(j);
322  t_Q(i, j) = t_kd(i, j) - t_P(i, j);
323 
324  const double alpha = t_w * vol;
325 
326  int rr = 0;
327  for (; rr != nbRows / 3; rr++) {
328  auto t_vec = get_t_vec(3 * rr);
329  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
330  const double a = alpha * g * t_row_base;
331 
332  for (int cc = 0; cc != nbCols; cc++) {
333  t_vec(i) += a * (t_Q(i, m) * t_col_diff_base(m));
334  ++t_vec;
335  ++t_col_diff_base;
336  }
337 
338  ++t_row_base;
339  }
340 
341  for (; rr < nbRowBaseFunctions; ++rr)
342  ++t_row_base;
343 
344  ++t_w;
345  ++t_normal;
346  }
347 
349  }
350 };
351 
352 struct Example {
353 
354  Example(MoFEM::Interface &m_field) : mField(m_field) {}
355 
357 
358 private:
359  MoFEM::Interface &mField;
360 
369 
370  boost::shared_ptr<FEMethod> domianLhsFEPtr;
371  boost::shared_ptr<FEMethod> domianRhsFEPtr;
372 };
373 
374 //! [Create common data]
377 
379 }
380 //! [Create common data]
381 
382 //! [Run problem]
385  CHKERR readMesh();
386  CHKERR setupProblem();
387  CHKERR createCommonData();
388  CHKERR boundaryCondition();
389  CHKERR assembleSystem();
390  CHKERR solveSystem();
391  CHKERR outputResults();
392  CHKERR checkResults();
394 }
395 //! [Run problem]
396 
397 //! [Read mesh]
400  auto simple = mField.getInterface<Simple>();
401  CHKERR simple->getOptions();
402  CHKERR simple->loadFile();
404 }
405 //! [Read mesh]
406 
407 //! [Set up problem]
410  Simple *simple = mField.getInterface<Simple>();
411  // Add field
412  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 3);
413  CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
414  CHKERR simple->addDataField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
415 
416  int order = 2;
417  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
418  CHKERR simple->setFieldOrder("U", order);
419  CHKERR simple->setFieldOrder("H", order);
420  CHKERR simple->setFieldOrder("HO_POSITIONS", 3);
421 
422  CHKERR simple->setUp();
424 }
425 //! [Set up problem]
426 
427 //! [Boundary condition]
430 
431  PetscBool is_restart = PETSC_FALSE;
432  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_restart", &is_restart,
433  PETSC_NULL);
434 
435  auto restart_vector = [&]() {
437  auto simple = mField.getInterface<Simple>();
438  auto dm = simple->getDM();
439  MOFEM_LOG("SW", Sev::inform)
440  << "reading vector in binary from vector.dat ...";
441  PetscViewer viewer;
442  PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_READ,
443  &viewer);
444  auto T = smartCreateDMVector(simple->getDM());
445  VecLoad(T, viewer);
446  CHKERR DMoFEMMeshToLocalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
448  };
449 
450  if (is_restart) {
451 
452  CHKERR restart_vector();
453 
454  } else {
455 
456  const double e_n = exp(-4 / pow(phi_1 - phi_0, 2));
457  const double u0 = u_max / e_n;
458 
459  FTensor::Tensor1<double, 3> t_k{0., 0., 1.};
461  t_A(i, m) = levi_civita(i, j, m) * t_k(j);
462 
463  auto get_phi = [&](const double x, const double y, const double z) {
464  FTensor::Tensor1<double, 3> t_r{x, y, 0.};
465  const double r = std::sqrt(t_r(i) * t_r(i));
466  return atan2(z, r);
467  };
468 
469  auto init_u_phi = [&](const double phi) {
470  if (phi > phi_0 && phi < phi_1) {
471  return u0 * exp(1. / ((phi - phi_0) * (phi - phi_1)));
472  } else {
473  return 0.;
474  }
475  };
476 
477  auto init_u = [&](const double x, const double y, const double z) {
478  FTensor::Tensor1<double, 3> t_u{0., 0., 0.};
479  const double u_phi = init_u_phi(get_phi(x, y, z));
480  if (u_phi > 0) {
481  FTensor::Tensor1<double, 3> t_r{x, y, 0.};
482  t_r.normalize();
483  t_u(i) = ((t_A(i, j) * t_r(j)) * u_phi);
484  }
485  return t_u;
486  };
487 
488  auto init_h = [&](const double x, const double y, const double z) {
489  const double a = std::sqrt(x * x + y * y + z * z);
490 
491  auto integral = [&](const double fi) {
492  const double u_phi = init_u_phi(fi);
493  const auto f = 2 * omega * sin(fi);
494  return a * u_phi * (f + (tan(fi) / a) * u_phi);
495  };
496 
497  auto montain = [&](const double lambda, const double fi) {
498  if (lambda > -M_PI && lambda < M_PI)
499  return h_hat * cos(fi) * exp(-pow(lambda / alpha_montain, 2)) *
500  exp(-pow((phi_2 - fi) / beta_montain, 2));
501  else
502  return 0.;
503  };
504 
505  const double fi = get_phi(x, y, z);
506  const double lambda = atan2(x, y);
507 
508  double h1 = 0;
509  if (fi > phi_0)
510  h1 = gauss_kronrod<double, 32>::integrate(
511  integral, phi_0, fi, 0, std::numeric_limits<float>::epsilon());
512 
513  return h0 + (montain(lambda, fi) - (h1 / g));
514  };
515 
516  auto set_domain_general = [&](auto &pipeline) {
517  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
518  pipeline.push_back(new OpGetHONormalsOnFace("HO_POSITIONS"));
519  pipeline.push_back(new OpCalculateHOCoords("HO_POSITIONS"));
520  pipeline.push_back(new OpSetHOWeightsOnFace());
521  };
522 
523  auto set_domain_rhs = [&](auto &pipeline) {
524  pipeline.push_back(new OpSourceU("U", init_u));
525  pipeline.push_back(new OpSourceH("H", init_h));
526  };
527 
528  auto set_domain_lhs = [&](auto &pipeline) {
529  pipeline.push_back(
530  new OpMassUU("U", "U", [](double, double, double) { return 1; }));
531  pipeline.push_back(
532  new OpMassHH("H", "H", [](double, double, double) { return 1; }));
533  };
534 
535  auto post_proc = [&]() {
537  auto simple = mField.getInterface<Simple>();
538  auto dm = simple->getDM();
539 
540  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
541  post_proc_fe->generateReferenceElementMesh();
542 
543  auto det_ptr = boost::make_shared<VectorDouble>();
544  auto jac_ptr = boost::make_shared<MatrixDouble>();
545  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
546 
547  post_proc_fe->getOpPtrVector().push_back(
548  new OpGetHONormalsOnFace("HO_POSITIONS"));
549  post_proc_fe->getOpPtrVector().push_back(
550  new OpCalculateHOCoords("HO_POSITIONS"));
551  post_proc_fe->getOpPtrVector().push_back(
553  post_proc_fe->getOpPtrVector().push_back(
554  new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
555  post_proc_fe->addFieldValuesPostProc("U");
556  post_proc_fe->addFieldValuesPostProc("H");
557  post_proc_fe->addFieldValuesPostProc("HO_POSITIONS");
558 
559  CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
560  CHKERR post_proc_fe->writeFile("out_init.h5m");
561 
563  };
564 
565  auto solve_init = [&]() {
567  auto simple = mField.getInterface<Simple>();
568  auto pipeline_mng = mField.getInterface<PipelineManager>();
569 
570  auto solver = pipeline_mng->createKSP();
571  CHKERR KSPSetFromOptions(solver);
572  PC pc;
573  CHKERR KSPGetPC(solver, &pc);
574  PetscBool is_pcfs = PETSC_FALSE;
575  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
576  if (is_pcfs == PETSC_TRUE) {
577  auto bc_mng = mField.getInterface<BcManager>();
578  auto name_prb = simple->getProblemName();
579  SmartPetscObj<IS> is_u;
580  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
581  name_prb, ROW, "U", 0, 3, is_u);
582  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
583  CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_ADDITIVE);
584  }
585 
586  CHKERR KSPSetUp(solver);
587 
588  auto dm = simple->getDM();
589  auto D = smartCreateDMVector(dm);
590  auto F = smartVectorDuplicate(D);
591 
592  CHKERR KSPSolve(solver, F, D);
593  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
594  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
595  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
596 
598  };
599 
600  auto pipeline_mng = mField.getInterface<PipelineManager>();
601 
602  auto integration_rule = [](int, int, int approx_order) {
603  return 2 * approx_order + 4;
604  };
605  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
606  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
607 
608  set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
609  set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
610  set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
611  set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
612 
613  CHKERR solve_init();
614  CHKERR post_proc();
615 
616  // Clear pipelines
617  pipeline_mng->getOpDomainRhsPipeline().clear();
618  pipeline_mng->getOpDomainLhsPipeline().clear();
619  }
620 
622 }
623 //! [Boundary condition]
624 
625 //! [Push operators to pipeline]
628 
629  // Push element from reference configuration to current configuration in 3d
630  // space
631  auto set_domain_general = [&](auto &pipeline) {
632  auto det_ptr = boost::make_shared<VectorDouble>();
633  auto jac_ptr = boost::make_shared<MatrixDouble>();
634  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
635  pipeline.push_back(new OpGetHONormalsOnFace("HO_POSITIONS"));
636  pipeline.push_back(new OpCalculateHOCoords("HO_POSITIONS"));
637  pipeline.push_back(new OpSetHOWeightsOnFace());
638  pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
639  pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
640  pipeline.push_back(new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
641  };
642 
643  auto set_domain_rhs = [&](auto &pipeline) {
644  auto dot_u_ptr = boost::make_shared<MatrixDouble>();
645  auto u_ptr = boost::make_shared<MatrixDouble>();
646  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
647  auto div_u_ptr = boost::make_shared<VectorDouble>();
648  auto dot_h_ptr = boost::make_shared<VectorDouble>();
649  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
650 
651  pipeline.push_back(new OpCalculateVectorFieldValuesDot<3>("U", dot_u_ptr));
652  pipeline.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
653 
654  pipeline.push_back(new OpCalculateVectorFieldValues<3>("U", u_ptr));
655  pipeline.push_back(
656  new OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
657  pipeline.push_back(
658  new OpCalculateDivergenceVectorFieldValues<3>("U", div_u_ptr));
659  pipeline.push_back(new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
660 
661  pipeline.push_back(new OpBaseTimesDotH(
662  "H", dot_h_ptr, [](double, double, double) { return 1.; }));
663  pipeline.push_back(new OpBaseTimesDivU(
664  "H", div_u_ptr, [](double, double, double) { return h0; }));
665  pipeline.push_back(new OpConvectiveH("H", u_ptr, grad_h_ptr));
666  pipeline.push_back(
667  new OpURhs("U", u_ptr, dot_u_ptr, grad_u_ptr, grad_h_ptr));
668  };
669 
670  auto set_domain_lhs = [&](auto &pipeline) {
671  auto u_ptr = boost::make_shared<MatrixDouble>();
672  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
673  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
674 
675  pipeline.push_back(new OpCalculateVectorFieldValues<3>("U", u_ptr));
676  pipeline.push_back(
677  new OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
678  pipeline.push_back(new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
679 
680  pipeline.push_back(new OpMassHH("H", "H", [&](double, double, double) {
681  return domianLhsFEPtr->ts_a;
682  }));
683  pipeline.push_back(new OpBaseDivU(
684  "H", "U", [](const double, const double, const double) { return h0; },
685  false, false));
686  pipeline.push_back(
687  new OpConvectiveH_dU("H", "U", grad_h_ptr, []() { return 1; }));
688  pipeline.push_back(
689  new OpConvectiveH_dGradH("H", "H", u_ptr, []() { return 1; }));
690  pipeline.push_back(new OpULhs_dU("U", "U", u_ptr, grad_u_ptr));
691  pipeline.push_back(new OpULhs_dH("U", "H"));
692  };
693 
694  auto *pipeline_mng = mField.getInterface<PipelineManager>();
695 
696  auto integration_rule = [](int, int, int approx_order) {
697  return 2 * approx_order + 4;
698  };
699  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
700  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
701 
702  set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
703  set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
704 
705  set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
706  set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
707 
708  domianLhsFEPtr = pipeline_mng->getDomainLhsFE();
709  domianRhsFEPtr = pipeline_mng->getDomainRhsFE();
710 
712 }
713 //! [Push operators to pipeline]
714 
715 /**
716  * @brief Monitor solution
717  *
718  * This functions is called by TS solver at the end of each step. It is used
719  * to output results to the hard drive.
720  */
721 struct Monitor : public FEMethod {
722  Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
723  : dM(dm), postProc(post_proc){};
726  constexpr int save_every_nth_step = 50;
727  if (ts_step % save_every_nth_step == 0) {
728  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
729  CHKERR postProc->writeFile(
730  "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
731  MOFEM_LOG("SW", Sev::verbose)
732  << "writing vector in binary to vector.dat ...";
733  PetscViewer viewer;
734  PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_WRITE,
735  &viewer);
736  VecView(ts_u, viewer);
737  PetscViewerDestroy(&viewer);
738  }
740  }
741 
742 private:
744  boost::shared_ptr<PostProcEle> postProc;
745 };
746 
747 //! [Solve]
750  auto *simple = mField.getInterface<Simple>();
751  auto *pipeline_mng = mField.getInterface<PipelineManager>();
752 
753  auto dm = simple->getDM();
754 
755  auto set_initial_step = [&](auto ts) {
757  int step = 0;
758  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-step", &step,
759  PETSC_NULL);
760  CHKERR TSSetStepNumber(ts, step);
762  };
763 
764  auto set_fieldsplit_preconditioner_ksp = [&](auto ksp) {
766  PC pc;
767  CHKERR KSPGetPC(ksp, &pc);
768  PetscBool is_pcfs = PETSC_FALSE;
769  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
770  if (is_pcfs == PETSC_TRUE) {
771  auto bc_mng = mField.getInterface<BcManager>();
772  auto name_prb = simple->getProblemName();
773  SmartPetscObj<IS> is_u;
774  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
775  name_prb, ROW, "U", 0, 3, is_u);
776  CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
777  }
779  };
780 
781  // Setup postprocessing
782  auto get_fe_post_proc = [&]() {
783  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
784  post_proc_fe->generateReferenceElementMesh();
785 
786  auto det_ptr = boost::make_shared<VectorDouble>();
787  auto jac_ptr = boost::make_shared<MatrixDouble>();
788  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
789 
790  post_proc_fe->getOpPtrVector().push_back(
791  new OpGetHONormalsOnFace("HO_POSITIONS"));
792  post_proc_fe->getOpPtrVector().push_back(
793  new OpCalculateHOCoords("HO_POSITIONS"));
794  post_proc_fe->getOpPtrVector().push_back(
796  post_proc_fe->getOpPtrVector().push_back(
797  new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
798  post_proc_fe->getOpPtrVector().push_back(
799  new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
800  post_proc_fe->addFieldValuesPostProc("U");
801  post_proc_fe->addFieldValuesPostProc("H");
802  post_proc_fe->addFieldValuesGradientPostProc("U");
803  post_proc_fe->addFieldValuesGradientPostProc("H");
804  post_proc_fe->addFieldValuesPostProc("HO_POSITIONS");
805  return post_proc_fe;
806  };
807 
808  auto set_fieldsplit_preconditioner_ts = [&](auto solver) {
810  SNES snes;
811  CHKERR TSGetSNES(solver, &snes);
812  KSP ksp;
813  CHKERR SNESGetKSP(snes, &ksp);
814  CHKERR set_fieldsplit_preconditioner_ksp(ksp);
816  };
817 
819  ts = pipeline_mng->createTSIM();
820 
821  boost::shared_ptr<FEMethod> null_fe;
822  auto monitor_ptr = boost::make_shared<Monitor>(dm, get_fe_post_proc());
823  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
824  null_fe, monitor_ptr);
825 
826  // Add monitor to time solver
827  double ftime = 1;
828  // CHKERR TSSetMaxTime(ts, ftime);
829  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
830 
831  auto T = smartCreateDMVector(simple->getDM());
832  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
833  SCATTER_FORWARD);
834  CHKERR TSSetSolution(ts, T);
835  CHKERR TSSetFromOptions(ts);
836  CHKERR set_fieldsplit_preconditioner_ts(ts);
837  CHKERR TSSetUp(ts);
838  CHKERR set_initial_step(ts);
839  CHKERR TSSolve(ts, NULL);
840  CHKERR TSGetTime(ts, &ftime);
841 
843 }
844 //! [Solve]
845 
846 //! [Postprocess results]
850 }
851 //! [Postprocess results]
852 
853 //! [Check]
857 }
858 //! [Check]
859 
860 static char help[] = "...\n\n";
861 
862 int main(int argc, char *argv[]) {
863 
864  // Initialisation of MoFEM/PETSc and MOAB data structures
865  const char param_file[] = "param_file.petsc";
866  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
867 
868  // Add logging channel for example
869  auto core_log = logging::core::get();
870  core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "SW"));
871  LogManager::setLog("SW");
872  MOFEM_LOG_TAG("SW", "example");
873 
874  try {
875 
876  //! [Register MoFEM discrete manager in PETSc]
877  DMType dm_name = "DMMOFEM";
878  CHKERR DMRegister_MoFEM(dm_name);
879  //! [Register MoFEM discrete manager in PETSc
880 
881  //! [Create MoAB]
882  moab::Core mb_instance; ///< mesh database
883  moab::Interface &moab = mb_instance; ///< mesh database interface
884  //! [Create MoAB]
885 
886  //! [Create MoFEM]
887  MoFEM::Core core(moab); ///< finite element database
888  MoFEM::Interface &m_field = core; ///< finite element database insterface
889  //! [Create MoFEM]
890 
891  //! [Example]
892  Example ex(m_field);
893  CHKERR ex.runProblem();
894  //! [Example]
895  }
896  CATCH_ERRORS;
897 
899 }
std::string param_file
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
constexpr double a
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double alpha
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
const double init_u
auto integration_rule
constexpr double lambda
auto init_h
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:478
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:541
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:309
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:340
const double T
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMMoFEM.cpp:991
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1095
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1949
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
const double r
rate factor
const double u0
inital vale on blocksets
int save_every_nth_step
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
constexpr int FE_DIM
constexpr double mu
int main(int argc, char *argv[])
constexpr double phi_0
static char help[]
[Check]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpMassHH
constexpr double phi_1
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpSourceU
constexpr double penalty
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBaseTimesDotH
constexpr double beta_montain
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMassUU
FTensor::Index< 'l', 3 > l
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpConvectiveTermLhsDu< 1, 1, 3 > OpConvectiveH_dU
constexpr double h0
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpConvectiveTermRhs< 1, 1, 3 > OpConvectiveH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< 3 > OpBaseDivU
OpBaseTimesDotH OpBaseTimesDivU
constexpr double u_max
FTensor::Index< 'j', 3 > j
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpConvectiveTermLhsDy< 1, 1, 3 > OpConvectiveH_dGradH
constexpr double omega
FTensor::Index< 'i', 3 > i
DataForcesAndSourcesCore::EntData EntData
constexpr double g
FTensor::Index< 'm', 3 > m
constexpr double h_hat
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpSourceH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixVectorTimesGrad< 1, 3, 3 > OpBaseGradH
constexpr double phi_2
constexpr double alpha_montain
PipelineManager::FaceEle DomainEle
[Example]
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
boost::shared_ptr< FEMethod > domianRhsFEPtr
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
boost::shared_ptr< FEMethod > domianLhsFEPtr
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Deprecated interface functions.
structure for User Loop Methods on finite elements
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:33
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Calculate HO coordinates at gauss points.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field valuse for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
MoFEMErrorCode postProcess()
function is run at the end of loop
OpULhs_dH(const std::string field_name_row, const std::string field_name_col)
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > uPtr
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > uGradPtr
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)
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data)
Class dedicated to integrate operator.
boost::shared_ptr< MatrixDouble > uDotPtr
boost::shared_ptr< MatrixDouble > uGradPtr
boost::shared_ptr< MatrixDouble > hGradPtr
boost::shared_ptr< MatrixDouble > uPtr
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)
Postprocess on face.