psurfaceboundary.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef PSURFACE_BOUNDARY_HH
#define PSURFACE_BOUNDARY_HH
 
#include "../../../common/gridfactory.hh"
 
#if HAVE_PSURFACE
#include <psurface/PSurface.h>
#include "psurface/AmiraMeshIO.h"
#if HAVE_PSURFACE_2_0
#include <psurface/Hdf5IO.h>
#endif
 
#if HAVE_AMIRAMESH
#include <amiramesh/AmiraMesh.h>
#endif
 
 
namespace Dune {
 
  template <int dim>
  class PSurfaceBoundary
  {
    dune_static_assert((dim==1 or dim==2), "PSurfaceBoundaries can only have dimensions 1 or 2!");
 
  public:
 
    class PSurfaceBoundarySegment : public Dune::BoundarySegment<dim+1>
    {
    public:
 
      PSurfaceBoundarySegment(const shared_ptr<PSurfaceBoundary<dim> >& psurfaceBoundary, int segment)
        : psurfaceBoundary_(psurfaceBoundary),
          segment_(segment)
      {}
 
      virtual Dune::FieldVector<double, dim+1> operator()(const Dune::FieldVector<double,dim>& local) const {
 
        Dune::FieldVector<double, dim+1> result;
 
        // Transform local to barycentric coordinates
        PSURFACE_NAMESPACE StaticVector<float,dim> barCoords;
 
        if (dim==2) {
          barCoords[0] = 1 - local[0] - local[1];
          barCoords[1] = local[0];
        } else {          // dim==1
          barCoords[0] = 1 - local[0];
        }
 
        PSURFACE_NAMESPACE StaticVector<float,dim+1> r;
 
        if (!psurfaceBoundary_->getPSurfaceObject()->positionMap(segment_, barCoords, r))
          DUNE_THROW(Dune::GridError, "psurface::positionMap returned error code");
 
        for (int i=0; i<dim+1; i++)
          result[i] = r[i];
 
        return result;
      }
 
      shared_ptr<PSurfaceBoundary<dim> > psurfaceBoundary_;
      int segment_;
    };
 
 
 
    PSurfaceBoundary(PSURFACE_NAMESPACE PSurface<dim,float>* psurface)
      : psurface_(psurface)
    {}
 
    PSURFACE_NAMESPACE PSurface<dim,float>* getPSurfaceObject()
    {
      return psurface_.get();
    }
 
    static shared_ptr<PSurfaceBoundary<dim> > read(const std::string& filename)
    {
      PSURFACE_NAMESPACE PSurface<dim,float>* newDomain;
 
#if HAVE_PSURFACE_2_0
      // Try to read the file as an hdf5 file
      if (filename.find(".h5")==filename.length()-3) {
        newDomain = psurface::Hdf5IO<float,dim>::read(filename);
        if (newDomain)
          return make_shared<PSurfaceBoundary<dim> >(newDomain);
      }
#endif
 
#if HAVE_AMIRAMESH
      std::auto_ptr<AmiraMesh> am(AmiraMesh::read(filename.c_str()));
 
      if (!am.get())
        DUNE_THROW(IOError, "An error has occured while reading " << filename);
 
      newDomain
        = (PSURFACE_NAMESPACE PSurface<dim,float>*) PSURFACE_NAMESPACE AmiraMeshIO<float>::readAmiraMesh(am.get(), filename.c_str());
 
      if (!newDomain)
        DUNE_THROW(IOError, "An error has occured while reading " << filename);
 
      return make_shared<PSurfaceBoundary<dim> >(newDomain);
#else
      DUNE_THROW(IOError, "The given file is not in a supported format!");
#endif
    }
 
  private:
 
    std::auto_ptr<PSURFACE_NAMESPACE PSurface<dim,float> > psurface_;
 
  };
 
}
 
#endif // #if HAVE_PSURFACE
#endif // #ifndef PSURFACE_BOUNDARY_HH