C++ MeshEditor::add_entity_point方法代码示例

//-----------------------------------------------------------------------------
void dolfin::p_refine(Mesh& refined_mesh, const Mesh& mesh)
{
  MeshEditor editor;
  if (mesh.geometry().degree() != 1)
  {
    dolfin_error("refine.cpp",
                 "increase polynomial degree of mesh",
                 "Currently only linear -> quadratic is supported");
  }

  const CellType::Type cell_type = mesh.type().cell_type();

  if (cell_type != CellType::Type::triangle
      and cell_type != CellType::Type::tetrahedron
      and cell_type != CellType::Type::interval)
  {
    dolfin_error("refine.cpp",
                 "increase polynomial degree of mesh",
                 "Unsupported cell type");
  }

  const std::size_t tdim = mesh.topology().dim();
  const std::size_t gdim = mesh.geometry().dim();

  editor.open(refined_mesh, cell_type, tdim, gdim, 2);

  // Copy over mesh
  editor.init_vertices_global(mesh.num_entities(0), mesh.num_entities_global(0));
  for (VertexIterator v(mesh); !v.end(); ++v)
    editor.add_vertex(v->index(), v->point());

  editor.init_cells_global(mesh.num_entities(tdim), mesh.num_entities_global(tdim));
  std::vector<std::size_t> verts(tdim + 1);
  for (CellIterator c(mesh); !c.end(); ++c)
  {
    std::copy(c->entities(0), c->entities(0) + tdim + 1, verts.begin());
    editor.add_cell(c->index(), verts);
  }

  // Initialise edges
  editor.init_entities();

  // Add points at centres of edges
  for (EdgeIterator e(refined_mesh); !e.end(); ++e)
    editor.add_entity_point(1, 0, e->index(), e->midpoint());

  editor.close();
}

//-----------------------------------------------------------------------------
UnitDiscMesh::UnitDiscMesh(MPI_Comm comm, std::size_t n,
                           std::size_t degree, std::size_t gdim)
  : Mesh(comm)
{
  dolfin_assert(n > 0);
  dolfin_assert(gdim == 2 or gdim == 3);
  dolfin_assert(degree == 1 or degree == 2);

  MeshEditor editor;
  editor.open(*this, 2, gdim, degree);
  editor.init_vertices_global(1 + 3*n*(n + 1),
                              1 + 3*n*(n + 1));

  std::size_t c = 0;
  editor.add_vertex(c, Point(0,0,0));
  ++c;

  for (std::size_t i = 1; i <= n; ++i)
    for (std::size_t j = 0; j < 6*i; ++j)
    {
      double r = (double)i/(double)n;
      double th = 2*M_PI*(double)j/(double)(6*i);
      double x = r*cos(th);
      double y = r*sin(th);
      editor.add_vertex(c, Point(x, y, 0));
      ++c;
    }

  editor.init_cells(6*n*n);

  c = 0;
  std::size_t base_i = 0;
  std::size_t row_i = 1;
  for (std::size_t i = 1; i <= n; ++i)
  {
    std::size_t base_m = base_i;
    base_i = 1 + 3*i*(i - 1);
    std::size_t row_m = row_i;
    row_i = 6*i;

    for (std::size_t k = 0; k != 6; ++k)
      for (std::size_t j = 0; j < (i*2 - 1); ++j)
      {
        std::size_t i0, i1, i2;
        if (j%2 == 0)
        {
          i0 = base_i + (k*i + j/2)%row_i;
          i1 = base_i + (k*i + j/2 + 1)%row_i;
          i2 = base_m + (k*(i-1) + j/2)%row_m;
        }
        else
        {
          i0 = base_m + (k*(i-1) + j/2)%row_m;
          i1 = base_m + (k*(i-1) + j/2 + 1)%row_m;
          i2 = base_i + (k*i + j/2 + 1)%row_i;
        }

        editor.add_cell(c, i0, i1, i2);
        ++c;
      }
  }

  // Initialise entities required for this degree polynomial mesh
  // and allocate space for the point coordinate data

  if (degree == 2)
  {
    editor.init_entities();

    for (EdgeIterator e(*this); !e.end(); ++e)
    {
      Point v0 = Vertex(*this, e->entities(0)[0]).point();
      Point v1 = Vertex(*this, e->entities(0)[1]).point();
      Point pt = e->midpoint();

      if (std::abs(v0.norm() - 1.0) < 1e-6 and
          std::abs(v1.norm() - 1.0) < 1e-6)
        pt *= v0.norm()/pt.norm();

      // Add Edge-based point
      editor.add_entity_point(1, 0, e->index(), pt);
    }
  }

  editor.close();
}