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 469 of file free_surface.cpp.

Constructor & Destructor Documentation

◆ FreeSurface()

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

Definition at line 471 of file free_surface.cpp.

471 : mField(m_field) {}

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode FreeSurface::assembleSystem ( )
private

[Data projection]

[Push operators to pip]

Examples
free_surface.cpp.

Definition at line 1582 of file free_surface.cpp.

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

◆ boundaryCondition()

MoFEMErrorCode FreeSurface::boundaryCondition ( )
private

[Set up problem]

[Boundary condition]

Examples
free_surface.cpp.

Definition at line 662 of file free_surface.cpp.

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

◆ 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 2456 of file free_surface.cpp.

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

◆ 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 952 of file free_surface.cpp.

952  {
953  MoFEMFunctionBegin;
954 
955  // FIXME: Functionality in this method should be improved (projection should
956  // be done field by field), generalised, and move to become core
957  // functionality.
958 
959  auto simple = mField.getInterface<Simple>();
960  auto pip_mng = mField.getInterface<PipelineManager>();
961  auto bc_mng = mField.getInterface<BcManager>();
962  auto bit_mng = mField.getInterface<BitRefManager>();
963  auto field_blas = mField.getInterface<FieldBlas>();
964 
965  // Store all existing elements pipelines, replace them by data projection
966  // pipelines, to put them back when projection is done.
967  auto fe_domain_lhs = pip_mng->getDomainLhsFE();
968  auto fe_domain_rhs = pip_mng->getDomainRhsFE();
969  auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
970  auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
971 
972  pip_mng->getDomainLhsFE().reset();
973  pip_mng->getDomainRhsFE().reset();
974  pip_mng->getBoundaryLhsFE().reset();
975  pip_mng->getBoundaryRhsFE().reset();
976 
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  int ii = 0;
1460  for (auto v : ts_solver_vecs) {
1461  MOFEM_LOG("FS", Sev::inform) << "Solve projection vector";
1462 
1463  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1464  SCATTER_REVERSE);
1465  CHKERR solve(sub_x);
1466 
1467  for (auto f : {"U", "P", "H", "G", "L"}) {
1468  MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1469  CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1470  }
1471 
1472  CHKERR DMoFEMMeshToLocalVector(subdm, sub_x, INSERT_VALUES,
1473  SCATTER_REVERSE);
1474  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1475  SCATTER_FORWARD);
1476 
1477  MOFEM_LOG("FS", Sev::inform) << "Norm V " << get_norm(v);
1478  }
1479 
1480  CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_X0",
1481  &ts_solver_vecs[0]);
1482  CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_Xdot",
1483  &ts_solver_vecs[1]);
1484  }
1485 
1486  for (auto f : {"U", "P", "H", "G", "L"}) {
1487  MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1488  CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1489  }
1490  CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
1491  }
1492 
1493  MoFEMFunctionReturn(0);
1494  };
1495 
1496  // postprocessing (only for debugging proposes)
1497  auto post_proc = [&](auto exe_test) {
1498  MoFEMFunctionBegin;
1499  auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
1500  auto &pip = post_proc_fe->getOpPtrVector();
1501  post_proc_fe->exeTestHook = exe_test;
1502 
1503  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
1504 
1505  CHKERR setParentDofs(
1506  post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1507 
1508  [&]() {
1509  boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1510  new DomianParentEle(mField));
1511  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1512  fe_parent->getOpPtrVector(), {H1});
1513  return fe_parent;
1514  },
1515 
1516  QUIET, Sev::noisy);
1517 
1518  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
1519 
1520  auto u_ptr = boost::make_shared<MatrixDouble>();
1521  auto p_ptr = boost::make_shared<VectorDouble>();
1522  auto h_ptr = boost::make_shared<VectorDouble>();
1523  auto g_ptr = boost::make_shared<VectorDouble>();
1524 
1525  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U",
1526  bit(get_skin_parent_bit()));
1527  pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1528  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P",
1529  bit(get_skin_parent_bit()));
1530  pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1531  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H",
1532  bit(get_skin_parent_bit()));
1533  pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1534  CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G",
1535  bit(get_skin_parent_bit()));
1536  pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1537 
1538  post_proc_fe->getOpPtrVector().push_back(
1539 
1540  new OpPPMap(post_proc_fe->getPostProcMesh(),
1541  post_proc_fe->getMapGaussPts(),
1542 
1543  {{"P", p_ptr}, {"H", h_ptr}, {"G", g_ptr}},
1544 
1545  {{"U", u_ptr}},
1546 
1547  {},
1548 
1549  {}
1550 
1551  )
1552 
1553  );
1554 
1555  auto dm = simple->getDM();
1556  CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
1557  CHKERR post_proc_fe->writeFile("out_projection.h5m");
1558 
1559  MoFEMFunctionReturn(0);
1560  };
1561 
1562  CHKERR solve_projection([](FEMethod *fe_ptr) {
1563  return get_fe_bit(fe_ptr).test(get_current_bit());
1564  });
1565 
1566  if constexpr (debug) {
1567  CHKERR post_proc([](FEMethod *fe_ptr) {
1568  return get_fe_bit(fe_ptr).test(get_current_bit());
1569  });
1570  }
1571 
1572  fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
1573  fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
1574  fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
1575  fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
1576 
1577  MoFEMFunctionReturn(0);
1578 }

◆ readMesh()

MoFEMErrorCode FreeSurface::readMesh ( )
private

[Run programme]

[Read mesh]

Examples
free_surface.cpp.

Definition at line 534 of file free_surface.cpp.

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

◆ refineMesh()

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

Definition at line 2578 of file free_surface.cpp.

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

◆ runProblem()

MoFEMErrorCode FreeSurface::runProblem ( )

[Run programme]

Examples
free_surface.cpp.

Definition at line 522 of file free_surface.cpp.

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

◆ 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 2811 of file free_surface.cpp.

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

◆ setupProblem()

MoFEMErrorCode FreeSurface::setupProblem ( )
private

[Read mesh]

[Set up problem]

Examples
free_surface.cpp.

Definition at line 550 of file free_surface.cpp.

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

◆ solveSystem()

MoFEMErrorCode FreeSurface::solveSystem ( )
private

[Solve]

Examples
free_surface.cpp.

Definition at line 2045 of file free_surface.cpp.

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

Friends And Related Function Documentation

◆ TSPrePostProc

friend struct TSPrePostProc
friend

Definition at line 518 of file free_surface.cpp.

Member Data Documentation

◆ mField

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

Definition at line 477 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 refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
rho_ave
double rho_ave
Definition: free_surface.cpp:178
FreeSurface::findEntitiesCrossedByPhaseInterface
std::vector< Range > findEntitiesCrossedByPhaseInterface(size_t overlap)
Find entities on refinement levels.
Definition: free_surface.cpp:2456
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
get_fe_bit
auto get_fe_bit
Definition: free_surface.cpp:314
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: free_surface.cpp:132
g
constexpr double g
Definition: shallow_wave.cpp:63
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
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:284
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
sdf_hertz.d
float d
Definition: sdf_hertz.py:5
CARTESIAN
@ CARTESIAN
Definition: definitions.h:115
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPSPACE
@ OPSPACE
operator do Work is execute on space data
Definition: ForcesAndSourcesCore.hpp:570
MoFEM::OpCalculateDivergenceVectorFieldValues
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:486
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
FreeSurfaceOps::OpWettingAngleRhs
Definition: FreeSurfaceOps.hpp:138
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
get_skin_child_bit
auto get_skin_child_bit
Definition: free_surface.cpp:153
MoFEM::OpRunParent
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Definition: MeshProjectionDataOperators.hpp:17
eps
double eps
Definition: free_surface.cpp:175
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition: PetscSmartObj.hpp:245
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:700
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
get_skin_parent_bit
auto get_skin_parent_bit
Definition: free_surface.cpp:152
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:460
mu_m
double mu_m
Definition: free_surface.cpp:162
FreeSurfaceOps::OpWettingAngleLhs
Definition: FreeSurfaceOps.hpp:246
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FreeSurface::assembleSystem
MoFEMErrorCode assembleSystem()
[Data projection]
Definition: free_surface.cpp:1582
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
md
double md
Definition: free_surface.cpp:174
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
MoFEM::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:527
OpBoundaryAssembleScalar
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
Definition: free_surface.cpp:121
W
double W
Definition: free_surface.cpp:166
FreeSurfaceOps::OpLhsH_dG
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1118
l
FTensor::Index< 'l', SPACE_DIM > l
Definition: free_surface.cpp:135
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:144
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:501
OpMixScalarTimesDiv
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
Definition: free_surface.cpp:125
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
MoFEM::CoreInterface::get_dofs
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
MoFEM::MeshsetsManager::getMeshset
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:705
h
double h
Definition: free_surface.cpp:169
FreeSurfaceOps::OpRhsH
Definition: FreeSurfaceOps.hpp:796
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
cylindrical
auto cylindrical
Definition: free_surface.cpp:187
OpDomainAssembleVector
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
Definition: free_surface.cpp:117
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:442
U_FIELD_DIM
constexpr int U_FIELD_DIM
Definition: free_surface.cpp:79
OpBoundaryMassL
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
Definition: free_surface.cpp:114
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
FreeSurfaceOps::OpLhsH_dU
Lhs for H dU.
Definition: FreeSurfaceOps.hpp:916
coord_type
int coord_type
Definition: free_surface.cpp:75
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1018
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:137
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
save_range
auto save_range
Definition: free_surface.cpp:324
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
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
rho_diff
double rho_diff
Definition: free_surface.cpp:179
get_current_bit
auto get_current_bit
dofs bit used to do calculations
Definition: free_surface.cpp:149
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:219
double
convert.type
type
Definition: convert.py:64
OpDomainMassH
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
Definition: free_surface.cpp:110
get_global_size
auto get_global_size
Definition: free_surface.cpp:318
FreeSurfaceOps::OpLhsH_dH
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:974
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
nb_levels
int nb_levels
Definition: free_surface.cpp:141
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
eta2
double eta2
Definition: free_surface.cpp:171
BoundaryParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
Definition: free_surface.cpp:92
OpDomainMassP
OpDomainMassH OpDomainMassP
Definition: free_surface.cpp:111
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
get_projection_bit
auto get_projection_bit
Definition: free_surface.cpp:154
eta
double eta
Definition: free_surface.cpp:170
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
bit
auto bit
Definition: free_surface.cpp:312
FreeSurfaceOps::OpRhsG
Definition: FreeSurfaceOps.hpp:1182
SPACE_DIM
constexpr int SPACE_DIM
Definition: free_surface.cpp:78
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:118
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:701
FreeSurfaceOps::OpLhsU_dG
Lhs for G dH.
Definition: FreeSurfaceOps.hpp:733
get_skin_projection_bit
auto get_skin_projection_bit
Definition: free_surface.cpp:155
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
order
int order
approximation order
Definition: free_surface.cpp:140
mu_diff
double mu_diff
Definition: free_surface.cpp:181
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
convert.n
n
Definition: convert.py:82
rho_p
double rho_p
Definition: free_surface.cpp:163
FreeSurfaceOps::OpLhsG_dG
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1322
kappa
double kappa
Definition: free_surface.cpp:183
rho_m
double rho_m
Definition: free_surface.cpp:161
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
tsPrePostProc
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
Definition: free_surface.cpp:467
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:368
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:217
CYLINDRICAL
@ CYLINDRICAL
Definition: definitions.h:117
lambda
double lambda
surface tension
Definition: free_surface.cpp:165
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
FreeSurface::setParentDofs
MoFEMErrorCode setParentDofs(boost::shared_ptr< FEMethod > fe_top, std::string field_name, ForcesAndSourcesCore::UserDataOperator::OpType op, BitRefLevel child_ent_bit, boost::function< boost::shared_ptr< ForcesAndSourcesCore >()> get_elem, int verbosity, LogManager::SeverityLevel sev)
Create hierarchy of elements run on parents levels.
Definition: free_surface.cpp:2811
FreeSurfaceOps::OpLhsU_dU
Lhs for U dU.
Definition: FreeSurfaceOps.hpp:473
mu_ave
double mu_ave
Definition: free_surface.cpp:180
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:108
refine_overlap
int refine_overlap
Definition: free_surface.cpp:142
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:345
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:1060
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
FreeSurface::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: free_surface.cpp:2045
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
get_start_bit
auto get_start_bit
Definition: free_surface.cpp:146
FreeSurface::readMesh
MoFEMErrorCode readMesh()
[Run programme]
Definition: free_surface.cpp:534
LAST_COORDINATE_SYSTEM
@ LAST_COORDINATE_SYSTEM
Definition: definitions.h:119
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
BoundaryEleOp
BoundaryEle::UserDataOperator BoundaryEleOp
Definition: free_surface.cpp:93
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:232
FreeSurfaceOps::OpLhsU_dH
Lhs for U dH.
Definition: FreeSurfaceOps.hpp:611
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
a0
double a0
Definition: free_surface.cpp:160
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:1046
FreeSurface::mField
MoFEM::Interface & mField
Definition: free_surface.cpp:477
QUIET
@ QUIET
Definition: definitions.h:208
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:18
mu_p
double mu_p
Definition: free_surface.cpp:164
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
FreeSurfaceOps::OpLhsG_dH
Lhs for H dH.
Definition: FreeSurfaceOps.hpp:1250
MoFEM::SmartPetscObj< DM >
FreeSurface::refineMesh
MoFEMErrorCode refineMesh(size_t overlap)
Definition: free_surface.cpp:2578
FreeSurfaceOps::OpNormalConstrainLhs
Definition: FreeSurfaceOps.hpp:195
FreeSurface::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: free_surface.cpp:550
OpDomainAssembleScalar
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
Definition: free_surface.cpp:119
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:590
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
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:89
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1289
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
OpDomainMassG
OpDomainMassH OpDomainMassG
Definition: free_surface.cpp:112
FreeSurface::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: free_surface.cpp:662
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
Generated by Doxygen 1.8.17 and hosted at University of Glasgow