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