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

Mesh build_mesh(const std::vector<unsigned int>& cells, const std::vector<double>& vertices, int dim)
{
    // vertices and cells are flattened
    unsigned int vlen = vertices.size() / dim;
    unsigned int clen = cells.size() / (dim + 1);

    Mesh mesh;
    
    MeshEditor editor;
    editor.open(mesh, dim, dim);
    editor.init_vertices(vlen);
    editor.init_cells(clen);
    if (dim==3)
    {
        for (int i=0; i<vlen; i++)
	    editor.add_vertex(i, vertices[3*i], vertices[3*i+1], vertices[3*i+2]);
    	for (int i=0; i<clen; i++)
	    editor.add_cell(i, cells[4*i], cells[4*i+1], cells[4*i+2], cells[4*i+3]);
    }
    else
    {
        for (int i=0; i<vlen; i++)
	    editor.add_vertex(i, vertices[2*i], vertices[2*i+1]);
    	for (int i=0; i<clen; i++)
	    editor.add_cell(i, cells[3*i], cells[3*i+1], cells[3*i+2]);
    }
    editor.close();

    return mesh;
}

//-----------------------------------------------------------------------------
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();
}

//-----------------------------------------------------------------------------
void DynamicMeshEditor::close(bool order)
{
  dolfin_assert(_mesh);
  dolfin_assert(_cell_type);

  // Open default mesh editor
  MeshEditor editor;
  editor.open(*_mesh, _cell_type->cell_type(), _tdim, _gdim);

  // Set number of vertices
  const std::size_t num_vertices = vertex_coordinates.size()/_gdim;
  editor.init_vertices(num_vertices);

  // Set number of cells
  const std::size_t vertices_per_cell = _cell_type->num_vertices(_gdim);
  const std::size_t num_cells = cell_vertices.size()/vertices_per_cell;
  editor.init_cells(num_cells);

  // Add vertices
  std::vector<double> p(_gdim);
  for (std::size_t v = 0; v < num_vertices; v++)
  {
    const std::size_t offset = v*_gdim;
    for (std::size_t i = 0; i < _gdim; i++)
      p[i] = vertex_coordinates[offset + i];
    editor.add_vertex(v, p);
  }

  // Add cells
  std::vector<std::size_t> vertices(vertices_per_cell);
  for (std::size_t c = 0; c < num_cells; c++)
  {
    const std::size_t offset = c*vertices_per_cell;
    for (std::size_t i = 0; i < vertices_per_cell; i++)
      vertices[i] = cell_vertices[offset + i];
    editor.add_cell(c, vertices);
  }

  // Close editor
  editor.close(order);

  // Clear data
  clear();
}

//-----------------------------------------------------------------------------
UnitTriangleMesh::UnitTriangleMesh() : Mesh()
{
  // Receive mesh according to parallel policy
  if (MPI::is_receiver(this->mpi_comm()))
  {
    MeshPartitioning::build_distributed_mesh(*this);
    return;
  }

  // Open mesh for editing
  MeshEditor editor;
  editor.open(*this, CellType::triangle, 2, 2);

  // Create vertices
  editor.init_vertices_global(3, 3);
  std::vector<double> x(2);
  x[0] = 0.0; x[1] = 0.0;
  editor.add_vertex(0, x);
  x[0] = 1.0; x[1] = 0.0;
  editor.add_vertex(1, x);
  x[0] = 0.0; x[1] = 1.0;
  editor.add_vertex(2, x);

  // Create cells
  editor.init_cells_global(1, 1);
  std::vector<std::size_t> cell_data(3);
  cell_data[0] = 0; cell_data[1] = 1; cell_data[2] = 2;
  editor.add_cell(0, cell_data);

  // Close mesh editor
  editor.close();

  // Broadcast mesh according to parallel policy
  if (MPI::is_broadcaster(this->mpi_comm()))
  {
    MeshPartitioning::build_distributed_mesh(*this);
    return;
  }
}

//-----------------------------------------------------------------------------
void ParallelRefinement::build_local(Mesh& new_mesh) const
{
  MeshEditor ed;
  const std::size_t tdim = _mesh.topology().dim();
  const std::size_t gdim = _mesh.geometry().dim();
  dolfin_assert(new_vertex_coordinates.size()%gdim == 0);
  const std::size_t num_vertices = new_vertex_coordinates.size()/gdim;

  const std::size_t num_cell_vertices = tdim + 1;
  dolfin_assert(new_cell_topology.size()%num_cell_vertices == 0);
  const std::size_t num_cells = new_cell_topology.size()/num_cell_vertices;

  ed.open(new_mesh, tdim, gdim);
  ed.init_vertices(num_vertices);
  std::size_t i = 0;
  for (auto p = new_vertex_coordinates.begin();
       p != new_vertex_coordinates.end(); p += gdim)
  {
    std::vector<double> vertex(p, p + gdim);
    ed.add_vertex(i, vertex);
    ++i;
  }

  ed.init_cells(num_cells);
  i = 0;
  std::vector<std::size_t> cell(num_cell_vertices);
  for (auto p = new_cell_topology.begin(); p != new_cell_topology.end();
       p += num_cell_vertices)
  {
    std::copy(p, p + num_cell_vertices, cell.begin());
    ed.add_cell(i, cell);
    ++i;
  }
  ed.close();

}

//-----------------------------------------------------------------------------
void IntervalMesh::build(std::size_t nx, double a, double b)
{
  // Receive mesh according to parallel policy
  if (MPI::is_receiver(this->mpi_comm()))
  {
    MeshPartitioning::build_distributed_mesh(*this);
    return;
  }

  if (std::abs(a - b) < DOLFIN_EPS)
  {
    dolfin_error("Interval.cpp",
                 "create interval",
                 "Length of interval is zero. Consider checking your dimensions");
  }

  if (b < a)
  {
    dolfin_error("Interval.cpp",
                 "create interval",
                 "Length of interval is negative. Consider checking the order of your arguments");
  }

  if (nx < 1)
  {
    dolfin_error("Interval.cpp",
                 "create interval",
                 "Number of points on interval is (%d), it must be at least 1", nx);
  }

  rename("mesh", "Mesh of the interval (a, b)");

  // Open mesh for editing
  MeshEditor editor;
  editor.open(*this, CellType::interval, 1, 1);

  // Create vertices and cells:
  editor.init_vertices_global((nx+1), (nx+1));
  editor.init_cells_global(nx, nx);

  // Create main vertices:
  for (std::size_t ix = 0; ix <= nx; ix++)
  {
    const std::vector<double>
      x(1, a + (static_cast<double>(ix)*(b - a)/static_cast<double>(nx)));
    editor.add_vertex(ix, x);
  }

  // Create intervals
  for (std::size_t ix = 0; ix < nx; ix++)
  {
    std::vector<std::size_t> cell(2);
    cell[0] = ix; cell[1] = ix + 1;
    editor.add_cell(ix, cell);
  }

  // Close mesh editor
  editor.close();

  // Broadcast mesh according to parallel policy
  if (MPI::is_broadcaster(this->mpi_comm()))
  {
    std::cout << "Building mesh (dist 0a)" << std::endl;
    MeshPartitioning::build_distributed_mesh(*this);
    std::cout << "Building mesh (dist 1a)" << std::endl;
    return;
  }
}

//-----------------------------------------------------------------------------
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();
}

//-----------------------------------------------------------------------------
void BoundaryComputation::compute_boundary(const Mesh& mesh,
                                           const std::string type,
                                           BoundaryMesh& boundary)
{
  // We iterate over all facets in the mesh and check if they are on
  // the boundary. A facet is on the boundary if it is connected to
  // exactly one cell.

  log(TRACE, "Computing boundary mesh.");

  bool exterior = true;
  bool interior = true;
  if (type == "exterior")
    interior = false;
  else if (type == "interior")
    exterior = false;
  else if (type != "local")
  {
    dolfin_error("BoundaryComputation.cpp",
                 "determine boundary mesh type",
                 "Unknown boundary type (%d)", type.c_str());
  }

  // Get my MPI process rank and number of MPI processes
  const std::size_t my_rank = MPI::rank(mesh.mpi_comm());
  const std::size_t num_processes = MPI::size(mesh.mpi_comm());

  // Open boundary mesh for editing
  const std::size_t D = mesh.topology().dim();
  MeshEditor editor;
  editor.open(boundary, mesh.type().facet_type(), D - 1, mesh.geometry().dim());

  // Generate facet - cell connectivity if not generated
  mesh.init(D - 1, D);

  // Temporary arrays for assignment of indices to vertices on the boundary
  std::map<std::size_t, std::size_t> boundary_vertices;

  // Map of index "owners" (process responsible for assigning global index)
  std::map< std::size_t, std::size_t > global_index_owner;

  // Shared vertices for full mesh
  // FIXME: const_cast
  const std::map<unsigned int, std::set<unsigned int>> &
    shared_vertices = const_cast<Mesh&>(mesh).topology().shared_entities(0);

  // Shared vertices for boundary mesh
  std::map<unsigned int, std::set<unsigned int>> shared_boundary_vertices;
  if (exterior)
  {
    // Extract shared vertices if vertex is identified as part of globally
    // exterior facet.
    std::vector<std::size_t> boundary_global_indices;
    for (std::map<unsigned int, std::set<unsigned int>>::const_iterator
        sv_it=shared_vertices.begin(); sv_it != shared_vertices.end(); ++sv_it)
    {
      std::size_t local_mesh_index = sv_it->first;
      Vertex v(mesh, local_mesh_index);

      for (FacetIterator f(v); !f.end(); ++f)
      {
        if (f->num_global_entities(D) == 1)
        {
          const std::size_t global_mesh_index
            = mesh.topology().global_indices(0)[local_mesh_index];
          shared_boundary_vertices[local_mesh_index] = sv_it->second;
          boundary_global_indices.push_back(global_mesh_index);
          break;
        }
      }
    }

    // Distribute all shared boundary vertices
    std::vector<std::vector<std::size_t>> boundary_global_indices_all;
    MPI::all_gather(mesh.mpi_comm(), boundary_global_indices,
                     boundary_global_indices_all);

    // Identify and clean up discrepancies between shared vertices of full mesh
    // and shared vertices of boundary mesh
    for (auto sbv_it = shared_boundary_vertices.begin();
         sbv_it != shared_boundary_vertices.end(); )
    {
      std::size_t local_mesh_index = sbv_it->first;
      const std::size_t global_mesh_index
        = mesh.topology().global_indices(0)[local_mesh_index];

      // Check if this vertex is identified as boundary vertex on
      // other processes sharing this vertex
      std::set<unsigned int> &other_processes = sbv_it->second;
      for (auto  op_it=other_processes.begin();
           op_it != other_processes.end(); )
      {
        // Check if vertex is identified as boundary vertex on process *op_it
        bool is_boundary_vertex
          = (std::find(boundary_global_indices_all[*op_it].begin(),
                      boundary_global_indices_all[*op_it].end(),
                      global_mesh_index)
             != boundary_global_indices_all[*op_it].end());

//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
bool LocalMeshCoarsening::coarsen_cell(Mesh& mesh, Mesh& coarse_mesh,
				      int cellid,
				      std::vector<int>& old2new_vertex,
				      std::vector<int>& old2new_cell,
				      bool coarsen_boundary)
{
  cout << "coarsen_cell: " << cellid << endl;
  cout << "num_cells: " << mesh.num_cells() << endl;

  const std::size_t num_vertices = mesh.size(0);
  const std::size_t num_cells = mesh.size(mesh.topology().dim());

  auto _mesh = reference_to_no_delete_pointer(mesh);

  // Initialise forbidden vertices
  MeshFunction<bool> vertex_forbidden(_mesh);
  vertex_forbidden.init(0);
  for (VertexIterator v(mesh); !v.end(); ++v)
    vertex_forbidden[v->index()] = false;

  // Initialise boundary vertices
  MeshFunction<bool> vertex_boundary(_mesh);
  vertex_boundary.init(0);
  for (VertexIterator v(mesh); !v.end(); ++v)
    vertex_boundary[v->index()] = false;

  BoundaryMesh boundary(mesh, "exterior");
  MeshFunction<std::size_t>& bnd_vertex_map = boundary.entity_map(0);
  for (VertexIterator v(boundary); !v.end(); ++v)
    vertex_boundary[bnd_vertex_map[v->index()]] = true;

  // If coarsen boundary is forbidden
  if (coarsen_boundary == false)
  {
    for (VertexIterator v(boundary); !v.end(); ++v)
      vertex_forbidden[bnd_vertex_map[v->index()]] = true;
  }

  // Initialise data for finding which vertex to remove
  bool _collapse_edge = false;
  const unsigned int* edge_vertex;
  std::size_t shortest_edge_index = 0;
  double lmin, l;
  std::size_t num_cells_to_remove = 0;

  // Get cell type
  const CellType& cell_type = mesh.type();
  const Cell cell(mesh, cellid);

  MeshEditor editor;
  editor.open(coarse_mesh, cell_type.cell_type(),
              mesh.topology().dim(), mesh.geometry().dim());

  MeshFunction<bool> cell_to_remove(_mesh);
  cell_to_remove.init(mesh.topology().dim());

  for (CellIterator ci(mesh); !ci.end(); ++ci)
    cell_to_remove[ci->index()] = false;

  MeshFunction<bool> cell_to_regenerate(_mesh);
  cell_to_regenerate.init(mesh.topology().dim());
  for (CellIterator ci(mesh); !ci.end(); ++ci)
    cell_to_regenerate[ci->index()] = false;

  // Find shortest edge of cell c
  _collapse_edge = false;
  lmin = 1.0e10*cell.diameter();
  for (EdgeIterator e(cell); !e.end(); ++e)
  {
    edge_vertex = e->entities(0);
    if (!vertex_forbidden[edge_vertex[0]] || !vertex_forbidden[edge_vertex[1]])
    {
      l = e->length();
      if ( lmin > l )
      {
        lmin = l;
        shortest_edge_index = e->index();
        _collapse_edge = true;
      }
    }
  }

  Edge shortest_edge(mesh, shortest_edge_index);

  // Decide which vertex to remove
  std::size_t vert2remove_idx = 0;

  // If at least one vertex should be removed
  if ( _collapse_edge == true )
  {
    edge_vertex = shortest_edge.entities(0);

    if(vertex_forbidden[edge_vertex[0]] &&
       vertex_forbidden[edge_vertex[1]])
    {
      // Both vertices are forbidden, cannot coarsen

      cout << "both vertices forbidden" << endl;

//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
void UniformMeshRefinement::refine(Mesh& refined_mesh,
                                   const Mesh& mesh)
{
  not_working_in_parallel("UniformMeshRefinement::refine");

  log(TRACE, "Refining simplicial mesh uniformly.");

  // Check that refined_mesh and mesh are not the same
  if (&refined_mesh == &mesh)
  {
    dolfin_error("UniformMeshRefinement.cpp",
                 "refine mesh",
                 "Refined_mesh and mesh point to the same object");
  }

  // Generate cell - edge connectivity if not generated
  mesh.init(mesh.topology().dim(), 1);

  // Generate edge - vertex connectivity if not generated
  mesh.init(1, 0);

  // Mesh needs to be ordered (so we can pick right combination of vertices/edges)
  if (!mesh.ordered())
    dolfin_error("UniformMeshRefinement.cpp",
                 "refine mesh",
                 "Mesh is not ordered according to the UFC numbering convention, consider calling mesh.order()");

  // Get cell type
  const CellType& cell_type = mesh.type();

  // Open new mesh for editing
  MeshEditor editor;
  editor.open(refined_mesh, cell_type.cell_type(),
              mesh.topology().dim(), mesh.geometry().dim());

  // Get size of mesh
  const std::size_t num_vertices = mesh.size(0);
  const std::size_t num_edges = mesh.size(1);
  const std::size_t num_cells = mesh.size(mesh.topology().dim());

  // Specify number of vertices and cells
  editor.init_vertices_global(num_vertices + num_edges, num_vertices + num_edges);
  editor.init_cells_global(ipow(2, mesh.topology().dim())*num_cells,
                           ipow(2, mesh.topology().dim())*num_cells);

  // Add old vertices
  std::size_t vertex = 0;
  for (VertexIterator v(mesh); !v.end(); ++v)
  {
    editor.add_vertex(vertex, v->point());
    vertex++;
  }

  // Add new vertices
  for (EdgeIterator e(mesh); !e.end(); ++e)
  {
    editor.add_vertex(vertex, e->midpoint());
    vertex++;
  }

  // Add cells
  std::size_t current_cell = 0;
  for (CellIterator c(mesh); !c.end(); ++c)
    cell_type.refine_cell(*c, editor, current_cell);

  // Close editor
  editor.close();

  // Make sure that mesh is ordered after refinement
  //refined_mesh.order();
}

int main(int argc, char *argv[]) {
  
  if (argc < 3 ) {
    std::cout<< "Usage: ./extract <checkpoint> <num>" << std::endl;
    return 0;
  }
  
  std::ostringstream _fname;
  _fname << argv[1] << "_";

  std::vector<std::string> chkp_files;

  for (int i = 0; i <  atoi(argv[2]); i++) {
    std::ostringstream tmp;
    tmp << _fname.str();
    tmp << i << ".chkp";
    chkp_files.push_back(tmp.str());
  }
  

  CellType::Type type;
  double t;
  unsigned int id, tdim, gdim, num_vertices, num_cells, num_entities;
  std::ifstream in;
  std::list<dVertex> vlist;
  std::list<dCell> dmesh;
  for(std::vector<std::string>::iterator it = chkp_files.begin(); 
      it != chkp_files.end(); ++it) { 

    in.open(it->c_str(), std::ifstream::binary);

    if(in.good()) {
      in.read((char *)&id, sizeof(unsigned int));
      in.read((char *)&t, sizeof(double));
      
      in.read((char *)&type, sizeof(CellType::Type));
      in.read((char *)&tdim, sizeof(unsigned int));
      in.read((char *)&gdim, sizeof(unsigned int));

      in.read((char *)&num_vertices, sizeof(unsigned int));
      in.read((char *)&num_cells, sizeof(unsigned int));
      in.read((char *)&num_entities, sizeof(unsigned int));

      double *coords = new double[gdim *num_vertices];
      in.read((char *)coords, (gdim * num_vertices) * sizeof(double));
      
      std::list<dVertex> _vlist;
      int vi = 0;
      for(unsigned int i = 0; i < gdim * num_vertices; i += gdim) {
	switch(gdim)
	{
	case 2:
	  _vlist.push_back(dVertex(vi++,  Point(coords[i], coords[i+1]))); break;
	case 3:
	  _vlist.push_back(dVertex(vi++,  Point(coords[i], coords[i+1], coords[i+2]))); break;
	}
      }
      delete[] coords;

      unsigned int *cells = new unsigned int[num_entities * num_cells];
      in.read((char *)cells, (num_entities * num_cells) * sizeof(unsigned int));

      unsigned int *mapping = new unsigned int[_vlist.size()];
      in.read((char *)mapping, _vlist.size() * sizeof(unsigned int));

      unsigned int *mp = &mapping[0];
      std::map<unsigned int, unsigned int> vmap;      
      for(std::list<dVertex>::iterator it = _vlist.begin(); 
	  it != _vlist.end(); ++it) 
	vmap[it->id] = *(mp++);
      delete[] mapping;	
      
      Array<unsigned int> v;
      for(unsigned int i = 0; i < num_entities * num_cells; i += num_entities) {
	v.clear();
	for(unsigned int j = 0; j < num_entities; j++)
	  v.push_back(vmap[cells[i+j]]);
	dmesh.push_back(dCell(v));
      }
      delete[] cells;

      unsigned int num_ghost;
      in.read((char *)&num_ghost, sizeof(unsigned int));
      unsigned int *ghosts = new unsigned int[2 * num_ghost];
      in.read((char *)ghosts, 2*num_ghost * sizeof(unsigned int));
      
      std::set<unsigned int> ghost_set;
      for (unsigned int i = 0; i < 2 * num_ghost; i += 2)
	ghost_set.insert(ghosts[i]);
      delete[] ghosts;

      for(std::list<dVertex>::iterator it = _vlist.begin(); 
	  it != _vlist.end(); ++it) 
	if(ghost_set.find(it->id) == ghost_set.end()) 
	  vlist.push_back(dVertex(vmap[it->id], it->p));
    }
    in.close();
  }

  Mesh mesh;
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
dolfin::Mesh MeshRenumbering::renumber_by_color(const Mesh& mesh,
                                 const std::vector<std::size_t> coloring_type)
{
  // Start timer
  Timer timer("Renumber mesh by color");

  // Get some some mesh
  const std::size_t tdim = mesh.topology().dim();
  const std::size_t gdim = mesh.geometry().dim();
  const std::size_t num_vertices = mesh.num_vertices();
  const std::size_t num_cells = mesh.num_cells();

  // Check that requested coloring is a cell coloring
  if (coloring_type[0] != tdim)
  {
    dolfin_error("MeshRenumbering.cpp",
                 "renumber mesh by color",
                 "Coloring is not a cell coloring: only cell colorings are supported");
  }

  // Compute renumbering
  std::vector<double> new_coordinates;
  std::vector<std::size_t> new_connections;
  MeshRenumbering::compute_renumbering(mesh, coloring_type, new_coordinates,
                                       new_connections);

  // Create new mesh
  Mesh new_mesh;

  // Create mesh editor
  MeshEditor editor;
  editor.open(new_mesh, mesh.type().cell_type(), tdim, gdim);
  editor.init_cells(num_cells);
  editor.init_vertices(num_vertices);

  // Add vertices
  dolfin_assert(new_coordinates.size() == num_vertices*gdim);
  for (std::size_t i = 0; i < num_vertices; ++i)
  {
    std::vector<double> x(gdim);
    for (std::size_t j = 0; j < gdim; ++j)
      x[j] = new_coordinates[i*gdim + j];
    editor.add_vertex(i, x);
  }

  cout << "Done adding vertices" << endl;

  // Add cells
  dolfin_assert(new_coordinates.size() == num_vertices*gdim);
  const std::size_t vertices_per_cell = mesh.type().num_entities(0);
  for (std::size_t i = 0; i < num_cells; ++i)
  {
    std::vector<std::size_t> c(vertices_per_cell);
    std::copy(new_connections.begin() + i*vertices_per_cell,
              new_connections.begin() + i*vertices_per_cell + vertices_per_cell,
              c.begin());
    editor.add_cell(i, c);
  }

  editor.close();

  cout << "Close editor" << endl;

  // Initialise coloring data
  typedef std::map<const std::vector<std::size_t>, std::pair<std::vector<std::size_t>,
           std::vector<std::vector<std::size_t> > > >::const_iterator ConstMeshColoringData;

  // Get old coloring
  ConstMeshColoringData mesh_coloring
    = mesh.topology().coloring.find(coloring_type);
  if (mesh_coloring == mesh.topology().coloring.end())
  {
    dolfin_error("MeshRenumbering.cpp",
                 "renumber mesh by color",
                 "Requested mesh coloring has not been computed");
  }

  // Get old coloring data
  const std::vector<std::size_t>& colors = mesh_coloring->second.first;
  const std::vector<std::vector<std::size_t> >&
    entities_of_color = mesh_coloring->second.second;
  dolfin_assert(colors.size() == num_cells);
  dolfin_assert(!entities_of_color.empty());
  const std::size_t num_colors = entities_of_color.size();

  // New coloring data
  dolfin_assert(new_mesh.topology().coloring.empty());
  std::vector<std::size_t> new_colors(colors.size());
  std::vector<std::vector<std::size_t> > new_entities_of_color(num_colors);

  std::size_t current_cell = 0;
  for (std::size_t color = 0; color < num_colors; color++)
  {
    // Get the array of cell indices of current color
    const std::vector<std::size_t>& colored_cells = entities_of_color[color];
    std::vector<std::size_t>& new_colored_cells = new_entities_of_color[color];

    // Update cell color data
    for (std::size_t i = 0; i < colored_cells.size(); i++)
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
void SubMesh::init(const Mesh& mesh,
                   const std::vector<std::size_t>& sub_domains,
                   std::size_t sub_domain)
{
  // Open mesh for editing
  MeshEditor editor;
  const std::size_t D = mesh.topology().dim();
  editor.open(*this, mesh.type().cell_type(), D,
              mesh.geometry().dim());

  // Build set of cells that are in sub-mesh
  std::vector<bool> parent_cell_in_subdomain(mesh.num_cells(), false);
  std::set<std::size_t> submesh_cells;
  for (CellIterator cell(mesh); !cell.end(); ++cell)
  {
    if (sub_domains[cell->index()] == sub_domain)
    {
      parent_cell_in_subdomain[cell->index()] = true;
      submesh_cells.insert(cell->index());
    }
  }

  // Map from parent vertex index to submesh vertex index
  std::map<std::size_t, std::size_t> parent_to_submesh_vertex_indices;

  // Map from submesh cell to parent cell
  std::vector<std::size_t> submesh_cell_parent_indices;
  submesh_cell_parent_indices.reserve(submesh_cells.size());

  // Vector from parent cell index to submesh cell index
  std::vector<std::size_t> parent_to_submesh_cell_indices(mesh.num_cells(), 0);

  // Add sub-mesh cells
  editor.init_cells_global(submesh_cells.size(), submesh_cells.size());
  std::size_t current_cell = 0;
  std::size_t current_vertex = 0;
  for (std::set<std::size_t>::iterator cell_it = submesh_cells.begin();
       cell_it != submesh_cells.end(); ++cell_it)
  {
    // Data structure to hold new vertex indices for cell
    std::vector<std::size_t> cell_vertices;

    // Create cell
    Cell cell(mesh, *cell_it);

    // Iterate over cell vertices
    for (VertexIterator vertex(cell); !vertex.end(); ++vertex)
    {
      const std::size_t parent_vertex_index = vertex->index();

      // Look for parent vertex in map
      std::map<std::size_t, std::size_t>::iterator vertex_it
        = parent_to_submesh_vertex_indices.find(parent_vertex_index);

      // If vertex has been inserted, get new index, otherwise
      // increment and insert
      std::size_t submesh_vertex_index = 0;
      if (vertex_it != parent_to_submesh_vertex_indices.end())
        submesh_vertex_index = vertex_it->second;
      else
      {
        submesh_vertex_index = current_vertex++;
        parent_to_submesh_vertex_indices[parent_vertex_index]
          = submesh_vertex_index;
      }

      // Add vertex to list of cell vertices (new indexing)
      cell_vertices.push_back(submesh_vertex_index);
    }

    // Add parent cell index to list
    submesh_cell_parent_indices.push_back(cell.index());

    // Store parent cell -> submesh cell indices
    parent_to_submesh_cell_indices[cell.index()] = current_cell;

    // Add cell to mesh
    editor.add_cell(current_cell++, cell_vertices);
  }

  // Vector to hold submesh vertex -> parent vertex
  std::vector<std::size_t> parent_vertex_indices;
  parent_vertex_indices.resize(parent_to_submesh_vertex_indices.size());

  // Initialise mesh editor
  editor.init_vertices_global(parent_to_submesh_vertex_indices.size(),
                              parent_to_submesh_vertex_indices.size());

  // Add vertices
  for (std::map<std::size_t, std::size_t>::iterator it
         = parent_to_submesh_vertex_indices.begin();
       it != parent_to_submesh_vertex_indices.end(); ++it)
  {
    Vertex vertex(mesh, it->first);
    if (MPI::size(mesh.mpi_comm()) > 1)
      error("SubMesh::init not working in parallel");

    // FIXME: Get global vertex index
    editor.add_vertex(it->second, vertex.point());
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
UnitHexMesh::UnitHexMesh(MPI_Comm comm, std::size_t nx, std::size_t ny,
                         std::size_t nz) : Mesh(comm)
{
  // Receive mesh according to parallel policy
  if (MPI::is_receiver(this->mpi_comm()))
  {
    MeshPartitioning::build_distributed_mesh(*this);
    return;
  }

  MeshEditor editor;
  editor.open(*this, CellType::hexahedron, 3, 3);

  // Create vertices and cells:
  editor.init_vertices_global((nx + 1)*(ny + 1)*(nz + 1),
                              (nx + 1)*(ny + 1)*(nz + 1));
  editor.init_cells_global(nx*ny*nz, nx*ny*nz);

  // Storage for vertices
  std::vector<double> x(3);

  const double a = 0.0;
  const double b = 1.0;
  const double c = 0.0;
  const double d = 1.0;
  const double e = 0.0;
  const double f = 1.0;

  // Create main vertices:
  std::size_t vertex = 0;
  for (std::size_t iz = 0; iz <= nz; iz++)
  {
    x[2] = e + ((static_cast<double>(iz))*(f - e)/static_cast<double>(nz));
    for (std::size_t iy = 0; iy <= ny; iy++)
    {
      x[1] = c + ((static_cast<double>(iy))*(d - c)/static_cast<double>(ny));
      for (std::size_t ix = 0; ix <= nx; ix++)
      {
        x[0] = a + ((static_cast<double>(ix))*(b - a)/static_cast<double>(nx));
        editor.add_vertex(vertex, x);
        vertex++;
      }
    }
  }

  // Create cuboids
  std::size_t cell = 0;
  std::vector<std::size_t> v(8);
  for (std::size_t iz = 0; iz < nz; iz++)
    for (std::size_t iy = 0; iy < ny; iy++)
      for (std::size_t ix = 0; ix < nx; ix++)
      {
        v[0] = (iz*(ny + 1) + iy)*(nx + 1) + ix;
        v[1] = v[0] + 1;
        v[2] = v[0] + (nx + 1);
        v[3] = v[1] + (nx + 1);
        v[4] = v[0] + (nx + 1)*(ny + 1);
        v[5] = v[1] + (nx + 1)*(ny + 1);
        v[6] = v[2] + (nx + 1)*(ny + 1);
        v[7] = v[3] + (nx + 1)*(ny + 1);
        editor.add_cell(cell, v);
        ++cell;
      }

  // Close mesh editor
  editor.close();

  // Broadcast mesh according to parallel policy
  if (MPI::is_broadcaster(this->mpi_comm()))
  {
    MeshPartitioning::build_distributed_mesh(*this);
    return;
  }


}

//-----------------------------------------------------------------------------
void RegularCutRefinement::refine_marked(Mesh& refined_mesh,
                                  const Mesh& mesh,
                                  const std::vector<int>& refinement_markers,
                                  const IndexSet& marked_edges)
{
  // Count the number of cells in refined mesh
  std::size_t num_cells = 0;

  // Data structure to hold a cell
  std::vector<std::size_t> cell_data(3);

  for (CellIterator cell(mesh); !cell.end(); ++cell)
  {
    const int marker = refinement_markers[cell->index()];
    switch (marker)
    {
    case no_refinement:
      num_cells += 1;
      break;
    case regular_refinement:
      num_cells += 4;
      break;
    case backtrack_bisection:
      num_cells += 2;
      break;
    case backtrack_bisection_refine:
      num_cells += 3;
      break;
    default:
      num_cells += 2;
    }
  }

  // Initialize mesh editor
  const std::size_t num_vertices = mesh.num_vertices() + marked_edges.size();
  MeshEditor editor;
  editor.open(refined_mesh, mesh.topology().dim(), mesh.geometry().dim());
  editor.init_vertices(num_vertices);
  editor.init_cells(num_cells);

  // Set vertex coordinates
  std::size_t current_vertex = 0;
  for (VertexIterator vertex(mesh); !vertex.end(); ++vertex)
  {
    editor.add_vertex(current_vertex, vertex->point());
    current_vertex++;
  }
  for (std::size_t i = 0; i < marked_edges.size(); i++)
  {
    Edge edge(mesh, marked_edges[i]);
    editor.add_vertex(current_vertex, edge.midpoint());
    current_vertex++;
  }

  // Get bisection data for old mesh
  const std::size_t D = mesh.topology().dim();
  const std::vector<std::size_t>*  bisection_twins = NULL;
  if (mesh.data().exists("bisection_twins", D))
    bisection_twins = &(mesh.data().array("bisection_twins", D));

  // Markers for bisected cells pointing to their bisection twins in
  // refined mesh
  std::vector<std::size_t>& refined_bisection_twins
    = refined_mesh.data().create_array("bisection_twins", D);
  refined_bisection_twins.resize(num_cells);
  for (std::size_t i = 0; i < num_cells; i++)
    refined_bisection_twins[i] = i;

  // Mapping from old to new unrefined cells (-1 means refined or not
  // yet processed)
  std::vector<int> unrefined_cells(mesh.num_cells());
  std::fill(unrefined_cells.begin(), unrefined_cells.end(), -1);

  // Iterate over all cells and add new cells
  std::size_t current_cell = 0;
  std::vector<std::vector<std::size_t> > cells(4, std::vector<std::size_t>(3));
  for (CellIterator cell(mesh); !cell.end(); ++cell)
  {
    // Get marker
    const int marker = refinement_markers[cell->index()];

    if (marker == no_refinement)
    {
      // No refinement: just copy cell to new mesh

      std::vector<std::size_t> vertices;
      for (VertexIterator vertex(*cell); !vertex.end(); ++vertex)
        vertices.push_back(vertex->index());
      editor.add_cell(current_cell++, vertices);

      // Store mapping to new cell index
      unrefined_cells[cell->index()] = current_cell - 1;

      // Remember unrefined bisection twins
      if (bisection_twins)
      {
        const std::size_t bisection_twin = (*bisection_twins)[cell->index()];
        const int twin_marker = refinement_markers[bisection_twin];
        dolfin_assert(twin_marker == no_refinement);
//.........这里部分代码省略.........