mesh_tools.c

一个用来实现偏微分方程中网格的计算库

  return std::distance(begin, end);
}
   


unsigned int MeshTools::n_p_levels (const MeshBase& mesh)
{
  parallel_only();
  
  unsigned int max_p_level = 0;

  // first my local elements
  MeshBase::const_element_iterator
    el     = mesh.local_elements_begin(),
    end_el = mesh.local_elements_end();

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

  // then any unpartitioned objects
  el     = mesh.unpartitioned_elements_begin();
  end_el = mesh.unpartitioned_elements_end();

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

  Parallel::max(max_p_level);
  return max_p_level + 1;
}



void MeshTools::find_nodal_neighbors(const MeshBase&, const Node& n, 
                                     std::vector<std::vector<const Elem*> >& nodes_to_elem_map, 
                                     std::vector<const Node*>& neighbors)
{
  unsigned int global_id = n.id();
  
  //Iterators to iterate through the elements that include this node
  std::vector<const Elem*>::const_iterator el     = nodes_to_elem_map[global_id].begin();
  std::vector<const Elem*>::const_iterator end_el = nodes_to_elem_map[global_id].end();
  
  unsigned int n_ed=0; // Number of edges on the element
  unsigned int ed=0;   // Current edge
  unsigned int l_n=0;  // Local node number
  unsigned int o_n=0;  // Other node on this edge
  
  //Assume we find a edge... then prove ourselves wrong...
  bool found_edge=true;
  
  Node * node_to_save = NULL;
  
  //Look through the elements that contain this node
  //find the local node id... then find the side that
  //node lives on in the element
  //next, look for the _other_ node on that side
  //That other node is a "nodal_neighbor"... save it
  for(;el != end_el;el++)
  {
    //We only care about active elements...
    if((*el)->active())
    {
      n_ed=(*el)->n_edges();
      
      //Find the local node id
      while(global_id != (*el)->node(l_n++)) { }
      l_n--; //Hmmm... take the last one back off
      
      while(ed<n_ed)
      {
        
        //Find the edge the node is on
        while(found_edge && !(*el)->is_node_on_edge(l_n,ed++))
        {
          //This only happens if all the edges have already been found
          if(ed>=n_ed)
            found_edge=false;
        }
        
        //Did we find one?
        if(found_edge)
        {
          ed--; //Take the last one back off again
          
          //Now find the other node on that edge
          while(!(*el)->is_node_on_edge(o_n++,ed) || global_id==(*el)->node(o_n-1)) { }
          o_n--;
          
          //We've found one!  Save it..
          node_to_save=(*el)->get_node(o_n);
          
          //Search to see if we've already found this one
          std::vector<const Node*>::const_iterator result = std::find(neighbors.begin(),neighbors.end(),node_to_save);
          
          //If we didn't find it and add it to the vector
          if(result == neighbors.end())
            neighbors.push_back(node_to_save);
        }
        
        //Reset to look for another
        o_n=0;
        
        //Keep looking for edges, node may be on more than one edge
        ed++;
      }
      
      //Reset to get ready for the next element
      l_n=ed=0;
      found_edge=true;
    }
  }
}

void MeshTools::find_hanging_nodes_and_parents(const MeshBase& mesh, std::map<unsigned int, std::vector<unsigned int> >& hanging_nodes)
{
  MeshBase::const_element_iterator it  = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end = mesh.active_local_elements_end();
  
  //Loop through all the elements
  for (; it != end; ++it)
  {
    //Save it off for easier access
    const Elem* elem = (*it);
    
    //Right now this only works for quad4's
    //libmesh_assert(elem->type() == libMeshEnums::QUAD4);
    if(elem->type() == libMeshEnums::QUAD4)
    {
      //Loop over the sides looking for sides that have hanging nodes
      //This code is inspired by compute_proj_constraints()
      for (unsigned int s=0; s<elem->n_sides(); s++)
      {
        //If not a boundary node
        if (elem->neighbor(s) != NULL)
        {
          // Get pointers to the element's neighbor.
          const Elem* neigh = elem->neighbor(s);
          
          //Is there a coarser element next to this one?
          if (neigh->level() < elem->level()) 
          {
            const Elem *ancestor = elem;
            while (neigh->level() < ancestor->level())
              ancestor = ancestor->parent();
            unsigned int s_neigh = neigh->which_neighbor_am_i(ancestor);
            libmesh_assert (s_neigh < neigh->n_neighbors());
            
            //Couple of helper uints...
            unsigned int node1=0;
            unsigned int node2=0;
            unsigned int hanging_node=0;
                
            bool found_in_neighbor = false;
                
            //Find the two vertices that make up this side
            while(!elem->is_node_on_side(node1++,s)) { }
            node1--;
                
                //Start looking for the second one with the next node
            node2=node1+1;
            
            //Find the other one
            while(!elem->is_node_on_side(node2++,s)) { }
            node2--;
  
                //Pull out their global ids:
            node1 = elem->node(node1);
            node2 = elem->node(node2);
            
            //Now find which node is present in the neighbor
            //FIXME This assumes a level one rule!
            //The _other_ one is the hanging node
            
            //First look for the first one
            //FIXME could be streamlined a bit
            for(unsigned int n=0;n<neigh->n_sides();n++)
            {
              if(neigh->node(n) == node1)
                found_in_neighbor=true;
            }
                
                
            if(!found_in_neighbor)
              hanging_node=node1;
            else //If it wasn't node1 then it must be node2!
              hanging_node=node2;
            
            //Reset these for reuse
            node1=0;
            node2=0;
            
            //Find the first node that makes up the side in the neighbor (these should be the parent nodes)
            while(!neigh->is_node_on_side(node1++,s_neigh)) { }
            node1--;
                
            node2=node1+1;
            
            //Find the second node...
            while(!neigh->is_node_on_side(node2++,s_neigh)) { }
            node2--;
            
            //Save them if we haven't already found the parents for this one
            if(hanging_nodes[hanging_node].size()<2)
            {
              hanging_nodes[hanging_node].push_back(neigh->node(node1));
              hanging_nodes[hanging_node].push_back(neigh->node(node2));
            }
          }
        }
      }
    }
  }
}



void MeshTools::correct_node_proc_ids
  (MeshBase &mesh,
   LocationMap<Node> &loc_map)
{
  // This function must be run on all processors at once
  parallel_only();

  // We'll need the new_nodes_map to answer other processors'
  // requests.  It should never be empty unless we don't have any
  // nodes.
  libmesh_assert(mesh.nodes_begin() == mesh.nodes_end() ||
                 !loc_map.empty());

  // Fix all nodes' processor ids.  Coarsening may have left us with
  // nodes which are no longer touched by any elements of the same
  // processor id, and for DofMap to work we need to fix that.

  // In the first pass, invalidate processor ids for nodes on active
  // elements.  We avoid touching subactive-only nodes.
  MeshBase::element_iterator       e_it  = mesh.active_elements_begin();
  const MeshBase::element_iterator e_end = mesh.active_elements_end();
  for (; e_it != e_end; ++e_it)
    {
      Elem *elem = *e_it;
      for (unsigned int n=0; n != elem->n_nodes(); ++n)
        {
          Node *node = elem->get_node(n);
          node->invalidate_processor_id();
        }
    }

  // In the second pass, find the lowest processor ids on active
  // elements touching each node, and set the node processor id.
  for (e_it = mesh.active_elements_begin(); e_it != e_end; ++e_it)
    {
      Elem *elem = *e_it;
      unsigned int proc_id = elem->processor_id();
      for (unsigned int n=0; n != elem->n_nodes(); ++n)
        {
          Node *node = elem->get_node(n);
          if (node->processor_id() == DofObject::invalid_processor_id ||
              node->processor_id() > proc_id)
            node->processor_id() = proc_id;
        }
    }

  // Those two passes will correct every node that touches a local
  // element, but we can't be sure about nodes touching remote
  // elements.  Fix those now.
  MeshCommunication().make_node_proc_ids_parallel_consistent
    (mesh, loc_map);
}



#ifdef DEBUG
void MeshTools::libmesh_assert_valid_node_pointers(const MeshBase &mesh)
{
  const MeshBase::const_element_iterator el_end =
    mesh.elements_end();
  for (MeshBase::const_element_iterator el = 
       mesh.elements_begin(); el != el_end; ++el)
    {
      const Elem* elem = *el;
      libmesh_assert (elem);
      while (elem)
	{
	  elem->libmesh_assert_valid_node_pointers();
          for (unsigned int n=0; n != elem->n_neighbors(); ++n)
            if (elem->neighbor(n) &&
                elem->neighbor(n) != remote_elem)
	      elem->neighbor(n)->libmesh_assert_valid_node_pointers();
	    
          libmesh_assert (elem->parent() != remote_elem);
	  elem = elem->parent();
	}
    }
}


void MeshTools::libmesh_assert_valid_remote_elems(const MeshBase &mesh)
{
  const MeshBase::const_element_iterator el_end =
    mesh.local_elements_end();
  for (MeshBase::const_element_iterator el = 
       mesh.local_elements_begin(); el != el_end; ++el)
    {
      const Elem* elem = *el;
      libmesh_assert (elem);
      for (unsigned int n=0; n != elem->n_neighbors(); ++n)
	libmesh_assert (elem->neighbor(n) != remote_elem);
#ifdef ENABLE_AMR
      const Elem* parent = elem->parent();
      if (parent)
	{
          libmesh_assert (parent != remote_elem);
          for (unsigned int c=0; c != elem->n_children(); ++c)
	    libmesh_assert (parent->child(c) != remote_elem);
	}
#endif
    }
}


void MeshTools::libmesh_assert_no_links_to_elem(const MeshBase &mesh,
                                                const Elem *bad_elem)
{
  const MeshBase::const_element_iterator el_end =
    mesh.elements_end();
  for (MeshBase::const_element_iterator el = 
       mesh.elements_begin(); el != el_end; ++el)
    {
      const Elem* elem = *el;
      libmesh_assert (elem);
      libmesh_assert (elem->parent() != bad_elem);
      for (unsigned int n=0; n != elem->n_neighbors(); ++n)
	libmesh_assert(elem->neighbor(n) != bad_elem);
#ifdef ENABLE_AMR
      if (elem->has_children())
        for (unsigned int c=0; c != elem->n_children(); ++c)
	  libmesh_assert(elem->child(c) != bad_elem);
#endif
    }
}



void MeshTools::libmesh_assert_valid_elem_ids(const MeshBase &mesh)
{
  unsigned int lastprocid = 0;
  unsigned int lastelemid = 0;

  const MeshBase::const_element_iterator el_end =
    mesh.active_elements_end();
  for (MeshBase::const_element_iterator el = 
       mesh.active_elements_begin(); el != el_end; ++el)
    {
      const Elem* elem = *el;
      libmesh_assert (elem);
      unsigned int elemprocid = elem->processor_id();
      unsigned int elemid = elem->id();

      libmesh_assert(elemid >= lastelemid);
      libmesh_assert(elemprocid >= lastprocid);

      lastelemid = elemid;
      lastprocid = elemprocid;
    }
}



void MeshTools::libmesh_assert_valid_node_procids(const MeshBase &mesh)
{
  parallel_only();
  if (libMesh::n_processors() == 1)
    return;

  libmesh_assert(Parallel::verify(mesh.max_node_id()));
  std::vector<bool> node_touched_by_me(mesh.max_node_id(), false);

  const MeshBase::const_element_iterator el_end =
    mesh.active_local_elements_end();
  for (MeshBase::const_element_iterator el = 
       mesh.active_local_elements_begin(); el != el_end; ++el)
    {
      const Elem* elem = *el;
      libmesh_assert (elem);

      for (unsigned int i=0; i != elem->n_nodes(); ++i)
        {
          const Node *node = elem->get_node(i);
          unsigned int nodeid = node->id();
          node_touched_by_me[nodeid] = true;
        }
    }
  std::vector<bool> node_touched_by_anyone(node_touched_by_me);
  Parallel::max(node_touched_by_anyone);

  const MeshBase::const_node_iterator nd_end = mesh.local_nodes_end();
  for (MeshBase::const_node_iterator nd = mesh.local_nodes_begin();
       nd != nd_end; ++nd)
    {
      const Node *node = *nd;
      libmesh_assert(node);

      unsigned int nodeid = node->id();
      libmesh_assert(!node_touched_by_anyone[nodeid] ||
                     node_touched_by_me[nodeid]);
    }
}



#ifdef ENABLE_AMR
void MeshTools::libmesh_assert_valid_refinement_flags(const MeshBase &mesh)
{
  parallel_only();
  if (libMesh::n_processors() == 1)
    return;

  std::vector<unsigned char> my_elem_h_state(mesh.max_elem_id(), 255);
  std::vector<unsigned char> my_elem_p_state(mesh.max_elem_id(), 255);

  const MeshBase::const_element_iterator el_end =
    mesh.elements_end();
  for (MeshBase::const_element_iterator el = 
       mesh.elements_begin(); el != el_end; ++el)
    {
      const Elem* elem = *el;
      libmesh_assert (elem);
      unsigned int elemid = elem->id();

      my_elem_h_state[elemid] =
        static_cast<unsigned char>(elem->refinement_flag());

      my_elem_p_state[elemid] =
        static_cast<unsigned char>(elem->p_refinement_flag());
    }
  std::vector<unsigned char> min_elem_h_state(my_elem_h_state);
  Parallel::min(min_elem_h_state);

  std::vector<unsigned char> min_elem_p_state(my_elem_p_state);
  Parallel::min(min_elem_p_state);

  for (unsigned int i=0; i!= mesh.max_elem_id(); ++i)
    {
      libmesh_assert(my_elem_h_state[i] == 255 ||
                     my_elem_h_state[i] == min_elem_h_state[i]);
      libmesh_assert(my_elem_p_state[i] == 255 ||
                     my_elem_p_state[i] == min_elem_p_state[i]);
    }
}
#else
void MeshTools::libmesh_assert_valid_refinement_flags(const MeshBase &)
{
}
#endif // ENABLE_AMR



void MeshTools::libmesh_assert_valid_neighbors(const MeshBase &mesh)
{
  const MeshBase::const_element_iterator el_end = mesh.elements_end();
  for (MeshBase::const_element_iterator el = mesh.elements_begin();
       el != el_end; ++el)
    {
      const Elem* elem = *el;
      libmesh_assert (elem);
      elem->libmesh_assert_valid_neighbors();
    }
}
#endif



void MeshTools::Private::globally_renumber_nodes_and_elements (MeshBase& mesh)
{
  MeshCommunication().assign_global_indices(mesh);
}



void MeshTools::Private::fix_broken_node_and_element_numbering (SerialMesh &mesh)
{
   // Nodes first
  for (unsigned int n=0; n<mesh._nodes.size(); n++)
    if (mesh._nodes[n] != NULL)
      mesh._nodes[n]->set_id() = n;

  // Elements next
  for (unsigned int e=0; e<mesh._elements.size(); e++)
    if (mesh._elements[e] != NULL)
      mesh._elements[e]->set_id() = e; 
}



void MeshTools::Private::fix_broken_node_and_element_numbering (ParallelMesh &mesh)
{
  // We need access to iterators for the underlying containers,
  // not the mapvector<> reimplementations.
  mapvector<Node*>::maptype &nodes = mesh._nodes;
  mapvector<Elem*>::maptype &elem  = mesh._elements;
  
  // Nodes first
  {
    mapvector<Node*>::maptype::iterator
      it  = nodes.begin(),
      end = nodes.end();

    for (; it != end; ++it)
      if (it->second != NULL)
	it->second->set_id() = it->first;
  }

  // Elements next
  {
    mapvector<Elem*>::maptype::iterator
      it  = elem.begin(),
      end = elem.end();

    for (; it != end; ++it)
      if (it->second != NULL)
	it->second->set_id() = it->first;
  }
}