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

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 TriangleCell::refine_cell(Cell& cell, MeshEditor& editor,
                               std::size_t& current_cell) const
{
  // Get vertices and edges
  const unsigned int* v = cell.entities(0);
  const unsigned int* e = cell.entities(1);
  dolfin_assert(v);
  dolfin_assert(e);

  // Get offset for new vertex indices
  const std::size_t offset = cell.mesh().num_vertices();

  // Compute indices for the six new vertices
  const std::size_t v0 = v[0];
  const std::size_t v1 = v[1];
  const std::size_t v2 = v[2];
  const std::size_t e0 = offset + e[find_edge(0, cell)];
  const std::size_t e1 = offset + e[find_edge(1, cell)];
  const std::size_t e2 = offset + e[find_edge(2, cell)];

  // Create four new cells
  std::vector<std::vector<std::size_t>> cells(4, std::vector<std::size_t>(3));
  cells[0][0] = v0; cells[0][1] = e2; cells[0][2] = e1;
  cells[1][0] = v1; cells[1][1] = e0; cells[1][2] = e2;
  cells[2][0] = v2; cells[2][1] = e1; cells[2][2] = e0;
  cells[3][0] = e0; cells[3][1] = e1; cells[3][2] = e2;

  // Add cells
  std::vector<std::vector<std::size_t>>::const_iterator _cell;
  for (_cell = cells.begin(); _cell != cells.end(); ++_cell)
    editor.add_cell(current_cell++, *_cell);
}

//-----------------------------------------------------------------------------
void IntervalCell::refine_cell(Cell& cell, MeshEditor& editor,
                               std::size_t& current_cell) const
{
  // Get vertices
  const unsigned int* v = cell.entities(0);
  dolfin_assert(v);

  // Get offset for new vertex indices
  const std::size_t offset = cell.mesh().num_vertices();

  // Compute indices for the three new vertices
  const std::size_t v0 = v[0];
  const std::size_t v1 = v[1];
  const std::size_t e0 = offset + cell.index();

  // Add the two new cells
  std::vector<std::size_t> new_cell(2);

  new_cell[0] = v0; new_cell[1] = e0;
  editor.add_cell(current_cell++, new_cell);

  new_cell[0] = e0; new_cell[1] = v1;
  editor.add_cell(current_cell++, new_cell);
}

//-----------------------------------------------------------------------------
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 LocalMeshCoarsening::collapse_edge(Mesh& mesh, Edge& edge,
                                       Vertex& vertex_to_remove,
                                       MeshFunction<bool>& cell_to_remove,
                                       std::vector<int>& old2new_vertex,
                                       std::vector<int>& old2new_cell,
                                       MeshEditor& editor,
                                       std::size_t& current_cell)
{
  const CellType& cell_type = mesh.type();
  std::vector<std::size_t> cell_vertices(cell_type.num_entities(0));

  std::size_t vert_slave = vertex_to_remove.index();
  std::size_t vert_master = 0;
  const unsigned int* edge_vertex = edge.entities(0);

  if ( edge_vertex[0] == vert_slave )
    vert_master = edge_vertex[1];
  else if ( edge_vertex[1] == vert_slave )
    vert_master = edge_vertex[0];
  else
    dolfin_error("LocalMeshCoarsening.cpp",
                 "collapse edge",
                 "The node to delete and edge to collapse are not compatible");

  for (CellIterator c(vertex_to_remove); !c.end(); ++c)
  {
    if ( cell_to_remove[*c] == false )
    {
      std::size_t cv_idx = 0;
      for (VertexIterator v(*c); !v.end(); ++v)
      {
        if ( v->index() == vert_slave )
          cell_vertices[cv_idx++] = old2new_vertex[vert_master];
        else
          cell_vertices[cv_idx++] = old2new_vertex[v->index()];
      }
      //cout << "adding new cell" << endl;
      editor.add_cell(current_cell++, cell_vertices);

      // Update cell map
      old2new_cell[c->index()] = current_cell - 1;
    }
  }

}

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

}

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


}

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

//.........这里部分代码省略.........
      global_indices[global_mesh_index] = current_index++;
  }

  // Send and receive requests from other processes
  std::vector<std::vector<std::size_t>> global_index_requests(num_processes);
  MPI::all_to_all(mesh.mpi_comm(), request_global_indices,
                  global_index_requests);

  // Find response to requests of global indices
  std::vector<std::vector<std::size_t>> respond_global_indices(num_processes);
  for (std::size_t i = 0; i < num_processes; i++)
  {
    const std::size_t N = global_index_requests[i].size();
    respond_global_indices[i].resize(N);

    for (std::size_t j = 0; j < N; j++)
      respond_global_indices[i][j]
        = global_indices[global_index_requests[i][j]];
  }

  // Scatter responses back to requesting processes
  std::vector<std::vector<std::size_t>> global_index_responses(num_processes);
  MPI::all_to_all(mesh.mpi_comm(), respond_global_indices,
                  global_index_responses);

  // Update global_indices
  for (std::size_t i = 0; i < num_processes; i++)
  {
    const std::size_t N = global_index_responses[i].size();
    // Check that responses are the same size as the requests made
    dolfin_assert(global_index_responses[i].size()
                  == request_global_indices[i].size());
    for (std::size_t j = 0; j < N; j++)
    {
      const std::size_t global_mesh_index = request_global_indices[i][j];
      const std::size_t global_boundary_index = global_index_responses[i][j];
      global_indices[global_mesh_index] = global_boundary_index;
    }
  }

  // Create vertices
  for (std::size_t local_boundary_index = 0;
       local_boundary_index < num_boundary_vertices; local_boundary_index++)
  {
    const std::size_t local_mesh_index = vertex_map[local_boundary_index];
    const std::size_t global_mesh_index
      = mesh.topology().global_indices(0)[local_mesh_index];
    const std::size_t global_boundary_index = global_indices[global_mesh_index];

    Vertex v(mesh, local_mesh_index);

    editor.add_vertex_global(local_boundary_index, global_boundary_index,
                             v.point());
  }

  // Find global index to start cell numbering from for current process
  std::vector<std::size_t> cell_distribution(num_processes);
  MPI::all_gather(mesh.mpi_comm(), num_boundary_cells, cell_distribution);
  std::size_t start_cell_index = 0;
  for (std::size_t i = 0; i < my_rank; i++)
    start_cell_index += cell_distribution[i];

  // Create cells (facets) and map between boundary mesh cells and facets parent
  MeshFunction<std::size_t>& cell_map = boundary.entity_map(D - 1);
  if (num_boundary_cells > 0)
  {
    cell_map.init(reference_to_no_delete_pointer(boundary), D - 1,
                  num_boundary_cells);
  }
  std::vector<std::size_t>
    cell(boundary.type().num_vertices(boundary.topology().dim()));
  std::size_t current_cell = 0;
  for (FacetIterator f(mesh); !f.end(); ++f)
  {
    if (boundary_facet[*f])
    {
      // Compute new vertex numbers for cell
      const unsigned int* vertices = f->entities(0);
      for (std::size_t i = 0; i < cell.size(); i++)
        cell[i] = boundary_vertices[vertices[i]];

      // Reorder vertices so facet is right-oriented w.r.t. facet
      // normal
      reorder(cell, *f);

      // Create mapping from boundary cell to mesh facet if requested
      if (!cell_map.empty())
        cell_map[current_cell] = f->index();

      // Add cell
      editor.add_cell(current_cell, start_cell_index+current_cell, cell);
      current_cell++;
    }
  }

  // Close mesh editor. Note the argument order=false to prevent
  // ordering from destroying the orientation of facets accomplished
  // by calling reorder() below.
  editor.close(false);
}

//.........这里部分代码省略.........
    // Update cell map (will be filled in with correct index)
    old2new_cell[cn->index()] = -1;
  }

  // Specify number of vertices and cells
  editor.init_vertices_global(num_vertices - 1, num_vertices - 1);
  editor.init_cells_global(num_cells - num_cells_to_remove,
                           num_cells - num_cells_to_remove);

  cout << "Number of cells in old mesh: " << num_cells << "; to remove: " <<
    num_cells_to_remove << endl;

  // Add old vertices
  std::size_t vertex = 0;
  for(VertexIterator v(mesh); !v.end(); ++v)
  {
    if(vertex_to_remove.index() == v->index())
      old2new_vertex[v->index()] = -1;
    else
    {
      //cout << "adding old vertex at: " << v->point() << endl;
      old2new_vertex[v->index()] = vertex;
      editor.add_vertex(vertex, v->point());
      vertex++;
    }
  }

  // Add old unrefined cells
  std::size_t cv_idx;
  std::size_t current_cell = 0;
  std::vector<std::size_t> cell_vertices(cell_type.num_entities(0));
  for (CellIterator c(mesh); !c.end(); ++c)
  {
    if (cell_to_remove[*c] == false && cell_to_regenerate[*c] == false)
    {
      cv_idx = 0;
      for (VertexIterator v(*c); !v.end(); ++v)
        cell_vertices[cv_idx++] = old2new_vertex[v->index()];
      //cout << "adding old cell" << endl;
      editor.add_cell(current_cell++, cell_vertices);

      // Update cell maps
      old2new_cell[c->index()] = current_cell - 1;
    }
  }

  // Add new cells.
  collapse_edge(mesh, shortest_edge, vertex_to_remove, cell_to_remove,
                old2new_vertex, old2new_cell, editor, current_cell);

  editor.close();

  // Set volume tolerance. This parameter determines a quality criterion
  // for the new mesh: higher value indicates a sharper criterion.
  double vol_tol = 1.0e-3;

  bool mesh_ok = true;

  Cell removed_cell(mesh, cellid);

  // Check mesh quality (volume)
  for (CellIterator c(removed_cell); !c.end(); ++c)
  {
    std::size_t id = c->index();
    int nid = old2new_cell[id];

    if(nid != -1)
    {
      Cell cn(coarse_mesh, nid);
      double qm = cn.volume() / cn.diameter();
      if(qm < vol_tol)
      {
        warning("Cell quality too low");
        cout << "qm: " << qm << endl;
        mesh_ok = false;
        return mesh_ok;
      }
    }
  }

  // Checking for inverted cells
  for (CellIterator c(removed_cell); !c.end(); ++c)
  {
    std::size_t id = c->index();
    int nid = old2new_cell[id];

    if(nid != -1)
    {
      Cell cn(coarse_mesh, nid);

      if(c->orientation() != cn.orientation())
      {
        cout << "cell orientation inverted" << endl;
        mesh_ok = false;
        return mesh_ok;
      }
    }
  }
  return mesh_ok;
}

//.........这里部分代码省略.........
  {
    for (std::size_t iy = 0; iy < ny; iy++)
    {
      x[1] = c +(static_cast<double>(iy) + 0.5)*(d - c)/static_cast<double>(ny);
      for (std::size_t ix = 0; ix < nx; ix++)
      {
        x[0] = a + (static_cast<double>(ix) + 0.5)*(b - a)/static_cast<double>(nx);
        editor.add_vertex(vertex, x);
        vertex++;
      }
    }
  }

  // Create triangles
  std::size_t cell = 0;
  if (diagonal == "crossed")
  {
    boost::multi_array<std::size_t, 2> cells(boost::extents[4][3]);
    for (std::size_t iy = 0; iy < ny; iy++)
    {
      for (std::size_t ix = 0; ix < nx; ix++)
      {
        const std::size_t v0 = iy*(nx + 1) + ix;
        const std::size_t v1 = v0 + 1;
        const std::size_t v2 = v0 + (nx + 1);
        const std::size_t v3 = v1 + (nx + 1);
        const std::size_t vmid = (nx + 1)*(ny + 1) + iy*nx + ix;

        // Note that v0 < v1 < v2 < v3 < vmid.
        cells[0][0] = v0; cells[0][1] = v1; cells[0][2] = vmid;
        cells[1][0] = v0; cells[1][1] = v2; cells[1][2] = vmid;
        cells[2][0] = v1; cells[2][1] = v3; cells[2][2] = vmid;
        cells[3][0] = v2; cells[3][1] = v3; cells[3][2] = vmid;

        // Add cells
        for (auto _cell = cells.begin(); _cell != cells.end(); ++_cell)
          editor.add_cell(cell++, *_cell);
      }
    }
  }
  else if (diagonal == "left" || diagonal == "right" || diagonal == "right/left" || diagonal == "left/right")
  {
    std::string local_diagonal = diagonal;
    boost::multi_array<std::size_t, 2> cells(boost::extents[2][3]);
    for (std::size_t iy = 0; iy < ny; iy++)
    {
      // Set up alternating diagonal
      if (diagonal == "right/left")
      {
        if (iy % 2)
          local_diagonal = "right";
        else
          local_diagonal = "left";
      }
      if (diagonal == "left/right")
      {
        if (iy % 2)
          local_diagonal = "left";
        else
          local_diagonal = "right";
      }

      for (std::size_t ix = 0; ix < nx; ix++)
      {
        const std::size_t v0 = iy*(nx + 1) + ix;
        const std::size_t v1 = v0 + 1;
        const std::size_t v2 = v0 + (nx + 1);
        const std::size_t v3 = v1 + (nx + 1);
        std::vector<std::size_t> cell_data;

        if(local_diagonal == "left")
        {
          cells[0][0] = v0; cells[0][1] = v1; cells[0][2] = v2;
          cells[1][0] = v1; cells[1][1] = v2; cells[1][2] = v3;
          if (diagonal == "right/left" || diagonal == "left/right")
            local_diagonal = "right";
        }
        else
        {
          cells[0][0] = v0; cells[0][1] = v1; cells[0][2] = v3;
          cells[1][0] = v0; cells[1][1] = v2; cells[1][2] = v3;
          if (diagonal == "right/left" || diagonal == "left/right")
            local_diagonal = "left";
        }
        editor.add_cell(cell++, cells[0]);
        editor.add_cell(cell++, cells[1]);
      }
    }
  }

  // Close mesh editor
  editor.close();

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

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