v0.13.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 
10 /* This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #include <MoFEM.hpp>
25 
26 using namespace MoFEM;
27 
28 static char help[] = "...\n\n";
29 
30 #include <BasicFiniteElements.hpp>
31 
32 constexpr int BASE_DIM = 1;
33 constexpr int SPACE_DIM = 2;
34 constexpr int U_FIELD_DIM = SPACE_DIM;
36  CARTESIAN; ///< select coordinate system <CARTESIAN, CYLINDRICAL>;
37 
38 template <int DIM> struct ElementsAndOps {};
39 
40 template <> struct ElementsAndOps<2> {
46 };
47 
48 template <> struct ElementsAndOps<3> {
54 };
55 
60 
62 
68 
74 
79 
82 
85 
90 
92 
93 // Physical parameters
94 constexpr double a0 = 0.98;
95 constexpr double rho_m = 0.998;
96 constexpr double mu_m = 0.0101;
97 constexpr double rho_p = 0.0012;
98 constexpr double mu_p = 0.000182;
99 constexpr double lambda = 7.4;
100 constexpr double W = 0.25;
101 
102 // Model parameters
103 constexpr double h = 0.02; // mesh size
104 constexpr double eta = h;
105 constexpr double eta2 = eta * eta;
106 
107 // Numerical parameteres
108 constexpr double md = 1e-2;
109 constexpr double eps = 1e-12;
110 constexpr double tol = std::numeric_limits<float>::epsilon();
111 
112 constexpr double rho_ave = (rho_p + rho_m) / 2;
113 constexpr double rho_diff = (rho_p - rho_m) / 2;
114 constexpr double mu_ave = (mu_p + mu_m) / 2;
115 constexpr double mu_diff = (mu_p - mu_m) / 2;
116 
117 const double kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
118 
120  return 2 * approx_order;
121 };
122 
123 auto cylindrical = [](const double r) {
124  // When we move to C++17 add if constexpr()
125  if constexpr (coord_type == CYLINDRICAL)
126  return 2 * M_PI * r;
127  else
128  return 1.;
129 };
130 
131 auto my_max = [](const double x) { return (x - 1 + std::abs(x + 1)) / 2; };
132 auto my_min = [](const double x) { return (x + 1 - std::abs(x - 1)) / 2; };
133 auto cut_off = [](const double h) {
134  return my_max(my_min(h));
135 };
136 auto d_cut_off = [](const double h) {
137  if (h >= -1 && h < 1)
138  return 1.;
139  else
140  return 0.;
141 };
142 
143 auto phase_function = [](const double h, const double diff, const double ave) {
144  return diff * cut_off(h) + ave;
145 };
146 
147 auto d_phase_function_h = [](const double h, const double diff) {
148  return diff * d_cut_off(h);
149 };
150 
151 auto get_f = [](const double h) { return 4 * W * h * (h * h - 1); };
152 auto get_f_dh = [](const double h) { return 4 * W * (3 * h * h - 1); };
153 
154 auto get_M0 = [](auto h) { return md; };
155 auto get_M0_dh = [](auto h) { return 0; };
156 
157 auto get_M2 = [](auto h_tmp) {
158  const double h = cut_off(h_tmp);
159  return md * (1 - h * h);
160 };
161 
162 auto get_M2_dh = [](auto h_tmp) {
163  const double h = cut_off(h_tmp);
164  return -md * 2 * h * d_cut_off(h_tmp);
165 };
166 
167 auto get_M3 = [](auto h_tmp) {
168  const double h = cut_off(h_tmp);
169  const double h2 = h * h;
170  const double h3 = h2 * h;
171  if (h >= 0)
172  return md * (2 * h3 - 3 * h2 + 1);
173  else
174  return md * (-2 * h3 - 3 * h2 + 1);
175 };
176 
177 auto get_M3_dh = [](auto h_tmp) {
178  const double h = cut_off(h_tmp);
179  if (h >= 0)
180  return md * (6 * h * (h - 1)) * d_cut_off(h_tmp);
181  else
182  return md * (-6 * h * (h + 1)) * d_cut_off(h_tmp);
183 };
184 
185 auto get_M = [](auto h) { return get_M0(h); };
186 auto get_M_dh = [](auto h) { return get_M0_dh(h); };
187 
188 auto get_D = [](const double A) {
190  t_D(i, j, k, l) = A * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
191  return t_D;
192 };
193 
194 auto kernel_oscillation = [](double r, double y, double) {
195  constexpr int n = 3;
196  constexpr double R = 0.0125;
197  constexpr double A = R * 0.2;
198  const double theta = atan2(r, y);
199  const double w = R + A * cos(n * theta);
200  const double d = std::sqrt(r * r + y * y);
201  return tanh((w - d) / (eta * std::sqrt(2)));
202 };
203 
204 auto kernel_eye = [](double r, double y, double) {
205  constexpr double y0 = 0.5;
206  constexpr double R = 0.4;
207  y -= y0;
208  const double d = std::sqrt(r * r + y * y);
209  return tanh((R - d) / (eta * std::sqrt(2)));
210 };
211 
212 auto init_h = [](double r, double y, double theta) {
213  return kernel_eye(r, y, theta);
214 };
215 
216 #include <FreeSurfaceOps.hpp>
217 using namespace FreeSurfaceOps;
218 
219 struct FreeSurface {
220 
221  FreeSurface(MoFEM::Interface &m_field) : mField(m_field) {}
222 
223  MoFEMErrorCode runProblem();
224 
225 private:
226  MoFEMErrorCode readMesh();
227  MoFEMErrorCode setupProblem();
228  MoFEMErrorCode boundaryCondition();
229  MoFEMErrorCode assembleSystem();
230  MoFEMErrorCode solveSystem();
231 
233 
234  boost::shared_ptr<FEMethod> domianLhsFEPtr;
235 };
236 
237 //! [Run programme]
240  CHKERR readMesh();
241  CHKERR setupProblem();
242  CHKERR boundaryCondition();
243  CHKERR assembleSystem();
244  CHKERR solveSystem();
246 }
247 //! [Run programme]
248 
249 //! [Read mesh]
252 
253  auto simple = mField.getInterface<Simple>();
254  CHKERR simple->getOptions();
255  CHKERR simple->loadFile();
256 
258 }
259 //! [Read mesh]
260 
261 //! [Set up problem]
264 
265  auto simple = mField.getInterface<Simple>();
266 
267  // Fields on domain
268 
269  // Velocity field
270  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, U_FIELD_DIM);
271  // Pressure field
272  CHKERR simple->addDomainField("P", H1, AINSWORTH_LEGENDRE_BASE, 1);
273  // Order/phase fild
274  CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
275  //Chemical potential
276  CHKERR simple->addDomainField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
277 
278  // Field on boundary
279  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, U_FIELD_DIM);
280  CHKERR simple->addBoundaryField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
281  // Lagrange multiplier which constrains slip conditions
282  CHKERR simple->addBoundaryField("L", H1, AINSWORTH_LEGENDRE_BASE, 1);
283 
284  constexpr int order = 3;
285  CHKERR simple->setFieldOrder("U", order);
286  CHKERR simple->setFieldOrder("P", order - 1);
287  CHKERR simple->setFieldOrder("H", order);
288  CHKERR simple->setFieldOrder("G", order);
289  CHKERR simple->setFieldOrder("L", order);
290  CHKERR simple->setUp();
291 
293 }
294 //! [Set up problem]
295 
296 //! [Boundary condition]
299 
300  auto simple = mField.getInterface<Simple>();
301  auto pipeline_mng = mField.getInterface<PipelineManager>();
302  auto bc_mng = mField.getInterface<BcManager>();
303  auto dm = simple->getDM();
304 
305  auto h_ptr = boost::make_shared<VectorDouble>();
306  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
307  auto g_ptr = boost::make_shared<VectorDouble>();
308  auto grad_g_ptr = boost::make_shared<MatrixDouble>();
309 
310  auto post_proc = [&]() {
312  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
313  post_proc_fe->generateReferenceElementMesh();
314 
315  auto det_ptr = boost::make_shared<VectorDouble>();
316  auto jac_ptr = boost::make_shared<MatrixDouble>();
317  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
318 
319  post_proc_fe->addFieldValuesPostProc("H");
320  post_proc_fe->addFieldValuesPostProc("G");
321  post_proc_fe->addFieldValuesGradientPostProc("G", 2);
322  post_proc_fe->addFieldValuesGradientPostProc("H", 2);
323 
324  CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
325  CHKERR post_proc_fe->writeFile("out_init.h5m");
326 
328  };
329 
330  auto solve_init = [&]() {
332 
333  auto set_generic = [&](auto &pipeline) {
334  auto det_ptr = boost::make_shared<VectorDouble>();
335  auto jac_ptr = boost::make_shared<MatrixDouble>();
336  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
337  pipeline.push_back(new OpSetHOWeightsOnFace());
338  pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
339  pipeline.push_back(
340  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
341  pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
342 
343  pipeline.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
344  pipeline.push_back(
345  new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
346 
347  pipeline.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
348  pipeline.push_back(
349  new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
350  };
351 
352  auto set_domain_rhs = [&](auto &pipeline) {
353  pipeline.push_back(new OpRhsH<true>("H", nullptr, nullptr, h_ptr,
354  grad_h_ptr, grad_g_ptr));
355  pipeline.push_back(new OpRhsG<true>("G", h_ptr, grad_h_ptr, g_ptr));
356  };
357 
358  auto set_domain_lhs = [&](auto &pipeline) {
359  pipeline.push_back(new OpLhsH_dH<true>("H", nullptr, h_ptr, grad_g_ptr));
360  pipeline.push_back(new OpLhsH_dG<true>("H", "G", h_ptr));
361  pipeline.push_back(new OpLhsG_dH<true>("G", "H", h_ptr));
362  pipeline.push_back(new OpLhsG_dG("G"));
363  };
364 
365  auto create_subdm = [&]() {
366  DM subdm;
367  CHKERR DMCreate(mField.get_comm(), &subdm);
368  CHKERR DMSetType(subdm, "DMMOFEM");
369  CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB");
370  CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName().c_str());
371  CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
372  CHKERR DMMoFEMAddSubFieldRow(subdm, "H");
373  CHKERR DMMoFEMAddSubFieldRow(subdm, "G");
374  CHKERR DMMoFEMAddSubFieldCol(subdm, "H");
375  CHKERR DMMoFEMAddSubFieldCol(subdm, "G");
376  CHKERR DMSetUp(subdm);
377  return SmartPetscObj<DM>(subdm);
378  };
379 
380  auto subdm = create_subdm();
381  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
382  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
383 
384  set_generic(pipeline_mng->getOpDomainRhsPipeline());
385  set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
386  set_generic(pipeline_mng->getOpDomainLhsPipeline());
387  set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
388 
389  auto D = smartCreateDMVector(subdm);
390  auto snes = pipeline_mng->createSNES(subdm);
391 
392  auto set_section_monitor = [&](auto solver) {
394  PetscViewerAndFormat *vf;
395  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
396  PETSC_VIEWER_DEFAULT, &vf);
397  CHKERR SNESMonitorSet(
398  solver,
399  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
400  void *))SNESMonitorFields,
401  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
402  auto section = mField.getInterface<ISManager>()->sectionCreate("SUB");
403  PetscInt num_fields;
404  CHKERR PetscSectionGetNumFields(section, &num_fields);
405  for (int f = 0; f < num_fields; ++f) {
406  const char *field_name;
407  CHKERR PetscSectionGetFieldName(section, f, &field_name);
408  MOFEM_LOG("FS", Sev::inform)
409  << "Field " << f << " " << std::string(field_name);
410  }
412  };
413 
414  CHKERR set_section_monitor(snes);
415 
416  CHKERR SNESSetFromOptions(snes);
417  CHKERR SNESSolve(snes, PETSC_NULL, D);
418 
419  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
420  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
421  CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
422 
424  };
425 
426  CHKERR solve_init();
427  CHKERR post_proc();
428 
429  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
430  "SYMETRY", "U", 0, 0);
431  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
432  "SYMETRY", "L", 0, 0);
433  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX",
434  "U", 0, SPACE_DIM);
435  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX",
436  "L", 0, 0);
437  // Clear pipelines
438  pipeline_mng->getOpDomainRhsPipeline().clear();
439  pipeline_mng->getOpDomainLhsPipeline().clear();
440 
442 }
443 //! [Boundary condition]
444 
445 //! [Push operators to pipeline]
448 
449  auto dot_u_ptr = boost::make_shared<MatrixDouble>();
450  auto u_ptr = boost::make_shared<MatrixDouble>();
451  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
452  auto dot_h_ptr = boost::make_shared<VectorDouble>();
453  auto h_ptr = boost::make_shared<VectorDouble>();
454  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
455  auto g_ptr = boost::make_shared<VectorDouble>();
456  auto grad_g_ptr = boost::make_shared<MatrixDouble>();
457  auto lambda_ptr = boost::make_shared<VectorDouble>();
458  auto p_ptr = boost::make_shared<VectorDouble>();
459  auto div_u_ptr = boost::make_shared<VectorDouble>();
460 
461  // Push element from reference configuration to current configuration in 3d
462  // space
463  auto set_domain_general = [&](auto &pipeline) {
464  auto det_ptr = boost::make_shared<VectorDouble>();
465  auto jac_ptr = boost::make_shared<MatrixDouble>();
466  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
467  pipeline.push_back(new OpSetHOWeightsOnFace());
468  pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
469  pipeline.push_back(
470  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
471  pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
472 
473  pipeline.push_back(
475  pipeline.push_back(
477  pipeline.push_back(
479  grad_u_ptr));
480 
481  pipeline.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
482  pipeline.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
483  pipeline.push_back(
484  new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
485  pipeline.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
486  pipeline.push_back(
487  new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
488  pipeline.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
489  pipeline.push_back(
491  "U", div_u_ptr));
492  };
493 
494  auto set_domain_rhs = [&](auto &pipeline) {
495  pipeline.push_back(new OpRhsU("U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
496  grad_h_ptr, g_ptr, p_ptr));
497  pipeline.push_back(new OpRhsH<false>("H", u_ptr, dot_h_ptr, h_ptr,
498  grad_h_ptr, grad_g_ptr));
499  pipeline.push_back(new OpRhsG<false>("G", h_ptr, grad_h_ptr, g_ptr));
500  pipeline.push_back(new OpBaseTimesScalarField(
501  "P", div_u_ptr, [](const double r, const double, const double) {
502  return cylindrical(r);
503  }));
504  pipeline.push_back(new OpBaseTimesScalarField(
505  "P", p_ptr, [](const double r, const double, const double) {
506  return eps * cylindrical(r) ;
507  }));
508  };
509 
510  auto set_domain_lhs = [&](auto &pipeline) {
511  pipeline.push_back(new OpLhsU_dU("U", u_ptr, grad_u_ptr, h_ptr));
512  pipeline.push_back(
513  new OpLhsU_dH("U", "H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
514  pipeline.push_back(new OpLhsU_dG("U", "G", grad_h_ptr));
515 
516  pipeline.push_back(new OpLhsH_dU("H", "U", grad_h_ptr));
517  pipeline.push_back(new OpLhsH_dH<false>("H", u_ptr, h_ptr, grad_g_ptr));
518  pipeline.push_back(new OpLhsH_dG<false>("H", "G", h_ptr));
519 
520  pipeline.push_back(new OpLhsG_dH<false>("G", "H", h_ptr));
521  pipeline.push_back(new OpLhsG_dG("G"));
522 
523  pipeline.push_back(new OpMixScalarTimesDiv(
524  "P", "U",
525  [](const double r, const double, const double) {
526  return cylindrical(r);
527  },
528  true, false));
529  pipeline.push_back(
530  new OpDomainMassP("P", "P", [](double r, double, double) {
531  return eps * cylindrical(r);
532  }));
533  };
534 
535  auto set_boundary_rhs = [&](auto &pipeline) {
536  pipeline.push_back(
538  pipeline.push_back(new OpCalculateScalarFieldValues("L", lambda_ptr));
539  pipeline.push_back(new OpNormalConstrainRhs("L", u_ptr));
540  pipeline.push_back(new OpNormalForcebRhs("U", lambda_ptr));
541  };
542 
543  auto set_boundary_lhs = [&](auto &pipeline) {
544  pipeline.push_back(new OpNormalConstrainLhs("L", "U"));
545  };
546 
547  auto *pipeline_mng = mField.getInterface<PipelineManager>();
548 
549  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
550  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
551  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
552  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
553 
554  set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
555  set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
556  set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
557  set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
558  set_boundary_rhs(pipeline_mng->getOpBoundaryRhsPipeline());
559  set_boundary_lhs(pipeline_mng->getOpBoundaryLhsPipeline());
560 
561  domianLhsFEPtr = pipeline_mng->getDomainLhsFE();
562 
564 }
565 //! [Push operators to pipeline]
566 
567 /**
568  * @brief Monitor solution
569  *
570  * This functions is called by TS solver at the end of each step. It is used
571  * to output results to the hard drive.
572  */
573 struct Monitor : public FEMethod {
575  SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc,
576  std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
577  p)
578  : dM(dm), postProc(post_proc), liftFE(p.first), liftVec(p.second) {}
581  constexpr int save_every_nth_step = 1;
582  if (ts_step % save_every_nth_step == 0) {
583  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc,
584  this->getCacheWeakPtr());
585  CHKERR postProc->writeFile(
586  "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
587 
588 
589  // MOFEM_LOG("FS", Sev::verbose)
590  // << "writing vector in binary to vector.dat ...";
591  // PetscViewer viewer;
592  // PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_WRITE,
593  // &viewer);
594  // VecView(ts_u, viewer);
595  // PetscViewerDestroy(&viewer);
596  }
597 
598  liftVec->resize(SPACE_DIM, false);
599  liftVec->clear();
600  CHKERR DMoFEMLoopFiniteElements(dM, "bFE", liftFE, this->getCacheWeakPtr());
601  MPI_Allreduce(MPI_IN_PLACE, &(*liftVec)[0], SPACE_DIM, MPI_DOUBLE, MPI_SUM,
602  MPI_COMM_WORLD);
603  MOFEM_LOG("FS", Sev::inform)
604  << "Step " << ts_step << " time " << ts_t
605  << " lift vec x: " << (*liftVec)[0] << " y: " << (*liftVec)[1];
606 
608  }
609 
610 private:
612  boost::shared_ptr<PostProcEle> postProc;
613  boost::shared_ptr<BoundaryEle> liftFE;
614  boost::shared_ptr<VectorDouble> liftVec;
615 };
616 
617 //! [Solve]
620 
621  auto *simple = mField.getInterface<Simple>();
622  auto *pipeline_mng = mField.getInterface<PipelineManager>();
623  auto dm = simple->getDM();
624 
625  auto get_fe_post_proc = [&]() {
626  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
627  post_proc_fe->generateReferenceElementMesh();
628  auto det_ptr = boost::make_shared<VectorDouble>();
629  auto jac_ptr = boost::make_shared<MatrixDouble>();
630  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
631 
632  post_proc_fe->getOpPtrVector().push_back(
633  new OpCalculateHOJacForFace(jac_ptr));
634  post_proc_fe->getOpPtrVector().push_back(
635  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
636  post_proc_fe->getOpPtrVector().push_back(
637  new OpSetInvJacH1ForFace(inv_jac_ptr));
638 
639  post_proc_fe->addFieldValuesPostProc("U");
640  post_proc_fe->addFieldValuesPostProc("H");
641  post_proc_fe->addFieldValuesPostProc("P");
642  post_proc_fe->addFieldValuesPostProc("G");
643  post_proc_fe->addFieldValuesGradientPostProc("U", 2);
644  post_proc_fe->addFieldValuesGradientPostProc("H", 2);
645  post_proc_fe->addFieldValuesGradientPostProc("G", 2);
646  return post_proc_fe;
647  };
648 
649  auto get_lift_fe = [&]() {
650  auto fe = boost::make_shared<BoundaryEle>(mField);
651  auto lift_ptr = boost::make_shared<VectorDouble>();
652  auto p_ptr = boost::make_shared<VectorDouble>();
653  auto ents_ptr = boost::make_shared<Range>();
654 
655  std::vector<const CubitMeshSets *> vec_ptr;
656  CHKERR mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
657  std::regex("LIFT"), vec_ptr);
658  for (auto m_ptr : vec_ptr) {
659  auto meshset = m_ptr->getMeshset();
660  Range ents;
661  CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
662  ents, true);
663  ents_ptr->merge(ents);
664  }
665 
666  MOFEM_LOG("FS", Sev::inform) << "Lift ents " << (*ents_ptr);
667 
668  fe->getOpPtrVector().push_back(
669  new OpCalculateScalarFieldValues("P", p_ptr));
670  fe->getOpPtrVector().push_back(
671  new OpCalculateLift("P", p_ptr, lift_ptr, ents_ptr));
672 
673  return std::make_pair(fe, lift_ptr);
674  };
675 
676  auto set_ts = [&](auto solver) {
678  SNES snes;
679  CHKERR TSGetSNES(solver, &snes);
680  KSP ksp;
681  CHKERR SNESGetKSP(snes, &ksp);
683  };
684 
685  auto ts = pipeline_mng->createTSIM();
686  CHKERR TSSetType(ts, TSALPHA);
687 
688  auto set_post_proc_monitor = [&](auto dm) {
690  boost::shared_ptr<FEMethod> null_fe;
691  auto monitor_ptr =
692  boost::make_shared<Monitor>(dm, get_fe_post_proc(), get_lift_fe());
693  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
694  null_fe, monitor_ptr);
696  };
697  CHKERR set_post_proc_monitor(dm);
698 
699  auto set_section_monitor = [&](auto solver) {
701  SNES snes;
702  CHKERR TSGetSNES(solver, &snes);
703  PetscViewerAndFormat *vf;
704  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
705  PETSC_VIEWER_DEFAULT, &vf);
706  CHKERR SNESMonitorSet(
707  snes,
708  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
709  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
710 
711  auto section = mField.getInterface<ISManager>()->sectionCreate(
712  simple->getProblemName());
713  PetscInt num_fields;
714  CHKERR PetscSectionGetNumFields(section, &num_fields);
715  for (int f = 0; f < num_fields; ++f) {
716  const char *field_name;
717  CHKERR PetscSectionGetFieldName(section, f, &field_name);
718  MOFEM_LOG("FS", Sev::inform)
719  << "Field " << f << " " << std::string(field_name);
720  }
721 
723  };
724 
725  // Add monitor to time solver
726  double ftime = 1;
727  // CHKERR TSSetMaxTime(ts, ftime);
728  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
729 
730  auto T = smartCreateDMVector(simple->getDM());
731  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
732  SCATTER_FORWARD);
733  CHKERR TSSetSolution(ts, T);
734  CHKERR TSSetFromOptions(ts);
735  CHKERR set_ts(ts);
736  CHKERR set_section_monitor(ts);
737  CHKERR TSSetUp(ts);
738  CHKERR TSSolve(ts, NULL);
739  CHKERR TSGetTime(ts, &ftime);
740 
742 }
743 
744 int main(int argc, char *argv[]) {
745 
746  // Initialisation of MoFEM/PETSc and MOAB data structures
747  const char param_file[] = "param_file.petsc";
748  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
749 
750  // Add logging channel for example
751  auto core_log = logging::core::get();
752  core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FS"));
753  LogManager::setLog("FS");
754  MOFEM_LOG_TAG("FS", "free surface");
755 
756  try {
757 
758  //! [Register MoFEM discrete manager in PETSc]
759  DMType dm_name = "DMMOFEM";
760  CHKERR DMRegister_MoFEM(dm_name);
761  //! [Register MoFEM discrete manager in PETSc
762 
763  //! [Create MoAB]
764  moab::Core mb_instance; ///< mesh database
765  moab::Interface &moab = mb_instance; ///< mesh database interface
766  //! [Create MoAB]
767 
768  //! [Create MoFEM]
769  MoFEM::Core core(moab); ///< finite element database
770  MoFEM::Interface &m_field = core; ///< finite element database insterface
771  //! [Create MoFEM]
772 
773  //! [FreeSurface]
774  FreeSurface ex(m_field);
775  CHKERR ex.runProblem();
776  //! [FreeSurface]
777  }
778  CATCH_ERRORS;
779 
781 }
static Index< 'n', 3 > n
static Index< 'p', 3 > p
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
CoordinateTypes
Coodinate system.
Definition: definitions.h:127
@ CYLINDRICAL
Definition: definitions.h:130
@ CARTESIAN
Definition: definitions.h:128
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto cylindrical
int main(int argc, char *argv[])
auto get_M2
EntitiesFieldData::EntData EntData
constexpr int U_FIELD_DIM
constexpr double mu_m
constexpr double rho_ave
FTensor::Index< 'j', SPACE_DIM > j
auto get_M2_dh
constexpr double rho_m
auto integration_rule
static char help[]
constexpr CoordinateTypes coord_type
select coordinate system <CARTESIAN, CYLINDRICAL>;
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, U_FIELD_DIM > OpDomainSourceU
constexpr double mu_ave
auto get_M_dh
auto kernel_eye
constexpr double mu_p
constexpr double eps
FTensor::Index< 'k', SPACE_DIM > k
constexpr double lambda
auto get_M3_dh
FTensor::Index< 'i', SPACE_DIM > i
constexpr double mu_diff
auto init_h
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
constexpr auto t_kd
auto d_phase_function_h
constexpr int SPACE_DIM
constexpr double eta
FTensor::Index< 'l', SPACE_DIM > l
auto kernel_oscillation
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, 1 > OpDomainMassH
auto get_M0
constexpr double a0
constexpr double tol
auto get_f_dh
auto get_M3
OpDomainMassH OpDomainMassP
constexpr int BASE_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1, 1 > OpBaseTimesScalarField
constexpr double W
auto get_M
auto phase_function
constexpr double h
constexpr double rho_p
constexpr double rho_diff
auto get_D
auto get_M0_dh
auto my_max
auto cut_off
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, 1 > OpDomainSourceH
constexpr double eta2
const double kappa
constexpr double md
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
auto get_f
auto d_cut_off
auto my_min
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMMoFEM.cpp:213
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMMoFEM.cpp:414
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:234
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:461
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:257
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
const double T
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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: DMMMoFEM.cpp:994
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
double w(const double g, const double t)
const double D
diffusivity
const double r
rate factor
double A
int save_every_nth_step
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
PipelineManager::FaceEle DomainEle
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
FreeSurface(MoFEM::Interface &m_field)
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEM::Interface & mField
boost::shared_ptr< FEMethod > domianLhsFEPtr
MoFEMErrorCode solveSystem()
[Solve]
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:33
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
Interface for managing meshsets containing materials and boundary conditions.
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field valuse for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< VectorDouble > liftVec
boost::shared_ptr< BoundaryEle > liftFE
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc, std::pair< boost::shared_ptr< BoundaryEle >, boost::shared_ptr< VectorDouble >> p)
Postprocess on face.
Post processing.