MeshTransformer.cpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2005 Stéphane Del Pino

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program 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 General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: MeshTransformer.cpp,v 1.4 2008/03/04 00:17:57 delpinux Exp $

#include <MeshTransformer.hpp>

#include <FieldOfScalarFunction.hpp>

#include <Structured3DMesh.hpp>
#include <MeshOfTetrahedra.hpp>
#include <MeshOfHexahedra.hpp>

#include <SpectralMesh.hpp>

#include <SurfaceMeshOfTriangles.hpp>
#include <SurfaceMeshOfQuadrangles.hpp>

template <typename CellType,
          typename OriginalCellType>
CellType MeshTransformer::
__getCell(VerticesSet& V,
          const OriginalCellType& originalCell)
{
  throw ErrorHandler(__FILE__,__LINE__,
                     "non-supported cell type",
                     ErrorHandler::unexpected);
}

template <>
Hexahedron MeshTransformer::
__getCell<Hexahedron,
          CartesianHexahedron>(VerticesSet& V,
                               const CartesianHexahedron& H)
{
  const size_t i0 = __input->vertexNumber(H(0));
  const size_t i1 = __input->vertexNumber(H(1));
  const size_t i2 = __input->vertexNumber(H(2));
  const size_t i3 = __input->vertexNumber(H(3));
  const size_t i4 = __input->vertexNumber(H(4));
  const size_t i5 = __input->vertexNumber(H(5));
  const size_t i6 = __input->vertexNumber(H(6));
  const size_t i7 = __input->vertexNumber(H(7));

  return Hexahedron(V[i0],V[i1],V[i2],V[i3],
                    V[i4],V[i5],V[i6],V[i7],
                    H.reference());
}

template <>
Tetrahedron MeshTransformer::
__getCell<Tetrahedron,
          Tetrahedron>(VerticesSet& V,
                       const Tetrahedron& T)
{
  const size_t i0 = __input->vertexNumber(T(0));
  const size_t i1 = __input->vertexNumber(T(1));
  const size_t i2 = __input->vertexNumber(T(2));
  const size_t i3 = __input->vertexNumber(T(3));

  return Tetrahedron(V[i0],V[i1],V[i2],V[i3],
                     T.reference());
}

template <>
Triangle MeshTransformer::
__getCell<Triangle,
          Triangle>(VerticesSet& V,
                    const Triangle& T)
{
  const size_t i0 = __input->vertexNumber(T(0));
  const size_t i1 = __input->vertexNumber(T(1));
  const size_t i2 = __input->vertexNumber(T(2));

  return Triangle(V[i0],V[i1],V[i2],
                  T.reference());
}

template <>
Quadrangle MeshTransformer::
__getCell<Quadrangle,
          Quadrangle>(VerticesSet& V,
                      const Quadrangle& Q)
{
  const size_t i0 = __input->vertexNumber(Q(0));
  const size_t i1 = __input->vertexNumber(Q(1));
  const size_t i2 = __input->vertexNumber(Q(2));
  const size_t i3 = __input->vertexNumber(Q(3));

  return Quadrangle(V[i0],V[i1],V[i2],V[i3],
                    Q.reference());
}

template <typename MeshType>
void MeshTransformer::
__transformSurface()
{
  const MeshType& m
    = static_cast<const MeshType&>(*__input);

  typedef typename MeshType::CellType OldCellType;

  typedef typename MeshType::Transformed NewMeshType;
  typedef typename NewMeshType::CellType NewCellType;

  ReferenceCounting<VerticesSet> v = new VerticesSet(m.numberOfVertices());

  VerticesSet& vertices = *v;

  if (__field->numberOfComponents() != 3) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "Transformation field require 3 components, you provided:"
                       +stringify(*__field),
                       ErrorHandler::normal);
  }

  const ScalarFunctionBase& f1 = *__field->function(0);
  const ScalarFunctionBase& f2 = *__field->function(1);
  const ScalarFunctionBase& f3 = *__field->function(2);

  for (size_t i=0; i<m.numberOfVertices(); ++i) {
    const Vertex& X = m.vertex(i);

    vertices[i][0] = f1(X);
    vertices[i][1] = f2(X);
    vertices[i][2] = f3(X);

    vertices[i].reference() = X.reference();
  }

  ReferenceCounting<Vector<NewCellType> > elements
    = new Vector<NewCellType>(m.numberOfCells());


  for (size_t i=0; i<m.numberOfCells(); ++i) {
    (*elements)[i]
      = __getCell<NewCellType,
                  OldCellType>(vertices, m.cell(i));
  }

  ReferenceCounting<VerticesCorrespondance> correspondances
    = new VerticesCorrespondance(m.numberOfVertices());
  __mesh = new NewMeshType(v, correspondances, elements);
}

template <typename MeshType>
void MeshTransformer::
__transformVolume()
{
  const MeshType& m
    = static_cast<const MeshType&>(*__input);

  typedef typename MeshType::CellType OldCellType;

  typedef typename MeshType::Transformed NewMeshType;
  typedef typename NewMeshType::CellType NewCellType;

  ReferenceCounting<VerticesSet> v = new VerticesSet(m.numberOfVertices());

  VerticesSet& vertices = *v;

  if (__field->numberOfComponents() != 3) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "Transformation field require 3 components, you provided:"
                       +stringify(*__field),
                       ErrorHandler::normal);
  }

  const ScalarFunctionBase& f1 = *__field->function(0);
  const ScalarFunctionBase& f2 = *__field->function(1);
  const ScalarFunctionBase& f3 = *__field->function(2);

  for (size_t i=0; i<m.numberOfVertices(); ++i) {
    const Vertex& X = m.vertex(i);

    vertices[i][0] = f1(X);
    vertices[i][1] = f2(X);
    vertices[i][2] = f3(X);

    vertices[i].reference() = X.reference();
  }

  ReferenceCounting<Vector<NewCellType> > elements
    = new Vector<NewCellType>(m.numberOfCells());


  for (size_t i=0; i<m.numberOfCells(); ++i) {
    (*elements)[i]
      = __getCell<NewCellType,
                  OldCellType>(vertices, m.cell(i));
  }

  ReferenceCounting<VerticesCorrespondance> correspondances
    = new VerticesCorrespondance(m.numberOfVertices());

  typedef
    typename NewMeshType::BorderMeshType
    NewBorderMeshType;

  ReferenceCounting<NewBorderMeshType> borderMesh=0;

  if (m.hasBorderMesh()) {
    typedef
      typename NewBorderMeshType::CellType
      NewFaceType;

    typedef
      typename MeshType::BorderMeshType
      OldBorderMeshType;

    typedef
      typename OldBorderMeshType::CellType
      OldFaceType;

    const OldBorderMeshType& s = *m.borderMesh();

    ReferenceCounting<Vector<NewFaceType> > boundaryElements
      = new Vector<NewFaceType>(s.numberOfCells());

    for (size_t i=0; i<s.numberOfCells(); ++i) {
      (*boundaryElements)[i]
        = __getCell<NewFaceType,
                    OldFaceType>(vertices,s.cell(i));
    }

    borderMesh = new NewBorderMeshType(v,correspondances,
                                       boundaryElements);
  }

  __mesh = new NewMeshType(v, correspondances,
                           elements, borderMesh);
}

void MeshTransformer::
transform()
{
  switch (__input->type()) {
  case Mesh::cartesianHexahedraMesh: {
    this->__transformVolume<Structured3DMesh>();
    break;
  }
  case Mesh::hexahedraMesh: {
    this->__transformVolume<MeshOfHexahedra>();
    break;
  }
  case Mesh::tetrahedraMesh: {
    this->__transformVolume<MeshOfTetrahedra>();
    break;
  }
  case Mesh::spectralMesh: {
    this->__transformVolume<SpectralMesh>();
    break;
  }
  case Mesh::surfaceMeshTriangles: {
    this->__transformSurface<SurfaceMeshOfTriangles>();
    break;
  }
  case Mesh::surfaceMeshQuadrangles: {
    this->__transformSurface<SurfaceMeshOfQuadrangles>();
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "not implemented",
                       ErrorHandler::unexpected);
  }
  }
}

---