v0.13.2
Loading...
Searching...
No Matches
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 * Implementation based on \cite Lovric2019-qn
10 */
11
12#include <MoFEM.hpp>
13#include <petsc/private/petscimpl.h>
14
15using namespace MoFEM;
16
17static char help[] = "...\n\n";
18
19#ifdef PYTHON_INIT_SURFACE
20#include <boost/python.hpp>
21#include <boost/python/def.hpp>
22namespace bp = boost::python;
23
24struct SurfacePython {
25 SurfacePython() = default;
26 virtual ~SurfacePython() = default;
27
28 MoFEMErrorCode surfaceInit(const std::string py_file) {
30 try {
31
32 // create main module
33 auto main_module = bp::import("__main__");
34 mainNamespace = main_module.attr("__dict__");
35 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
36 // create a reference to python function
37 surfaceFun = mainNamespace["surface"];
38
39 } catch (bp::error_already_set const &) {
40 // print all other errors to stderr
41 PyErr_Print();
43 }
45 };
46
47 MoFEMErrorCode evalSurface(
48
49 double x, double y, double z, double eta, double &s
50
51 ) {
53 try {
54
55 // call python function
56 s = bp::extract<double>(surfaceFun(x, y, z, eta));
57
58 } catch (bp::error_already_set const &) {
59 // print all other errors to stderr
60 PyErr_Print();
62 }
64 }
65
66private:
67 bp::object mainNamespace;
68 bp::object surfaceFun;
69};
70
71static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
72
73#endif
74
75constexpr int BASE_DIM = 1;
76constexpr int SPACE_DIM = 2;
77constexpr int U_FIELD_DIM = SPACE_DIM;
79 EXECUTABLE_COORD_TYPE; ///< select coordinate system <CARTESIAN,
80 ///< CYLINDRICAL>;
81
82constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
83constexpr IntegrationType I =
84 IntegrationType::GAUSS; //< selected integration type
85
86template <int DIM>
88
96using SideOp = SideEle::UserDataOperator;
97
101
103
107
116
123
126
127// Flux is applied by Lagrange Multiplie
131
136
138
139// mesh refinement
140int order = 3; ///< approximation order
141int nb_levels = 4; //< number of refinement levels
142int refine_overlap = 4; //< mesh overlap while refine
143
144constexpr bool debug = true;
145
146auto get_start_bit = []() {
147 return nb_levels + 1;
148}; //< first refinement level for computational mesh
149auto get_current_bit = []() {
150 return 2 * get_start_bit() + 1;
151}; ///< dofs bit used to do calculations
152auto get_skin_parent_bit = []() { return 2 * get_start_bit() + 2; };
153auto get_skin_child_bit = []() { return 2 * get_start_bit() + 3; };
154auto get_projection_bit = []() { return 2 * get_start_bit() + 4; };
155auto get_skin_projection_bit = []() { return 2 * get_start_bit() + 5; };
156
157// FIXME: Set parameters from command line
158
159// Physical parameters
160constexpr double a0 = 0; // 980;
161constexpr double rho_m = 0.998;
162constexpr double mu_m = 0.010101;
163constexpr double rho_p = 0.0012;
164constexpr double mu_p = 0.000182;
165constexpr double lambda = 73 / 4; ///< surface tension
166constexpr double W = 0.25;
167
168// Model parameters
169constexpr double h = 0.0015 / 20; // mesh size
170constexpr double eta = h;
171constexpr double eta2 = eta * eta;
172
173// Numerical parameters
174constexpr double md = 1e-2;
175constexpr double eps = 1e-12;
176constexpr double tol = std::numeric_limits<float>::epsilon();
177
178constexpr double rho_ave = (rho_p + rho_m) / 2;
179constexpr double rho_diff = (rho_p - rho_m) / 2;
180constexpr double mu_ave = (mu_p + mu_m) / 2;
181constexpr double mu_diff = (mu_p - mu_m) / 2;
182
183const double kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
184
185auto integration_rule = [](int, int, int) { return 2 * order + 1; };
186
187auto cylindrical = [](const double r) {
188 // When we move to C++17 add if constexpr()
189 if constexpr (coord_type == CYLINDRICAL)
190 return 2 * M_PI * r;
191 else
192 return 1.;
193};
194
195auto wetting_angle_sub_stepping = [](auto ts_step) {
196 constexpr int sub_stepping = 16;
197 return std::min(1., static_cast<double>(ts_step) / sub_stepping);
198};
199
200// cut off function
201
202auto my_max = [](const double x) { return (x - 1 + std::abs(x + 1)) / 2; };
203auto my_min = [](const double x) { return (x + 1 - std::abs(x - 1)) / 2; };
204auto cut_off = [](const double h) { return my_max(my_min(h)); };
205auto d_cut_off = [](const double h) {
206 if (h >= -1 && h < 1)
207 return 1.;
208 else
209 return 0.;
210};
211
212auto phase_function = [](const double h, const double diff, const double ave) {
213 return diff * cut_off(h) + ave;
214};
215
216auto d_phase_function_h = [](const double h, const double diff) {
217 return diff * d_cut_off(h);
218};
219
220// order function (free energy)
221
222auto get_f = [](const double h) { return 4 * W * h * (h * h - 1); };
223auto get_f_dh = [](const double h) { return 4 * W * (3 * h * h - 1); };
224
225auto get_M0 = [](auto h) { return md; };
226auto get_M0_dh = [](auto h) { return 0; };
227
228auto get_M2 = [](auto h_tmp) {
229 const double h = cut_off(h_tmp);
230 return md * (1 - h * h);
231};
232
233auto get_M2_dh = [](auto h_tmp) {
234 const double h = cut_off(h_tmp);
235 return -md * 2 * h * d_cut_off(h_tmp);
236};
237
238auto get_M3 = [](auto h_tmp) {
239 const double h = cut_off(h_tmp);
240 const double h2 = h * h;
241 const double h3 = h2 * h;
242 if (h >= 0)
243 return md * (2 * h3 - 3 * h2 + 1);
244 else
245 return md * (-2 * h3 - 3 * h2 + 1);
246};
247
248auto get_M3_dh = [](auto h_tmp) {
249 const double h = cut_off(h_tmp);
250 if (h >= 0)
251 return md * (6 * h * (h - 1)) * d_cut_off(h_tmp);
252 else
253 return md * (-6 * h * (h + 1)) * d_cut_off(h_tmp);
254};
255
256auto get_M = [](auto h) { return get_M0(h); };
257auto get_M_dh = [](auto h) { return get_M0_dh(h); };
258
259auto get_D = [](const double A) {
261 t_D(i, j, k, l) = A * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
262 return t_D;
263};
264
265// some examples of initialisation functions
266
267auto kernel_oscillation = [](double r, double y, double) {
268 constexpr int n = 3;
269 constexpr double R = 0.0125;
270 constexpr double A = R * 0.2;
271 const double theta = atan2(r, y);
272 const double w = R + A * cos(n * theta);
273 const double d = std::sqrt(r * r + y * y);
274 return tanh((w - d) / (eta * std::sqrt(2)));
275};
276
277auto kernel_eye = [](double r, double y, double) {
278 constexpr double y0 = 0.5;
279 constexpr double R = 0.4;
280 y -= y0;
281 const double d = std::sqrt(r * r + y * y);
282 return tanh((R - d) / (eta * std::sqrt(2)));
283};
284
285auto capillary_tube = [](double x, double y, double z) {
286 constexpr double water_height = 0.;
287 return tanh((water_height - y) / (eta * std::sqrt(2)));
288 ;
289};
290
291auto bubble_device = [](double x, double y, double z) {
292 return -tanh((-0.039 - x) / (eta * std::sqrt(2)));
293};
294
295/**
296 * @brief Initialisation function
297 *
298 * @note If UMs are compiled with Python to initialise phase field "H"
299 * surface.py function is used, which has to be present in execution folder.
300 *
301 */
302auto init_h = [](double r, double y, double theta) {
303#ifdef PYTHON_INIT_SURFACE
304 double s = 1;
305 if (auto ptr = surfacePythonWeakPtr.lock()) {
306 CHK_THROW_MESSAGE(ptr->evalSurface(r, y, theta, eta, s),
307 "error eval python");
308 }
309 return s;
310#else
311 return bubble_device(r, y, theta);
312 // return capillary_tube(r, y, theta);
313 // return kernel_eye(r, y, theta);
314#endif
315};
316
317auto wetting_angle = [](double water_level) { return water_level; };
318
319auto bit = [](auto b) { return BitRefLevel().set(b); };
320auto marker = [](auto b) { return BitRefLevel().set(BITREFLEVEL_SIZE - b); };
321auto get_fe_bit = [](FEMethod *fe_ptr) {
322 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
323};
324
325auto get_global_size = [](int l_size) {
326 int g_size;
327 MPI_Allreduce(&l_size, &g_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
328 return g_size;
329};
330
331auto save_range = [](moab::Interface &moab, const std::string name,
332 const Range r) {
334 if (get_global_size(r.size())) {
335 auto out_meshset = get_temp_meshset_ptr(moab);
336 CHKERR moab.add_entities(*out_meshset, r);
337 CHKERR moab.write_file(name.c_str(), "MOAB", "PARALLEL=WRITE_PART",
338 out_meshset->get_ptr(), 1);
339 }
341};
342
343/**
344 * @brief get entities of dofs in the problem - used for debugging
345 *
346 */
347auto get_dofs_ents_by_field_name = [](auto dm, auto field_name) {
348 auto prb_ptr = getProblemPtr(dm);
349 std::vector<EntityHandle> ents_vec;
350
351 MoFEM::Interface *m_field_ptr;
352 CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
353
354 auto bit_number = m_field_ptr->get_field_bit_number(field_name);
355 auto dofs = prb_ptr->numeredRowDofsPtr;
356 auto lo_it = dofs->lower_bound(FieldEntity::getLoBitNumberUId(bit_number));
357 auto hi_it = dofs->upper_bound(FieldEntity::getHiBitNumberUId(bit_number));
358 ents_vec.reserve(std::distance(lo_it, hi_it));
359
360 for (; lo_it != hi_it; ++lo_it) {
361 ents_vec.push_back((*lo_it)->getEnt());
362 }
363
364 std::sort(ents_vec.begin(), ents_vec.end());
365 auto it = std::unique(ents_vec.begin(), ents_vec.end());
366 Range r;
367 r.insert_list(ents_vec.begin(), it);
368 return r;
369};
370
371/**
372 * @brief get entities of dofs in the problem - used for debugging
373 *
374 */
375auto get_dofs_ents_all = [](auto dm) {
376 auto prb_ptr = getProblemPtr(dm);
377 std::vector<EntityHandle> ents_vec;
378
379 MoFEM::Interface *m_field_ptr;
380 CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
381
382 auto dofs = prb_ptr->numeredRowDofsPtr;
383 ents_vec.reserve(dofs->size());
384
385 for (auto d : *dofs) {
386 ents_vec.push_back(d->getEnt());
387 }
388
389 std::sort(ents_vec.begin(), ents_vec.end());
390 auto it = std::unique(ents_vec.begin(), ents_vec.end());
391 Range r;
392 r.insert_list(ents_vec.begin(), it);
393 return r;
394};
395
396#include <FreeSurfaceOps.hpp>
397using namespace FreeSurfaceOps;
398
399struct FreeSurface;
400
401enum FR { F, R }; // F - forward, and reverse
402
403/**
404 * @brief Set of functions called by PETSc solver used to refine and update
405 * mesh.
406 *
407 * @note Currently theta method is only handled by this code.
408 *
409 */
411
412 TSPrePostProc() = default;
413 virtual ~TSPrePostProc() = default;
414
415 /**
416 * @brief Used to setup TS solver
417 *
418 * @param ts
419 * @return MoFEMErrorCode
420 */
421 MoFEMErrorCode tsSetUp(TS ts);
422
423 SmartPetscObj<VecScatter> getScatter(Vec x, Vec y, enum FR fr);
425
429
430private:
431 /**
432 * @brief Pre process time step
433 *
434 * Refine mesh and update fields
435 *
436 * @param ts
437 * @return MoFEMErrorCode
438 */
439 static MoFEMErrorCode tsPreProc(TS ts);
440
441 /**
442 * @brief Post process time step.
443 *
444 * Currently that function do not make anything major
445 *
446 * @param ts
447 * @return MoFEMErrorCode
448 */
449 static MoFEMErrorCode tsPostProc(TS ts);
450
452
453 static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t,
454 Vec f,
455 void *ctx); //< Wrapper for SNES Rhs
456 static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
457 PetscReal a, Mat A, Mat B,
458 void *ctx); ///< Wrapper for SNES Lhs
459 static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u,
460 void *ctx); ///< Wrapper for TS monitor
461 static MoFEMErrorCode pcSetup(PC pc);
462 static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x);
463
464 SmartPetscObj<Vec> globRes; //< global residual
465 SmartPetscObj<Mat> subB; //< sub problem tangent matrix
466 SmartPetscObj<KSP> subKSP; //< sub problem KSP solver
467
468 boost::shared_ptr<SnesCtx>
469 snesCtxPtr; //< infernal data (context) for MoFEM SNES fuctions
470 boost::shared_ptr<TsCtx>
471 tsCtxPtr; //< internal data (context) for MoFEM TS functions.
472};
473
474static boost::weak_ptr<TSPrePostProc> tsPrePostProc;
475
477
478 FreeSurface(MoFEM::Interface &m_field) : mField(m_field) {}
479
481
483
485
486private:
493
494 /// @brief Find entities on refinement levels
495 /// @param overlap level of overlap
496 /// @return
497 std::vector<Range> findEntitiesCrossedByPhaseInterface(size_t overlap);
498
499 /// @brief
500 /// @param ents
501 /// @param level
502 /// @param mask
503 /// @return
505
506 /// @brief
507 /// @param overlap
508 /// @return
509 MoFEMErrorCode refineMesh(size_t overlap);
510
511 /// @brief Create hierarchy of elements run on parents levels
512 /// @param fe_top pipeline element to which hierarchy is attached
513 /// @param op type of operator OPSPACE, OPROW, OPCOL or OPROWCOL
514 /// @param get_elema lambda function to create element on hierarchy
515 /// @param verbosity verbosity level
516 /// @param sev severity level
517 /// @return
519 boost::shared_ptr<FEMethod> fe_top, std::string field_name,
520 ForcesAndSourcesCore::UserDataOperator::OpType op,
521 BitRefLevel child_ent_bit,
522 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
523 int verbosity, LogManager::SeverityLevel sev);
524
525 friend struct TSPrePostProc;
526};
527
528//! [Run programme]
537}
538//! [Run programme]
539
540//! [Read mesh]
543 MOFEM_LOG("FS", Sev::inform)
544 << "Read mesh for problem in " << EXECUTABLE_COORD_TYPE;
546
547 simple->getParentAdjacencies() = true;
548 simple->getBitRefLevel() = BitRefLevel();
549
550 CHKERR simple->getOptions();
551 CHKERR simple->loadFile();
552
554}
555//! [Read mesh]
556
557//! [Set up problem]
560
561 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
562 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-nb_levels", &nb_levels,
563 PETSC_NULL);
564 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-refine_overlap", &refine_overlap,
565 PETSC_NULL);
566
567 MOFEM_LOG("FS", Sev::inform) << "order = " << order;
568 MOFEM_LOG("FS", Sev::inform) << "nb_levels = " << nb_levels;
569 nb_levels += 1;
570
572 auto bit_mng = mField.getInterface<BitRefManager>();
573
574 // Fields on domain
575
576 // Velocity field
577 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, U_FIELD_DIM);
578 // Pressure field
579 CHKERR simple->addDomainField("P", H1, AINSWORTH_LEGENDRE_BASE, 1);
580 // Order/phase fild
581 CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
582 // Chemical potential
583 CHKERR simple->addDomainField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
584
585 // Field on boundary
586 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE,
588 CHKERR simple->addBoundaryField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
589 CHKERR simple->addBoundaryField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
590 // Lagrange multiplier which constrains slip conditions
591 CHKERR simple->addBoundaryField("L", H1, AINSWORTH_LEGENDRE_BASE, 1);
592
593 CHKERR simple->setFieldOrder("U", order);
594 CHKERR simple->setFieldOrder("P", order - 1);
595 CHKERR simple->setFieldOrder("H", order);
596 CHKERR simple->setFieldOrder("G", order);
597 CHKERR simple->setFieldOrder("L", order);
598
599 // Initialise bit ref levels
600 auto set_problem_bit = [&]() {
602 // Set bits to build adjacencies between parents and children. That is
603 // used by simple interface.
604 simple->getBitAdjEnt() = BitRefLevel().set();
605 simple->getBitAdjParent() = BitRefLevel().set();
606 simple->getBitRefLevel() = BitRefLevel().set();
607 simple->getBitRefLevelMask() = BitRefLevel().set();
609 };
610
611 CHKERR set_problem_bit();
612
613 CHKERR simple->setUp();
614
616}
617//! [Set up problem]
618
619//! [Boundary condition]
622
623#ifdef PYTHON_INIT_SURFACE
624 auto get_py_surface_init = []() {
625 auto py_surf_init = boost::make_shared<SurfacePython>();
626 CHKERR py_surf_init->surfaceInit("surface.py");
627 surfacePythonWeakPtr = py_surf_init;
628 return py_surf_init;
629 };
630 auto py_surf_init = get_py_surface_init();
631#endif
632
634 auto pip_mng = mField.getInterface<PipelineManager>();
635 auto bc_mng = mField.getInterface<BcManager>();
636 auto bit_mng = mField.getInterface<BitRefManager>();
637 auto dm = simple->getDM();
638
640
641 auto reset_bits = [&]() {
643 BitRefLevel start_mask;
644 for (auto s = 0; s != get_start_bit(); ++s)
645 start_mask[s] = true;
646 // reset bit ref levels
647 CHKERR bit_mng->lambdaBitRefLevel(
648 [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
649 Range level0;
650 CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
651 CHKERR bit_mng->setNthBitRefLevel(level0, get_current_bit(), true);
652 CHKERR bit_mng->setNthBitRefLevel(level0, get_projection_bit(), true);
654 };
655
656 auto add_parent_field = [&](auto fe, auto op, auto field) {
657 return setParentDofs(
658 fe, field, op, bit(get_skin_parent_bit()),
659
660 [&]() {
661 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
663 return fe_parent;
664 },
665
666 QUIET, Sev::noisy);
667 };
668
669 auto h_ptr = boost::make_shared<VectorDouble>();
670 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
671 auto g_ptr = boost::make_shared<VectorDouble>();
672 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
673
674 auto set_generic = [&](auto fe) {
676 auto &pip = fe->getOpPtrVector();
677
679
681 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
682
683 [&]() {
684 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
687 fe_parent->getOpPtrVector(), {H1});
688 return fe_parent;
689 },
690
691 QUIET, Sev::noisy);
692
693 CHKERR add_parent_field(fe, UDO::OPROW, "H");
694 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
695 pip.push_back(
696 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
697
698 CHKERR add_parent_field(fe, UDO::OPROW, "G");
699 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
700 pip.push_back(
701 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
702
704 };
705
706 auto post_proc = [&](auto exe_test) {
708 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
709 post_proc_fe->exeTestHook = exe_test;
710
711 CHKERR set_generic(post_proc_fe);
712
714
715 post_proc_fe->getOpPtrVector().push_back(
716
717 new OpPPMap(post_proc_fe->getPostProcMesh(),
718 post_proc_fe->getMapGaussPts(),
719
720 {{"H", h_ptr}, {"G", g_ptr}},
721
722 {{"GRAD_H", grad_h_ptr}, {"GRAD_G", grad_g_ptr}},
723
724 {},
725
726 {}
727
728 )
729
730 );
731
732 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
733 CHKERR post_proc_fe->writeFile("out_init.h5m");
734
736 };
737
738 auto solve_init = [&](auto exe_test) {
740
741 pip_mng->getOpDomainRhsPipeline().clear();
742 pip_mng->getOpDomainLhsPipeline().clear();
743
744 auto set_domain_rhs = [&](auto fe) {
746 CHKERR set_generic(fe);
747 auto &pip = fe->getOpPtrVector();
748
749 CHKERR add_parent_field(fe, UDO::OPROW, "H");
750 pip.push_back(new OpRhsH<true>("H", nullptr, nullptr, h_ptr, grad_h_ptr,
751 grad_g_ptr));
752 CHKERR add_parent_field(fe, UDO::OPROW, "G");
753 pip.push_back(new OpRhsG<true>("G", h_ptr, grad_h_ptr, g_ptr));
755 };
756
757 auto set_domain_lhs = [&](auto fe) {
759
760 CHKERR set_generic(fe);
761 auto &pip = fe->getOpPtrVector();
762
763 CHKERR add_parent_field(fe, UDO::OPROW, "H");
764 CHKERR add_parent_field(fe, UDO::OPCOL, "H");
765 pip.push_back(new OpLhsH_dH<true>("H", nullptr, h_ptr, grad_g_ptr));
766
767 CHKERR add_parent_field(fe, UDO::OPCOL, "G");
768 pip.push_back(new OpLhsH_dG<true>("H", "G", h_ptr));
769
770 CHKERR add_parent_field(fe, UDO::OPROW, "G");
771 pip.push_back(new OpLhsG_dG("G"));
772
773 CHKERR add_parent_field(fe, UDO::OPCOL, "H");
774 pip.push_back(new OpLhsG_dH<true>("G", "H", h_ptr));
776 };
777
778 auto create_subdm = [&]() {
779 auto level_ents_ptr = boost::make_shared<Range>();
780 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
781 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
782
783 DM subdm;
784 CHKERR DMCreate(mField.get_comm(), &subdm);
785 CHKERR DMSetType(subdm, "DMMOFEM");
786 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_INIT");
787 CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
788 CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
789 CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_TRUE);
790
791 for (auto f : {"H", "G"}) {
792 CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
793 CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
794 }
795 CHKERR DMSetUp(subdm);
796
797 if constexpr (debug) {
798 if (mField.get_comm_size() == 1) {
799 auto dm_ents = get_dofs_ents_all(SmartPetscObj<DM>(subdm, true));
800 CHKERR save_range(mField.get_moab(), "sub_dm_init_ents_verts.h5m",
801 dm_ents.subset_by_type(MBVERTEX));
802 dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
803 CHKERR save_range(mField.get_moab(), "sub_dm_init_ents.h5m", dm_ents);
804 }
805 }
806
807 return SmartPetscObj<DM>(subdm);
808 };
809
810 auto subdm = create_subdm();
811
812 pip_mng->getDomainRhsFE().reset();
813 pip_mng->getDomainLhsFE().reset();
814 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
815 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
816 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
817 pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
818
819 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
820 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
821
822 auto D = createDMVector(subdm);
823 auto snes = pip_mng->createSNES(subdm);
824 auto snes_ctx_ptr = getDMSnesCtx(subdm);
825
826 auto set_section_monitor = [&](auto solver) {
828 CHKERR SNESMonitorSet(solver,
829 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
831 (void *)(snes_ctx_ptr.get()), nullptr);
832
834 };
835
836 auto print_section_field = [&]() {
838
839 auto section =
840 mField.getInterface<ISManager>()->sectionCreate("SUB_INIT");
841 PetscInt num_fields;
842 CHKERR PetscSectionGetNumFields(section, &num_fields);
843 for (int f = 0; f < num_fields; ++f) {
844 const char *field_name;
845 CHKERR PetscSectionGetFieldName(section, f, &field_name);
846 MOFEM_LOG("FS", Sev::inform)
847 << "Field " << f << " " << std::string(field_name);
848 }
849
851 };
852
853 CHKERR set_section_monitor(snes);
854 CHKERR print_section_field();
855
856 for (auto f : {"U", "P", "H", "G", "L"}) {
857 CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
858 }
859
860 CHKERR SNESSetFromOptions(snes);
861 CHKERR SNESSolve(snes, PETSC_NULL, D);
862
863 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
864 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
865 CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
866
868 };
869
870 CHKERR reset_bits();
871 CHKERR solve_init(
872 [](FEMethod *fe_ptr) { return get_fe_bit(fe_ptr).test(0); });
873 CHKERR refineMesh(refine_overlap);
874 for (auto f : {"U", "P", "H", "G", "L"}) {
875 CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
876 }
877 CHKERR solve_init([](FEMethod *fe_ptr) {
878 return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
879 });
880 CHKERR post_proc([](FEMethod *fe_ptr) {
881 return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
882 });
883
884 pip_mng->getOpDomainRhsPipeline().clear();
885 pip_mng->getOpDomainLhsPipeline().clear();
886
887 // Remove DOFs where boundary conditions are set
888 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
889 "U", 0, 0);
890 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
891 "L", 0, 0);
892 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
893 "U", 1, 1);
894 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
895 "L", 1, 1);
896 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "U",
897 0, SPACE_DIM);
898 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "L",
899 0, 0);
900 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "ZERO",
901 "L", 0, 0);
903}
904//! [Boundary condition]
905
906//! [Data projection]
909
910 // FIXME: Functionality in this method should be improved (projection should
911 // be done field by field), generalised, and move to become core
912 // functionality.
913
915 auto pip_mng = mField.getInterface<PipelineManager>();
916 auto bc_mng = mField.getInterface<BcManager>();
917 auto bit_mng = mField.getInterface<BitRefManager>();
918 auto field_blas = mField.getInterface<FieldBlas>();
919
920 // Store all existing elements pipelines, replace them by data projection
921 // pipelines, to put them back when projection is done.
922 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
923 auto fe_domain_rhs = pip_mng->getDomainRhsFE();
924 auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
925 auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
926
927 pip_mng->getDomainLhsFE().reset();
928 pip_mng->getDomainRhsFE().reset();
929 pip_mng->getBoundaryLhsFE().reset();
930 pip_mng->getBoundaryRhsFE().reset();
931
933
934 // extract field data for domain parent element
935 auto add_parent_field_domain = [&](auto fe, auto op, auto field, auto bit) {
936 return setParentDofs(
937 fe, field, op, bit,
938
939 [&]() {
940 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
942 return fe_parent;
943 },
944
945 QUIET, Sev::noisy);
946 };
947
948 // extract field data for boundary parent element
949 auto add_parent_field_bdy = [&](auto fe, auto op, auto field, auto bit) {
950 return setParentDofs(
951 fe, field, op, bit,
952
953 [&]() {
954 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
956 return fe_parent;
957 },
958
959 QUIET, Sev::noisy);
960 };
961
962 // run this element on element with given bit, otherwise run other nested
963 // element
964 auto create_run_parent_op = [&](auto parent_fe_ptr, auto this_fe_ptr,
965 auto fe_bit) {
966 auto parent_mask = fe_bit;
967 parent_mask.flip();
968 return new OpRunParent(parent_fe_ptr, BitRefLevel().set(), parent_mask,
969 this_fe_ptr, fe_bit, BitRefLevel().set(), QUIET,
970 Sev::inform);
971 };
972
973 // create hierarchy of nested elements
974 auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr, auto fe_bit) {
975 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
976 for (int l = 0; l < nb_levels; ++l)
977 parents_elems_ptr_vec.emplace_back(
978 boost::make_shared<DomianParentEle>(mField));
979 for (auto l = 1; l < nb_levels; ++l) {
980 parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
981 create_run_parent_op(parents_elems_ptr_vec[l], this_fe_ptr, fe_bit));
982 }
983 return parents_elems_ptr_vec[0];
984 };
985
986 // solve projection problem
987 auto solve_projection = [&](auto exe_test) {
989
990 auto set_domain_rhs = [&](auto fe) {
992 auto &pip = fe->getOpPtrVector();
993
994 auto u_ptr = boost::make_shared<MatrixDouble>();
995 auto p_ptr = boost::make_shared<VectorDouble>();
996 auto h_ptr = boost::make_shared<VectorDouble>();
997 auto g_ptr = boost::make_shared<VectorDouble>();
998
999 auto eval_fe_ptr = boost::make_shared<DomianParentEle>(mField);
1000 {
1001 auto &pip = eval_fe_ptr->getOpPtrVector();
1004 eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
1005
1006 [&]() {
1007 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1008 new DomianParentEle(mField));
1009 return fe_parent;
1010 },
1011
1012 QUIET, Sev::noisy);
1013 // That can be done much smarter, by block, field by field. For
1014 // simplicity is like that.
1015 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "U",
1017 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1018 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "P",
1020 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1021 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "H",
1023 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1024 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "G",
1026 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1027 }
1028 auto parent_eval_fe_ptr =
1029 get_parents_vel_fe_ptr(eval_fe_ptr, bit(get_projection_bit()));
1030 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1032
1033 auto assemble_fe_ptr = boost::make_shared<DomianParentEle>(mField);
1034 {
1035 auto &pip = assemble_fe_ptr->getOpPtrVector();
1038 assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1039
1040 [&]() {
1041 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1042 new DomianParentEle(mField));
1043 return fe_parent;
1044 },
1045
1046 QUIET, Sev::noisy);
1047 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "U",
1049 pip.push_back(new OpDomainAssembleVector("U", u_ptr));
1050 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "P",
1052 pip.push_back(new OpDomainAssembleScalar("P", p_ptr));
1053 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "H",
1055 pip.push_back(new OpDomainAssembleScalar("H", h_ptr));
1056 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "G",
1058 pip.push_back(new OpDomainAssembleScalar("G", g_ptr));
1059 }
1060 auto parent_assemble_fe_ptr =
1061 get_parents_vel_fe_ptr(assemble_fe_ptr, bit(get_current_bit()));
1062 pip.push_back(create_run_parent_op(
1063 parent_assemble_fe_ptr, assemble_fe_ptr, bit(get_current_bit())));
1064
1066 };
1067
1068 auto set_domain_lhs = [&](auto fe) {
1070
1071 auto &pip = fe->getOpPtrVector();
1072
1074
1076 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1077
1078 [&]() {
1079 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1080 new DomianParentEle(mField));
1081 return fe_parent;
1082 },
1083
1084 QUIET, Sev::noisy);
1085
1086 // That can be done much smarter, by block, field by field. For simplicity
1087 // is like that.
1088 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U",
1090 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U",
1092 pip.push_back(new OpDomainMassU("U", "U"));
1093 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P",
1095 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P",
1097 pip.push_back(new OpDomainMassP("P", "P"));
1098 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H",
1100 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H",
1102 pip.push_back(new OpDomainMassH("H", "H"));
1103 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G",
1105 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G",
1107 pip.push_back(new OpDomainMassG("G", "G"));
1108
1110 };
1111
1112 auto set_bdy_rhs = [&](auto fe) {
1114 auto &pip = fe->getOpPtrVector();
1115
1116 auto l_ptr = boost::make_shared<VectorDouble>();
1117
1118 auto eval_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
1119 {
1120 auto &pip = eval_fe_ptr->getOpPtrVector();
1122 eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
1123
1124 [&]() {
1125 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1127 return fe_parent;
1128 },
1129
1130 QUIET, Sev::noisy);
1131 // That can be done much smarter, by block, field by field. For
1132 // simplicity is like that.
1133 CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPROW, "L",
1135 pip.push_back(new OpCalculateScalarFieldValues("L", l_ptr));
1136 }
1137 auto parent_eval_fe_ptr =
1138 get_parents_vel_fe_ptr(eval_fe_ptr, bit(get_projection_bit()));
1139 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1141
1142 auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
1143 {
1144 auto &pip = assemble_fe_ptr->getOpPtrVector();
1146 assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1147
1148 [&]() {
1149 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1151 return fe_parent;
1152 },
1153
1154 QUIET, Sev::noisy);
1155
1156 struct OpLSize : public BoundaryEleOp {
1157 OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
1158 : BoundaryEleOp(NOSPACE, DomainEleOp::OPSPACE), lPtr(l_ptr) {}
1159 MoFEMErrorCode doWork(int, EntityType, EntData &) {
1161 if (lPtr->size() != getGaussPts().size2()) {
1162 lPtr->resize(getGaussPts().size2());
1163 lPtr->clear();
1164 }
1166 }
1167
1168 private:
1169 boost::shared_ptr<VectorDouble> lPtr;
1170 };
1171
1172 pip.push_back(new OpLSize(l_ptr));
1173
1174 CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW, "L",
1176 pip.push_back(new OpBoundaryAssembleScalar("L", l_ptr));
1177 }
1178 auto parent_assemble_fe_ptr =
1179 get_parents_vel_fe_ptr(assemble_fe_ptr, bit(get_current_bit()));
1180 pip.push_back(create_run_parent_op(
1181 parent_assemble_fe_ptr, assemble_fe_ptr, bit(get_current_bit())));
1182
1184 };
1185
1186 auto set_bdy_lhs = [&](auto fe) {
1188
1189 auto &pip = fe->getOpPtrVector();
1190
1192 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1193
1194 [&]() {
1195 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1197 return fe_parent;
1198 },
1199
1200 QUIET, Sev::noisy);
1201
1202 // That can be done much smarter, by block, field by field. For simplicity
1203 // is like that.
1204 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L",
1206 CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "L",
1208 pip.push_back(new OpBoundaryMassL("L", "L"));
1209
1211 };
1212
1213 auto create_subdm = [&]() -> SmartPetscObj<DM> {
1214 auto level_ents_ptr = boost::make_shared<Range>();
1215 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
1216 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
1217
1218 auto get_prj_ents = [&]() {
1219 Range prj_mesh;
1220 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_projection_bit()),
1221 BitRefLevel().set(),
1222 SPACE_DIM, prj_mesh);
1223 auto common_ents = intersect(prj_mesh, *level_ents_ptr);
1224 prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
1225 .subset_by_dimension(SPACE_DIM);
1226
1227 return prj_mesh;
1228 };
1229
1230 auto prj_ents = get_prj_ents();
1231
1232 if (get_global_size(prj_ents.size())) {
1233
1234 auto rebuild = [&]() {
1235 auto prb_mng = mField.getInterface<ProblemsManager>();
1237
1238 std::vector<std::string> fields{"U", "P", "H", "G", "L"};
1239 std::map<std::string, boost::shared_ptr<Range>> range_maps{
1240
1241 {"U", level_ents_ptr},
1242 {"P", level_ents_ptr},
1243 {"H", level_ents_ptr},
1244 {"G", level_ents_ptr},
1245 {"L", level_ents_ptr}
1246
1247 };
1248
1249 CHKERR prb_mng->buildSubProblem("SUB_SOLVER", fields, fields,
1250 simple->getProblemName(), PETSC_TRUE,
1251 &range_maps, &range_maps);
1252
1253 // partition problem
1254 CHKERR prb_mng->partitionFiniteElements("SUB_SOLVER", true, 0,
1256 // set ghost nodes
1257 CHKERR prb_mng->partitionGhostDofsOnDistributedMesh("SUB_SOLVER");
1259 };
1260
1261 MOFEM_LOG("FS", Sev::verbose) << "Create projection problem";
1262
1263 CHKERR rebuild();
1264
1265 auto dm = simple->getDM();
1266 DM subdm;
1267 CHKERR DMCreate(mField.get_comm(), &subdm);
1268 CHKERR DMSetType(subdm, "DMMOFEM");
1269 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
1270 return SmartPetscObj<DM>(subdm);
1271 }
1272
1273 MOFEM_LOG("FS", Sev::inform) << "Nothing to project";
1274
1275 return SmartPetscObj<DM>();
1276 };
1277
1278 auto create_dummy_dm = [&]() {
1279 auto dummy_dm = createDM(mField.get_comm(), "DMMOFEM");
1281 simple->getProblemName().c_str(),
1283 "create dummy dm");
1284 return dummy_dm;
1285 };
1286
1287 auto subdm = create_subdm();
1288 if (subdm) {
1289
1290 pip_mng->getDomainRhsFE().reset();
1291 pip_mng->getDomainLhsFE().reset();
1292 pip_mng->getBoundaryRhsFE().reset();
1293 pip_mng->getBoundaryLhsFE().reset();
1294 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
1295 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
1296 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
1297 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
1298 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1299 pip_mng->getDomainRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
1300 return get_fe_bit(fe_ptr).test(nb_levels - 1);
1301 };
1302 pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
1303 pip_mng->getBoundaryRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
1304 return get_fe_bit(fe_ptr).test(nb_levels - 1);
1305 };
1306
1307 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1308 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1309 CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
1310 CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
1311
1312 auto D = createDMVector(subdm);
1313 auto F = vectorDuplicate(D);
1314
1315 auto ksp = pip_mng->createKSP(subdm);
1316 CHKERR KSPSetFromOptions(ksp);
1317 CHKERR KSPSetUp(ksp);
1318
1319 // get vector norm
1320 auto get_norm = [&](auto x) {
1321 double nrm;
1322 CHKERR VecNorm(x, NORM_2, &nrm);
1323 return nrm;
1324 };
1325
1326 /**
1327 * @brief Zero DOFs, used by FieldBlas
1328 */
1329 auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
1331 for (auto &v : ent_ptr->getEntFieldData()) {
1332 v = 0;
1333 }
1335 };
1336
1337 /**
1338 * @brief cut-off values at nodes, i.e. abs("H") <= 1
1339 *
1340 */
1341 auto cut_off_dofs = [&]() {
1343
1344 Range current_verts;
1345 CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit(get_current_bit()),
1346 BitRefLevel().set(),
1347 MBVERTEX, current_verts);
1348
1349 auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
1351 for (auto &h : ent_ptr->getEntFieldData()) {
1352 h = cut_off(h);
1353 }
1355 };
1356
1357 auto field_blas = mField.getInterface<FieldBlas>();
1358 CHKERR field_blas->fieldLambdaOnEntities(cut_off_verts, "H",
1359 &current_verts);
1361 };
1362
1363 auto solve = [&](auto S) {
1365 CHKERR VecZeroEntries(S);
1366 CHKERR VecZeroEntries(F);
1367 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
1368 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
1369 CHKERR KSPSolve(ksp, F, S);
1370 CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
1371 CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
1373 };
1374
1375 MOFEM_LOG("FS", Sev::inform) << "Solve projection";
1376 CHKERR solve(D);
1377
1378 auto glob_x = createDMVector(simple->getDM());
1379 auto sub_x = createDMVector(subdm);
1380 auto dummy_dm = create_dummy_dm();
1381
1382 /**
1383 * @brief get TSTheta data operators
1384 */
1385 auto apply_restrict = [&]() {
1386 auto get_is = [](auto v) {
1387 IS iy;
1388 auto create = [&]() {
1390 int n, ystart;
1391 CHKERR VecGetLocalSize(v, &n);
1392 CHKERR VecGetOwnershipRange(v, &ystart, NULL);
1393 CHKERR ISCreateStride(PETSC_COMM_SELF, n, ystart, 1, &iy);
1395 };
1396 CHK_THROW_MESSAGE(create(), "create is");
1397 return SmartPetscObj<IS>(iy);
1398 };
1399
1400 auto iy = get_is(glob_x);
1401 auto s = createVecScatter(glob_x, PETSC_NULL, glob_x, iy);
1402
1404 DMSubDomainRestrict(simple->getDM(), s, PETSC_NULL, dummy_dm),
1405 "restrict");
1406 Vec X0, Xdot;
1407 CHK_THROW_MESSAGE(DMGetNamedGlobalVector(dummy_dm, "TSTheta_X0", &X0),
1408 "get X0");
1410 DMGetNamedGlobalVector(dummy_dm, "TSTheta_Xdot", &Xdot),
1411 "get Xdot");
1412
1413 auto forward_ghost = [](auto g) {
1415 CHKERR VecGhostUpdateBegin(g, INSERT_VALUES, SCATTER_FORWARD);
1416 CHKERR VecGhostUpdateEnd(g, INSERT_VALUES, SCATTER_FORWARD);
1418 };
1419
1420 CHK_THROW_MESSAGE(forward_ghost(X0), "");
1421 CHK_THROW_MESSAGE(forward_ghost(Xdot), "");
1422
1423 if constexpr (debug) {
1424 MOFEM_LOG("FS", Sev::inform)
1425 << "Reverse restrict: X0 " << get_norm(X0) << " Xdot "
1426 << get_norm(Xdot);
1427 }
1428
1429 return std::vector<Vec>{X0, Xdot};
1430 };
1431
1432 auto ts_solver_vecs = apply_restrict();
1433
1434 if (ts_solver_vecs.size()) {
1435
1436 for (auto v : ts_solver_vecs) {
1437 MOFEM_LOG("FS", Sev::inform) << "Solve projection vector";
1438
1439 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1440 SCATTER_REVERSE);
1441 CHKERR solve(sub_x);
1442
1443 for (auto f : {"U", "P", "H", "G", "L"}) {
1444 MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1445 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1446 }
1447 CHKERR DMoFEMMeshToLocalVector(subdm, sub_x, INSERT_VALUES,
1448 SCATTER_REVERSE);
1449 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1450 SCATTER_FORWARD);
1451
1452 MOFEM_LOG("FS", Sev::inform) << "Norm V " << get_norm(v);
1453 }
1454
1455 CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_X0",
1456 &ts_solver_vecs[0]);
1457 CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_Xdot",
1458 &ts_solver_vecs[1]);
1459 }
1460
1461 for (auto f : {"U", "P", "H", "G", "L"}) {
1462 MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1463 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1464 }
1465 CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
1466 CHKERR cut_off_dofs();
1467 }
1468
1470 };
1471
1472 // postprocessing (only for debugging proposes)
1473 auto post_proc = [&](auto exe_test) {
1475 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
1476 auto &pip = post_proc_fe->getOpPtrVector();
1477 post_proc_fe->exeTestHook = exe_test;
1478
1480
1482 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1483
1484 [&]() {
1485 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1486 new DomianParentEle(mField));
1488 fe_parent->getOpPtrVector(), {H1});
1489 return fe_parent;
1490 },
1491
1492 QUIET, Sev::noisy);
1493
1495
1496 auto u_ptr = boost::make_shared<MatrixDouble>();
1497 auto p_ptr = boost::make_shared<VectorDouble>();
1498 auto h_ptr = boost::make_shared<VectorDouble>();
1499 auto g_ptr = boost::make_shared<VectorDouble>();
1500
1501 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U",
1503 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1504 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P",
1506 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1507 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H",
1509 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1510 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G",
1512 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1513
1514 post_proc_fe->getOpPtrVector().push_back(
1515
1516 new OpPPMap(post_proc_fe->getPostProcMesh(),
1517 post_proc_fe->getMapGaussPts(),
1518
1519 {{"P", p_ptr}, {"H", h_ptr}, {"G", g_ptr}},
1520
1521 {{"U", u_ptr}},
1522
1523 {},
1524
1525 {}
1526
1527 )
1528
1529 );
1530
1531 auto dm = simple->getDM();
1532 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
1533 CHKERR post_proc_fe->writeFile("out_projection.h5m");
1534
1536 };
1537
1538 CHKERR solve_projection([](FEMethod *fe_ptr) {
1539 return get_fe_bit(fe_ptr).test(get_current_bit());
1540 });
1541
1542 if constexpr (debug) {
1543 CHKERR post_proc([](FEMethod *fe_ptr) {
1544 return get_fe_bit(fe_ptr).test(get_current_bit());
1545 });
1546 }
1547
1548 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
1549 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
1550 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
1551 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
1552
1554}
1555//! [Data projection]
1556
1557//! [Push operators to pip]
1560 auto simple = mField.getInterface<Simple>();
1561
1563
1564 auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
1565 return setParentDofs(
1566 fe, field, op, bit(get_skin_parent_bit()),
1567
1568 [&]() {
1569 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1570 new DomianParentEle(mField));
1571 return fe_parent;
1572 },
1573
1574 QUIET, Sev::noisy);
1575 };
1576
1577 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
1578 return setParentDofs(
1579 fe, field, op, bit(get_skin_parent_bit()),
1580
1581 [&]() {
1582 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1584 return fe_parent;
1585 },
1586
1587 QUIET, Sev::noisy);
1588 };
1589
1590 auto test_bit_child = [](FEMethod *fe_ptr) {
1591 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1592 get_start_bit() + nb_levels - 1);
1593 };
1594
1595 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
1596 auto u_ptr = boost::make_shared<MatrixDouble>();
1597 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
1598 auto dot_h_ptr = boost::make_shared<VectorDouble>();
1599 auto h_ptr = boost::make_shared<VectorDouble>();
1600 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
1601 auto g_ptr = boost::make_shared<VectorDouble>();
1602 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
1603 auto lambda_ptr = boost::make_shared<VectorDouble>();
1604 auto p_ptr = boost::make_shared<VectorDouble>();
1605 auto div_u_ptr = boost::make_shared<VectorDouble>();
1606
1607 // Push element from reference configuration to current configuration in 3d
1608 // space
1609 auto set_domain_general = [&](auto fe) {
1611 auto &pip = fe->getOpPtrVector();
1612
1614
1616 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1617
1618 [&]() {
1619 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1620 new DomianParentEle(mField));
1622 fe_parent->getOpPtrVector(), {H1});
1623 return fe_parent;
1624 },
1625
1626 QUIET, Sev::noisy);
1627
1628 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1629 pip.push_back(
1631 pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
1633 "U", grad_u_ptr));
1634 pip.push_back(
1636 "U", div_u_ptr));
1637
1638 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1639 pip.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
1640 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1641 pip.push_back(
1642 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1643
1644 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1645 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1646 pip.push_back(
1647 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
1648
1649 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1650 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1652 };
1653
1654 auto set_domain_rhs = [&](auto fe) {
1656 auto &pip = fe->getOpPtrVector();
1657
1658 CHKERR set_domain_general(fe);
1659
1660 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1661 pip.push_back(new OpRhsU("U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
1662 grad_h_ptr, g_ptr, p_ptr));
1663
1664 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1665 pip.push_back(new OpRhsH<false>("H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
1666 grad_g_ptr));
1667
1668 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1669 pip.push_back(new OpRhsG<false>("G", h_ptr, grad_h_ptr, g_ptr));
1670
1671 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1672 pip.push_back(new OpDomainAssembleScalar(
1673 "P", div_u_ptr, [](const double r, const double, const double) {
1674 return cylindrical(r);
1675 }));
1676 pip.push_back(new OpDomainAssembleScalar(
1677 "P", p_ptr, [](const double r, const double, const double) {
1678 return eps * cylindrical(r);
1679 }));
1681 };
1682
1683 auto set_domain_lhs = [&](auto fe) {
1685 auto &pip = fe->getOpPtrVector();
1686
1687 CHKERR set_domain_general(fe);
1688
1689 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1690 {
1691 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1692 pip.push_back(new OpLhsU_dU("U", u_ptr, grad_u_ptr, h_ptr));
1693 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1694 pip.push_back(
1695 new OpLhsU_dH("U", "H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
1696
1697 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1698 pip.push_back(new OpLhsU_dG("U", "G", grad_h_ptr));
1699 }
1700
1701 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1702 {
1703 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1704 pip.push_back(new OpLhsH_dU("H", "U", grad_h_ptr));
1705 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1706 pip.push_back(new OpLhsH_dH<false>("H", u_ptr, h_ptr, grad_g_ptr));
1707 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1708 pip.push_back(new OpLhsH_dG<false>("H", "G", h_ptr));
1709 }
1710
1711 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1712 {
1713 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1714 pip.push_back(new OpLhsG_dH<false>("G", "H", h_ptr));
1715 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1716 pip.push_back(new OpLhsG_dG("G"));
1717 }
1718
1719 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1720 {
1721 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1722 pip.push_back(new OpMixScalarTimesDiv(
1723 "P", "U",
1724 [](const double r, const double, const double) {
1725 return cylindrical(r);
1726 },
1727 true, false));
1728 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P");
1729 pip.push_back(new OpDomainMassP("P", "P", [](double r, double, double) {
1730 return eps * cylindrical(r);
1731 }));
1732 }
1733
1735 };
1736
1737 auto get_block_name = [](auto name) {
1738 return boost::format("%s(.*)") % "WETTING_ANGLE";
1739 };
1740
1741 auto get_blocks = [&](auto &&name) {
1742 return mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
1743 std::regex(name.str()));
1744 };
1745
1746 auto set_boundary_rhs = [&](auto fe) {
1748 auto &pip = fe->getOpPtrVector();
1749
1751 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1752
1753 [&]() {
1754 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1756 return fe_parent;
1757 },
1758
1759 QUIET, Sev::noisy);
1760
1761 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
1762 pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
1763
1764 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
1765 pip.push_back(new OpCalculateScalarFieldValues("L", lambda_ptr));
1766 pip.push_back(new OpNormalConstrainRhs("L", u_ptr));
1767
1769 pip, mField, "L", {}, "INFLUX", Sev::inform);
1770
1771 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
1772 pip.push_back(new OpNormalForceRhs("U", lambda_ptr));
1773
1774 // push operators to the side element which is called from op_bdy_side
1775 auto op_bdy_side =
1776 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
1777 op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
1778
1780 op_bdy_side->getOpPtrVector(), {H1});
1781
1783 op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
1785
1786 [&]() {
1787 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1788 new DomianParentEle(mField));
1790 fe_parent->getOpPtrVector(), {H1});
1791 return fe_parent;
1792 },
1793
1794 QUIET, Sev::noisy);
1795
1796 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1797 "H");
1798 op_bdy_side->getOpPtrVector().push_back(
1799 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1800 // push bdy side op
1801 pip.push_back(op_bdy_side);
1802
1803 // push operators for rhs wetting angle
1804 for (auto &b : get_blocks(get_block_name("WETTING_ANGLE"))) {
1805 Range force_edges;
1806 std::vector<double> attr_vec;
1807 CHKERR b->getMeshsetIdEntitiesByDimension(
1808 mField.get_moab(), SPACE_DIM - 1, force_edges, true);
1809 b->getAttributes(attr_vec);
1810 if (attr_vec.size() != 1)
1811 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Should be one attribute");
1812 MOFEM_LOG("FS", Sev::inform) << "Wetting angle: " << attr_vec.front();
1813 // need to find the attributes and pass to operator
1814 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
1815 pip.push_back(new OpWettingAngleRhs(
1816 "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
1817 attr_vec.front()));
1818 }
1819
1821 };
1822
1823 auto set_boundary_lhs = [&](auto fe) {
1825 auto &pip = fe->getOpPtrVector();
1826
1828 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1829
1830 [&]() {
1831 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1833 return fe_parent;
1834 },
1835
1836 QUIET, Sev::noisy);
1837
1838 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
1839 CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "U");
1840 pip.push_back(new OpNormalConstrainLhs("L", "U"));
1841
1842 auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
1843 auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
1844
1845 // push operators to the side element which is called from op_bdy_side
1846 auto op_bdy_side =
1847 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
1848 op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
1849
1851 op_bdy_side->getOpPtrVector(), {H1});
1852
1854 op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
1856
1857 [&]() {
1858 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1859 new DomianParentEle(mField));
1861 fe_parent->getOpPtrVector(), {H1});
1862 return fe_parent;
1863 },
1864
1865 QUIET, Sev::noisy);
1866
1867 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1868 "H");
1869 op_bdy_side->getOpPtrVector().push_back(
1870 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1871 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
1872 "H");
1873 op_bdy_side->getOpPtrVector().push_back(
1874 new OpLoopSideGetDataForSideEle("H", col_ind_ptr, col_diff_base_ptr));
1875
1876 // push bdy side op
1877 pip.push_back(op_bdy_side);
1878
1879 // push operators for lhs wetting angle
1880 for (auto &b : get_blocks(get_block_name("WETTING_ANGLE"))) {
1881 Range force_edges;
1882 std::vector<double> attr_vec;
1883 CHKERR b->getMeshsetIdEntitiesByDimension(
1884 mField.get_moab(), SPACE_DIM - 1, force_edges, true);
1885 b->getAttributes(attr_vec);
1886 if (attr_vec.size() != 1)
1887 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Should be one attribute");
1888 MOFEM_LOG("FS", Sev::inform)
1889 << "wetting angle edges size " << force_edges.size();
1890
1891 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
1892 pip.push_back(new OpWettingAngleLhs(
1893 "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
1894 boost::make_shared<Range>(force_edges), attr_vec.front()));
1895 }
1896
1898 };
1899
1900 auto *pip_mng = mField.getInterface<PipelineManager>();
1901
1902 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1903 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1904 CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
1905 CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
1906
1907 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
1908 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
1909 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
1910 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
1911
1912 pip_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
1913 pip_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
1914 pip_mng->getBoundaryLhsFE()->exeTestHook = test_bit_child;
1915 pip_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
1916
1918}
1919//! [Push operators to pip]
1920
1921/**
1922 * @brief Monitor solution
1923 *
1924 * This functions is cplled by TS kso at the end of each step. It is used
1925 */
1926struct Monitor : public FEMethod {
1928 SmartPetscObj<DM> dm, boost::shared_ptr<moab::Core> post_proc_mesh,
1929 boost::shared_ptr<PostProcEleDomainCont> post_proc,
1930 boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
1931 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
1932 p)
1933 : dM(dm), postProcMesh(post_proc_mesh), postProc(post_proc),
1934 postProcEdge(post_proc_edge), liftFE(p.first), liftVec(p.second) {}
1937
1938 MOFEM_LOG("FS", Sev::verbose) << "Monitor";
1939 constexpr int save_every_nth_step = 1;
1940 if (ts_step % save_every_nth_step == 0) {
1941 MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc";
1942 MoFEM::Interface *m_field_ptr;
1943 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
1944 auto post_proc_begin =
1945 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
1946 postProcMesh);
1947 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
1948 *m_field_ptr, postProcMesh);
1949 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
1951 this->getCacheWeakPtr());
1953 this->getCacheWeakPtr());
1954 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end->getFEMethod());
1955 CHKERR post_proc_end->writeFile(
1956 "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
1957 MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc done";
1958 }
1959
1960 liftVec->resize(SPACE_DIM, false);
1961 liftVec->clear();
1963 MPI_Allreduce(MPI_IN_PLACE, &(*liftVec)[0], SPACE_DIM, MPI_DOUBLE, MPI_SUM,
1964 MPI_COMM_WORLD);
1965 MOFEM_LOG("FS", Sev::inform)
1966 << "Step " << ts_step << " time " << ts_t
1967 << " lift vec x: " << (*liftVec)[0] << " y: " << (*liftVec)[1];
1968
1970 }
1971
1972private:
1974 boost::shared_ptr<moab::Core> postProcMesh;
1975 boost::shared_ptr<PostProcEleDomainCont> postProc;
1976 boost::shared_ptr<PostProcEleBdyCont> postProcEdge;
1977 boost::shared_ptr<BoundaryEle> liftFE;
1978 boost::shared_ptr<VectorDouble> liftVec;
1979};
1980
1981//! [Solve]
1984
1986
1987 auto simple = mField.getInterface<Simple>();
1988 auto pip_mng = mField.getInterface<PipelineManager>();
1989
1990 auto create_solver_dm = [&](auto dm) -> SmartPetscObj<DM> {
1991 DM subdm;
1992
1993 auto setup_subdm = [&](auto dm) {
1995 auto simple = mField.getInterface<Simple>();
1996 auto bit_mng = mField.getInterface<BitRefManager>();
1997 auto dm = simple->getDM();
1998 auto level_ents_ptr = boost::make_shared<Range>();
1999 CHKERR bit_mng->getEntitiesByRefLevel(
2000 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
2001 CHKERR DMCreate(mField.get_comm(), &subdm);
2002 CHKERR DMSetType(subdm, "DMMOFEM");
2003 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
2004 CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
2005 CHKERR DMMoFEMAddElement(subdm, simple->getBoundaryFEName());
2006 CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
2007 CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_FALSE);
2008 for (auto f : {"U", "P", "H", "G", "L"}) {
2009 CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
2010 CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
2011 }
2012 CHKERR DMSetUp(subdm);
2014 };
2015
2016 CHK_THROW_MESSAGE(setup_subdm(dm), "create subdm");
2017
2018 return SmartPetscObj<DM>(subdm);
2019 };
2020
2021 auto get_fe_post_proc = [&](auto post_proc_mesh) {
2022 auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
2023 return setParentDofs(
2024 fe, field, op, bit(get_skin_parent_bit()),
2025
2026 [&]() {
2027 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2028 new DomianParentEle(mField));
2029 return fe_parent;
2030 },
2031
2032 QUIET, Sev::noisy);
2033 };
2034
2035 auto post_proc_fe =
2036 boost::make_shared<PostProcEleDomainCont>(mField, post_proc_mesh);
2037 post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
2038 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2039 get_start_bit() + nb_levels - 1);
2040 };
2041
2042 auto u_ptr = boost::make_shared<MatrixDouble>();
2043 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2044 auto h_ptr = boost::make_shared<VectorDouble>();
2045 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2046 auto p_ptr = boost::make_shared<VectorDouble>();
2047 auto g_ptr = boost::make_shared<VectorDouble>();
2048 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2049
2051 post_proc_fe->getOpPtrVector(), {H1});
2052
2054 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2055
2056 [&]() {
2057 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2058 new DomianParentEle(mField));
2060 fe_parent->getOpPtrVector(), {H1});
2061 return fe_parent;
2062 },
2063
2064 QUIET, Sev::noisy);
2065
2066 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U");
2067 post_proc_fe->getOpPtrVector().push_back(
2069 post_proc_fe->getOpPtrVector().push_back(
2071 grad_u_ptr));
2072
2073 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H");
2074 post_proc_fe->getOpPtrVector().push_back(
2075 new OpCalculateScalarFieldValues("H", h_ptr));
2076 post_proc_fe->getOpPtrVector().push_back(
2077 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
2078
2079 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P");
2080 post_proc_fe->getOpPtrVector().push_back(
2081 new OpCalculateScalarFieldValues("P", p_ptr));
2082
2083 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G");
2084 post_proc_fe->getOpPtrVector().push_back(
2085 new OpCalculateScalarFieldValues("G", g_ptr));
2086 post_proc_fe->getOpPtrVector().push_back(
2087 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
2088
2090
2091 post_proc_fe->getOpPtrVector().push_back(
2092
2093 new OpPPMap(
2094 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2095
2096 {{"H", h_ptr}, {"P", p_ptr}, {"G", g_ptr}},
2097
2098 {{"U", u_ptr}, {"H_GRAD", grad_h_ptr}, {"G_GRAD", grad_g_ptr}},
2099
2100 {{"GRAD_U", grad_u_ptr}},
2101
2102 {}
2103
2104 )
2105
2106 );
2107
2108 return post_proc_fe;
2109 };
2110
2111 auto get_bdy_post_proc_fe = [&](auto post_proc_mesh) {
2112 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2113 return setParentDofs(
2114 fe, field, op, bit(get_skin_parent_bit()),
2115
2116 [&]() {
2117 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2118 new BoundaryParentEle(mField));
2119 return fe_parent;
2120 },
2121
2122 QUIET, Sev::noisy);
2123 };
2124
2125 auto post_proc_fe =
2126 boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
2127 post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
2128 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2129 get_start_bit() + nb_levels - 1);
2130 };
2131
2132 CHKERR setParentDofs(
2133 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2134
2135 [&]() {
2136 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2137 new BoundaryParentEle(mField));
2138 return fe_parent;
2139 },
2140
2141 QUIET, Sev::noisy);
2142
2143 auto u_ptr = boost::make_shared<MatrixDouble>();
2144 auto p_ptr = boost::make_shared<VectorDouble>();
2145 auto lambda_ptr = boost::make_shared<VectorDouble>();
2146
2147 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "U");
2148 post_proc_fe->getOpPtrVector().push_back(
2150
2151 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "L");
2152 post_proc_fe->getOpPtrVector().push_back(
2153 new OpCalculateScalarFieldValues("L", lambda_ptr));
2154
2155 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "P");
2156 post_proc_fe->getOpPtrVector().push_back(
2157 new OpCalculateScalarFieldValues("P", p_ptr));
2158
2160
2161 post_proc_fe->getOpPtrVector().push_back(
2162
2163 new OpPPMap(post_proc_fe->getPostProcMesh(),
2164 post_proc_fe->getMapGaussPts(),
2165
2166 OpPPMap::DataMapVec{{"L", lambda_ptr}, {"P", p_ptr}},
2167
2168 OpPPMap::DataMapMat{{"U", u_ptr}},
2169
2171
2173
2174 )
2175
2176 );
2177
2178 return post_proc_fe;
2179 };
2180
2181 auto get_lift_fe = [&]() {
2182 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2183 return setParentDofs(
2184 fe, field, op, bit(get_skin_parent_bit()),
2185
2186 [&]() {
2187 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2188 new BoundaryParentEle(mField));
2189 return fe_parent;
2190 },
2191
2192 QUIET, Sev::noisy);
2193 };
2194
2195 auto fe = boost::make_shared<BoundaryEle>(mField);
2196 fe->exeTestHook = [](FEMethod *fe_ptr) {
2197 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2198 get_start_bit() + nb_levels - 1);
2199 };
2200
2201 auto lift_ptr = boost::make_shared<VectorDouble>();
2202 auto p_ptr = boost::make_shared<VectorDouble>();
2203 auto ents_ptr = boost::make_shared<Range>();
2204
2205 CHKERR setParentDofs(
2206 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2207
2208 [&]() {
2209 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2210 new BoundaryParentEle(mField));
2211 return fe_parent;
2212 },
2213
2214 QUIET, Sev::noisy);
2215
2216 std::vector<const CubitMeshSets *> vec_ptr;
2217 CHKERR mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
2218 std::regex("LIFT"), vec_ptr);
2219 for (auto m_ptr : vec_ptr) {
2220 auto meshset = m_ptr->getMeshset();
2221 Range ents;
2222 CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
2223 ents, true);
2224 ents_ptr->merge(ents);
2225 }
2226
2227 MOFEM_LOG("FS", Sev::noisy) << "Lift ents " << (*ents_ptr);
2228
2229 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "P");
2230 fe->getOpPtrVector().push_back(
2231 new OpCalculateScalarFieldValues("P", p_ptr));
2232 fe->getOpPtrVector().push_back(
2233 new OpCalculateLift("P", p_ptr, lift_ptr, ents_ptr));
2234
2235 return std::make_pair(fe, lift_ptr);
2236 };
2237
2238 auto set_post_proc_monitor = [&](auto ts) {
2240 DM dm;
2241 CHKERR TSGetDM(ts, &dm);
2242 boost::shared_ptr<FEMethod> null_fe;
2243 auto post_proc_mesh = boost::make_shared<moab::Core>();
2244 auto monitor_ptr = boost::make_shared<Monitor>(
2245 SmartPetscObj<DM>(dm, true), post_proc_mesh,
2246 get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
2247 get_lift_fe());
2248 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
2249 null_fe, monitor_ptr);
2251 };
2252
2253 auto dm = simple->getDM();
2254 auto ts = createTS(mField.get_comm());
2255 CHKERR TSSetDM(ts, dm);
2256
2257 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
2258 tsPrePostProc = ts_pre_post_proc;
2259
2260 if (auto ptr = tsPrePostProc.lock()) {
2261
2262 ptr->fsRawPtr = this;
2263 ptr->solverSubDM = create_solver_dm(simple->getDM());
2264 ptr->globSol = createDMVector(dm);
2265 CHKERR DMoFEMMeshToLocalVector(dm, ptr->globSol, INSERT_VALUES,
2266 SCATTER_FORWARD);
2267 CHKERR VecAssemblyBegin(ptr->globSol);
2268 CHKERR VecAssemblyEnd(ptr->globSol);
2269
2270 auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
2271
2272 CHKERR set_post_proc_monitor(sub_ts);
2273
2274 // Add monitor to time solver
2275 double ftime = 1;
2276 // CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
2277 // CHKERR TSSetSolution(ts, ptr->globSol);
2278 CHKERR TSSetFromOptions(ts);
2279 CHKERR ptr->tsSetUp(ts);
2280 CHKERR TSSetUp(ts);
2281
2282 auto print_fields_in_section = [&]() {
2284 auto section = mField.getInterface<ISManager>()->sectionCreate(
2285 simple->getProblemName());
2286 PetscInt num_fields;
2287 CHKERR PetscSectionGetNumFields(section, &num_fields);
2288 for (int f = 0; f < num_fields; ++f) {
2289 const char *field_name;
2290 CHKERR PetscSectionGetFieldName(section, f, &field_name);
2291 MOFEM_LOG("FS", Sev::inform)
2292 << "Field " << f << " " << std::string(field_name);
2293 }
2295 };
2296
2297 CHKERR print_fields_in_section();
2298
2299 CHKERR TSSolve(ts, ptr->globSol);
2300 }
2301
2303}
2304
2305int main(int argc, char *argv[]) {
2306
2307#ifdef PYTHON_INIT_SURFACE
2308 Py_Initialize();
2309#endif
2310
2311 // Initialisation of MoFEM/PETSc and MOAB data structures
2312 const char param_file[] = "param_file.petsc";
2313 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
2314
2315 // Add logging channel for example
2316 auto core_log = logging::core::get();
2317 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FS"));
2318 LogManager::setLog("FS");
2319 MOFEM_LOG_TAG("FS", "free surface");
2320
2321 try {
2322
2323 //! [Register MoFEM discrete manager in PETSc]
2324 DMType dm_name = "DMMOFEM";
2325 CHKERR DMRegister_MoFEM(dm_name);
2326 //! [Register MoFEM discrete manager in PETSc
2327
2328 //! [Create MoAB]
2329 moab::Core mb_instance; ///< mesh database
2330 moab::Interface &moab = mb_instance; ///< mesh database interface
2331 //! [Create MoAB]
2332
2333 //! [Create MoFEM]
2334 MoFEM::Core core(moab); ///< finite element database
2335 MoFEM::Interface &m_field = core; ///< finite element database interface
2336 //! [Create MoFEM]
2337
2338 //! [FreeSurface]
2339 FreeSurface ex(m_field);
2340 CHKERR ex.runProblem();
2341 //! [FreeSurface]
2342 }
2344
2346
2347#ifdef PYTHON_INIT_SURFACE
2348 if (Py_FinalizeEx() < 0) {
2349 exit(120);
2350 }
2351#endif
2352}
2353
2354std::vector<Range>
2356
2357 auto &moab = mField.get_moab();
2358 auto bit_mng = mField.getInterface<BitRefManager>();
2359 auto comm_mng = mField.getInterface<CommInterface>();
2360
2361 Range vertices;
2362 CHK_THROW_MESSAGE(bit_mng->getEntitiesByTypeAndRefLevel(
2363 bit(0), BitRefLevel().set(), MBVERTEX, vertices),
2364 "can not get vertices on bit 0");
2365
2366 auto &dofs_mi = mField.get_dofs()->get<Unique_mi_tag>();
2367 auto field_bit_number = mField.get_field_bit_number("H");
2368
2369 Range plus_range, minus_range;
2370 std::vector<EntityHandle> plus, minus;
2371
2372 // get vertices on level 0 on plus and minus side
2373 for (auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
2374
2375 const auto f = p->first;
2376 const auto s = p->second;
2377
2378 // Lowest Dof UId for given field (field bit number) on entity f
2379 const auto lo_uid = DofEntity::getLoFieldEntityUId(field_bit_number, f);
2380 const auto hi_uid = DofEntity::getHiFieldEntityUId(field_bit_number, s);
2381 auto it = dofs_mi.lower_bound(lo_uid);
2382 const auto hi_it = dofs_mi.upper_bound(hi_uid);
2383
2384 plus.clear();
2385 minus.clear();
2386 plus.reserve(std::distance(it, hi_it));
2387 minus.reserve(std::distance(it, hi_it));
2388
2389 for (; it != hi_it; ++it) {
2390 const auto v = (*it)->getFieldData();
2391 if (v > 0)
2392 plus.push_back((*it)->getEnt());
2393 else
2394 minus.push_back((*it)->getEnt());
2395 }
2396
2397 plus_range.insert_list(plus.begin(), plus.end());
2398 minus_range.insert_list(minus.begin(), minus.end());
2399 }
2400
2401 MOFEM_LOG_CHANNEL("SYNC");
2402 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
2403 << "Plus range " << plus_range << endl;
2404 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
2405 << "Minus range " << minus_range << endl;
2407
2408 auto get_elems = [&](auto &ents, auto bit, auto mask) {
2409 Range adj;
2410 CHK_MOAB_THROW(moab.get_adjacencies(ents, SPACE_DIM, false, adj,
2411 moab::Interface::UNION),
2412 "can not get adjacencies");
2413 CHK_THROW_MESSAGE(bit_mng->filterEntitiesByRefLevel(bit, mask, adj),
2414 "can not filter elements with bit 0");
2415 return adj;
2416 };
2417
2418 CHKERR comm_mng->synchroniseEntities(plus_range);
2419 CHKERR comm_mng->synchroniseEntities(minus_range);
2420
2421 std::vector<Range> ele_plus(nb_levels), ele_minus(nb_levels);
2422 ele_plus[0] = get_elems(plus_range, bit(0), BitRefLevel().set());
2423 ele_minus[0] = get_elems(minus_range, bit(0), BitRefLevel().set());
2424 auto common = intersect(ele_plus[0], ele_minus[0]);
2425 ele_plus[0] = subtract(ele_plus[0], common);
2426 ele_minus[0] = subtract(ele_minus[0], common);
2427
2428 auto get_children = [&](auto &p, auto &c) {
2430 CHKERR bit_mng->updateRangeByChildren(p, c);
2431 c = c.subset_by_dimension(SPACE_DIM);
2433 };
2434
2435 for (auto l = 1; l != nb_levels; ++l) {
2436 CHK_THROW_MESSAGE(get_children(ele_plus[l - 1], ele_plus[l]),
2437 "get children");
2438 CHK_THROW_MESSAGE(get_children(ele_minus[l - 1], ele_minus[l]),
2439 "get children");
2440 }
2441
2442 auto get_level = [&](auto &p, auto &m, auto z, auto bit, auto mask) {
2443 Range l;
2445 bit_mng->getEntitiesByDimAndRefLevel(bit, mask, SPACE_DIM, l),
2446 "can not get vertices on bit");
2447 l = subtract(l, p);
2448 l = subtract(l, m);
2449 for (auto f = 0; f != z; ++f) {
2450 Range conn;
2451 CHK_MOAB_THROW(moab.get_connectivity(l, conn, true), "");
2452 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(conn);
2453 l = get_elems(conn, bit, mask);
2454 }
2455 return l;
2456 };
2457
2458 std::vector<Range> vec_levels(nb_levels);
2459 for (auto l = nb_levels - 1; l >= 0; --l) {
2460 vec_levels[l] = get_level(ele_plus[l], ele_minus[l], 2 * overlap, bit(l),
2461 BitRefLevel().set());
2462 }
2463
2464 if constexpr (debug) {
2465 if (mField.get_comm_size() == 1) {
2466 for (auto l = 0; l != nb_levels; ++l) {
2467 std::string name = (boost::format("out_r%d.h5m") % l).str();
2468 CHK_THROW_MESSAGE(save_range(mField.get_moab(), name, vec_levels[l]),
2469 "save mesh");
2470 }
2471 }
2472 }
2473
2474 return vec_levels;
2475}
2476
2479
2480 auto simple = mField.getInterface<Simple>();
2481 auto bit_mng = mField.getInterface<BitRefManager>();
2482 auto prb_mng = mField.getInterface<ProblemsManager>();
2483
2484 BitRefLevel start_mask;
2485 for (auto s = 0; s != get_start_bit(); ++s)
2486 start_mask[s] = true;
2487
2488 // store prev_level
2489 Range prev_level;
2490 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_current_bit()),
2491 BitRefLevel().set(), prev_level);
2492 Range prev_level_skin;
2493 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_parent_bit()),
2494 BitRefLevel().set(), prev_level_skin);
2495 // reset bit ref levels
2496 CHKERR bit_mng->lambdaBitRefLevel(
2497 [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
2498 CHKERR bit_mng->setNthBitRefLevel(prev_level, get_projection_bit(), true);
2499 CHKERR bit_mng->setNthBitRefLevel(prev_level_skin, get_skin_projection_bit(),
2500 true);
2501
2502 // set refinement levels
2503 auto set_levels = [&](auto &&
2504 vec_levels /*entities are refined on each level*/) {
2506
2507 // start with zero level, which is the coarsest mesh
2508 Range level0;
2509 CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
2510 CHKERR bit_mng->setNthBitRefLevel(level0, get_start_bit(), true);
2511
2512 // get lower dimension entities
2513 auto get_adj = [&](auto ents) {
2514 Range conn;
2515 CHK_MOAB_THROW(mField.get_moab().get_connectivity(ents, conn, true),
2516 "get conn");
2517 for (auto d = 1; d != SPACE_DIM; ++d) {
2518 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
2519 ents.subset_by_dimension(SPACE_DIM), d, false, conn,
2520 moab::Interface::UNION),
2521 "get adj");
2522 }
2523 ents.merge(conn);
2524 return ents;
2525 };
2526
2527 // set bit levels
2528 for (auto l = 1; l != nb_levels; ++l) {
2529 Range level_prev;
2530 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_start_bit() + l - 1),
2531 BitRefLevel().set(),
2532 SPACE_DIM, level_prev);
2533 Range parents;
2534 CHKERR bit_mng->updateRangeByParent(vec_levels[l], parents);
2535 // subtract entities from previous level, which are refined, so should be
2536 // not there
2537 level_prev = subtract(level_prev, parents);
2538 // and instead add their children
2539 level_prev.merge(vec_levels[l]);
2540 // set bit to each level
2541 CHKERR bit_mng->setNthBitRefLevel(level_prev, get_start_bit() + l, true);
2542 }
2543
2544 // set bit levels to lower dimension entities
2545 for (auto l = 1; l != nb_levels; ++l) {
2546 Range level;
2547 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2548 bit(get_start_bit() + l), BitRefLevel().set(), SPACE_DIM, level);
2549 level = get_adj(level);
2550 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(level);
2551 CHKERR bit_mng->setNthBitRefLevel(level, get_start_bit() + l, true);
2552 }
2553
2555 };
2556
2557 // resolve skin between refined levels
2558 auto set_skins = [&]() {
2560
2561 moab::Skinner skinner(&mField.get_moab());
2562 ParallelComm *pcomm =
2563 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2564
2565 // get skin of bit level
2566 auto get_bit_skin = [&](BitRefLevel bit, BitRefLevel mask) {
2567 Range bit_ents;
2570 bit, mask, SPACE_DIM, bit_ents),
2571 "can't get bit level");
2572 Range bit_skin;
2573 CHK_MOAB_THROW(skinner.find_skin(0, bit_ents, false, bit_skin),
2574 "can't get skin");
2575 return bit_skin;
2576 };
2577
2578 auto get_level_skin = [&]() {
2579 Range skin;
2580 BitRefLevel bit_prev;
2581 for (auto l = 1; l != nb_levels; ++l) {
2582 auto skin_level_mesh = get_bit_skin(bit(l), BitRefLevel().set());
2583 // filter (remove) all entities which are on partitions boarders
2584 CHKERR pcomm->filter_pstatus(skin_level_mesh,
2585 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2586 PSTATUS_NOT, -1, nullptr);
2587 auto skin_level =
2588 get_bit_skin(bit(get_start_bit() + l), BitRefLevel().set());
2589 skin_level = subtract(skin_level,
2590 skin_level_mesh); // get only internal skins, not
2591 // on the body boundary
2592 // get lower dimension adjacencies (FIXME: add edges if 3D)
2593 Range skin_level_verts;
2594 CHKERR mField.get_moab().get_connectivity(skin_level, skin_level_verts,
2595 true);
2596 skin_level.merge(skin_level_verts);
2597
2598 // remove previous level
2599 bit_prev.set(l - 1);
2600 Range level_prev;
2601 CHKERR bit_mng->getEntitiesByRefLevel(bit_prev, BitRefLevel().set(),
2602 level_prev);
2603 skin.merge(subtract(skin_level, level_prev));
2604 }
2605
2606 return skin;
2607 };
2608
2609 auto resolve_shared = [&](auto &&skin) {
2610 Range tmp_skin = skin;
2611
2612 map<int, Range> map_procs;
2613 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2614 tmp_skin, &map_procs);
2615
2616 Range from_other_procs; // entities which also exist on other processors
2617 for (auto &m : map_procs) {
2618 if (m.first != mField.get_comm_rank()) {
2619 from_other_procs.merge(m.second);
2620 }
2621 }
2622
2623 auto common = intersect(
2624 skin, from_other_procs); // entities which are on internal skin
2625 skin.merge(from_other_procs);
2626
2627 // entities which are on processors boards, and several processors are not
2628 // true skin.
2629 if (!common.empty()) {
2630 // skin is internal exist on other procs
2631 skin = subtract(skin, common);
2632 }
2633
2634 return skin;
2635 };
2636
2637 // get parents of entities
2638 auto get_parent_level_skin = [&](auto skin) {
2639 Range skin_parents;
2640 CHKERR bit_mng->updateRangeByParent(
2641 skin.subset_by_dimension(SPACE_DIM - 1), skin_parents);
2642 Range skin_parent_verts;
2643 CHKERR mField.get_moab().get_connectivity(skin_parents, skin_parent_verts,
2644 true);
2645 skin_parents.merge(skin_parent_verts);
2646 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2647 skin_parents);
2648 return skin_parents;
2649 };
2650
2651 auto child_skin = resolve_shared(get_level_skin());
2652 auto parent_skin = get_parent_level_skin(child_skin);
2653
2654 child_skin = subtract(child_skin, parent_skin);
2655 CHKERR bit_mng->setNthBitRefLevel(child_skin, get_skin_child_bit(), true);
2656 CHKERR bit_mng->setNthBitRefLevel(parent_skin, get_skin_parent_bit(), true);
2657
2659 };
2660
2661 // take last level, remove childs on boarder, and set bit
2662 auto set_current = [&]() {
2664 Range last_level;
2665 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_start_bit() + nb_levels - 1),
2666 BitRefLevel().set(), last_level);
2667 Range skin_child;
2668 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_child_bit()),
2669 BitRefLevel().set(), skin_child);
2670
2671 last_level = subtract(last_level, skin_child);
2672 CHKERR bit_mng->setNthBitRefLevel(last_level, get_current_bit(), true);
2674 };
2675
2676 // set bits to levels
2677 CHKERR set_levels(findEntitiesCrossedByPhaseInterface(overlap));
2678 // set bits to skin
2679 CHKERR set_skins();
2680 // set current level bit
2681 CHKERR set_current();
2682
2683 if constexpr (debug) {
2684 if (mField.get_comm_size() == 1) {
2685 for (auto l = 0; l != nb_levels; ++l) {
2686 std::string name = (boost::format("out_level%d.h5m") % l).str();
2687 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_start_bit() + l),
2688 BitRefLevel().set(), name.c_str(), "MOAB",
2689 "PARALLEL=WRITE_PART");
2690 }
2691 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_current_bit()),
2692 BitRefLevel().set(), "current_bit.h5m",
2693 "MOAB", "PARALLEL=WRITE_PART");
2694 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_projection_bit()),
2695 BitRefLevel().set(), "projection_bit.h5m",
2696 "MOAB", "PARALLEL=WRITE_PART");
2697
2698 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_child_bit()),
2699 BitRefLevel().set(), "skin_child_bit.h5m",
2700 "MOAB", "PARALLEL=WRITE_PART");
2701 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_parent_bit()),
2702 BitRefLevel().set(), "skin_parent_bit.h5m",
2703 "MOAB", "PARALLEL=WRITE_PART");
2704 }
2705 }
2706
2708};
2709
2711 boost::shared_ptr<FEMethod> fe_top, std::string field_name,
2712 ForcesAndSourcesCore::UserDataOperator::OpType op,
2713 BitRefLevel child_ent_bit,
2714 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
2715 int verbosity, LogManager::SeverityLevel sev) {
2717
2718 /**
2719 * @brief Collect data from parent elements to child
2720 */
2721 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
2722 add_parent_level =
2723 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
2724 // Evaluate if not last parent element
2725 if (level > 0) {
2726
2727 // Create domain parent FE
2728 auto fe_ptr_current = get_elem();
2729
2730 // Call next level
2731 add_parent_level(
2732 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
2733 fe_ptr_current),
2734 level - 1);
2735
2736 // Add data to curent fe level
2737 if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
2738
2739 // Only base
2740 parent_fe_pt->getOpPtrVector().push_back(
2741
2743
2744 H1, op, fe_ptr_current,
2745
2746 BitRefLevel().set(), BitRefLevel().set(),
2747
2748 child_ent_bit, BitRefLevel().set(),
2749
2750 verbosity, sev));
2751
2752 } else {
2753
2754 // Filed data
2755 parent_fe_pt->getOpPtrVector().push_back(
2756
2758
2759 field_name, op, fe_ptr_current,
2760
2761 BitRefLevel().set(), BitRefLevel().set(),
2762
2763 child_ent_bit, BitRefLevel().set(),
2764
2765 verbosity, sev));
2766 }
2767 }
2768 };
2769
2770 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
2771 nb_levels);
2772
2774}
2775
2778
2779 if (auto ptr = tsPrePostProc.lock()) {
2780 MOFEM_LOG("FS", Sev::inform) << "Run step pre proc";
2781
2782 auto &m_field = ptr->fsRawPtr->mField;
2783 auto simple = m_field.getInterface<Simple>();
2784 auto bit_mng = m_field.getInterface<BitRefManager>();
2785 auto bc_mng = m_field.getInterface<BcManager>();
2786 auto field_blas = m_field.getInterface<FieldBlas>();
2787 auto opt = m_field.getInterface<OperatorsTester>();
2788 auto pip_mng = m_field.getInterface<PipelineManager>();
2789 auto prb_mng = m_field.getInterface<ProblemsManager>();
2790
2791 // get vector norm
2792 auto get_norm = [&](auto x) {
2793 double nrm;
2794 CHKERR VecNorm(x, NORM_2, &nrm);
2795 return nrm;
2796 };
2797
2798 // refine problem and project data, including theta data
2799 auto refine_problem = [&]() {
2801 MOFEM_LOG("FS", Sev::inform) << "Refine problem";
2802 CHKERR ptr->fsRawPtr->refineMesh(refine_overlap);
2803 CHKERR ptr->fsRawPtr->projectData();
2805 };
2806
2807 // set new jacobin operator, since problem and thus tangent matrix size has
2808 // changed
2809 auto set_jacobian_operators = [&]() {
2811 ptr->subB = createDMMatrix(ptr->solverSubDM);
2812 CHKERR KSPReset(ptr->subKSP);
2814 };
2815
2816 // set new solution
2817 auto set_solution = [&]() {
2819 MOFEM_LOG("FS", Sev::inform) << "Set solution";
2820
2821 PetscObjectState state;
2822
2823 // Record the state, and set it again. This is to fool PETSc that solution
2824 // vector is not updated. Otherwise PETSc will treat every step as a first
2825 // step.
2826
2827 // globSol is updated as result mesh refinement - this is not really set
2828 // a new solution.
2829
2830 CHKERR PetscObjectStateGet(getPetscObject(ptr->globSol.get()), &state);
2831 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), ptr->globSol,
2832 INSERT_VALUES, SCATTER_FORWARD);
2833 CHKERR PetscObjectStateSet(getPetscObject(ptr->globSol.get()), state);
2834 MOFEM_LOG("FS", Sev::verbose)
2835 << "Set solution, vector norm " << get_norm(ptr->globSol);
2837 };
2838
2839 PetscBool is_theta;
2840 PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
2841 if (is_theta) {
2842
2843 CHKERR refine_problem(); // refine problem
2844 CHKERR set_jacobian_operators(); // set new jacobian
2845 CHKERR set_solution(); // set solution
2846
2847 } else {
2848 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
2849 "Sorry, only TSTheta handling is implemented");
2850 }
2851
2852 // Need barriers, somme functions in TS solver need are called collectively
2853 // and requite the same state of variables
2854 PetscBarrier((PetscObject)ts);
2855
2856 MOFEM_LOG_CHANNEL("SYNC");
2857 MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PreProc done";
2858 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
2859 }
2860
2862}
2863
2865 if (auto ptr = tsPrePostProc.lock()) {
2866 auto &m_field = ptr->fsRawPtr->mField;
2867 MOFEM_LOG_CHANNEL("SYNC");
2868 MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PostProc done";
2869 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
2870 }
2871 return 0;
2872}
2873
2874MoFEMErrorCode TSPrePostProc::tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t,
2875 Vec f, void *ctx) {
2877 if (auto ptr = tsPrePostProc.lock()) {
2878 auto sub_u = ptr->getSubVector();
2879 auto sub_u_t = vectorDuplicate(sub_u);
2880 auto sub_f = vectorDuplicate(sub_u);
2881 auto scatter = ptr->getScatter(sub_u, u, R);
2882
2883 auto apply_scatter_and_update = [&](auto x, auto sub_x) {
2885 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
2886 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
2887 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
2888 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
2890 };
2891
2892 CHKERR apply_scatter_and_update(u, sub_u);
2893 CHKERR apply_scatter_and_update(u_t, sub_u_t);
2894
2895 CHKERR TsSetIFunction(ts, t, sub_u, sub_u_t, sub_f, ptr->tsCtxPtr.get());
2896 CHKERR VecScatterBegin(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
2897 CHKERR VecScatterEnd(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
2898 }
2900}
2901
2902MoFEMErrorCode TSPrePostProc::tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
2903 PetscReal a, Mat A, Mat B,
2904 void *ctx) {
2906 if (auto ptr = tsPrePostProc.lock()) {
2907 auto sub_u = ptr->getSubVector();
2908 auto sub_u_t = vectorDuplicate(sub_u);
2909 auto scatter = ptr->getScatter(sub_u, u, R);
2910
2911 auto apply_scatter_and_update = [&](auto x, auto sub_x) {
2913 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
2914 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
2915 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
2916 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
2918 };
2919
2920 CHKERR apply_scatter_and_update(u, sub_u);
2921 CHKERR apply_scatter_and_update(u_t, sub_u_t);
2922
2923 CHKERR TsSetIJacobian(ts, t, sub_u, sub_u_t, a, ptr->subB, ptr->subB,
2924 ptr->tsCtxPtr.get());
2925 }
2927}
2928
2929MoFEMErrorCode TSPrePostProc::tsMonitor(TS ts, PetscInt step, PetscReal t,
2930 Vec u, void *ctx) {
2932 if (auto ptr = tsPrePostProc.lock()) {
2933 auto get_norm = [&](auto x) {
2934 double nrm;
2935 CHKERR VecNorm(x, NORM_2, &nrm);
2936 return nrm;
2937 };
2938
2939 auto sub_u = ptr->getSubVector();
2940 auto scatter = ptr->getScatter(sub_u, u, R);
2941 CHKERR VecScatterBegin(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
2942 CHKERR VecScatterEnd(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
2943 CHKERR VecGhostUpdateBegin(sub_u, INSERT_VALUES, SCATTER_FORWARD);
2944 CHKERR VecGhostUpdateEnd(sub_u, INSERT_VALUES, SCATTER_FORWARD);
2945
2946 MOFEM_LOG("FS", Sev::verbose)
2947 << "u norm " << get_norm(u) << " u sub nom " << get_norm(sub_u);
2948
2949 CHKERR TsMonitorSet(ts, step, t, sub_u, ptr->tsCtxPtr.get());
2950 }
2952}
2953
2956 if (auto ptr = tsPrePostProc.lock()) {
2957 MOFEM_LOG("FS", Sev::verbose) << "SetUP sub PC";
2958 ptr->subKSP = createKSP(ptr->fsRawPtr->mField.get_comm());
2959 CHKERR KSPSetFromOptions(ptr->subKSP);
2960 }
2962};
2963
2964MoFEMErrorCode TSPrePostProc::pcApply(PC pc, Vec pc_f, Vec pc_x) {
2966 if (auto ptr = tsPrePostProc.lock()) {
2967 auto sub_x = ptr->getSubVector();
2968 auto sub_f = vectorDuplicate(sub_x);
2969 auto scatter = ptr->getScatter(sub_x, pc_x, R);
2970 CHKERR VecScatterBegin(scatter, pc_f, sub_f, INSERT_VALUES,
2971 SCATTER_REVERSE);
2972 CHKERR VecScatterEnd(scatter, pc_f, sub_f, INSERT_VALUES, SCATTER_REVERSE);
2973 CHKERR KSPSetOperators(ptr->subKSP, ptr->subB, ptr->subB);
2974 MOFEM_LOG("FS", Sev::verbose) << "PCShell solve";
2975 CHKERR KSPSolve(ptr->subKSP, sub_f, sub_x);
2976 MOFEM_LOG("FS", Sev::verbose) << "PCShell solve <- done";
2977 CHKERR VecScatterBegin(scatter, sub_x, pc_x, INSERT_VALUES,
2978 SCATTER_FORWARD);
2979 CHKERR VecScatterEnd(scatter, sub_x, pc_x, INSERT_VALUES, SCATTER_FORWARD);
2980 }
2982};
2983
2985 if (auto ptr = tsPrePostProc.lock()) {
2986 auto prb_ptr = ptr->fsRawPtr->mField.get_problem("SUB_SOLVER");
2987 if (auto sub_data = prb_ptr->getSubData()) {
2988 auto is = sub_data->getSmartColIs();
2989 VecScatter s;
2990 if (fr == R) {
2991 CHK_THROW_MESSAGE(VecScatterCreate(x, PETSC_NULL, y, is, &s),
2992 "crate scatter");
2993 } else {
2994 CHK_THROW_MESSAGE(VecScatterCreate(x, is, y, PETSC_NULL, &s),
2995 "crate scatter");
2996 }
2997 return SmartPetscObj<VecScatter>(s);
2998 }
2999 }
3002}
3003
3006}
3007
3009
3010 auto &m_field = fsRawPtr->mField;
3011 auto simple = m_field.getInterface<Simple>();
3012
3014
3015 auto dm = simple->getDM();
3016
3018 CHKERR TSSetIFunction(ts, globRes, tsSetIFunction, nullptr);
3019 CHKERR TSSetIJacobian(ts, PETSC_NULL, PETSC_NULL, tsSetIJacobian, nullptr);
3020 CHKERR TSMonitorSet(ts, tsMonitor, fsRawPtr, PETSC_NULL);
3021
3022 SNES snes;
3023 CHKERR TSGetSNES(ts, &snes);
3024 auto snes_ctx_ptr = getDMSnesCtx(dm);
3025
3026 auto set_section_monitor = [&](auto snes) {
3028 CHKERR SNESMonitorSet(snes,
3029 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
3030 void *))MoFEMSNESMonitorFields,
3031 (void *)(snes_ctx_ptr.get()), nullptr);
3033 };
3034
3035 CHKERR set_section_monitor(snes);
3036
3037 auto ksp = createKSP(m_field.get_comm());
3038 CHKERR KSPSetType(ksp, KSPPREONLY); // Run KSP internally in ShellPC
3039 auto sub_pc = createPC(fsRawPtr->mField.get_comm());
3040 CHKERR PCSetType(sub_pc, PCSHELL);
3041 CHKERR PCShellSetSetUp(sub_pc, pcSetup);
3042 CHKERR PCShellSetApply(sub_pc, pcApply);
3043 CHKERR KSPSetPC(ksp, sub_pc);
3044 CHKERR SNESSetKSP(snes, ksp);
3045
3048
3049 CHKERR TSSetPreStep(ts, tsPreProc);
3050 CHKERR TSSetPostStep(ts, tsPostProc);
3051
3053}
static Index< 'p', 3 > p
std::string param_file
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
BoundaryEle::UserDataOperator BoundaryEleOp
Kronecker Delta class symmetric.
@ QUIET
Definition: definitions.h:208
#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
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
CoordinateTypes
Coodinate system.
Definition: definitions.h:114
@ CYLINDRICAL
Definition: definitions.h:117
#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
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: elastic.cpp:66
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
Definition: elastic.cpp:67
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
auto cylindrical
auto marker
auto bubble_device
auto save_range
auto get_M2
constexpr int U_FIELD_DIM
constexpr double mu_m
constexpr double rho_ave
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
FTensor::Index< 'j', SPACE_DIM > j
auto get_M2_dh
constexpr double rho_m
@ R
@ F
auto integration_rule
static char help[]
constexpr CoordinateTypes coord_type
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
surface tension
OpDomainMassH OpDomainMassG
auto get_M3_dh
FTensor::Index< 'i', SPACE_DIM > i
constexpr double mu_diff
auto init_h
Initialisation function.
constexpr auto t_kd
auto d_phase_function_h
auto get_fe_bit
constexpr int SPACE_DIM
constexpr double eta
FTensor::Index< 'l', SPACE_DIM > l
auto kernel_oscillation
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
auto get_skin_parent_bit
auto get_M0
constexpr bool debug
constexpr double a0
constexpr double tol
auto get_f_dh
auto get_M3
auto get_start_bit
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
OpDomainMassH OpDomainMassP
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
constexpr int BASE_DIM
int nb_levels
auto wetting_angle
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
constexpr double W
auto get_M
auto get_dofs_ents_all
get entities of dofs in the problem - used for debugging
auto phase_function
constexpr IntegrationType I
auto get_skin_projection_bit
auto get_global_size
constexpr double h
auto get_dofs_ents_by_field_name
get entities of dofs in the problem - used for debugging
ElementsAndOps< SPACE_DIM >::DomianParentEle DomianParentEle
constexpr AssemblyType A
constexpr double rho_p
constexpr double rho_diff
auto get_D
auto wetting_angle_sub_stepping
int refine_overlap
auto get_M0_dh
auto my_max
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
auto cut_off
auto get_skin_child_bit
constexpr double eta2
int order
approximation order
SideEle::UserDataOperator SideOp
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
auto get_current_bit
dofs bit used to do calculations
const double kappa
auto bit
constexpr double md
auto capillary_tube
auto get_f
auto get_projection_bit
auto d_cut_off
auto my_min
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:219
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:442
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:542
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.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: DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:400
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:275
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:988
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:532
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
MoFEMErrorCode getEntitiesByDimAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const int dim, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
const double c
speed of light (cm/ns)
double D
const double v
phase velocity of light in medium (cm/ns)
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomianParentEle DomianParentEle
Definition: level_set.cpp:25
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: Common.hpp:10
auto createKSP(MPI_Comm comm)
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1042
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobina in TS solver.
Definition: TsCtx.cpp:131
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:221
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:424
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:226
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
Definition: TsCtx.cpp:25
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVecScatter(Vec x, IS ix, Vec y, IS iy)
Create a Vec Scatter object.
auto createTS(MPI_Comm comm)
auto createPC(MPI_Comm comm)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
PetscObject getPetscObject(T obj)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1626
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1031
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:976
int save_every_nth_step
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
constexpr double g
Range findParentsToRefine(Range ents, BitRefLevel level, BitRefLevel mask)
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
FreeSurface(MoFEM::Interface &m_field)
MoFEMErrorCode projectData()
[Boundary condition]
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode setParentDofs(boost::shared_ptr< FEMethod > fe_top, std::string field_name, ForcesAndSourcesCore::UserDataOperator::OpType op, BitRefLevel child_ent_bit, boost::function< boost::shared_ptr< ForcesAndSourcesCore >()> get_elem, int verbosity, LogManager::SeverityLevel sev)
Create hierarchy of elements run on parents levels.
std::vector< Range > findEntitiesCrossedByPhaseInterface(size_t overlap)
Find entities on refinement levels.
MoFEMErrorCode assembleSystem()
[Data projection]
MoFEM::Interface & mField
MoFEMErrorCode refineMesh(size_t overlap)
MoFEMErrorCode makeRefProblem()
MoFEMErrorCode solveSystem()
[Solve]
Add operators pushing bases from local to physical configuration.
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
Managing BitRefLevels.
Managing BitRefLevels.
virtual int get_comm_size() const =0
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
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.
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Basic algebra on fields.
Definition: FieldBlas.hpp:21
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
@ OPSPACE
operator do Work is execute on space data
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition: Natural.hpp:65
Operator to project base functions from parent entity to child.
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 values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Element used to execute operators on side of the element.
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
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
intrusive_ptr for managing petsc objects
PetscReal ts_t
time
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.
boost::shared_ptr< moab::Core > postProcMesh
SmartPetscObj< DM > dM
MoFEMErrorCode postProcess()
function is run at the end of loop
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< moab::Core > post_proc_mesh, boost::shared_ptr< PostProcEleDomainCont > post_proc, boost::shared_ptr< PostProcEleBdyCont > post_proc_edge, std::pair< boost::shared_ptr< BoundaryEle >, boost::shared_ptr< VectorDouble > > p)
boost::shared_ptr< VectorDouble > liftVec
boost::shared_ptr< BoundaryEle > liftFE
boost::shared_ptr< PostProcEle > postProc
boost::shared_ptr< PostProcEleBdyCont > postProcEdge
boost::shared_ptr< PostProcEleDomainCont > postProc
Set of functions called by PETSc solver used to refine and update mesh.
boost::shared_ptr< TsCtx > tsCtxPtr
virtual ~TSPrePostProc()=default
SmartPetscObj< Vec > globRes
SmartPetscObj< Mat > subB
SmartPetscObj< Vec > globSol
SmartPetscObj< Vec > getSubVector()
static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Wrapper for TS monitor.
static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x)
SmartPetscObj< VecScatter > getScatter(Vec x, Vec y, enum FR fr)
SmartPetscObj< DM > solverSubDM
static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Wrapper for SNES Lhs.
FreeSurface * fsRawPtr
boost::shared_ptr< SnesCtx > snesCtxPtr
static MoFEMErrorCode tsPreProc(TS ts)
Pre process time step.
TSPrePostProc()=default
static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec f, void *ctx)
static MoFEMErrorCode pcSetup(PC pc)
static MoFEMErrorCode tsPostProc(TS ts)
Post process time step.
static MoFEMErrorCode tsPreStage(TS ts)
MoFEMErrorCode tsSetUp(TS ts)
Used to setup TS solver.
SmartPetscObj< KSP > subKSP