Poster of Linux kernelThe best gift for a Linux geek
InfElemBuilder

InfElemBuilder

Section: C Library Functions (3) Updated: Thu Apr 7 2011
Local index Up
 

NAME

InfElemBuilder -  

SYNOPSIS


#include <inf_elem_builder.h>  

Public Types


typedef std::pair< bool, double > InfElemOriginValue
 

Public Member Functions


InfElemBuilder (MeshBase &mesh)

const Point build_inf_elem (const bool be_verbose=false)

const Point build_inf_elem (const InfElemOriginValue &origin_x, const InfElemOriginValue &origin_y, const InfElemOriginValue &origin_z, const bool x_sym=false, const bool y_sym=false, const bool z_sym=false, const bool be_verbose=false, std::vector< const Node * > *inner_boundary_nodes=NULL)
 

Private Member Functions


void build_inf_elem (const Point &origin, const bool x_sym=false, const bool y_sym=false, const bool z_sym=false, const bool be_verbose=false, std::set< std::pair< unsigned int, unsigned int > > *inner_faces=NULL)
 

Private Attributes


MeshBase & _mesh
 

Detailed Description

This class is used to build infinite elements on top of an existing mesh. It only makes sense to use this if LIBMESH_ENABLE_INFINITE_ELEMENTS is true.

Author:

original code by Daniel Dreyer, converted to InfElemBuilder class by John W. Peterson, 2004.

Definition at line 49 of file inf_elem_builder.h.  

Member Typedef Documentation

 

typedef std::pair<bool, double> InfElemBuilder::InfElemOriginValueUseful typedef

Definition at line 60 of file inf_elem_builder.h.  

Constructor & Destructor Documentation

 

InfElemBuilder::InfElemBuilder (MeshBase &mesh) [inline]Constructor.

Definition at line 55 of file inf_elem_builder.h.

: _mesh(mesh) {}
 

Member Function Documentation

 

const Point InfElemBuilder::build_inf_elem (const boolbe_verbose = false)Build infinite elements atop a volume-based mesh, determine origin automatically. Also returns the origin as a const Point to make it more obvious that the origin should not change after the infinite elements have been built. When symmetry planes are present, use the version with optional symmetry switches. The flag be_verbose enables some diagnostic output.

Definition at line 43 of file inf_elem_builder.C.

References _mesh, MeshTools::bounding_box(), MeshBase::prepare_for_use(), libMesh::processor_id(), and TypeVector< T >::write_unformatted().

Referenced by build_inf_elem().

{
  // 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;
}
 

void InfElemBuilder::build_inf_elem (const Point &origin, const boolx_sym = false, const booly_sym = false, const boolz_sym = false, const boolbe_verbose = false, std::set< std::pair< unsigned int, unsigned int > > *inner_faces = NULL) [private]Build infinite elements atop a volume-based mesh. Actual implementation.

Definition at line 278 of file inf_elem_builder.C.

References _mesh, MeshBase::active_elements_begin(), MeshBase::active_elements_end(), MeshBase::add_elem(), MeshBase::add_point(), Elem::build_side(), libMeshEnums::EDGE2, libMeshEnums::EDGE3, MeshBase::elem(), MeshBase::find_neighbors(), DofObject::id(), MeshBase::is_serial(), ParallelMesh::libmesh_assert_valid_parallel_ids(), MeshBase::max_elem_id(), MeshBase::max_node_id(), MeshBase::n_elem(), Elem::n_neighbors(), Elem::n_vertices(), Elem::neighbor(), MeshBase::node(), MeshBase::point(), DofObject::processor_id(), libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, DofObject::set_id(), Elem::set_node(), TypeVector< T >::size(), libMeshEnums::TRI3, and libMeshEnums::TRI6.

{
  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.
  {
    MeshBase::element_iterator       it  = this->_mesh.active_elements_begin();
    const MeshBase::element_iterator end = this->_mesh.active_elements_end();
  
    for(; it != end; ++it)
      {
        Elem* elem = *it;
        
        for (unsigned int s=0; s<elem->n_neighbors(); s++)
          {
            // check if elem(e) is on the boundary
            if (elem->neighbor(s) == NULL)
              {
                // note that it is safe to use the Elem::side() method, 
                // which gives a non-full-ordered element 
                AutoPtr<Elem> side(elem->build_side(s));

                // bool flags for symmetry detection
                bool sym_side=false;            
                bool on_x_sym=true;
                bool on_y_sym=true;
                bool on_z_sym=true;
                
                
                // Loop over the nodes to check whether they are on the symmetry planes,
                // and therefore sufficient to use a non-full-ordered side element
                for(unsigned int n=0; n<side->n_nodes(); n++)
                  {
                    const Point dist_from_origin = this->_mesh.point(side->node(n)) - origin;

                    if(x_sym)
                      if( std::abs(dist_from_origin(0)) > 1.e-3 )
                        on_x_sym=false;

                    if(y_sym)
                      if( std::abs(dist_from_origin(1)) > 1.e-3 )
                        on_y_sym=false;

                    if(z_sym)
                      if( std::abs(dist_from_origin(2)) > 1.e-3 )
                        on_z_sym=false;
              
                    //        if(x_sym)
                    //          if( std::abs(dist_from_origin(0)) > 1.e-6 )
                    //            on_x_sym=false;

                    //        if(y_sym)
                    //          if( std::abs(dist_from_origin(1)) > 1.e-6 )
                    //            on_y_sym=false;

                    //        if(z_sym)
                    //          if( std::abs(dist_from_origin(2)) > 1.e-6 )
                    //            on_z_sym=false;

                    //find the node most distant from origin
              
                    Real r = dist_from_origin.size();
                    if (r > max_r)
                      {
                        max_r = r;
                        max_r_node=side->node(n);
                      }

                  }

                sym_side = (x_sym && on_x_sym) || (y_sym && on_y_sym) || (z_sym && on_z_sym);
                                
                if (!sym_side)
                  faces.insert( std::make_pair(elem->id(), s) );                                        
                                            
              } // neighbor(s) == NULL
          } // sides
      } // elems
  }



  
                        

  //  If a boundary side has one node on the outer boundary,
  //  all points of this side are on the outer boundary.
  //  Start with the node most distant from origin, which has
  //  to be on the outer boundary, then recursively find all
  //  sides and nodes connected to it. Found sides are moved
  //  from faces to ofaces, nodes are collected in onodes.  
  //  Here, the search is done iteratively, because, depending on
  //  the mesh, a very high level of recursion might be necessary.
  onodes.insert(max_r_node);    


  {
    std::set< std::pair<unsigned int,unsigned int> >::iterator face_it = faces.begin();
    unsigned int facesfound=0;
    do {
                        
      std::pair<unsigned int, unsigned int> p;
      p = *face_it;

      // This has to be a full-ordered side element,
      // since we need the correct n_nodes,
      AutoPtr<Elem> side(this->_mesh.elem(p.first)->build_side(p.second));
                
      bool found=false;
      for(unsigned int sn=0; sn<side->n_nodes(); sn++)
        if(onodes.count(side->node(sn)))
          {
            found=true;
            break;
          }
        
        
      // If a new oface is found, include its nodes in onodes
      if(found)         
        {
          for(unsigned int sn=0; sn<side->n_nodes(); sn++)
            onodes.insert(side->node(sn));
                                
          ofaces.insert(p);
          face_it++;                    // iteration is done here
          faces.erase(p);
                
          facesfound++;
        }
    
      else
        face_it++;                      // iteration is done here


      // If at least one new oface was found in this cycle, 
      // do another search cycle.
      if(facesfound>0 && face_it == faces.end())
        {       
          facesfound = 0;
          face_it    = faces.begin();
        }
    
    } while(face_it != faces.end());
  }
  
        
#ifdef DEBUG
  if (be_verbose)
    std::cout << '  found ' 
              << faces.size() 
              << ' inner and ' 
              << ofaces.size() 
              << ' outer boundary faces' 
              << std::endl;
#endif  
        
  // When the user provided a non-null pointer to
  // inner_faces, that implies he wants to have
  // this std::set.  For now, simply copy the data.
  if (inner_faces != NULL)
    *inner_faces = faces;

  // free memory, clear our local variable, no need
  // for it any more.
  faces.clear();


  // outer_nodes maps onodes to their duplicates
  std::map<unsigned int, Node *> outer_nodes;

  // We may need to pick our own object ids in parallel
  unsigned int old_max_node_id = _mesh.max_node_id();
  unsigned int old_max_elem_id = _mesh.max_elem_id();

  // for each boundary node, add an outer_node with 
  // double distance from origin.
  std::set<unsigned int>::iterator on_it = onodes.begin();
  for( ; on_it != onodes.end(); ++on_it)
    {
      Point p = (Point(this->_mesh.point(*on_it)) * 2.) - origin;
      if (_mesh.is_serial())
        {
          // Add with a default id in serial
          outer_nodes[*on_it]=this->_mesh.add_point(p);
        }
      else
        {
          // Pick a unique id in parallel
          Node &bnode = _mesh.node(*on_it);
          unsigned int new_id = bnode.id() + old_max_node_id;
          outer_nodes[*on_it] = 
            this->_mesh.add_point(p, new_id,
                                  bnode.processor_id());
        }
    }


#ifdef DEBUG
  // for verbose, remember n_elem
  unsigned int n_conventional_elem = this->_mesh.n_elem();
#endif


  // build Elems based on boundary side type
  std::set< std::pair<unsigned int,unsigned int> >::iterator face_it = ofaces.begin();
  for( ; face_it != ofaces.end(); ++face_it)
    {
      // Shortcut to the pair being iterated over
      std::pair<unsigned int,unsigned int> p = *face_it;
        
      // build a full-ordered side element to get the base nodes
      AutoPtr<Elem> side(this->_mesh.elem(p.first)->build_side(p.second)); 

      // create cell depending on side type, assign nodes,
      // use braces to force scope.
      bool is_higher_order_elem = false;
      {
        Elem* el;
        switch(side->type())
        {
          // 3D infinite elements
          // TRIs                                       
          case TRI3:    
            el=new InfPrism6;
            break;
                                                        
          case TRI6: 
            el=new InfPrism12; 
            is_higher_order_elem = true;
            break;
                                                        
          // QUADs                                      
          case QUAD4: 
            el=new InfHex8;
            break;
                                                        
          case QUAD8: 
            el=new InfHex16;
            is_higher_order_elem = true;
            break;
                                                        
          case QUAD9: 
            el=new InfHex18;

            // the method of assigning nodes (which follows below)
            // omits in the case of QUAD9 the bubble node; therefore
            // we assign these first by hand here.
            el->set_node(16) = side->get_node(8);
            el->set_node(17) = outer_nodes[side->node(8)];
            is_higher_order_elem=true;
            break;

          // 2D infinite elements                
          case EDGE2:
            el=new InfQuad4;
            break;

          case EDGE3:
            el=new InfQuad6;
            el->set_node(4) = side->get_node(2);
            break;

          // 1D infinite elements not supported
          default: 
            std::cout << 'InfElemBuilder::build_inf_elem(Point, bool, bool, bool, bool): '
                      << 'invalid face element '
                      << std::endl;
            continue;
        }

        // In parallel, assign unique ids to the new element
        if (!_mesh.is_serial())
          {
            Elem *belem = _mesh.elem(p.first);
            el->processor_id() = belem->processor_id();
            // We'd better not have elements with more than 6 sides
            el->set_id (belem->id() * 6 + p.second + old_max_elem_id);
          }

        // assign vertices to the new infinite element
        const unsigned int n_base_vertices = side->n_vertices();
        for(unsigned int i=0; i<n_base_vertices; i++)
          {
            el->set_node(i                ) = side->get_node(i);
            el->set_node(i+n_base_vertices) = outer_nodes[side->node(i)];
          }


        // when this is a higher order element, 
        // assign also the nodes in between
        if (is_higher_order_elem)
          {
            // n_safe_base_nodes is the number of nodes in 
side // that may be safely assigned using below for loop. // Actually, n_safe_base_nodes is _identical_ with el->n_vertices(), // since for QUAD9, the 9th node was already assigned above const unsigned int n_safe_base_nodes = el->n_vertices(); for(unsigned int i=n_base_vertices; i<n_safe_base_nodes; i++) { el->set_node(i+n_base_vertices) = side->get_node(i); el->set_node(i+n_safe_base_nodes) = outer_nodes[side->node(i)]; } } // add infinite element to mesh this->_mesh.add_elem(el); } // el goes out of scope } // for #ifdef DEBUG ParallelMesh *pmesh = dynamic_cast<ParallelMesh *>(&this->_mesh); if (pmesh) pmesh->libmesh_assert_valid_parallel_ids(); if (be_verbose) std::cout << ' added ' << this->_mesh.n_elem() - n_conventional_elem << ' infinite elements and ' << onodes.size() << ' nodes to the mesh' << std::endl << std::endl; #endif STOP_LOG('build_inf_elem()', 'InfElemBuilder'); }
 

const Point InfElemBuilder::build_inf_elem (const InfElemOriginValue &origin_x, const InfElemOriginValue &origin_y, const InfElemOriginValue &origin_z, const boolx_sym = false, const booly_sym = false, const boolz_sym = false, const boolbe_verbose = false, std::vector< const Node * > *inner_boundary_nodes = NULL)Returns:

the origin of the infinite elements. Builds infinite elements atop a volume-based mesh. Finds all faces on the outer boundary and build infinite elements on them. Using the InfElemOriginValue the user can prescribe only selected origin coordinates. The remaining coordinates are computed from the center of the bounding box of the mesh.

During the search for faces on which infinite elements are built, interior faces that are not on symmetry planes are found, too. When an (optional) pointer to inner_boundary_nodes is provided, then this vector will be filled with the nodes that lie on the inner boundary.

Faces which lie in at least one symmetry plane are skipped. The three optional booleans x_sym, y_sym, z_sym indicate symmetry planes (through the origin, obviously) perpendicular to the x, y and z direction, respectively. The flag be_verbose enables some diagnostic output.

Definition at line 84 of file inf_elem_builder.C.

References _mesh, MeshTools::bounding_box(), build_inf_elem(), Elem::build_side(), MeshBase::elem(), MeshBase::node(), MeshBase::prepare_for_use(), MeshBase::print_info(), and TypeVector< T >::write_unformatted().

{
  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 
// 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
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; }
 

Member Data Documentation

 

MeshBase& InfElemBuilder::_mesh [private]Reference to the mesh we're building infinite elements for.

Definition at line 122 of file inf_elem_builder.h.

Referenced by build_inf_elem().

 

Author

Generated automatically by Doxygen for libMesh from the source code.


 

Index

NAME
SYNOPSIS
Public Types
Public Member Functions
Private Member Functions
Private Attributes
Detailed Description
Member Typedef Documentation
typedef std::pair<bool, double> InfElemBuilder::InfElemOriginValueUseful typedef
Constructor & Destructor Documentation
InfElemBuilder::InfElemBuilder (MeshBase &mesh) [inline]Constructor.
Member Function Documentation
const Point InfElemBuilder::build_inf_elem (const boolbe_verbose = false)Build infinite elements atop a volume-based mesh, determine origin automatically. Also returns the origin as a const Point to make it more obvious that the origin should not change after the infinite elements have been built. When symmetry planes are present, use the version with optional symmetry switches. The flag be_verbose enables some diagnostic output.
void InfElemBuilder::build_inf_elem (const Point &origin, const boolx_sym = false, const booly_sym = false, const boolz_sym = false, const boolbe_verbose = false, std::set< std::pair< unsigned int, unsigned int > > *inner_faces = NULL) [private]Build infinite elements atop a volume-based mesh. Actual implementation.
const Point InfElemBuilder::build_inf_elem (const InfElemOriginValue &origin_x, const InfElemOriginValue &origin_y, const InfElemOriginValue &origin_z, const boolx_sym = false, const booly_sym = false, const boolz_sym = false, const boolbe_verbose = false, std::vector< const Node * > *inner_boundary_nodes = NULL)Returns:
Member Data Documentation
MeshBase& InfElemBuilder::_mesh [private]Reference to the mesh we're building infinite elements for.
Author

This document was created by man2html, using the manual pages.
Time: 21:48:29 GMT, April 16, 2011