inf_elem_builder.c

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

// $Id: inf_elem_builder.C 2808 2008-04-25 18:41:43Z 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

#include "libmesh_config.h"

#ifdef ENABLE_INFINITE_ELEMENTS

// C++ includes

// Local includes
#include "inf_elem_builder.h"
#include "libmesh_logging.h"
#include "mesh_tools.h"
#include "face_inf_quad4.h"
#include "face_inf_quad6.h"
#include "cell_inf_prism6.h"
#include "cell_inf_prism12.h"
#include "cell_inf_hex8.h"
#include "cell_inf_hex16.h"
#include "cell_inf_hex18.h"
#include "mesh_base.h"

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

const Point InfElemBuilder::build_inf_elem(bool be_verbose)
{
  // determine origin automatically,
  // works only if the mesh has no symmetry planes.
  const MeshTools::BoundingBox b_box = MeshTools::bounding_box(_mesh);
  Point origin = (b_box.first + b_box.second) / 2.;

  if (be_verbose && libMesh::processor_id() == 0)
    {
#ifdef DEBUG
      std::cout << " Determined origin for Infinite Elements:" 
		<< std::endl
		<< "  ";
      origin.write_unformatted(std::cout);
      std::cout << std::endl;
#endif
    }

  // Call the protected implementation function with the
  // automatically determined origin.
  this->build_inf_elem(origin, false, false, false, be_verbose);

  // when finished with building the Ifems,
  // it remains to prepare the mesh for use:
  // find neighbors (again), partition (if needed)...
  this->_mesh.prepare_for_use ();

  return origin;
}












const Point InfElemBuilder::build_inf_elem (const InfElemOriginValue& origin_x,
					    const InfElemOriginValue& origin_y,
					    const InfElemOriginValue& origin_z,
					    const bool x_sym,
					    const bool y_sym,
					    const bool z_sym,
					    const bool be_verbose,
					    std::vector<const Node*>* inner_boundary_nodes)
{
  START_LOG("build_inf_elem()", "InfElemBuilder");

  // first determine the origin of the 
  // infinite elements.  For this, the
  // origin defaults to the given values,
  // and may be overridden when the user
  // provided values
  Point origin(origin_x.second, origin_y.second, origin_z.second);

  // when only _one_ of the origin coordinates is _not_
  // given, we have to determine it on our own
  if ( !origin_x.first || !origin_y.first || !origin_z.first)
    {
      // determine origin
      const MeshTools::BoundingBox b_box = MeshTools::bounding_box(_mesh);
      const Point auto_origin = (b_box.first+b_box.second)/2.;

      // override default values, if necessary
      if (!origin_x.first)
	  origin(0) = auto_origin(0);
      if (!origin_y.first)
	  origin(1) = auto_origin(1);
      if (!origin_z.first)
	  origin(2) = auto_origin(2);

      if (be_verbose)
        {
	  std::cout << " Origin for Infinite Elements:" << std::endl;

	  if (!origin_x.first)
	      std::cout << "  determined x-coordinate" << std::endl;
	  if (!origin_y.first)
	      std::cout << "  determined y-coordinate" << std::endl;
	  if (!origin_z.first)
	      std::cout << "  determined z-coordinate" << std::endl;

	  std::cout << "  coordinates: ";
	  origin.write_unformatted(std::cout);
	  std::cout << std::endl;
	}
    }

  else if (be_verbose)

    {
      std::cout << " Origin for Infinite Elements:" << std::endl;
      std::cout << "  coordinates: ";
      origin.write_unformatted(std::cout);
      std::cout << std::endl;
    }



  // Now that we have the origin, check if the user provided an \p
  // inner_boundary_nodes.  If so, we pass a std::set to the actual
  // implementation of the build_inf_elem(), so that we can convert
  // this to the Node* vector
  if (inner_boundary_nodes != NULL)
    {
      // note that the std::set that we will get
      // from build_inf_elem() uses the index of
      // the element in this->_elements vector,
      // and the second entry is the side index
      // for this element.  Therefore, we do _not_
      // need to renumber nodes and elements
      // prior to building the infinite elements.
      //
      // However, note that this method here uses
      // node id's... Do we need to renumber?


      // Form the list of faces of elements which finally
      // will tell us which nodes should receive boundary
      // conditions (to form the std::vector<const Node*>)
      std::set< std::pair<unsigned int,
	                  unsigned int> > inner_faces;


      // build infinite elements
      this->build_inf_elem(origin, 
			   x_sym, y_sym, z_sym, 
			   be_verbose,
			   &inner_faces);

      if (be_verbose)
        {
	  this->_mesh.print_info();
	  std::cout << "Data pre-processing:" << std::endl
		    << " convert the <int,int> list to a Node* list..."
		    << std::endl;
	}

      // First use a std::vector<unsigned int> that holds
      // the global node numbers.  Then sort this vector,
      // so that it can be made unique (no multiple occurence
      // of a node), and then finally insert the Node* in
      // the vector inner_boundary_nodes.
      //
      // Reserve memory for the vector<unsigned int> with 
      // 4 times the size of the number of elements in the 
      // std::set. This is a good bet for Quad4 face elements.  
      // For higher-order elements, this probably _has_ to lead
      // to additional allocations...
      // Practice has to show how this affects performance.
      std::vector<unsigned int> inner_boundary_node_numbers;
      inner_boundary_node_numbers.reserve(4*inner_faces.size());

      // Now transform the set of pairs to a list of (possibly
      // duplicate) global node numbers.
      std::set< std::pair<unsigned int,unsigned int> >::iterator face_it = inner_faces.begin();
      const std::set< std::pair<unsigned int,unsigned int> >::iterator face_end = inner_faces.end();
      for(; face_it!=face_end; ++face_it)
        {
	  std::pair<unsigned int,unsigned int> p = *face_it;
	
	  // build a full-ordered side element to get _all_ the base nodes
	  AutoPtr<Elem> side( this->_mesh.elem(p.first)->build_side(p.second) ); 

	  // insert all the node numbers in inner_boundary_node_numbers
	  for (unsigned int n=0; n< side->n_nodes(); n++)
	      inner_boundary_node_numbers.push_back(side->node(n));
	}


      // inner_boundary_node_numbers now still holds multiple entries of
      // node numbers.  So first sort, then unique the vector.
      // Note that \p std::unique only puts the new ones in
      // front, while to leftovers are @e not deleted.  Instead,
      // it returns a pointer to the end of the unique range.
      //TODO:[BSK] int_ibn_size_before is not the same type as unique_size!
      const unsigned int ibn_size_before = inner_boundary_node_numbers.size();
      std::sort (inner_boundary_node_numbers.begin(), inner_boundary_node_numbers.end());
      std::vector<unsigned int>::iterator unique_end = 
	  std::unique (inner_boundary_node_numbers.begin(), inner_boundary_node_numbers.end());

      const int unique_size = std::distance(inner_boundary_node_numbers.begin(), unique_end);
      libmesh_assert (unique_size <= static_cast<int>(ibn_size_before));

      // Finally, create const Node* in the inner_boundary_nodes
      // vector.  Reserve, not resize (otherwise, the push_back
      // would append the interesting nodes, while NULL-nodes
      // live in the resize'd area...
      inner_boundary_nodes->reserve (unique_size);
      inner_boundary_nodes->clear();


      std::vector<unsigned int>::iterator pos_it = inner_boundary_node_numbers.begin();
      for (; pos_it != unique_end; ++pos_it)
        {
	  const Node& node = this->_mesh.node(*pos_it);
	  inner_boundary_nodes->push_back(&node);
	}

      if (be_verbose)
	  std::cout << "  finished identifying " << unique_size 
		    << " target nodes." << std::endl;  
    }

  else

    {
      // There are no inner boundary nodes, so simply build the infinite elements
      this->build_inf_elem(origin, x_sym, y_sym, z_sym, be_verbose);
    }


  STOP_LOG("build_inf_elem()", "InfElemBuilder");

  // when finished with building the Ifems,
  // it remains to prepare the mesh for use:
  // find neighbors again, partition (if needed)...
  this->_mesh.prepare_for_use ();

  return origin;
}









// The actual implementation of building elements.
void InfElemBuilder::build_inf_elem(const Point& origin, 
				    const bool x_sym, 
				    const bool y_sym, 
				    const bool z_sym,
				    const bool be_verbose,
				    std::set< std::pair<unsigned int,
				    unsigned int> >* inner_faces)
{
  if (be_verbose)
    {
#ifdef DEBUG
      std::cout << " Building Infinite Elements:" << std::endl;
      std::cout << "  updating element neighbor tables..." << std::endl;
#else
      std::cout << " Verbose mode disabled in non-debug mode." << std::endl;
#endif
    }


  // update element neighbors
  this->_mesh.find_neighbors();	

  START_LOG("build_inf_elem()", "InfElemBuilder");

  // A set for storing element number, side number pairs.
  // pair.first == element number, pair.second == side number
  std::set< std::pair<unsigned int,unsigned int> > faces;
  std::set< std::pair<unsigned int,unsigned int> > ofaces;

  // A set for storing node numbers on the outer faces.
  std::set<unsigned int> onodes;

  // The distance to the farthest point in the mesh from the origin
  Real max_r=0.;

  // The index of the farthest point in the mesh from the origin
  unsigned int max_r_node;

#ifdef DEBUG
  if (be_verbose)
    {
      std::cout << "  collecting boundary sides";
      if (x_sym || y_sym || z_sym)
	std::cout << ", skipping sides in symmetry planes..." << std::endl;
      else
	std::cout << "..." << std::endl;
    }
#endif

  // Iterate through all elements and sides, collect indices of all active
  // boundary sides in the faces set. Skip sides which lie in symmetry planes.
  // Later, sides of the inner boundary will be sorted out.
  {