#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< dof_id_type, 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 52 of file inf_elem_builder.h.

Member Typedef Documentation

typedef std::pair<bool, double> libMesh::InfElemBuilder::InfElemOriginValue

Useful typedef

Definition at line 64 of file inf_elem_builder.h.

Constructor & Destructor Documentation

libMesh::InfElemBuilder::InfElemBuilder ( MeshBase & mesh ) [inline, explicit]

Constructor.

Definition at line 59 of file inf_elem_builder.h.

: _mesh(mesh) {}

Member Function Documentation

const Point libMesh::InfElemBuilder::build_inf_elem ( const bool be_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 44 of file inf_elem_builder.C.

References _mesh, libMesh::MeshTools::bounding_box(), libMesh::out, libMesh::MeshBase::prepare_for_use(), libMesh::ParallelObject::processor_id(), and libMesh::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 && _mesh.processor_id() == 0)
    {
#ifdef DEBUG
      libMesh::out << " Determined origin for Infinite Elements:"
                   << std::endl
                   << "  ";
      origin.write_unformatted(libMesh::out);
      libMesh::out << 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 (/*skip_renumber =*/ false);

  return origin;
}

const Point libMesh::InfElemBuilder::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`
	)

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 85 of file inf_elem_builder.C.

References _mesh, libMesh::MeshTools::bounding_box(), build_inf_elem(), libMesh::Elem::build_side(), libMesh::MeshBase::elem(), libMesh::MeshBase::node(), libMesh::out, libMesh::MeshBase::prepare_for_use(), libMesh::MeshBase::print_info(), side, libMesh::START_LOG(), and libMesh::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)
        {
          libMesh::out << " Origin for Infinite Elements:" << std::endl;

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

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

  else if (be_verbose)

    {
      libMesh::out << " Origin for Infinite Elements:" << std::endl;
      libMesh::out << "  coordinates: ";
      origin.write_unformatted(libMesh::out);
      libMesh::out << 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<dof_id_type,
        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();
          libMesh::out << "Data pre-processing:" << std::endl
                       << " convert the <int,int> list to a Node* list..."
                       << std::endl;
        }

      // First use a std::vector<dof_id_type> 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<> 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<dof_id_type> 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<dof_id_type,unsigned int> >::iterator face_it = inner_faces.begin();
      const std::set< std::pair<dof_id_type,unsigned int> >::iterator face_end = inner_faces.end();
      for(; face_it!=face_end; ++face_it)
        {
          std::pair<dof_id_type,unsigned int> p = *face_it;

          // build a full-ordered side element to get _all_ the base nodes
          UniquePtr<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!
#ifndef NDEBUG
      const std::size_t ibn_size_before = inner_boundary_node_numbers.size();
#endif
      std::sort (inner_boundary_node_numbers.begin(), inner_boundary_node_numbers.end());
      std::vector<dof_id_type>::iterator unique_end =
        std::unique (inner_boundary_node_numbers.begin(), inner_boundary_node_numbers.end());

      std::size_t unique_size = std::distance(inner_boundary_node_numbers.begin(), unique_end);
      libmesh_assert_less_equal (unique_size, 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<dof_id_type>::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)
        libMesh::out << "  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 (/*skip_renumber =*/ false);

  return origin;
}

void libMesh::InfElemBuilder::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< dof_id_type, unsigned int > > *	inner_faces = `NULL`
	)		`[private]`

Build infinite elements atop a volume-based mesh. Actual implementation.

Definition at line 281 of file inf_elem_builder.C.

{
  if (be_verbose)
    {
#ifdef DEBUG
      libMesh::out << " Building Infinite Elements:" << std::endl;
      libMesh::out << "  updating element neighbor tables..." << std::endl;
#else
      libMesh::out << " 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<dof_id_type,unsigned int> > faces;
  std::set< std::pair<dof_id_type,unsigned int> > ofaces;

  // A set for storing node numbers on the outer faces.
  std::set<dof_id_type> 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
  int max_r_node = -1;

#ifdef DEBUG
  if (be_verbose)
    {
      libMesh::out << "  collecting boundary sides";
      if (x_sym || y_sym || z_sym)
        libMesh::out << ", skipping sides in symmetry planes..." << std::endl;
      else
        libMesh::out << "..." << 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
                UniquePtr<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.
  if (max_r_node > 0)
    onodes.insert(max_r_node);


  {
    std::set< std::pair<dof_id_type,unsigned int> >::iterator face_it = faces.begin();
    unsigned int facesfound=0;
    while (face_it != faces.end()) {

      std::pair<dof_id_type, unsigned int> p;
      p = *face_it;

      // This has to be a full-ordered side element,
      // since we need the correct n_nodes,
      UniquePtr<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();
        }

    }
  }


#ifdef DEBUG
  if (be_verbose)
    libMesh::out << "  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<dof_id_type, Node *> outer_nodes;

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

  // for each boundary node, add an outer_node with
  // double distance from origin.
  std::set<dof_id_type>::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);
          dof_id_type 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
  dof_id_type n_conventional_elem = this->_mesh.n_elem();
#endif


  // build Elems based on boundary side type
  std::set< std::pair<dof_id_type,unsigned int> >::iterator face_it = ofaces.begin();
  for( ; face_it != ofaces.end(); ++face_it)
    {
      // Shortcut to the pair being iterated over
      std::pair<dof_id_type,unsigned int> p = *face_it;

      // build a full-ordered side element to get the base nodes
      UniquePtr<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:
          libMesh::out << "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 \p 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);
    } // for


#ifdef DEBUG
  _mesh.libmesh_assert_valid_parallel_ids();

  if (be_verbose)
    libMesh::out << "  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");
}

Member Data Documentation

MeshBase& libMesh::InfElemBuilder::_mesh [private]

Reference to the mesh we're building infinite elements for.

Definition at line 126 of file inf_elem_builder.h.

Referenced by build_inf_elem().

The documentation for this class was generated from the following files:

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation