C++ MeshBase::n_elem方法代码示例

void CentroidPartitioner::compute_centroids (MeshBase& mesh)
{
  _elem_centroids.clear();
  _elem_centroids.reserve(mesh.n_elem());

  //   elem_iterator it(mesh.elements_begin());
  //   const elem_iterator it_end(mesh.elements_end());

  MeshBase::element_iterator       it     = mesh.elements_begin();
  const MeshBase::element_iterator it_end = mesh.elements_end();

  for (; it != it_end; ++it)
    {
      Elem* elem = *it;

      _elem_centroids.push_back(std::make_pair(elem->centroid(), elem));
    }
}

// ------------------------------------------------------------
// SFCPartitioner implementation
void SFCPartitioner::_do_partition (MeshBase & mesh,
                                    const unsigned int n)
{

  libmesh_assert_greater (n, 0);

  // Check for an easy return
  if (n == 1)
    {
      this->single_partition (mesh);
      return;
    }

  // What to do if the sfcurves library IS NOT present
#ifndef LIBMESH_HAVE_SFCURVES

  libmesh_here();
  libMesh::err << "ERROR: The library has been built without"    << std::endl
               << "Space Filling Curve support.  Using a linear" << std::endl
               << "partitioner instead!" << std::endl;

  LinearPartitioner lp;

  lp.partition (mesh, n);

  // What to do if the sfcurves library IS present
#else

  LOG_SCOPE("sfc_partition()", "SFCPartitioner");

  const dof_id_type n_active_elem = mesh.n_active_elem();
  const dof_id_type n_elem        = mesh.n_elem();

  // the forward_map maps the active element id
  // into a contiguous block of indices
  std::vector<dof_id_type>
    forward_map (n_elem, DofObject::invalid_id);

  // the reverse_map maps the contiguous ids back
  // to active elements
  std::vector<Elem *> reverse_map (n_active_elem, libmesh_nullptr);

  int size = static_cast<int>(n_active_elem);
  std::vector<double> x      (size);
  std::vector<double> y      (size);
  std::vector<double> z      (size);
  std::vector<int>    table  (size);


  // We need to map the active element ids into a
  // contiguous range.
  {
    MeshBase::element_iterator       elem_it  = mesh.active_elements_begin();
    const MeshBase::element_iterator elem_end = mesh.active_elements_end();

    dof_id_type el_num = 0;

    for (; elem_it != elem_end; ++elem_it)
      {
        libmesh_assert_less ((*elem_it)->id(), forward_map.size());
        libmesh_assert_less (el_num, reverse_map.size());

        forward_map[(*elem_it)->id()] = el_num;
        reverse_map[el_num]           = *elem_it;
        el_num++;
      }
    libmesh_assert_equal_to (el_num, n_active_elem);
  }


  // Get the centroid for each active element
  {
    //     const_active_elem_iterator       elem_it (mesh.const_elements_begin());
    //     const const_active_elem_iterator elem_end(mesh.const_elements_end());

    MeshBase::element_iterator       elem_it  = mesh.active_elements_begin();
    const MeshBase::element_iterator elem_end = mesh.active_elements_end();

    for (; elem_it != elem_end; ++elem_it)
      {
        const Elem * elem = *elem_it;

        libmesh_assert_less (elem->id(), forward_map.size());

        const Point p = elem->centroid();

        x[forward_map[elem->id()]] = p(0);
        y[forward_map[elem->id()]] = p(1);
        z[forward_map[elem->id()]] = p(2);
      }
  }

  // build the space-filling curve
  if (_sfc_type == "Hilbert")
    Sfc::hilbert (&x[0], &y[0], &z[0], &size, &table[0]);

  else if (_sfc_type == "Morton")
    Sfc::morton  (&x[0], &y[0], &z[0], &size, &table[0]);
//.........这里部分代码省略.........

//---------------------------------------------------------
// CentroidPartitioner methods
void CentroidPartitioner::_do_partition (MeshBase& mesh,
                                         const unsigned int n)
{
  // Check for an easy return
  if (n == 1)
    {
      this->single_partition (mesh);
      return;
    }


  // Possibly reconstruct centroids
  if (mesh.n_elem() != _elem_centroids.size())
    this->compute_centroids (mesh);



  switch (this->sort_method())
    {
    case X:
      {
        std::sort(_elem_centroids.begin(),
                  _elem_centroids.end(),
                  CentroidPartitioner::sort_x);

        break;
      }


    case Y:
      {
        std::sort(_elem_centroids.begin(),
                  _elem_centroids.end(),
                  CentroidPartitioner::sort_y);

        break;

      }


    case Z:
      {
        std::sort(_elem_centroids.begin(),
                  _elem_centroids.end(),
                  CentroidPartitioner::sort_z);

        break;
      }


    case RADIAL:
      {
        std::sort(_elem_centroids.begin(),
                  _elem_centroids.end(),
                  CentroidPartitioner::sort_radial);

        break;
      }
    default:
      libmesh_error();
    }


  // Make sure the user has not handed us an
  // invalid number of partitions.
  libmesh_assert_greater (n, 0);

  // the number of elements, e.g. 1000
  const dof_id_type n_elem      = mesh.n_elem();
  // the number of elements per processor, e.g 400
  const dof_id_type target_size = n_elem / n;

  // Make sure the mesh hasn't changed since the
  // last time we computed the centroids.
  libmesh_assert_equal_to (mesh.n_elem(), _elem_centroids.size());

  for (dof_id_type i=0; i<n_elem; i++)
    {
      Elem* elem = _elem_centroids[i].second;

      elem->processor_id() =
        std::min (libmesh_cast_int<processor_id_type>(i / target_size),
                  libmesh_cast_int<processor_id_type>(n-1));
    }
}

void MeshTools::Subdivision::all_subdivision(MeshBase & mesh)
{
  std::vector<Elem *> new_elements;
  new_elements.reserve(mesh.n_elem());
  const bool mesh_has_boundary_data =
    (mesh.get_boundary_info().n_boundary_ids() > 0);

  std::vector<Elem *> new_boundary_elements;
  std::vector<short int> new_boundary_sides;
  std::vector<boundary_id_type> new_boundary_ids;

  MeshBase::const_element_iterator       el     = mesh.elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.elements_end();
  for (; el != end_el; ++el)
    {
      const Elem * elem = *el;
      libmesh_assert_equal_to(elem->type(), TRI3);

      Elem * tri = new Tri3Subdivision;
      tri->set_id(elem->id());
      tri->subdomain_id() = elem->subdomain_id();
      tri->set_node(0) = (*el)->get_node(0);
      tri->set_node(1) = (*el)->get_node(1);
      tri->set_node(2) = (*el)->get_node(2);

      if (mesh_has_boundary_data)
        {
          for (unsigned short side = 0; side < elem->n_sides(); ++side)
            {
              const boundary_id_type boundary_id =
                mesh.get_boundary_info().boundary_id(elem, side);
              if (boundary_id != BoundaryInfo::invalid_id)
                {
                  // add the boundary id to the list of new boundary ids
                  new_boundary_ids.push_back(boundary_id);
                  new_boundary_elements.push_back(tri);
                  new_boundary_sides.push_back(side);
                }
            }

          // remove the original element from the BoundaryInfo structure
          mesh.get_boundary_info().remove(elem);
        }

      new_elements.push_back(tri);
      mesh.insert_elem(tri);
    }
  mesh.prepare_for_use();

  if (mesh_has_boundary_data)
    {
      // If the old mesh had boundary data, the new mesh better have some too.
      libmesh_assert_greater(new_boundary_elements.size(), 0);

      // We should also be sure that the lengths of the new boundary data vectors
      // are all the same.
      libmesh_assert_equal_to(new_boundary_sides.size(), new_boundary_elements.size());
      libmesh_assert_equal_to(new_boundary_sides.size(), new_boundary_ids.size());

      // Add the new boundary info to the mesh.
      for (unsigned int s = 0; s < new_boundary_elements.size(); ++s)
        mesh.get_boundary_info().add_side(new_boundary_elements[s],
                                          new_boundary_sides[s],
                                          new_boundary_ids[s]);
    }

  mesh.prepare_for_use();
}

void MeshTools::Subdivision::add_boundary_ghosts(MeshBase & mesh)
{
  static const Real tol = 1e-5;

  // add the mirrored ghost elements (without using iterators, because the mesh is modified in the course)
  std::vector<Tri3Subdivision *> ghost_elems;
  std::vector<Node *> ghost_nodes;
  const unsigned int n_elem = mesh.n_elem();
  for (unsigned int eid = 0; eid < n_elem; ++eid)
    {
      Elem * elem = mesh.elem(eid);
      libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);

      // If the triangle happens to be in a corner (two boundary
      // edges), we perform a counter-clockwise loop by mirroring the
      // previous triangle until we come back to the original
      // triangle.  This prevents degenerated triangles in the mesh
      // corners and guarantees that the node in the middle of the
      // loop is of valence=6.
      for (unsigned int i = 0; i < elem->n_sides(); ++i)
        {
          libmesh_assert_not_equal_to(elem->neighbor(i), elem);

          if (elem->neighbor(i) == libmesh_nullptr &&
              elem->neighbor(next[i]) == libmesh_nullptr)
            {
              Elem * nelem = elem;
              unsigned int k = i;
              for (unsigned int l=0;l<4;l++)
                {
                  // this is the vertex to be mirrored
                  Point point = nelem->point(k) + nelem->point(next[k]) - nelem->point(prev[k]);

                  // Check if the proposed vertex doesn't coincide
                  // with one of the existing vertices.  This is
                  // necessary because for some triangulations, it can
                  // happen that two mirrored ghost vertices coincide,
                  // which would then lead to a zero size ghost
                  // element below.
                  Node * node = libmesh_nullptr;
                  for (unsigned int j = 0; j < ghost_nodes.size(); ++j)
                    {
                      if ((*ghost_nodes[j] - point).size() < tol * (elem->point(k) - point).size())
                        {
                          node = ghost_nodes[j];
                          break;
                        }
                    }

                  // add the new vertex only if no other is nearby
                  if (node == libmesh_nullptr)
                    {
                      node = mesh.add_point(point);
                      ghost_nodes.push_back(node);
                    }

                  Tri3Subdivision * newelem = new Tri3Subdivision();

                  // add the first new ghost element to the list just as in the non-corner case
                  if (l == 0)
                    ghost_elems.push_back(newelem);

                  newelem->set_node(0) = nelem->get_node(next[k]);
                  newelem->set_node(1) = nelem->get_node(k);
                  newelem->set_node(2) = node;
                  newelem->set_neighbor(0, nelem);
                  newelem->set_ghost(true);
                  if (l>0)
                    newelem->set_neighbor(2, libmesh_nullptr);
                  nelem->set_neighbor(k, newelem);

                  mesh.add_elem(newelem);
                  mesh.get_boundary_info().add_node(nelem->get_node(k), 1);
                  mesh.get_boundary_info().add_node(nelem->get_node(next[k]), 1);
                  mesh.get_boundary_info().add_node(nelem->get_node(prev[k]), 1);
                  mesh.get_boundary_info().add_node(node, 1);

                  nelem = newelem;
                  k = 2 ;
                }

              Tri3Subdivision * newelem = new Tri3Subdivision();

              newelem->set_node(0) = elem->get_node(next[i]);
              newelem->set_node(1) = nelem->get_node(2);
              newelem->set_node(2) = elem->get_node(prev[i]);
              newelem->set_neighbor(0, nelem);
              nelem->set_neighbor(2, newelem);
              newelem->set_ghost(true);
              newelem->set_neighbor(2, elem);
              elem->set_neighbor(next[i],newelem);

              mesh.add_elem(newelem);

              break;
            }
        }

      for (unsigned int i = 0; i < elem->n_sides(); ++i)
        {
//.........这里部分代码省略.........