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

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

}

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

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

//-----------------------------------------------------------------------------
void MeshPartitioning::build_mesh(Mesh& mesh,
              const std::vector<std::size_t>& global_cell_indices,
              const boost::multi_array<std::size_t, 2>& cell_global_vertices,
              const std::vector<std::size_t>& vertex_indices,
              const boost::multi_array<double, 2>& vertex_coordinates,
              const std::map<std::size_t, std::size_t>& vertex_global_to_local,
              std::size_t tdim, std::size_t gdim, std::size_t num_global_cells,
              std::size_t num_global_vertices)
{
  Timer timer("PARALLEL 3: Build mesh (from local mesh data)");

  // Get number of processes and process number
  const std::size_t num_processes = MPI::num_processes();
  const std::size_t process_number = MPI::process_number();

  // Open mesh for editing
  mesh.clear();
  MeshEditor editor;
  editor.open(mesh, tdim, gdim);

  // Add vertices
  editor.init_vertices(vertex_coordinates.size());
  Point point(gdim);
  dolfin_assert(vertex_indices.size() == vertex_coordinates.size());
  for (std::size_t i = 0; i < vertex_coordinates.size(); ++i)
  {
    for (std::size_t j = 0; j < gdim; ++j)
      point[j] = vertex_coordinates[i][j];
    editor.add_vertex_global(i, vertex_indices[i], point);
  }

  // Add cells
  editor.init_cells(cell_global_vertices.size());
  const std::size_t num_cell_vertices = tdim + 1;
  std::vector<std::size_t> cell(num_cell_vertices);
  for (std::size_t i = 0; i < cell_global_vertices.size(); ++i)
  {
    for (std::size_t j = 0; j < num_cell_vertices; ++j)
    {
      // Get local cell vertex
      std::map<std::size_t, std::size_t>::const_iterator iter
          = vertex_global_to_local.find(cell_global_vertices[i][j]);
      dolfin_assert(iter != vertex_global_to_local.end());
      cell[j] = iter->second;
    }
    editor.add_cell(i, global_cell_indices[i], cell);
  }

  // Close mesh: Note that this must be done after creating the global
  // vertex map or otherwise the ordering in mesh.close() will be wrong
  // (based on local numbers).
  editor.close();

  // Set global number of cells and vertices
  mesh.topology().init_global(0, num_global_vertices);
  mesh.topology().init_global(tdim,  num_global_cells);

  // Construct boundary mesh
  BoundaryMesh bmesh(mesh, "exterior");

  const MeshFunction<std::size_t>& boundary_vertex_map = bmesh.entity_map(0);
  const std::size_t boundary_size = boundary_vertex_map.size();

  // Build sorted array of global boundary vertex indices (global
  // numbering)
  std::vector<std::size_t> global_vertex_send(boundary_size);
  for (std::size_t i = 0; i < boundary_size; ++i)
    global_vertex_send[i] = vertex_indices[boundary_vertex_map[i]];
  std::sort(global_vertex_send.begin(), global_vertex_send.end());

  // Receive buffer
  std::vector<std::size_t> global_vertex_recv;

  // Create shared_vertices data structure: mapping from shared vertices
  // to list of neighboring processes
  std::map<unsigned int, std::set<unsigned int> >& shared_vertices
        = mesh.topology().shared_entities(0);
  shared_vertices.clear();

  // FIXME: Remove computation from inside communication loop

  // Build shared vertex to sharing processes map
  for (std::size_t i = 1; i < num_processes; ++i)
  {
    // We send data to process p - i (i steps to the left)
    const int p = (process_number - i + num_processes) % num_processes;

    // We receive data from process p + i (i steps to the right)
    const int q = (process_number + i) % num_processes;

    // Send and receive
    MPI::send_recv(global_vertex_send, p, global_vertex_recv, q);

    // Compute intersection of global indices
    std::vector<std::size_t> intersection(std::min(global_vertex_send.size(),
                                                   global_vertex_recv.size()));
    std::vector<std::size_t>::iterator intersection_end
      = std::set_intersection(global_vertex_send.begin(),
                              global_vertex_send.end(),
//.........这里部分代码省略.........