C++ MeshEditor类代码示例

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

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

//-----------------------------------------------------------------------------
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 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 RectangleMesh::build(const Point& p0, const Point& p1,
                          std::size_t nx, std::size_t ny,
                          std::string diagonal)
{
  // Receive mesh according to parallel policy
  if (MPI::is_receiver(this->mpi_comm()))
  {
    MeshPartitioning::build_distributed_mesh(*this);
    return;
  }

  // Check options
  if (diagonal != "left" && diagonal != "right" && diagonal != "right/left"
          && diagonal != "left/right" && diagonal != "crossed")
  {
    dolfin_error("RectangleMesh.cpp",
                 "create rectangle",
                 "Unknown mesh diagonal definition: allowed options are \"left\", \"right\", \"left/right\", \"right/left\" and \"crossed\"");
  }

  // Extract minimum and maximum coordinates
  const double x0 = std::min(p0.x(), p1.x());
  const double x1 = std::max(p0.x(), p1.x());
  const double y0 = std::min(p0.y(), p1.y());
  const double y1 = std::max(p0.y(), p1.y());

  const double a = x0;
  const double b = x1;
  const double c = y0;
  const double d = y1;

  if (std::abs(x0 - x1) < DOLFIN_EPS || std::abs(y0 - y1) < DOLFIN_EPS)
  {
    dolfin_error("Rectangle.cpp",
                 "create rectangle",
                 "Rectangle seems to have zero width, height or depth. Consider checking your dimensions");
  }

  if (nx < 1 || ny < 1)
  {
    dolfin_error("RectangleMesh.cpp",
                 "create rectangle",
                 "Rectangle has non-positive number of vertices in some dimension: number of vertices must be at least 1 in each dimension");
  }

  rename("mesh", "Mesh of the unit square (a,b) x (c,d)");
  // Open mesh for editing
  MeshEditor editor;
  editor.open(*this, CellType::triangle, 2, 2);

  // Create vertices and cells:
  if (diagonal == "crossed")
  {
    editor.init_vertices_global((nx + 1)*(ny + 1) + nx*ny,
                                  (nx + 1)*(ny + 1) + nx*ny);
    editor.init_cells_global(4*nx*ny, 4*nx*ny);
  }
  else
  {
    editor.init_vertices_global((nx + 1)*(ny + 1), (nx + 1)*(ny + 1));
    editor.init_cells_global(2*nx*ny, 2*nx*ny);
  }

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

  // Create main vertices:
  std::size_t vertex = 0;
  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 midpoint vertices if the mesh type is crossed
  if (diagonal == "crossed")
  {
    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]);
//.........这里部分代码省略.........

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

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

//-----------------------------------------------------------------------------
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++)
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
    using namespace Eigen;
    using EigMat = SparseMatrix<double>;
    using EigVec = VectorXd;
    
    EigMat A(6,6);
    for(int i=0;i<6;++i)
        A.insert(i, i) = 1.0/i;
    

    SimplicialLLT<EigMat> solver(A);

    EigVec b(6);
    b << 1,2,3,4,5,6;
    EigVec X = solver.solve(b);
    
    cout << X << endl;


    ArgExtracter ae(argc, argv);
	
    init_glut(argc, argv);
    init_gl();
    
    me.init();
    
    register_console_funcs(&me);
    
    glutMainLoop();
    return 0;
}