C++ MeshBase类代码示例

void UCDIO::write_interior_elems(std::ostream& out_stream, const MeshBase& mesh)
{
  std::string type[] =
    { "edge",  "edge",  "edge",
      "tri",   "tri",
      "quad",  "quad",  "quad",
      "tet",   "tet",
      "hex",   "hex",   "hex",
      "prism", "prism", "prism",
      "pyramid" };

  MeshBase::const_element_iterator it  = mesh.elements_begin();
  const MeshBase::const_element_iterator end = mesh.elements_end();

  unsigned int e=1; // 1-based element number for UCD

  // Write element information
  for (; it != end; ++it)
    {
      libmesh_assert (out_stream.good());

      // PB: I believe these are the only supported ElemTypes.
      const ElemType etype = (*it)->type();
      if( (etype != TRI3) && (etype != QUAD4) &&
          (etype != TET4) && (etype != HEX8) &&
          (etype != PRISM6) && (etype != PYRAMID5) )
        libmesh_error_msg("Error: Unsupported ElemType for UCDIO.");

      out_stream << e++ << " 0 " << type[etype] << "\t";
      // (*it)->write_ucd_connectivity(out_stream);
      (*it)->write_connectivity(out_stream, UCD);
    }

  return;
}

void TopologyMap::fill(const MeshBase& mesh)
{
  // Populate the nodes map
  MeshBase::const_element_iterator
    it = mesh.elements_begin(),
    end = mesh.elements_end();
  for (; it != end; ++it)
    {
      const Elem* elem = *it;

      // We only need to add nodes which might be added during mesh
      // refinement; this means they need to be child nodes.
      if (!elem->has_children())
        continue;

      for (unsigned int c = 0; c != elem->n_children(); ++c)
        {
          if (elem->child(c)->is_remote())
            continue;

          for (unsigned int n = 0; n != elem->n_nodes_in_child(c); ++n)
            {
              const std::vector<std::pair<dof_id_type, dof_id_type> >
                bracketing_nodes = elem->bracketing_nodes(c,n);

              this->add_node(*elem->child(c)->get_node(n),
                             bracketing_nodes);
            }
        }
    }
}

void MeshTools::Subdivision::tag_boundary_ghosts(MeshBase & mesh)
{
  MeshBase::element_iterator       el     = mesh.elements_begin();
  const MeshBase::element_iterator end_el = mesh.elements_end();
  for (; el != end_el; ++el)
    {
      Elem * elem = *el;
      libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);

      Tri3Subdivision * sd_elem = static_cast<Tri3Subdivision *>(elem);
      for (unsigned int i = 0; i < elem->n_sides(); ++i)
        {
          if (elem->neighbor(i) == libmesh_nullptr)
            {
              sd_elem->set_ghost(true);
              // set all other neighbors to ghosts as well
              if (elem->neighbor(next[i]))
                {
                  Tri3Subdivision * nb = static_cast<Tri3Subdivision *>(elem->neighbor(next[i]));
                  nb->set_ghost(true);
                }
              if (elem->neighbor(prev[i]))
                {
                  Tri3Subdivision * nb = static_cast<Tri3Subdivision *>(elem->neighbor(prev[i]));
                  nb->set_ghost(true);
                }
            }
        }
    }
}

void findElementsIntersectedByPlane(const Plane & plane, const MeshBase & mesh, std::vector<const Elem *> & intersected_elems)
{
  // Loop over all elements to find elements intersected by the plane
  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;
    bool intersected = false;

    // Check whether the first node of this element is below or above the plane
    const Node node0 = *elem->get_node(0);
    bool node0_above_plane = plane.above_surface(node0);

    // Loop over the rest of the nodes and check if any node is on the other side of the plane
    for (unsigned int i = 1; i < elem->n_nodes(); ++i)
    {
      const Node node = *elem->get_node(i);

      bool node_above_plane = plane.above_surface(node);
      if (node0_above_plane != node_above_plane)
        intersected = true;
    }

    if (intersected)
      intersected_elems.push_back(elem);
  }
}

void SubdomainPartitioner::_do_partition (MeshBase & mesh,
                                          const unsigned int n)
{
  libmesh_assert_greater (n, 0);

  // Check for an easy return.  If the user has not specified any
  // entries in the chunks vector, we just do a single partitioning.
  if ((n == 1) || chunks.empty())
    {
      this->single_partition (mesh);
      return;
    }

  // Now actually do the partitioning.
  LOG_SCOPE ("_do_partition()", "SubdomainPartitioner");

  // For each chunk, construct an iterator range for the set of
  // subdomains in question, and pass it to the internal Partitioner.
  for (std::size_t c=0; c<chunks.size(); ++c)
    {
      MeshBase::element_iterator
        it = mesh.active_subdomain_set_elements_begin(chunks[c]),
        end = mesh.active_subdomain_set_elements_end(chunks[c]);

      _internal_partitioner->partition_range(mesh, it, end, n);
    }
}

void Partitioner::repartition (MeshBase & mesh,
                               const unsigned int n)
{
  // we cannot partition into more pieces than we have
  // active elements!
  const unsigned int n_parts =
    static_cast<unsigned int>
    (std::min(mesh.n_active_elem(), static_cast<dof_id_type>(n)));

  // Set the number of partitions in the mesh
  mesh.set_n_partitions()=n_parts;

  if (n_parts == 1)
    {
      this->single_partition (mesh);
      return;
    }

  // First assign a temporary partitioning to any unpartitioned elements
  Partitioner::partition_unpartitioned_elements(mesh, n_parts);

  // Call the partitioning function
  this->_do_repartition(mesh,n_parts);

  // Set the parent's processor ids
  Partitioner::set_parent_processor_ids(mesh);

  // Set the node's processor ids
  Partitioner::set_node_processor_ids(mesh);
}

void SFCPartitioner::_do_partition (MeshBase & mesh,
                                    const unsigned int n)
{
  this->partition_range(mesh,
                        mesh.active_elements_begin(),
                        mesh.active_elements_end(),
                        n);
}

void LocationMap<Elem>::fill(MeshBase& mesh)
{
  // Populate the elem map
  MeshBase::element_iterator       it  = mesh.active_elements_begin(),
                                   end = mesh.active_elements_end();
  for (; it != end; ++it)
    this->insert(**it);
}

void MetisPartitioner::_do_partition (MeshBase & mesh,
                                      const unsigned int n_pieces)
{
  this->partition_range(mesh,
                        mesh.active_elements_begin(),
                        mesh.active_elements_end(),
                        n_pieces);
}

void MeshBase::addAttribs(const MeshBase& other)
{
    int num = other.numAttribs();
    for (int i = 0; i < num; i++)
    {
        const AttribSpec& spec = other.attribSpec(i);
        addAttrib(spec.type, spec.format, spec.length);
    }
}

void Partitioner::partition (MeshBase & mesh,
                             const unsigned int n)
{
  libmesh_parallel_only(mesh.comm());

  // BSK - temporary fix while redistribution is integrated 6/26/2008
  // Uncomment this to not repartition in parallel
  //   if (!mesh.is_serial())
  //     return;

  // we cannot partition into more pieces than we have
  // active elements!
  const unsigned int n_parts =
    static_cast<unsigned int>
    (std::min(mesh.n_active_elem(), static_cast<dof_id_type>(n)));

  // Set the number of partitions in the mesh
  mesh.set_n_partitions()=n_parts;

  if (n_parts == 1)
    {
      this->single_partition (mesh);
      return;
    }

  // First assign a temporary partitioning to any unpartitioned elements
  Partitioner::partition_unpartitioned_elements(mesh, n_parts);

  // Call the partitioning function
  this->_do_partition(mesh,n_parts);

  // Set the parent's processor ids
  Partitioner::set_parent_processor_ids(mesh);

  // Redistribute elements if necessary, before setting node processor
  // ids, to make sure those will be set consistently
  mesh.redistribute();

#ifdef DEBUG
  MeshTools::libmesh_assert_valid_remote_elems(mesh);

  // Messed up elem processor_id()s can leave us without the child
  // elements we need to restrict vectors on a distributed mesh
  MeshTools::libmesh_assert_valid_procids<Elem>(mesh);
#endif

  // Set the node's processor ids
  Partitioner::set_node_processor_ids(mesh);

#ifdef DEBUG
  MeshTools::libmesh_assert_valid_procids<Elem>(mesh);
#endif

  // Give derived Mesh classes a chance to update any cached data to
  // reflect the new partitioning
  mesh.update_post_partitioning();
}

void MeshData::translate (const MeshBase & out_mesh,
                          std::vector<Number> & values,
                          std::vector<std::string> & names) const
{
  libmesh_assert (_active || _compatibility_mode);

  START_LOG("translate()", "MeshData");

  const unsigned int n_comp = this->n_val_per_node();

  // transfer our nodal data to a vector
  // that may be written concurrently
  // with the \p out_mesh.
  {
    // reserve memory for the nodal data
    values.reserve(n_comp*out_mesh.n_nodes());

    // iterate over the mesh's nodes
    MeshBase::const_node_iterator       nodes_it  = out_mesh.nodes_begin();
    const MeshBase::const_node_iterator nodes_end = out_mesh.nodes_end();

    // Do not use the \p get_data() method, but the operator()
    // method, since this returns by default a zero value,
    // when there is no nodal data.
    for (; nodes_it != nodes_end; ++nodes_it)
      {
        const Node * node = *nodes_it;

        for (unsigned int c= 0; c<n_comp; c++)
          values.push_back(this->operator()(node, c));
      }
  }



  // Now we have the data, nicely stored in \p values.
  // It remains to give names to the data, then write to
  // file.
  {
    names.reserve(n_comp);

    // this naming scheme only works up to n_comp=100
    // (at least for gmv-accepted variable names)
    libmesh_assert_less (n_comp, 100);

    for (unsigned int n=0; n<n_comp; n++)
      {
        std::ostringstream name_buf;
        name_buf << "bc_" << n;
        names.push_back(name_buf.str());
      }
  }

  STOP_LOG("translate()", "MeshData");
}

void
MoosePreconditioner::copyVarValues(MeshBase & mesh,
                                   const unsigned int from_system,
                                   const unsigned int from_var,
                                   const NumericVector<Number> & from_vector,
                                   const unsigned int to_system,
                                   const unsigned int to_var,
                                   NumericVector<Number> & to_vector)
{
  {
    MeshBase::node_iterator it = mesh.local_nodes_begin();
    MeshBase::node_iterator it_end = mesh.local_nodes_end();

    for (; it != it_end; ++it)
    {
      Node * node = *it;

      unsigned int n_comp = node->n_comp(from_system, from_var);

      mooseAssert(node->n_comp(from_system, from_var) == node->n_comp(to_system, to_var),
                  "Number of components does not match in each system");

      for (unsigned int i = 0; i < n_comp; i++)
      {
        dof_id_type from_dof = node->dof_number(from_system, from_var, i);
        dof_id_type to_dof = node->dof_number(to_system, to_var, i);

        to_vector.set(to_dof, from_vector(from_dof));
      }
    }
  }
  {
    MeshBase::element_iterator it = mesh.local_elements_begin();
    MeshBase::element_iterator it_end = mesh.local_elements_end();

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

      unsigned int n_comp = elem->n_comp(from_system, from_var);

      mooseAssert(elem->n_comp(from_system, from_var) == elem->n_comp(to_system, to_var),
                  "Number of components does not match in each system");

      for (unsigned int i = 0; i < n_comp; i++)
      {
        dof_id_type from_dof = elem->dof_number(from_system, from_var, i);
        dof_id_type to_dof = elem->dof_number(to_system, to_var, i);

        to_vector.set(to_dof, from_vector(from_dof));
      }
    }
  }
}

unsigned int CheckpointIO::n_active_levels_on_processor(const MeshBase & mesh) const
{
  unsigned int max_level = 0;

  MeshBase::const_element_iterator el = mesh.active_elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.active_elements_end();

  for( ; el != end_el; ++el)
    max_level = std::max((*el)->level(), max_level);

  return max_level + 1;
}

void integrate_function (const MeshBase &mesh)
{
#if defined(LIBMESH_HAVE_TRIANGLE) && defined(LIBMESH_HAVE_TETGEN)
    MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();

    std::vector<Real> vertex_distance;

    QComposite<QGauss> qrule (mesh.mesh_dimension(), FIRST);
    //QGauss qrule (mesh.mesh_dimension(), FIRST);

    UniquePtr<FEBase> fe (FEBase::build (mesh.mesh_dimension(), FEType (FIRST, LAGRANGE)));

    Real int_val=0.;

    const std::vector<Point> &q_points = fe->get_xyz();
    const std::vector<Real>  &JxW      = fe->get_JxW();

    for (; el!=end_el; ++el)
    {
        const Elem *elem = *el;

        vertex_distance.clear();

        for (unsigned int v=0; v<elem->n_vertices(); v++)
            vertex_distance.push_back (distance(elem->point(v)));

        qrule.init (*elem, vertex_distance);

        fe->reinit (elem,
                    &(qrule.get_points()),
                    &(qrule.get_weights()));


        // TODO:  would it be valuable to have the composite quadrature rule sort
        // from smallest to largest JxW value to help prevent
        // ... large + small + large + large + small ...
        // type truncation errors?
        for (unsigned int qp=0; qp<q_points.size(); qp++)
            int_val += JxW[qp] * integrand(q_points[qp]);
    }

    mesh.comm().sum (int_val);

    std::cout  << "\n***********************************\n"
               << " int_val   = " << int_val << std::endl
               << " exact_val = " <<  1*(2*2 - radius*radius*pi) + 10.*(radius*radius*pi)
               << "\n***********************************\n"
               << std::endl;
#else
    libmesh_ignore(mesh);
#endif
}