v0.13.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>
26using namespace MoFEM;
27
28#include <boost/math/quadrature/gauss_kronrod.hpp>
29using namespace boost::math::quadrature;
30
31template <int DIM> struct ElementsAndOps {};
32
33template <> struct ElementsAndOps<2> {
37};
38
39constexpr 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
77constexpr double omega = 7.292 * 1e-5;
78constexpr double g = 9.80616;
79constexpr double mu = 1e4;
80constexpr double h0 = 1e4;
81
82constexpr double h_hat = 120;
83constexpr double u_max = 80;
84constexpr double phi_0 = M_PI / 7;
85constexpr double phi_1 = M_PI / 2 - phi_0;
86constexpr double phi_2 = M_PI / 4;
87constexpr double alpha_montain = 1. / 3.;
88constexpr double beta_montain = 1. / 15.;
89
90constexpr double penalty = 1;
91
96
97struct 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)
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
176private:
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
195 EntitiesFieldData::EntData &col_data) {
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
277private:
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
291 EntitiesFieldData::EntData &col_data) {
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
352struct Example {
353
354 Example(MoFEM::Interface &m_field) : mField(m_field) {}
355
357
358private:
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]
394}
395//! [Run problem]
396
397//! [Read mesh]
401 CHKERR simple->getOptions();
402 CHKERR simple->loadFile();
404}
405//! [Read mesh]
406
407//! [Set up problem]
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 = [&]() {
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 = [&]() {
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 = [&]() {
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();
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(
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 */
721struct 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) {
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
742private:
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();
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(
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
860static char help[] = "...\n\n";
861
862int main(int argc, char *argv[]) {
863
864 // Initialisation of MoFEM/PETSc and MOAB data structures
865 const char param_file[] = "param_file.petsc";
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 }
897
899}
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
constexpr double a
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
@ 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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double alpha
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
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:482
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:545
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
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:1015
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)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
const double r
rate factor
const double u0
inital vale on blocksets
int save_every_nth_step
static double phi
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
constexpr int FE_DIM
constexpr double mu
int main(int argc, char *argv[])
constexpr double phi_0
EntitiesFieldData::EntData EntData
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
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
[Example]
Definition: plastic.cpp:130
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:137
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
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
structure for User Loop Methods on finite elements
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:33
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
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:35
PetscInt ts_step
time step
Vec & ts_u
state vector
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
SmartPetscObj< DM > dM
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
OpULhs_dH(const std::string field_name_row, const std::string field_name_col)
boost::shared_ptr< MatrixDouble > uPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::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(EntitiesFieldData::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.