C++ MeshEntity::dim方法代码示例

//-----------------------------------------------------------------------------
double QuadrilateralCell::volume(const MeshEntity& cell) const
{
  if (cell.dim() != 2)
  {
    dolfin_error("QuadrilateralCell.cpp",
                 "compute volume (area) of cell",
                 "Illegal mesh entity");
  }

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

  // Get the coordinates of the four vertices
  const unsigned int* vertices = cell.entities(0);
  const Point p0 = geometry.point(vertices[0]);
  const Point p1 = geometry.point(vertices[1]);
  const Point p2 = geometry.point(vertices[2]);
  const Point p3 = geometry.point(vertices[3]);

  if (geometry.dim() == 2)
  {
    const Point c = (p0 - p2).cross(p1 - p3);
    return 0.5 * c.norm();
  }
  else
    dolfin_error("QuadrilateralCell.cpp",
                 "compute volume of quadrilateral",
                 "Only know how to compute volume in R^2");

  // FIXME: could work in R^3 but need to check co-planarity

  return 0.0;
}

//-----------------------------------------------------------------------------
double HexahedronCell::volume(const MeshEntity& cell) const
{
  if (cell.dim() != 2)
  {
    dolfin_error("HexahedronCell.cpp",
                 "compute volume (area) of cell",
                 "Illegal mesh entity");
  }

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

  // Get the coordinates of the four vertices
  const unsigned int* vertices = cell.entities(0);
  const Point p0 = geometry.point(vertices[0]);
  const Point p1 = geometry.point(vertices[1]);
  const Point p2 = geometry.point(vertices[2]);
  const Point p3 = geometry.point(vertices[3]);
  const Point p4 = geometry.point(vertices[4]);
  const Point p5 = geometry.point(vertices[5]);

  dolfin_error("HexahedronCell.cpp",
               "compute volume of hexahedron",
               "Not implemented");

  return 0.0;
}

//-----------------------------------------------------------------------------
double IntervalCell::volume(const MeshEntity& interval) const
{
  // Check that we get an interval
  if (interval.dim() != 1)
  {
    dolfin_error("IntervalCell.cpp",
                 "compute volume (length) of interval cell",
                 "Illegal mesh entity, not an interval");
  }

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

  // Get the coordinates of the two vertices
  const unsigned int* vertices = interval.entities(0);
  const double* x0 = geometry.x(vertices[0]);
  const double* x1 = geometry.x(vertices[1]);

  // Compute length of interval (line segment)
  double sum = 0.0;
  for (std::size_t i = 0; i < geometry.dim(); ++i)
  {
    const double dx = x1[i] - x0[i];
    sum += dx*dx;
  }

  return std::sqrt(sum);
}

//-----------------------------------------------------------------------------
double IntervalCell::diameter(const MeshEntity& interval) const
{
  // Check that we get an interval
  if (interval.dim() != 1)
  {
    dolfin_error("IntervalCell.cpp",
                 "compute diameter of interval cell",
                 "Illegal mesh entity, not an interval");
  }

  // Diameter is same as volume for interval (line segment)
  return volume(interval);
}

//-----------------------------------------------------------------------------
double IntervalCell::circumradius(const MeshEntity& interval) const
{
  // Check that we get an interval
  if (interval.dim() != 1)
  {
    dolfin_error("IntervalCell.cpp",
                 "compute diameter of interval cell",
                 "Illegal mesh entity, not an interval");
  }

  // Circumradius is half the volume for an interval (line segment)
  return volume(interval)/2.0;
}

//-----------------------------------------------------------------------------
double TriangleCell::volume(const MeshEntity& triangle) const
{
  // Check that we get a triangle
  if (triangle.dim() != 2)
  {
    dolfin_error("TriangleCell.cpp",
                 "compute volume (area) of triangle cell",
                 "Illegal mesh entity, not a triangle");
  }

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

  // Get the coordinates of the three vertices
  const unsigned int* vertices = triangle.entities(0);
  const double* x0 = geometry.x(vertices[0]);
  const double* x1 = geometry.x(vertices[1]);
  const double* x2 = geometry.x(vertices[2]);

  if (geometry.dim() == 2)
  {
    // Compute area of triangle embedded in R^2
    double v2 = (x0[0]*x1[1] + x0[1]*x2[0] + x1[0]*x2[1])
      - (x2[0]*x1[1] + x2[1]*x0[0] + x1[0]*x0[1]);

    // Formula for volume from http://mathworld.wolfram.com
    return 0.5 * std::abs(v2);
  }
  else if (geometry.dim() == 3)
  {
    // Compute area of triangle embedded in R^3
    const double v0 = (x0[1]*x1[2] + x0[2]*x2[1] + x1[1]*x2[2])
      - (x2[1]*x1[2] + x2[2]*x0[1] + x1[1]*x0[2]);
    const double v1 = (x0[2]*x1[0] + x0[0]*x2[2] + x1[2]*x2[0])
      - (x2[2]*x1[0] + x2[0]*x0[2] + x1[2]*x0[0]);
    const double v2 = (x0[0]*x1[1] + x0[1]*x2[0] + x1[0]*x2[1])
      - (x2[0]*x1[1] + x2[1]*x0[0] + x1[0]*x0[1]);

    // Formula for volume from http://mathworld.wolfram.com
    return  0.5*sqrt(v0*v0 + v1*v1 + v2*v2);
  }
  else
  {
    dolfin_error("TriangleCell.cpp",
                 "compute volume of triangle",
                 "Only know how to compute volume when embedded in R^2 or R^3");
  }

  return 0.0;
}

//-----------------------------------------------------------------------------
double HexahedronCell::diameter(const MeshEntity& cell) const
{
  // Check that we get a cell
  if (cell.dim() != 2)
  {
    dolfin_error("HexahedronCell.cpp",
                 "compute diameter of hexahedron cell",
                 "Illegal mesh entity");
  }

  dolfin_error("HexahedronCell.cpp",
               "compute diameter of hexahedron cell",
               "Don't know how to compute diameter");

  dolfin_not_implemented();
  return 0.0;
}

//-----------------------------------------------------------------------------
double IntervalCell::volume(const MeshEntity& interval) const
{
  // Check that we get an interval
  if (interval.dim() != 1)
  {
    dolfin_error("IntervalCell.cpp",
                 "compute volume (length) of interval cell",
                 "Illegal mesh entity, not an interval");
  }

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

  // Get the coordinates of the two vertices
  const unsigned int* vertices = interval.entities(0);
  const Point x0 = geometry.point(vertices[0]);
  const Point x1 = geometry.point(vertices[1]);

  return x1.distance(x0);
}

//-----------------------------------------------------------------------------
bool CollisionDetection::collides(const MeshEntity& entity,
                                  const Point& point)
{
  switch (entity.dim())
  {
  case 0:
    dolfin_not_implemented();
    break;
  case 1:
    return collides_interval_point(entity, point);
  case 2:
    return collides_triangle_point(entity, point);
  case 3:
    return collides_tetrahedron_point(entity, point);
  default:
    dolfin_error("CollisionDetection.cpp",
		 "collides entity with point",
		 "Unknown dimension of entity");
  }

  return false;
}

//-----------------------------------------------------------------------------
double TriangleCell::diameter(const MeshEntity& triangle) const
{
  // Check that we get a triangle
  if (triangle.dim() != 2)
  {
    dolfin_error("TriangleCell.cpp",
                 "compute diameter of triangle cell",
                 "Illegal mesh entity, not a triangle");
  }

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

  // Only know how to compute the diameter when embedded in R^2 or R^3
  if (geometry.dim() != 2 && geometry.dim() != 3)
    dolfin_error("TriangleCell.cpp",
                 "compute diameter of triangle",
                 "Only know how to compute diameter when embedded in R^2 or R^3");

  // Get the coordinates of the three vertices
  const unsigned int* vertices = triangle.entities(0);
  const Point p0 = geometry.point(vertices[0]);
  const Point p1 = geometry.point(vertices[1]);
  const Point p2 = geometry.point(vertices[2]);

  // FIXME: Assuming 3D coordinates, could be more efficient if
  // FIXME: if we assumed 2D coordinates in 2D

  // Compute side lengths
  const double a  = p1.distance(p2);
  const double b  = p0.distance(p2);
  const double c  = p0.distance(p1);

  // Formula for diameter (2*circumradius) from http://mathworld.wolfram.com
  return 0.5*a*b*c / volume(triangle);
}

//-----------------------------------------------------------------------------
bool CollisionPredicates::collides(const MeshEntity& entity_0,
				   const MeshEntity& entity_1)
{
  // Intersection is only implemented for simplex meshes
  if (!entity_0.mesh().type().is_simplex() ||
      !entity_1.mesh().type().is_simplex())
  {
    dolfin_error("Cell.cpp",
		 "intersect cell and point",
		 "intersection is only implemented for simplex meshes");
  }

  // Get data
  const MeshGeometry& g0 = entity_0.mesh().geometry();
  const MeshGeometry& g1 = entity_1.mesh().geometry();
  const unsigned int* v0 = entity_0.entities(0);
  const unsigned int* v1 = entity_1.entities(0);
  const std::size_t d0 = entity_0.dim();
  const std::size_t d1 = entity_1.dim();
  const std::size_t gdim = g0.dim();
  dolfin_assert(gdim == g1.dim());

  // Pick correct specialized implementation
  if (d0 == 1 && d1 == 1)
  {
    return collides_segment_segment(g0.point(v0[0]),
                                    g0.point(v0[1]),
                                    g1.point(v1[0]),
                                    g1.point(v1[1]),
                                    gdim);
  }

  if (d0 == 1 && d1 == 2)
  {
    return collides_triangle_segment(g1.point(v1[0]),
				     g1.point(v1[1]),
				     g1.point(v1[2]),
				     g0.point(v0[0]),
				     g0.point(v0[1]),
				     gdim);
  }

  if (d0 == 2 && d1 == 1)
  {
    return collides_triangle_segment(g0.point(v0[0]),
                                     g0.point(v0[1]),
                                     g0.point(v0[2]),
                                     g1.point(v1[0]),
                                     g1.point(v1[1]),
                                     gdim);
  }

  if (d0 == 2 && d1 == 2)
  {
    return collides_triangle_triangle(g0.point(v0[0]),
                                      g0.point(v0[1]),
                                      g0.point(v0[2]),
                                      g1.point(v1[0]),
                                      g1.point(v1[1]),
                                      g1.point(v1[2]),
                                      gdim);
  }

  if (d0 == 2 && d1 == 3)
  {
    return collides_tetrahedron_triangle_3d(g1.point(v1[0]),
                                            g1.point(v1[1]),
                                            g1.point(v1[2]),
                                            g1.point(v1[3]),
                                            g0.point(v0[0]),
                                            g0.point(v0[1]),
                                            g0.point(v0[2]));
  }

  if (d0 == 3 && d1 == 2)
  {
    return collides_tetrahedron_triangle_3d(g0.point(v0[0]),
                                            g0.point(v0[1]),
                                            g0.point(v0[2]),
                                            g0.point(v0[3]),
                                            g1.point(v1[0]),
                                            g1.point(v1[1]),
                                            g1.point(v1[2]));
  }

  if (d0 == 3 && d1 == 3)
  {
    return collides_tetrahedron_tetrahedron_3d(g0.point(v0[0]),
                                               g0.point(v0[1]),
                                               g0.point(v0[2]),
                                               g0.point(v0[3]),
                                               g1.point(v1[0]),
                                               g1.point(v1[1]),
                                               g1.point(v1[2]),
                                               g1.point(v1[3]));
  }

  dolfin_error("CollisionPredicates.cpp",
               "compute entity-entity collision",
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
bool
CollisionDetection::collides(const MeshEntity& entity_0,
			     const MeshEntity& entity_1)
{
  switch (entity_0.dim())
  {
  case 0:
    // Collision with PointCell
    dolfin_not_implemented();
    break;
  case 1:
    // Collision with interval
    switch (entity_1.dim())
    {
    case 0:
      dolfin_not_implemented();
      break;
    case 1:
      return collides_interval_interval(entity_1, entity_0);
      break;
    case 2:
      dolfin_not_implemented();
      break;
    case 3:
      dolfin_not_implemented();
      break;
    default:
      dolfin_error("CollisionDetection.cpp",
                   "collides entity_0 with entity_1",
                   "Unknown dimension of entity_1 in IntervalCell collision");
    }
    break;
  case 2:
    // Collision with triangle
    switch (entity_1.dim())
    {
    case 0:
      dolfin_not_implemented();
      break;
    case 1:
      dolfin_not_implemented();
      break;
    case 2:
      return collides_triangle_triangle(entity_0, entity_1);
    case 3:
      return collides_tetrahedron_triangle(entity_1, entity_0);
    default:
      dolfin_error("CollisionDetection.cpp",
		   "collides entity_0 with entity_1",
		   "Unknown dimension of entity_1 in TriangleCell collision");
    }
    break;
  case 3:
    // Collision with tetrahedron
    switch (entity_1.dim())
    {
    case 0:
     dolfin_not_implemented();
      break;
    case 1:
      dolfin_not_implemented();
      break;
    case 2:
      return collides_tetrahedron_triangle(entity_0, entity_1);
      break;
    case 3:
      return collides_tetrahedron_tetrahedron(entity_0, entity_1);
      break;
    default:
      dolfin_error("CollisionDetection.cpp",
		   "collides entity_0 with entity_1",
		   "Unknown dimension of entity_1 in TetrahedronCell collision");
    }
    break;
  default:
    dolfin_error("CollisionDetection.cpp",
		 "collides entity_0 with entity_1",
		 "Unknown dimension of entity_0");
  }

  return false;
}