v0.14.0
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
75int coord_type = EXECUTABLE_COORD_TYPE;
76
77constexpr int BASE_DIM = 1;
78constexpr int SPACE_DIM = 2;
79constexpr int U_FIELD_DIM = SPACE_DIM;
80
81constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
82constexpr IntegrationType I =
83 IntegrationType::GAUSS; //< selected integration type
84
85template <int DIM>
87
95using SideOp = SideEle::UserDataOperator;
96
100
102
106
115
122
123template <CoordinateTypes COORD_TYPE>
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
160double a0 = 980;
161double rho_m = 0.998;
162double mu_m = 0.010101 * 1e2;
163double rho_p = 0.0012;
164double mu_p = 0.000182 * 1e2;
165double lambda = 73; ///< surface tension
166double W = 0.25;
167
168// Model parameters
169double h = 0.01 / 8; // mesh size
170double eta = h;
171double eta2 = eta * eta;
172
173// Numerical parameters
174double md = 1e-2; // mobility
175double eps = 1e-12;
176double tol = std::numeric_limits<float>::epsilon();
177
178double rho_ave = (rho_p + rho_m) / 2;
179double rho_diff = (rho_p - rho_m) / 2;
180double mu_ave = (mu_p + mu_m) / 2;
181double mu_diff = (mu_p - mu_m) / 2;
182
183double 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 if (coord_type == CYLINDRICAL)
189 return 2 * M_PI * r;
190 else
191 return 1.;
192};
193
194auto wetting_angle_sub_stepping = [](auto ts_step) {
195 constexpr int sub_stepping = 16;
196 return std::min(1., static_cast<double>(ts_step) / sub_stepping);
197};
198
199// cut off function
200
201auto my_max = [](const double x) { return (x - 1 + std::abs(x + 1)) / 2; };
202auto my_min = [](const double x) { return (x + 1 - std::abs(x - 1)) / 2; };
203auto cut_off = [](const double h) { return my_max(my_min(h)); };
204auto d_cut_off = [](const double h) {
205 if (h >= -1 && h < 1)
206 return 1.;
207 else
208 return 0.;
209};
210
211auto phase_function = [](const double h, const double diff, const double ave) {
212 return diff * h + ave;
213};
214
215auto d_phase_function_h = [](const double h, const double diff) {
216 return diff;
217};
218
219// order function (free energy)
220
221auto get_f = [](const double h) { return 4 * W * h * (h * h - 1); };
222auto get_f_dh = [](const double h) { return 4 * W * (3 * h * h - 1); };
223
224auto get_M0 = [](auto h) { return md; };
225auto get_M0_dh = [](auto h) { return 0; };
226
227auto get_M2 = [](auto h) {
228 return md * (1 - h * h);
229};
230
231auto get_M2_dh = [](auto h) { return -md * 2 * h; };
232
233auto get_M3 = [](auto h) {
234 const double h2 = h * h;
235 const double h3 = h2 * h;
236 if (h >= 0)
237 return md * (2 * h3 - 3 * h2 + 1);
238 else
239 return md * (-2 * h3 - 3 * h2 + 1);
240};
241
242auto get_M3_dh = [](auto h) {
243 if (h >= 0)
244 return md * (6 * h * (h - 1));
245 else
246 return md * (-6 * h * (h + 1));
247};
248
249auto get_M = [](auto h) { return get_M0(h); };
250auto get_M_dh = [](auto h) { return get_M0_dh(h); };
251
252auto get_D = [](const double A) {
254 t_D(i, j, k, l) = A * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
255 return t_D;
256};
257
258// some examples of initialisation functions
259
260auto kernel_oscillation = [](double r, double y, double) {
261 constexpr int n = 3;
262 constexpr double R = 0.0125;
263 constexpr double A = R * 0.2;
264 const double theta = atan2(r, y);
265 const double w = R + A * cos(n * theta);
266 const double d = std::sqrt(r * r + y * y);
267 return tanh((w - d) / (eta * std::sqrt(2)));
268};
269
270auto kernel_eye = [](double r, double y, double) {
271 constexpr double y0 = 0.5;
272 constexpr double R = 0.4;
273 y -= y0;
274 const double d = std::sqrt(r * r + y * y);
275 return tanh((R - d) / (eta * std::sqrt(2)));
276};
277
278auto capillary_tube = [](double x, double y, double z) {
279 constexpr double water_height = 0.;
280 return tanh((water_height - y) / (eta * std::sqrt(2)));
281 ;
282};
283
284auto bubble_device = [](double x, double y, double z) {
285 return -tanh((-0.039 - x) / (eta * std::sqrt(2)));
286};
287
288/**
289 * @brief Initialisation function
290 *
291 * @note If UMs are compiled with Python to initialise phase field "H"
292 * surface.py function is used, which has to be present in execution folder.
293 *
294 */
295auto init_h = [](double r, double y, double theta) {
296#ifdef PYTHON_INIT_SURFACE
297 double s = 1;
298 if (auto ptr = surfacePythonWeakPtr.lock()) {
299 CHK_THROW_MESSAGE(ptr->evalSurface(r, y, theta, eta, s),
300 "error eval python");
301 }
302 return s;
303#else
304 return bubble_device(r, y, theta);
305 // return capillary_tube(r, y, theta);
306 // return kernel_eye(r, y, theta);
307#endif
308};
309
310auto wetting_angle = [](double water_level) { return water_level; };
311
312auto bit = [](auto b) { return BitRefLevel().set(b); };
313auto marker = [](auto b) { return BitRefLevel().set(BITREFLEVEL_SIZE - b); };
314auto get_fe_bit = [](FEMethod *fe_ptr) {
315 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
316};
317
318auto get_global_size = [](int l_size) {
319 int g_size;
320 MPI_Allreduce(&l_size, &g_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
321 return g_size;
322};
323
324auto save_range = [](moab::Interface &moab, const std::string name,
325 const Range r) {
327 if (get_global_size(r.size())) {
328 auto out_meshset = get_temp_meshset_ptr(moab);
329 CHKERR moab.add_entities(*out_meshset, r);
330 CHKERR moab.write_file(name.c_str(), "MOAB", "PARALLEL=WRITE_PART",
331 out_meshset->get_ptr(), 1);
332 }
334};
335
336/**
337 * @brief get entities of dofs in the problem - used for debugging
338 *
339 */
340auto get_dofs_ents_by_field_name = [](auto dm, auto field_name) {
341 auto prb_ptr = getProblemPtr(dm);
342 std::vector<EntityHandle> ents_vec;
343
344 MoFEM::Interface *m_field_ptr;
345 CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
346
347 auto bit_number = m_field_ptr->get_field_bit_number(field_name);
348 auto dofs = prb_ptr->numeredRowDofsPtr;
349 auto lo_it = dofs->lower_bound(FieldEntity::getLoBitNumberUId(bit_number));
350 auto hi_it = dofs->upper_bound(FieldEntity::getHiBitNumberUId(bit_number));
351 ents_vec.reserve(std::distance(lo_it, hi_it));
352
353 for (; lo_it != hi_it; ++lo_it) {
354 ents_vec.push_back((*lo_it)->getEnt());
355 }
356
357 std::sort(ents_vec.begin(), ents_vec.end());
358 auto it = std::unique(ents_vec.begin(), ents_vec.end());
359 Range r;
360 r.insert_list(ents_vec.begin(), it);
361 return r;
362};
363
364/**
365 * @brief get entities of dofs in the problem - used for debugging
366 *
367 */
368auto get_dofs_ents_all = [](auto dm) {
369 auto prb_ptr = getProblemPtr(dm);
370 std::vector<EntityHandle> ents_vec;
371
372 MoFEM::Interface *m_field_ptr;
373 CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
374
375 auto dofs = prb_ptr->numeredRowDofsPtr;
376 ents_vec.reserve(dofs->size());
377
378 for (auto d : *dofs) {
379 ents_vec.push_back(d->getEnt());
380 }
381
382 std::sort(ents_vec.begin(), ents_vec.end());
383 auto it = std::unique(ents_vec.begin(), ents_vec.end());
384 Range r;
385 r.insert_list(ents_vec.begin(), it);
386 return r;
387};
388
389#include <FreeSurfaceOps.hpp>
390using namespace FreeSurfaceOps;
391
392struct FreeSurface;
393
394enum FR { F, R }; // F - forward, and reverse
395
396/**
397 * @brief Set of functions called by PETSc solver used to refine and update
398 * mesh.
399 *
400 * @note Currently theta method is only handled by this code.
401 *
402 */
403struct TSPrePostProc {
404
405 TSPrePostProc() = default;
406 virtual ~TSPrePostProc() = default;
407
408 /**
409 * @brief Used to setup TS solver
410 *
411 * @param ts
412 * @return MoFEMErrorCode
413 */
415
416 SmartPetscObj<VecScatter> getScatter(Vec x, Vec y, enum FR fr);
418
422
423private:
424 /**
425 * @brief Pre process time step
426 *
427 * Refine mesh and update fields
428 *
429 * @param ts
430 * @return MoFEMErrorCode
431 */
432 static MoFEMErrorCode tsPreProc(TS ts);
433
434 /**
435 * @brief Post process time step.
436 *
437 * Currently that function do not make anything major
438 *
439 * @param ts
440 * @return MoFEMErrorCode
441 */
442 static MoFEMErrorCode tsPostProc(TS ts);
443
445
446 static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t,
447 Vec f,
448 void *ctx); //< Wrapper for SNES Rhs
449 static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
450 PetscReal a, Mat A, Mat B,
451 void *ctx); ///< Wrapper for SNES Lhs
452 static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u,
453 void *ctx); ///< Wrapper for TS monitor
454 static MoFEMErrorCode pcSetup(PC pc);
455 static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x);
456
457 SmartPetscObj<Vec> globRes; //< global residual
458 SmartPetscObj<Mat> subB; //< sub problem tangent matrix
459 SmartPetscObj<KSP> subKSP; //< sub problem KSP solver
460
461 boost::shared_ptr<SnesCtx>
462 snesCtxPtr; //< infernal data (context) for MoFEM SNES fuctions
463 boost::shared_ptr<TsCtx>
464 tsCtxPtr; //< internal data (context) for MoFEM TS functions.
465};
466
467static boost::weak_ptr<TSPrePostProc> tsPrePostProc;
468
470
471 FreeSurface(MoFEM::Interface &m_field) : mField(m_field) {}
472
474
476
478
479private:
486
487 /// @brief Find entities on refinement levels
488 /// @param overlap level of overlap
489 /// @return
490 std::vector<Range> findEntitiesCrossedByPhaseInterface(size_t overlap);
491
492 /// @brief
493 /// @param ents
494 /// @param level
495 /// @param mask
496 /// @return
498
499 /// @brief
500 /// @param overlap
501 /// @return
502 MoFEMErrorCode refineMesh(size_t overlap);
503
504 /// @brief Create hierarchy of elements run on parents levels
505 /// @param fe_top pipeline element to which hierarchy is attached
506 /// @param op type of operator OPSPACE, OPROW, OPCOL or OPROWCOL
507 /// @param get_elema lambda function to create element on hierarchy
508 /// @param verbosity verbosity level
509 /// @param sev severity level
510 /// @return
512 boost::shared_ptr<FEMethod> fe_top, std::string field_name,
513 ForcesAndSourcesCore::UserDataOperator::OpType op,
514 BitRefLevel child_ent_bit,
515 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
516 int verbosity, LogManager::SeverityLevel sev);
517
518 friend struct TSPrePostProc;
519};
520
521//! [Run programme]
530}
531//! [Run programme]
532
533//! [Read mesh]
536 MOFEM_LOG("FS", Sev::inform) << "Read mesh for problem";
538
539 simple->getParentAdjacencies() = true;
540 simple->getBitRefLevel() = BitRefLevel();
541
542 CHKERR simple->getOptions();
543 CHKERR simple->loadFile();
544
546}
547//! [Read mesh]
548
549//! [Set up problem]
552
553 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
554 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-nb_levels", &nb_levels,
555 PETSC_NULL);
556 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-refine_overlap", &refine_overlap,
557 PETSC_NULL);
558
559 const char *coord_type_names[] = {"cartesian", "polar", "cylindrical",
560 "spherical"};
561 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-coords", coord_type_names,
562 LAST_COORDINATE_SYSTEM, &coord_type, PETSC_NULL);
563
564 MOFEM_LOG("FS", Sev::inform) << "Approximation order = " << order;
565 MOFEM_LOG("FS", Sev::inform)
566 << "Number of refinement levels nb_levels = " << nb_levels;
567 nb_levels += 1;
568
570 auto bit_mng = mField.getInterface<BitRefManager>();
571
572 CHKERR PetscOptionsGetScalar(PETSC_NULL, "Acceleration", "-a0", &a0, PETSC_NULL);
573 CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Minus\" phase density rho_m",
574 "-rho_m", &rho_m, PETSC_NULL);
575 CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Minus\" phase viscosity", "-mu_m",
576 &mu_m, PETSC_NULL);
577 CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Plus\" phase density", "-rho_p",
578 &rho_p, PETSC_NULL);
579 CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Plus\" phase viscosity", "-mu_p",
580 &mu_p, PETSC_NULL);
581 CHKERR PetscOptionsGetScalar(PETSC_NULL, "Surface tension", "-lambda",
582 &lambda, PETSC_NULL);
583 CHKERR PetscOptionsGetScalar(PETSC_NULL,
584 "Height of the well in energy functional", "-W",
585 &W, PETSC_NULL);
586 rho_ave = (rho_p + rho_m) / 2;
587 rho_diff = (rho_p - rho_m) / 2;
588 mu_ave = (mu_p + mu_m) / 2;
589 mu_diff = (mu_p - mu_m) / 2;
590
591 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-h", &h, PETSC_NULL);
592 eta = h;
593 eta2 = eta * eta;
594 kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
595
596 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-md", &eta, PETSC_NULL);
597
598 MOFEM_LOG("FS", Sev::inform) << "Acceleration a0 = " << a0;
599 MOFEM_LOG("FS", Sev::inform) << "\"Minus\" phase density rho_m = " << rho_m;
600 MOFEM_LOG("FS", Sev::inform) << "\"Minus\" phase viscosity mu_m = " << mu_m;
601 MOFEM_LOG("FS", Sev::inform) << "\"Plus\" phase density rho_p = " << rho_p;
602 MOFEM_LOG("FS", Sev::inform) << "\"Plus\" phase viscosity mu_p = " << mu_p;
603 MOFEM_LOG("FS", Sev::inform) << "Surface tension lambda = " << lambda;
604 MOFEM_LOG("FS", Sev::inform)
605 << "Height of the well in energy functional W = " << W;
606 MOFEM_LOG("FS", Sev::inform) << "Characteristic mesh size h = " << h;
607 MOFEM_LOG("FS", Sev::inform) << "Mobility md = " << md;
608
609 MOFEM_LOG("FS", Sev::inform) << "Average density rho_ave = " << rho_ave;
610 MOFEM_LOG("FS", Sev::inform) << "Difference density rho_diff = " << rho_diff;
611 MOFEM_LOG("FS", Sev::inform) << "Average viscosity mu_ave = " << mu_ave;
612 MOFEM_LOG("FS", Sev::inform) << "Difference viscosity mu_diff = " << mu_diff;
613 MOFEM_LOG("FS", Sev::inform) << "kappa = " << kappa;
614
615
616 // Fields on domain
617
618 // Velocity field
619 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, U_FIELD_DIM);
620 // Pressure field
621 CHKERR simple->addDomainField("P", H1, AINSWORTH_LEGENDRE_BASE, 1);
622 // Order/phase fild
623 CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
624 // Chemical potential
625 CHKERR simple->addDomainField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
626
627 // Field on boundary
628 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE,
630 CHKERR simple->addBoundaryField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
631 CHKERR simple->addBoundaryField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
632 // Lagrange multiplier which constrains slip conditions
633 CHKERR simple->addBoundaryField("L", H1, AINSWORTH_LEGENDRE_BASE, 1);
634
635 CHKERR simple->setFieldOrder("U", order);
636 CHKERR simple->setFieldOrder("P", order - 1);
637 CHKERR simple->setFieldOrder("H", order);
638 CHKERR simple->setFieldOrder("G", order);
639 CHKERR simple->setFieldOrder("L", order);
640
641 // Initialise bit ref levels
642 auto set_problem_bit = [&]() {
644 // Set bits to build adjacencies between parents and children. That is
645 // used by simple interface.
646 simple->getBitAdjEnt() = BitRefLevel().set();
647 simple->getBitAdjParent() = BitRefLevel().set();
648 simple->getBitRefLevel() = BitRefLevel().set();
649 simple->getBitRefLevelMask() = BitRefLevel().set();
651 };
652
653 CHKERR set_problem_bit();
654
655 CHKERR simple->setUp();
656
658}
659//! [Set up problem]
660
661//! [Boundary condition]
664
665#ifdef PYTHON_INIT_SURFACE
666 auto get_py_surface_init = []() {
667 auto py_surf_init = boost::make_shared<SurfacePython>();
668 CHKERR py_surf_init->surfaceInit("surface.py");
669 surfacePythonWeakPtr = py_surf_init;
670 return py_surf_init;
671 };
672 auto py_surf_init = get_py_surface_init();
673#endif
674
676 auto pip_mng = mField.getInterface<PipelineManager>();
677 auto bc_mng = mField.getInterface<BcManager>();
678 auto bit_mng = mField.getInterface<BitRefManager>();
679 auto dm = simple->getDM();
680
682
683 auto reset_bits = [&]() {
685 BitRefLevel start_mask;
686 for (auto s = 0; s != get_start_bit(); ++s)
687 start_mask[s] = true;
688 // reset bit ref levels
689 CHKERR bit_mng->lambdaBitRefLevel(
690 [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
691 Range level0;
692 CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
693 CHKERR bit_mng->setNthBitRefLevel(level0, get_current_bit(), true);
694 CHKERR bit_mng->setNthBitRefLevel(level0, get_projection_bit(), true);
696 };
697
698 auto add_parent_field = [&](auto fe, auto op, auto field) {
699 return setParentDofs(
700 fe, field, op, bit(get_skin_parent_bit()),
701
702 [&]() {
703 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
705 return fe_parent;
706 },
707
708 QUIET, Sev::noisy);
709 };
710
711 auto h_ptr = boost::make_shared<VectorDouble>();
712 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
713 auto g_ptr = boost::make_shared<VectorDouble>();
714 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
715
716 auto set_generic = [&](auto fe) {
718 auto &pip = fe->getOpPtrVector();
719
721
723 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
724
725 [&]() {
726 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
729 fe_parent->getOpPtrVector(), {H1});
730 return fe_parent;
731 },
732
733 QUIET, Sev::noisy);
734
735 CHKERR add_parent_field(fe, UDO::OPROW, "H");
736 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
737 pip.push_back(
738 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
739
740 CHKERR add_parent_field(fe, UDO::OPROW, "G");
741 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
742 pip.push_back(
743 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
744
746 };
747
748 auto post_proc = [&](auto exe_test) {
750 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
751 post_proc_fe->exeTestHook = exe_test;
752
753 CHKERR set_generic(post_proc_fe);
754
756
757 post_proc_fe->getOpPtrVector().push_back(
758
759 new OpPPMap(post_proc_fe->getPostProcMesh(),
760 post_proc_fe->getMapGaussPts(),
761
762 {{"H", h_ptr}, {"G", g_ptr}},
763
764 {{"GRAD_H", grad_h_ptr}, {"GRAD_G", grad_g_ptr}},
765
766 {},
767
768 {}
769
770 )
771
772 );
773
774 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
775 CHKERR post_proc_fe->writeFile("out_init.h5m");
776
778 };
779
780 auto solve_init = [&](auto exe_test) {
782
783 pip_mng->getOpDomainRhsPipeline().clear();
784 pip_mng->getOpDomainLhsPipeline().clear();
785
786 auto set_domain_rhs = [&](auto fe) {
788 CHKERR set_generic(fe);
789 auto &pip = fe->getOpPtrVector();
790
791 CHKERR add_parent_field(fe, UDO::OPROW, "H");
792 pip.push_back(new OpRhsH<true>("H", nullptr, nullptr, h_ptr, grad_h_ptr,
793 grad_g_ptr));
794 CHKERR add_parent_field(fe, UDO::OPROW, "G");
795 pip.push_back(new OpRhsG<true>("G", h_ptr, grad_h_ptr, g_ptr));
797 };
798
799 auto set_domain_lhs = [&](auto fe) {
801
802 CHKERR set_generic(fe);
803 auto &pip = fe->getOpPtrVector();
804
805 CHKERR add_parent_field(fe, UDO::OPROW, "H");
806 CHKERR add_parent_field(fe, UDO::OPCOL, "H");
807 pip.push_back(new OpLhsH_dH<true>("H", nullptr, h_ptr, grad_g_ptr));
808
809 CHKERR add_parent_field(fe, UDO::OPCOL, "G");
810 pip.push_back(new OpLhsH_dG<true>("H", "G", h_ptr));
811
812 CHKERR add_parent_field(fe, UDO::OPROW, "G");
813 pip.push_back(new OpLhsG_dG("G"));
814
815 CHKERR add_parent_field(fe, UDO::OPCOL, "H");
816 pip.push_back(new OpLhsG_dH<true>("G", "H", h_ptr));
818 };
819
820 auto create_subdm = [&]() {
821 auto level_ents_ptr = boost::make_shared<Range>();
822 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
823 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
824
825 DM subdm;
826 CHKERR DMCreate(mField.get_comm(), &subdm);
827 CHKERR DMSetType(subdm, "DMMOFEM");
828 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_INIT");
829 CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
830 CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
831 CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_TRUE);
832
833 for (auto f : {"H", "G"}) {
834 CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
835 CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
836 }
837 CHKERR DMSetUp(subdm);
838
839 if constexpr (debug) {
840 if (mField.get_comm_size() == 1) {
841 auto dm_ents = get_dofs_ents_all(SmartPetscObj<DM>(subdm, true));
842 CHKERR save_range(mField.get_moab(), "sub_dm_init_ents_verts.h5m",
843 dm_ents.subset_by_type(MBVERTEX));
844 dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
845 CHKERR save_range(mField.get_moab(), "sub_dm_init_ents.h5m", dm_ents);
846 }
847 }
848
849 return SmartPetscObj<DM>(subdm);
850 };
851
852 auto subdm = create_subdm();
853
854 pip_mng->getDomainRhsFE().reset();
855 pip_mng->getDomainLhsFE().reset();
856 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
857 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
858 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
859 pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
860
861 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
862 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
863
864 auto D = createDMVector(subdm);
865 auto snes = pip_mng->createSNES(subdm);
866 auto snes_ctx_ptr = getDMSnesCtx(subdm);
867
868 auto set_section_monitor = [&](auto solver) {
870 CHKERR SNESMonitorSet(solver,
871 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
873 (void *)(snes_ctx_ptr.get()), nullptr);
874
876 };
877
878 auto print_section_field = [&]() {
880
881 auto section =
882 mField.getInterface<ISManager>()->sectionCreate("SUB_INIT");
883 PetscInt num_fields;
884 CHKERR PetscSectionGetNumFields(section, &num_fields);
885 for (int f = 0; f < num_fields; ++f) {
886 const char *field_name;
887 CHKERR PetscSectionGetFieldName(section, f, &field_name);
888 MOFEM_LOG("FS", Sev::inform)
889 << "Field " << f << " " << std::string(field_name);
890 }
891
893 };
894
895 CHKERR set_section_monitor(snes);
896 CHKERR print_section_field();
897
898 for (auto f : {"U", "P", "H", "G", "L"}) {
899 CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
900 }
901
902 CHKERR SNESSetFromOptions(snes);
903 CHKERR SNESSolve(snes, PETSC_NULL, D);
904
905 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
906 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
907 CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
908
910 };
911
912 CHKERR reset_bits();
913 CHKERR solve_init(
914 [](FEMethod *fe_ptr) { return get_fe_bit(fe_ptr).test(0); });
915 CHKERR refineMesh(refine_overlap);
916
917 for (auto f : {"U", "P", "H", "G", "L"}) {
918 CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
919 }
920 CHKERR solve_init([](FEMethod *fe_ptr) {
921 return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
922 });
923
924 CHKERR post_proc([](FEMethod *fe_ptr) {
925 return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
926 });
927
928 pip_mng->getOpDomainRhsPipeline().clear();
929 pip_mng->getOpDomainLhsPipeline().clear();
930
931 // Remove DOFs where boundary conditions are set
932 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
933 "U", 0, 0);
934 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
935 "L", 0, 0);
936 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
937 "U", 1, 1);
938 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
939 "L", 1, 1);
940 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "U",
941 0, SPACE_DIM);
942 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "L",
943 0, 0);
944 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "ZERO",
945 "L", 0, 0);
946
948}
949//! [Boundary condition]
950
951//! [Data projection]
954
955 // FIXME: Functionality in this method should be improved (projection should
956 // be done field by field), generalised, and move to become core
957 // functionality.
958
960 auto pip_mng = mField.getInterface<PipelineManager>();
961 auto bc_mng = mField.getInterface<BcManager>();
962 auto bit_mng = mField.getInterface<BitRefManager>();
963 auto field_blas = mField.getInterface<FieldBlas>();
964
965 // Store all existing elements pipelines, replace them by data projection
966 // pipelines, to put them back when projection is done.
967 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
968 auto fe_domain_rhs = pip_mng->getDomainRhsFE();
969 auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
970 auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
971
972 pip_mng->getDomainLhsFE().reset();
973 pip_mng->getDomainRhsFE().reset();
974 pip_mng->getBoundaryLhsFE().reset();
975 pip_mng->getBoundaryRhsFE().reset();
976
978
979 // extract field data for domain parent element
980 auto add_parent_field_domain = [&](auto fe, auto op, auto field, auto bit) {
981 return setParentDofs(
982 fe, field, op, bit,
983
984 [&]() {
985 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
987 return fe_parent;
988 },
989
990 QUIET, Sev::noisy);
991 };
992
993 // extract field data for boundary parent element
994 auto add_parent_field_bdy = [&](auto fe, auto op, auto field, auto bit) {
995 return setParentDofs(
996 fe, field, op, bit,
997
998 [&]() {
999 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1001 return fe_parent;
1002 },
1003
1004 QUIET, Sev::noisy);
1005 };
1006
1007 // run this element on element with given bit, otherwise run other nested
1008 // element
1009 auto create_run_parent_op = [&](auto parent_fe_ptr, auto this_fe_ptr,
1010 auto fe_bit) {
1011 auto parent_mask = fe_bit;
1012 parent_mask.flip();
1013 return new OpRunParent(parent_fe_ptr, BitRefLevel().set(), parent_mask,
1014 this_fe_ptr, fe_bit, BitRefLevel().set(), QUIET,
1015 Sev::inform);
1016 };
1017
1018 // create hierarchy of nested elements
1019 auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr, auto fe_bit) {
1020 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
1021 for (int l = 0; l < nb_levels; ++l)
1022 parents_elems_ptr_vec.emplace_back(
1023 boost::make_shared<DomianParentEle>(mField));
1024 for (auto l = 1; l < nb_levels; ++l) {
1025 parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
1026 create_run_parent_op(parents_elems_ptr_vec[l], this_fe_ptr, fe_bit));
1027 }
1028 return parents_elems_ptr_vec[0];
1029 };
1030
1031 // solve projection problem
1032 auto solve_projection = [&](auto exe_test) {
1034
1035 auto set_domain_rhs = [&](auto fe) {
1037 auto &pip = fe->getOpPtrVector();
1038
1039 auto u_ptr = boost::make_shared<MatrixDouble>();
1040 auto p_ptr = boost::make_shared<VectorDouble>();
1041 auto h_ptr = boost::make_shared<VectorDouble>();
1042 auto g_ptr = boost::make_shared<VectorDouble>();
1043
1044 auto eval_fe_ptr = boost::make_shared<DomianParentEle>(mField);
1045 {
1046 auto &pip = eval_fe_ptr->getOpPtrVector();
1049 eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
1050
1051 [&]() {
1052 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1053 new DomianParentEle(mField));
1054 return fe_parent;
1055 },
1056
1057 QUIET, Sev::noisy);
1058 // That can be done much smarter, by block, field by field. For
1059 // simplicity is like that.
1060 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "U",
1062 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1063 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "P",
1065 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1066 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "H",
1068 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1069 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "G",
1071 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1072 }
1073 auto parent_eval_fe_ptr =
1074 get_parents_vel_fe_ptr(eval_fe_ptr, bit(get_projection_bit()));
1075 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1077
1078 auto assemble_fe_ptr = boost::make_shared<DomianParentEle>(mField);
1079 {
1080 auto &pip = assemble_fe_ptr->getOpPtrVector();
1083 assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1084
1085 [&]() {
1086 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1087 new DomianParentEle(mField));
1088 return fe_parent;
1089 },
1090
1091 QUIET, Sev::noisy);
1092 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "U",
1094 pip.push_back(new OpDomainAssembleVector("U", u_ptr));
1095 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "P",
1097 pip.push_back(new OpDomainAssembleScalar("P", p_ptr));
1098 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "H",
1100 pip.push_back(new OpDomainAssembleScalar("H", h_ptr));
1101 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "G",
1103 pip.push_back(new OpDomainAssembleScalar("G", g_ptr));
1104 }
1105 auto parent_assemble_fe_ptr =
1106 get_parents_vel_fe_ptr(assemble_fe_ptr, bit(get_current_bit()));
1107 pip.push_back(create_run_parent_op(
1108 parent_assemble_fe_ptr, assemble_fe_ptr, bit(get_current_bit())));
1109
1111 };
1112
1113 auto set_domain_lhs = [&](auto fe) {
1115
1116 auto &pip = fe->getOpPtrVector();
1117
1119
1121 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1122
1123 [&]() {
1124 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1125 new DomianParentEle(mField));
1126 return fe_parent;
1127 },
1128
1129 QUIET, Sev::noisy);
1130
1131 // That can be done much smarter, by block, field by field. For simplicity
1132 // is like that.
1133 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U",
1135 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U",
1137 pip.push_back(new OpDomainMassU("U", "U"));
1138 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P",
1140 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P",
1142 pip.push_back(new OpDomainMassP("P", "P"));
1143 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H",
1145 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H",
1147 pip.push_back(new OpDomainMassH("H", "H"));
1148 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G",
1150 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G",
1152 pip.push_back(new OpDomainMassG("G", "G"));
1153
1155 };
1156
1157 auto set_bdy_rhs = [&](auto fe) {
1159 auto &pip = fe->getOpPtrVector();
1160
1161 auto l_ptr = boost::make_shared<VectorDouble>();
1162
1163 auto eval_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
1164 {
1165 auto &pip = eval_fe_ptr->getOpPtrVector();
1167 eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
1168
1169 [&]() {
1170 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1172 return fe_parent;
1173 },
1174
1175 QUIET, Sev::noisy);
1176 // That can be done much smarter, by block, field by field. For
1177 // simplicity is like that.
1178 CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPROW, "L",
1180 pip.push_back(new OpCalculateScalarFieldValues("L", l_ptr));
1181 }
1182 auto parent_eval_fe_ptr =
1183 get_parents_vel_fe_ptr(eval_fe_ptr, bit(get_projection_bit()));
1184 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1186
1187 auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
1188 {
1189 auto &pip = assemble_fe_ptr->getOpPtrVector();
1191 assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1192
1193 [&]() {
1194 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1196 return fe_parent;
1197 },
1198
1199 QUIET, Sev::noisy);
1200
1201 struct OpLSize : public BoundaryEleOp {
1202 OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
1203 : BoundaryEleOp(NOSPACE, DomainEleOp::OPSPACE), lPtr(l_ptr) {}
1204 MoFEMErrorCode doWork(int, EntityType, EntData &) {
1206 if (lPtr->size() != getGaussPts().size2()) {
1207 lPtr->resize(getGaussPts().size2());
1208 lPtr->clear();
1209 }
1211 }
1212
1213 private:
1214 boost::shared_ptr<VectorDouble> lPtr;
1215 };
1216
1217 pip.push_back(new OpLSize(l_ptr));
1218
1219 CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW, "L",
1221 pip.push_back(new OpBoundaryAssembleScalar("L", l_ptr));
1222 }
1223 auto parent_assemble_fe_ptr =
1224 get_parents_vel_fe_ptr(assemble_fe_ptr, bit(get_current_bit()));
1225 pip.push_back(create_run_parent_op(
1226 parent_assemble_fe_ptr, assemble_fe_ptr, bit(get_current_bit())));
1227
1229 };
1230
1231 auto set_bdy_lhs = [&](auto fe) {
1233
1234 auto &pip = fe->getOpPtrVector();
1235
1237 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1238
1239 [&]() {
1240 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1242 return fe_parent;
1243 },
1244
1245 QUIET, Sev::noisy);
1246
1247 // That can be done much smarter, by block, field by field. For simplicity
1248 // is like that.
1249 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L",
1251 CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "L",
1253 pip.push_back(new OpBoundaryMassL("L", "L"));
1254
1256 };
1257
1258 auto create_subdm = [&]() -> SmartPetscObj<DM> {
1259 auto level_ents_ptr = boost::make_shared<Range>();
1260 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
1261 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
1262
1263 auto get_prj_ents = [&]() {
1264 Range prj_mesh;
1265 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_projection_bit()),
1266 BitRefLevel().set(),
1267 SPACE_DIM, prj_mesh);
1268 auto common_ents = intersect(prj_mesh, *level_ents_ptr);
1269 prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
1270 .subset_by_dimension(SPACE_DIM);
1271
1272 return prj_mesh;
1273 };
1274
1275 auto prj_ents = get_prj_ents();
1276
1277 if (get_global_size(prj_ents.size())) {
1278
1279 auto rebuild = [&]() {
1280 auto prb_mng = mField.getInterface<ProblemsManager>();
1282
1283 std::vector<std::string> fields{"U", "P", "H", "G", "L"};
1284 std::map<std::string, boost::shared_ptr<Range>> range_maps{
1285
1286 {"U", level_ents_ptr},
1287 {"P", level_ents_ptr},
1288 {"H", level_ents_ptr},
1289 {"G", level_ents_ptr},
1290 {"L", level_ents_ptr}
1291
1292 };
1293
1294 CHKERR prb_mng->buildSubProblem("SUB_SOLVER", fields, fields,
1295 simple->getProblemName(), PETSC_TRUE,
1296 &range_maps, &range_maps);
1297
1298 // partition problem
1299 CHKERR prb_mng->partitionFiniteElements("SUB_SOLVER", true, 0,
1301 // set ghost nodes
1302 CHKERR prb_mng->partitionGhostDofsOnDistributedMesh("SUB_SOLVER");
1303
1305 };
1306
1307 MOFEM_LOG("FS", Sev::verbose) << "Create projection problem";
1308
1309 CHKERR rebuild();
1310
1311 auto dm = simple->getDM();
1312 DM subdm;
1313 CHKERR DMCreate(mField.get_comm(), &subdm);
1314 CHKERR DMSetType(subdm, "DMMOFEM");
1315 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
1316 return SmartPetscObj<DM>(subdm);
1317 }
1318
1319 MOFEM_LOG("FS", Sev::inform) << "Nothing to project";
1320
1321 return SmartPetscObj<DM>();
1322 };
1323
1324 auto create_dummy_dm = [&]() {
1325 auto dummy_dm = createDM(mField.get_comm(), "DMMOFEM");
1327 simple->getProblemName().c_str(),
1329 "create dummy dm");
1330 return dummy_dm;
1331 };
1332
1333 auto subdm = create_subdm();
1334 if (subdm) {
1335
1336 pip_mng->getDomainRhsFE().reset();
1337 pip_mng->getDomainLhsFE().reset();
1338 pip_mng->getBoundaryRhsFE().reset();
1339 pip_mng->getBoundaryLhsFE().reset();
1340 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
1341 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
1342 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
1343 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
1344 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1345 pip_mng->getDomainRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
1346 return get_fe_bit(fe_ptr).test(nb_levels - 1);
1347 };
1348 pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
1349 pip_mng->getBoundaryRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
1350 return get_fe_bit(fe_ptr).test(nb_levels - 1);
1351 };
1352
1353 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1354 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1355 CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
1356 CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
1357
1358 auto D = createDMVector(subdm);
1359 auto F = vectorDuplicate(D);
1360
1361 auto ksp = pip_mng->createKSP(subdm);
1362 CHKERR KSPSetFromOptions(ksp);
1363 CHKERR KSPSetUp(ksp);
1364
1365 // get vector norm
1366 auto get_norm = [&](auto x) {
1367 double nrm;
1368 CHKERR VecNorm(x, NORM_2, &nrm);
1369 return nrm;
1370 };
1371
1372 /**
1373 * @brief Zero DOFs, used by FieldBlas
1374 */
1375 auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
1377 for (auto &v : ent_ptr->getEntFieldData()) {
1378 v = 0;
1379 }
1381 };
1382
1383 auto solve = [&](auto S) {
1385 CHKERR VecZeroEntries(S);
1386 CHKERR VecZeroEntries(F);
1387 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
1388 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
1389 CHKERR KSPSolve(ksp, F, S);
1390 CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
1391 CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
1392
1393
1394
1396 };
1397
1398 MOFEM_LOG("FS", Sev::inform) << "Solve projection";
1399 CHKERR solve(D);
1400
1401 auto glob_x = createDMVector(simple->getDM());
1402 auto sub_x = createDMVector(subdm);
1403 auto dummy_dm = create_dummy_dm();
1404
1405 /**
1406 * @brief get TSTheta data operators
1407 */
1408 auto apply_restrict = [&]() {
1409 auto get_is = [](auto v) {
1410 IS iy;
1411 auto create = [&]() {
1413 int n, ystart;
1414 CHKERR VecGetLocalSize(v, &n);
1415 CHKERR VecGetOwnershipRange(v, &ystart, NULL);
1416 CHKERR ISCreateStride(PETSC_COMM_SELF, n, ystart, 1, &iy);
1418 };
1419 CHK_THROW_MESSAGE(create(), "create is");
1420 return SmartPetscObj<IS>(iy);
1421 };
1422
1423 auto iy = get_is(glob_x);
1424 auto s = createVecScatter(glob_x, PETSC_NULL, glob_x, iy);
1425
1427 DMSubDomainRestrict(simple->getDM(), s, PETSC_NULL, dummy_dm),
1428 "restrict");
1429 Vec X0, Xdot;
1430 CHK_THROW_MESSAGE(DMGetNamedGlobalVector(dummy_dm, "TSTheta_X0", &X0),
1431 "get X0");
1433 DMGetNamedGlobalVector(dummy_dm, "TSTheta_Xdot", &Xdot),
1434 "get Xdot");
1435
1436 auto forward_ghost = [](auto g) {
1438 CHKERR VecGhostUpdateBegin(g, INSERT_VALUES, SCATTER_FORWARD);
1439 CHKERR VecGhostUpdateEnd(g, INSERT_VALUES, SCATTER_FORWARD);
1441 };
1442
1443 CHK_THROW_MESSAGE(forward_ghost(X0), "");
1444 CHK_THROW_MESSAGE(forward_ghost(Xdot), "");
1445
1446 if constexpr (debug) {
1447 MOFEM_LOG("FS", Sev::inform)
1448 << "Reverse restrict: X0 " << get_norm(X0) << " Xdot "
1449 << get_norm(Xdot);
1450 }
1451
1452 return std::vector<Vec>{X0, Xdot};
1453 };
1454
1455 auto ts_solver_vecs = apply_restrict();
1456
1457 if (ts_solver_vecs.size()) {
1458
1459 int ii = 0;
1460 for (auto v : ts_solver_vecs) {
1461 MOFEM_LOG("FS", Sev::inform) << "Solve projection vector";
1462
1463 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1464 SCATTER_REVERSE);
1465 CHKERR solve(sub_x);
1466
1467 for (auto f : {"U", "P", "H", "G", "L"}) {
1468 MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1469 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1470 }
1471
1472 CHKERR DMoFEMMeshToLocalVector(subdm, sub_x, INSERT_VALUES,
1473 SCATTER_REVERSE);
1474 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1475 SCATTER_FORWARD);
1476
1477 MOFEM_LOG("FS", Sev::inform) << "Norm V " << get_norm(v);
1478 }
1479
1480 CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_X0",
1481 &ts_solver_vecs[0]);
1482 CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_Xdot",
1483 &ts_solver_vecs[1]);
1484 }
1485
1486 for (auto f : {"U", "P", "H", "G", "L"}) {
1487 MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1488 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1489 }
1490 CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
1491 }
1492
1494 };
1495
1496 // postprocessing (only for debugging proposes)
1497 auto post_proc = [&](auto exe_test) {
1499 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
1500 auto &pip = post_proc_fe->getOpPtrVector();
1501 post_proc_fe->exeTestHook = exe_test;
1502
1504
1506 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1507
1508 [&]() {
1509 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1510 new DomianParentEle(mField));
1512 fe_parent->getOpPtrVector(), {H1});
1513 return fe_parent;
1514 },
1515
1516 QUIET, Sev::noisy);
1517
1519
1520 auto u_ptr = boost::make_shared<MatrixDouble>();
1521 auto p_ptr = boost::make_shared<VectorDouble>();
1522 auto h_ptr = boost::make_shared<VectorDouble>();
1523 auto g_ptr = boost::make_shared<VectorDouble>();
1524
1525 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U",
1527 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1528 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P",
1530 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1531 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H",
1533 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1534 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G",
1536 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1537
1538 post_proc_fe->getOpPtrVector().push_back(
1539
1540 new OpPPMap(post_proc_fe->getPostProcMesh(),
1541 post_proc_fe->getMapGaussPts(),
1542
1543 {{"P", p_ptr}, {"H", h_ptr}, {"G", g_ptr}},
1544
1545 {{"U", u_ptr}},
1546
1547 {},
1548
1549 {}
1550
1551 )
1552
1553 );
1554
1555 auto dm = simple->getDM();
1556 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
1557 CHKERR post_proc_fe->writeFile("out_projection.h5m");
1558
1560 };
1561
1562 CHKERR solve_projection([](FEMethod *fe_ptr) {
1563 return get_fe_bit(fe_ptr).test(get_current_bit());
1564 });
1565
1566 if constexpr (debug) {
1567 CHKERR post_proc([](FEMethod *fe_ptr) {
1568 return get_fe_bit(fe_ptr).test(get_current_bit());
1569 });
1570 }
1571
1572 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
1573 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
1574 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
1575 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
1576
1578}
1579//! [Data projection]
1580
1581//! [Push operators to pip]
1584 auto simple = mField.getInterface<Simple>();
1585
1587
1588 auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
1589 return setParentDofs(
1590 fe, field, op, bit(get_skin_parent_bit()),
1591
1592 [&]() {
1593 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1594 new DomianParentEle(mField));
1595 return fe_parent;
1596 },
1597
1598 QUIET, Sev::noisy);
1599 };
1600
1601 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
1602 return setParentDofs(
1603 fe, field, op, bit(get_skin_parent_bit()),
1604
1605 [&]() {
1606 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1608 return fe_parent;
1609 },
1610
1611 QUIET, Sev::noisy);
1612 };
1613
1614 auto test_bit_child = [](FEMethod *fe_ptr) {
1615 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1616 get_start_bit() + nb_levels - 1);
1617 };
1618
1619 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
1620 auto u_ptr = boost::make_shared<MatrixDouble>();
1621 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
1622 auto dot_h_ptr = boost::make_shared<VectorDouble>();
1623 auto h_ptr = boost::make_shared<VectorDouble>();
1624 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
1625 auto g_ptr = boost::make_shared<VectorDouble>();
1626 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
1627 auto lambda_ptr = boost::make_shared<VectorDouble>();
1628 auto p_ptr = boost::make_shared<VectorDouble>();
1629 auto div_u_ptr = boost::make_shared<VectorDouble>();
1630
1631 // Push element from reference configuration to current configuration in 3d
1632 // space
1633 auto set_domain_general = [&](auto fe) {
1635 auto &pip = fe->getOpPtrVector();
1636
1638
1640 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1641
1642 [&]() {
1643 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1644 new DomianParentEle(mField));
1646 fe_parent->getOpPtrVector(), {H1});
1647 return fe_parent;
1648 },
1649
1650 QUIET, Sev::noisy);
1651
1652 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1653 pip.push_back(
1655 pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
1657 "U", grad_u_ptr));
1658
1659 switch (coord_type) {
1660 case CARTESIAN:
1661 pip.push_back(
1663 "U", div_u_ptr));
1664 break;
1665 case CYLINDRICAL:
1666 pip.push_back(
1668 "U", div_u_ptr));
1669 break;
1670 default:
1671 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
1672 "Coordinate system not implemented");
1673 }
1674
1675 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1676 pip.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
1677 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1678 pip.push_back(
1679 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1680
1681 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1682 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1683 pip.push_back(
1684 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
1685
1686 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1687 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1689 };
1690
1691 auto set_domain_rhs = [&](auto fe) {
1693 auto &pip = fe->getOpPtrVector();
1694
1695 CHKERR set_domain_general(fe);
1696
1697 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1698 pip.push_back(new OpRhsU("U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
1699 grad_h_ptr, g_ptr, p_ptr));
1700
1701 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1702 pip.push_back(new OpRhsH<false>("H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
1703 grad_g_ptr));
1704
1705 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1706 pip.push_back(new OpRhsG<false>("G", h_ptr, grad_h_ptr, g_ptr));
1707
1708 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1709 pip.push_back(new OpDomainAssembleScalar(
1710 "P", div_u_ptr, [](const double r, const double, const double) {
1711 return cylindrical(r);
1712 }));
1713 pip.push_back(new OpDomainAssembleScalar(
1714 "P", p_ptr, [](const double r, const double, const double) {
1715 return eps * cylindrical(r);
1716 }));
1718 };
1719
1720 auto set_domain_lhs = [&](auto fe) {
1722 auto &pip = fe->getOpPtrVector();
1723
1724 CHKERR set_domain_general(fe);
1725
1726 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1727 {
1728 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1729 pip.push_back(new OpLhsU_dU("U", u_ptr, grad_u_ptr, h_ptr));
1730 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1731 pip.push_back(
1732 new OpLhsU_dH("U", "H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
1733
1734 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1735 pip.push_back(new OpLhsU_dG("U", "G", grad_h_ptr));
1736 }
1737
1738 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1739 {
1740 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1741 pip.push_back(new OpLhsH_dU("H", "U", grad_h_ptr));
1742 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1743 pip.push_back(new OpLhsH_dH<false>("H", u_ptr, h_ptr, grad_g_ptr));
1744 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1745 pip.push_back(new OpLhsH_dG<false>("H", "G", h_ptr));
1746 }
1747
1748 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1749 {
1750 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1751 pip.push_back(new OpLhsG_dH<false>("G", "H", h_ptr));
1752 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1753 pip.push_back(new OpLhsG_dG("G"));
1754 }
1755
1756 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1757 {
1758 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1759
1760 switch (coord_type) {
1761 case CARTESIAN:
1762 pip.push_back(new OpMixScalarTimesDiv<CARTESIAN>(
1763 "P", "U",
1764 [](const double r, const double, const double) {
1765 return cylindrical(r);
1766 },
1767 true, false));
1768 break;
1769 case CYLINDRICAL:
1770 pip.push_back(new OpMixScalarTimesDiv<CYLINDRICAL>(
1771 "P", "U",
1772 [](const double r, const double, const double) {
1773 return cylindrical(r);
1774 },
1775 true, false));
1776 break;
1777 default:
1778 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
1779 "Coordinate system not implemented");
1780 }
1781
1782 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P");
1783 pip.push_back(new OpDomainMassP("P", "P", [](double r, double, double) {
1784 return eps * cylindrical(r);
1785 }));
1786 }
1787
1789 };
1790
1791 auto get_block_name = [](auto name) {
1792 return boost::format("%s(.*)") % "WETTING_ANGLE";
1793 };
1794
1795 auto get_blocks = [&](auto &&name) {
1796 return mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
1797 std::regex(name.str()));
1798 };
1799
1800 auto set_boundary_rhs = [&](auto fe) {
1802 auto &pip = fe->getOpPtrVector();
1803
1805 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1806
1807 [&]() {
1808 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1810 return fe_parent;
1811 },
1812
1813 QUIET, Sev::noisy);
1814
1815 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
1816 pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
1817
1818 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
1819 pip.push_back(new OpCalculateScalarFieldValues("L", lambda_ptr));
1820 pip.push_back(new OpNormalConstrainRhs("L", u_ptr));
1821
1823 pip, mField, "L", {}, "INFLUX",
1824 [](double r, double, double) { return cylindrical(r); }, Sev::inform);
1825
1826 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
1827 pip.push_back(new OpNormalForceRhs("U", lambda_ptr));
1828
1829 auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
1830 if (wetting_block.size()) {
1831 // push operators to the side element which is called from op_bdy_side
1832 auto op_bdy_side =
1833 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
1834 op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
1835
1837 op_bdy_side->getOpPtrVector(), {H1});
1838
1840 op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
1842
1843 [&]() {
1844 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1845 new DomianParentEle(mField));
1847 fe_parent->getOpPtrVector(), {H1});
1848 return fe_parent;
1849 },
1850
1851 QUIET, Sev::noisy);
1852
1853 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1854 "H");
1855 op_bdy_side->getOpPtrVector().push_back(
1856 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1857 // push bdy side op
1858 pip.push_back(op_bdy_side);
1859
1860 // push operators for rhs wetting angle
1861 for (auto &b : wetting_block) {
1862 Range force_edges;
1863 std::vector<double> attr_vec;
1864 CHKERR b->getMeshsetIdEntitiesByDimension(
1865 mField.get_moab(), SPACE_DIM - 1, force_edges, true);
1866 b->getAttributes(attr_vec);
1867 if (attr_vec.size() != 1)
1868 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1869 "Should be one attribute");
1870 MOFEM_LOG("FS", Sev::inform) << "Wetting angle: " << attr_vec.front();
1871 // need to find the attributes and pass to operator
1872 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
1873 pip.push_back(new OpWettingAngleRhs(
1874 "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
1875 attr_vec.front()));
1876 }
1877 }
1878
1880 };
1881
1882 auto set_boundary_lhs = [&](auto fe) {
1884 auto &pip = fe->getOpPtrVector();
1885
1887 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1888
1889 [&]() {
1890 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1892 return fe_parent;
1893 },
1894
1895 QUIET, Sev::noisy);
1896
1897 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
1898 CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "U");
1899 pip.push_back(new OpNormalConstrainLhs("L", "U"));
1900
1901 auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
1902 if (wetting_block.size()) {
1903 auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
1904 auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
1905
1906 // push operators to the side element which is called from op_bdy_side
1907 auto op_bdy_side =
1908 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
1909 op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
1910
1912 op_bdy_side->getOpPtrVector(), {H1});
1913
1915 op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
1917
1918 [&]() {
1919 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1920 new DomianParentEle(mField));
1922 fe_parent->getOpPtrVector(), {H1});
1923 return fe_parent;
1924 },
1925
1926 QUIET, Sev::noisy);
1927
1928 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1929 "H");
1930 op_bdy_side->getOpPtrVector().push_back(
1931 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1932 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
1933 "H");
1934 op_bdy_side->getOpPtrVector().push_back(
1935 new OpLoopSideGetDataForSideEle("H", col_ind_ptr, col_diff_base_ptr));
1936
1937 // push bdy side op
1938 pip.push_back(op_bdy_side);
1939
1940 // push operators for lhs wetting angle
1941 for (auto &b : wetting_block) {
1942 Range force_edges;
1943 std::vector<double> attr_vec;
1944 CHKERR b->getMeshsetIdEntitiesByDimension(
1945 mField.get_moab(), SPACE_DIM - 1, force_edges, true);
1946 b->getAttributes(attr_vec);
1947 if (attr_vec.size() != 1)
1948 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1949 "Should be one attribute");
1950 MOFEM_LOG("FS", Sev::inform)
1951 << "wetting angle edges size " << force_edges.size();
1952
1953 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
1954 pip.push_back(new OpWettingAngleLhs(
1955 "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
1956 boost::make_shared<Range>(force_edges), attr_vec.front()));
1957 }
1958 }
1959
1961 };
1962
1963 auto *pip_mng = mField.getInterface<PipelineManager>();
1964
1965 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1966 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1967 CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
1968 CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
1969
1970 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
1971 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
1972 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
1973 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
1974
1975 pip_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
1976 pip_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
1977 pip_mng->getBoundaryLhsFE()->exeTestHook = test_bit_child;
1978 pip_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
1979
1981}
1982//! [Push operators to pip]
1983
1984/**
1985 * @brief Monitor solution
1986 *
1987 * This functions is called by TS kso at the end of each step. It is used
1988 */
1989struct Monitor : public FEMethod {
1991 SmartPetscObj<DM> dm, boost::shared_ptr<moab::Core> post_proc_mesh,
1992 boost::shared_ptr<PostProcEleDomainCont> post_proc,
1993 boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
1994 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
1995 p)
1996 : dM(dm), postProcMesh(post_proc_mesh), postProc(post_proc),
1997 postProcEdge(post_proc_edge), liftFE(p.first), liftVec(p.second) {}
2000
2001 MOFEM_LOG("FS", Sev::verbose) << "Monitor";
2002 constexpr int save_every_nth_step = 1;
2003 if (ts_step % save_every_nth_step == 0) {
2004 MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc";
2005 MoFEM::Interface *m_field_ptr;
2006 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
2007 auto post_proc_begin =
2008 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
2009 postProcMesh);
2010 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
2011 *m_field_ptr, postProcMesh);
2012 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
2014 this->getCacheWeakPtr());
2016 this->getCacheWeakPtr());
2017 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end->getFEMethod());
2018 CHKERR post_proc_end->writeFile(
2019 "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
2020 MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc done";
2021 }
2022
2023 liftVec->resize(SPACE_DIM, false);
2024 liftVec->clear();
2026 MPI_Allreduce(MPI_IN_PLACE, &(*liftVec)[0], SPACE_DIM, MPI_DOUBLE, MPI_SUM,
2027 MPI_COMM_WORLD);
2028 MOFEM_LOG("FS", Sev::inform)
2029 << "Step " << ts_step << " time " << ts_t
2030 << " lift vec x: " << (*liftVec)[0] << " y: " << (*liftVec)[1];
2031
2033 }
2034
2035private:
2037 boost::shared_ptr<moab::Core> postProcMesh;
2038 boost::shared_ptr<PostProcEleDomainCont> postProc;
2039 boost::shared_ptr<PostProcEleBdyCont> postProcEdge;
2040 boost::shared_ptr<BoundaryEle> liftFE;
2041 boost::shared_ptr<VectorDouble> liftVec;
2042};
2043
2044//! [Solve]
2047
2049
2050 auto simple = mField.getInterface<Simple>();
2051 auto pip_mng = mField.getInterface<PipelineManager>();
2052
2053 auto create_solver_dm = [&](auto dm) -> SmartPetscObj<DM> {
2054 DM subdm;
2055
2056 auto setup_subdm = [&](auto dm) {
2058 auto simple = mField.getInterface<Simple>();
2059 auto bit_mng = mField.getInterface<BitRefManager>();
2060 auto dm = simple->getDM();
2061 auto level_ents_ptr = boost::make_shared<Range>();
2062 CHKERR bit_mng->getEntitiesByRefLevel(
2063 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
2064 CHKERR DMCreate(mField.get_comm(), &subdm);
2065 CHKERR DMSetType(subdm, "DMMOFEM");
2066 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
2067 CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
2068 CHKERR DMMoFEMAddElement(subdm, simple->getBoundaryFEName());
2069 CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
2070 CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_FALSE);
2071 for (auto f : {"U", "P", "H", "G", "L"}) {
2072 CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
2073 CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
2074 }
2075 CHKERR DMSetUp(subdm);
2077 };
2078
2079 CHK_THROW_MESSAGE(setup_subdm(dm), "create subdm");
2080
2081 return SmartPetscObj<DM>(subdm);
2082 };
2083
2084 auto get_fe_post_proc = [&](auto post_proc_mesh) {
2085 auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
2086 return setParentDofs(
2087 fe, field, op, bit(get_skin_parent_bit()),
2088
2089 [&]() {
2090 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2091 new DomianParentEle(mField));
2092 return fe_parent;
2093 },
2094
2095 QUIET, Sev::noisy);
2096 };
2097
2098 auto post_proc_fe =
2099 boost::make_shared<PostProcEleDomainCont>(mField, post_proc_mesh);
2100 post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
2101 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2102 get_start_bit() + nb_levels - 1);
2103 };
2104
2105 auto u_ptr = boost::make_shared<MatrixDouble>();
2106 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2107 auto h_ptr = boost::make_shared<VectorDouble>();
2108 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2109 auto p_ptr = boost::make_shared<VectorDouble>();
2110 auto g_ptr = boost::make_shared<VectorDouble>();
2111 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2112
2114 post_proc_fe->getOpPtrVector(), {H1});
2115
2117 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2118
2119 [&]() {
2120 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2121 new DomianParentEle(mField));
2123 fe_parent->getOpPtrVector(), {H1});
2124 return fe_parent;
2125 },
2126
2127 QUIET, Sev::noisy);
2128
2129 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U");
2130 post_proc_fe->getOpPtrVector().push_back(
2132 post_proc_fe->getOpPtrVector().push_back(
2134 grad_u_ptr));
2135
2136 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H");
2137 post_proc_fe->getOpPtrVector().push_back(
2138 new OpCalculateScalarFieldValues("H", h_ptr));
2139 post_proc_fe->getOpPtrVector().push_back(
2140 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
2141
2142 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P");
2143 post_proc_fe->getOpPtrVector().push_back(
2144 new OpCalculateScalarFieldValues("P", p_ptr));
2145
2146 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G");
2147 post_proc_fe->getOpPtrVector().push_back(
2148 new OpCalculateScalarFieldValues("G", g_ptr));
2149 post_proc_fe->getOpPtrVector().push_back(
2150 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
2151
2153
2154 post_proc_fe->getOpPtrVector().push_back(
2155
2156 new OpPPMap(
2157 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2158
2159 {{"H", h_ptr}, {"P", p_ptr}, {"G", g_ptr}},
2160
2161 {{"U", u_ptr}, {"H_GRAD", grad_h_ptr}, {"G_GRAD", grad_g_ptr}},
2162
2163 {{"GRAD_U", grad_u_ptr}},
2164
2165 {}
2166
2167 )
2168
2169 );
2170
2171 return post_proc_fe;
2172 };
2173
2174 auto get_bdy_post_proc_fe = [&](auto post_proc_mesh) {
2175 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2176 return setParentDofs(
2177 fe, field, op, bit(get_skin_parent_bit()),
2178
2179 [&]() {
2180 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2181 new BoundaryParentEle(mField));
2182 return fe_parent;
2183 },
2184
2185 QUIET, Sev::noisy);
2186 };
2187
2188 auto post_proc_fe =
2189 boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
2190 post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
2191 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2192 get_start_bit() + nb_levels - 1);
2193 };
2194
2195 CHKERR setParentDofs(
2196 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2197
2198 [&]() {
2199 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2200 new BoundaryParentEle(mField));
2201 return fe_parent;
2202 },
2203
2204 QUIET, Sev::noisy);
2205
2206 struct OpGetNormal : public BoundaryEleOp {
2207
2208 OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
2209 boost::shared_ptr<MatrixDouble> n_ptr)
2210 : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), ptrL(l_ptr),
2211 ptrNormal(n_ptr) {}
2212
2213 MoFEMErrorCode doWork(int side, EntityType type,
2216 auto t_l = getFTensor0FromVec(*ptrL);
2217 auto t_n_fe = getFTensor1NormalsAtGaussPts();
2218 ptrNormal->resize(SPACE_DIM, getGaussPts().size2(), false);
2219 auto t_n = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
2220 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
2221 t_n(i) = t_n_fe(i) * t_l / std::sqrt(t_n_fe(i) * t_n_fe(i));
2222 ++t_n_fe;
2223 ++t_l;
2224 ++t_n;
2225 }
2227 };
2228
2229 protected:
2230 boost::shared_ptr<VectorDouble> ptrL;
2231 boost::shared_ptr<MatrixDouble> ptrNormal;
2232 };
2233
2234 auto u_ptr = boost::make_shared<MatrixDouble>();
2235 auto p_ptr = boost::make_shared<VectorDouble>();
2236 auto lambda_ptr = boost::make_shared<VectorDouble>();
2237 auto normal_l_ptr = boost::make_shared<MatrixDouble>();
2238
2239 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "U");
2240 post_proc_fe->getOpPtrVector().push_back(
2242
2243 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "L");
2244 post_proc_fe->getOpPtrVector().push_back(
2245 new OpCalculateScalarFieldValues("L", lambda_ptr));
2246 post_proc_fe->getOpPtrVector().push_back(
2247 new OpGetNormal(lambda_ptr, normal_l_ptr));
2248
2249 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "P");
2250 post_proc_fe->getOpPtrVector().push_back(
2251 new OpCalculateScalarFieldValues("P", p_ptr));
2252
2253 auto op_ptr = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
2254 op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
2258 auto op_ptr = static_cast<BoundaryEleOp *>(base_op_ptr);
2260 };
2261
2263
2264 post_proc_fe->getOpPtrVector().push_back(
2265
2266 new OpPPMap(post_proc_fe->getPostProcMesh(),
2267 post_proc_fe->getMapGaussPts(),
2268
2269 OpPPMap::DataMapVec{{"P", p_ptr}},
2270
2271 OpPPMap::DataMapMat{{"U", u_ptr}, {"L", normal_l_ptr}},
2272
2274
2276
2277 )
2278
2279 );
2280
2281 return post_proc_fe;
2282 };
2283
2284 auto get_lift_fe = [&]() {
2285 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2286 return setParentDofs(
2287 fe, field, op, bit(get_skin_parent_bit()),
2288
2289 [&]() {
2290 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2291 new BoundaryParentEle(mField));
2292 return fe_parent;
2293 },
2294
2295 QUIET, Sev::noisy);
2296 };
2297
2298 auto fe = boost::make_shared<BoundaryEle>(mField);
2299 fe->exeTestHook = [](FEMethod *fe_ptr) {
2300 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2301 get_start_bit() + nb_levels - 1);
2302 };
2303
2304 auto lift_ptr = boost::make_shared<VectorDouble>();
2305 auto p_ptr = boost::make_shared<VectorDouble>();
2306 auto ents_ptr = boost::make_shared<Range>();
2307
2308 CHKERR setParentDofs(
2309 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2310
2311 [&]() {
2312 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2313 new BoundaryParentEle(mField));
2314 return fe_parent;
2315 },
2316
2317 QUIET, Sev::noisy);
2318
2319 std::vector<const CubitMeshSets *> vec_ptr;
2320 CHKERR mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
2321 std::regex("LIFT"), vec_ptr);
2322 for (auto m_ptr : vec_ptr) {
2323 auto meshset = m_ptr->getMeshset();
2324 Range ents;
2325 CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
2326 ents, true);
2327 ents_ptr->merge(ents);
2328 }
2329
2330 MOFEM_LOG("FS", Sev::noisy) << "Lift ents " << (*ents_ptr);
2331
2332 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "P");
2333 fe->getOpPtrVector().push_back(
2334 new OpCalculateScalarFieldValues("L", p_ptr));
2335 fe->getOpPtrVector().push_back(
2336 new OpCalculateLift("L", p_ptr, lift_ptr, ents_ptr));
2337
2338 return std::make_pair(fe, lift_ptr);
2339 };
2340
2341 auto set_post_proc_monitor = [&](auto ts) {
2343 DM dm;
2344 CHKERR TSGetDM(ts, &dm);
2345 boost::shared_ptr<FEMethod> null_fe;
2346 auto post_proc_mesh = boost::make_shared<moab::Core>();
2347 auto monitor_ptr = boost::make_shared<Monitor>(
2348 SmartPetscObj<DM>(dm, true), post_proc_mesh,
2349 get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
2350 get_lift_fe());
2351 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
2352 null_fe, monitor_ptr);
2354 };
2355
2356 auto dm = simple->getDM();
2357 auto ts = createTS(mField.get_comm());
2358 CHKERR TSSetDM(ts, dm);
2359
2360 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
2361 tsPrePostProc = ts_pre_post_proc;
2362
2363 if (auto ptr = tsPrePostProc.lock()) {
2364
2365 ptr->fsRawPtr = this;
2366 ptr->solverSubDM = create_solver_dm(simple->getDM());
2367 ptr->globSol = createDMVector(dm);
2368 CHKERR DMoFEMMeshToLocalVector(dm, ptr->globSol, INSERT_VALUES,
2369 SCATTER_FORWARD);
2370 CHKERR VecAssemblyBegin(ptr->globSol);
2371 CHKERR VecAssemblyEnd(ptr->globSol);
2372
2373 auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
2374
2375 CHKERR set_post_proc_monitor(sub_ts);
2376
2377 // Add monitor to time solver
2378 double ftime = 1;
2379 CHKERR TSSetFromOptions(ts);
2380 CHKERR ptr->tsSetUp(ts);
2381 CHKERR TSSetUp(ts);
2382
2383 auto print_fields_in_section = [&]() {
2385 auto section = mField.getInterface<ISManager>()->sectionCreate(
2386 simple->getProblemName());
2387 PetscInt num_fields;
2388 CHKERR PetscSectionGetNumFields(section, &num_fields);
2389 for (int f = 0; f < num_fields; ++f) {
2390 const char *field_name;
2391 CHKERR PetscSectionGetFieldName(section, f, &field_name);
2392 MOFEM_LOG("FS", Sev::inform)
2393 << "Field " << f << " " << std::string(field_name);
2394 }
2396 };
2397
2398 CHKERR print_fields_in_section();
2399
2400 CHKERR TSSolve(ts, ptr->globSol);
2401 }
2402
2404}
2405
2406int main(int argc, char *argv[]) {
2407
2408#ifdef PYTHON_INIT_SURFACE
2409 Py_Initialize();
2410#endif
2411
2412 // Initialisation of MoFEM/PETSc and MOAB data structures
2413 const char param_file[] = "param_file.petsc";
2414 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
2415
2416 // Add logging channel for example
2417 auto core_log = logging::core::get();
2418 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FS"));
2419 LogManager::setLog("FS");
2420 MOFEM_LOG_TAG("FS", "free surface");
2421
2422 try {
2423
2424 //! [Register MoFEM discrete manager in PETSc]
2425 DMType dm_name = "DMMOFEM";
2426 CHKERR DMRegister_MoFEM(dm_name);
2427 //! [Register MoFEM discrete manager in PETSc
2428
2429 //! [Create MoAB]
2430 moab::Core mb_instance; ///< mesh database
2431 moab::Interface &moab = mb_instance; ///< mesh database interface
2432 //! [Create MoAB]
2433
2434 //! [Create MoFEM]
2435 MoFEM::Core core(moab); ///< finite element database
2436 MoFEM::Interface &m_field = core; ///< finite element database interface
2437 //! [Create MoFEM]
2438
2439 //! [FreeSurface]
2440 FreeSurface ex(m_field);
2441 CHKERR ex.runProblem();
2442 //! [FreeSurface]
2443 }
2445
2447
2448#ifdef PYTHON_INIT_SURFACE
2449 if (Py_FinalizeEx() < 0) {
2450 exit(120);
2451 }
2452#endif
2453}
2454
2455std::vector<Range>
2457
2458 auto &moab = mField.get_moab();
2459 auto bit_mng = mField.getInterface<BitRefManager>();
2460 auto comm_mng = mField.getInterface<CommInterface>();
2461
2462 Range vertices;
2463 CHK_THROW_MESSAGE(bit_mng->getEntitiesByTypeAndRefLevel(
2464 bit(0), BitRefLevel().set(), MBVERTEX, vertices),
2465 "can not get vertices on bit 0");
2466
2467 auto &dofs_mi = mField.get_dofs()->get<Unique_mi_tag>();
2468 auto field_bit_number = mField.get_field_bit_number("H");
2469
2470 Range plus_range, minus_range;
2471 std::vector<EntityHandle> plus, minus;
2472
2473 // get vertices on level 0 on plus and minus side
2474 for (auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
2475
2476 const auto f = p->first;
2477 const auto s = p->second;
2478
2479 // Lowest Dof UId for given field (field bit number) on entity f
2480 const auto lo_uid = DofEntity::getLoFieldEntityUId(field_bit_number, f);
2481 const auto hi_uid = DofEntity::getHiFieldEntityUId(field_bit_number, s);
2482 auto it = dofs_mi.lower_bound(lo_uid);
2483 const auto hi_it = dofs_mi.upper_bound(hi_uid);
2484
2485 plus.clear();
2486 minus.clear();
2487 plus.reserve(std::distance(it, hi_it));
2488 minus.reserve(std::distance(it, hi_it));
2489
2490 for (; it != hi_it; ++it) {
2491 const auto v = (*it)->getFieldData();
2492 if (v > 0)
2493 plus.push_back((*it)->getEnt());
2494 else
2495 minus.push_back((*it)->getEnt());
2496 }
2497
2498 plus_range.insert_list(plus.begin(), plus.end());
2499 minus_range.insert_list(minus.begin(), minus.end());
2500 }
2501
2502 MOFEM_LOG_CHANNEL("SYNC");
2503 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
2504 << "Plus range " << plus_range << endl;
2505 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
2506 << "Minus range " << minus_range << endl;
2508
2509 auto get_elems = [&](auto &ents, auto bit, auto mask) {
2510 Range adj;
2511 CHK_MOAB_THROW(moab.get_adjacencies(ents, SPACE_DIM, false, adj,
2512 moab::Interface::UNION),
2513 "can not get adjacencies");
2514 CHK_THROW_MESSAGE(bit_mng->filterEntitiesByRefLevel(bit, mask, adj),
2515 "can not filter elements with bit 0");
2516 return adj;
2517 };
2518
2519 CHKERR comm_mng->synchroniseEntities(plus_range);
2520 CHKERR comm_mng->synchroniseEntities(minus_range);
2521
2522 std::vector<Range> ele_plus(nb_levels), ele_minus(nb_levels);
2523 ele_plus[0] = get_elems(plus_range, bit(0), BitRefLevel().set());
2524 ele_minus[0] = get_elems(minus_range, bit(0), BitRefLevel().set());
2525 auto common = intersect(ele_plus[0], ele_minus[0]);
2526 ele_plus[0] = subtract(ele_plus[0], common);
2527 ele_minus[0] = subtract(ele_minus[0], common);
2528
2529 auto get_children = [&](auto &p, auto &c) {
2531 CHKERR bit_mng->updateRangeByChildren(p, c);
2532 c = c.subset_by_dimension(SPACE_DIM);
2534 };
2535
2536 for (auto l = 1; l != nb_levels; ++l) {
2537 CHK_THROW_MESSAGE(get_children(ele_plus[l - 1], ele_plus[l]),
2538 "get children");
2539 CHK_THROW_MESSAGE(get_children(ele_minus[l - 1], ele_minus[l]),
2540 "get children");
2541 }
2542
2543 auto get_level = [&](auto &p, auto &m, auto z, auto bit, auto mask) {
2544 Range l;
2546 bit_mng->getEntitiesByDimAndRefLevel(bit, mask, SPACE_DIM, l),
2547 "can not get vertices on bit");
2548 l = subtract(l, p);
2549 l = subtract(l, m);
2550 for (auto f = 0; f != z; ++f) {
2551 Range conn;
2552 CHK_MOAB_THROW(moab.get_connectivity(l, conn, true), "");
2553 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(conn);
2554 l = get_elems(conn, bit, mask);
2555 }
2556 return l;
2557 };
2558
2559 std::vector<Range> vec_levels(nb_levels);
2560 for (auto l = nb_levels - 1; l >= 0; --l) {
2561 vec_levels[l] = get_level(ele_plus[l], ele_minus[l], 2 * overlap, bit(l),
2562 BitRefLevel().set());
2563 }
2564
2565 if constexpr (debug) {
2566 if (mField.get_comm_size() == 1) {
2567 for (auto l = 0; l != nb_levels; ++l) {
2568 std::string name = (boost::format("out_r%d.h5m") % l).str();
2569 CHK_THROW_MESSAGE(save_range(mField.get_moab(), name, vec_levels[l]),
2570 "save mesh");
2571 }
2572 }
2573 }
2574
2575 return vec_levels;
2576}
2577
2580
2581 auto simple = mField.getInterface<Simple>();
2582 auto bit_mng = mField.getInterface<BitRefManager>();
2583 auto prb_mng = mField.getInterface<ProblemsManager>();
2584
2585 BitRefLevel start_mask;
2586 for (auto s = 0; s != get_start_bit(); ++s)
2587 start_mask[s] = true;
2588
2589 // store prev_level
2590 Range prev_level;
2591 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_current_bit()),
2592 BitRefLevel().set(), prev_level);
2593 Range prev_level_skin;
2594 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_parent_bit()),
2595 BitRefLevel().set(), prev_level_skin);
2596 // reset bit ref levels
2597 CHKERR bit_mng->lambdaBitRefLevel(
2598 [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
2599 CHKERR bit_mng->setNthBitRefLevel(prev_level, get_projection_bit(), true);
2600 CHKERR bit_mng->setNthBitRefLevel(prev_level_skin, get_skin_projection_bit(),
2601 true);
2602
2603 // set refinement levels
2604 auto set_levels = [&](auto &&
2605 vec_levels /*entities are refined on each level*/) {
2607
2608 // start with zero level, which is the coarsest mesh
2609 Range level0;
2610 CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
2611 CHKERR bit_mng->setNthBitRefLevel(level0, get_start_bit(), true);
2612
2613 // get lower dimension entities
2614 auto get_adj = [&](auto ents) {
2615 Range conn;
2616 CHK_MOAB_THROW(mField.get_moab().get_connectivity(ents, conn, true),
2617 "get conn");
2618 for (auto d = 1; d != SPACE_DIM; ++d) {
2619 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
2620 ents.subset_by_dimension(SPACE_DIM), d, false, conn,
2621 moab::Interface::UNION),
2622 "get adj");
2623 }
2624 ents.merge(conn);
2625 return ents;
2626 };
2627
2628 // set bit levels
2629 for (auto l = 1; l != nb_levels; ++l) {
2630 Range level_prev;
2631 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_start_bit() + l - 1),
2632 BitRefLevel().set(),
2633 SPACE_DIM, level_prev);
2634 Range parents;
2635 CHKERR bit_mng->updateRangeByParent(vec_levels[l], parents);
2636 // subtract entities from previous level, which are refined, so should be
2637 // not there
2638 level_prev = subtract(level_prev, parents);
2639 // and instead add their children
2640 level_prev.merge(vec_levels[l]);
2641 // set bit to each level
2642 CHKERR bit_mng->setNthBitRefLevel(level_prev, get_start_bit() + l, true);
2643 }
2644
2645 // set bit levels to lower dimension entities
2646 for (auto l = 1; l != nb_levels; ++l) {
2647 Range level;
2648 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2649 bit(get_start_bit() + l), BitRefLevel().set(), SPACE_DIM, level);
2650 level = get_adj(level);
2651 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(level);
2652 CHKERR bit_mng->setNthBitRefLevel(level, get_start_bit() + l, true);
2653 }
2654
2656 };
2657
2658 // resolve skin between refined levels
2659 auto set_skins = [&]() {
2661
2662 moab::Skinner skinner(&mField.get_moab());
2663 ParallelComm *pcomm =
2664 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2665
2666 // get skin of bit level
2667 auto get_bit_skin = [&](BitRefLevel bit, BitRefLevel mask) {
2668 Range bit_ents;
2671 bit, mask, SPACE_DIM, bit_ents),
2672 "can't get bit level");
2673 Range bit_skin;
2674 CHK_MOAB_THROW(skinner.find_skin(0, bit_ents, false, bit_skin),
2675 "can't get skin");
2676 return bit_skin;
2677 };
2678
2679 auto get_level_skin = [&]() {
2680 Range skin;
2681 BitRefLevel bit_prev;
2682 for (auto l = 1; l != nb_levels; ++l) {
2683 auto skin_level_mesh = get_bit_skin(bit(l), BitRefLevel().set());
2684 // filter (remove) all entities which are on partitions boarders
2685 CHKERR pcomm->filter_pstatus(skin_level_mesh,
2686 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2687 PSTATUS_NOT, -1, nullptr);
2688 auto skin_level =
2689 get_bit_skin(bit(get_start_bit() + l), BitRefLevel().set());
2690 skin_level = subtract(skin_level,
2691 skin_level_mesh); // get only internal skins, not
2692 // on the body boundary
2693 // get lower dimension adjacencies (FIXME: add edges if 3D)
2694 Range skin_level_verts;
2695 CHKERR mField.get_moab().get_connectivity(skin_level, skin_level_verts,
2696 true);
2697 skin_level.merge(skin_level_verts);
2698
2699 // remove previous level
2700 bit_prev.set(l - 1);
2701 Range level_prev;
2702 CHKERR bit_mng->getEntitiesByRefLevel(bit_prev, BitRefLevel().set(),
2703 level_prev);
2704 skin.merge(subtract(skin_level, level_prev));
2705 }
2706
2707 return skin;
2708 };
2709
2710 auto resolve_shared = [&](auto &&skin) {
2711 Range tmp_skin = skin;
2712
2713 map<int, Range> map_procs;
2714 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2715 tmp_skin, &map_procs);
2716
2717 Range from_other_procs; // entities which also exist on other processors
2718 for (auto &m : map_procs) {
2719 if (m.first != mField.get_comm_rank()) {
2720 from_other_procs.merge(m.second);
2721 }
2722 }
2723
2724 auto common = intersect(
2725 skin, from_other_procs); // entities which are on internal skin
2726 skin.merge(from_other_procs);
2727
2728 // entities which are on processors boards, and several processors are not
2729 // true skin.
2730 if (!common.empty()) {
2731 // skin is internal exist on other procs
2732 skin = subtract(skin, common);
2733 }
2734
2735 return skin;
2736 };
2737
2738 // get parents of entities
2739 auto get_parent_level_skin = [&](auto skin) {
2740 Range skin_parents;
2741 CHKERR bit_mng->updateRangeByParent(
2742 skin.subset_by_dimension(SPACE_DIM - 1), skin_parents);
2743 Range skin_parent_verts;
2744 CHKERR mField.get_moab().get_connectivity(skin_parents, skin_parent_verts,
2745 true);
2746 skin_parents.merge(skin_parent_verts);
2747 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2748 skin_parents);
2749 return skin_parents;
2750 };
2751
2752 auto child_skin = resolve_shared(get_level_skin());
2753 auto parent_skin = get_parent_level_skin(child_skin);
2754
2755 child_skin = subtract(child_skin, parent_skin);
2756 CHKERR bit_mng->setNthBitRefLevel(child_skin, get_skin_child_bit(), true);
2757 CHKERR bit_mng->setNthBitRefLevel(parent_skin, get_skin_parent_bit(), true);
2758
2760 };
2761
2762 // take last level, remove childs on boarder, and set bit
2763 auto set_current = [&]() {
2765 Range last_level;
2766 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_start_bit() + nb_levels - 1),
2767 BitRefLevel().set(), last_level);
2768 Range skin_child;
2769 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_child_bit()),
2770 BitRefLevel().set(), skin_child);
2771
2772 last_level = subtract(last_level, skin_child);
2773 CHKERR bit_mng->setNthBitRefLevel(last_level, get_current_bit(), true);
2775 };
2776
2777 // set bits to levels
2778 CHKERR set_levels(findEntitiesCrossedByPhaseInterface(overlap));
2779 // set bits to skin
2780 CHKERR set_skins();
2781 // set current level bit
2782 CHKERR set_current();
2783
2784 if constexpr (debug) {
2785 if (mField.get_comm_size() == 1) {
2786 for (auto l = 0; l != nb_levels; ++l) {
2787 std::string name = (boost::format("out_level%d.h5m") % l).str();
2788 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_start_bit() + l),
2789 BitRefLevel().set(), name.c_str(), "MOAB",
2790 "PARALLEL=WRITE_PART");
2791 }
2792 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_current_bit()),
2793 BitRefLevel().set(), "current_bit.h5m",
2794 "MOAB", "PARALLEL=WRITE_PART");
2795 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_projection_bit()),
2796 BitRefLevel().set(), "projection_bit.h5m",
2797 "MOAB", "PARALLEL=WRITE_PART");
2798
2799 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_child_bit()),
2800 BitRefLevel().set(), "skin_child_bit.h5m",
2801 "MOAB", "PARALLEL=WRITE_PART");
2802 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_parent_bit()),
2803 BitRefLevel().set(), "skin_parent_bit.h5m",
2804 "MOAB", "PARALLEL=WRITE_PART");
2805 }
2806 }
2807
2809};
2810
2812 boost::shared_ptr<FEMethod> fe_top, std::string field_name,
2813 ForcesAndSourcesCore::UserDataOperator::OpType op,
2814 BitRefLevel child_ent_bit,
2815 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
2816 int verbosity, LogManager::SeverityLevel sev) {
2818
2819 /**
2820 * @brief Collect data from parent elements to child
2821 */
2822 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
2823 add_parent_level =
2824 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
2825 // Evaluate if not last parent element
2826 if (level > 0) {
2827
2828 // Create domain parent FE
2829 auto fe_ptr_current = get_elem();
2830
2831 // Call next level
2832 add_parent_level(
2833 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
2834 fe_ptr_current),
2835 level - 1);
2836
2837 // Add data to curent fe level
2838 if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
2839
2840 // Only base
2841 parent_fe_pt->getOpPtrVector().push_back(
2842
2844
2845 H1, op, fe_ptr_current,
2846
2847 BitRefLevel().set(), BitRefLevel().set(),
2848
2849 child_ent_bit, BitRefLevel().set(),
2850
2851 verbosity, sev));
2852
2853 } else {
2854
2855 // Filed data
2856 parent_fe_pt->getOpPtrVector().push_back(
2857
2859
2860 field_name, op, fe_ptr_current,
2861
2862 BitRefLevel().set(), BitRefLevel().set(),
2863
2864 child_ent_bit, BitRefLevel().set(),
2865
2866 verbosity, sev));
2867 }
2868 }
2869 };
2870
2871 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
2872 nb_levels);
2873
2875}
2876
2879
2880 if (auto ptr = tsPrePostProc.lock()) {
2881
2882 /**
2883 * @brief cut-off values at nodes, i.e. abs("H") <= 1
2884 *
2885 */
2886 auto cut_off_dofs = [&]() {
2888
2889 auto &m_field = ptr->fsRawPtr->mField;
2890
2891 Range current_verts;
2892 auto bit_mng = m_field.getInterface<BitRefManager>();
2893 CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
2894 bit(get_current_bit()), BitRefLevel().set(), MBVERTEX, current_verts);
2895
2896 auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
2898 for (auto &h : ent_ptr->getEntFieldData()) {
2899 h = cut_off(h);
2900 }
2902 };
2903
2904 auto field_blas = m_field.getInterface<FieldBlas>();
2905 CHKERR field_blas->fieldLambdaOnEntities(cut_off_verts, "H",
2906 &current_verts);
2908 };
2909
2910 CHKERR cut_off_dofs();
2911 }
2912
2913 if (auto ptr = tsPrePostProc.lock()) {
2914 MOFEM_LOG("FS", Sev::inform) << "Run step pre proc";
2915
2916 auto &m_field = ptr->fsRawPtr->mField;
2917 auto simple = m_field.getInterface<Simple>();
2918 auto bit_mng = m_field.getInterface<BitRefManager>();
2919 auto bc_mng = m_field.getInterface<BcManager>();
2920 auto field_blas = m_field.getInterface<FieldBlas>();
2921 auto opt = m_field.getInterface<OperatorsTester>();
2922 auto pip_mng = m_field.getInterface<PipelineManager>();
2923 auto prb_mng = m_field.getInterface<ProblemsManager>();
2924
2925 // get vector norm
2926 auto get_norm = [&](auto x) {
2927 double nrm;
2928 CHKERR VecNorm(x, NORM_2, &nrm);
2929 return nrm;
2930 };
2931
2932 // refine problem and project data, including theta data
2933 auto refine_problem = [&]() {
2935 MOFEM_LOG("FS", Sev::inform) << "Refine problem";
2936 CHKERR ptr->fsRawPtr->refineMesh(refine_overlap);
2937 CHKERR ptr->fsRawPtr->projectData();
2939 };
2940
2941 // set new jacobin operator, since problem and thus tangent matrix size has
2942 // changed
2943 auto set_jacobian_operators = [&]() {
2945 ptr->subB = createDMMatrix(ptr->solverSubDM);
2946 CHKERR KSPReset(ptr->subKSP);
2948 };
2949
2950 // set new solution
2951 auto set_solution = [&]() {
2953 MOFEM_LOG("FS", Sev::inform) << "Set solution";
2954
2955 PetscObjectState state;
2956
2957 // Record the state, and set it again. This is to fool PETSc that solution
2958 // vector is not updated. Otherwise PETSc will treat every step as a first
2959 // step.
2960
2961 // globSol is updated as result mesh refinement - this is not really set
2962 // a new solution.
2963
2964 CHKERR PetscObjectStateGet(getPetscObject(ptr->globSol.get()), &state);
2965 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), ptr->globSol,
2966 INSERT_VALUES, SCATTER_FORWARD);
2967 CHKERR PetscObjectStateSet(getPetscObject(ptr->globSol.get()), state);
2968 MOFEM_LOG("FS", Sev::verbose)
2969 << "Set solution, vector norm " << get_norm(ptr->globSol);
2971 };
2972
2973 PetscBool is_theta;
2974 PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
2975 if (is_theta) {
2976
2977 CHKERR refine_problem(); // refine problem
2978 CHKERR set_jacobian_operators(); // set new jacobian
2979 CHKERR set_solution(); // set solution
2980
2981 } else {
2982 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
2983 "Sorry, only TSTheta handling is implemented");
2984 }
2985
2986 // Need barriers, somme functions in TS solver need are called collectively
2987 // and requite the same state of variables
2988 PetscBarrier((PetscObject)ts);
2989
2990 MOFEM_LOG_CHANNEL("SYNC");
2991 MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PreProc done";
2992 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
2993 }
2994
2996}
2997
2999 if (auto ptr = tsPrePostProc.lock()) {
3000 auto &m_field = ptr->fsRawPtr->mField;
3001 MOFEM_LOG_CHANNEL("SYNC");
3002 MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PostProc done";
3003 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
3004 }
3005 return 0;
3006}
3007
3008MoFEMErrorCode TSPrePostProc::tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t,
3009 Vec f, void *ctx) {
3011 if (auto ptr = tsPrePostProc.lock()) {
3012 auto sub_u = ptr->getSubVector();
3013 auto sub_u_t = vectorDuplicate(sub_u);
3014 auto sub_f = vectorDuplicate(sub_u);
3015 auto scatter = ptr->getScatter(sub_u, u, R);
3016
3017 auto apply_scatter_and_update = [&](auto x, auto sub_x) {
3019 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3020 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3021 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3022 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3024 };
3025
3026 CHKERR apply_scatter_and_update(u, sub_u);
3027 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3028
3029 CHKERR TsSetIFunction(ts, t, sub_u, sub_u_t, sub_f, ptr->tsCtxPtr.get());
3030 CHKERR VecScatterBegin(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3031 CHKERR VecScatterEnd(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3032 }
3034}
3035
3036MoFEMErrorCode TSPrePostProc::tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
3037 PetscReal a, Mat A, Mat B,
3038 void *ctx) {
3040 if (auto ptr = tsPrePostProc.lock()) {
3041 auto sub_u = ptr->getSubVector();
3042 auto sub_u_t = vectorDuplicate(sub_u);
3043 auto scatter = ptr->getScatter(sub_u, u, R);
3044
3045 auto apply_scatter_and_update = [&](auto x, auto sub_x) {
3047 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3048 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3049 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3050 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3052 };
3053
3054 CHKERR apply_scatter_and_update(u, sub_u);
3055 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3056
3057 CHKERR TsSetIJacobian(ts, t, sub_u, sub_u_t, a, ptr->subB, ptr->subB,
3058 ptr->tsCtxPtr.get());
3059 }
3061}
3062
3063MoFEMErrorCode TSPrePostProc::tsMonitor(TS ts, PetscInt step, PetscReal t,
3064 Vec u, void *ctx) {
3066 if (auto ptr = tsPrePostProc.lock()) {
3067 auto get_norm = [&](auto x) {
3068 double nrm;
3069 CHKERR VecNorm(x, NORM_2, &nrm);
3070 return nrm;
3071 };
3072
3073 auto sub_u = ptr->getSubVector();
3074 auto scatter = ptr->getScatter(sub_u, u, R);
3075 CHKERR VecScatterBegin(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3076 CHKERR VecScatterEnd(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3077 CHKERR VecGhostUpdateBegin(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3078 CHKERR VecGhostUpdateEnd(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3079
3080 MOFEM_LOG("FS", Sev::verbose)
3081 << "u norm " << get_norm(u) << " u sub nom " << get_norm(sub_u);
3082
3083 CHKERR TsMonitorSet(ts, step, t, sub_u, ptr->tsCtxPtr.get());
3084 }
3086}
3087
3090 if (auto ptr = tsPrePostProc.lock()) {
3091 MOFEM_LOG("FS", Sev::verbose) << "SetUP sub PC";
3092 ptr->subKSP = createKSP(ptr->fsRawPtr->mField.get_comm());
3093 CHKERR KSPSetFromOptions(ptr->subKSP);
3094 }
3096};
3097
3098MoFEMErrorCode TSPrePostProc::pcApply(PC pc, Vec pc_f, Vec pc_x) {
3100 if (auto ptr = tsPrePostProc.lock()) {
3101 auto sub_x = ptr->getSubVector();
3102 auto sub_f = vectorDuplicate(sub_x);
3103 auto scatter = ptr->getScatter(sub_x, pc_x, R);
3104 CHKERR VecScatterBegin(scatter, pc_f, sub_f, INSERT_VALUES,
3105 SCATTER_REVERSE);
3106 CHKERR VecScatterEnd(scatter, pc_f, sub_f, INSERT_VALUES, SCATTER_REVERSE);
3107 CHKERR KSPSetOperators(ptr->subKSP, ptr->subB, ptr->subB);
3108 MOFEM_LOG("FS", Sev::verbose) << "PCShell solve";
3109 CHKERR KSPSolve(ptr->subKSP, sub_f, sub_x);
3110 MOFEM_LOG("FS", Sev::verbose) << "PCShell solve <- done";
3111 CHKERR VecScatterBegin(scatter, sub_x, pc_x, INSERT_VALUES,
3112 SCATTER_FORWARD);
3113 CHKERR VecScatterEnd(scatter, sub_x, pc_x, INSERT_VALUES, SCATTER_FORWARD);
3114 }
3116};
3117
3119 if (auto ptr = tsPrePostProc.lock()) {
3120 auto prb_ptr = ptr->fsRawPtr->mField.get_problem("SUB_SOLVER");
3121 if (auto sub_data = prb_ptr->getSubData()) {
3122 auto is = sub_data->getSmartColIs();
3123 VecScatter s;
3124 if (fr == R) {
3125 CHK_THROW_MESSAGE(VecScatterCreate(x, PETSC_NULL, y, is, &s),
3126 "crate scatter");
3127 } else {
3128 CHK_THROW_MESSAGE(VecScatterCreate(x, is, y, PETSC_NULL, &s),
3129 "crate scatter");
3130 }
3131 return SmartPetscObj<VecScatter>(s);
3132 }
3133 }
3136}
3137
3140}
3141
3143
3144 auto &m_field = fsRawPtr->mField;
3145 auto simple = m_field.getInterface<Simple>();
3146
3148
3149 auto dm = simple->getDM();
3150
3152 CHKERR TSSetIFunction(ts, globRes, tsSetIFunction, nullptr);
3153 CHKERR TSSetIJacobian(ts, PETSC_NULL, PETSC_NULL, tsSetIJacobian, nullptr);
3154 CHKERR TSMonitorSet(ts, tsMonitor, fsRawPtr, PETSC_NULL);
3155
3156 SNES snes;
3157 CHKERR TSGetSNES(ts, &snes);
3158 auto snes_ctx_ptr = getDMSnesCtx(dm);
3159
3160 auto set_section_monitor = [&](auto snes) {
3162 CHKERR SNESMonitorSet(snes,
3163 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
3164 void *))MoFEMSNESMonitorFields,
3165 (void *)(snes_ctx_ptr.get()), nullptr);
3167 };
3168
3169 CHKERR set_section_monitor(snes);
3170
3171 auto ksp = createKSP(m_field.get_comm());
3172 CHKERR KSPSetType(ksp, KSPPREONLY); // Run KSP internally in ShellPC
3173 auto sub_pc = createPC(fsRawPtr->mField.get_comm());
3174 CHKERR PCSetType(sub_pc, PCSHELL);
3175 CHKERR PCShellSetSetUp(sub_pc, pcSetup);
3176 CHKERR PCShellSetApply(sub_pc, pcApply);
3177 CHKERR KSPSetPC(ksp, sub_pc);
3178 CHKERR SNESSetKSP(snes, ksp);
3179
3182
3183 CHKERR TSSetPreStep(ts, tsPreProc);
3184 CHKERR TSSetPostStep(ts, tsPostProc);
3185
3187}
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
constexpr int SPACE_DIM
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
@ LAST_COORDINATE_SYSTEM
Definition: definitions.h:119
@ CYLINDRICAL
Definition: definitions.h:117
@ CARTESIAN
Definition: definitions.h:115
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
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
double kappa
auto marker
auto bubble_device
double mu_diff
auto save_range
auto get_M2
constexpr int U_FIELD_DIM
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
FTensor::Index< 'j', SPACE_DIM > j
auto get_M2_dh
double tol
@ R
@ F
auto integration_rule
static char help[]
auto get_M_dh
auto kernel_eye
double mu_m
double lambda
surface tension
FTensor::Index< 'k', SPACE_DIM > k
OpDomainMassH OpDomainMassG
auto get_M3_dh
FTensor::Index< 'i', SPACE_DIM > i
double rho_diff
auto init_h
Initialisation function.
constexpr auto t_kd
auto d_phase_function_h
auto get_fe_bit
constexpr int SPACE_DIM
double rho_ave
FTensor::Index< 'l', SPACE_DIM > l
int coord_type
double W
auto kernel_oscillation
double eta2
auto get_skin_parent_bit
auto get_M0
constexpr bool debug
auto get_f_dh
auto get_M3
double rho_p
double mu_p
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
double h
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
auto get_M
double mu_ave
double eps
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
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
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
int order
approximation order
SideEle::UserDataOperator SideOp
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
double a0
double rho_m
auto get_current_bit
dofs bit used to do calculations
BoundaryEle::UserDataOperator BoundaryEleOp
double eta
auto bit
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.
FTensor::Index< 'i', SPACE_DIM > i
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 PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
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)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscObject getPetscObject(T obj)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1767
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
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
constexpr double g
MoFEM::Interface & mField
Definition: plastic.cpp:229
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:25
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
base operator to do operations at Gauss Pt. level
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
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
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
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.
[Push operators to pipeline]
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
Example * fsRawPtr
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
constexpr CoordinateTypes COORD_TYPE