v0.13.1
free_surface.cpp
Go to the documentation of this file.
1/**
2 * \file free_surface.cpp
3 * \example free_surface.cpp
4 *
5 * Using PipelineManager interface calculate the divergence of base functions,
6 * and integral of flux on the boundary. Since the h-div space is used, volume
7 * integral and boundary integral should give the same result.
8 */
9
10#include <MoFEM.hpp>
11
12using namespace MoFEM;
13
14static char help[] = "...\n\n";
15
16#include <BasicFiniteElements.hpp>
17
18constexpr int BASE_DIM = 1;
19constexpr int SPACE_DIM = 2;
20constexpr int U_FIELD_DIM = SPACE_DIM;
22 CARTESIAN; ///< select coordinate system <CARTESIAN, CYLINDRICAL>;
23
24template <int DIM> struct ElementsAndOps {};
25
26template <> struct ElementsAndOps<2> {
32};
33
40
42
45
50
56
61
64
67
72
74
75// mesh refinement
76constexpr int order = 3; ///< approximation order
77
78// Physical parameters
79constexpr double a0 = 0.98;
80constexpr double rho_m = 0.998;
81constexpr double mu_m = 0.0101;
82constexpr double rho_p = 0.0012;
83constexpr double mu_p = 0.000182;
84constexpr double lambda = 7.4;
85constexpr double W = 0.25;
86
87template <int T> constexpr int powof2() {
88 if constexpr (T == 0)
89 return 1;
90 else
91 return powof2<T - 1>() * 2;
92};
93
94// Model parameters
95constexpr double h = 0.025; // mesh size
96constexpr double eta = h;
97constexpr double eta2 = eta * eta;
98
99// Numerical parameteres
100constexpr double md = 1e-2;
101constexpr double eps = 1e-12;
102constexpr double tol = std::numeric_limits<float>::epsilon();
103
104constexpr double rho_ave = (rho_p + rho_m) / 2;
105constexpr double rho_diff = (rho_p - rho_m) / 2;
106constexpr double mu_ave = (mu_p + mu_m) / 2;
107constexpr double mu_diff = (mu_p - mu_m) / 2;
108
109const double kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
110
111auto integration_rule = [](int, int, int) { return 2 * order; };
112
113auto cylindrical = [](const double r) {
114 // When we move to C++17 add if constexpr()
115 if constexpr (coord_type == CYLINDRICAL)
116 return 2 * M_PI * r;
117 else
118 return 1.;
119};
120
121auto my_max = [](const double x) { return (x - 1 + std::abs(x + 1)) / 2; };
122auto my_min = [](const double x) { return (x + 1 - std::abs(x - 1)) / 2; };
123auto cut_off = [](const double h) { return my_max(my_min(h)); };
124auto d_cut_off = [](const double h) {
125 if (h >= -1 && h < 1)
126 return 1.;
127 else
128 return 0.;
129};
130
131auto phase_function = [](const double h, const double diff, const double ave) {
132 return diff * cut_off(h) + ave;
133};
134
135auto d_phase_function_h = [](const double h, const double diff) {
136 return diff * d_cut_off(h);
137};
138
139auto get_f = [](const double h) { return 4 * W * h * (h * h - 1); };
140auto get_f_dh = [](const double h) { return 4 * W * (3 * h * h - 1); };
141
142auto get_M0 = [](auto h) { return md; };
143auto get_M0_dh = [](auto h) { return 0; };
144
145auto get_M2 = [](auto h_tmp) {
146 const double h = cut_off(h_tmp);
147 return md * (1 - h * h);
148};
149
150auto get_M2_dh = [](auto h_tmp) {
151 const double h = cut_off(h_tmp);
152 return -md * 2 * h * d_cut_off(h_tmp);
153};
154
155auto get_M3 = [](auto h_tmp) {
156 const double h = cut_off(h_tmp);
157 const double h2 = h * h;
158 const double h3 = h2 * h;
159 if (h >= 0)
160 return md * (2 * h3 - 3 * h2 + 1);
161 else
162 return md * (-2 * h3 - 3 * h2 + 1);
163};
164
165auto get_M3_dh = [](auto h_tmp) {
166 const double h = cut_off(h_tmp);
167 if (h >= 0)
168 return md * (6 * h * (h - 1)) * d_cut_off(h_tmp);
169 else
170 return md * (-6 * h * (h + 1)) * d_cut_off(h_tmp);
171};
172
173auto get_M = [](auto h) { return get_M0(h); };
174auto get_M_dh = [](auto h) { return get_M0_dh(h); };
175
176auto get_D = [](const double A) {
178 t_D(i, j, k, l) = A * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
179 return t_D;
180};
181
182auto kernel_oscillation = [](double r, double y, double) {
183 constexpr int n = 3;
184 constexpr double R = 0.0125;
185 constexpr double A = R * 0.2;
186 const double theta = atan2(r, y);
187 const double w = R + A * cos(n * theta);
188 const double d = std::sqrt(r * r + y * y);
189 return tanh((w - d) / (eta * std::sqrt(2)));
190};
191
192auto kernel_eye = [](double r, double y, double) {
193 constexpr double y0 = 0.5;
194 constexpr double R = 0.4;
195 y -= y0;
196 const double d = std::sqrt(r * r + y * y);
197 return tanh((R - d) / (eta * std::sqrt(2)));
198};
199
200auto init_h = [](double r, double y, double theta) {
201 return kernel_eye(r, y, theta);
202};
203
204auto bit = [](auto b) { return BitRefLevel().set(b); };
205auto marker = [](auto b) { return BitRefLevel().set(BITREFLEVEL_SIZE - b); };
206
207auto save_range = [](moab::Interface &moab, const std::string name,
208 const Range r) {
210 EntityHandle out_meshset;
211 CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
212 CHKERR moab.add_entities(out_meshset, r);
213 CHKERR moab.write_file(name.c_str(), "VTK", "", &out_meshset, 1);
214 CHKERR moab.delete_entities(&out_meshset, 1);
216};
217
218auto get_dofs_ents = [](auto dm) {
219 auto prb_ptr = getProblemPtr(dm);
220 std::vector<EntityHandle> ents_vec;
221 ents_vec.reserve(prb_ptr->numeredRowDofsPtr->size());
222 for (auto dof : *prb_ptr->numeredRowDofsPtr) {
223 ents_vec.push_back(dof->getEnt());
224 }
225 std::sort(ents_vec.begin(), ents_vec.end());
226 auto it = std::unique(ents_vec.begin(), ents_vec.end());
227 Range r;
228 r.insert_list(ents_vec.begin(), it);
229 return r;
230};
231
232template <typename PARENT> struct ExtractParentType { using Prent = PARENT; };
233
234#include <FreeSurfaceOps.hpp>
235using namespace FreeSurfaceOps;
236
238
239 FreeSurface(MoFEM::Interface &m_field) : mField(m_field) {}
240
242
244
246
247private:
253
254 boost::shared_ptr<FEMethod> domianLhsFEPtr;
255};
256
257//! [Run programme]
266}
267//! [Run programme]
268
269//! [Read mesh]
272
274
275 CHKERR simple->getOptions();
276 CHKERR simple->loadFile();
277
279}
280//! [Read mesh]
281
282//! [Set up problem]
285
287
288 // Fields on domain
289
290 // Velocity field
291 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, U_FIELD_DIM);
292 // Pressure field
293 CHKERR simple->addDomainField("P", H1, AINSWORTH_LEGENDRE_BASE, 1);
294 // Order/phase fild
295 CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
296 // Chemical potential
297 CHKERR simple->addDomainField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
298
299 // Field on boundary
300 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE,
302 CHKERR simple->addBoundaryField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
303 // Lagrange multiplier which constrains slip conditions
304 CHKERR simple->addBoundaryField("L", H1, AINSWORTH_LEGENDRE_BASE, 1);
305
306 CHKERR simple->setFieldOrder("U", order);
307 CHKERR simple->setFieldOrder("P", order - 1);
308 CHKERR simple->setFieldOrder("H", order);
309 CHKERR simple->setFieldOrder("G", order);
310 CHKERR simple->setFieldOrder("L", order);
311
312 CHKERR simple->setUp();
313
315}
316//! [Set up problem]
317
318//! [Boundary condition]
321
323 auto pipeline_mng = mField.getInterface<PipelineManager>();
324 auto bc_mng = mField.getInterface<BcManager>();
325 auto dm = simple->getDM();
326
327 auto h_ptr = boost::make_shared<VectorDouble>();
328 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
329 auto g_ptr = boost::make_shared<VectorDouble>();
330 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
331
332 auto set_generic = [&](auto &pipeline, auto &fe) {
333 auto det_ptr = boost::make_shared<VectorDouble>();
334 auto jac_ptr = boost::make_shared<MatrixDouble>();
335 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
336 pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
337 pipeline.push_back(
338 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
339 pipeline.push_back(
340 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
341
342 pipeline.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
343 pipeline.push_back(
344 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
345
346 pipeline.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
347 pipeline.push_back(
348 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
349 };
350
351 auto post_proc = [&]() {
353 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
354
355 set_generic(post_proc_fe->getOpPtrVector(), post_proc_fe);
356
358
359 post_proc_fe->getOpPtrVector().push_back(
360
361 new OpPPMap(
362 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
363
364 {{"H", h_ptr}, {"G", g_ptr}},
365
366 {{"GRAD_H", grad_h_ptr}, {"GRAD_G", grad_g_ptr}},
367
368 {},
369
370 {}
371
372 )
373
374 );
375
376 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
377 CHKERR post_proc_fe->writeFile("out_init.h5m");
378
380 };
381
382 auto solve_init = [&]() {
384
385 auto set_domain_rhs = [&](auto &pipeline, auto &fe) {
386 set_generic(pipeline, fe);
387 pipeline.push_back(new OpRhsH<true>("H", nullptr, nullptr, h_ptr,
388 grad_h_ptr, grad_g_ptr));
389 pipeline.push_back(new OpRhsG<true>("G", h_ptr, grad_h_ptr, g_ptr));
390 };
391
392 auto set_domain_lhs = [&](auto &pipeline, auto &fe) {
393 set_generic(pipeline, fe);
394 pipeline.push_back(new OpLhsH_dH<true>("H", nullptr, h_ptr, grad_g_ptr));
395 pipeline.push_back(new OpLhsH_dG<true>("H", "G", h_ptr));
396 pipeline.push_back(new OpLhsG_dG("G"));
397 pipeline.push_back(new OpLhsG_dH<true>("G", "H", h_ptr));
398 };
399
400 auto create_subdm = [&]() {
401 DM subdm;
402 CHKERR DMCreate(mField.get_comm(), &subdm);
403 CHKERR DMSetType(subdm, "DMMOFEM");
404 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB");
405 CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName().c_str());
406 CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
407 CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_TRUE);
408 CHKERR DMMoFEMAddSubFieldRow(subdm, "H");
409 CHKERR DMMoFEMAddSubFieldRow(subdm, "G");
410 CHKERR DMMoFEMAddSubFieldCol(subdm, "H");
411 CHKERR DMMoFEMAddSubFieldCol(subdm, "G");
412 CHKERR DMSetUp(subdm);
413 return SmartPetscObj<DM>(subdm);
414 };
415
416 auto subdm = create_subdm();
417
418 auto prb_ents = get_dofs_ents(subdm);
419
420 pipeline_mng->getDomainRhsFE().reset();
421 pipeline_mng->getDomainLhsFE().reset();
422 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
423 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
424
425 set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline(),
426 pipeline_mng->getDomainRhsFE());
427 set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline(),
428 pipeline_mng->getDomainLhsFE());
429
430 auto D = smartCreateDMVector(subdm);
431 auto snes = pipeline_mng->createSNES(subdm);
432 auto snes_ctx_ptr = smartGetDMSnesCtx(subdm);
433
434 auto set_section_monitor = [&](auto solver) {
436 CHKERR SNESMonitorSet(snes,
437 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
439 (void *)(snes_ctx_ptr.get()), nullptr);
440 auto section = mField.getInterface<ISManager>()->sectionCreate("SUB");
441 PetscInt num_fields;
442 CHKERR PetscSectionGetNumFields(section, &num_fields);
443 for (int f = 0; f < num_fields; ++f) {
444 const char *field_name;
445 CHKERR PetscSectionGetFieldName(section, f, &field_name);
446 MOFEM_LOG("FS", Sev::inform)
447 << "Field " << f << " " << std::string(field_name);
448 }
450 };
451
452 CHKERR set_section_monitor(snes);
453
454 CHKERR SNESSetFromOptions(snes);
455 CHKERR SNESSolve(snes, PETSC_NULL, D);
456
457 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
458 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
459 CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
460
462 };
463
464 CHKERR solve_init();
465 CHKERR post_proc();
466
467 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYMETRY",
468 "U", 0, 0);
469 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYMETRY",
470 "L", 0, 0);
471 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "U",
472 0, SPACE_DIM);
473 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "L",
474 0, 0);
475 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "ZERO",
476 "L", 0, 0);
477
478 // Clear pipelines
479 pipeline_mng->getOpDomainRhsPipeline().clear();
480 pipeline_mng->getOpDomainLhsPipeline().clear();
481
483}
484//! [Boundary condition]
485
486//! [Push operators to pipeline]
489
490 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
491 auto u_ptr = boost::make_shared<MatrixDouble>();
492 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
493 auto dot_h_ptr = boost::make_shared<VectorDouble>();
494 auto h_ptr = boost::make_shared<VectorDouble>();
495 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
496 auto g_ptr = boost::make_shared<VectorDouble>();
497 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
498 auto lambda_ptr = boost::make_shared<VectorDouble>();
499 auto p_ptr = boost::make_shared<VectorDouble>();
500 auto div_u_ptr = boost::make_shared<VectorDouble>();
501
502 // Push element from reference configuration to current configuration in 3d
503 // space
504 auto set_domain_general = [&](auto &pipeline, auto &fe) {
505 auto det_ptr = boost::make_shared<VectorDouble>();
506 auto jac_ptr = boost::make_shared<MatrixDouble>();
507 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
508 pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
509 pipeline.push_back(
510 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
511 pipeline.push_back(
512 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
513 pipeline.push_back(new OpSetHOWeightsOnFace());
514
515 pipeline.push_back(
517 pipeline.push_back(
519 pipeline.push_back(
521 grad_u_ptr));
522 pipeline.push_back(
524 "U", div_u_ptr));
525
526 pipeline.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
527 pipeline.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
528 pipeline.push_back(
529 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
530
531 pipeline.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
532 pipeline.push_back(
533 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
534
535 pipeline.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
536 };
537
538 auto set_domain_rhs = [&](auto &pipeline, auto &fe) {
539 set_domain_general(pipeline, fe);
540 pipeline.push_back(new OpRhsU("U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
541 grad_h_ptr, g_ptr, p_ptr));
542 pipeline.push_back(new OpRhsH<false>("H", u_ptr, dot_h_ptr, h_ptr,
543 grad_h_ptr, grad_g_ptr));
544 pipeline.push_back(new OpRhsG<false>("G", h_ptr, grad_h_ptr, g_ptr));
545 pipeline.push_back(new OpBaseTimesScalar(
546 "P", div_u_ptr, [](const double r, const double, const double) {
547 return cylindrical(r);
548 }));
549 pipeline.push_back(new OpBaseTimesScalar(
550 "P", p_ptr, [](const double r, const double, const double) {
551 return eps * cylindrical(r);
552 }));
553 };
554
555 auto set_domain_lhs = [&](auto &pipeline, auto &fe) {
556 set_domain_general(pipeline, fe);
557 pipeline.push_back(new OpLhsU_dU("U", u_ptr, grad_u_ptr, h_ptr));
558 pipeline.push_back(
559 new OpLhsU_dH("U", "H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
560 pipeline.push_back(new OpLhsU_dG("U", "G", grad_h_ptr));
561
562 pipeline.push_back(new OpLhsH_dU("H", "U", grad_h_ptr));
563 pipeline.push_back(new OpLhsH_dH<false>("H", u_ptr, h_ptr, grad_g_ptr));
564 pipeline.push_back(new OpLhsH_dG<false>("H", "G", h_ptr));
565
566 pipeline.push_back(new OpLhsG_dH<false>("G", "H", h_ptr));
567 pipeline.push_back(new OpLhsG_dG("G"));
568
569 pipeline.push_back(new OpMixScalarTimesDiv(
570 "P", "U",
571 [](const double r, const double, const double) {
572 return cylindrical(r);
573 },
574 true, false));
575 pipeline.push_back(
576 new OpDomainMassP("P", "P", [](double r, double, double) {
577 return eps * cylindrical(r);
578 }));
579 };
580
581 auto set_boundary_rhs = [&](auto &pipeline, auto &fe) {
582 pipeline.push_back(
584 pipeline.push_back(new OpCalculateScalarFieldValues("L", lambda_ptr));
585 pipeline.push_back(new OpNormalConstrainRhs("L", u_ptr));
586 pipeline.push_back(new OpNormalForceRhs("U", lambda_ptr));
587 };
588
589 auto set_boundary_lhs = [&](auto &pipeline, auto &fe) {
590 pipeline.push_back(new OpNormalConstrainLhs("L", "U"));
591 };
592
593 auto *pipeline_mng = mField.getInterface<PipelineManager>();
594
595 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
596 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
597
598 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
599 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
600
601 set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline(),
602 pipeline_mng->getDomainRhsFE());
603 set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline(),
604 pipeline_mng->getDomainLhsFE());
605
606 set_boundary_rhs(pipeline_mng->getOpBoundaryRhsPipeline(),
607 pipeline_mng->getBoundaryRhsFE());
608 set_boundary_lhs(pipeline_mng->getOpBoundaryLhsPipeline(),
609 pipeline_mng->getBoundaryLhsFE());
610
611 domianLhsFEPtr = pipeline_mng->getDomainLhsFE();
612
614}
615//! [Push operators to pipeline]
616
617/**
618 * @brief Monitor solution
619 *
620 * This functions is called by TS solver at the end of each step. It is used
621 * to output results to the hard drive.
622 */
623struct Monitor : public FEMethod {
625 SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc,
626 boost::shared_ptr<PostProcEdgeEle> post_proc_edge,
627 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
628 p)
629 : dM(dm), postProc(post_proc), postProcEdge(post_proc_edge),
630 liftFE(p.first), liftVec(p.second) {}
633 constexpr int save_every_nth_step = 1;
634 if (ts_step % save_every_nth_step == 0) {
636 this->getCacheWeakPtr());
637 CHKERR postProc->writeFile(
638 "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
639
641 this->getCacheWeakPtr());
642 CHKERR postProcEdge->writeFile(
643 "out_step_bdy_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
644 }
645
646 liftVec->resize(SPACE_DIM, false);
647 liftVec->clear();
649 MPI_Allreduce(MPI_IN_PLACE, &(*liftVec)[0], SPACE_DIM, MPI_DOUBLE, MPI_SUM,
650 MPI_COMM_WORLD);
651 MOFEM_LOG("FS", Sev::inform)
652 << "Step " << ts_step << " time " << ts_t
653 << " lift vec x: " << (*liftVec)[0] << " y: " << (*liftVec)[1];
654
656 }
657
658private:
660 boost::shared_ptr<PostProcEle> postProc;
661 boost::shared_ptr<PostProcEdgeEle> postProcEdge;
662 boost::shared_ptr<BoundaryEle> liftFE;
663 boost::shared_ptr<VectorDouble> liftVec;
664};
665
666//! [Solve]
669
670 auto *simple = mField.getInterface<Simple>();
671 auto *pipeline_mng = mField.getInterface<PipelineManager>();
672 auto dm = simple->getDM();
673
674 auto get_fe_post_proc = [&]() {
675 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
676
677 auto det_ptr = boost::make_shared<VectorDouble>();
678 auto jac_ptr = boost::make_shared<MatrixDouble>();
679 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
680
681 auto u_ptr = boost::make_shared<MatrixDouble>();
682 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
683 auto h_ptr = boost::make_shared<VectorDouble>();
684 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
685 auto p_ptr = boost::make_shared<VectorDouble>();
686 auto g_ptr = boost::make_shared<VectorDouble>();
687 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
688
689 post_proc_fe->getOpPtrVector().push_back(
690 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
691 post_proc_fe->getOpPtrVector().push_back(
692 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
693 post_proc_fe->getOpPtrVector().push_back(
694 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
695
696 post_proc_fe->getOpPtrVector().push_back(
698 post_proc_fe->getOpPtrVector().push_back(
700 grad_u_ptr));
701
702 post_proc_fe->getOpPtrVector().push_back(
703 new OpCalculateScalarFieldValues("H", h_ptr));
704 post_proc_fe->getOpPtrVector().push_back(
705 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
706
707 post_proc_fe->getOpPtrVector().push_back(
708 new OpCalculateScalarFieldValues("P", p_ptr));
709 post_proc_fe->getOpPtrVector().push_back(
710 new OpCalculateScalarFieldValues("G", g_ptr));
711 post_proc_fe->getOpPtrVector().push_back(
712 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
713
715
716 post_proc_fe->getOpPtrVector().push_back(
717
718 new OpPPMap(
719 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
720
721 {{"H", h_ptr}, {"P", p_ptr}, {"G", g_ptr}},
722
723 {{"U", u_ptr}, {"H_GRAD", grad_h_ptr}, {"G_GRAD", grad_g_ptr}},
724
725 {{"GRAD_U", grad_u_ptr}},
726
727 {}
728
729 )
730
731 );
732
733 return post_proc_fe;
734 };
735
736 auto get_bdy_post_proc_fe = [&]() {
737 auto post_proc_fe = boost::make_shared<PostProcEdgeEle>(mField);
738
739 auto u_ptr = boost::make_shared<MatrixDouble>();
740 auto p_ptr = boost::make_shared<VectorDouble>();
741 auto lambda_ptr = boost::make_shared<VectorDouble>();
742
743 post_proc_fe->getOpPtrVector().push_back(
745 post_proc_fe->getOpPtrVector().push_back(
746 new OpCalculateScalarFieldValues("L", lambda_ptr));
747 post_proc_fe->getOpPtrVector().push_back(
748 new OpCalculateScalarFieldValues("P", p_ptr));
749
751
752 post_proc_fe->getOpPtrVector().push_back(
753
754 new OpPPMap(post_proc_fe->getPostProcMesh(),
755 post_proc_fe->getMapGaussPts(),
756
757 OpPPMap::DataMapVec{{"L", lambda_ptr}, {"P", p_ptr}},
758
759 OpPPMap::DataMapMat{{"U", u_ptr}},
760
762
764
765 )
766
767 );
768
769 return post_proc_fe;
770 };
771
772 auto get_lift_fe = [&]() {
773 auto fe = boost::make_shared<BoundaryEle>(mField);
774 auto lift_ptr = boost::make_shared<VectorDouble>();
775 auto p_ptr = boost::make_shared<VectorDouble>();
776 auto ents_ptr = boost::make_shared<Range>();
777
778 std::vector<const CubitMeshSets *> vec_ptr;
779 CHKERR mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
780 std::regex("LIFT"), vec_ptr);
781 for (auto m_ptr : vec_ptr) {
782 auto meshset = m_ptr->getMeshset();
783 Range ents;
784 CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
785 ents, true);
786 ents_ptr->merge(ents);
787 }
788
789 MOFEM_LOG("FS", Sev::noisy) << "Lift ents " << (*ents_ptr);
790
791 fe->getOpPtrVector().push_back(
792 new OpCalculateScalarFieldValues("P", p_ptr));
793 fe->getOpPtrVector().push_back(
794 new OpCalculateLift("P", p_ptr, lift_ptr, ents_ptr));
795
796 return std::make_pair(fe, lift_ptr);
797 };
798
799 auto set_ts = [&](auto solver) {
801 SNES snes;
802 CHKERR TSGetSNES(solver, &snes);
803 KSP ksp;
804 CHKERR SNESGetKSP(snes, &ksp);
806 };
807
808 auto ts = pipeline_mng->createTSIM();
809 CHKERR TSSetType(ts, TSALPHA);
810
811 auto set_post_proc_monitor = [&](auto dm) {
813 boost::shared_ptr<FEMethod> null_fe;
814 auto monitor_ptr = boost::make_shared<Monitor>(
815 dm, get_fe_post_proc(), get_bdy_post_proc_fe(), get_lift_fe());
816 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
817 null_fe, monitor_ptr);
819 };
820 CHKERR set_post_proc_monitor(dm);
821
822 // Add monitor to time solver
823 double ftime = 1;
824 // CHKERR TSSetMaxTime(ts, ftime);
825 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
826
827 auto T = smartCreateDMVector(simple->getDM());
828 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
829 SCATTER_FORWARD);
830 CHKERR TSSetSolution(ts, T);
831 CHKERR TSSetFromOptions(ts);
832 CHKERR set_ts(ts);
833 CHKERR TSSetUp(ts);
834
835 auto print_fields_in_section = [&]() {
837
838 auto section = mField.getInterface<ISManager>()->sectionCreate(
839 simple->getProblemName());
840 PetscInt num_fields;
841 CHKERR PetscSectionGetNumFields(section, &num_fields);
842 for (int f = 0; f < num_fields; ++f) {
843 const char *field_name;
844 CHKERR PetscSectionGetFieldName(section, f, &field_name);
845 MOFEM_LOG("FS", Sev::inform)
846 << "Field " << f << " " << std::string(field_name);
847 }
849 };
850
851 CHKERR print_fields_in_section();
852
853 CHKERR TSSolve(ts, NULL);
854 CHKERR TSGetTime(ts, &ftime);
855
857}
858
859int main(int argc, char *argv[]) {
860
861 // Initialisation of MoFEM/PETSc and MOAB data structures
862 const char param_file[] = "param_file.petsc";
864
865 // Add logging channel for example
866 auto core_log = logging::core::get();
867 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FS"));
868 LogManager::setLog("FS");
869 MOFEM_LOG_TAG("FS", "free surface");
870
871 try {
872
873 //! [Register MoFEM discrete manager in PETSc]
874 DMType dm_name = "DMMOFEM";
875 CHKERR DMRegister_MoFEM(dm_name);
876 //! [Register MoFEM discrete manager in PETSc
877
878 //! [Create MoAB]
879 moab::Core mb_instance; ///< mesh database
880 moab::Interface &moab = mb_instance; ///< mesh database interface
881 //! [Create MoAB]
882
883 //! [Create MoFEM]
884 MoFEM::Core core(moab); ///< finite element database
885 MoFEM::Interface &m_field = core; ///< finite element database insterface
886 //! [Create MoFEM]
887
888 //! [FreeSurface]
889 FreeSurface ex(m_field);
890 CHKERR ex.runProblem();
891 //! [FreeSurface]
892 }
894
896}
static Index< 'p', 3 > p
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
CoordinateTypes
Coodinate system.
Definition: definitions.h:114
@ CYLINDRICAL
Definition: definitions.h:117
@ CARTESIAN
Definition: definitions.h:115
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'n', SPACE_DIM > n
auto cylindrical
auto marker
int main(int argc, char *argv[])
auto save_range
auto get_M2
EntitiesFieldData::EntData EntData
constexpr int U_FIELD_DIM
constexpr double mu_m
constexpr double rho_ave
FTensor::Index< 'j', SPACE_DIM > j
auto get_M2_dh
constexpr double rho_m
auto integration_rule
static char help[]
constexpr CoordinateTypes coord_type
select coordinate system <CARTESIAN, CYLINDRICAL>;
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, U_FIELD_DIM > OpDomainSourceU
constexpr double mu_ave
auto get_M_dh
auto kernel_eye
constexpr double mu_p
constexpr double eps
FTensor::Index< 'k', SPACE_DIM > k
constexpr double lambda
auto get_M3_dh
FTensor::Index< 'i', SPACE_DIM > i
constexpr double mu_diff
auto init_h
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
constexpr auto t_kd
auto d_phase_function_h
constexpr int SPACE_DIM
constexpr double eta
FTensor::Index< 'l', SPACE_DIM > l
auto kernel_oscillation
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, 1 > OpDomainMassH
auto get_M0
constexpr double a0
constexpr double tol
auto get_f_dh
auto get_M3
OpDomainMassH OpDomainMassP
auto get_dofs_ents
constexpr int BASE_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1, 1 > OpBaseTimesScalar
constexpr double W
auto get_M
auto phase_function
constexpr double h
constexpr int powof2()
ElementsAndOps< SPACE_DIM >::DomianParentEle DomianParentEle
constexpr double rho_p
constexpr double rho_diff
auto get_D
auto get_M0_dh
auto my_max
auto cut_off
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, 1 > OpDomainSourceH
constexpr double eta2
constexpr int order
approximation order
const double kappa
auto bit
constexpr double md
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
auto get_f
auto d_cut_off
auto my_min
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMMoFEM.cpp:201
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMMoFEM.cpp:403
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:223
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:450
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:470
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:246
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
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:533
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
const double T
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
auto f
Definition: HenckyOps.hpp:5
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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:1003
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMMoFEM.cpp:385
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:224
CoreTmp< 0 > Core
Definition: Core.hpp:1086
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:941
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:987
double w(const double g, const double t)
const double D
diffusivity
const double r
rate factor
double A
int save_every_nth_step
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
FreeSurface(MoFEM::Interface &m_field)
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEM::Interface & mField
boost::shared_ptr< FEMethod > domianLhsFEPtr
MoFEMErrorCode makeRefProblem()
MoFEMErrorCode solveSystem()
[Solve]
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:23
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Interface for managing meshsets containing materials and boundary conditions.
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar 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.
Post post-proc data at points from hash maps.
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
PetscReal ts_t
time
PetscInt ts_step
time step
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.
boost::shared_ptr< PostProcEdgeEle > postProcEdge
SmartPetscObj< DM > dM
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< VectorDouble > liftVec
boost::shared_ptr< BoundaryEle > liftFE
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc, boost::shared_ptr< PostProcEdgeEle > post_proc_edge, std::pair< boost::shared_ptr< BoundaryEle >, boost::shared_ptr< VectorDouble > > p)
boost::shared_ptr< PostProcEle > postProc
Postprocess on edge.