C++ Facet::entities方法代码示例

//-----------------------------------------------------------------------------
double HexahedronCell::facet_area(const Cell& cell, std::size_t facet) const
{
  // Create facet from the mesh and local facet number
  const Facet f(cell.mesh(), cell.entities(1)[facet]);

  // Get global index of vertices on the facet
  const std::size_t v0 = f.entities(0)[0];
  const std::size_t v1 = f.entities(0)[1];
  const std::size_t v2 = f.entities(0)[2];
  const std::size_t v3 = f.entities(0)[3];

  // Get mesh geometry
  const MeshGeometry& geometry = cell.mesh().geometry();

  // Need to check points are co-planar

  const Point p0 = geometry.point(v0);
  const Point p1 = geometry.point(v1);
  const Point p2 = geometry.point(v2);
  const Point p3 = geometry.point(v3);

  dolfin_not_implemented();

  return 0.0;
}

//-----------------------------------------------------------------------------
bool DirichletBC::on_facet(const double* coordinates, const Facet& facet) const
{
  // Check if the coordinates are on the same line as the line segment
  if (facet.dim() == 1)
  {
    // Create points
    Point p(coordinates[0], coordinates[1]);
    const Point v0 = Vertex(facet.mesh(), facet.entities(0)[0]).point();
    const Point v1 = Vertex(facet.mesh(), facet.entities(0)[1]).point();

    // Create vectors
    const Point v01 = v1 - v0;
    const Point vp0 = v0 - p;
    const Point vp1 = v1 - p;

    // Check if the length of the sum of the two line segments vp0 and
    // vp1 is equal to the total length of the facet
    if ( std::abs(v01.norm() - vp0.norm() - vp1.norm()) < DOLFIN_EPS )
      return true;
    else
      return false;
  }
  // Check if the coordinates are in the same plane as the triangular
  // facet
  else if (facet.dim() == 2)
  {
    // Create points
    const Point p(coordinates[0], coordinates[1], coordinates[2]);
    const Point v0 = Vertex(facet.mesh(), facet.entities(0)[0]).point();
    const Point v1 = Vertex(facet.mesh(), facet.entities(0)[1]).point();
    const Point v2 = Vertex(facet.mesh(), facet.entities(0)[2]).point();

    // Create vectors
    const Point v01 = v1 - v0;
    const Point v02 = v2 - v0;
    const Point vp0 = v0 - p;
    const Point vp1 = v1 - p;
    const Point vp2 = v2 - p;

    // Check if the sum of the area of the sub triangles is equal to
    // the total area of the facet
    if (std::abs(v01.cross(v02).norm() - vp0.cross(vp1).norm()
                 - vp1.cross(vp2).norm() - vp2.cross(vp0).norm()) < DOLFIN_EPS)
    {
      return true;
    }
    else
      return false;
  }

  dolfin_error("DirichletBC.cpp",
               "determine if given point is on facet",
               "Not implemented for given facet dimension");

  return false;
}

//-----------------------------------------------------------------------------
bool NewDirichletBC::onFacet(real* coordinates, Facet& facet)
{
  // Check if the coordinates are on the same line as the line segment
  if ( facet.dim() == 1 )
  {
    // Create points
    Point p(coordinates[0], coordinates[1]);
    Point v0 = Vertex(facet.mesh(), facet.entities(0)[0]).point();
    Point v1 = Vertex(facet.mesh(), facet.entities(0)[1]).point();

    // Create vectors
    Point v01 = v1 - v0;
    Point vp0 = v0 - p;
    Point vp1 = v1 - p;

    // Check if the length of the sum of the two line segments vp0 and vp1 is
    // equal to the total length of the facet
    if ( std::abs(v01.norm() - vp0.norm() - vp1.norm()) < DOLFIN_EPS )
      return true;
    else
      return false;
  }
  // Check if the coordinates are in the same plane as the triangular facet
  else if ( facet.dim() == 2 )
  {
    // Create points
    Point p(coordinates[0], coordinates[1], coordinates[2]);
    Point v0 = Vertex(facet.mesh(), facet.entities(0)[0]).point();
    Point v1 = Vertex(facet.mesh(), facet.entities(0)[1]).point();
    Point v2 = Vertex(facet.mesh(), facet.entities(0)[2]).point();

    // Create vectors
    Point v01 = v1 - v0;
    Point v02 = v2 - v0;
    Point vp0 = v0 - p;
    Point vp1 = v1 - p;
    Point vp2 = v2 - p;

    // Check if the sum of the area of the sub triangles is equal to the total
    // area of the facet
    if ( std::abs(v01.cross(v02).norm() - vp0.cross(vp1).norm() - vp1.cross(vp2).norm()
        - vp2.cross(vp0).norm()) < DOLFIN_EPS )
      return true;
    else
      return false;
  }

  error("Unable to determine if given point is on facet (not implemented for given facet dimension).");

  return false;
}

//-----------------------------------------------------------------------------
double TriangleCell::facet_area(const Cell& cell, std::size_t facet) const
{
  // Create facet from the mesh and local facet number
  const Facet f(cell.mesh(), cell.entities(1)[facet]);

  // Get global index of vertices on the facet
  const std::size_t v0 = f.entities(0)[0];
  const std::size_t v1 = f.entities(0)[1];

  // Get mesh geometry
  const MeshGeometry& geometry = cell.mesh().geometry();

  // Get the coordinates of the two vertices
  const double* p0 = geometry.x(v0);
  const double* p1 = geometry.x(v1);

  // Compute distance between vertices
  double d = 0.0;
  for (std::size_t i = 0; i < geometry.dim(); i++)
  {
    const double dp = p0[i] - p1[i];
    d += dp*dp;
  }

  return std::sqrt(d);
}

//-----------------------------------------------------------------------------
void NewDirichletBC::computeBCTopological(std::map<uint, real>& boundary_values,
                                       Facet& facet,
                                       BoundaryCondition::LocalData& data)
{
  // Get cell to which facet belongs (there may be two, but pick first)
  Cell cell(_mesh, facet.entities(facet.dim() + 1)[0]);
  UFCCell ufc_cell(cell);
  
  // Get local index of facet with respect to the cell
  const uint local_facet = cell.index(facet);
  
  // Interpolate function on cell
  g.interpolate(data.w, ufc_cell, *data.finite_element, cell, local_facet);
  
  // Tabulate dofs on cell
  data.dof_map->tabulate_dofs(data.cell_dofs, ufc_cell);

  // Tabulate which dofs are on the facet
  data.dof_map->tabulate_facet_dofs(data.facet_dofs, local_facet);
  
  // Pick values for facet
  for (uint i = 0; i < data.dof_map->num_facet_dofs(); i++)
  {
    const uint dof = data.offset + data.cell_dofs[data.facet_dofs[i]];
    const real value = data.w[data.facet_dofs[i]];
    boundary_values[dof] = value;
  }
}

//-----------------------------------------------------------------------------
double QuadrilateralCell::facet_area(const Cell& cell, std::size_t facet) const
{
  // Create facet from the mesh and local facet number
  const Facet f(cell.mesh(), cell.entities(1)[facet]);

  // Get global index of vertices on the facet
  const std::size_t v0 = f.entities(0)[0];
  const std::size_t v1 = f.entities(0)[1];

  // Get mesh geometry
  const MeshGeometry& geometry = cell.mesh().geometry();

  const Point p0 = geometry.point(v0);
  const Point p1 = geometry.point(v1);

  return (p0 - p1).norm();
}

//-----------------------------------------------------------------------------
double TriangleCell::facet_area(const Cell& cell, std::size_t facet) const
{
  // Create facet from the mesh and local facet number
  const Facet f(cell.mesh(), cell.entities(1)[facet]);

  // Get global index of vertices on the facet
  const std::size_t v0 = f.entities(0)[0];
  const std::size_t v1 = f.entities(0)[1];

  // Get mesh geometry
  const MeshGeometry& geometry = cell.mesh().geometry();

  // Get the coordinates of the two vertices
  const Point p0 = geometry.point(v0);
  const Point p1 = geometry.point(v1);

  return p1.distance(p0);
}

//-----------------------------------------------------------------------------
void BoundaryComputation::reorder(std::vector<std::size_t>& vertices,
                                  const Facet& facet)
{
  // Get mesh
  const Mesh& mesh = facet.mesh();

  // Get the vertex opposite to the facet (the one we remove)
  std::size_t vertex = 0;
  const Cell cell(mesh, facet.entities(mesh.topology().dim())[0]);
  for (std::size_t i = 0; i < cell.num_entities(0); i++)
  {
    bool not_in_facet = true;
    vertex = cell.entities(0)[i];
    for (std::size_t j = 0; j < facet.num_entities(0); j++)
    {
      if (vertex == facet.entities(0)[j])
      {
        not_in_facet = false;
        break;
      }
    }
    if (not_in_facet)
      break;
  }
  const Point p = mesh.geometry().point(vertex);

  // Check orientation
  switch (mesh.type().cell_type())
  {
  case CellType::interval:
    // Do nothing
    break;
  case CellType::triangle:
    {
      dolfin_assert(facet.num_entities(0) == 2);

      const Point p0 = mesh.geometry().point(facet.entities(0)[0]);
      const Point p1 = mesh.geometry().point(facet.entities(0)[1]);
      const Point v = p1 - p0;
      const Point n(v.y(), -v.x());

      if (n.dot(p0 - p) < 0.0)
      {
        const std::size_t tmp = vertices[0];
        vertices[0] = vertices[1];
        vertices[1] = tmp;
      }
    }
    break;
  case CellType::tetrahedron:
    {
      dolfin_assert(facet.num_entities(0) == 3);

      const Point p0 = mesh.geometry().point(facet.entities(0)[0]);
      const Point p1 = mesh.geometry().point(facet.entities(0)[1]);
      const Point p2 = mesh.geometry().point(facet.entities(0)[2]);
      const Point v1 = p1 - p0;
      const Point v2 = p2 - p0;
      const Point n  = v1.cross(v2);

      if (n.dot(p0 - p) < 0.0)
      {
        const std::size_t tmp = vertices[0];
        vertices[0] = vertices[1];
        vertices[1] = tmp;
      }
    }
    break;
  default:
    {
      dolfin_error("BoundaryComputation.cpp",
                   "reorder cell for extraction of mesh boundary",
                   "Unknown cell type (%d)",
                   mesh.type().cell_type());
    }
  }
}