v0.15.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#undef PYTHON_INIT_SURFACE
20
21#ifdef PYTHON_INIT_SURFACE
22#include <boost/python.hpp>
23#include <boost/python/def.hpp>
24namespace bp = boost::python;
25
26struct SurfacePython {
27 SurfacePython() = default;
28 virtual ~SurfacePython() = default;
29
30 MoFEMErrorCode surfaceInit(const std::string py_file) {
32 try {
33
34 // create main module
35 auto main_module = bp::import("__main__");
36 mainNamespace = main_module.attr("__dict__");
37 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
38 // create a reference to python function
39 surfaceFun = mainNamespace["surface"];
40
41 } catch (bp::error_already_set const &) {
42 // print all other errors to stderr
43 PyErr_Print();
45 }
47 };
48
49 MoFEMErrorCode evalSurface(
50
51 double x, double y, double z, double eta, double &s
52
53 ) {
55 try {
56
57 // call python function
58 s = bp::extract<double>(surfaceFun(x, y, z, eta));
59
60 } catch (bp::error_already_set const &) {
61 // print all other errors to stderr
62 PyErr_Print();
64 }
66 }
67
68private:
69 bp::object mainNamespace;
70 bp::object surfaceFun;
71};
72
73static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
74
75#endif
76
77int coord_type = EXECUTABLE_COORD_TYPE;
78
79constexpr int BASE_DIM = 1;
80constexpr int SPACE_DIM = 2;
81constexpr int U_FIELD_DIM = SPACE_DIM;
82
83constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
84constexpr IntegrationType I =
85 IntegrationType::GAUSS; //< selected integration type
86
87template <int DIM>
89
92using DomainEleOp = DomainEle::UserDataOperator;
95using BoundaryEleOp = BoundaryEle::UserDataOperator;
97using SideOp = SideEle::UserDataOperator;
98
102
104
108
117
124
125template <CoordinateTypes COORD_TYPE>
128
129// Flux is applied by Lagrange Multiplie
133
138
140
141// mesh refinement
142int order = 3; ///< approximation order
143int nb_levels = 4; //< number of refinement levels
144int refine_overlap = 4; //< mesh overlap while refine
145
146constexpr bool debug = true;
147
148auto get_start_bit = []() {
149 return nb_levels + 1;
150}; //< first refinement level for computational mesh
151auto get_current_bit = []() {
152 return 2 * get_start_bit() + 1;
153}; ///< dofs bit used to do calculations
154auto get_skin_parent_bit = []() { return 2 * get_start_bit() + 2; };
155auto get_skin_child_bit = []() { return 2 * get_start_bit() + 3; };
156auto get_projection_bit = []() { return 2 * get_start_bit() + 4; };
157auto get_skin_projection_bit = []() { return 2 * get_start_bit() + 5; };
158
159// FIXME: Set parameters from command line
160
161// Physical parameters
162double a0 = 980;
163double rho_m = 0.998;
164double mu_m = 0.010101 * 1e2;
165double rho_p = 0.0012;
166double mu_p = 0.000182 * 1e2;
167double lambda = 73; ///< surface tension
168double W = 0.25;
169
170// Model parameters
171double h = 0.03; // mesh size
172double eta = h;
173double eta2 = eta * eta;
174
175// Numerical parameters
176double md = 1e-2; // mobility
177double eps = 1e-12;
178double tol = std::numeric_limits<float>::epsilon();
179
180double rho_ave = (rho_p + rho_m) / 2;
181double rho_diff = (rho_p - rho_m) / 2;
182double mu_ave = (mu_p + mu_m) / 2;
183double mu_diff = (mu_p - mu_m) / 2;
184
185double kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
186
187auto integration_rule = [](int, int, int) { return 2 * order + 1; };
188
189auto cylindrical = [](const double r) {
190 if (coord_type == CYLINDRICAL)
191 return 2 * M_PI * r;
192 else
193 return 1.;
194};
195
196auto wetting_angle_sub_stepping = [](auto ts_step) {
197 constexpr int sub_stepping = 16;
198 return std::min(1., static_cast<double>(ts_step) / sub_stepping);
199};
200
201// cut off function
202
203auto my_max = [](const double x) { return (x - 1 + std::abs(x + 1)) / 2; };
204auto my_min = [](const double x) { return (x + 1 - std::abs(x - 1)) / 2; };
205auto cut_off = [](const double h) { return my_max(my_min(h)); };
206auto d_cut_off = [](const double h) {
207 if (h >= -1 && h < 1)
208 return 1.;
209 else
210 return 0.;
211};
212
213auto phase_function = [](const double h, const double diff, const double ave) {
214 return diff * h + ave;
215};
216
217auto d_phase_function_h = [](const double h, const double diff) {
218 return diff;
219};
220
221// order function (free energy)
222
223auto get_f = [](const double h) { return 4 * W * h * (h * h - 1); };
224auto get_f_dh = [](const double h) { return 4 * W * (3 * h * h - 1); };
225
226auto get_M0 = [](auto h) { return md; };
227auto get_M0_dh = [](auto h) { return 0; };
228
229auto get_M2 = [](auto h) {
230 return md * (1 - h * h);
231};
232
233auto get_M2_dh = [](auto h) { return -md * 2 * h; };
234
235auto get_M3 = [](auto h) {
236 const double h2 = h * h;
237 const double h3 = h2 * h;
238 if (h >= 0)
239 return md * (2 * h3 - 3 * h2 + 1);
240 else
241 return md * (-2 * h3 - 3 * h2 + 1);
242};
243
244auto get_M3_dh = [](auto h) {
245 if (h >= 0)
246 return md * (6 * h * (h - 1));
247 else
248 return md * (-6 * h * (h + 1));
249};
250
251auto get_M = [](auto h) { return get_M0(h); };
252auto get_M_dh = [](auto h) { return get_M0_dh(h); };
253
254auto get_D = [](const double A) {
256 t_D(i, j, k, l) = A * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
257 return t_D;
258};
259
260// some examples of initialisation functions
261
262auto kernel_oscillation = [](double r, double y, double) {
263 constexpr int n = 3;
264 constexpr double R = 0.0125;
265 constexpr double A = R * 0.2;
266 const double theta = atan2(r, y);
267 const double w = R + A * cos(n * theta);
268 const double d = std::sqrt(r * r + y * y);
269 return tanh((w - d) / (eta * std::sqrt(2)));
270};
271
272auto kernel_eye = [](double r, double y, double) {
273 constexpr double y0 = 0.5;
274 constexpr double R = 0.4;
275 y -= y0;
276 const double d = std::sqrt(r * r + y * y);
277 return tanh((R - d) / (eta * std::sqrt(2)));
278};
279
280auto capillary_tube = [](double x, double y, double z) {
281 constexpr double water_height = 0.;
282 return tanh((water_height - y) / (eta * std::sqrt(2)));
283 ;
284};
285
286auto bubble_device = [](double x, double y, double z) {
287 return -tanh((-0.039 - x) / (eta * std::sqrt(2)));
288};
289
290/**
291 * @brief Initialisation function
292 *
293 * @note If UMs are compiled with Python to initialise phase field "H"
294 * surface.py function is used, which has to be present in execution folder.
295 *
296 */
297auto init_h = [](double r, double y, double theta) {
298#ifdef PYTHON_INIT_SURFACE
299 double s = 1;
300 if (auto ptr = surfacePythonWeakPtr.lock()) {
301 CHK_THROW_MESSAGE(ptr->evalSurface(r, y, theta, eta, s),
302 "error eval python");
303 }
304 return s;
305#else
306 // return bubble_device(r, y, theta);
307 return capillary_tube(r, y, theta);
308 // return kernel_eye(r, y, theta);
309#endif
310};
311
312auto wetting_angle = [](double water_level) { return water_level; };
313
314auto bit = [](auto b) { return BitRefLevel().set(b); };
315auto marker = [](auto b) { return BitRefLevel().set(BITREFLEVEL_SIZE - b); };
316auto get_fe_bit = [](FEMethod *fe_ptr) {
317 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
318};
319
320auto get_global_size = [](int l_size) {
321 int g_size;
322 MPI_Allreduce(&l_size, &g_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
323 return g_size;
324};
325
326auto save_range = [](moab::Interface &moab, const std::string name,
327 const Range r) {
329 if (get_global_size(r.size())) {
330 auto out_meshset = get_temp_meshset_ptr(moab);
331 CHKERR moab.add_entities(*out_meshset, r);
332 CHKERR moab.write_file(name.c_str(), "MOAB", "PARALLEL=WRITE_PART",
333 out_meshset->get_ptr(), 1);
334 }
336};
337
338/**
339 * @brief get entities of dofs in the problem - used for debugging
340 *
341 */
342auto get_dofs_ents_by_field_name = [](auto dm, auto field_name) {
343 auto prb_ptr = getProblemPtr(dm);
344 std::vector<EntityHandle> ents_vec;
345
346 MoFEM::Interface *m_field_ptr;
347 CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
348
349 auto bit_number = m_field_ptr->get_field_bit_number(field_name);
350 auto dofs = prb_ptr->numeredRowDofsPtr;
351 auto lo_it = dofs->lower_bound(FieldEntity::getLoBitNumberUId(bit_number));
352 auto hi_it = dofs->upper_bound(FieldEntity::getHiBitNumberUId(bit_number));
353 ents_vec.reserve(std::distance(lo_it, hi_it));
354
355 for (; lo_it != hi_it; ++lo_it) {
356 ents_vec.push_back((*lo_it)->getEnt());
357 }
358
359 std::sort(ents_vec.begin(), ents_vec.end());
360 auto it = std::unique(ents_vec.begin(), ents_vec.end());
361 Range r;
362 r.insert_list(ents_vec.begin(), it);
363 return r;
364};
365
366/**
367 * @brief get entities of dofs in the problem - used for debugging
368 *
369 */
370auto get_dofs_ents_all = [](auto dm) {
371 auto prb_ptr = getProblemPtr(dm);
372 std::vector<EntityHandle> ents_vec;
373
374 MoFEM::Interface *m_field_ptr;
375 CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
376
377 auto dofs = prb_ptr->numeredRowDofsPtr;
378 ents_vec.reserve(dofs->size());
379
380 for (auto d : *dofs) {
381 ents_vec.push_back(d->getEnt());
382 }
383
384 std::sort(ents_vec.begin(), ents_vec.end());
385 auto it = std::unique(ents_vec.begin(), ents_vec.end());
386 Range r;
387 r.insert_list(ents_vec.begin(), it);
388 return r;
389};
390
391#include <FreeSurfaceOps.hpp>
392using namespace FreeSurfaceOps;
393
394struct FreeSurface;
395
396enum FR { F, R }; // F - forward, and reverse
397
398/**
399 * @brief Set of functions called by PETSc solver used to refine and update
400 * mesh.
401 *
402 * @note Currently theta method is only handled by this code.
403 *
404 */
405struct TSPrePostProc {
406
407 TSPrePostProc() = default;
408 virtual ~TSPrePostProc() = default;
409
410 /**
411 * @brief Used to setup TS solver
412 *
413 * @param ts
414 * @return MoFEMErrorCode
415 */
417
418 SmartPetscObj<VecScatter> getScatter(Vec x, Vec y, enum FR fr);
420
424
425private:
426 /**
427 * @brief Pre process time step
428 *
429 * Refine mesh and update fields
430 *
431 * @param ts
432 * @return MoFEMErrorCode
433 */
434 static MoFEMErrorCode tsPreProc(TS ts);
435
436 /**
437 * @brief Post process time step.
438 *
439 * Currently that function do not make anything major
440 *
441 * @param ts
442 * @return MoFEMErrorCode
443 */
444 static MoFEMErrorCode tsPostProc(TS ts);
445
447
448 static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t,
449 Vec f,
450 void *ctx); //< Wrapper for SNES Rhs
451 static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
452 PetscReal a, Mat A, Mat B,
453 void *ctx); ///< Wrapper for SNES Lhs
454 static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u,
455 void *ctx); ///< Wrapper for TS monitor
456 static MoFEMErrorCode pcSetup(PC pc);
457 static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x);
458
459 SmartPetscObj<Vec> globRes; //< global residual
460 SmartPetscObj<Mat> subB; //< sub problem tangent matrix
461 SmartPetscObj<KSP> subKSP; //< sub problem KSP solver
462
463 boost::shared_ptr<SnesCtx>
464 snesCtxPtr; //< infernal data (context) for MoFEM SNES fuctions
465 boost::shared_ptr<TsCtx>
466 tsCtxPtr; //< internal data (context) for MoFEM TS functions.
467};
468
469static boost::weak_ptr<TSPrePostProc> tsPrePostProc;
470
472
473 FreeSurface(MoFEM::Interface &m_field) : mField(m_field) {}
474
476
478
480
481private:
488
489 /// @brief Find entities on refinement levels
490 /// @param overlap level of overlap
491 /// @return
492 std::vector<Range> findEntitiesCrossedByPhaseInterface(size_t overlap);
493
494 /// @brief
495 /// @param ents
496 /// @param level
497 /// @param mask
498 /// @return
500
501 /// @brief
502 /// @param overlap
503 /// @return
504 MoFEMErrorCode refineMesh(size_t overlap);
505
506 /// @brief Create hierarchy of elements run on parents levels
507 /// @param fe_top pipeline element to which hierarchy is attached
508 /// @param op type of operator OPSPACE, OPROW, OPCOL or OPROWCOL
509 /// @param get_elema lambda function to create element on hierarchy
510 /// @param verbosity verbosity level
511 /// @param sev severity level
512 /// @return
514 boost::shared_ptr<FEMethod> fe_top, std::string field_name,
515 ForcesAndSourcesCore::UserDataOperator::OpType op,
516 BitRefLevel child_ent_bit,
517 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
518 int verbosity, LogManager::SeverityLevel sev);
519
520 friend struct TSPrePostProc;
521};
522
523//! [Run programme]
533//! [Run programme]
534
535//! [Read mesh]
538 MOFEM_LOG("FS", Sev::inform) << "Read mesh for problem";
540
542 simple->getBitRefLevel() = BitRefLevel();
543
544 CHKERR simple->getOptions();
545 CHKERR simple->loadFile();
546
548}
549//! [Read mesh]
550
551//! [Set up problem]
554
555 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
556 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-nb_levels", &nb_levels,
557 PETSC_NULLPTR);
558 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-refine_overlap", &refine_overlap,
559 PETSC_NULLPTR);
560
561 const char *coord_type_names[] = {"cartesian", "polar", "cylindrical",
562 "spherical"};
563 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-coords", coord_type_names,
564 LAST_COORDINATE_SYSTEM, &coord_type, PETSC_NULLPTR);
565
566 MOFEM_LOG("FS", Sev::inform) << "Approximation order = " << order;
567 MOFEM_LOG("FS", Sev::inform)
568 << "Number of refinement levels nb_levels = " << nb_levels;
569 nb_levels += 1;
570
572
573 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "Acceleration", "-a0", &a0, PETSC_NULLPTR);
574 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "\"Minus\" phase density rho_m",
575 "-rho_m", &rho_m, PETSC_NULLPTR);
576 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "\"Minus\" phase viscosity", "-mu_m",
577 &mu_m, PETSC_NULLPTR);
578 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "\"Plus\" phase density", "-rho_p",
579 &rho_p, PETSC_NULLPTR);
580 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "\"Plus\" phase viscosity", "-mu_p",
581 &mu_p, PETSC_NULLPTR);
582 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "Surface tension", "-lambda",
583 &lambda, PETSC_NULLPTR);
584 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR,
585 "Height of the well in energy functional", "-W",
586 &W, PETSC_NULLPTR);
587 rho_ave = (rho_p + rho_m) / 2;
588 rho_diff = (rho_p - rho_m) / 2;
589 mu_ave = (mu_p + mu_m) / 2;
590 mu_diff = (mu_p - mu_m) / 2;
591
592 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-h", &h, PETSC_NULLPTR);
593 eta = h;
594 eta2 = eta * eta;
595 kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
596
597 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-md", &eta, PETSC_NULLPTR);
598
599 MOFEM_LOG("FS", Sev::inform) << "Acceleration a0 = " << a0;
600 MOFEM_LOG("FS", Sev::inform) << "\"Minus\" phase density rho_m = " << rho_m;
601 MOFEM_LOG("FS", Sev::inform) << "\"Minus\" phase viscosity mu_m = " << mu_m;
602 MOFEM_LOG("FS", Sev::inform) << "\"Plus\" phase density rho_p = " << rho_p;
603 MOFEM_LOG("FS", Sev::inform) << "\"Plus\" phase viscosity mu_p = " << mu_p;
604 MOFEM_LOG("FS", Sev::inform) << "Surface tension lambda = " << lambda;
605 MOFEM_LOG("FS", Sev::inform)
606 << "Height of the well in energy functional W = " << W;
607 MOFEM_LOG("FS", Sev::inform) << "Characteristic mesh size h = " << h;
608 MOFEM_LOG("FS", Sev::inform) << "Mobility md = " << md;
609
610 MOFEM_LOG("FS", Sev::inform) << "Average density rho_ave = " << rho_ave;
611 MOFEM_LOG("FS", Sev::inform) << "Difference density rho_diff = " << rho_diff;
612 MOFEM_LOG("FS", Sev::inform) << "Average viscosity mu_ave = " << mu_ave;
613 MOFEM_LOG("FS", Sev::inform) << "Difference viscosity mu_diff = " << mu_diff;
614 MOFEM_LOG("FS", Sev::inform) << "kappa = " << kappa;
615
616
617 // Fields on domain
618
619 // Velocity field
620 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, U_FIELD_DIM);
621 // Pressure field
622 CHKERR simple->addDomainField("P", H1, AINSWORTH_LEGENDRE_BASE, 1);
623 // Order/phase fild
624 CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
625 // Chemical potential
626 CHKERR simple->addDomainField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
627
628 // Field on boundary
629 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE,
631 CHKERR simple->addBoundaryField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
632 CHKERR simple->addBoundaryField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
633 // Lagrange multiplier which constrains slip conditions
634 CHKERR simple->addBoundaryField("L", H1, AINSWORTH_LEGENDRE_BASE, 1);
635
636 CHKERR simple->setFieldOrder("U", order);
637 CHKERR simple->setFieldOrder("P", order - 1);
638 CHKERR simple->setFieldOrder("H", order);
639 CHKERR simple->setFieldOrder("G", order);
640 CHKERR simple->setFieldOrder("L", order);
641
642 // Initialise bit ref levels
643 auto set_problem_bit = [&]() {
645 // Set bits to build adjacencies between parents and children. That is
646 // used by simple interface.
647 simple->getBitAdjEnt() = BitRefLevel().set();
648 simple->getBitAdjParent() = BitRefLevel().set();
649 simple->getBitRefLevel() = BitRefLevel().set();
650 simple->getBitRefLevelMask() = BitRefLevel().set();
652 };
653
654 CHKERR set_problem_bit();
655
656 CHKERR simple->setUp();
657
659}
660//! [Set up problem]
661
662//! [Boundary condition]
665
666#ifdef PYTHON_INIT_SURFACE
667 auto get_py_surface_init = []() {
668 auto py_surf_init = boost::make_shared<SurfacePython>();
669 CHKERR py_surf_init->surfaceInit("surface.py");
670 surfacePythonWeakPtr = py_surf_init;
671 return py_surf_init;
672 };
673 auto py_surf_init = get_py_surface_init();
674#endif
675
677 auto pip_mng = mField.getInterface<PipelineManager>();
678 auto bc_mng = mField.getInterface<BcManager>();
679 auto bit_mng = mField.getInterface<BitRefManager>();
680 auto dm = simple->getDM();
681
682 using UDO = ForcesAndSourcesCore::UserDataOperator;
683
684 auto reset_bits = [&]() {
686 BitRefLevel start_mask;
687 for (auto s = 0; s != get_start_bit(); ++s)
688 start_mask[s] = true;
689 // reset bit ref levels
690 CHKERR bit_mng->lambdaBitRefLevel(
691 [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
692 Range level0;
693 CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
694 CHKERR bit_mng->setNthBitRefLevel(level0, get_current_bit(), true);
695 CHKERR bit_mng->setNthBitRefLevel(level0, get_projection_bit(), true);
697 };
698
699 auto add_parent_field = [&](auto fe, auto op, auto field) {
700 return setParentDofs(
701 fe, field, op, bit(get_skin_parent_bit()),
702
703 [&]() {
704 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
706 return fe_parent;
707 },
708
709 QUIET, Sev::noisy);
710 };
711
712 auto h_ptr = boost::make_shared<VectorDouble>();
713 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
714 auto g_ptr = boost::make_shared<VectorDouble>();
715 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
716
717 auto set_generic = [&](auto fe) {
719 auto &pip = fe->getOpPtrVector();
720
722
724 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
725
726 [&]() {
727 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
730 fe_parent->getOpPtrVector(), {H1});
731 return fe_parent;
732 },
733
734 QUIET, Sev::noisy);
735
736 CHKERR add_parent_field(fe, UDO::OPROW, "H");
737 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
738 pip.push_back(
739 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
740
741 CHKERR add_parent_field(fe, UDO::OPROW, "G");
742 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
743 pip.push_back(
744 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
745
747 };
748
749 auto post_proc = [&](auto exe_test) {
751 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
752 post_proc_fe->exeTestHook = exe_test;
753
754 CHKERR set_generic(post_proc_fe);
755
757
758 post_proc_fe->getOpPtrVector().push_back(
759
760 new OpPPMap(post_proc_fe->getPostProcMesh(),
761 post_proc_fe->getMapGaussPts(),
762
763 {{"H", h_ptr}, {"G", g_ptr}},
764
765 {{"GRAD_H", grad_h_ptr}, {"GRAD_G", grad_g_ptr}},
766
767 {},
768
769 {}
770
771 )
772
773 );
774
775 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
776 CHKERR post_proc_fe->writeFile("out_init.h5m");
777
779 };
780
781 auto solve_init = [&](auto exe_test) {
783
784 pip_mng->getOpDomainRhsPipeline().clear();
785 pip_mng->getOpDomainLhsPipeline().clear();
786
787 auto set_domain_rhs = [&](auto fe) {
789 CHKERR set_generic(fe);
790 auto &pip = fe->getOpPtrVector();
791
792 CHKERR add_parent_field(fe, UDO::OPROW, "H");
793 pip.push_back(new OpRhsH<true>("H", nullptr, nullptr, h_ptr, grad_h_ptr,
794 grad_g_ptr));
795 CHKERR add_parent_field(fe, UDO::OPROW, "G");
796 pip.push_back(new OpRhsG<true>("G", h_ptr, grad_h_ptr, g_ptr));
798 };
799
800 auto set_domain_lhs = [&](auto fe) {
802
803 CHKERR set_generic(fe);
804 auto &pip = fe->getOpPtrVector();
805
806 CHKERR add_parent_field(fe, UDO::OPROW, "H");
807 CHKERR add_parent_field(fe, UDO::OPCOL, "H");
808 pip.push_back(new OpLhsH_dH<true>("H", nullptr, h_ptr, grad_g_ptr));
809
810 CHKERR add_parent_field(fe, UDO::OPCOL, "G");
811 pip.push_back(new OpLhsH_dG<true>("H", "G", h_ptr));
812
813 CHKERR add_parent_field(fe, UDO::OPROW, "G");
814 pip.push_back(new OpLhsG_dG("G"));
815
816 CHKERR add_parent_field(fe, UDO::OPCOL, "H");
817 pip.push_back(new OpLhsG_dH<true>("G", "H", h_ptr));
819 };
820
821 auto create_subdm = [&]() {
822 auto level_ents_ptr = boost::make_shared<Range>();
823 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
824 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
825
826 DM subdm;
827 CHKERR DMCreate(mField.get_comm(), &subdm);
828 CHKERR DMSetType(subdm, "DMMOFEM");
829 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_INIT");
830 CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
831 CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
832 CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_TRUE);
833
834 for (auto f : {"H", "G"}) {
835 CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
836 CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
837 }
838 CHKERR DMSetUp(subdm);
839
840 if constexpr (debug) {
841 if (mField.get_comm_size() == 1) {
842 auto dm_ents = get_dofs_ents_all(SmartPetscObj<DM>(subdm, true));
843 CHKERR save_range(mField.get_moab(), "sub_dm_init_ents_verts.h5m",
844 dm_ents.subset_by_type(MBVERTEX));
845 dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
846 CHKERR save_range(mField.get_moab(), "sub_dm_init_ents.h5m", dm_ents);
847 }
848 }
849
850 return SmartPetscObj<DM>(subdm);
851 };
852
853 auto subdm = create_subdm();
854
855 pip_mng->getDomainRhsFE().reset();
856 pip_mng->getDomainLhsFE().reset();
857 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
858 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
859 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
860 pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
861
862 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
863 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
864
865 auto D = createDMVector(subdm);
866 auto snes = pip_mng->createSNES(subdm);
867 auto snes_ctx_ptr = getDMSnesCtx(subdm);
868
869 auto set_section_monitor = [&](auto solver) {
871 CHKERR SNESMonitorSet(solver,
872 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
874 (void *)(snes_ctx_ptr.get()), nullptr);
875
877 };
878
879 auto print_section_field = [&]() {
881
882 auto section =
883 mField.getInterface<ISManager>()->sectionCreate("SUB_INIT");
884 PetscInt num_fields;
885 CHKERR PetscSectionGetNumFields(section, &num_fields);
886 for (int f = 0; f < num_fields; ++f) {
887 const char *field_name;
888 CHKERR PetscSectionGetFieldName(section, f, &field_name);
889 MOFEM_LOG("FS", Sev::inform)
890 << "Field " << f << " " << std::string(field_name);
891 }
892
894 };
895
896 CHKERR set_section_monitor(snes);
897 CHKERR print_section_field();
898
899 for (auto f : {"U", "P", "H", "G", "L"}) {
900 CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
901 }
902
903 CHKERR SNESSetFromOptions(snes);
904 CHKERR SNESSolve(snes, PETSC_NULLPTR, D);
905
906 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
907 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
908 CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
909
911 };
912
913 CHKERR reset_bits();
914 CHKERR solve_init(
915 [](FEMethod *fe_ptr) { return get_fe_bit(fe_ptr).test(0); });
916 CHKERR refineMesh(refine_overlap);
917
918 for (auto f : {"U", "P", "H", "G", "L"}) {
919 CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
920 }
921 CHKERR solve_init([](FEMethod *fe_ptr) {
922 return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
923 });
924
925 CHKERR post_proc([](FEMethod *fe_ptr) {
926 return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
927 });
928
929 pip_mng->getOpDomainRhsPipeline().clear();
930 pip_mng->getOpDomainLhsPipeline().clear();
931
932 // Remove DOFs where boundary conditions are set
933 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
934 "U", 0, 0);
935 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
936 "L", 0, 0);
937 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
938 "U", 1, 1);
939 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
940 "L", 1, 1);
941 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "U",
942 0, SPACE_DIM);
943 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "L",
944 0, 0);
945 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "ZERO",
946 "L", 0, 0);
947
949}
950//! [Boundary condition]
951
952//! [Data projection]
955
956 // FIXME: Functionality in this method should be improved (projection should
957 // be done field by field), generalised, and move to become core
958 // functionality.
959
961 auto pip_mng = mField.getInterface<PipelineManager>();
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
977 using UDO = ForcesAndSourcesCore::UserDataOperator;
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_NULLPTR, glob_x, iy);
1425
1427 DMSubDomainRestrict(simple->getDM(), s, PETSC_NULLPTR, 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 for (auto v : ts_solver_vecs) {
1460 MOFEM_LOG("FS", Sev::inform) << "Solve projection vector";
1461
1462 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1463 SCATTER_REVERSE);
1464 CHKERR solve(sub_x);
1465
1466 for (auto f : {"U", "P", "H", "G", "L"}) {
1467 MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1468 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1469 }
1470
1471 CHKERR DMoFEMMeshToLocalVector(subdm, sub_x, INSERT_VALUES,
1472 SCATTER_REVERSE);
1473 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1474 SCATTER_FORWARD);
1475
1476 MOFEM_LOG("FS", Sev::inform) << "Norm V " << get_norm(v);
1477 }
1478
1479 CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_X0",
1480 &ts_solver_vecs[0]);
1481 CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_Xdot",
1482 &ts_solver_vecs[1]);
1483 }
1484
1485 for (auto f : {"U", "P", "H", "G", "L"}) {
1486 MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1487 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1488 }
1489 CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
1490 }
1491
1493 };
1494
1495 // postprocessing (only for debugging proposes)
1496 auto post_proc = [&](auto exe_test) {
1498 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
1499 auto &pip = post_proc_fe->getOpPtrVector();
1500 post_proc_fe->exeTestHook = exe_test;
1501
1503
1505 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1506
1507 [&]() {
1508 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1509 new DomianParentEle(mField));
1511 fe_parent->getOpPtrVector(), {H1});
1512 return fe_parent;
1513 },
1514
1515 QUIET, Sev::noisy);
1516
1518
1519 auto u_ptr = boost::make_shared<MatrixDouble>();
1520 auto p_ptr = boost::make_shared<VectorDouble>();
1521 auto h_ptr = boost::make_shared<VectorDouble>();
1522 auto g_ptr = boost::make_shared<VectorDouble>();
1523
1524 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U",
1526 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1527 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P",
1529 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1530 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H",
1532 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1533 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G",
1535 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1536
1537 post_proc_fe->getOpPtrVector().push_back(
1538
1539 new OpPPMap(post_proc_fe->getPostProcMesh(),
1540 post_proc_fe->getMapGaussPts(),
1541
1542 {{"P", p_ptr}, {"H", h_ptr}, {"G", g_ptr}},
1543
1544 {{"U", u_ptr}},
1545
1546 {},
1547
1548 {}
1549
1550 )
1551
1552 );
1553
1554 auto dm = simple->getDM();
1555 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
1556 CHKERR post_proc_fe->writeFile("out_projection.h5m");
1557
1559 };
1560
1561 CHKERR solve_projection([](FEMethod *fe_ptr) {
1562 return get_fe_bit(fe_ptr).test(get_current_bit());
1563 });
1564
1565 if constexpr (debug) {
1566 CHKERR post_proc([](FEMethod *fe_ptr) {
1567 return get_fe_bit(fe_ptr).test(get_current_bit());
1568 });
1569 }
1570
1571 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
1572 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
1573 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
1574 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
1575
1577}
1578//! [Data projection]
1579
1580//! [Push operators to pip]
1583 auto simple = mField.getInterface<Simple>();
1584
1585 using UDO = ForcesAndSourcesCore::UserDataOperator;
1586
1587 auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
1588 return setParentDofs(
1589 fe, field, op, bit(get_skin_parent_bit()),
1590
1591 [&]() {
1592 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1593 new DomianParentEle(mField));
1594 return fe_parent;
1595 },
1596
1597 QUIET, Sev::noisy);
1598 };
1599
1600 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
1601 return setParentDofs(
1602 fe, field, op, bit(get_skin_parent_bit()),
1603
1604 [&]() {
1605 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1607 return fe_parent;
1608 },
1609
1610 QUIET, Sev::noisy);
1611 };
1612
1613 auto test_bit_child = [](FEMethod *fe_ptr) {
1614 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1615 get_start_bit() + nb_levels - 1);
1616 };
1617
1618 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
1619 auto u_ptr = boost::make_shared<MatrixDouble>();
1620 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
1621 auto dot_h_ptr = boost::make_shared<VectorDouble>();
1622 auto h_ptr = boost::make_shared<VectorDouble>();
1623 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
1624 auto g_ptr = boost::make_shared<VectorDouble>();
1625 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
1626 auto lambda_ptr = boost::make_shared<VectorDouble>();
1627 auto p_ptr = boost::make_shared<VectorDouble>();
1628 auto div_u_ptr = boost::make_shared<VectorDouble>();
1629
1630 // Push element from reference configuration to current configuration in 3d
1631 // space
1632 auto set_domain_general = [&](auto fe) {
1634 auto &pip = fe->getOpPtrVector();
1635
1637
1639 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1640
1641 [&]() {
1642 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1643 new DomianParentEle(mField));
1645 fe_parent->getOpPtrVector(), {H1});
1646 return fe_parent;
1647 },
1648
1649 QUIET, Sev::noisy);
1650
1651 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1652 pip.push_back(
1654 pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
1656 "U", grad_u_ptr));
1657
1658 switch (coord_type) {
1659 case CARTESIAN:
1660 pip.push_back(
1662 "U", div_u_ptr));
1663 break;
1664 case CYLINDRICAL:
1665 pip.push_back(
1667 "U", div_u_ptr));
1668 break;
1669 default:
1670 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
1671 "Coordinate system not implemented");
1672 }
1673
1674 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1675 pip.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
1676 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1677 pip.push_back(
1678 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1679
1680 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1681 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1682 pip.push_back(
1683 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
1684
1685 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1686 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1688 };
1689
1690 auto set_domain_rhs = [&](auto fe) {
1692 auto &pip = fe->getOpPtrVector();
1693
1694 CHKERR set_domain_general(fe);
1695
1696 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1697 pip.push_back(new OpRhsU("U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
1698 grad_h_ptr, g_ptr, p_ptr));
1699
1700 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1701 pip.push_back(new OpRhsH<false>("H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
1702 grad_g_ptr));
1703
1704 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1705 pip.push_back(new OpRhsG<false>("G", h_ptr, grad_h_ptr, g_ptr));
1706
1707 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1708 pip.push_back(new OpDomainAssembleScalar(
1709 "P", div_u_ptr, [](const double r, const double, const double) {
1710 return cylindrical(r);
1711 }));
1712 pip.push_back(new OpDomainAssembleScalar(
1713 "P", p_ptr, [](const double r, const double, const double) {
1714 return eps * cylindrical(r);
1715 }));
1717 };
1718
1719 auto set_domain_lhs = [&](auto fe) {
1721 auto &pip = fe->getOpPtrVector();
1722
1723 CHKERR set_domain_general(fe);
1724
1725 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1726 {
1727 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1728 pip.push_back(new OpLhsU_dU("U", u_ptr, grad_u_ptr, h_ptr));
1729 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1730 pip.push_back(
1731 new OpLhsU_dH("U", "H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
1732
1733 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1734 pip.push_back(new OpLhsU_dG("U", "G", grad_h_ptr));
1735 }
1736
1737 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1738 {
1739 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1740 pip.push_back(new OpLhsH_dU("H", "U", grad_h_ptr));
1741 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1742 pip.push_back(new OpLhsH_dH<false>("H", u_ptr, h_ptr, grad_g_ptr));
1743 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1744 pip.push_back(new OpLhsH_dG<false>("H", "G", h_ptr));
1745 }
1746
1747 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1748 {
1749 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1750 pip.push_back(new OpLhsG_dH<false>("G", "H", h_ptr));
1751 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1752 pip.push_back(new OpLhsG_dG("G"));
1753 }
1754
1755 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1756 {
1757 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1758
1759 switch (coord_type) {
1760 case CARTESIAN:
1761 pip.push_back(new OpMixScalarTimesDiv<CARTESIAN>(
1762 "P", "U",
1763 [](const double r, const double, const double) {
1764 return cylindrical(r);
1765 },
1766 true, false));
1767 break;
1768 case CYLINDRICAL:
1769 pip.push_back(new OpMixScalarTimesDiv<CYLINDRICAL>(
1770 "P", "U",
1771 [](const double r, const double, const double) {
1772 return cylindrical(r);
1773 },
1774 true, false));
1775 break;
1776 default:
1777 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
1778 "Coordinate system not implemented");
1779 }
1780
1781 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P");
1782 pip.push_back(new OpDomainMassP("P", "P", [](double r, double, double) {
1783 return eps * cylindrical(r);
1784 }));
1785 }
1786
1788 };
1789
1790 auto get_block_name = [](auto name) {
1791 return boost::format("%s(.*)") % "WETTING_ANGLE";
1792 };
1793
1794 auto get_blocks = [&](auto &&name) {
1795 return mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
1796 std::regex(name.str()));
1797 };
1798
1799 auto set_boundary_rhs = [&](auto fe) {
1801 auto &pip = fe->getOpPtrVector();
1802
1804 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1805
1806 [&]() {
1807 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1809 return fe_parent;
1810 },
1811
1812 QUIET, Sev::noisy);
1813
1814 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
1815 pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
1816
1817 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
1818 pip.push_back(new OpCalculateScalarFieldValues("L", lambda_ptr));
1819 pip.push_back(new OpNormalConstrainRhs("L", u_ptr));
1820
1822 pip, mField, "L", {}, "INFLUX",
1823 [](double r, double, double) { return cylindrical(r); }, Sev::inform);
1824
1825 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
1826 pip.push_back(new OpNormalForceRhs("U", lambda_ptr));
1827
1828 auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
1829 if (wetting_block.size()) {
1830 // push operators to the side element which is called from op_bdy_side
1831 auto op_bdy_side =
1832 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
1833 op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
1834
1836 op_bdy_side->getOpPtrVector(), {H1});
1837
1839 op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
1841
1842 [&]() {
1843 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1844 new DomianParentEle(mField));
1846 fe_parent->getOpPtrVector(), {H1});
1847 return fe_parent;
1848 },
1849
1850 QUIET, Sev::noisy);
1851
1852 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1853 "H");
1854 op_bdy_side->getOpPtrVector().push_back(
1855 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1856 // push bdy side op
1857 pip.push_back(op_bdy_side);
1858
1859 // push operators for rhs wetting angle
1860 for (auto &b : wetting_block) {
1861 Range force_edges;
1862 std::vector<double> attr_vec;
1863 CHKERR b->getMeshsetIdEntitiesByDimension(
1864 mField.get_moab(), SPACE_DIM - 1, force_edges, true);
1865 b->getAttributes(attr_vec);
1866 if (attr_vec.size() != 1)
1867 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1868 "Should be one attribute");
1869 MOFEM_LOG("FS", Sev::inform) << "Wetting angle: " << attr_vec.front();
1870 // need to find the attributes and pass to operator
1871 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
1872 pip.push_back(new OpWettingAngleRhs(
1873 "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
1874 attr_vec.front()));
1875 }
1876 }
1877
1879 };
1880
1881 auto set_boundary_lhs = [&](auto fe) {
1883 auto &pip = fe->getOpPtrVector();
1884
1886 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1887
1888 [&]() {
1889 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1891 return fe_parent;
1892 },
1893
1894 QUIET, Sev::noisy);
1895
1896 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
1897 CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "U");
1898 pip.push_back(new OpNormalConstrainLhs("L", "U"));
1899
1900 auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
1901 if (wetting_block.size()) {
1902 auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
1903 auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
1904
1905 // push operators to the side element which is called from op_bdy_side
1906 auto op_bdy_side =
1907 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
1908 op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
1909
1911 op_bdy_side->getOpPtrVector(), {H1});
1912
1914 op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
1916
1917 [&]() {
1918 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1919 new DomianParentEle(mField));
1921 fe_parent->getOpPtrVector(), {H1});
1922 return fe_parent;
1923 },
1924
1925 QUIET, Sev::noisy);
1926
1927 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1928 "H");
1929 op_bdy_side->getOpPtrVector().push_back(
1930 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1931 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
1932 "H");
1933 op_bdy_side->getOpPtrVector().push_back(
1934 new OpLoopSideGetDataForSideEle("H", col_ind_ptr, col_diff_base_ptr));
1935
1936 // push bdy side op
1937 pip.push_back(op_bdy_side);
1938
1939 // push operators for lhs wetting angle
1940 for (auto &b : wetting_block) {
1941 Range force_edges;
1942 std::vector<double> attr_vec;
1943 CHKERR b->getMeshsetIdEntitiesByDimension(
1944 mField.get_moab(), SPACE_DIM - 1, force_edges, true);
1945 b->getAttributes(attr_vec);
1946 if (attr_vec.size() != 1)
1947 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1948 "Should be one attribute");
1949 MOFEM_LOG("FS", Sev::inform)
1950 << "wetting angle edges size " << force_edges.size();
1951
1952 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
1953 pip.push_back(new OpWettingAngleLhs(
1954 "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
1955 boost::make_shared<Range>(force_edges), attr_vec.front()));
1956 }
1957 }
1958
1960 };
1961
1962 auto *pip_mng = mField.getInterface<PipelineManager>();
1963
1964 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1965 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1966 CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
1967 CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
1968
1969 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
1970 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
1971 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
1972 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
1973
1974 pip_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
1975 pip_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
1976 pip_mng->getBoundaryLhsFE()->exeTestHook = test_bit_child;
1977 pip_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
1978
1980}
1981//! [Push operators to pip]
1982
1983/**
1984 * @brief Monitor solution
1985 *
1986 * This functions is called by TS kso at the end of each step. It is used
1987 */
1988struct Monitor : public FEMethod {
1990 SmartPetscObj<DM> dm, boost::shared_ptr<moab::Core> post_proc_mesh,
1991 boost::shared_ptr<PostProcEleDomainCont> post_proc,
1992 boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
1993 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
1994 p)
1995 : dM(dm), postProcMesh(post_proc_mesh), postProc(post_proc),
1996 postProcEdge(post_proc_edge), liftFE(p.first), liftVec(p.second) {}
1999
2000 MOFEM_LOG("FS", Sev::verbose) << "Monitor";
2001 constexpr int save_every_nth_step = 1;
2002 if (ts_step % save_every_nth_step == 0) {
2003 MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc";
2004 MoFEM::Interface *m_field_ptr;
2005 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
2006 auto post_proc_begin =
2007 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
2008 postProcMesh);
2009 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
2010 *m_field_ptr, postProcMesh);
2011 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
2013 this->getCacheWeakPtr());
2015 this->getCacheWeakPtr());
2016 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end->getFEMethod());
2017 CHKERR post_proc_end->writeFile(
2018 "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
2019 MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc done";
2020 }
2021
2022 liftVec->resize(SPACE_DIM, false);
2023 liftVec->clear();
2025 MPI_Allreduce(MPI_IN_PLACE, &(*liftVec)[0], SPACE_DIM, MPI_DOUBLE, MPI_SUM,
2026 MPI_COMM_WORLD);
2027 MOFEM_LOG("FS", Sev::inform)
2028 << "Step " << ts_step << " time " << ts_t
2029 << " lift vec x: " << (*liftVec)[0] << " y: " << (*liftVec)[1];
2030
2032 }
2033
2034private:
2036 boost::shared_ptr<moab::Core> postProcMesh;
2037 boost::shared_ptr<PostProcEleDomainCont> postProc;
2038 boost::shared_ptr<PostProcEleBdyCont> postProcEdge;
2039 boost::shared_ptr<BoundaryEle> liftFE;
2040 boost::shared_ptr<VectorDouble> liftVec;
2041};
2042
2043//! [Solve]
2046
2047 using UDO = ForcesAndSourcesCore::UserDataOperator;
2048
2049 auto simple = mField.getInterface<Simple>();
2050 auto pip_mng = mField.getInterface<PipelineManager>();
2051
2052 auto create_solver_dm = [&](auto dm) -> SmartPetscObj<DM> {
2053 DM subdm;
2054
2055 auto setup_subdm = [&](auto dm) {
2057 auto simple = mField.getInterface<Simple>();
2058 auto bit_mng = mField.getInterface<BitRefManager>();
2059 auto dm = simple->getDM();
2060 auto level_ents_ptr = boost::make_shared<Range>();
2061 CHKERR bit_mng->getEntitiesByRefLevel(
2062 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
2063 CHKERR DMCreate(mField.get_comm(), &subdm);
2064 CHKERR DMSetType(subdm, "DMMOFEM");
2065 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
2066 CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
2067 CHKERR DMMoFEMAddElement(subdm, simple->getBoundaryFEName());
2068 CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
2069 CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_FALSE);
2070 for (auto f : {"U", "P", "H", "G", "L"}) {
2071 CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
2072 CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
2073 }
2074 CHKERR DMSetUp(subdm);
2076 };
2077
2078 CHK_THROW_MESSAGE(setup_subdm(dm), "create subdm");
2079
2080 return SmartPetscObj<DM>(subdm);
2081 };
2082
2083 auto get_fe_post_proc = [&](auto post_proc_mesh) {
2084 auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
2085 return setParentDofs(
2086 fe, field, op, bit(get_skin_parent_bit()),
2087
2088 [&]() {
2089 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2090 new DomianParentEle(mField));
2091 return fe_parent;
2092 },
2093
2094 QUIET, Sev::noisy);
2095 };
2096
2097 auto post_proc_fe =
2098 boost::make_shared<PostProcEleDomainCont>(mField, post_proc_mesh);
2099 post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
2100 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2101 get_start_bit() + nb_levels - 1);
2102 };
2103
2104 auto u_ptr = boost::make_shared<MatrixDouble>();
2105 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2106 auto h_ptr = boost::make_shared<VectorDouble>();
2107 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2108 auto p_ptr = boost::make_shared<VectorDouble>();
2109 auto g_ptr = boost::make_shared<VectorDouble>();
2110 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2111
2113 post_proc_fe->getOpPtrVector(), {H1});
2114
2116 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2117
2118 [&]() {
2119 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2120 new DomianParentEle(mField));
2122 fe_parent->getOpPtrVector(), {H1});
2123 return fe_parent;
2124 },
2125
2126 QUIET, Sev::noisy);
2127
2128 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U");
2129 post_proc_fe->getOpPtrVector().push_back(
2131 post_proc_fe->getOpPtrVector().push_back(
2133 grad_u_ptr));
2134
2135 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H");
2136 post_proc_fe->getOpPtrVector().push_back(
2137 new OpCalculateScalarFieldValues("H", h_ptr));
2138 post_proc_fe->getOpPtrVector().push_back(
2139 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
2140
2141 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P");
2142 post_proc_fe->getOpPtrVector().push_back(
2143 new OpCalculateScalarFieldValues("P", p_ptr));
2144
2145 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G");
2146 post_proc_fe->getOpPtrVector().push_back(
2147 new OpCalculateScalarFieldValues("G", g_ptr));
2148 post_proc_fe->getOpPtrVector().push_back(
2149 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
2150
2152
2153 post_proc_fe->getOpPtrVector().push_back(
2154
2155 new OpPPMap(
2156 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2157
2158 {{"H", h_ptr}, {"P", p_ptr}, {"G", g_ptr}},
2159
2160 {{"U", u_ptr}, {"H_GRAD", grad_h_ptr}, {"G_GRAD", grad_g_ptr}},
2161
2162 {{"GRAD_U", grad_u_ptr}},
2163
2164 {}
2165
2166 )
2167
2168 );
2169
2170 return post_proc_fe;
2171 };
2172
2173 auto get_bdy_post_proc_fe = [&](auto post_proc_mesh) {
2174 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2175 return setParentDofs(
2176 fe, field, op, bit(get_skin_parent_bit()),
2177
2178 [&]() {
2179 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2180 new BoundaryParentEle(mField));
2181 return fe_parent;
2182 },
2183
2184 QUIET, Sev::noisy);
2185 };
2186
2187 auto post_proc_fe =
2188 boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
2189 post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
2190 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2191 get_start_bit() + nb_levels - 1);
2192 };
2193
2194 CHKERR setParentDofs(
2195 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2196
2197 [&]() {
2198 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2199 new BoundaryParentEle(mField));
2200 return fe_parent;
2201 },
2202
2203 QUIET, Sev::noisy);
2204
2205 struct OpGetNormal : public BoundaryEleOp {
2206
2207 OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
2208 boost::shared_ptr<MatrixDouble> n_ptr)
2209 : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), ptrL(l_ptr),
2210 ptrNormal(n_ptr) {}
2211
2212 MoFEMErrorCode doWork(int side, EntityType type,
2215 auto t_l = getFTensor0FromVec(*ptrL);
2216 auto t_n_fe = getFTensor1NormalsAtGaussPts();
2217 ptrNormal->resize(SPACE_DIM, getGaussPts().size2(), false);
2218 auto t_n = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
2219 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
2220 t_n(i) = t_n_fe(i) * t_l / std::sqrt(t_n_fe(i) * t_n_fe(i));
2221 ++t_n_fe;
2222 ++t_l;
2223 ++t_n;
2224 }
2226 };
2227
2228 protected:
2229 boost::shared_ptr<VectorDouble> ptrL;
2230 boost::shared_ptr<MatrixDouble> ptrNormal;
2231 };
2232
2233 auto u_ptr = boost::make_shared<MatrixDouble>();
2234 auto p_ptr = boost::make_shared<VectorDouble>();
2235 auto lambda_ptr = boost::make_shared<VectorDouble>();
2236 auto normal_l_ptr = boost::make_shared<MatrixDouble>();
2237
2238 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "U");
2239 post_proc_fe->getOpPtrVector().push_back(
2241
2242 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "L");
2243 post_proc_fe->getOpPtrVector().push_back(
2244 new OpCalculateScalarFieldValues("L", lambda_ptr));
2245 post_proc_fe->getOpPtrVector().push_back(
2246 new OpGetNormal(lambda_ptr, normal_l_ptr));
2247
2248 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "P");
2249 post_proc_fe->getOpPtrVector().push_back(
2250 new OpCalculateScalarFieldValues("P", p_ptr));
2251
2252 auto op_ptr = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
2253 op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
2254 EntityType type,
2258 };
2259
2261
2262 post_proc_fe->getOpPtrVector().push_back(
2263
2264 new OpPPMap(post_proc_fe->getPostProcMesh(),
2265 post_proc_fe->getMapGaussPts(),
2266
2267 OpPPMap::DataMapVec{{"P", p_ptr}},
2268
2269 OpPPMap::DataMapMat{{"U", u_ptr}, {"L", normal_l_ptr}},
2270
2272
2274
2275 )
2276
2277 );
2278
2279 return post_proc_fe;
2280 };
2281
2282 auto get_lift_fe = [&]() {
2283 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2284 return setParentDofs(
2285 fe, field, op, bit(get_skin_parent_bit()),
2286
2287 [&]() {
2288 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2289 new BoundaryParentEle(mField));
2290 return fe_parent;
2291 },
2292
2293 QUIET, Sev::noisy);
2294 };
2295
2296 auto fe = boost::make_shared<BoundaryEle>(mField);
2297 fe->exeTestHook = [](FEMethod *fe_ptr) {
2298 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2299 get_start_bit() + nb_levels - 1);
2300 };
2301
2302 auto lift_ptr = boost::make_shared<VectorDouble>();
2303 auto p_ptr = boost::make_shared<VectorDouble>();
2304 auto ents_ptr = boost::make_shared<Range>();
2305
2306 CHKERR setParentDofs(
2307 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2308
2309 [&]() {
2310 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2311 new BoundaryParentEle(mField));
2312 return fe_parent;
2313 },
2314
2315 QUIET, Sev::noisy);
2316
2317 std::vector<const CubitMeshSets *> vec_ptr;
2318 CHKERR mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
2319 std::regex("LIFT"), vec_ptr);
2320 for (auto m_ptr : vec_ptr) {
2321 auto meshset = m_ptr->getMeshset();
2322 Range ents;
2323 CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
2324 ents, true);
2325 ents_ptr->merge(ents);
2326 }
2327
2328 MOFEM_LOG("FS", Sev::noisy) << "Lift ents " << (*ents_ptr);
2329
2330 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "P");
2331 fe->getOpPtrVector().push_back(
2332 new OpCalculateScalarFieldValues("L", p_ptr));
2333 fe->getOpPtrVector().push_back(
2334 new OpCalculateLift("L", p_ptr, lift_ptr, ents_ptr));
2335
2336 return std::make_pair(fe, lift_ptr);
2337 };
2338
2339 auto set_post_proc_monitor = [&](auto ts) {
2341 DM dm;
2342 CHKERR TSGetDM(ts, &dm);
2343 boost::shared_ptr<FEMethod> null_fe;
2344 auto post_proc_mesh = boost::make_shared<moab::Core>();
2345 auto monitor_ptr = boost::make_shared<Monitor>(
2346 SmartPetscObj<DM>(dm, true), post_proc_mesh,
2347 get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
2348 get_lift_fe());
2349 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
2350 null_fe, monitor_ptr);
2352 };
2353
2354 auto dm = simple->getDM();
2355 auto ts = createTS(mField.get_comm());
2356 CHKERR TSSetDM(ts, dm);
2357
2358 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
2359 tsPrePostProc = ts_pre_post_proc;
2360
2361 if (auto ptr = tsPrePostProc.lock()) {
2362
2363 ptr->fsRawPtr = this;
2364 ptr->solverSubDM = create_solver_dm(simple->getDM());
2365 ptr->globSol = createDMVector(dm);
2366 CHKERR DMoFEMMeshToLocalVector(dm, ptr->globSol, INSERT_VALUES,
2367 SCATTER_FORWARD);
2368 CHKERR VecAssemblyBegin(ptr->globSol);
2369 CHKERR VecAssemblyEnd(ptr->globSol);
2370
2371 auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
2372
2373 CHKERR set_post_proc_monitor(sub_ts);
2374
2375 // Add monitor to time solver
2376 CHKERR TSSetFromOptions(ts);
2377 CHKERR ptr->tsSetUp(ts);
2378 CHKERR TSSetUp(ts);
2379
2380 auto print_fields_in_section = [&]() {
2382 auto section = mField.getInterface<ISManager>()->sectionCreate(
2383 simple->getProblemName());
2384 PetscInt num_fields;
2385 CHKERR PetscSectionGetNumFields(section, &num_fields);
2386 for (int f = 0; f < num_fields; ++f) {
2387 const char *field_name;
2388 CHKERR PetscSectionGetFieldName(section, f, &field_name);
2389 MOFEM_LOG("FS", Sev::inform)
2390 << "Field " << f << " " << std::string(field_name);
2391 }
2393 };
2394
2395 CHKERR print_fields_in_section();
2396
2397 CHKERR TSSolve(ts, ptr->globSol);
2398 }
2399
2401}
2402
2403int main(int argc, char *argv[]) {
2404
2405#ifdef PYTHON_INIT_SURFACE
2406 Py_Initialize();
2407#endif
2408
2409 // Initialisation of MoFEM/PETSc and MOAB data structures
2410 const char param_file[] = "param_file.petsc";
2411 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
2412
2413 // Add logging channel for example
2414 auto core_log = logging::core::get();
2415 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FS"));
2416 LogManager::setLog("FS");
2417 MOFEM_LOG_TAG("FS", "free surface");
2418
2419 try {
2420
2421 //! [Register MoFEM discrete manager in PETSc]
2422 DMType dm_name = "DMMOFEM";
2423 CHKERR DMRegister_MoFEM(dm_name);
2424 //! [Register MoFEM discrete manager in PETSc
2425
2426 //! [Create MoAB]
2427 moab::Core mb_instance; ///< mesh database
2428 moab::Interface &moab = mb_instance; ///< mesh database interface
2429 //! [Create MoAB]
2430
2431 //! [Create MoFEM]
2432 MoFEM::Core core(moab); ///< finite element database
2433 MoFEM::Interface &m_field = core; ///< finite element database interface
2434 //! [Create MoFEM]
2435
2436 //! [FreeSurface]
2437 FreeSurface ex(m_field);
2438 CHKERR ex.runProblem();
2439 //! [FreeSurface]
2440 }
2442
2444
2445#ifdef PYTHON_INIT_SURFACE
2446 if (Py_FinalizeEx() < 0) {
2447 exit(120);
2448 }
2449#endif
2450}
2451
2452std::vector<Range>
2454
2455 auto &moab = mField.get_moab();
2456 auto bit_mng = mField.getInterface<BitRefManager>();
2457 auto comm_mng = mField.getInterface<CommInterface>();
2458
2459 Range vertices;
2460 CHK_THROW_MESSAGE(bit_mng->getEntitiesByTypeAndRefLevel(
2461 bit(0), BitRefLevel().set(), MBVERTEX, vertices),
2462 "can not get vertices on bit 0");
2463
2464 auto &dofs_mi = mField.get_dofs()->get<Unique_mi_tag>();
2465 auto field_bit_number = mField.get_field_bit_number("H");
2466
2467 Range plus_range, minus_range;
2468 std::vector<EntityHandle> plus, minus;
2469
2470 // get vertices on level 0 on plus and minus side
2471 for (auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
2472
2473 const auto f = p->first;
2474 const auto s = p->second;
2475
2476 // Lowest Dof UId for given field (field bit number) on entity f
2477 const auto lo_uid = DofEntity::getLoFieldEntityUId(field_bit_number, f);
2478 const auto hi_uid = DofEntity::getHiFieldEntityUId(field_bit_number, s);
2479 auto it = dofs_mi.lower_bound(lo_uid);
2480 const auto hi_it = dofs_mi.upper_bound(hi_uid);
2481
2482 plus.clear();
2483 minus.clear();
2484 plus.reserve(std::distance(it, hi_it));
2485 minus.reserve(std::distance(it, hi_it));
2486
2487 for (; it != hi_it; ++it) {
2488 const auto v = (*it)->getFieldData();
2489 if (v > 0)
2490 plus.push_back((*it)->getEnt());
2491 else
2492 minus.push_back((*it)->getEnt());
2493 }
2494
2495 plus_range.insert_list(plus.begin(), plus.end());
2496 minus_range.insert_list(minus.begin(), minus.end());
2497 }
2498
2499 MOFEM_LOG_CHANNEL("SYNC");
2500 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
2501 << "Plus range " << plus_range << endl;
2502 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
2503 << "Minus range " << minus_range << endl;
2505
2506 auto get_elems = [&](auto &ents, auto bit, auto mask) {
2507 Range adj;
2508 CHK_MOAB_THROW(moab.get_adjacencies(ents, SPACE_DIM, false, adj,
2509 moab::Interface::UNION),
2510 "can not get adjacencies");
2511 CHK_THROW_MESSAGE(bit_mng->filterEntitiesByRefLevel(bit, mask, adj),
2512 "can not filter elements with bit 0");
2513 return adj;
2514 };
2515
2516 CHKERR comm_mng->synchroniseEntities(plus_range);
2517 CHKERR comm_mng->synchroniseEntities(minus_range);
2518
2519 std::vector<Range> ele_plus(nb_levels), ele_minus(nb_levels);
2520 ele_plus[0] = get_elems(plus_range, bit(0), BitRefLevel().set());
2521 ele_minus[0] = get_elems(minus_range, bit(0), BitRefLevel().set());
2522 auto common = intersect(ele_plus[0], ele_minus[0]);
2523 ele_plus[0] = subtract(ele_plus[0], common);
2524 ele_minus[0] = subtract(ele_minus[0], common);
2525
2526 auto get_children = [&](auto &p, auto &c) {
2528 CHKERR bit_mng->updateRangeByChildren(p, c);
2529 c = c.subset_by_dimension(SPACE_DIM);
2531 };
2532
2533 for (auto l = 1; l != nb_levels; ++l) {
2534 CHK_THROW_MESSAGE(get_children(ele_plus[l - 1], ele_plus[l]),
2535 "get children");
2536 CHK_THROW_MESSAGE(get_children(ele_minus[l - 1], ele_minus[l]),
2537 "get children");
2538 }
2539
2540 auto get_level = [&](auto &p, auto &m, auto z, auto bit, auto mask) {
2541 Range l;
2543 bit_mng->getEntitiesByDimAndRefLevel(bit, mask, SPACE_DIM, l),
2544 "can not get vertices on bit");
2545 l = subtract(l, p);
2546 l = subtract(l, m);
2547 for (auto f = 0; f != z; ++f) {
2548 Range conn;
2549 CHK_MOAB_THROW(moab.get_connectivity(l, conn, true), "");
2550 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(conn);
2551 l = get_elems(conn, bit, mask);
2552 }
2553 return l;
2554 };
2555
2556 std::vector<Range> vec_levels(nb_levels);
2557 for (auto l = nb_levels - 1; l >= 0; --l) {
2558 vec_levels[l] = get_level(ele_plus[l], ele_minus[l], 2 * overlap, bit(l),
2559 BitRefLevel().set());
2560 }
2561
2562 if constexpr (debug) {
2563 if (mField.get_comm_size() == 1) {
2564 for (auto l = 0; l != nb_levels; ++l) {
2565 std::string name = (boost::format("out_r%d.h5m") % l).str();
2566 CHK_THROW_MESSAGE(save_range(mField.get_moab(), name, vec_levels[l]),
2567 "save mesh");
2568 }
2569 }
2570 }
2571
2572 return vec_levels;
2573}
2574
2577
2578 auto bit_mng = mField.getInterface<BitRefManager>();
2579
2580 BitRefLevel start_mask;
2581 for (auto s = 0; s != get_start_bit(); ++s)
2582 start_mask[s] = true;
2583
2584 // store prev_level
2585 Range prev_level;
2586 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_current_bit()),
2587 BitRefLevel().set(), prev_level);
2588 Range prev_level_skin;
2589 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_parent_bit()),
2590 BitRefLevel().set(), prev_level_skin);
2591 // reset bit ref levels
2592 CHKERR bit_mng->lambdaBitRefLevel(
2593 [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
2594 CHKERR bit_mng->setNthBitRefLevel(prev_level, get_projection_bit(), true);
2595 CHKERR bit_mng->setNthBitRefLevel(prev_level_skin, get_skin_projection_bit(),
2596 true);
2597
2598 // set refinement levels
2599 auto set_levels = [&](auto &&
2600 vec_levels /*entities are refined on each level*/) {
2602
2603 // start with zero level, which is the coarsest mesh
2604 Range level0;
2605 CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
2606 CHKERR bit_mng->setNthBitRefLevel(level0, get_start_bit(), true);
2607
2608 // get lower dimension entities
2609 auto get_adj = [&](auto ents) {
2610 Range conn;
2611 CHK_MOAB_THROW(mField.get_moab().get_connectivity(ents, conn, true),
2612 "get conn");
2613 for (auto d = 1; d != SPACE_DIM; ++d) {
2614 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
2615 ents.subset_by_dimension(SPACE_DIM), d, false, conn,
2616 moab::Interface::UNION),
2617 "get adj");
2618 }
2619 ents.merge(conn);
2620 return ents;
2621 };
2622
2623 // set bit levels
2624 for (auto l = 1; l != nb_levels; ++l) {
2625 Range level_prev;
2626 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_start_bit() + l - 1),
2627 BitRefLevel().set(),
2628 SPACE_DIM, level_prev);
2629 Range parents;
2630 CHKERR bit_mng->updateRangeByParent(vec_levels[l], parents);
2631 // subtract entities from previous level, which are refined, so should be
2632 // not there
2633 level_prev = subtract(level_prev, parents);
2634 // and instead add their children
2635 level_prev.merge(vec_levels[l]);
2636 // set bit to each level
2637 CHKERR bit_mng->setNthBitRefLevel(level_prev, get_start_bit() + l, true);
2638 }
2639
2640 // set bit levels to lower dimension entities
2641 for (auto l = 1; l != nb_levels; ++l) {
2642 Range level;
2643 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2644 bit(get_start_bit() + l), BitRefLevel().set(), SPACE_DIM, level);
2645 level = get_adj(level);
2646 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(level);
2647 CHKERR bit_mng->setNthBitRefLevel(level, get_start_bit() + l, true);
2648 }
2649
2651 };
2652
2653 // resolve skin between refined levels
2654 auto set_skins = [&]() {
2656
2657 moab::Skinner skinner(&mField.get_moab());
2658 ParallelComm *pcomm =
2659 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2660
2661 // get skin of bit level
2662 auto get_bit_skin = [&](BitRefLevel bit, BitRefLevel mask) {
2663 Range bit_ents;
2666 bit, mask, SPACE_DIM, bit_ents),
2667 "can't get bit level");
2668 Range bit_skin;
2669 CHK_MOAB_THROW(skinner.find_skin(0, bit_ents, false, bit_skin),
2670 "can't get skin");
2671 return bit_skin;
2672 };
2673
2674 auto get_level_skin = [&]() {
2675 Range skin;
2676 BitRefLevel bit_prev;
2677 for (auto l = 1; l != nb_levels; ++l) {
2678 auto skin_level_mesh = get_bit_skin(bit(l), BitRefLevel().set());
2679 // filter (remove) all entities which are on partitions boarders
2680 CHKERR pcomm->filter_pstatus(skin_level_mesh,
2681 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2682 PSTATUS_NOT, -1, nullptr);
2683 auto skin_level =
2684 get_bit_skin(bit(get_start_bit() + l), BitRefLevel().set());
2685 skin_level = subtract(skin_level,
2686 skin_level_mesh); // get only internal skins, not
2687 // on the body boundary
2688 // get lower dimension adjacencies (FIXME: add edges if 3D)
2689 Range skin_level_verts;
2690 CHKERR mField.get_moab().get_connectivity(skin_level, skin_level_verts,
2691 true);
2692 skin_level.merge(skin_level_verts);
2693
2694 // remove previous level
2695 bit_prev.set(l - 1);
2696 Range level_prev;
2697 CHKERR bit_mng->getEntitiesByRefLevel(bit_prev, BitRefLevel().set(),
2698 level_prev);
2699 skin.merge(subtract(skin_level, level_prev));
2700 }
2701
2702 return skin;
2703 };
2704
2705 auto resolve_shared = [&](auto &&skin) {
2706 Range tmp_skin = skin;
2707
2708 map<int, Range> map_procs;
2709 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2710 tmp_skin, &map_procs);
2711
2712 Range from_other_procs; // entities which also exist on other processors
2713 for (auto &m : map_procs) {
2714 if (m.first != mField.get_comm_rank()) {
2715 from_other_procs.merge(m.second);
2716 }
2717 }
2718
2719 auto common = intersect(
2720 skin, from_other_procs); // entities which are on internal skin
2721 skin.merge(from_other_procs);
2722
2723 // entities which are on processors boards, and several processors are not
2724 // true skin.
2725 if (!common.empty()) {
2726 // skin is internal exist on other procs
2727 skin = subtract(skin, common);
2728 }
2729
2730 return skin;
2731 };
2732
2733 // get parents of entities
2734 auto get_parent_level_skin = [&](auto skin) {
2735 Range skin_parents;
2736 CHKERR bit_mng->updateRangeByParent(
2737 skin.subset_by_dimension(SPACE_DIM - 1), skin_parents);
2738 Range skin_parent_verts;
2739 CHKERR mField.get_moab().get_connectivity(skin_parents, skin_parent_verts,
2740 true);
2741 skin_parents.merge(skin_parent_verts);
2742 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2743 skin_parents);
2744 return skin_parents;
2745 };
2746
2747 auto child_skin = resolve_shared(get_level_skin());
2748 auto parent_skin = get_parent_level_skin(child_skin);
2749
2750 child_skin = subtract(child_skin, parent_skin);
2751 CHKERR bit_mng->setNthBitRefLevel(child_skin, get_skin_child_bit(), true);
2752 CHKERR bit_mng->setNthBitRefLevel(parent_skin, get_skin_parent_bit(), true);
2753
2755 };
2756
2757 // take last level, remove childs on boarder, and set bit
2758 auto set_current = [&]() {
2760 Range last_level;
2761 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_start_bit() + nb_levels - 1),
2762 BitRefLevel().set(), last_level);
2763 Range skin_child;
2764 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_child_bit()),
2765 BitRefLevel().set(), skin_child);
2766
2767 last_level = subtract(last_level, skin_child);
2768 CHKERR bit_mng->setNthBitRefLevel(last_level, get_current_bit(), true);
2770 };
2771
2772 // set bits to levels
2773 CHKERR set_levels(findEntitiesCrossedByPhaseInterface(overlap));
2774 // set bits to skin
2775 CHKERR set_skins();
2776 // set current level bit
2777 CHKERR set_current();
2778
2779 if constexpr (debug) {
2780 if (mField.get_comm_size() == 1) {
2781 for (auto l = 0; l != nb_levels; ++l) {
2782 std::string name = (boost::format("out_level%d.h5m") % l).str();
2783 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_start_bit() + l),
2784 BitRefLevel().set(), name.c_str(), "MOAB",
2785 "PARALLEL=WRITE_PART");
2786 }
2787 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_current_bit()),
2788 BitRefLevel().set(), "current_bit.h5m",
2789 "MOAB", "PARALLEL=WRITE_PART");
2790 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_projection_bit()),
2791 BitRefLevel().set(), "projection_bit.h5m",
2792 "MOAB", "PARALLEL=WRITE_PART");
2793
2794 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_child_bit()),
2795 BitRefLevel().set(), "skin_child_bit.h5m",
2796 "MOAB", "PARALLEL=WRITE_PART");
2797 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_parent_bit()),
2798 BitRefLevel().set(), "skin_parent_bit.h5m",
2799 "MOAB", "PARALLEL=WRITE_PART");
2800 }
2801 }
2802
2804};
2805
2807 boost::shared_ptr<FEMethod> fe_top, std::string field_name,
2808 ForcesAndSourcesCore::UserDataOperator::OpType op,
2809 BitRefLevel child_ent_bit,
2810 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
2811 int verbosity, LogManager::SeverityLevel sev) {
2813
2814 /**
2815 * @brief Collect data from parent elements to child
2816 */
2817 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
2818 add_parent_level =
2819 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
2820 // Evaluate if not last parent element
2821 if (level > 0) {
2822
2823 // Create domain parent FE
2824 auto fe_ptr_current = get_elem();
2825
2826 // Call next level
2827 add_parent_level(
2828 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
2829 fe_ptr_current),
2830 level - 1);
2831
2832 // Add data to curent fe level
2833 if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
2834
2835 // Only base
2836 parent_fe_pt->getOpPtrVector().push_back(
2837
2839
2840 H1, op, fe_ptr_current,
2841
2842 BitRefLevel().set(), BitRefLevel().set(),
2843
2844 child_ent_bit, BitRefLevel().set(),
2845
2846 verbosity, sev));
2847
2848 } else {
2849
2850 // Filed data
2851 parent_fe_pt->getOpPtrVector().push_back(
2852
2854
2855 field_name, op, fe_ptr_current,
2856
2857 BitRefLevel().set(), BitRefLevel().set(),
2858
2859 child_ent_bit, BitRefLevel().set(),
2860
2861 verbosity, sev));
2862 }
2863 }
2864 };
2865
2866 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
2867 nb_levels);
2868
2870}
2871
2874
2875 if (auto ptr = tsPrePostProc.lock()) {
2876
2877 /**
2878 * @brief cut-off values at nodes, i.e. abs("H") <= 1
2879 *
2880 */
2881 auto cut_off_dofs = [&]() {
2883
2884 auto &m_field = ptr->fsRawPtr->mField;
2885
2886 Range current_verts;
2887 auto bit_mng = m_field.getInterface<BitRefManager>();
2889 bit(get_current_bit()), BitRefLevel().set(), MBVERTEX, current_verts);
2890
2891 auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
2893 for (auto &h : ent_ptr->getEntFieldData()) {
2894 h = cut_off(h);
2895 }
2897 };
2898
2899 auto field_blas = m_field.getInterface<FieldBlas>();
2900 CHKERR field_blas->fieldLambdaOnEntities(cut_off_verts, "H",
2901 &current_verts);
2903 };
2904
2905 CHKERR cut_off_dofs();
2906 }
2907
2908 if (auto ptr = tsPrePostProc.lock()) {
2909 MOFEM_LOG("FS", Sev::inform) << "Run step pre proc";
2910
2911 auto &m_field = ptr->fsRawPtr->mField;
2912 auto simple = m_field.getInterface<Simple>();
2913
2914 // get vector norm
2915 auto get_norm = [&](auto x) {
2916 double nrm;
2917 CHKERR VecNorm(x, NORM_2, &nrm);
2918 return nrm;
2919 };
2920
2921 // refine problem and project data, including theta data
2922 auto refine_problem = [&]() {
2924 MOFEM_LOG("FS", Sev::inform) << "Refine problem";
2925 CHKERR ptr->fsRawPtr->refineMesh(refine_overlap);
2926 CHKERR ptr->fsRawPtr->projectData();
2928 };
2929
2930 // set new jacobin operator, since problem and thus tangent matrix size has
2931 // changed
2932 auto set_jacobian_operators = [&]() {
2934 ptr->subB = createDMMatrix(ptr->solverSubDM);
2935 CHKERR KSPReset(ptr->subKSP);
2937 };
2938
2939 // set new solution
2940 auto set_solution = [&]() {
2942 MOFEM_LOG("FS", Sev::inform) << "Set solution";
2943
2944 PetscObjectState state;
2945
2946 // Record the state, and set it again. This is to fool PETSc that solution
2947 // vector is not updated. Otherwise PETSc will treat every step as a first
2948 // step.
2949
2950 // globSol is updated as result mesh refinement - this is not really set
2951 // a new solution.
2952
2953 CHKERR PetscObjectStateGet(getPetscObject(ptr->globSol.get()), &state);
2954 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), ptr->globSol,
2955 INSERT_VALUES, SCATTER_FORWARD);
2956 CHKERR PetscObjectStateSet(getPetscObject(ptr->globSol.get()), state);
2957 MOFEM_LOG("FS", Sev::verbose)
2958 << "Set solution, vector norm " << get_norm(ptr->globSol);
2960 };
2961
2962 PetscBool is_theta;
2963 PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
2964 if (is_theta) {
2965
2966 CHKERR refine_problem(); // refine problem
2967 CHKERR set_jacobian_operators(); // set new jacobian
2968 CHKERR set_solution(); // set solution
2969
2970 } else {
2971 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
2972 "Sorry, only TSTheta handling is implemented");
2973 }
2974
2975 // Need barriers, somme functions in TS solver need are called collectively
2976 // and requite the same state of variables
2977 PetscBarrier((PetscObject)ts);
2978
2979 MOFEM_LOG_CHANNEL("SYNC");
2980 MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PreProc done";
2981 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
2982 }
2983
2985}
2986
2988 if (auto ptr = tsPrePostProc.lock()) {
2989 auto &m_field = ptr->fsRawPtr->mField;
2990 MOFEM_LOG_CHANNEL("SYNC");
2991 MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PostProc done";
2992 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
2993 }
2994 return 0;
2995}
2996
2997MoFEMErrorCode TSPrePostProc::tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t,
2998 Vec f, void *ctx) {
3000 if (auto ptr = tsPrePostProc.lock()) {
3001 auto sub_u = ptr->getSubVector();
3002 auto sub_u_t = vectorDuplicate(sub_u);
3003 auto sub_f = vectorDuplicate(sub_u);
3004 auto scatter = ptr->getScatter(sub_u, u, R);
3005
3006 auto apply_scatter_and_update = [&](auto x, auto sub_x) {
3008 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3009 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3010 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3011 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3013 };
3014
3015 CHKERR apply_scatter_and_update(u, sub_u);
3016 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3017
3018 CHKERR TsSetIFunction(ts, t, sub_u, sub_u_t, sub_f, ptr->tsCtxPtr.get());
3019 CHKERR VecScatterBegin(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3020 CHKERR VecScatterEnd(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3021 }
3023}
3024
3025MoFEMErrorCode TSPrePostProc::tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
3026 PetscReal a, Mat A, Mat B,
3027 void *ctx) {
3029 if (auto ptr = tsPrePostProc.lock()) {
3030 auto sub_u = ptr->getSubVector();
3031 auto sub_u_t = vectorDuplicate(sub_u);
3032 auto scatter = ptr->getScatter(sub_u, u, R);
3033
3034 auto apply_scatter_and_update = [&](auto x, auto sub_x) {
3036 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3037 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3038 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3039 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3041 };
3042
3043 CHKERR apply_scatter_and_update(u, sub_u);
3044 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3045
3046 CHKERR TsSetIJacobian(ts, t, sub_u, sub_u_t, a, ptr->subB, ptr->subB,
3047 ptr->tsCtxPtr.get());
3048 }
3050}
3051
3052MoFEMErrorCode TSPrePostProc::tsMonitor(TS ts, PetscInt step, PetscReal t,
3053 Vec u, void *ctx) {
3055 if (auto ptr = tsPrePostProc.lock()) {
3056 auto get_norm = [&](auto x) {
3057 double nrm;
3058 CHKERR VecNorm(x, NORM_2, &nrm);
3059 return nrm;
3060 };
3061
3062 auto sub_u = ptr->getSubVector();
3063 auto scatter = ptr->getScatter(sub_u, u, R);
3064 CHKERR VecScatterBegin(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3065 CHKERR VecScatterEnd(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3066 CHKERR VecGhostUpdateBegin(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3067 CHKERR VecGhostUpdateEnd(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3068
3069 MOFEM_LOG("FS", Sev::verbose)
3070 << "u norm " << get_norm(u) << " u sub nom " << get_norm(sub_u);
3071
3072 CHKERR TsMonitorSet(ts, step, t, sub_u, ptr->tsCtxPtr.get());
3073 }
3075}
3076
3079 if (auto ptr = tsPrePostProc.lock()) {
3080 MOFEM_LOG("FS", Sev::verbose) << "SetUP sub PC";
3081 ptr->subKSP = createKSP(ptr->fsRawPtr->mField.get_comm());
3082 CHKERR KSPSetFromOptions(ptr->subKSP);
3083 }
3085};
3086
3087MoFEMErrorCode TSPrePostProc::pcApply(PC pc, Vec pc_f, Vec pc_x) {
3089 if (auto ptr = tsPrePostProc.lock()) {
3090 auto sub_x = ptr->getSubVector();
3091 auto sub_f = vectorDuplicate(sub_x);
3092 auto scatter = ptr->getScatter(sub_x, pc_x, R);
3093 CHKERR VecScatterBegin(scatter, pc_f, sub_f, INSERT_VALUES,
3094 SCATTER_REVERSE);
3095 CHKERR VecScatterEnd(scatter, pc_f, sub_f, INSERT_VALUES, SCATTER_REVERSE);
3096 CHKERR KSPSetOperators(ptr->subKSP, ptr->subB, ptr->subB);
3097 MOFEM_LOG("FS", Sev::verbose) << "PCShell solve";
3098 CHKERR KSPSolve(ptr->subKSP, sub_f, sub_x);
3099 MOFEM_LOG("FS", Sev::verbose) << "PCShell solve <- done";
3100 CHKERR VecScatterBegin(scatter, sub_x, pc_x, INSERT_VALUES,
3101 SCATTER_FORWARD);
3102 CHKERR VecScatterEnd(scatter, sub_x, pc_x, INSERT_VALUES, SCATTER_FORWARD);
3103 }
3105};
3106
3108 if (auto ptr = tsPrePostProc.lock()) {
3109 auto prb_ptr = ptr->fsRawPtr->mField.get_problem("SUB_SOLVER");
3110 if (auto sub_data = prb_ptr->getSubData()) {
3111 auto is = sub_data->getSmartColIs();
3112 VecScatter s;
3113 if (fr == R) {
3114 CHK_THROW_MESSAGE(VecScatterCreate(x, PETSC_NULLPTR, y, is, &s),
3115 "crate scatter");
3116 } else {
3117 CHK_THROW_MESSAGE(VecScatterCreate(x, is, y, PETSC_NULLPTR, &s),
3118 "crate scatter");
3119 }
3120 return SmartPetscObj<VecScatter>(s);
3121 }
3122 }
3125}
3126
3130
3132
3133 auto &m_field = fsRawPtr->mField;
3134 auto simple = m_field.getInterface<Simple>();
3135
3137
3138 auto dm = simple->getDM();
3139
3141 CHKERR TSSetIFunction(ts, globRes, tsSetIFunction, nullptr);
3142 CHKERR TSSetIJacobian(ts, PETSC_NULLPTR, PETSC_NULLPTR, tsSetIJacobian, nullptr);
3143 CHKERR TSMonitorSet(ts, tsMonitor, fsRawPtr, PETSC_NULLPTR);
3144
3145 SNES snes;
3146 CHKERR TSGetSNES(ts, &snes);
3147 auto snes_ctx_ptr = getDMSnesCtx(dm);
3148
3149 auto set_section_monitor = [&](auto snes) {
3151 CHKERR SNESMonitorSet(snes,
3152 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
3153 void *))MoFEMSNESMonitorFields,
3154 (void *)(snes_ctx_ptr.get()), nullptr);
3156 };
3157
3158 CHKERR set_section_monitor(snes);
3159
3160 auto ksp = createKSP(m_field.get_comm());
3161 CHKERR KSPSetType(ksp, KSPPREONLY); // Run KSP internally in ShellPC
3162 auto sub_pc = createPC(fsRawPtr->mField.get_comm());
3163 CHKERR PCSetType(sub_pc, PCSHELL);
3164 CHKERR PCShellSetSetUp(sub_pc, pcSetup);
3165 CHKERR PCShellSetApply(sub_pc, pcApply);
3166 CHKERR KSPSetPC(ksp, sub_pc);
3167 CHKERR SNESSetKSP(snes, ksp);
3168
3171
3172 CHKERR TSSetPreStep(ts, tsPreProc);
3173 CHKERR TSSetPostStep(ts, tsPostProc);
3174
3176}
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
int main()
constexpr int SPACE_DIM
constexpr double a
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
BoundaryEle::UserDataOperator BoundaryEleOp
Kronecker Delta class symmetric.
@ QUIET
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
#define BITREFLEVEL_SIZE
max number of refinements
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ H1
continuous field
Definition definitions.h:85
@ NOSPACE
Definition definitions.h:83
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ 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
@ CYLINDRICAL
@ CARTESIAN
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
auto cylindrical
double kappa
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
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
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
@ 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
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
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
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
int coord_type
double W
auto kernel_oscillation
double eta2
auto get_skin_parent_bit
auto get_M0
auto get_f_dh
auto get_M3
double rho_p
double mu_p
auto get_start_bit
OpDomainMassH OpDomainMassP
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
constexpr int BASE_DIM
int nb_levels
auto wetting_angle
double h
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
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
auto get_D
auto wetting_angle_sub_stepping
int refine_overlap
auto get_M0_dh
auto my_max
auto cut_off
auto get_skin_child_bit
int order
approximation order
SideEle::UserDataOperator SideOp
double a0
double rho_m
auto get_current_bit
dofs bit used to do calculations
double eta
auto bit
double md
auto capillary_tube
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
auto get_f
auto get_projection_bit
auto d_cut_off
auto my_min
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:215
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
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:114
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:546
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:514
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
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:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition DMMoFEM.cpp:410
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:280
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1059
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:536
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
MoFEMErrorCode getEntitiesByTypeAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const EntityType type, 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.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
MoFEMErrorCode buildSubProblem(const std::string out_name, const std::vector< std::string > &fields_row, const std::vector< std::string > &fields_col, const std::string main_problem, const bool square_matrix=true, const map< std::string, boost::shared_ptr< Range > > *entityMapRow=nullptr, const map< std::string, boost::shared_ptr< Range > > *entityMapCol=nullptr, int verb=VERBOSE)
build sub problem
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)
const double n
refractive index of diffusive medium
PipelineManager::ElementsAndOpsByDim< FE_DIM >::DomianParentEle DomianParentEle
Definition level_set.cpp:38
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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:1046
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition TsCtx.cpp:169
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition TsCtx.cpp:263
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1144
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition DMMoFEM.cpp:434
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
static const bool debug
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:592
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:56
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.
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
Definition Schur.cpp:1296
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
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.
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.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition DMMoFEM.hpp:1130
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:1047
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:58
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
constexpr double g
FTensor::Index< 'm', 3 > m
MoFEM::Interface & mField
Definition plastic.cpp:226
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:29
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:118
base operator to do operations at Gauss Pt. level
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
MoFEMErrorCode fieldLambdaOnEntities(OneFieldFunctionOnEntities lambda, const std::string field_name, Range *ents_ptr=nullptr)
field lambda
Definition FieldBlas.cpp:50
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
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.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
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 field rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor field 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< MatrixDouble > > DataMapMat
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Operator to execute finite element instance on parent element. This operator is typically used to pro...
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
bool & getParentAdjacencies()
Get the addParentAdjacencies.
Definition Simple.hpp:514
intrusive_ptr for managing petsc objects
PetscReal ts_t
time
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference 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
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