Home
Publications
Presentations
Support
Download
Examples
Developers
Class Docs
Namespaces
Classes
Files
Mailing Lists
Gallery
Wiki

mesh_tools.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



// C++ includes
#include <limits>
#include <numeric> // for std::accumulate
#include <set>

// Local includes
#include "libmesh/elem.h"
#include "libmesh/elem_range.h"
#include "libmesh/location_maps.h"
#include "libmesh/mesh_base.h"
#include "libmesh/mesh_communication.h"
#include "libmesh/mesh_tools.h"
#include "libmesh/node_range.h"
#include "libmesh/parallel.h"
#include "libmesh/parallel_mesh.h"
#include "libmesh/serial_mesh.h"
#include "libmesh/sphere.h"
#include "libmesh/threads.h"

#ifdef DEBUG
#  include "libmesh/remote_elem.h"
#endif



// ------------------------------------------------------------
// anonymous namespace for helper classes
namespace {

  using namespace libMesh;

  class SumElemWeight
  {
  public:
    SumElemWeight () :
      _weight(0)
    {}

    SumElemWeight (SumElemWeight &, Threads::split) :
      _weight(0)
    {}

    void operator()(const ConstElemRange &range)
    {
      for (ConstElemRange::const_iterator it = range.begin(); it !=range.end(); ++it)
        _weight += (*it)->n_nodes();
    }

    dof_id_type weight() const
    { return _weight; }

// If we don't have threads we never need a join, and icpc yells a
// warning if it sees an anonymous function that's never used
#ifdef LIBMESH_HAVE_TBB_API
    void join (const SumElemWeight &other)
    { _weight += other.weight(); }
#endif

  private:
    dof_id_type _weight;
  };


  class FindBBox
  {
  public:
    FindBBox () :
      _vmin(LIBMESH_DIM,  std::numeric_limits<Real>::max()),
      _vmax(LIBMESH_DIM, -std::numeric_limits<Real>::max())
    {}

    FindBBox (FindBBox &other, Threads::split) :
      _vmin(other._vmin),
      _vmax(other._vmax)
    {}

    std::vector<Real> & min() { return _vmin; }
    std::vector<Real> & max() { return _vmax; }

    void operator()(const ConstNodeRange &range)
    {
      for (ConstNodeRange::const_iterator it = range.begin(); it != range.end(); ++it)
        {
          const Node *node = *it;
          libmesh_assert(node);

          for (unsigned int i=0; i<LIBMESH_DIM; i++)
            {
              _vmin[i] = std::min(_vmin[i], (*node)(i));
              _vmax[i] = std::max(_vmax[i], (*node)(i));
            }
        }
    }

    void operator()(const ConstElemRange &range)
    {
      for (ConstElemRange::const_iterator it = range.begin(); it != range.end(); ++it)
        {
          const Elem *elem = *it;
          libmesh_assert(elem);

          for (unsigned int n=0; n<elem->n_nodes(); n++)
            {
              const Point &point = elem->point(n);

              for (unsigned int i=0; i<LIBMESH_DIM; i++)
                {
                  _vmin[i] = std::min(_vmin[i], point(i));
                  _vmax[i] = std::max(_vmax[i], point(i));
                }
            }
        }
    }

// If we don't have threads we never need a join, and icpc yells a
// warning if it sees an anonymous function that's never used
#ifdef LIBMESH_HAVE_TBB_API
    void join (const FindBBox &other)
    {
      for (unsigned int i=0; i<LIBMESH_DIM; i++)
        {
          _vmin[i] = std::min(_vmin[i], other._vmin[i]);
          _vmax[i] = std::max(_vmax[i], other._vmax[i]);
        }
    }
#endif

    MeshTools::BoundingBox bbox () const
    {
      Point pmin(_vmin[0]
#if LIBMESH_DIM > 1
                 , _vmin[1]
#endif
#if LIBMESH_DIM > 2
                 , _vmin[2]
#endif
                );
      Point pmax(_vmax[0]
#if LIBMESH_DIM > 1
                 , _vmax[1]
#endif
#if LIBMESH_DIM > 2
                 , _vmax[2]
#endif
                );

      const MeshTools::BoundingBox ret_val(pmin, pmax);

      return ret_val;
    }

  private:
    std::vector<Real> _vmin;
    std::vector<Real> _vmax;
  };

#ifdef DEBUG
  void assert_semiverify_dofobj(const DofObject *d)
    {
      if (d)
        {
          const unsigned int n_sys = d->n_systems();
         
          std::vector<unsigned int> n_vars (n_sys, 0);
          for (unsigned int s = 0; s != n_sys; ++s)
            n_vars[s] = d->n_vars(s);

          const unsigned int tot_n_vars =
            std::accumulate(n_vars.begin(), n_vars.end(), 0);

          std::vector<unsigned int> n_comp (tot_n_vars, 0);
          std::vector<dof_id_type> first_dof (tot_n_vars, 0);
         
          for (unsigned int s = 0, i=0; s != n_sys; ++s)
            for (unsigned int v = 0; v != n_vars[s]; ++v, ++i)
              {
                n_comp[i] = d->n_comp(s,v);
                first_dof[i] = n_comp[i] ? d->dof_number(s,v,0) : libMesh::invalid_uint;
              }

          libmesh_assert(CommWorld.semiverify(&n_sys));
          libmesh_assert(CommWorld.semiverify(&n_vars));
          libmesh_assert(CommWorld.semiverify(&n_comp));
          libmesh_assert(CommWorld.semiverify(&first_dof));
        }
      else
        {
          const unsigned int* p_ui = NULL;
          const std::vector<unsigned int>* p_vui = NULL;

          libmesh_assert(CommWorld.semiverify(p_ui));
          libmesh_assert(CommWorld.semiverify(p_vui));
          libmesh_assert(CommWorld.semiverify(p_vui));
          libmesh_assert(CommWorld.semiverify(p_vui));
        }
    }
#endif // DEBUG

}


namespace libMesh
{
// Small helper function to make intersect more readable.
bool is_between(Real min, Real check, Real max)
{
  return min <= check && check <= max;
}

bool MeshTools::BoundingBox::intersect (const BoundingBox & other_box) const
{
  // Make local variables first to make thiings more clear in a moment
  const Real& my_min_x = this->first(0);
  const Real& my_max_x = this->second(0);
  const Real& other_min_x = other_box.first(0);
  const Real& other_max_x = other_box.second(0);

  const bool x_int = is_between(my_min_x, other_min_x, my_max_x) || is_between(my_min_x, other_max_x, my_max_x) ||
                     is_between(other_min_x, my_min_x, other_max_x) || is_between(other_min_x, my_max_x, other_max_x);

  bool intersection_true = x_int;

#if LIBMESH_DIM > 1
  const Real& my_min_y = this->first(1);
  const Real& my_max_y = this->second(1);
  const Real& other_min_y = other_box.first(1);
  const Real& other_max_y = other_box.second(1);

  const bool y_int = is_between(my_min_y, other_min_y, my_max_y) || is_between(my_min_y, other_max_y, my_max_y) ||
                     is_between(other_min_y, my_min_y, other_max_y) || is_between(other_min_y, my_max_y, other_max_y);

  intersection_true = intersection_true && y_int;
#endif

#if LIBMESH_DIM > 2
  const Real& my_min_z = this->first(2);
  const Real& my_max_z = this->second(2);
  const Real& other_min_z = other_box.first(2);
  const Real& other_max_z = other_box.second(2);

  const bool z_int = is_between(my_min_z, other_min_z, my_max_z) || is_between(my_min_z, other_max_z, my_max_z) ||
                     is_between(other_min_z, my_min_z, other_max_z) || is_between(other_min_z, my_max_z, other_max_z);

  intersection_true = intersection_true && z_int;
#endif

  return intersection_true;
}

bool MeshTools::BoundingBox::contains_point (const Point & p) const
{
  // Make local variables first to make thiings more clear in a moment
  Real my_min_x = this->first(0);
  Real my_max_x = this->second(0);
  bool x_int = is_between(my_min_x, p(0), my_max_x);

  bool intersection_true = x_int;

#if LIBMESH_DIM > 1
  Real my_min_y = this->first(1);
  Real my_max_y = this->second(1);
  bool y_int = is_between(my_min_y, p(1), my_max_y);

  intersection_true = intersection_true && y_int;
#endif


#if LIBMESH_DIM > 2
  Real my_min_z = this->first(2);
  Real my_max_z = this->second(2);
  bool z_int = is_between(my_min_z, p(2), my_max_z);

  intersection_true = intersection_true && z_int;
#endif

  return intersection_true;
}

// ------------------------------------------------------------
// MeshTools functions
dof_id_type MeshTools::total_weight(const MeshBase& mesh)
{
  if (!mesh.is_serial())
    {
      parallel_only();
      dof_id_type weight = MeshTools::weight (mesh, libMesh::processor_id());
      CommWorld.sum(weight);
      dof_id_type unpartitioned_weight =
        MeshTools::weight (mesh, DofObject::invalid_processor_id);
      return weight + unpartitioned_weight;
    }

  SumElemWeight sew;

  Threads::parallel_reduce (ConstElemRange (mesh.elements_begin(),
                                            mesh.elements_end()),
                            sew);
  return sew.weight();

}



dof_id_type MeshTools::weight(const MeshBase& mesh, const processor_id_type pid)
{
  SumElemWeight sew;

  Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(pid),
                                            mesh.pid_elements_end(pid)),
                            sew);
  return sew.weight();
}



void MeshTools::build_nodes_to_elem_map (const MeshBase& mesh,
                                         std::vector<std::vector<dof_id_type> >& nodes_to_elem_map)
{
  nodes_to_elem_map.resize (mesh.n_nodes());

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

  for (; el != end; ++el)
    for (unsigned int n=0; n<(*el)->n_nodes(); n++)
      {
        libmesh_assert_less ((*el)->node(n), nodes_to_elem_map.size());
        libmesh_assert_less ((*el)->id(), mesh.n_elem());

        nodes_to_elem_map[(*el)->node(n)].push_back((*el)->id());
      }
}



void MeshTools::build_nodes_to_elem_map (const MeshBase& mesh,
                                         std::vector<std::vector<const Elem*> >& nodes_to_elem_map)
{
  nodes_to_elem_map.resize (mesh.n_nodes());

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

  for (; el != end; ++el)
    for (unsigned int n=0; n<(*el)->n_nodes(); n++)
      {
        libmesh_assert_less ((*el)->node(n), nodes_to_elem_map.size());

        nodes_to_elem_map[(*el)->node(n)].push_back(*el);
      }
}



void MeshTools::find_boundary_nodes (const MeshBase& mesh,
                                     std::vector<bool>& on_boundary)
{
  // Resize the vector which holds boundary nodes and fill with false.
  on_boundary.resize(mesh.n_nodes());
  std::fill(on_boundary.begin(),
            on_boundary.end(),
            false);

  // Loop over elements, find those on boundary, and
  // mark them as true in on_boundary.
  MeshBase::const_element_iterator       el  = mesh.active_elements_begin();
  const MeshBase::const_element_iterator end = mesh.active_elements_end();

  for (; el != end; ++el)
    for (unsigned int s=0; s<(*el)->n_neighbors(); s++)
      if ((*el)->neighbor(s) == NULL) // on the boundary
        {
          const AutoPtr<Elem> side((*el)->build_side(s));

          for (unsigned int n=0; n<side->n_nodes(); n++)
            on_boundary[side->node(n)] = true;
        }
}



MeshTools::BoundingBox
MeshTools::bounding_box(const MeshBase& mesh)
{
  // This function must be run on all processors at once
  parallel_only();

  FindBBox find_bbox;

  Threads::parallel_reduce (ConstNodeRange (mesh.local_nodes_begin(),
                                            mesh.local_nodes_end()),
                            find_bbox);

  // and the unpartitioned nodes
  Threads::parallel_reduce (ConstNodeRange (mesh.pid_nodes_begin(DofObject::invalid_processor_id),
                                            mesh.pid_nodes_end(DofObject::invalid_processor_id)),
                            find_bbox);

  // Compare the bounding boxes across processors
  CommWorld.min(find_bbox.min());
  CommWorld.max(find_bbox.max());

  return find_bbox.bbox();
}



Sphere
MeshTools::bounding_sphere(const MeshBase& mesh)
{
  BoundingBox bbox = bounding_box(mesh);

  const Real  diag = (bbox.second - bbox.first).size();
  const Point cent = (bbox.second + bbox.first)/2;

  return Sphere (cent, .5*diag);
}



MeshTools::BoundingBox
MeshTools::processor_bounding_box (const MeshBase& mesh,
                                   const processor_id_type pid)
{
  libmesh_assert_less (pid, libMesh::n_processors());

  FindBBox find_bbox;

  Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(pid),
                                            mesh.pid_elements_end(pid)),
                            find_bbox);

  return find_bbox.bbox();
}



Sphere
MeshTools::processor_bounding_sphere (const MeshBase& mesh,
                                      const processor_id_type pid)
{
  BoundingBox bbox = processor_bounding_box(mesh,pid);

  const Real  diag = (bbox.second - bbox.first).size();
  const Point cent = (bbox.second + bbox.first)/2;

  return Sphere (cent, .5*diag);
}



MeshTools::BoundingBox
MeshTools::subdomain_bounding_box (const MeshBase& mesh,
                                   const subdomain_id_type sid)
{
  libmesh_assert_not_equal_to (mesh.n_nodes(), 0);

  Point min( 1.e30,  1.e30,  1.e30);
  Point max(-1.e30, -1.e30, -1.e30);

  for (unsigned int e=0; e<mesh.n_elem(); e++)
    if (mesh.elem(e)->subdomain_id() == sid)
      for (unsigned int n=0; n<mesh.elem(e)->n_nodes(); n++)
        for (unsigned int i=0; i<mesh.spatial_dimension(); i++)
          {
            min(i) = std::min(min(i), mesh.point(mesh.elem(e)->node(n))(i));
            max(i) = std::max(max(i), mesh.point(mesh.elem(e)->node(n))(i));
          }

  return BoundingBox (min, max);
}



Sphere
MeshTools::subdomain_bounding_sphere (const MeshBase& mesh,
                                      const subdomain_id_type sid)
{
  BoundingBox bbox = subdomain_bounding_box(mesh,sid);

  const Real  diag = (bbox.second - bbox.first).size();
  const Point cent = (bbox.second + bbox.first)/2;

  return Sphere (cent, .5*diag);
}



void MeshTools::elem_types (const MeshBase& mesh,
                            std::vector<ElemType>& et)
{
  MeshBase::const_element_iterator       el  = mesh.elements_begin();
  const MeshBase::const_element_iterator end = mesh.elements_end();

  // Automatically get the first type
  et.push_back((*el)->type());  ++el;

  // Loop over the rest of the elements.
  // If the current element type isn't in the
  // vector, insert it.
  for (; el != end; ++el)
    if (!std::count(et.begin(), et.end(), (*el)->type()))
      et.push_back((*el)->type());
}



dof_id_type MeshTools::n_elem_of_type (const MeshBase& mesh,
                                       const ElemType type)
{
  return static_cast<dof_id_type>(std::distance(mesh.type_elements_begin(type),
                                                mesh.type_elements_end  (type)));
}



dof_id_type MeshTools::n_active_elem_of_type (const MeshBase& mesh,
                                              const ElemType type)
{
  return static_cast<dof_id_type>(std::distance(mesh.active_type_elements_begin(type),
                                                mesh.active_type_elements_end  (type)));
}

dof_id_type MeshTools::n_non_subactive_elem_of_type_at_level(const MeshBase& mesh,
                                                             const ElemType type,
                                                             const unsigned int level)
{
  dof_id_type cnt = 0;
  // iterate over the elements of the specified type
  MeshBase::const_element_iterator el = mesh.type_elements_begin(type);
  const MeshBase::const_element_iterator end = mesh.type_elements_end(type);

  for(; el!=end; ++el)
    if( ((*el)->level() == level) && !(*el)->subactive())
      cnt++;

  return cnt;
}


unsigned int MeshTools::n_active_local_levels(const MeshBase& mesh)
{
  unsigned int max_level = 0;

  MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();

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

  return max_level + 1;
}



unsigned int MeshTools::n_active_levels(const MeshBase& mesh)
{
  parallel_only();

  unsigned int nl = MeshTools::n_active_local_levels(mesh);

  MeshBase::const_element_iterator el =
    mesh.unpartitioned_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.unpartitioned_elements_end();

  for( ; el != end_el; ++el)
    if ((*el)->active())
      nl = std::max((*el)->level() + 1, nl);

  CommWorld.max(nl);
  return nl;
}



unsigned int MeshTools::n_local_levels(const MeshBase& mesh)
{
  unsigned int max_level = 0;

  MeshBase::const_element_iterator el = mesh.local_elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.local_elements_end();

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

  return max_level + 1;
}



unsigned int MeshTools::n_levels(const MeshBase& mesh)
{
  parallel_only();

  unsigned int nl = MeshTools::n_local_levels(mesh);

  MeshBase::const_element_iterator el =
    mesh.unpartitioned_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.unpartitioned_elements_end();

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

  CommWorld.max(nl);
  return nl;
}



void MeshTools::get_not_subactive_node_ids(const MeshBase& mesh,
    std::set<dof_id_type>& not_subactive_node_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);
    if(!elem->subactive())
      for (unsigned int n=0; n<elem->n_nodes(); ++n)
        not_subactive_node_ids.insert(elem->node(n));
  }
}



dof_id_type MeshTools::n_elem (const MeshBase::const_element_iterator &begin,
                               const MeshBase::const_element_iterator &end)
{
  return std::distance(begin, end);
}



dof_id_type MeshTools::n_nodes (const MeshBase::const_node_iterator &begin,
                                const MeshBase::const_node_iterator &end)
{
  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);

  CommWorld.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)
{
  dof_id_type 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<dof_id_type, std::vector<dof_id_type> >& 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_equal_to (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_less (s_neigh, neigh->n_neighbors());

            //Couple of helper uints...
            unsigned int local_node1=0;
            unsigned int local_node2=0;

            bool found_in_neighbor = false;

            //Find the two vertices that make up this side
            while(!elem->is_node_on_side(local_node1++,s)) { }
            local_node1--;

                //Start looking for the second one with the next node
            local_node2=local_node1+1;

            //Find the other one
            while(!elem->is_node_on_side(local_node2++,s)) { }
            local_node2--;

                //Pull out their global ids:
            dof_id_type node1 = elem->node(local_node1);
            dof_id_type node2 = elem->node(local_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;
            }

            dof_id_type hanging_node=0;

            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
            local_node1=0;
            local_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(local_node1++,s_neigh)) { }
            local_node1--;

            local_node2=local_node1+1;

            //Find the second node...
            while(!neigh->is_node_on_side(local_node2++,s_neigh)) { }
            local_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(local_node1));
              hanging_nodes[hanging_node].push_back(neigh->node(local_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;
      processor_id_type 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_equal_n_systems (const MeshBase &mesh)
{
  MeshBase::const_element_iterator el =
    mesh.elements_begin();
  const MeshBase::const_element_iterator el_end =
    mesh.elements_end();
  if (el == el_end)
    return;

  const unsigned int n_sys = (*el)->n_systems();

  for (; el != el_end; ++el)
    {
      const Elem *elem = *el;
      libmesh_assert_equal_to (elem->n_systems(), n_sys);
    }

  MeshBase::const_node_iterator node_it =
    mesh.nodes_begin();
  const MeshBase::const_node_iterator node_end =
    mesh.nodes_end();

  if (node_it == node_end)
    return;

  for (; node_it != node_end; ++node_it)
    {
      const Node *node = *node_it;
      libmesh_assert_equal_to (node->n_systems(), n_sys);
    }
}



void MeshTools::libmesh_assert_old_dof_objects (const MeshBase &mesh)
{
#ifdef LIBMESH_ENABLE_AMR
  MeshBase::const_element_iterator el =
    mesh.elements_begin();
  const MeshBase::const_element_iterator el_end =
    mesh.elements_end();

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

      if (elem->refinement_flag() == Elem::JUST_REFINED ||
          elem->refinement_flag() == Elem::INACTIVE)
        continue;

      if (elem->has_dofs())
        libmesh_assert(elem->old_dof_object);

      for (unsigned int n=0; n != elem->n_nodes(); ++n)
        {
          const Node *node = elem->get_node(n);
          if (node->has_dofs())
            libmesh_assert(elem->get_node(n)->old_dof_object);
        }
    }
#endif // LIBMESH_ENABLE_AMR
}



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_not_equal_to (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_not_equal_to (elem->neighbor(n), remote_elem);
#ifdef LIBMESH_ENABLE_AMR
      const Elem* parent = elem->parent();
      if (parent)
        {
          libmesh_assert_not_equal_to (parent, remote_elem);
          for (unsigned int c=0; c != elem->n_children(); ++c)
            libmesh_assert_not_equal_to (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_not_equal_to (elem->parent(), bad_elem);
      for (unsigned int n=0; n != elem->n_neighbors(); ++n)
        libmesh_assert_not_equal_to (elem->neighbor(n), bad_elem);
#ifdef LIBMESH_ENABLE_AMR
      if (elem->has_children())
        for (unsigned int c=0; c != elem->n_children(); ++c)
          libmesh_assert_not_equal_to (elem->child(c), bad_elem);
#endif
    }
}



void MeshTools::libmesh_assert_valid_elem_ids(const MeshBase &mesh)
{
  processor_id_type lastprocid = 0;
  dof_id_type 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);
      processor_id_type elemprocid = elem->processor_id();
      dof_id_type elemid = elem->id();

      libmesh_assert_greater_equal (elemid, lastelemid);
      libmesh_assert_greater_equal (elemprocid, lastprocid);

      lastelemid = elemid;
      lastprocid = elemprocid;
    }
}



void MeshTools::libmesh_assert_valid_amr_elem_ids(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);

      const Elem* parent = elem->parent();

      if (parent)
        {
          libmesh_assert_greater_equal (elem->id(), parent->id());
          libmesh_assert_greater_equal (elem->processor_id(), parent->processor_id());
        }
    }
}



void MeshTools::libmesh_assert_connected_nodes (const MeshBase &mesh)
{
  std::set<const Node*> used_nodes;

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

      for (unsigned int n=0; n<elem->n_nodes(); ++n)
        used_nodes.insert(elem->get_node(n));
    }

  const MeshBase::const_node_iterator node_end = mesh.nodes_end();

  for (MeshBase::const_node_iterator node_it = mesh.nodes_begin();
       node_it != node_end; ++node_it)
    {
      Node *node = *node_it;
      libmesh_assert(node);
      libmesh_assert(used_nodes.count(node));
    }
}



namespace MeshTools {

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

  parallel_only();

  dof_id_type pmax_elem_id = mesh.max_elem_id();
  CommWorld.max(pmax_elem_id);

  for (dof_id_type i=0; i != pmax_elem_id; ++i)
    assert_semiverify_dofobj(mesh.query_elem(i));

  dof_id_type pmax_node_id = mesh.max_node_id();
  CommWorld.max(pmax_node_id);

  for (dof_id_type i=0; i != pmax_node_id; ++i)
    assert_semiverify_dofobj(mesh.query_node_ptr(i));
}

template <>
void libmesh_assert_valid_procids<Elem>(const MeshBase& mesh)
{
  if (libMesh::n_processors() == 1)
    return;

  parallel_only();

  // We want this test to be valid even when called even after nodes
  // have been added asynchonously but before they're renumbered
  dof_id_type parallel_max_elem_id = mesh.max_elem_id();
  CommWorld.max(parallel_max_elem_id);

  // Check processor ids for consistency between processors

  for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
    {
      const Elem *elem = mesh.query_elem(i);

      processor_id_type min_id =
        elem ? elem->processor_id() :
               std::numeric_limits<processor_id_type>::max();
      CommWorld.min(min_id);

      processor_id_type max_id =
        elem ? elem->processor_id() :
               std::numeric_limits<processor_id_type>::min();
      CommWorld.max(max_id);

      if (elem)
        {
          libmesh_assert_equal_to (min_id, elem->processor_id());
          libmesh_assert_equal_to (max_id, elem->processor_id());
        }

      if (min_id == libMesh::processor_id())
        libmesh_assert(elem);
    }

  // If we're adaptively refining, check processor ids for consistency
  // between parents and children.
#ifdef LIBMESH_ENABLE_AMR

  // Ancestor elements we won't worry about, but subactive and active
  // elements ought to have parents with consistent processor ids

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

      if (!elem->active() && !elem->subactive())
        continue;

      const Elem *parent = elem->parent();

      if (parent)
        {
          libmesh_assert(parent->has_children());
          processor_id_type parent_procid = parent->processor_id();
          bool matching_child_id = false;
          for (unsigned int c = 0; c != parent->n_children(); ++c)
            {
              const Elem* child = parent->child(c);
              libmesh_assert(child);

              // If we've got a remote_elem then we don't know whether
              // it's responsible for the parent's processor id; all
              // we can do is assume it is and let its processor fail
              // an assert if there's something wrong.
              if (child == remote_elem ||
                  child->processor_id() == parent_procid)
                matching_child_id = true;
            }
          libmesh_assert(matching_child_id);
        }
    }
#endif
}



template <>
void libmesh_assert_valid_procids<Node>(const MeshBase& mesh)
{
  if (libMesh::n_processors() == 1)
    return;

  parallel_only();

  // We want this test to be valid even when called even after nodes
  // have been added asynchonously but before they're renumbered
  dof_id_type parallel_max_node_id = mesh.max_node_id();
  CommWorld.max(parallel_max_node_id);

  // Check processor ids for consistency between processors

  for (dof_id_type i=0; i != parallel_max_node_id; ++i)
    {
      const Node *node = mesh.query_node_ptr(i);

      processor_id_type min_id =
        node ? node->processor_id() :
               std::numeric_limits<processor_id_type>::max();
      CommWorld.min(min_id);

      processor_id_type max_id =
        node ? node->processor_id() :
               std::numeric_limits<processor_id_type>::min();
      CommWorld.max(max_id);

      if (node)
        {
          libmesh_assert_equal_to (min_id, node->processor_id());
          libmesh_assert_equal_to (max_id, node->processor_id());
        }

      if (min_id == libMesh::processor_id())
        libmesh_assert(node);
    }

  std::vector<bool> node_touched_by_me(parallel_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);
          dof_id_type nodeid = node->id();
          node_touched_by_me[nodeid] = true;
        }
    }
  std::vector<bool> node_touched_by_anyone(node_touched_by_me);
  CommWorld.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);

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

} // namespace MeshTools



#ifdef LIBMESH_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);
      dof_id_type 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);
  CommWorld.min(min_elem_h_state);

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

  for (dof_id_type 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 // LIBMESH_ENABLE_AMR



#ifdef LIBMESH_ENABLE_AMR
void MeshTools::libmesh_assert_valid_refinement_tree(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);
      if (elem->has_children())
        for (unsigned int n=0; n != elem->n_children(); ++n)
          {
            libmesh_assert(elem->child(n));
            if (elem->child(n) != remote_elem)
              libmesh_assert_equal_to (elem->child(n)->parent(), elem);
          }
      if (elem->active())
        {
          libmesh_assert(!elem->ancestor());
          libmesh_assert(!elem->subactive());
        }
      else if (elem->ancestor())
        {
          libmesh_assert(!elem->subactive());
        }
      else
        libmesh_assert(elem->subactive());
    }
}
#else
void MeshTools::libmesh_assert_valid_refinement_tree(const MeshBase &)
{
}
#endif // LIBMESH_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);
}

} // namespace libMesh

---

Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:47 UTC

Hosted By: