elem.c

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

// $Id: elem.C 2897 2008-06-30 14:25:32Z benkirk $

// 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 <algorithm> // for std::sort
#include <iterator>  // for std::ostream_iterator

// Local includes
#include "elem.h"
#include "fe_type.h"
#include "fe_interface.h"
#include "edge_edge2.h"
#include "edge_edge3.h"
#include "edge_edge4.h"
#include "edge_inf_edge2.h"
#include "face_tri3.h"
#include "face_tri6.h"
#include "face_quad4.h"
#include "face_quad8.h"
#include "face_quad9.h"
#include "face_inf_quad4.h"
#include "face_inf_quad6.h"
#include "cell_tet4.h"
#include "cell_tet10.h"
#include "cell_hex8.h"
#include "cell_hex20.h"
#include "cell_hex27.h"
#include "cell_inf_hex8.h"
#include "cell_inf_hex16.h"
#include "cell_inf_hex18.h"
#include "cell_prism6.h"
#include "cell_prism15.h"
#include "cell_prism18.h"
#include "cell_inf_prism6.h"
#include "cell_inf_prism12.h"
#include "cell_pyramid5.h"
#include "fe_base.h"
#include "quadrature_gauss.h"
#include "remote_elem.h"
#include "mesh_base.h"

// Initialize static member variables
const unsigned int Elem::_bp1 = 65449;
const unsigned int Elem::_bp2 = 48661;
const unsigned int Elem::type_to_n_nodes_map [] =
  {
    2,  // EDGE2
    3,  // EDGE3
    4,  // EDGE4
		 
    3,  // TRI3
    6,  // TRI6
		 
    4,  // QUAD4
    8,  // QUAD8
    9,  // QUAD9
    
    4,  // TET4
    10, // TET10
		 
    8,  // HEX8
    20, // HEX20
    27, // HEX27
		 
    6,  // PRISM6
    15, // PRISM15
    18, // PRISM18
    
    5,  // PYRAMID5
    
    2,  // INFEDGE2
    
    4,  // INFQUAD4
    6,  // INFQUAD6
		 
    8,  // INFHEX8
    16, // INFHEX16
    18, // INFHEX18
		 
    6,  // INFPRISM6
    16, // INFPRISM12

    1,  // NODEELEM
  };



// ------------------------------------------------------------
// Elem class member funcions
AutoPtr<Elem> Elem::build(const ElemType type,
			  Elem* p)
{
  Elem* elem = NULL;
 
  switch (type)
    {
      // 1D elements 
    case EDGE2:
      {
	elem = new Edge2(p);
	break;
      }
    case EDGE3:
      {
	elem = new Edge3(p);
	break;
      }
    case EDGE4:
      {
        elem = new Edge4(p);
        break;
      }


      
      // 2D elements
    case TRI3:
      {
	elem = new Tri3(p);
	break;
      }
    case TRI6:
      {
	elem = new Tri6(p);
	break;
      }
    case QUAD4:
      {
	elem = new Quad4(p);
	break;
      }
    case QUAD8:
      {
	elem = new Quad8(p);
	break;
      }
    case QUAD9:
      {
	elem = new Quad9(p);
	break;
      }
   

      // 3D elements
    case TET4:
      {
	elem = new Tet4(p);
	break;
      }
    case TET10:
      {
	elem = new Tet10(p);
	break;
      }
    case HEX8:
      {
	elem = new Hex8(p);
	break;
      }
    case HEX20:
      {
	elem = new Hex20(p);
	break;
      }
    case HEX27:
      {
	elem = new Hex27(p);
	break;
      }
    case PRISM6:
      {
	elem = new Prism6(p);
	break;
      }
    case PRISM15:
      {
	elem = new Prism15(p);
	break;
      }
    case PRISM18:
      {
	elem = new Prism18(p);
	break;
      }
    case PYRAMID5:
      {
	elem = new Pyramid5(p);
	break;
      }



#ifdef ENABLE_INFINITE_ELEMENTS

      // 1D infinite elements
    case INFEDGE2:
      {
	elem = new InfEdge2(p);
	break;
      }


      // 2D infinite elements
    case INFQUAD4:
      {
	elem = new InfQuad4(p);
	break;
      }
    case INFQUAD6:
      {
	elem = new InfQuad6(p);
	break;
      }

   
    // 3D infinite elements
    case INFHEX8:
      {
	elem = new InfHex8(p);
	break;
      }
    case INFHEX16:
      {
	elem = new InfHex16(p);
	break;
      }
    case INFHEX18:
      {
	elem = new InfHex18(p);
	break;
      }
    case INFPRISM6:
      {
	elem = new InfPrism6(p);
	break;
      }
    case INFPRISM12:
      {
	elem = new InfPrism12(p);
	break;
      }

#endif
           
    default:
      {
	std::cerr << "ERROR: Undefined element type!." << std::endl;
	libmesh_error();
      }
    }
  

  AutoPtr<Elem> ap(elem);
  return ap;
}



Point Elem::centroid() const
{
  Point cp;

  for (unsigned int n=0; n<this->n_vertices(); n++)
    cp.add (this->point(n));

  return (cp /= static_cast<Real>(this->n_vertices()));    
}



Real Elem::hmin() const
{
  Real h_min=1.e30;

  for (unsigned int n_outer=0; n_outer<this->n_vertices(); n_outer++)
    for (unsigned int n_inner=n_outer+1; n_inner<this->n_vertices(); n_inner++)
      {
	const Point diff = (this->point(n_outer) - this->point(n_inner));
	
	h_min = std::min(h_min,diff.size());
      }

  return h_min;
}



Real Elem::hmax() const
{
  Real h_max=0;

  for (unsigned int n_outer=0; n_outer<this->n_vertices(); n_outer++)
    for (unsigned int n_inner=n_outer+1; n_inner<this->n_vertices(); n_inner++)
      {
	const Point diff = (this->point(n_outer) - this->point(n_inner));
	
	h_max = std::max(h_max,diff.size());
      }

  return h_max;
}



Real Elem::length(const unsigned int n1, 
		  const unsigned int n2) const
{
  libmesh_assert ( n1 < this->n_vertices() );
  libmesh_assert ( n2 < this->n_vertices() );

  return (this->point(n1) - this->point(n2)).size();
}



bool Elem::operator == (const DofObject& rhs) const
{

    // Cast rhs to an Elem*
    const Elem* rhs_elem = dynamic_cast<const Elem*>(&rhs);

    // If we cannot cast to an Elem*, rhs must be a Node
    if(rhs_elem == static_cast<const Elem*>(NULL))
        return false;

//   libmesh_assert (n_nodes());
//   libmesh_assert (rhs.n_nodes());

//   // Elements can only be equal if they
//   // contain the same number of nodes.
//   if (this->n_nodes() == rhs.n_nodes())
//     {
//       // Create a set that contains our global
//       // node numbers and those of our neighbor.
//       // If the set is the same size as the number
//       // of nodes in both elements then they must
//       // be connected to the same nodes.     
//       std::set<unsigned int> nodes_set;

//       for (unsigned int n=0; n<this->n_nodes(); n++)
//         {
//           nodes_set.insert(this->node(n));
//           nodes_set.insert(rhs.node(n));
//         }

//       // If this passes the elements are connected
//       // to the same global nodes
//       if (nodes_set.size() == this->n_nodes())
//         return true;
//     }

//   // If we get here it is because the elements either
//   // do not have the same number of nodes or they are
//   // connected to different nodes.  Either way they
//   // are not the same element
//   return false;
  
  // Useful typedefs
  typedef std::vector<unsigned int>::iterator iterator;

  
  // Elements can only be equal if they
  // contain the same number of nodes.
  // However, we will only test the vertices,
  // which is sufficient & cheaper
  if (this->n_nodes() == rhs_elem->n_nodes())
    {
      // The number of nodes in the element
      const unsigned int nn = this->n_nodes();
      
      // Create a vector that contains our global
      // node numbers and those of our neighbor.
      // If the sorted, unique vector is the same size
      // as the number of nodes in both elements then
      // they must be connected to the same nodes.
      //
      // The vector will be no larger than 2*n_nodes(),
      // so we might as well reserve the space.
      std::vector<unsigned int> common_nodes;
      common_nodes.reserve (2*nn);

      // Add the global indices of the nodes
      for (unsigned int n=0; n<nn; n++)
	{
	  common_nodes.push_back (this->node(n));
	  common_nodes.push_back (rhs_elem->node(n));
	}

      // Sort the vector and find out how long
      // the sorted vector is.
      std::sort (common_nodes.begin(), common_nodes.end());
      
      iterator new_end = std::unique (common_nodes.begin(),
				      common_nodes.end());
      
      const int new_size = std::distance (common_nodes.begin(),
					  new_end);
      
      // If this passes the elements are connected
      // to the same global vertex nodes
      if (new_size == static_cast<int>(nn))
	return true;
    }

  // If we get here it is because the elements either
  // do not have the same number of nodes or they are
  // connected to different nodes.  Either way they
  // are not the same element
  return false;
}



bool Elem::contains_vertex_of(const Elem *e) const
{
  // Our vertices are the first numbered nodes
  for (unsigned int n = 0; n != e->n_vertices(); ++n)
    if (this->contains_point(e->point(n)))
      return true;
  return false;
}



void Elem::find_point_neighbors(std::set<const Elem *> &neighbor_set) const
{
  neighbor_set.clear();
  neighbor_set.insert(this);

  unsigned int old_size;
  do
    {
      old_size = neighbor_set.size();

      // Loop over all the elements in the patch
      std::set<const Elem*>::const_iterator       it  = neighbor_set.begin();
      const std::set<const Elem*>::const_iterator end = neighbor_set.end();

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

          for (unsigned int s=0; s<elem->n_sides(); s++)
            {
              const Elem* neighbor = elem->neighbor(s);
              if (neighbor &&
                  neighbor != remote_elem)    // we have a real neighbor on this side
                {
                  if (neighbor->active())                // ... if it is active
                    {
                      if (this->contains_vertex_of(neighbor) // ... and touches us
                          || neighbor->contains_vertex_of(this))  
                        neighbor_set.insert (neighbor);  // ... then add it
                    }
#ifdef ENABLE_AMR
                  else                                 // ... the neighbor is *not* active,
                    {                                  // ... so add *all* neighboring
                                                       // active children
                      std::vector<const Elem*> active_neighbor_children;
  
                      neighbor->active_family_tree_by_neighbor
                        (active_neighbor_children, elem);

                      std::vector<const Elem*>::const_iterator
                        child_it = active_neighbor_children.begin();
                      const std::vector<const Elem*>::const_iterator
                        child_end = active_neighbor_children.end();
                      for (; child_it != child_end; ++child_it)