v0.14.0
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 
15 using namespace MoFEM;
16 
17 static 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>
24 namespace bp = boost::python;
25 
26 struct 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 
68 private:
69  bp::object mainNamespace;
70  bp::object surfaceFun;
71 };
72 
73 static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
74 
75 #endif
76 
77 int coord_type = EXECUTABLE_COORD_TYPE;
78 
79 constexpr int BASE_DIM = 1;
80 constexpr int SPACE_DIM = 2;
81 constexpr int U_FIELD_DIM = SPACE_DIM;
82 
83 constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
84 constexpr IntegrationType I =
85  IntegrationType::GAUSS; //< selected integration type
86 
87 template <int DIM>
89 
98 
102 
104 
106 using AssemblyBoundaryEleOp =
108 
117 
119  A>::LinearForm<I>::OpBaseTimesVector<BASE_DIM, SPACE_DIM, 1>;
121  A>::LinearForm<I>::OpBaseTimesScalar<BASE_DIM>;
123  A>::LinearForm<I>::OpBaseTimesScalar<BASE_DIM>;
124 
125 template <CoordinateTypes COORD_TYPE>
128 
129 // Flux is applied by Lagrange Multiplie
131 using OpFluidFlux =
133 
138 
140 
141 // mesh refinement
142 int order = 3; ///< approximation order
143 int nb_levels = 4; //< number of refinement levels
144 int refine_overlap = 4; //< mesh overlap while refine
145 
146 constexpr bool debug = true;
147 
148 auto get_start_bit = []() {
149  return nb_levels + 1;
150 }; //< first refinement level for computational mesh
151 auto get_current_bit = []() {
152  return 2 * get_start_bit() + 1;
153 }; ///< dofs bit used to do calculations
154 auto get_skin_parent_bit = []() { return 2 * get_start_bit() + 2; };
155 auto get_skin_child_bit = []() { return 2 * get_start_bit() + 3; };
156 auto get_projection_bit = []() { return 2 * get_start_bit() + 4; };
157 auto get_skin_projection_bit = []() { return 2 * get_start_bit() + 5; };
158 
159 // FIXME: Set parameters from command line
160 
161 // Physical parameters
162 double a0 = 980;
163 double rho_m = 0.998;
164 double mu_m = 0.010101 * 1e2;
165 double rho_p = 0.0012;
166 double mu_p = 0.000182 * 1e2;
167 double lambda = 73; ///< surface tension
168 double W = 0.25;
169 
170 // Model parameters
171 double h = 0.03; // mesh size
172 double eta = h;
173 double eta2 = eta * eta;
174 
175 // Numerical parameters
176 double md = 1e-2; // mobility
177 double eps = 1e-12;
178 double tol = std::numeric_limits<float>::epsilon();
179 
180 double rho_ave = (rho_p + rho_m) / 2;
181 double rho_diff = (rho_p - rho_m) / 2;
182 double mu_ave = (mu_p + mu_m) / 2;
183 double mu_diff = (mu_p - mu_m) / 2;
184 
185 double kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
186 
187 auto integration_rule = [](int, int, int) { return 2 * order + 1; };
188 
189 auto cylindrical = [](const double r) {
190  if (coord_type == CYLINDRICAL)
191  return 2 * M_PI * r;
192  else
193  return 1.;
194 };
195 
196 auto 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 
203 auto my_max = [](const double x) { return (x - 1 + std::abs(x + 1)) / 2; };
204 auto my_min = [](const double x) { return (x + 1 - std::abs(x - 1)) / 2; };
205 auto cut_off = [](const double h) { return my_max(my_min(h)); };
206 auto d_cut_off = [](const double h) {
207  if (h >= -1 && h < 1)
208  return 1.;
209  else
210  return 0.;
211 };
212 
213 auto phase_function = [](const double h, const double diff, const double ave) {
214  return diff * h + ave;
215 };
216 
217 auto d_phase_function_h = [](const double h, const double diff) {
218  return diff;
219 };
220 
221 // order function (free energy)
222 
223 auto get_f = [](const double h) { return 4 * W * h * (h * h - 1); };
224 auto get_f_dh = [](const double h) { return 4 * W * (3 * h * h - 1); };
225 
226 auto get_M0 = [](auto h) { return md; };
227 auto get_M0_dh = [](auto h) { return 0; };
228 
229 auto get_M2 = [](auto h) {
230  return md * (1 - h * h);
231 };
232 
233 auto get_M2_dh = [](auto h) { return -md * 2 * h; };
234 
235 auto 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 
244 auto 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 
251 auto get_M = [](auto h) { return get_M0(h); };
252 auto get_M_dh = [](auto h) { return get_M0_dh(h); };
253 
254 auto 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 
262 auto 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 
272 auto 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 
280 auto 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 
286 auto 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  */
297 auto 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 
312 auto wetting_angle = [](double water_level) { return water_level; };
313 
314 auto bit = [](auto b) { return BitRefLevel().set(b); };
315 auto marker = [](auto b) { return BitRefLevel().set(BITREFLEVEL_SIZE - b); };
316 auto get_fe_bit = [](FEMethod *fe_ptr) {
317  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
318 };
319 
320 auto 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 
326 auto 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  */
342 auto 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  */
370 auto 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>
392 using namespace FreeSurfaceOps;
393 
394 struct FreeSurface;
395 
396 enum 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  */
405 struct 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  */
416  MoFEMErrorCode tsSetUp(TS ts);
417 
418  SmartPetscObj<VecScatter> getScatter(Vec x, Vec y, enum FR fr);
419  SmartPetscObj<Vec> getSubVector();
420 
424 
425 private:
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 
446  static MoFEMErrorCode tsPreStage(TS ts);
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 
469 static boost::weak_ptr<TSPrePostProc> tsPrePostProc;
470 
471 struct FreeSurface {
472 
473  FreeSurface(MoFEM::Interface &m_field) : mField(m_field) {}
474 
475  MoFEMErrorCode runProblem();
476 
477  MoFEMErrorCode makeRefProblem();
478 
480 
481 private:
482  MoFEMErrorCode readMesh();
483  MoFEMErrorCode setupProblem();
484  MoFEMErrorCode boundaryCondition();
485  MoFEMErrorCode projectData();
486  MoFEMErrorCode assembleSystem();
487  MoFEMErrorCode solveSystem();
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
499  Range findParentsToRefine(Range ents, BitRefLevel level, BitRefLevel mask);
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
513  MoFEMErrorCode setParentDofs(
514  boost::shared_ptr<FEMethod> fe_top, std::string field_name,
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]
526  CHKERR readMesh();
527  CHKERR setupProblem();
528  CHKERR boundaryCondition();
529  CHKERR assembleSystem();
530  CHKERR solveSystem();
532 }
533 //! [Run programme]
534 
535 //! [Read mesh]
538  MOFEM_LOG("FS", Sev::inform) << "Read mesh for problem";
539  auto simple = mField.getInterface<Simple>();
540 
541  simple->getParentAdjacencies() = true;
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_NULL, "", "-order", &order, PETSC_NULL);
556  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-nb_levels", &nb_levels,
557  PETSC_NULL);
558  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-refine_overlap", &refine_overlap,
559  PETSC_NULL);
560 
561  const char *coord_type_names[] = {"cartesian", "polar", "cylindrical",
562  "spherical"};
563  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-coords", coord_type_names,
564  LAST_COORDINATE_SYSTEM, &coord_type, PETSC_NULL);
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 
571  auto simple = mField.getInterface<Simple>();
572 
573  CHKERR PetscOptionsGetScalar(PETSC_NULL, "Acceleration", "-a0", &a0, PETSC_NULL);
574  CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Minus\" phase density rho_m",
575  "-rho_m", &rho_m, PETSC_NULL);
576  CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Minus\" phase viscosity", "-mu_m",
577  &mu_m, PETSC_NULL);
578  CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Plus\" phase density", "-rho_p",
579  &rho_p, PETSC_NULL);
580  CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Plus\" phase viscosity", "-mu_p",
581  &mu_p, PETSC_NULL);
582  CHKERR PetscOptionsGetScalar(PETSC_NULL, "Surface tension", "-lambda",
583  &lambda, PETSC_NULL);
584  CHKERR PetscOptionsGetScalar(PETSC_NULL,
585  "Height of the well in energy functional", "-W",
586  &W, PETSC_NULL);
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_NULL, "", "-h", &h, PETSC_NULL);
593  eta = h;
594  eta2 = eta * eta;
595  kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
596 
597  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-md", &eta, PETSC_NULL);
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,
630  U_FIELD_DIM);
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 
676  auto simple = mField.getInterface<Simple>();
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 
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(
705  new DomianParentEle(mField));
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 
723  CHKERR setParentDofs(
724  fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
725 
726  [&]() {
727  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
728  new DomianParentEle(mField));
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,
873  void *))MoFEMSNESMonitorFields,
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_NULL, 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 
960  auto simple = mField.getInterface<Simple>();
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 
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(
986  new DomianParentEle(mField));
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(
1000  new BoundaryParentEle(mField));
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();
1048  CHKERR setParentDofs(
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,
1076  bit(get_projection_bit())));
1077 
1078  auto assemble_fe_ptr = boost::make_shared<DomianParentEle>(mField);
1079  {
1080  auto &pip = assemble_fe_ptr->getOpPtrVector();
1082  CHKERR setParentDofs(
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 
1120  CHKERR setParentDofs(
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();
1166  CHKERR setParentDofs(
1167  eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
1168 
1169  [&]() {
1170  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1171  new BoundaryParentEle(mField));
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,
1185  bit(get_projection_bit())));
1186 
1187  auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
1188  {
1189  auto &pip = assemble_fe_ptr->getOpPtrVector();
1190  CHKERR setParentDofs(
1191  assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1192 
1193  [&]() {
1194  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1195  new BoundaryParentEle(mField));
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 
1236  CHKERR setParentDofs(
1237  fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1238 
1239  [&]() {
1240  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1241  new BoundaryParentEle(mField));
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,
1300  mField.get_comm_size());
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");
1326  CHK_THROW_MESSAGE(DMMoFEMCreateMoFEM(dummy_dm, &mField,
1327  simple->getProblemName().c_str(),
1328  BitRefLevel(), BitRefLevel()),
1329  "create dummy dm");
1330  return dummy_dm;
1331  };
1332 
1333  auto subdm = create_subdm();
1334  if (subdm) {
1335 
1336  pip_mng->getDomainRhsFE().reset();
1337  pip_mng->getDomainLhsFE().reset();
1338  pip_mng->getBoundaryRhsFE().reset();
1339  pip_mng->getBoundaryLhsFE().reset();
1340  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
1341  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
1342  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
1343  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
1344  pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1345  pip_mng->getDomainRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
1346  return get_fe_bit(fe_ptr).test(nb_levels - 1);
1347  };
1348  pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
1349  pip_mng->getBoundaryRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
1350  return get_fe_bit(fe_ptr).test(nb_levels - 1);
1351  };
1352 
1353  CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1354  CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1355  CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
1356  CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
1357 
1358  auto D = createDMVector(subdm);
1359  auto F = vectorDuplicate(D);
1360 
1361  auto ksp = pip_mng->createKSP(subdm);
1362  CHKERR KSPSetFromOptions(ksp);
1363  CHKERR KSPSetUp(ksp);
1364 
1365  // get vector norm
1366  auto get_norm = [&](auto x) {
1367  double nrm;
1368  CHKERR VecNorm(x, NORM_2, &nrm);
1369  return nrm;
1370  };
1371 
1372  /**
1373  * @brief Zero DOFs, used by FieldBlas
1374  */
1375  auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
1377  for (auto &v : ent_ptr->getEntFieldData()) {
1378  v = 0;
1379  }
1381  };
1382 
1383  auto solve = [&](auto S) {
1385  CHKERR VecZeroEntries(S);
1386  CHKERR VecZeroEntries(F);
1387  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
1388  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
1389  CHKERR KSPSolve(ksp, F, S);
1390  CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
1391  CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
1392 
1393 
1394 
1396  };
1397 
1398  MOFEM_LOG("FS", Sev::inform) << "Solve projection";
1399  CHKERR solve(D);
1400 
1401  auto glob_x = createDMVector(simple->getDM());
1402  auto sub_x = createDMVector(subdm);
1403  auto dummy_dm = create_dummy_dm();
1404 
1405  /**
1406  * @brief get TSTheta data operators
1407  */
1408  auto apply_restrict = [&]() {
1409  auto get_is = [](auto v) {
1410  IS iy;
1411  auto create = [&]() {
1413  int n, ystart;
1414  CHKERR VecGetLocalSize(v, &n);
1415  CHKERR VecGetOwnershipRange(v, &ystart, NULL);
1416  CHKERR ISCreateStride(PETSC_COMM_SELF, n, ystart, 1, &iy);
1418  };
1419  CHK_THROW_MESSAGE(create(), "create is");
1420  return SmartPetscObj<IS>(iy);
1421  };
1422 
1423  auto iy = get_is(glob_x);
1424  auto s = createVecScatter(glob_x, PETSC_NULL, glob_x, iy);
1425 
1427  DMSubDomainRestrict(simple->getDM(), s, PETSC_NULL, dummy_dm),
1428  "restrict");
1429  Vec X0, Xdot;
1430  CHK_THROW_MESSAGE(DMGetNamedGlobalVector(dummy_dm, "TSTheta_X0", &X0),
1431  "get X0");
1433  DMGetNamedGlobalVector(dummy_dm, "TSTheta_Xdot", &Xdot),
1434  "get Xdot");
1435 
1436  auto forward_ghost = [](auto g) {
1438  CHKERR VecGhostUpdateBegin(g, INSERT_VALUES, SCATTER_FORWARD);
1439  CHKERR VecGhostUpdateEnd(g, INSERT_VALUES, SCATTER_FORWARD);
1441  };
1442 
1443  CHK_THROW_MESSAGE(forward_ghost(X0), "");
1444  CHK_THROW_MESSAGE(forward_ghost(Xdot), "");
1445 
1446  if constexpr (debug) {
1447  MOFEM_LOG("FS", Sev::inform)
1448  << "Reverse restrict: X0 " << get_norm(X0) << " Xdot "
1449  << get_norm(Xdot);
1450  }
1451 
1452  return std::vector<Vec>{X0, Xdot};
1453  };
1454 
1455  auto ts_solver_vecs = apply_restrict();
1456 
1457  if (ts_solver_vecs.size()) {
1458 
1459  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 
1504  CHKERR setParentDofs(
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 
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(
1606  new BoundaryParentEle(mField));
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(
1653  new OpCalculateVectorFieldValuesDot<U_FIELD_DIM>("U", dot_u_ptr));
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(
1808  new BoundaryParentEle(mField));
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(
1890  new BoundaryParentEle(mField));
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  */
1988 struct 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 
2034 private:
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 
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 
2403 int 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  }
2441  CATCH_ERRORS;
2442 
2444 
2445 #ifdef PYTHON_INIT_SURFACE
2446  if (Py_FinalizeEx() < 0) {
2447  exit(120);
2448  }
2449 #endif
2450 }
2451 
2452 std::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;
2504  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
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,
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 
2838  new OpAddParentEntData(
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 
2853  new OpAddParentEntData(
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>();
2888  CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
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 
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 
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 
3052 MoFEMErrorCode 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 
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_NULL, y, is, &s),
3115  "crate scatter");
3116  } else {
3117  CHK_THROW_MESSAGE(VecScatterCreate(x, is, y, PETSC_NULL, &s),
3118  "crate scatter");
3119  }
3120  return SmartPetscObj<VecScatter>(s);
3121  }
3122  }
3123  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "No prb pinter");
3124  return SmartPetscObj<VecScatter>();
3125 }
3126 
3128  return createDMVector(solverSubDM);
3129 }
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_NULL, PETSC_NULL, tsSetIJacobian, nullptr);
3143  CHKERR TSMonitorSet(ts, tsMonitor, fsRawPtr, PETSC_NULL);
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 }
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
get_D
auto get_D
Definition: free_surface.cpp:254
rho_ave
double rho_ave
Definition: free_surface.cpp:180
TSPrePostProc::solverSubDM
SmartPetscObj< DM > solverSubDM
Definition: free_surface.cpp:421
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:712
FreeSurface::findEntitiesCrossedByPhaseInterface
std::vector< Range > findEntitiesCrossedByPhaseInterface(size_t overlap)
Find entities on refinement levels.
Definition: free_surface.cpp:2453
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
get_fe_bit
auto get_fe_bit
Definition: free_surface.cpp:316
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: free_surface.cpp:134
g
constexpr double g
Definition: shallow_wave.cpp:63
FreeSurface::projectData
MoFEMErrorCode projectData()
[Boundary condition]
Definition: free_surface.cpp:953
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::debug
const static bool debug
Definition: ConstrainMatrixCtx.cpp:9
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
FreeSurfaceOps::OpNormalForceRhs
Definition: FreeSurfaceOps.hpp:92
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:24
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
kernel_oscillation
auto kernel_oscillation
Definition: free_surface.cpp:262
Monitor::liftFE
boost::shared_ptr< BoundaryEle > liftFE
Definition: free_surface.cpp:2039
CARTESIAN
@ CARTESIAN
Definition: definitions.h:128
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
get_M_dh
auto get_M_dh
Definition: free_surface.cpp:252
MoFEM::OpCalculateDivergenceVectorFieldValues
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:488
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
FreeSurface::FreeSurface
FreeSurface(MoFEM::Interface &m_field)
Definition: free_surface.cpp:473
FreeSurfaceOps::OpWettingAngleRhs
Definition: FreeSurfaceOps.hpp:138
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:373
get_skin_child_bit
auto get_skin_child_bit
Definition: free_surface.cpp:155
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
cut_off
auto cut_off
Definition: free_surface.cpp:205
get_M
auto get_M
Definition: free_surface.cpp:251
marker
auto marker
Definition: free_surface.cpp:315
MoFEM::OpRunParent
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Definition: MeshProjectionDataOperators.hpp:17
eps
double eps
Definition: free_surface.cpp:177
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition: PetscSmartObj.hpp:249
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::TsSetIFunction
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
Definition: TsCtx.cpp:56
get_skin_parent_bit
auto get_skin_parent_bit
Definition: free_surface.cpp:154
k
FTensor::Index< 'k', SPACE_DIM > k
Definition: free_surface.cpp:136
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
init_h
auto init_h
Initialisation function.
Definition: free_surface.cpp:297
get_f_dh
auto get_f_dh
Definition: free_surface.cpp:224
mu_m
double mu_m
Definition: free_surface.cpp:164
FreeSurfaceOps::OpWettingAngleLhs
Definition: FreeSurfaceOps.hpp:246
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FreeSurface::assembleSystem
MoFEMErrorCode assembleSystem()
[Data projection]
Definition: free_surface.cpp:1581
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468
MoFEM::getPetscObject
PetscObject getPetscObject(T obj)
Definition: PetscSmartObj.hpp:9
d_phase_function_h
auto d_phase_function_h
Definition: free_surface.cpp:217
my_min
auto my_min
Definition: free_surface.cpp:204
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: dynamic_first_order_con_law.cpp:52
md
double md
Definition: free_surface.cpp:176
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
MoFEM::ForcesAndSourcesCore::UserDataOperator::OpType
OpType
Controls loop over entities on element.
Definition: ForcesAndSourcesCore.hpp:566
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
FreeSurface::runProblem
MoFEMErrorCode runProblem()
[Run programme]
Definition: free_surface.cpp:524
MoFEM.hpp
MoFEM::DMoFEMMeshToLocalVector
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:523
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
get_f
auto get_f
Definition: free_surface.cpp:223
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:261
OpBoundaryAssembleScalar
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
Definition: free_surface.cpp:123
W
double W
Definition: free_surface.cpp:168
FreeSurfaceOps::OpLhsH_dG
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1117
l
FTensor::Index< 'l', SPACE_DIM > l
Definition: free_surface.cpp:137
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
test_bit_child
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
Definition: hanging_node_approx.cpp:173
Monitor::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: free_surface.cpp:1997
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
TSPrePostProc::getScatter
SmartPetscObj< VecScatter > getScatter(Vec x, Vec y, enum FR fr)
Definition: free_surface.cpp:3107
get_M3
auto get_M3
Definition: free_surface.cpp:235
Monitor::postProc
boost::shared_ptr< PostProcEleDomainCont > postProc
Definition: free_surface.cpp:2037
debug
constexpr bool debug
Definition: free_surface.cpp:146
sdf.r
int r
Definition: sdf.py:8
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1056
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
TSPrePostProc::subB
SmartPetscObj< Mat > subB
Definition: free_surface.cpp:460
capillary_tube
auto capillary_tube
Definition: free_surface.cpp:280
OpMixScalarTimesDiv
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
Definition: free_surface.cpp:127
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1319
MoFEM::CoreInterface::get_dofs
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
MoFEM::MeshsetsManager::getMeshset
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:719
h
double h
Definition: free_surface.cpp:171
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:263
FreeSurfaceOps::OpRhsH
Definition: FreeSurfaceOps.hpp:795
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:597
EntData
EntitiesFieldData::EntData EntData
Definition: free_surface.cpp:103
cylindrical
auto cylindrical
Definition: free_surface.cpp:189
OpDomainAssembleVector
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
Definition: free_surface.cpp:119
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:438
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:677
tol
double tol
Definition: free_surface.cpp:178
TSPrePostProc::subKSP
SmartPetscObj< KSP > subKSP
Definition: free_surface.cpp:461
U_FIELD_DIM
constexpr int U_FIELD_DIM
Definition: free_surface.cpp:81
OpBoundaryMassL
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
Definition: free_surface.cpp:116
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
TSPrePostProc
Set of functions called by PETSc solver used to refine and update mesh.
Definition: dynamic_first_order_con_law.cpp:383
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
kernel_eye
auto kernel_eye
Definition: free_surface.cpp:272
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FreeSurfaceOps::OpLhsH_dU
Lhs for H dU.
Definition: FreeSurfaceOps.hpp:915
coord_type
int coord_type
Definition: free_surface.cpp:77
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
TSPrePostProc::tsMonitor
static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Wrapper for TS monitor.
Definition: free_surface.cpp:3052
BASE_DIM
constexpr int BASE_DIM
Definition: free_surface.cpp:79
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
wetting_angle
auto wetting_angle
Definition: free_surface.cpp:312
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
FreeSurfaceOps::OpLoopSideGetDataForSideEle
Definition: FreeSurfaceOps.hpp:576
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
BoundaryEleOp
save_range
auto save_range
Definition: free_surface.cpp:326
R
@ R
Definition: free_surface.cpp:396
MoFEM::BitRefManager::getEntitiesByDimAndRefLevel
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
Definition: BitRefManager.cpp:822
rho_diff
double rho_diff
Definition: free_surface.cpp:181
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1886
get_current_bit
auto get_current_bit
dofs bit used to do calculations
Definition: free_surface.cpp:151
j
FTensor::Index< 'j', SPACE_DIM > j
Definition: free_surface.cpp:135
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::DMMoFEMCreateSubDM
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
double
convert.type
type
Definition: convert.py:64
OpDomainMassH
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
Definition: free_surface.cpp:112
get_global_size
auto get_global_size
Definition: free_surface.cpp:320
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
FreeSurfaceOps::OpLhsH_dH
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:973
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
nb_levels
int nb_levels
Definition: free_surface.cpp:143
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
eta2
double eta2
Definition: free_surface.cpp:173
TSPrePostProc::tsSetIFunction
static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec f, void *ctx)
Definition: free_surface.cpp:2997
TSPrePostProc::tsSetUp
MoFEMErrorCode tsSetUp(TS ts)
Used to setup TS solver.
Definition: free_surface.cpp:3131
MoFEM::PostProcBrokenMeshInMoabBaseContImpl
Definition: PostProcBrokenMeshInMoabBase.hpp:890
BoundaryParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
Definition: free_surface.cpp:94
MoFEM::FieldEntity::getLoBitNumberUId
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:222
OpDomainMassP
OpDomainMassH OpDomainMassP
Definition: free_surface.cpp:113
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
get_projection_bit
auto get_projection_bit
Definition: free_surface.cpp:156
get_M0_dh
auto get_M0_dh
Definition: free_surface.cpp:227
TSPrePostProc::tsPostProc
static MoFEMErrorCode tsPostProc(TS ts)
Post process time step.
Definition: free_surface.cpp:2987
MoFEM::createPC
auto createPC(MPI_Comm comm)
Definition: PetscSmartObj.hpp:267
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
help
static char help[]
Definition: free_surface.cpp:17
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
eta
double eta
Definition: free_surface.cpp:172
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
main
int main(int argc, char *argv[])
Definition: free_surface.cpp:2403
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
bit
auto bit
Definition: free_surface.cpp:314
SideOp
SideEle::UserDataOperator SideOp
Definition: free_surface.cpp:97
FreeSurfaceOps::OpRhsG
Definition: FreeSurfaceOps.hpp:1181
phase_function
auto phase_function
Definition: free_surface.cpp:213
wetting_angle_sub_stepping
auto wetting_angle_sub_stepping
Definition: free_surface.cpp:196
SPACE_DIM
constexpr int SPACE_DIM
Definition: free_surface.cpp:80
MoFEM::solve
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
Definition: Schur.cpp:1296
Monitor::postProcMesh
boost::shared_ptr< moab::Core > postProcMesh
Definition: free_surface.cpp:2036
MoFEM::TsSetIJacobian
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
Monitor::dM
SmartPetscObj< DM > dM
Definition: dynamic_first_order_con_law.cpp:873
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
MoFEM::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
MoFEM::DMMoFEMCreateMoFEM
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
MoFEM::ForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: ForcesAndSourcesCore.hpp:482
TSPrePostProc::snesCtxPtr
boost::shared_ptr< SnesCtx > snesCtxPtr
Definition: free_surface.cpp:464
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:139
FreeSurfaceOps::OpLhsU_dG
Lhs for G dH.
Definition: FreeSurfaceOps.hpp:732
COORD_TYPE
constexpr CoordinateTypes COORD_TYPE
Definition: tensor_divergence_operator.cpp:36
BiLinearForm
get_skin_projection_bit
auto get_skin_projection_bit
Definition: free_surface.cpp:157
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
order
int order
approximation order
Definition: free_surface.cpp:142
FTensor::Index< 'i', SPACE_DIM >
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
mu_diff
double mu_diff
Definition: free_surface.cpp:183
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:417
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1561
convert.n
n
Definition: convert.py:82
rho_p
double rho_p
Definition: free_surface.cpp:165
FreeSurfaceOps::OpLhsG_dG
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1321
get_M3_dh
auto get_M3_dh
Definition: free_surface.cpp:244
get_M2_dh
auto get_M2_dh
Definition: free_surface.cpp:233
kappa
double kappa
Definition: free_surface.cpp:185
MoFEM::TSMethod::ts_step
PetscInt ts_step
time step number
Definition: LoopMethods.hpp:163
FreeSurfaceOps.hpp
rho_m
double rho_m
Definition: free_surface.cpp:163
FreeSurface
Definition: free_surface.cpp:471
integration_rule
auto integration_rule
Definition: free_surface.cpp:187
tsPrePostProc
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
Definition: free_surface.cpp:469
TSPrePostProc::tsPreProc
static MoFEMErrorCode tsPreProc(TS ts)
Pre process time step.
Definition: free_surface.cpp:2872
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ElementsAndOps
Definition: child_and_parent.cpp:18
Range
FreeSurfaceOps
Definition: FreeSurfaceOps.hpp:1
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
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
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
get_dofs_ents_all
auto get_dofs_ents_all
get entities of dofs in the problem - used for debugging
Definition: free_surface.cpp:370
Monitor::Monitor
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)
Definition: free_surface.cpp:1989
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
TSPrePostProc::tsCtxPtr
boost::shared_ptr< TsCtx > tsCtxPtr
Definition: free_surface.cpp:466
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CYLINDRICAL
@ CYLINDRICAL
Definition: definitions.h:130
lambda
double lambda
surface tension
Definition: free_surface.cpp:167
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::TSMethod::ts_t
PetscReal ts_t
time
Definition: LoopMethods.hpp:166
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
FreeSurface::setParentDofs
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.
Definition: free_surface.cpp:2806
FreeSurfaceOps::OpLhsU_dU
Lhs for U dU.
Definition: FreeSurfaceOps.hpp:473
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
mu_ave
double mu_ave
Definition: free_surface.cpp:182
FreeSurfaceOps::OpNormalConstrainRhs
Definition: FreeSurfaceOps.hpp:50
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:414
TSPrePostProc::tsSetIJacobian
static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Wrapper for SNES Lhs.
Definition: free_surface.cpp:3025
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
OpDomainMassU
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
Definition: free_surface.cpp:110
refine_overlap
int refine_overlap
Definition: free_surface.cpp:144
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
BITREFLEVEL_SIZE
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:232
Monitor
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:774
FTensor::Ddg
Definition: Ddg_value.hpp:7
MoFEM::createVecScatter
auto createVecScatter(Vec x, IS ix, Vec y, IS iy)
Create a Vec Scatter object.
Definition: PetscSmartObj.hpp:361
MoFEM::DMMoFEMTSSetMonitor
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:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
FreeSurface::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: free_surface.cpp:2044
MoFEM::getProblemPtr
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:1044
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:64
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
get_start_bit
auto get_start_bit
Definition: free_surface.cpp:148
FreeSurface::readMesh
MoFEMErrorCode readMesh()
[Run programme]
Definition: free_surface.cpp:536
LAST_COORDINATE_SYSTEM
@ LAST_COORDINATE_SYSTEM
Definition: definitions.h:132
d_cut_off
auto d_cut_off
Definition: free_surface.cpp:206
get_M2
auto get_M2
Definition: free_surface.cpp:229
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
Monitor::postProc
boost::shared_ptr< PostProcEle > postProc
Definition: dynamic_first_order_con_law.cpp:874
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
FreeSurfaceOps::OpRhsU
Rhs for U.
Definition: FreeSurfaceOps.hpp:356
FreeSurfaceOps::OpCalculateLift
Definition: FreeSurfaceOps.hpp:3
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
BoundaryEleOp
BoundaryEle::UserDataOperator BoundaryEleOp
Definition: free_surface.cpp:95
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:239
MoFEM::FieldEntity::getHiBitNumberUId
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:228
TSPrePostProc::globSol
SmartPetscObj< Vec > globSol
Definition: free_surface.cpp:422
TSPrePostProc::getSubVector
SmartPetscObj< Vec > getSubVector()
Definition: free_surface.cpp:3127
FreeSurfaceOps::OpLhsU_dH
Lhs for U dH.
Definition: FreeSurfaceOps.hpp:611
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
a0
double a0
Definition: free_surface.cpp:162
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
Example::mField
MoFEM::Interface & mField
Definition: plastic.cpp:226
TSPrePostProc::pcSetup
static MoFEMErrorCode pcSetup(PC pc)
Definition: free_surface.cpp:3077
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1127
MoFEM::FEMethod::numeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Definition: LoopMethods.hpp:387
BoundaryParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
Definition: child_and_parent.cpp:41
FreeSurface::mField
MoFEM::Interface & mField
Definition: free_surface.cpp:479
QUIET
@ QUIET
Definition: definitions.h:221
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:18
mu_p
double mu_p
Definition: free_surface.cpp:166
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
FreeSurfaceOps::OpLhsG_dH
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1249
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::BasicMethod::getCacheWeakPtr
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
Definition: LoopMethods.hpp:348
Monitor::postProcEdge
boost::shared_ptr< PostProcEleBdyCont > postProcEdge
Definition: free_surface.cpp:2038
MoFEM::SmartPetscObj< VecScatter >
FreeSurface::refineMesh
MoFEMErrorCode refineMesh(size_t overlap)
Definition: free_surface.cpp:2575
Monitor::liftVec
boost::shared_ptr< VectorDouble > liftVec
Definition: free_surface.cpp:2040
get_M0
auto get_M0
Definition: free_surface.cpp:226
FreeSurfaceOps::OpNormalConstrainLhs
Definition: FreeSurfaceOps.hpp:195
FreeSurface::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: free_surface.cpp:552
OpDomainAssembleScalar
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
Definition: free_surface.cpp:121
I
constexpr IntegrationType I
Definition: free_surface.cpp:84
MoFEM::DMoFEMLoopFiniteElements
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:586
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
TSPrePostProc::globRes
SmartPetscObj< Vec > globRes
Definition: free_surface.cpp:459
MoFEM::OpAddParentEntData
Operator to project base functions from parent entity to child.
Definition: MeshProjectionDataOperators.hpp:66
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
bubble_device
auto bubble_device
Definition: free_surface.cpp:286
DomianParentEle
ElementsAndOps< SPACE_DIM >::DomianParentEle DomianParentEle
Definition: free_surface.cpp:91
get_dofs_ents_by_field_name
auto get_dofs_ents_by_field_name
get entities of dofs in the problem - used for debugging
Definition: free_surface.cpp:342
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
my_max
auto my_max
Definition: free_surface.cpp:203
DomianParentEle
PipelineManager::ElementsAndOpsByDim< FE_DIM >::DomianParentEle DomianParentEle
Definition: level_set.cpp:39
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1291
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
A
constexpr AssemblyType A
Definition: free_surface.cpp:83
TSPrePostProc::fsRawPtr
Example * fsRawPtr
Definition: dynamic_first_order_con_law.cpp:397
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:711
OpDomainMassG
OpDomainMassH OpDomainMassG
Definition: free_surface.cpp:114
TSPrePostProc::fsRawPtr
FreeSurface * fsRawPtr
Definition: free_surface.cpp:423
TSPrePostProc::pcApply
static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x)
Definition: free_surface.cpp:3087
FreeSurface::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: free_surface.cpp:663
F
@ F
Definition: free_surface.cpp:396
FR
FR
Definition: free_surface.cpp:396
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:709