viennagrid/algorithm/centroid.hpp

Go to the documentation of this file.

#ifndef VIENNAGRID_ALGORITHM_CENTROID_HPP
#define VIENNAGRID_ALGORITHM_CENTROID_HPP

/* =======================================================================
   Copyright (c) 2011-2014, Institute for Microelectronics,
                            Institute for Analysis and Scientific Computing,
                            TU Wien.

                            -----------------
                     ViennaGrid - The Vienna Grid Library
                            -----------------

   License:      MIT (X11), see file LICENSE in the base directory
======================================================================= */

#include <iostream>
#include <sstream>
#include <string>
#include <stdexcept>

#include "viennagrid/forwards.hpp"
#include "viennagrid/topology/all.hpp"
#include "viennagrid/algorithm/volume.hpp"    //for mesh/segment centroid
#include "viennagrid/accessor.hpp"

namespace viennagrid
{
  namespace detail
  {
    //
    // Calculation of centroid
    //
    template <typename PointAccessorT, typename ElementT>
    typename PointAccessorT::value_type
    centroid(PointAccessorT const accessor, ElementT const & cell, viennagrid::triangle_tag)
    {
      typedef typename PointAccessorT::value_type PointType;
      typedef typename viennagrid::result_of::coord<PointType>::type    CoordType;

      typedef typename viennagrid::result_of::const_element_range<ElementT, vertex_tag>::type         VertexOnCellRange;
      typedef typename viennagrid::result_of::iterator<VertexOnCellRange>::type            VertexOnCellIterator;

      PointType p0(0.0, 0.0);

      VertexOnCellRange vertices = viennagrid::elements<vertex_tag>(cell);
      for (VertexOnCellIterator vocit = vertices.begin();
           vocit != vertices.end();
           ++vocit)
        p0 += accessor(*vocit);

      p0 /= static_cast<CoordType>(vertices.size());

      return p0;
    }

    //tetrahedron can reuse the algorithm defined for a triangle
    template <typename PointAccessorT, typename ElementT>
    typename PointAccessorT::value_type
    centroid(PointAccessorT const accessor, ElementT const & cell, viennagrid::tetrahedron_tag)
    {
      return centroid(accessor, cell, viennagrid::triangle_tag());
    }

    //
    // Note: This works for rectangles only, but not for general quadrilateral
    //
    template <typename PointAccessorT, typename ElementT>
    typename PointAccessorT::value_type
    centroid(PointAccessorT const accessor, ElementT const & cell, viennagrid::quadrilateral_tag)
    {
      return centroid(accessor, cell, viennagrid::triangle_tag());
    }

    //
    // Note: This works for cuboids only, but not for general hexahedra
    //
    template <typename PointAccessorT, typename ElementT>
    typename PointAccessorT::value_type
    centroid(PointAccessorT const accessor, ElementT const & cell, viennagrid::hexahedron_tag)
    {
      return centroid(accessor, cell, viennagrid::triangle_tag());
    }


    //a line can reuse the algorithm defined for a triangle
    template <typename PointAccessorT, typename ElementT>
    typename PointAccessorT::value_type
    centroid(PointAccessorT const accessor, ElementT const & cell, viennagrid::simplex_tag<1>)
    {
      return centroid(accessor, cell, viennagrid::triangle_tag());
    }


    template <typename PointAccessorT, typename ElementT>
    typename PointAccessorT::value_type
    centroid(PointAccessorT const accessor, ElementT const & cell, viennagrid::hypercube_tag<1>)
    {
      return centroid(accessor, cell, viennagrid::triangle_tag());
    }


    //a point is degenerate and returns its location
    template <typename PointAccessorT, typename ElementT>
    typename PointAccessorT::value_type
    centroid(PointAccessorT const accessor, ElementT const & cell, viennagrid::vertex_tag)
    {
      return accessor(cell);
    }



    template <typename ElementTOrTag, typename MeshSegmentHandleType, typename PointAccessorT>
    typename viennagrid::result_of::point<MeshSegmentHandleType>::type
    centroid_mesh(MeshSegmentHandleType const & mesh_obj, PointAccessorT const point_accessor)
    {
      typedef typename viennagrid::result_of::element_tag<ElementTOrTag>::type            ElementTag;
      typedef typename viennagrid::result_of::point<MeshSegmentHandleType>::type          PointType;
      typedef typename viennagrid::result_of::const_element_range<MeshSegmentHandleType,
                                                                ElementTag>::type         CellRange;
      typedef typename viennagrid::result_of::iterator<CellRange>::type                   CellIterator;

      PointType result = 0;
      double volume = 0;

      CellRange cells = viennagrid::elements<ElementTag>(mesh_obj);
      for (CellIterator cit = cells.begin(); cit != cells.end(); ++cit)
      {
        double vol_cell = viennagrid::volume( point_accessor, *cit );
        result += vol_cell * centroid( point_accessor, *cit);
        volume += vol_cell;
      }

      return result / volume;
    }

  } //namespace detail



  template <typename PointAccessorT, typename ElementT>
  typename PointAccessorT::value_type
  centroid(PointAccessorT const accessor, ElementT const & element)
  {
    return detail::centroid( accessor, element, typename ElementT::tag());
  }

  template <typename ElementTag, typename WrappedConfigT>
  typename viennagrid::result_of::point< viennagrid::element<ElementTag,WrappedConfigT> >::type
  centroid(viennagrid::element<ElementTag,WrappedConfigT> const & element)
  {
    return detail::centroid( default_point_accessor(element), element, ElementTag());
  }



  template<typename ElementTOrTagT, typename WrappedConfigT, typename PointAccessorT>
  typename viennagrid::result_of::point< mesh<WrappedConfigT> >::type
  centroid(mesh<WrappedConfigT> const & mesh_obj, PointAccessorT const point_accessor)
  {
    return detail::centroid_mesh<ElementTOrTagT>(mesh_obj, point_accessor);
  }

  template<typename WrappedConfigT, typename PointAccessorT>
  typename viennagrid::result_of::point< mesh<WrappedConfigT> >::type
  centroid(mesh<WrappedConfigT> const & mesh_obj, PointAccessorT const point_accessor)
  {
    typedef typename viennagrid::result_of::cell_tag< mesh<WrappedConfigT> >::type CellTag;
    return detail::centroid_mesh<CellTag>(mesh_obj, point_accessor);
  }


  template<typename ElementTOrTagT, typename WrappedConfigT>
  typename viennagrid::result_of::point< mesh<WrappedConfigT> >::type
  centroid(mesh<WrappedConfigT> const & mesh_obj)
  {
    return centroid<ElementTOrTagT>(mesh_obj, default_point_accessor(mesh_obj));
  }

  template<typename WrappedConfigT>
  typename viennagrid::result_of::point< mesh<WrappedConfigT> >::type
  centroid(mesh<WrappedConfigT> const & mesh_obj)
  {
    typedef typename viennagrid::result_of::cell_tag< mesh<WrappedConfigT> >::type CellTag;
    return centroid<CellTag>(mesh_obj, default_point_accessor(mesh_obj));
  }


  template<typename ElementTOrTagT, typename SegmentationT, typename PointAccessorT>
  typename viennagrid::result_of::point< segment_handle<SegmentationT> >::type
  centroid(segment_handle<SegmentationT> const & segment, PointAccessorT const point_accessor)
  {
    return detail::centroid_mesh<ElementTOrTagT>(segment, point_accessor);
  }

  template<typename SegmentationT, typename PointAccessorT>
  typename viennagrid::result_of::point< segment_handle<SegmentationT> >::type
  centroid(segment_handle<SegmentationT> const & segment, PointAccessorT const point_accessor)
  {
    typedef typename viennagrid::result_of::cell_tag< segment_handle<SegmentationT> >::type CellTag;
    return detail::centroid_mesh<CellTag>(segment, point_accessor);
  }


  template<typename ElementTOrTagT, typename SegmentationT>
  typename viennagrid::result_of::point< segment_handle<SegmentationT> >::type
  centroid(segment_handle<SegmentationT> const & segment)
  {
    return centroid<ElementTOrTagT>(segment, default_point_accessor(segment));
  }

  template<typename SegmentationT>
  typename viennagrid::result_of::point< segment_handle<SegmentationT> >::type
  centroid(segment_handle<SegmentationT> const & segment)
  {
    typedef typename viennagrid::result_of::cell_tag< segment_handle<SegmentationT> >::type CellTag;
    return centroid<CellTag>(segment, default_point_accessor(segment));
  }
} //namespace viennagrid
#endif