mesh_tools.c

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

// $Id: mesh_tools.C 2970 2008-08-11 00:17:08Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// 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 <set>

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

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



// ------------------------------------------------------------
// anonymous namespace for helper classes
namespace {
  /**
   * SumElemWeight(Range) sums the number of nodes per element 
   * for each element in the provided range. The join() method
   * defines how to combine the reduction operation from two
   * distinct instances of this class which may be executed on 
   * separate threads.
   */
  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();
    }

    unsigned int weight() const 
    { return _weight; }
    
    void join (const SumElemWeight &other)
    { _weight += other.weight(); }

  private:
    unsigned int _weight;
  };


  /**
   * FindBBox(Range) computes the bounding box for the objects
   * in the specified range.  This class may be split and subranges
   * can be executed on separate threads.  The join() method
   * defines how the results from two separate threads are combined.
   */
  class FindBBox
  {
  public:
    FindBBox () :
      _vmin(3,  std::numeric_limits<Real>::max()),
      _vmax(3, -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 != NULL);
	  
	  for (unsigned int i=0; i<3; 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 != NULL);

	  for (unsigned int n=0; n<elem->n_nodes(); n++)
	    for (unsigned int i=0; i<3; i++)
	      {
		_vmin[i] = std::min(_vmin[i], elem->point(n)(i));
		_vmax[i] = std::max(_vmax[i], elem->point(n)(i));
	      }      
        }
    }

    void join (const FindBBox &other)
    {
      for (unsigned int i=0; i<3; i++)
	{
	  _vmin[i] = std::min(_vmin[i], other._vmin[i]);
	  _vmax[i] = std::max(_vmax[i], other._vmax[i]);
	}      
    }

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

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

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



// ------------------------------------------------------------
// MeshTools functions
unsigned int MeshTools::total_weight(const MeshBase& mesh)
{
  if (!mesh.is_serial())
    {
      parallel_only();
      unsigned int weight = MeshTools::weight (mesh, libMesh::processor_id());
      Parallel::sum(weight);
      unsigned int 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();

}



unsigned int MeshTools::weight(const MeshBase& mesh, const unsigned int 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<unsigned int> >& 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 ((*el)->node(n) < nodes_to_elem_map.size());
	libmesh_assert ((*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 ((*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
  Parallel::min(find_bbox.min());
  Parallel::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 unsigned int pid)
{
  libmesh_assert (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 unsigned int 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 unsigned int sid)
{
  libmesh_assert (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));
	  }      
  
  const BoundingBox ret_val(min, max);
  
  return ret_val;  
}



Sphere
MeshTools::subdomain_bounding_sphere (const MeshBase& mesh,
				      const unsigned int 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());
}



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



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

unsigned int MeshTools::n_non_subactive_elem_of_type_at_level(const MeshBase& mesh,
                                                const ElemType type,
                                                const unsigned int level)
{
  unsigned int 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);

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

  Parallel::max(nl);
  return nl;
}



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



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



unsigned int MeshTools::n_nodes (const MeshBase::const_node_iterator &begin,
				 const MeshBase::const_node_iterator &end)
{