v0.14.0
Public Member Functions | Public Attributes | Private Member Functions | Friends | List of all members
FreeSurface Struct Reference
Collaboration diagram for FreeSurface:
[legend]

Public Member Functions

 FreeSurface (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 [Run programme] More...
 
MoFEMErrorCode makeRefProblem ()
 

Public Attributes

MoFEM::Interface & mField
 

Private Member Functions

MoFEMErrorCode readMesh ()
 [Run programme] More...
 
MoFEMErrorCode setupProblem ()
 [Read mesh] More...
 
MoFEMErrorCode boundaryCondition ()
 [Set up problem] More...
 
MoFEMErrorCode projectData ()
 [Boundary condition] More...
 
MoFEMErrorCode assembleSystem ()
 [Data projection] More...
 
MoFEMErrorCode solveSystem ()
 [Solve] More...
 
std::vector< Range > findEntitiesCrossedByPhaseInterface (size_t overlap)
 Find entities on refinement levels. More...
 
Range findParentsToRefine (Range ents, BitRefLevel level, BitRefLevel mask)
 
MoFEMErrorCode refineMesh (size_t overlap)
 
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. More...
 

Friends

struct TSPrePostProc
 

Detailed Description

Examples
free_surface.cpp.

Definition at line 471 of file free_surface.cpp.

Constructor & Destructor Documentation

◆ FreeSurface()

FreeSurface::FreeSurface ( MoFEM::Interface &  m_field)
inline

Definition at line 473 of file free_surface.cpp.

473 : mField(m_field) {}

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode FreeSurface::assembleSystem ( )
private

[Data projection]

[Push operators to pip]

Examples
free_surface.cpp.

Definition at line 1581 of file free_surface.cpp.

1581  {
1582  MoFEMFunctionBegin;
1583  auto simple = mField.getInterface<Simple>();
1584 
1585  using UDO = ForcesAndSourcesCore::UserDataOperator;
1586 
1587  auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
1588  return setParentDofs(
1589  fe, field, op, bit(get_skin_parent_bit()),
1590 
1591  [&]() {
1592  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1593  new DomianParentEle(mField));
1594  return fe_parent;
1595  },
1596 
1597  QUIET, Sev::noisy);
1598  };
1599 
1600  auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
1601  return setParentDofs(
1602  fe, field, op, bit(get_skin_parent_bit()),
1603 
1604  [&]() {
1605  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1606  new BoundaryParentEle(mField));
1607  return fe_parent;
1608  },
1609 
1610  QUIET, Sev::noisy);
1611  };
1612 
1613  auto test_bit_child = [](FEMethod *fe_ptr) {
1614  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
1615  get_start_bit() + nb_levels - 1);
1616  };
1617 
1618  auto dot_u_ptr = boost::make_shared<MatrixDouble>();
1619  auto u_ptr = boost::make_shared<MatrixDouble>();
1620  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
1621  auto dot_h_ptr = boost::make_shared<VectorDouble>();
1622  auto h_ptr = boost::make_shared<VectorDouble>();
1623  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
1624  auto g_ptr = boost::make_shared<VectorDouble>();
1625  auto grad_g_ptr = boost::make_shared<MatrixDouble>();
1626  auto lambda_ptr = boost::make_shared<VectorDouble>();
1627  auto p_ptr = boost::make_shared<VectorDouble>();
1628  auto div_u_ptr = boost::make_shared<VectorDouble>();
1629 
1630  // Push element from reference configuration to current configuration in 3d
1631  // space
1632  auto set_domain_general = [&](auto fe) {
1633  MoFEMFunctionBegin;
1634  auto &pip = fe->getOpPtrVector();
1635 
1636  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
1637 
1638  CHKERR setParentDofs(
1639  fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1640 
1641  [&]() {
1642  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1643  new DomianParentEle(mField));
1644  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1645  fe_parent->getOpPtrVector(), {H1});
1646  return fe_parent;
1647  },
1648 
1649  QUIET, Sev::noisy);
1650 
1651  CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1652  pip.push_back(
1653  new OpCalculateVectorFieldValuesDot<U_FIELD_DIM>("U", dot_u_ptr));
1654  pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
1655  pip.push_back(new OpCalculateVectorFieldGradient<U_FIELD_DIM, SPACE_DIM>(
1656  "U", grad_u_ptr));
1657 
1658  switch (coord_type) {
1659  case CARTESIAN:
1660  pip.push_back(
1661  new OpCalculateDivergenceVectorFieldValues<SPACE_DIM, CARTESIAN>(
1662  "U", div_u_ptr));
1663  break;
1664  case CYLINDRICAL:
1665  pip.push_back(
1666  new OpCalculateDivergenceVectorFieldValues<SPACE_DIM, CYLINDRICAL>(
1667  "U", div_u_ptr));
1668  break;
1669  default:
1670  SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
1671  "Coordinate system not implemented");
1672  }
1673 
1674  CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1675  pip.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
1676  pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1677  pip.push_back(
1678  new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1679 
1680  CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1681  pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1682  pip.push_back(
1683  new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
1684 
1685  CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1686  pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1687  MoFEMFunctionReturn(0);
1688  };
1689 
1690  auto set_domain_rhs = [&](auto fe) {
1691  MoFEMFunctionBegin;
1692  auto &pip = fe->getOpPtrVector();
1693 
1694  CHKERR set_domain_general(fe);
1695 
1696  CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1697  pip.push_back(new OpRhsU("U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
1698  grad_h_ptr, g_ptr, p_ptr));
1699 
1700  CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1701  pip.push_back(new OpRhsH<false>("H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
1702  grad_g_ptr));
1703 
1704  CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1705  pip.push_back(new OpRhsG<false>("G", h_ptr, grad_h_ptr, g_ptr));
1706 
1707  CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1708  pip.push_back(new OpDomainAssembleScalar(
1709  "P", div_u_ptr, [](const double r, const double, const double) {
1710  return cylindrical(r);
1711  }));
1712  pip.push_back(new OpDomainAssembleScalar(
1713  "P", p_ptr, [](const double r, const double, const double) {
1714  return eps * cylindrical(r);
1715  }));
1716  MoFEMFunctionReturn(0);
1717  };
1718 
1719  auto set_domain_lhs = [&](auto fe) {
1720  MoFEMFunctionBegin;
1721  auto &pip = fe->getOpPtrVector();
1722 
1723  CHKERR set_domain_general(fe);
1724 
1725  CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
1726  {
1727  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1728  pip.push_back(new OpLhsU_dU("U", u_ptr, grad_u_ptr, h_ptr));
1729  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1730  pip.push_back(
1731  new OpLhsU_dH("U", "H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
1732 
1733  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1734  pip.push_back(new OpLhsU_dG("U", "G", grad_h_ptr));
1735  }
1736 
1737  CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
1738  {
1739  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1740  pip.push_back(new OpLhsH_dU("H", "U", grad_h_ptr));
1741  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1742  pip.push_back(new OpLhsH_dH<false>("H", u_ptr, h_ptr, grad_g_ptr));
1743  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1744  pip.push_back(new OpLhsH_dG<false>("H", "G", h_ptr));
1745  }
1746 
1747  CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
1748  {
1749  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
1750  pip.push_back(new OpLhsG_dH<false>("G", "H", h_ptr));
1751  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
1752  pip.push_back(new OpLhsG_dG("G"));
1753  }
1754 
1755  CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
1756  {
1757  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
1758 
1759  switch (coord_type) {
1760  case CARTESIAN:
1761  pip.push_back(new OpMixScalarTimesDiv<CARTESIAN>(
1762  "P", "U",
1763  [](const double r, const double, const double) {
1764  return cylindrical(r);
1765  },
1766  true, false));
1767  break;
1768  case CYLINDRICAL:
1769  pip.push_back(new OpMixScalarTimesDiv<CYLINDRICAL>(
1770  "P", "U",
1771  [](const double r, const double, const double) {
1772  return cylindrical(r);
1773  },
1774  true, false));
1775  break;
1776  default:
1777  SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
1778  "Coordinate system not implemented");
1779  }
1780 
1781  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P");
1782  pip.push_back(new OpDomainMassP("P", "P", [](double r, double, double) {
1783  return eps * cylindrical(r);
1784  }));
1785  }
1786 
1787  MoFEMFunctionReturn(0);
1788  };
1789 
1790  auto get_block_name = [](auto name) {
1791  return boost::format("%s(.*)") % "WETTING_ANGLE";
1792  };
1793 
1794  auto get_blocks = [&](auto &&name) {
1795  return mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
1796  std::regex(name.str()));
1797  };
1798 
1799  auto set_boundary_rhs = [&](auto fe) {
1800  MoFEMFunctionBegin;
1801  auto &pip = fe->getOpPtrVector();
1802 
1803  CHKERR setParentDofs(
1804  fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1805 
1806  [&]() {
1807  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1808  new BoundaryParentEle(mField));
1809  return fe_parent;
1810  },
1811 
1812  QUIET, Sev::noisy);
1813 
1814  CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
1815  pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
1816 
1817  CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
1818  pip.push_back(new OpCalculateScalarFieldValues("L", lambda_ptr));
1819  pip.push_back(new OpNormalConstrainRhs("L", u_ptr));
1820 
1821  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpFluidFlux>::add(
1822  pip, mField, "L", {}, "INFLUX",
1823  [](double r, double, double) { return cylindrical(r); }, Sev::inform);
1824 
1825  CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
1826  pip.push_back(new OpNormalForceRhs("U", lambda_ptr));
1827 
1828  auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
1829  if (wetting_block.size()) {
1830  // push operators to the side element which is called from op_bdy_side
1831  auto op_bdy_side =
1832  new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
1833  op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
1834 
1835  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1836  op_bdy_side->getOpPtrVector(), {H1});
1837 
1838  CHKERR setParentDofs(
1839  op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
1840  bit(get_skin_parent_bit()),
1841 
1842  [&]() {
1843  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1844  new DomianParentEle(mField));
1845  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1846  fe_parent->getOpPtrVector(), {H1});
1847  return fe_parent;
1848  },
1849 
1850  QUIET, Sev::noisy);
1851 
1852  CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1853  "H");
1854  op_bdy_side->getOpPtrVector().push_back(
1855  new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1856  // push bdy side op
1857  pip.push_back(op_bdy_side);
1858 
1859  // push operators for rhs wetting angle
1860  for (auto &b : wetting_block) {
1861  Range force_edges;
1862  std::vector<double> attr_vec;
1863  CHKERR b->getMeshsetIdEntitiesByDimension(
1864  mField.get_moab(), SPACE_DIM - 1, force_edges, true);
1865  b->getAttributes(attr_vec);
1866  if (attr_vec.size() != 1)
1867  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1868  "Should be one attribute");
1869  MOFEM_LOG("FS", Sev::inform) << "Wetting angle: " << attr_vec.front();
1870  // need to find the attributes and pass to operator
1871  CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
1872  pip.push_back(new OpWettingAngleRhs(
1873  "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
1874  attr_vec.front()));
1875  }
1876  }
1877 
1878  MoFEMFunctionReturn(0);
1879  };
1880 
1881  auto set_boundary_lhs = [&](auto fe) {
1882  MoFEMFunctionBegin;
1883  auto &pip = fe->getOpPtrVector();
1884 
1885  CHKERR setParentDofs(
1886  fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1887 
1888  [&]() {
1889  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1890  new BoundaryParentEle(mField));
1891  return fe_parent;
1892  },
1893 
1894  QUIET, Sev::noisy);
1895 
1896  CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
1897  CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "U");
1898  pip.push_back(new OpNormalConstrainLhs("L", "U"));
1899 
1900  auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
1901  if (wetting_block.size()) {
1902  auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
1903  auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
1904 
1905  // push operators to the side element which is called from op_bdy_side
1906  auto op_bdy_side =
1907  new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
1908  op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
1909 
1910  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1911  op_bdy_side->getOpPtrVector(), {H1});
1912 
1913  CHKERR setParentDofs(
1914  op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
1915  bit(get_skin_parent_bit()),
1916 
1917  [&]() {
1918  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1919  new DomianParentEle(mField));
1920  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1921  fe_parent->getOpPtrVector(), {H1});
1922  return fe_parent;
1923  },
1924 
1925  QUIET, Sev::noisy);
1926 
1927  CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1928  "H");
1929  op_bdy_side->getOpPtrVector().push_back(
1930  new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1931  CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
1932  "H");
1933  op_bdy_side->getOpPtrVector().push_back(
1934  new OpLoopSideGetDataForSideEle("H", col_ind_ptr, col_diff_base_ptr));
1935 
1936  // push bdy side op
1937  pip.push_back(op_bdy_side);
1938 
1939  // push operators for lhs wetting angle
1940  for (auto &b : wetting_block) {
1941  Range force_edges;
1942  std::vector<double> attr_vec;
1943  CHKERR b->getMeshsetIdEntitiesByDimension(
1944  mField.get_moab(), SPACE_DIM - 1, force_edges, true);
1945  b->getAttributes(attr_vec);
1946  if (attr_vec.size() != 1)
1947  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1948  "Should be one attribute");
1949  MOFEM_LOG("FS", Sev::inform)
1950  << "wetting angle edges size " << force_edges.size();
1951 
1952  CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
1953  pip.push_back(new OpWettingAngleLhs(
1954  "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
1955  boost::make_shared<Range>(force_edges), attr_vec.front()));
1956  }
1957  }
1958 
1959  MoFEMFunctionReturn(0);
1960  };
1961 
1962  auto *pip_mng = mField.getInterface<PipelineManager>();
1963 
1964  CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1965  CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1966  CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
1967  CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
1968 
1969  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
1970  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
1971  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
1972  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
1973 
1974  pip_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
1975  pip_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
1976  pip_mng->getBoundaryLhsFE()->exeTestHook = test_bit_child;
1977  pip_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
1978 
1979  MoFEMFunctionReturn(0);
1980 }

◆ boundaryCondition()

MoFEMErrorCode FreeSurface::boundaryCondition ( )
private

[Set up problem]

[Boundary condition]

Examples
free_surface.cpp.

Definition at line 663 of file free_surface.cpp.

663  {
664  MoFEMFunctionBegin;
665 
666 #ifdef PYTHON_INIT_SURFACE
667  auto get_py_surface_init = []() {
668  auto py_surf_init = boost::make_shared<SurfacePython>();
669  CHKERR py_surf_init->surfaceInit("surface.py");
670  surfacePythonWeakPtr = py_surf_init;
671  return py_surf_init;
672  };
673  auto py_surf_init = get_py_surface_init();
674 #endif
675 
676  auto simple = mField.getInterface<Simple>();
677  auto pip_mng = mField.getInterface<PipelineManager>();
678  auto bc_mng = mField.getInterface<BcManager>();
679  auto bit_mng = mField.getInterface<BitRefManager>();
680  auto dm = simple->getDM();
681 
682  using UDO = ForcesAndSourcesCore::UserDataOperator;
683 
684  auto reset_bits = [&]() {
685  MoFEMFunctionBegin;
686  BitRefLevel start_mask;
687  for (auto s = 0; s != get_start_bit(); ++s)
688  start_mask[s] = true;
689  // reset bit ref levels
690  CHKERR bit_mng->lambdaBitRefLevel(
691  [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
692  Range level0;
693  CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
694  CHKERR bit_mng->setNthBitRefLevel(level0, get_current_bit(), true);
695  CHKERR bit_mng->setNthBitRefLevel(level0, get_projection_bit(), true);
696  MoFEMFunctionReturn(0);
697  };
698 
699  auto add_parent_field = [&](auto fe, auto op, auto field) {
700  return setParentDofs(
701  fe, field, op, bit(get_skin_parent_bit()),
702 
703  [&]() {
704  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
705  new DomianParentEle(mField));
706  return fe_parent;
707  },
708 
709  QUIET, Sev::noisy);
710  };
711 
712  auto h_ptr = boost::make_shared<VectorDouble>();
713  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
714  auto g_ptr = boost::make_shared<VectorDouble>();
715  auto grad_g_ptr = boost::make_shared<MatrixDouble>();
716 
717  auto set_generic = [&](auto fe) {
718  MoFEMFunctionBegin;
719  auto &pip = fe->getOpPtrVector();
720 
721  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
722 
723  CHKERR setParentDofs(
724  fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
725 
726  [&]() {
727  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
728  new DomianParentEle(mField));
729  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
730  fe_parent->getOpPtrVector(), {H1});
731  return fe_parent;
732  },
733 
734  QUIET, Sev::noisy);
735 
736  CHKERR add_parent_field(fe, UDO::OPROW, "H");
737  pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
738  pip.push_back(
739  new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
740 
741  CHKERR add_parent_field(fe, UDO::OPROW, "G");
742  pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
743  pip.push_back(
744  new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
745 
746  MoFEMFunctionReturn(0);
747  };
748 
749  auto post_proc = [&](auto exe_test) {
750  MoFEMFunctionBegin;
751  auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
752  post_proc_fe->exeTestHook = exe_test;
753 
754  CHKERR set_generic(post_proc_fe);
755 
756  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
757 
758  post_proc_fe->getOpPtrVector().push_back(
759 
760  new OpPPMap(post_proc_fe->getPostProcMesh(),
761  post_proc_fe->getMapGaussPts(),
762 
763  {{"H", h_ptr}, {"G", g_ptr}},
764 
765  {{"GRAD_H", grad_h_ptr}, {"GRAD_G", grad_g_ptr}},
766 
767  {},
768 
769  {}
770 
771  )
772 
773  );
774 
775  CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
776  CHKERR post_proc_fe->writeFile("out_init.h5m");
777 
778  MoFEMFunctionReturn(0);
779  };
780 
781  auto solve_init = [&](auto exe_test) {
782  MoFEMFunctionBegin;
783 
784  pip_mng->getOpDomainRhsPipeline().clear();
785  pip_mng->getOpDomainLhsPipeline().clear();
786 
787  auto set_domain_rhs = [&](auto fe) {
788  MoFEMFunctionBegin;
789  CHKERR set_generic(fe);
790  auto &pip = fe->getOpPtrVector();
791 
792  CHKERR add_parent_field(fe, UDO::OPROW, "H");
793  pip.push_back(new OpRhsH<true>("H", nullptr, nullptr, h_ptr, grad_h_ptr,
794  grad_g_ptr));
795  CHKERR add_parent_field(fe, UDO::OPROW, "G");
796  pip.push_back(new OpRhsG<true>("G", h_ptr, grad_h_ptr, g_ptr));
797  MoFEMFunctionReturn(0);
798  };
799 
800  auto set_domain_lhs = [&](auto fe) {
801  MoFEMFunctionBegin;
802 
803  CHKERR set_generic(fe);
804  auto &pip = fe->getOpPtrVector();
805 
806  CHKERR add_parent_field(fe, UDO::OPROW, "H");
807  CHKERR add_parent_field(fe, UDO::OPCOL, "H");
808  pip.push_back(new OpLhsH_dH<true>("H", nullptr, h_ptr, grad_g_ptr));
809 
810  CHKERR add_parent_field(fe, UDO::OPCOL, "G");
811  pip.push_back(new OpLhsH_dG<true>("H", "G", h_ptr));
812 
813  CHKERR add_parent_field(fe, UDO::OPROW, "G");
814  pip.push_back(new OpLhsG_dG("G"));
815 
816  CHKERR add_parent_field(fe, UDO::OPCOL, "H");
817  pip.push_back(new OpLhsG_dH<true>("G", "H", h_ptr));
818  MoFEMFunctionReturn(0);
819  };
820 
821  auto create_subdm = [&]() {
822  auto level_ents_ptr = boost::make_shared<Range>();
823  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
824  bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
825 
826  DM subdm;
827  CHKERR DMCreate(mField.get_comm(), &subdm);
828  CHKERR DMSetType(subdm, "DMMOFEM");
829  CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_INIT");
830  CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
831  CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
832  CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_TRUE);
833 
834  for (auto f : {"H", "G"}) {
835  CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
836  CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
837  }
838  CHKERR DMSetUp(subdm);
839 
840  if constexpr (debug) {
841  if (mField.get_comm_size() == 1) {
842  auto dm_ents = get_dofs_ents_all(SmartPetscObj<DM>(subdm, true));
843  CHKERR save_range(mField.get_moab(), "sub_dm_init_ents_verts.h5m",
844  dm_ents.subset_by_type(MBVERTEX));
845  dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
846  CHKERR save_range(mField.get_moab(), "sub_dm_init_ents.h5m", dm_ents);
847  }
848  }
849 
850  return SmartPetscObj<DM>(subdm);
851  };
852 
853  auto subdm = create_subdm();
854 
855  pip_mng->getDomainRhsFE().reset();
856  pip_mng->getDomainLhsFE().reset();
857  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
858  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
859  pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
860  pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
861 
862  CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
863  CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
864 
865  auto D = createDMVector(subdm);
866  auto snes = pip_mng->createSNES(subdm);
867  auto snes_ctx_ptr = getDMSnesCtx(subdm);
868 
869  auto set_section_monitor = [&](auto solver) {
870  MoFEMFunctionBegin;
871  CHKERR SNESMonitorSet(solver,
872  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
873  void *))MoFEMSNESMonitorFields,
874  (void *)(snes_ctx_ptr.get()), nullptr);
875 
876  MoFEMFunctionReturn(0);
877  };
878 
879  auto print_section_field = [&]() {
880  MoFEMFunctionBegin;
881 
882  auto section =
883  mField.getInterface<ISManager>()->sectionCreate("SUB_INIT");
884  PetscInt num_fields;
885  CHKERR PetscSectionGetNumFields(section, &num_fields);
886  for (int f = 0; f < num_fields; ++f) {
887  const char *field_name;
888  CHKERR PetscSectionGetFieldName(section, f, &field_name);
889  MOFEM_LOG("FS", Sev::inform)
890  << "Field " << f << " " << std::string(field_name);
891  }
892 
893  MoFEMFunctionReturn(0);
894  };
895 
896  CHKERR set_section_monitor(snes);
897  CHKERR print_section_field();
898 
899  for (auto f : {"U", "P", "H", "G", "L"}) {
900  CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
901  }
902 
903  CHKERR SNESSetFromOptions(snes);
904  CHKERR SNESSolve(snes, PETSC_NULL, D);
905 
906  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
907  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
908  CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
909 
910  MoFEMFunctionReturn(0);
911  };
912 
913  CHKERR reset_bits();
914  CHKERR solve_init(
915  [](FEMethod *fe_ptr) { return get_fe_bit(fe_ptr).test(0); });
916  CHKERR refineMesh(refine_overlap);
917 
918  for (auto f : {"U", "P", "H", "G", "L"}) {
919  CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
920  }
921  CHKERR solve_init([](FEMethod *fe_ptr) {
922  return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
923  });
924 
925  CHKERR post_proc([](FEMethod *fe_ptr) {
926  return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
927  });
928 
929  pip_mng->getOpDomainRhsPipeline().clear();
930  pip_mng->getOpDomainLhsPipeline().clear();
931 
932  // Remove DOFs where boundary conditions are set
933  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
934  "U", 0, 0);
935  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
936  "L", 0, 0);
937  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
938  "U", 1, 1);
939  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
940  "L", 1, 1);
941  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "U",
942  0, SPACE_DIM);
943  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "L",
944  0, 0);
945  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "ZERO",
946  "L", 0, 0);
947 
948  MoFEMFunctionReturn(0);
949 }

◆ findEntitiesCrossedByPhaseInterface()

std::vector< Range > FreeSurface::findEntitiesCrossedByPhaseInterface ( size_t  overlap)
private

Find entities on refinement levels.

Parameters
overlaplevel of overlap
Returns
Examples
free_surface.cpp.

Definition at line 2453 of file free_surface.cpp.

2453  {
2454 
2455  auto &moab = mField.get_moab();
2456  auto bit_mng = mField.getInterface<BitRefManager>();
2457  auto comm_mng = mField.getInterface<CommInterface>();
2458 
2459  Range vertices;
2460  CHK_THROW_MESSAGE(bit_mng->getEntitiesByTypeAndRefLevel(
2461  bit(0), BitRefLevel().set(), MBVERTEX, vertices),
2462  "can not get vertices on bit 0");
2463 
2464  auto &dofs_mi = mField.get_dofs()->get<Unique_mi_tag>();
2465  auto field_bit_number = mField.get_field_bit_number("H");
2466 
2467  Range plus_range, minus_range;
2468  std::vector<EntityHandle> plus, minus;
2469 
2470  // get vertices on level 0 on plus and minus side
2471  for (auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
2472 
2473  const auto f = p->first;
2474  const auto s = p->second;
2475 
2476  // Lowest Dof UId for given field (field bit number) on entity f
2477  const auto lo_uid = DofEntity::getLoFieldEntityUId(field_bit_number, f);
2478  const auto hi_uid = DofEntity::getHiFieldEntityUId(field_bit_number, s);
2479  auto it = dofs_mi.lower_bound(lo_uid);
2480  const auto hi_it = dofs_mi.upper_bound(hi_uid);
2481 
2482  plus.clear();
2483  minus.clear();
2484  plus.reserve(std::distance(it, hi_it));
2485  minus.reserve(std::distance(it, hi_it));
2486 
2487  for (; it != hi_it; ++it) {
2488  const auto v = (*it)->getFieldData();
2489  if (v > 0)
2490  plus.push_back((*it)->getEnt());
2491  else
2492  minus.push_back((*it)->getEnt());
2493  }
2494 
2495  plus_range.insert_list(plus.begin(), plus.end());
2496  minus_range.insert_list(minus.begin(), minus.end());
2497  }
2498 
2499  MOFEM_LOG_CHANNEL("SYNC");
2500  MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
2501  << "Plus range " << plus_range << endl;
2502  MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
2503  << "Minus range " << minus_range << endl;
2504  MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), Sev::noisy);
2505 
2506  auto get_elems = [&](auto &ents, auto bit, auto mask) {
2507  Range adj;
2508  CHK_MOAB_THROW(moab.get_adjacencies(ents, SPACE_DIM, false, adj,
2509  moab::Interface::UNION),
2510  "can not get adjacencies");
2511  CHK_THROW_MESSAGE(bit_mng->filterEntitiesByRefLevel(bit, mask, adj),
2512  "can not filter elements with bit 0");
2513  return adj;
2514  };
2515 
2516  CHKERR comm_mng->synchroniseEntities(plus_range);
2517  CHKERR comm_mng->synchroniseEntities(minus_range);
2518 
2519  std::vector<Range> ele_plus(nb_levels), ele_minus(nb_levels);
2520  ele_plus[0] = get_elems(plus_range, bit(0), BitRefLevel().set());
2521  ele_minus[0] = get_elems(minus_range, bit(0), BitRefLevel().set());
2522  auto common = intersect(ele_plus[0], ele_minus[0]);
2523  ele_plus[0] = subtract(ele_plus[0], common);
2524  ele_minus[0] = subtract(ele_minus[0], common);
2525 
2526  auto get_children = [&](auto &p, auto &c) {
2527  MoFEMFunctionBegin;
2528  CHKERR bit_mng->updateRangeByChildren(p, c);
2529  c = c.subset_by_dimension(SPACE_DIM);
2530  MoFEMFunctionReturn(0);
2531  };
2532 
2533  for (auto l = 1; l != nb_levels; ++l) {
2534  CHK_THROW_MESSAGE(get_children(ele_plus[l - 1], ele_plus[l]),
2535  "get children");
2536  CHK_THROW_MESSAGE(get_children(ele_minus[l - 1], ele_minus[l]),
2537  "get children");
2538  }
2539 
2540  auto get_level = [&](auto &p, auto &m, auto z, auto bit, auto mask) {
2541  Range l;
2542  CHK_THROW_MESSAGE(
2543  bit_mng->getEntitiesByDimAndRefLevel(bit, mask, SPACE_DIM, l),
2544  "can not get vertices on bit");
2545  l = subtract(l, p);
2546  l = subtract(l, m);
2547  for (auto f = 0; f != z; ++f) {
2548  Range conn;
2549  CHK_MOAB_THROW(moab.get_connectivity(l, conn, true), "");
2550  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(conn);
2551  l = get_elems(conn, bit, mask);
2552  }
2553  return l;
2554  };
2555 
2556  std::vector<Range> vec_levels(nb_levels);
2557  for (auto l = nb_levels - 1; l >= 0; --l) {
2558  vec_levels[l] = get_level(ele_plus[l], ele_minus[l], 2 * overlap, bit(l),
2559  BitRefLevel().set());
2560  }
2561 
2562  if constexpr (debug) {
2563  if (mField.get_comm_size() == 1) {
2564  for (auto l = 0; l != nb_levels; ++l) {
2565  std::string name = (boost::format("out_r%d.h5m") % l).str();
2566  CHK_THROW_MESSAGE(save_range(mField.get_moab(), name, vec_levels[l]),
2567  "save mesh");
2568  }
2569  }
2570  }
2571 
2572  return vec_levels;
2573 }

◆ findParentsToRefine()

Range FreeSurface::findParentsToRefine ( Range  ents,
BitRefLevel  level,
BitRefLevel  mask 
)
private
Parameters
ents
level
mask
Returns

◆ makeRefProblem()

MoFEMErrorCode FreeSurface::makeRefProblem ( )

◆ projectData()

MoFEMErrorCode FreeSurface::projectData ( )
private

[Boundary condition]

[Data projection]

Zero DOFs, used by FieldBlas

get TSTheta data operators

Examples
free_surface.cpp.

Definition at line 953 of file free_surface.cpp.

953  {
954  MoFEMFunctionBegin;
955 
956  // FIXME: Functionality in this method should be improved (projection should
957  // be done field by field), generalised, and move to become core
958  // functionality.
959 
960  auto simple = mField.getInterface<Simple>();
961  auto pip_mng = mField.getInterface<PipelineManager>();
962  auto bit_mng = mField.getInterface<BitRefManager>();
963  auto field_blas = mField.getInterface<FieldBlas>();
964 
965  // Store all existing elements pipelines, replace them by data projection
966  // pipelines, to put them back when projection is done.
967  auto fe_domain_lhs = pip_mng->getDomainLhsFE();
968  auto fe_domain_rhs = pip_mng->getDomainRhsFE();
969  auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
970  auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
971 
972  pip_mng->getDomainLhsFE().reset();
973  pip_mng->getDomainRhsFE().reset();
974  pip_mng->getBoundaryLhsFE().reset();
975  pip_mng->getBoundaryRhsFE().reset();
976 
977  using UDO = ForcesAndSourcesCore::UserDataOperator;
978 
979  // extract field data for domain parent element
980  auto add_parent_field_domain = [&](auto fe, auto op, auto field, auto bit) {
981  return setParentDofs(
982  fe, field, op, bit,
983 
984  [&]() {
985  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
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) {
1033  MoFEMFunctionBegin;
1034 
1035  auto set_domain_rhs = [&](auto fe) {
1036  MoFEMFunctionBegin;
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();
1047  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
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",
1061  bit(get_skin_projection_bit()));
1062  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1063  CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "P",
1064  bit(get_skin_projection_bit()));
1065  pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1066  CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "H",
1067  bit(get_skin_projection_bit()));
1068  pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1069  CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "G",
1070  bit(get_skin_projection_bit()));
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();
1081  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
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",
1093  bit(get_skin_parent_bit()));
1094  pip.push_back(new OpDomainAssembleVector("U", u_ptr));
1095  CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "P",
1096  bit(get_skin_parent_bit()));
1097  pip.push_back(new OpDomainAssembleScalar("P", p_ptr));
1098  CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "H",
1099  bit(get_skin_parent_bit()));
1100  pip.push_back(new OpDomainAssembleScalar("H", h_ptr));
1101  CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "G",
1102  bit(get_skin_parent_bit()));
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 
1110  MoFEMFunctionReturn(0);
1111  };
1112 
1113  auto set_domain_lhs = [&](auto fe) {
1114  MoFEMFunctionBegin;
1115 
1116  auto &pip = fe->getOpPtrVector();
1117 
1118  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
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",
1134  bit(get_skin_parent_bit()));
1135  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U",
1136  bit(get_skin_parent_bit()));
1137  pip.push_back(new OpDomainMassU("U", "U"));
1138  CHKERR add_parent_field_domain(fe, UDO::OPROW, "P",
1139  bit(get_skin_parent_bit()));
1140  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P",
1141  bit(get_skin_parent_bit()));
1142  pip.push_back(new OpDomainMassP("P", "P"));
1143  CHKERR add_parent_field_domain(fe, UDO::OPROW, "H",
1144  bit(get_skin_parent_bit()));
1145  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H",
1146  bit(get_skin_parent_bit()));
1147  pip.push_back(new OpDomainMassH("H", "H"));
1148  CHKERR add_parent_field_domain(fe, UDO::OPROW, "G",
1149  bit(get_skin_parent_bit()));
1150  CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G",
1151  bit(get_skin_parent_bit()));
1152  pip.push_back(new OpDomainMassG("G", "G"));
1153 
1154  MoFEMFunctionReturn(0);
1155  };
1156 
1157  auto set_bdy_rhs = [&](auto fe) {
1158  MoFEMFunctionBegin;
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",
1179  bit(get_skin_projection_bit()));
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 &) {
1205  MoFEMFunctionBegin;
1206  if (lPtr->size() != getGaussPts().size2()) {
1207  lPtr->resize(getGaussPts().size2());
1208  lPtr->clear();
1209  }
1210  MoFEMFunctionReturn(0);
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",
1220  bit(get_skin_parent_bit()));
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 
1228  MoFEMFunctionReturn(0);
1229  };
1230 
1231  auto set_bdy_lhs = [&](auto fe) {
1232  MoFEMFunctionBegin;
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",
1250  bit(get_skin_parent_bit()));
1251  CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "L",
1252  bit(get_skin_parent_bit()));
1253  pip.push_back(new OpBoundaryMassL("L", "L"));
1254 
1255  MoFEMFunctionReturn(0);
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>();
1281  MoFEMFunctionBegin;
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 
1304  MoFEMFunctionReturn(0);
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) {
1376  MoFEMFunctionBegin;
1377  for (auto &v : ent_ptr->getEntFieldData()) {
1378  v = 0;
1379  }
1380  MoFEMFunctionReturn(0);
1381  };
1382 
1383  auto solve = [&](auto S) {
1384  MoFEMFunctionBegin;
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 
1395  MoFEMFunctionReturn(0);
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 = [&]() {
1412  MoFEMFunctionBegin;
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);
1417  MoFEMFunctionReturn(0);
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 
1426  CHK_THROW_MESSAGE(
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");
1432  CHK_THROW_MESSAGE(
1433  DMGetNamedGlobalVector(dummy_dm, "TSTheta_Xdot", &Xdot),
1434  "get Xdot");
1435 
1436  auto forward_ghost = [](auto g) {
1437  MoFEMFunctionBegin;
1438  CHKERR VecGhostUpdateBegin(g, INSERT_VALUES, SCATTER_FORWARD);
1439  CHKERR VecGhostUpdateEnd(g, INSERT_VALUES, SCATTER_FORWARD);
1440  MoFEMFunctionReturn(0);
1441  };
1442 
1443  CHK_THROW_MESSAGE(forward_ghost(X0), "");
1444  CHK_THROW_MESSAGE(forward_ghost(Xdot), "");
1445 
1446  if constexpr (debug) {
1447  MOFEM_LOG("FS", Sev::inform)
1448  << "Reverse restrict: X0 " << get_norm(X0) << " Xdot "
1449  << get_norm(Xdot);
1450  }
1451 
1452  return std::vector<Vec>{X0, Xdot};
1453  };
1454 
1455  auto ts_solver_vecs = apply_restrict();
1456 
1457  if (ts_solver_vecs.size()) {
1458 
1459  for (auto v : ts_solver_vecs) {
1460  MOFEM_LOG("FS", Sev::inform) << "Solve projection vector";
1461 
1462  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1463  SCATTER_REVERSE);
1464  CHKERR solve(sub_x);
1465 
1466  for (auto f : {"U", "P", "H", "G", "L"}) {
1467  MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1468  CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1469  }
1470 
1471  CHKERR DMoFEMMeshToLocalVector(subdm, sub_x, INSERT_VALUES,
1472  SCATTER_REVERSE);
1473  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1474  SCATTER_FORWARD);
1475 
1476  MOFEM_LOG("FS", Sev::inform) << "Norm V " << get_norm(v);
1477  }
1478 
1479  CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_X0",
1480  &ts_solver_vecs[0]);
1481  CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_Xdot",
1482  &ts_solver_vecs[1]);
1483  }
1484 
1485  for (auto f : {"U", "P", "H", "G", "L"}) {
1486  MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1487  CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1488  }
1489  CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
1490  }
1491 
1492  MoFEMFunctionReturn(0);
1493  };
1494 
1495  // postprocessing (only for debugging proposes)
1496  auto post_proc = [&](auto exe_test) {
1497  MoFEMFunctionBegin;
1498  auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
1499  auto &pip = post_proc_fe->getOpPtrVector();
1500  post_proc_fe->exeTestHook = exe_test;
1501 
1502  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
1503 
1504  CHKERR setParentDofs(
1505  post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1506 
1507  [&]() {
1508  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1509  new DomianParentEle(mField));
1510  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1511  fe_parent->getOpPtrVector(), {H1});
1512  return fe_parent;
1513  },
1514 
1515  QUIET, Sev::noisy);
1516 
1517  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
1518 
1519  auto u_ptr = boost::make_shared<MatrixDouble>();
1520  auto p_ptr = boost::make_shared<VectorDouble>();
1521  auto h_ptr = boost::make_shared<VectorDouble>();
1522  auto g_ptr = boost::make_shared<VectorDouble>();
1523 
1524  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U",
1525  bit(get_skin_parent_bit()));
1526  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1527  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P",
1528  bit(get_skin_parent_bit()));
1529  pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1530  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H",
1531  bit(get_skin_parent_bit()));
1532  pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1533  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G",
1534  bit(get_skin_parent_bit()));
1535  pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1536 
1537  post_proc_fe->getOpPtrVector().push_back(
1538 
1539  new OpPPMap(post_proc_fe->getPostProcMesh(),
1540  post_proc_fe->getMapGaussPts(),
1541 
1542  {{"P", p_ptr}, {"H", h_ptr}, {"G", g_ptr}},
1543 
1544  {{"U", u_ptr}},
1545 
1546  {},
1547 
1548  {}
1549 
1550  )
1551 
1552  );
1553 
1554  auto dm = simple->getDM();
1555  CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
1556  CHKERR post_proc_fe->writeFile("out_projection.h5m");
1557 
1558  MoFEMFunctionReturn(0);
1559  };
1560 
1561  CHKERR solve_projection([](FEMethod *fe_ptr) {
1562  return get_fe_bit(fe_ptr).test(get_current_bit());
1563  });
1564 
1565  if constexpr (debug) {
1566  CHKERR post_proc([](FEMethod *fe_ptr) {
1567  return get_fe_bit(fe_ptr).test(get_current_bit());
1568  });
1569  }
1570 
1571  fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
1572  fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
1573  fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
1574  fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
1575 
1576  MoFEMFunctionReturn(0);
1577 }

◆ readMesh()

MoFEMErrorCode FreeSurface::readMesh ( )
private

[Run programme]

[Read mesh]

Examples
free_surface.cpp.

Definition at line 536 of file free_surface.cpp.

536  {
537  MoFEMFunctionBegin;
538  MOFEM_LOG("FS", Sev::inform) << "Read mesh for problem";
539  auto simple = mField.getInterface<Simple>();
540 
541  simple->getParentAdjacencies() = true;
542  simple->getBitRefLevel() = BitRefLevel();
543 
544  CHKERR simple->getOptions();
545  CHKERR simple->loadFile();
546 
547  MoFEMFunctionReturn(0);
548 }

◆ refineMesh()

MoFEMErrorCode FreeSurface::refineMesh ( size_t  overlap)
private
Parameters
overlap
Returns
Examples
free_surface.cpp.

Definition at line 2575 of file free_surface.cpp.

2575  {
2576  MoFEMFunctionBegin;
2577 
2578  auto bit_mng = mField.getInterface<BitRefManager>();
2579 
2580  BitRefLevel start_mask;
2581  for (auto s = 0; s != get_start_bit(); ++s)
2582  start_mask[s] = true;
2583 
2584  // store prev_level
2585  Range prev_level;
2586  CHKERR bit_mng->getEntitiesByRefLevel(bit(get_current_bit()),
2587  BitRefLevel().set(), prev_level);
2588  Range prev_level_skin;
2589  CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_parent_bit()),
2590  BitRefLevel().set(), prev_level_skin);
2591  // reset bit ref levels
2592  CHKERR bit_mng->lambdaBitRefLevel(
2593  [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
2594  CHKERR bit_mng->setNthBitRefLevel(prev_level, get_projection_bit(), true);
2595  CHKERR bit_mng->setNthBitRefLevel(prev_level_skin, get_skin_projection_bit(),
2596  true);
2597 
2598  // set refinement levels
2599  auto set_levels = [&](auto &&
2600  vec_levels /*entities are refined on each level*/) {
2601  MoFEMFunctionBegin;
2602 
2603  // start with zero level, which is the coarsest mesh
2604  Range level0;
2605  CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
2606  CHKERR bit_mng->setNthBitRefLevel(level0, get_start_bit(), true);
2607 
2608  // get lower dimension entities
2609  auto get_adj = [&](auto ents) {
2610  Range conn;
2611  CHK_MOAB_THROW(mField.get_moab().get_connectivity(ents, conn, true),
2612  "get conn");
2613  for (auto d = 1; d != SPACE_DIM; ++d) {
2614  CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
2615  ents.subset_by_dimension(SPACE_DIM), d, false, conn,
2616  moab::Interface::UNION),
2617  "get adj");
2618  }
2619  ents.merge(conn);
2620  return ents;
2621  };
2622 
2623  // set bit levels
2624  for (auto l = 1; l != nb_levels; ++l) {
2625  Range level_prev;
2626  CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_start_bit() + l - 1),
2627  BitRefLevel().set(),
2628  SPACE_DIM, level_prev);
2629  Range parents;
2630  CHKERR bit_mng->updateRangeByParent(vec_levels[l], parents);
2631  // subtract entities from previous level, which are refined, so should be
2632  // not there
2633  level_prev = subtract(level_prev, parents);
2634  // and instead add their children
2635  level_prev.merge(vec_levels[l]);
2636  // set bit to each level
2637  CHKERR bit_mng->setNthBitRefLevel(level_prev, get_start_bit() + l, true);
2638  }
2639 
2640  // set bit levels to lower dimension entities
2641  for (auto l = 1; l != nb_levels; ++l) {
2642  Range level;
2643  CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2644  bit(get_start_bit() + l), BitRefLevel().set(), SPACE_DIM, level);
2645  level = get_adj(level);
2646  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(level);
2647  CHKERR bit_mng->setNthBitRefLevel(level, get_start_bit() + l, true);
2648  }
2649 
2650  MoFEMFunctionReturn(0);
2651  };
2652 
2653  // resolve skin between refined levels
2654  auto set_skins = [&]() {
2655  MoFEMFunctionBegin;
2656 
2657  moab::Skinner skinner(&mField.get_moab());
2658  ParallelComm *pcomm =
2659  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2660 
2661  // get skin of bit level
2662  auto get_bit_skin = [&](BitRefLevel bit, BitRefLevel mask) {
2663  Range bit_ents;
2664  CHK_THROW_MESSAGE(
2665  mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
2666  bit, mask, SPACE_DIM, bit_ents),
2667  "can't get bit level");
2668  Range bit_skin;
2669  CHK_MOAB_THROW(skinner.find_skin(0, bit_ents, false, bit_skin),
2670  "can't get skin");
2671  return bit_skin;
2672  };
2673 
2674  auto get_level_skin = [&]() {
2675  Range skin;
2676  BitRefLevel bit_prev;
2677  for (auto l = 1; l != nb_levels; ++l) {
2678  auto skin_level_mesh = get_bit_skin(bit(l), BitRefLevel().set());
2679  // filter (remove) all entities which are on partitions boarders
2680  CHKERR pcomm->filter_pstatus(skin_level_mesh,
2681  PSTATUS_SHARED | PSTATUS_MULTISHARED,
2682  PSTATUS_NOT, -1, nullptr);
2683  auto skin_level =
2684  get_bit_skin(bit(get_start_bit() + l), BitRefLevel().set());
2685  skin_level = subtract(skin_level,
2686  skin_level_mesh); // get only internal skins, not
2687  // on the body boundary
2688  // get lower dimension adjacencies (FIXME: add edges if 3D)
2689  Range skin_level_verts;
2690  CHKERR mField.get_moab().get_connectivity(skin_level, skin_level_verts,
2691  true);
2692  skin_level.merge(skin_level_verts);
2693 
2694  // remove previous level
2695  bit_prev.set(l - 1);
2696  Range level_prev;
2697  CHKERR bit_mng->getEntitiesByRefLevel(bit_prev, BitRefLevel().set(),
2698  level_prev);
2699  skin.merge(subtract(skin_level, level_prev));
2700  }
2701 
2702  return skin;
2703  };
2704 
2705  auto resolve_shared = [&](auto &&skin) {
2706  Range tmp_skin = skin;
2707 
2708  map<int, Range> map_procs;
2709  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2710  tmp_skin, &map_procs);
2711 
2712  Range from_other_procs; // entities which also exist on other processors
2713  for (auto &m : map_procs) {
2714  if (m.first != mField.get_comm_rank()) {
2715  from_other_procs.merge(m.second);
2716  }
2717  }
2718 
2719  auto common = intersect(
2720  skin, from_other_procs); // entities which are on internal skin
2721  skin.merge(from_other_procs);
2722 
2723  // entities which are on processors boards, and several processors are not
2724  // true skin.
2725  if (!common.empty()) {
2726  // skin is internal exist on other procs
2727  skin = subtract(skin, common);
2728  }
2729 
2730  return skin;
2731  };
2732 
2733  // get parents of entities
2734  auto get_parent_level_skin = [&](auto skin) {
2735  Range skin_parents;
2736  CHKERR bit_mng->updateRangeByParent(
2737  skin.subset_by_dimension(SPACE_DIM - 1), skin_parents);
2738  Range skin_parent_verts;
2739  CHKERR mField.get_moab().get_connectivity(skin_parents, skin_parent_verts,
2740  true);
2741  skin_parents.merge(skin_parent_verts);
2742  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
2743  skin_parents);
2744  return skin_parents;
2745  };
2746 
2747  auto child_skin = resolve_shared(get_level_skin());
2748  auto parent_skin = get_parent_level_skin(child_skin);
2749 
2750  child_skin = subtract(child_skin, parent_skin);
2751  CHKERR bit_mng->setNthBitRefLevel(child_skin, get_skin_child_bit(), true);
2752  CHKERR bit_mng->setNthBitRefLevel(parent_skin, get_skin_parent_bit(), true);
2753 
2754  MoFEMFunctionReturn(0);
2755  };
2756 
2757  // take last level, remove childs on boarder, and set bit
2758  auto set_current = [&]() {
2759  MoFEMFunctionBegin;
2760  Range last_level;
2761  CHKERR bit_mng->getEntitiesByRefLevel(bit(get_start_bit() + nb_levels - 1),
2762  BitRefLevel().set(), last_level);
2763  Range skin_child;
2764  CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_child_bit()),
2765  BitRefLevel().set(), skin_child);
2766 
2767  last_level = subtract(last_level, skin_child);
2768  CHKERR bit_mng->setNthBitRefLevel(last_level, get_current_bit(), true);
2769  MoFEMFunctionReturn(0);
2770  };
2771 
2772  // set bits to levels
2773  CHKERR set_levels(findEntitiesCrossedByPhaseInterface(overlap));
2774  // set bits to skin
2775  CHKERR set_skins();
2776  // set current level bit
2777  CHKERR set_current();
2778 
2779  if constexpr (debug) {
2780  if (mField.get_comm_size() == 1) {
2781  for (auto l = 0; l != nb_levels; ++l) {
2782  std::string name = (boost::format("out_level%d.h5m") % l).str();
2783  CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_start_bit() + l),
2784  BitRefLevel().set(), name.c_str(), "MOAB",
2785  "PARALLEL=WRITE_PART");
2786  }
2787  CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_current_bit()),
2788  BitRefLevel().set(), "current_bit.h5m",
2789  "MOAB", "PARALLEL=WRITE_PART");
2790  CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_projection_bit()),
2791  BitRefLevel().set(), "projection_bit.h5m",
2792  "MOAB", "PARALLEL=WRITE_PART");
2793 
2794  CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_child_bit()),
2795  BitRefLevel().set(), "skin_child_bit.h5m",
2796  "MOAB", "PARALLEL=WRITE_PART");
2797  CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_parent_bit()),
2798  BitRefLevel().set(), "skin_parent_bit.h5m",
2799  "MOAB", "PARALLEL=WRITE_PART");
2800  }
2801  }
2802 
2803  MoFEMFunctionReturn(0);
2804 };

◆ runProblem()

MoFEMErrorCode FreeSurface::runProblem ( )

[Run programme]

Examples
free_surface.cpp.

Definition at line 524 of file free_surface.cpp.

524  {
525  MoFEMFunctionBegin;
526  CHKERR readMesh();
527  CHKERR setupProblem();
528  CHKERR boundaryCondition();
529  CHKERR assembleSystem();
530  CHKERR solveSystem();
531  MoFEMFunctionReturn(0);
532 }

◆ setParentDofs()

MoFEMErrorCode FreeSurface::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 
)
private

Create hierarchy of elements run on parents levels.

Parameters
fe_toppipeline element to which hierarchy is attached
optype of operator OPSPACE, OPROW, OPCOL or OPROWCOL
get_elemalambda function to create element on hierarchy
verbosityverbosity level
sevseverity level
Returns

Collect data from parent elements to child

Examples
free_surface.cpp.

Definition at line 2806 of file free_surface.cpp.

2811  {
2812  MoFEMFunctionBegin;
2813 
2814  /**
2815  * @brief Collect data from parent elements to child
2816  */
2817  boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
2818  add_parent_level =
2819  [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
2820  // Evaluate if not last parent element
2821  if (level > 0) {
2822 
2823  // Create domain parent FE
2824  auto fe_ptr_current = get_elem();
2825 
2826  // Call next level
2827  add_parent_level(
2828  boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
2829  fe_ptr_current),
2830  level - 1);
2831 
2832  // Add data to curent fe level
2833  if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
2834 
2835  // Only base
2836  parent_fe_pt->getOpPtrVector().push_back(
2837 
2838  new OpAddParentEntData(
2839 
2840  H1, op, fe_ptr_current,
2841 
2842  BitRefLevel().set(), BitRefLevel().set(),
2843 
2844  child_ent_bit, BitRefLevel().set(),
2845 
2846  verbosity, sev));
2847 
2848  } else {
2849 
2850  // Filed data
2851  parent_fe_pt->getOpPtrVector().push_back(
2852 
2853  new OpAddParentEntData(
2854 
2855  field_name, op, fe_ptr_current,
2856 
2857  BitRefLevel().set(), BitRefLevel().set(),
2858 
2859  child_ent_bit, BitRefLevel().set(),
2860 
2861  verbosity, sev));
2862  }
2863  }
2864  };
2865 
2866  add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
2867  nb_levels);
2868 
2869  MoFEMFunctionReturn(0);
2870 }

◆ setupProblem()

MoFEMErrorCode FreeSurface::setupProblem ( )
private

[Read mesh]

[Set up problem]

Examples
free_surface.cpp.

Definition at line 552 of file free_surface.cpp.

552  {
553  MoFEMFunctionBegin;
554 
555  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
556  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-nb_levels", &nb_levels,
557  PETSC_NULL);
558  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-refine_overlap", &refine_overlap,
559  PETSC_NULL);
560 
561  const char *coord_type_names[] = {"cartesian", "polar", "cylindrical",
562  "spherical"};
563  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-coords", coord_type_names,
564  LAST_COORDINATE_SYSTEM, &coord_type, PETSC_NULL);
565 
566  MOFEM_LOG("FS", Sev::inform) << "Approximation order = " << order;
567  MOFEM_LOG("FS", Sev::inform)
568  << "Number of refinement levels nb_levels = " << nb_levels;
569  nb_levels += 1;
570 
571  auto simple = mField.getInterface<Simple>();
572 
573  CHKERR PetscOptionsGetScalar(PETSC_NULL, "Acceleration", "-a0", &a0, PETSC_NULL);
574  CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Minus\" phase density rho_m",
575  "-rho_m", &rho_m, PETSC_NULL);
576  CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Minus\" phase viscosity", "-mu_m",
577  &mu_m, PETSC_NULL);
578  CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Plus\" phase density", "-rho_p",
579  &rho_p, PETSC_NULL);
580  CHKERR PetscOptionsGetScalar(PETSC_NULL, "\"Plus\" phase viscosity", "-mu_p",
581  &mu_p, PETSC_NULL);
582  CHKERR PetscOptionsGetScalar(PETSC_NULL, "Surface tension", "-lambda",
583  &lambda, PETSC_NULL);
584  CHKERR PetscOptionsGetScalar(PETSC_NULL,
585  "Height of the well in energy functional", "-W",
586  &W, PETSC_NULL);
587  rho_ave = (rho_p + rho_m) / 2;
588  rho_diff = (rho_p - rho_m) / 2;
589  mu_ave = (mu_p + mu_m) / 2;
590  mu_diff = (mu_p - mu_m) / 2;
591 
592  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-h", &h, PETSC_NULL);
593  eta = h;
594  eta2 = eta * eta;
595  kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
596 
597  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-md", &eta, PETSC_NULL);
598 
599  MOFEM_LOG("FS", Sev::inform) << "Acceleration a0 = " << a0;
600  MOFEM_LOG("FS", Sev::inform) << "\"Minus\" phase density rho_m = " << rho_m;
601  MOFEM_LOG("FS", Sev::inform) << "\"Minus\" phase viscosity mu_m = " << mu_m;
602  MOFEM_LOG("FS", Sev::inform) << "\"Plus\" phase density rho_p = " << rho_p;
603  MOFEM_LOG("FS", Sev::inform) << "\"Plus\" phase viscosity mu_p = " << mu_p;
604  MOFEM_LOG("FS", Sev::inform) << "Surface tension lambda = " << lambda;
605  MOFEM_LOG("FS", Sev::inform)
606  << "Height of the well in energy functional W = " << W;
607  MOFEM_LOG("FS", Sev::inform) << "Characteristic mesh size h = " << h;
608  MOFEM_LOG("FS", Sev::inform) << "Mobility md = " << md;
609 
610  MOFEM_LOG("FS", Sev::inform) << "Average density rho_ave = " << rho_ave;
611  MOFEM_LOG("FS", Sev::inform) << "Difference density rho_diff = " << rho_diff;
612  MOFEM_LOG("FS", Sev::inform) << "Average viscosity mu_ave = " << mu_ave;
613  MOFEM_LOG("FS", Sev::inform) << "Difference viscosity mu_diff = " << mu_diff;
614  MOFEM_LOG("FS", Sev::inform) << "kappa = " << kappa;
615 
616 
617  // Fields on domain
618 
619  // Velocity field
620  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, U_FIELD_DIM);
621  // Pressure field
622  CHKERR simple->addDomainField("P", H1, AINSWORTH_LEGENDRE_BASE, 1);
623  // Order/phase fild
624  CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
625  // Chemical potential
626  CHKERR simple->addDomainField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
627 
628  // Field on boundary
629  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE,
630  U_FIELD_DIM);
631  CHKERR simple->addBoundaryField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
632  CHKERR simple->addBoundaryField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
633  // Lagrange multiplier which constrains slip conditions
634  CHKERR simple->addBoundaryField("L", H1, AINSWORTH_LEGENDRE_BASE, 1);
635 
636  CHKERR simple->setFieldOrder("U", order);
637  CHKERR simple->setFieldOrder("P", order - 1);
638  CHKERR simple->setFieldOrder("H", order);
639  CHKERR simple->setFieldOrder("G", order);
640  CHKERR simple->setFieldOrder("L", order);
641 
642  // Initialise bit ref levels
643  auto set_problem_bit = [&]() {
644  MoFEMFunctionBegin;
645  // Set bits to build adjacencies between parents and children. That is
646  // used by simple interface.
647  simple->getBitAdjEnt() = BitRefLevel().set();
648  simple->getBitAdjParent() = BitRefLevel().set();
649  simple->getBitRefLevel() = BitRefLevel().set();
650  simple->getBitRefLevelMask() = BitRefLevel().set();
651  MoFEMFunctionReturn(0);
652  };
653 
654  CHKERR set_problem_bit();
655 
656  CHKERR simple->setUp();
657 
658  MoFEMFunctionReturn(0);
659 }

◆ solveSystem()

MoFEMErrorCode FreeSurface::solveSystem ( )
private

[Solve]

Examples
free_surface.cpp.

Definition at line 2044 of file free_surface.cpp.

2044  {
2045  MoFEMFunctionBegin;
2046 
2047  using UDO = ForcesAndSourcesCore::UserDataOperator;
2048 
2049  auto simple = mField.getInterface<Simple>();
2050  auto pip_mng = mField.getInterface<PipelineManager>();
2051 
2052  auto create_solver_dm = [&](auto dm) -> SmartPetscObj<DM> {
2053  DM subdm;
2054 
2055  auto setup_subdm = [&](auto dm) {
2056  MoFEMFunctionBegin;
2057  auto simple = mField.getInterface<Simple>();
2058  auto bit_mng = mField.getInterface<BitRefManager>();
2059  auto dm = simple->getDM();
2060  auto level_ents_ptr = boost::make_shared<Range>();
2061  CHKERR bit_mng->getEntitiesByRefLevel(
2062  bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
2063  CHKERR DMCreate(mField.get_comm(), &subdm);
2064  CHKERR DMSetType(subdm, "DMMOFEM");
2065  CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
2066  CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
2067  CHKERR DMMoFEMAddElement(subdm, simple->getBoundaryFEName());
2068  CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
2069  CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_FALSE);
2070  for (auto f : {"U", "P", "H", "G", "L"}) {
2071  CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
2072  CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
2073  }
2074  CHKERR DMSetUp(subdm);
2075  MoFEMFunctionReturn(0);
2076  };
2077 
2078  CHK_THROW_MESSAGE(setup_subdm(dm), "create subdm");
2079 
2080  return SmartPetscObj<DM>(subdm);
2081  };
2082 
2083  auto get_fe_post_proc = [&](auto post_proc_mesh) {
2084  auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
2085  return setParentDofs(
2086  fe, field, op, bit(get_skin_parent_bit()),
2087 
2088  [&]() {
2089  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2090  new DomianParentEle(mField));
2091  return fe_parent;
2092  },
2093 
2094  QUIET, Sev::noisy);
2095  };
2096 
2097  auto post_proc_fe =
2098  boost::make_shared<PostProcEleDomainCont>(mField, post_proc_mesh);
2099  post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
2100  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2101  get_start_bit() + nb_levels - 1);
2102  };
2103 
2104  auto u_ptr = boost::make_shared<MatrixDouble>();
2105  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2106  auto h_ptr = boost::make_shared<VectorDouble>();
2107  auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2108  auto p_ptr = boost::make_shared<VectorDouble>();
2109  auto g_ptr = boost::make_shared<VectorDouble>();
2110  auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2111 
2112  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2113  post_proc_fe->getOpPtrVector(), {H1});
2114 
2115  CHKERR setParentDofs(
2116  post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2117 
2118  [&]() {
2119  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2120  new DomianParentEle(mField));
2121  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2122  fe_parent->getOpPtrVector(), {H1});
2123  return fe_parent;
2124  },
2125 
2126  QUIET, Sev::noisy);
2127 
2128  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U");
2129  post_proc_fe->getOpPtrVector().push_back(
2130  new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
2131  post_proc_fe->getOpPtrVector().push_back(
2132  new OpCalculateVectorFieldGradient<U_FIELD_DIM, SPACE_DIM>("U",
2133  grad_u_ptr));
2134 
2135  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H");
2136  post_proc_fe->getOpPtrVector().push_back(
2137  new OpCalculateScalarFieldValues("H", h_ptr));
2138  post_proc_fe->getOpPtrVector().push_back(
2139  new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
2140 
2141  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P");
2142  post_proc_fe->getOpPtrVector().push_back(
2143  new OpCalculateScalarFieldValues("P", p_ptr));
2144 
2145  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G");
2146  post_proc_fe->getOpPtrVector().push_back(
2147  new OpCalculateScalarFieldValues("G", g_ptr));
2148  post_proc_fe->getOpPtrVector().push_back(
2149  new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
2150 
2151  using OpPPMap = OpPostProcMapInMoab<2, 2>;
2152 
2153  post_proc_fe->getOpPtrVector().push_back(
2154 
2155  new OpPPMap(
2156  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2157 
2158  {{"H", h_ptr}, {"P", p_ptr}, {"G", g_ptr}},
2159 
2160  {{"U", u_ptr}, {"H_GRAD", grad_h_ptr}, {"G_GRAD", grad_g_ptr}},
2161 
2162  {{"GRAD_U", grad_u_ptr}},
2163 
2164  {}
2165 
2166  )
2167 
2168  );
2169 
2170  return post_proc_fe;
2171  };
2172 
2173  auto get_bdy_post_proc_fe = [&](auto post_proc_mesh) {
2174  auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2175  return setParentDofs(
2176  fe, field, op, bit(get_skin_parent_bit()),
2177 
2178  [&]() {
2179  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2180  new BoundaryParentEle(mField));
2181  return fe_parent;
2182  },
2183 
2184  QUIET, Sev::noisy);
2185  };
2186 
2187  auto post_proc_fe =
2188  boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
2189  post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
2190  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2191  get_start_bit() + nb_levels - 1);
2192  };
2193 
2194  CHKERR setParentDofs(
2195  post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2196 
2197  [&]() {
2198  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2199  new BoundaryParentEle(mField));
2200  return fe_parent;
2201  },
2202 
2203  QUIET, Sev::noisy);
2204 
2205  struct OpGetNormal : public BoundaryEleOp {
2206 
2207  OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
2208  boost::shared_ptr<MatrixDouble> n_ptr)
2209  : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), ptrL(l_ptr),
2210  ptrNormal(n_ptr) {}
2211 
2212  MoFEMErrorCode doWork(int side, EntityType type,
2213  EntitiesFieldData::EntData &data) {
2214  MoFEMFunctionBegin;
2215  auto t_l = getFTensor0FromVec(*ptrL);
2216  auto t_n_fe = getFTensor1NormalsAtGaussPts();
2217  ptrNormal->resize(SPACE_DIM, getGaussPts().size2(), false);
2218  auto t_n = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
2219  for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
2220  t_n(i) = t_n_fe(i) * t_l / std::sqrt(t_n_fe(i) * t_n_fe(i));
2221  ++t_n_fe;
2222  ++t_l;
2223  ++t_n;
2224  }
2225  MoFEMFunctionReturn(0);
2226  };
2227 
2228  protected:
2229  boost::shared_ptr<VectorDouble> ptrL;
2230  boost::shared_ptr<MatrixDouble> ptrNormal;
2231  };
2232 
2233  auto u_ptr = boost::make_shared<MatrixDouble>();
2234  auto p_ptr = boost::make_shared<VectorDouble>();
2235  auto lambda_ptr = boost::make_shared<VectorDouble>();
2236  auto normal_l_ptr = boost::make_shared<MatrixDouble>();
2237 
2238  CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "U");
2239  post_proc_fe->getOpPtrVector().push_back(
2240  new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
2241 
2242  CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "L");
2243  post_proc_fe->getOpPtrVector().push_back(
2244  new OpCalculateScalarFieldValues("L", lambda_ptr));
2245  post_proc_fe->getOpPtrVector().push_back(
2246  new OpGetNormal(lambda_ptr, normal_l_ptr));
2247 
2248  CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "P");
2249  post_proc_fe->getOpPtrVector().push_back(
2250  new OpCalculateScalarFieldValues("P", p_ptr));
2251 
2252  auto op_ptr = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
2253  op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
2254  EntityType type,
2255  EntitiesFieldData::EntData &data) {
2256  MoFEMFunctionBegin;
2257  MoFEMFunctionReturn(0);
2258  };
2259 
2260  using OpPPMap = OpPostProcMapInMoab<2, 2>;
2261 
2262  post_proc_fe->getOpPtrVector().push_back(
2263 
2264  new OpPPMap(post_proc_fe->getPostProcMesh(),
2265  post_proc_fe->getMapGaussPts(),
2266 
2267  OpPPMap::DataMapVec{{"P", p_ptr}},
2268 
2269  OpPPMap::DataMapMat{{"U", u_ptr}, {"L", normal_l_ptr}},
2270 
2271  OpPPMap::DataMapMat(),
2272 
2273  OpPPMap::DataMapMat()
2274 
2275  )
2276 
2277  );
2278 
2279  return post_proc_fe;
2280  };
2281 
2282  auto get_lift_fe = [&]() {
2283  auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2284  return setParentDofs(
2285  fe, field, op, bit(get_skin_parent_bit()),
2286 
2287  [&]() {
2288  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2289  new BoundaryParentEle(mField));
2290  return fe_parent;
2291  },
2292 
2293  QUIET, Sev::noisy);
2294  };
2295 
2296  auto fe = boost::make_shared<BoundaryEle>(mField);
2297  fe->exeTestHook = [](FEMethod *fe_ptr) {
2298  return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2299  get_start_bit() + nb_levels - 1);
2300  };
2301 
2302  auto lift_ptr = boost::make_shared<VectorDouble>();
2303  auto p_ptr = boost::make_shared<VectorDouble>();
2304  auto ents_ptr = boost::make_shared<Range>();
2305 
2306  CHKERR setParentDofs(
2307  fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2308 
2309  [&]() {
2310  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2311  new BoundaryParentEle(mField));
2312  return fe_parent;
2313  },
2314 
2315  QUIET, Sev::noisy);
2316 
2317  std::vector<const CubitMeshSets *> vec_ptr;
2318  CHKERR mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
2319  std::regex("LIFT"), vec_ptr);
2320  for (auto m_ptr : vec_ptr) {
2321  auto meshset = m_ptr->getMeshset();
2322  Range ents;
2323  CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
2324  ents, true);
2325  ents_ptr->merge(ents);
2326  }
2327 
2328  MOFEM_LOG("FS", Sev::noisy) << "Lift ents " << (*ents_ptr);
2329 
2330  CHKERR add_parent_field_bdy(fe, UDO::OPROW, "P");
2331  fe->getOpPtrVector().push_back(
2332  new OpCalculateScalarFieldValues("L", p_ptr));
2333  fe->getOpPtrVector().push_back(
2334  new OpCalculateLift("L", p_ptr, lift_ptr, ents_ptr));
2335 
2336  return std::make_pair(fe, lift_ptr);
2337  };
2338 
2339  auto set_post_proc_monitor = [&](auto ts) {
2340  MoFEMFunctionBegin;
2341  DM dm;
2342  CHKERR TSGetDM(ts, &dm);
2343  boost::shared_ptr<FEMethod> null_fe;
2344  auto post_proc_mesh = boost::make_shared<moab::Core>();
2345  auto monitor_ptr = boost::make_shared<Monitor>(
2346  SmartPetscObj<DM>(dm, true), post_proc_mesh,
2347  get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
2348  get_lift_fe());
2349  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
2350  null_fe, monitor_ptr);
2351  MoFEMFunctionReturn(0);
2352  };
2353 
2354  auto dm = simple->getDM();
2355  auto ts = createTS(mField.get_comm());
2356  CHKERR TSSetDM(ts, dm);
2357 
2358  auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
2359  tsPrePostProc = ts_pre_post_proc;
2360 
2361  if (auto ptr = tsPrePostProc.lock()) {
2362 
2363  ptr->fsRawPtr = this;
2364  ptr->solverSubDM = create_solver_dm(simple->getDM());
2365  ptr->globSol = createDMVector(dm);
2366  CHKERR DMoFEMMeshToLocalVector(dm, ptr->globSol, INSERT_VALUES,
2367  SCATTER_FORWARD);
2368  CHKERR VecAssemblyBegin(ptr->globSol);
2369  CHKERR VecAssemblyEnd(ptr->globSol);
2370 
2371  auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
2372 
2373  CHKERR set_post_proc_monitor(sub_ts);
2374 
2375  // Add monitor to time solver
2376  CHKERR TSSetFromOptions(ts);
2377  CHKERR ptr->tsSetUp(ts);
2378  CHKERR TSSetUp(ts);
2379 
2380  auto print_fields_in_section = [&]() {
2381  MoFEMFunctionBegin;
2382  auto section = mField.getInterface<ISManager>()->sectionCreate(
2383  simple->getProblemName());
2384  PetscInt num_fields;
2385  CHKERR PetscSectionGetNumFields(section, &num_fields);
2386  for (int f = 0; f < num_fields; ++f) {
2387  const char *field_name;
2388  CHKERR PetscSectionGetFieldName(section, f, &field_name);
2389  MOFEM_LOG("FS", Sev::inform)
2390  << "Field " << f << " " << std::string(field_name);
2391  }
2392  MoFEMFunctionReturn(0);
2393  };
2394 
2395  CHKERR print_fields_in_section();
2396 
2397  CHKERR TSSolve(ts, ptr->globSol);
2398  }
2399 
2400  MoFEMFunctionReturn(0);
2401 }

Friends And Related Function Documentation

◆ TSPrePostProc

friend struct TSPrePostProc
friend

Definition at line 520 of file free_surface.cpp.

Member Data Documentation

◆ mField

MoFEM::Interface& FreeSurface::mField
Examples
free_surface.cpp.

Definition at line 479 of file free_surface.cpp.

The documentation for this struct was generated from the following file:
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
rho_ave
double rho_ave
Definition: free_surface.cpp:180
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:712
FreeSurface::findEntitiesCrossedByPhaseInterface
std::vector< Range > findEntitiesCrossedByPhaseInterface(size_t overlap)
Find entities on refinement levels.
Definition: free_surface.cpp:2453
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
get_fe_bit
auto get_fe_bit
Definition: free_surface.cpp:316
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: free_surface.cpp:134
g
constexpr double g
Definition: shallow_wave.cpp:63
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
H1
@ H1
continuous field
Definition: definitions.h:85
FreeSurfaceOps::OpNormalForceRhs
Definition: FreeSurfaceOps.hpp:92
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:24
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
CARTESIAN
@ CARTESIAN
Definition: definitions.h:128
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
MoFEM::OpCalculateDivergenceVectorFieldValues
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:488
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
FreeSurfaceOps::OpWettingAngleRhs
Definition: FreeSurfaceOps.hpp:138
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:373
get_skin_child_bit
auto get_skin_child_bit
Definition: free_surface.cpp:155
MoFEM::OpRunParent
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Definition: MeshProjectionDataOperators.hpp:17
eps
double eps
Definition: free_surface.cpp:177
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition: PetscSmartObj.hpp:249
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
get_skin_parent_bit
auto get_skin_parent_bit
Definition: free_surface.cpp:154
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
mu_m
double mu_m
Definition: free_surface.cpp:164
FreeSurfaceOps::OpWettingAngleLhs
Definition: FreeSurfaceOps.hpp:246
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FreeSurface::assembleSystem
MoFEMErrorCode assembleSystem()
[Data projection]
Definition: free_surface.cpp:1581
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468
MoFEM::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:176
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
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
OpBoundaryAssembleScalar
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
Definition: free_surface.cpp:123
W
double W
Definition: free_surface.cpp:168
FreeSurfaceOps::OpLhsH_dG
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1117
l
FTensor::Index< 'l', SPACE_DIM > l
Definition: free_surface.cpp:137
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
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
debug
constexpr bool debug
Definition: free_surface.cpp:146
sdf.r
int r
Definition: sdf.py:8
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
OpMixScalarTimesDiv
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
Definition: free_surface.cpp:127
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1319
MoFEM::CoreInterface::get_dofs
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
MoFEM::MeshsetsManager::getMeshset
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:719
h
double h
Definition: free_surface.cpp:171
FreeSurfaceOps::OpRhsH
Definition: FreeSurfaceOps.hpp:795
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:597
cylindrical
auto cylindrical
Definition: free_surface.cpp:189
OpDomainAssembleVector
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
Definition: free_surface.cpp:119
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:438
U_FIELD_DIM
constexpr int U_FIELD_DIM
Definition: free_surface.cpp:81
OpBoundaryMassL
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
Definition: free_surface.cpp:116
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FreeSurfaceOps::OpLhsH_dU
Lhs for H dU.
Definition: FreeSurfaceOps.hpp:915
coord_type
int coord_type
Definition: free_surface.cpp:77
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
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
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
BoundaryEleOp
save_range
auto save_range
Definition: free_surface.cpp:326
MoFEM::BitRefManager::getEntitiesByDimAndRefLevel
MoFEMErrorCode getEntitiesByDimAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const int dim, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
Definition: BitRefManager.cpp:822
rho_diff
double rho_diff
Definition: free_surface.cpp:181
get_current_bit
auto get_current_bit
dofs bit used to do calculations
Definition: free_surface.cpp:151
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:215
double
convert.type
type
Definition: convert.py:64
OpDomainMassH
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
Definition: free_surface.cpp:112
get_global_size
auto get_global_size
Definition: free_surface.cpp:320
FreeSurfaceOps::OpLhsH_dH
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:973
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
nb_levels
int nb_levels
Definition: free_surface.cpp:143
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
eta2
double eta2
Definition: free_surface.cpp:173
BoundaryParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
Definition: free_surface.cpp:94
OpDomainMassP
OpDomainMassH OpDomainMassP
Definition: free_surface.cpp:113
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
get_projection_bit
auto get_projection_bit
Definition: free_surface.cpp:156
eta
double eta
Definition: free_surface.cpp:172
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
bit
auto bit
Definition: free_surface.cpp:314
FreeSurfaceOps::OpRhsG
Definition: FreeSurfaceOps.hpp:1181
SPACE_DIM
constexpr int SPACE_DIM
Definition: free_surface.cpp:80
MoFEM::solve
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
Definition: Schur.cpp:1296
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
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
FreeSurfaceOps::OpLhsU_dG
Lhs for G dH.
Definition: FreeSurfaceOps.hpp:732
get_skin_projection_bit
auto get_skin_projection_bit
Definition: free_surface.cpp:157
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
order
int order
approximation order
Definition: free_surface.cpp:142
mu_diff
double mu_diff
Definition: free_surface.cpp:183
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:417
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1561
convert.n
n
Definition: convert.py:82
rho_p
double rho_p
Definition: free_surface.cpp:165
FreeSurfaceOps::OpLhsG_dG
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1321
kappa
double kappa
Definition: free_surface.cpp:185
rho_m
double rho_m
Definition: free_surface.cpp:163
integration_rule
auto integration_rule
Definition: free_surface.cpp:187
tsPrePostProc
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
Definition: free_surface.cpp:469
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
DomainEleOp
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
get_dofs_ents_all
auto get_dofs_ents_all
get entities of dofs in the problem - used for debugging
Definition: free_surface.cpp:370
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CYLINDRICAL
@ CYLINDRICAL
Definition: definitions.h:130
lambda
double lambda
surface tension
Definition: free_surface.cpp:167
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
FreeSurface::setParentDofs
MoFEMErrorCode setParentDofs(boost::shared_ptr< FEMethod > fe_top, std::string field_name, ForcesAndSourcesCore::UserDataOperator::OpType op, BitRefLevel child_ent_bit, boost::function< boost::shared_ptr< ForcesAndSourcesCore >()> get_elem, int verbosity, LogManager::SeverityLevel sev)
Create hierarchy of elements run on parents levels.
Definition: free_surface.cpp:2806
FreeSurfaceOps::OpLhsU_dU
Lhs for U dU.
Definition: FreeSurfaceOps.hpp:473
mu_ave
double mu_ave
Definition: free_surface.cpp:182
FreeSurfaceOps::OpNormalConstrainRhs
Definition: FreeSurfaceOps.hpp:50
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
OpDomainMassU
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
Definition: free_surface.cpp:110
refine_overlap
int refine_overlap
Definition: free_surface.cpp:144
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
MoFEM::createVecScatter
auto createVecScatter(Vec x, IS ix, Vec y, IS iy)
Create a Vec Scatter object.
Definition: PetscSmartObj.hpp:361
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
FreeSurface::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: free_surface.cpp:2044
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:64
get_start_bit
auto get_start_bit
Definition: free_surface.cpp:148
FreeSurface::readMesh
MoFEMErrorCode readMesh()
[Run programme]
Definition: free_surface.cpp:536
LAST_COORDINATE_SYSTEM
@ LAST_COORDINATE_SYSTEM
Definition: definitions.h:132
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
FreeSurfaceOps::OpRhsU
Rhs for U.
Definition: FreeSurfaceOps.hpp:356
FreeSurfaceOps::OpCalculateLift
Definition: FreeSurfaceOps.hpp:3
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
BoundaryEleOp
BoundaryEle::UserDataOperator BoundaryEleOp
Definition: free_surface.cpp:95
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:239
FreeSurfaceOps::OpLhsU_dH
Lhs for U dH.
Definition: FreeSurfaceOps.hpp:611
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
a0
double a0
Definition: free_surface.cpp:162
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1127
FreeSurface::mField
MoFEM::Interface & mField
Definition: free_surface.cpp:479
QUIET
@ QUIET
Definition: definitions.h:221
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:18
mu_p
double mu_p
Definition: free_surface.cpp:166
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
FreeSurfaceOps::OpLhsG_dH
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1249
MoFEM::SmartPetscObj< DM >
FreeSurface::refineMesh
MoFEMErrorCode refineMesh(size_t overlap)
Definition: free_surface.cpp:2575
FreeSurfaceOps::OpNormalConstrainLhs
Definition: FreeSurfaceOps.hpp:195
FreeSurface::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: free_surface.cpp:552
OpDomainAssembleScalar
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
Definition: free_surface.cpp:121
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
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
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
DomianParentEle
ElementsAndOps< SPACE_DIM >::DomianParentEle DomianParentEle
Definition: free_surface.cpp:91
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1291
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:711
OpDomainMassG
OpDomainMassH OpDomainMassG
Definition: free_surface.cpp:114
FreeSurface::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: free_surface.cpp:663
F
@ F
Definition: free_surface.cpp:396
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:709
Generated by Doxygen 1.8.17 and hosted at University of Glasgow