SpectralMesh.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 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: SpectralMesh.hpp,v 1.9 2008/03/04 00:17:57 delpinux Exp $

// This class allow 3D Structured Mesh management.

#ifndef  SPECTRAL_MESH_HPP
#define  SPECTRAL_MESH_HPP

#include <Mesh.hpp>

#include <TinyVector.hpp>

#include <CartesianHexahedron.hpp>

#include <SurfaceMeshOfQuadrangles.hpp>
#include <FacesSet.hpp>

#include <Connectivity.hpp>

#include <Structured3DMeshShape.hpp>

#include <Index.hpp>

class MeshOfHexahedra;

class SpectralMesh
  : public Mesh
{
public:
  typedef CartesianHexahedron CellType;
  typedef Quadrangle FaceType;

  typedef MeshOfHexahedra Transformed;
  typedef SurfaceMeshOfQuadrangles BorderMeshType;

  enum {
    family = Mesh::volume
  };
private:
  const TinyVector<3,size_t> __degrees;

  const Structured3DMeshShape __verticesShape;

  const Array3DShape __cellShape;

  Vector<CartesianHexahedron> __cells;

  ReferenceCounting<BorderMeshType>
  __borderMesh; 
  ReferenceCounting<FacesSet<FaceType> >
  __facesSet;  
  Connectivity<SpectralMesh> __connectivity; 
  SpectralMesh(const SpectralMesh&);

public:
  std::string typeName() const
  {
    return "spectral mesh";
  }

  const TinyVector<3,size_t>& degrees() const
  {
    return __degrees;
  }

  size_t degree(const size_t& i) const
  {
    return __degrees[i];
  }

  size_t dofNumber(const size_t& i,
                   const size_t& j,
                   const size_t& k) const
  {
    return __cellShape(i, j, k);
  }

  void buildEdges();

  bool hasBorderMesh() const
  {
    return (__borderMesh != 0);
  }

  ConstReferenceCounting<BorderMeshType> borderMesh() const
  {
    return __borderMesh;
  }

  ConstReferenceCounting<Mesh> borderBaseMesh() const
  {
    return static_cast<const BorderMeshType*>(__borderMesh);
  }

  void buildFaces();

  inline bool hasFaces() const
  {
    return __facesSet != 0;
  }

  const FaceType& face(const size_t& i) const
  {
    return (*__facesSet)[i];
  }

  size_t faceNumber(const FaceType& f) const
  {
    return (*__facesSet).number(f);
  }

  inline const size_t& numberOfFaces() const
  {
    return (*__facesSet).numberOfFaces();
  }

  inline const size_t& numberOfCells() const
  {
    return __cells.size();
  }

  size_t cellNumber(const CartesianHexahedron& h) const
  {
    return __cells.number(h);
  }

  typedef Mesh::T_iterator<SpectralMesh, CartesianHexahedron> iterator;
  typedef Mesh::T_iterator<const SpectralMesh, const CartesianHexahedron> const_iterator;

  inline SpectralMesh::const_iterator find(const double& x,
                                           const double& y,
                                           const double& z) const;

  inline SpectralMesh::const_iterator find(const TinyVector<3>& X) const
  {
    return find(X[0],X[1],X[2]);
  }

  inline Vertex& vertex(const size_t& i,
                        const size_t& j,
                        const size_t& k);

  inline const Vertex& vertex(const size_t& i,
                              const size_t& j,
                              const size_t& k) const;

  inline CartesianHexahedron& cell(const size_t& i,
                                   const size_t& j,
                                   const size_t& k);

  inline const CartesianHexahedron& cell(const size_t& i,
                                         const size_t& j,
                                         const size_t& k) const;

  inline CartesianHexahedron& cell(const size_t& i)
  {
    return __cells[i];
  }

  const Connectivity<SpectralMesh>& connectivity() const
  {
    return __connectivity;
  }

  Connectivity<SpectralMesh>& connectivity()
  {
    return __connectivity;
  }

  inline const CartesianHexahedron& cell(const size_t& i) const
  {
    return __cells[i];
  }

  const Structured3DMeshShape& shape()const
  {
    return __verticesShape;
  }

  inline Vertex& vertex(const size_t& i)
  {
    return Mesh::vertex(i);
  }

  inline const Vertex& vertex(const size_t& i) const
  {
    return Mesh::vertex(i);
  }

  inline Vertex& vertex(const Index& I);

  inline const Vertex& vertex(const Index& I) const;

  inline CartesianHexahedron& cell(const Index& I);

  inline const CartesianHexahedron& cell(const Index& I) const;

  inline Index cellIndex(const Vertex& V) const;

  inline Index cellIndex(const TinyVector<3>& V) const;

  inline Index vertexIndex(const Vertex& V) const;

  inline bool inside(const real_t& x, const real_t& y, const real_t& z) const
  {
    return ((x>=__verticesShape.a(0))&&
            (x<=__verticesShape.b(0))&&
            (y>=__verticesShape.a(1))&&
            (y<=__verticesShape.b(1))&&
            (z>=__verticesShape.a(2))&&
            (z<=__verticesShape.b(2)));
  }

  inline bool inside(const TinyVector<3>& p) const
  {
    return this->inside(p[0], p[1], p[2]);
  }

  SpectralMesh(const Structured3DMeshShape& s3dM,
               ReferenceCounting<VerticesCorrespondance> correspondance);

};

inline Vertex& SpectralMesh::vertex(const size_t& i,
                                    const size_t& j,
                                    const size_t& k)
{
  return (*__verticesSet)[__verticesShape(i, j, k)];
}

inline const Vertex& SpectralMesh::vertex(const size_t& i,
                                          const size_t& j,
                                          const size_t& k) const
{
  return (*__verticesSet)[__verticesShape(i, j, k)];
}

inline CartesianHexahedron& SpectralMesh::cell(const size_t& i,
                                               const size_t& j,
                                               const size_t& k)
{
  return (__cells[__cellShape(i, j, k)]);
}

inline const CartesianHexahedron& SpectralMesh::cell(const size_t& i,
                                                     const size_t& j,
                                                     const size_t& k) const
{
  return (__cells[__cellShape(i, j, k)]);
}
inline Vertex& SpectralMesh::vertex(const Index& I)
{
  return (*__verticesSet)[__verticesShape(I[0], I[1], I[2])];
}

inline const Vertex& SpectralMesh::vertex(const Index& I) const
{
  return (*__verticesSet)[__verticesShape(I[0], I[1], I[2])];
}

inline CartesianHexahedron& SpectralMesh::cell(const Index& I)
{
  return (__cells[__cellShape(I[0], I[1], I[2])]);
}

inline const CartesianHexahedron& SpectralMesh::cell(const Index& I) const
{
  return (__cells[__cellShape(I[0], I[1], I[2])]);
}

inline Index SpectralMesh::cellIndex(const TinyVector<3>& V) const
{
  const real_t& x = V[0];
  const real_t& y = V[1];
  const real_t& z = V[2];

  const real_t& x0 = __verticesShape.a(0);
  const real_t& y0 = __verticesShape.a(1);
  const real_t& z0 = __verticesShape.a(2);

  int i = int(std::floor((x - x0)/__verticesShape.hx()));
  int j = int(std::floor((y - y0)/__verticesShape.hy()));
  int k = int(std::floor((z - z0)/__verticesShape.hz()));

  i = (i<0) ? 0 : i;
  j = (j<0) ? 0 : j;
  k = (k<0) ? 0 : k;

  int nx_1 = __cellShape.nx() - 1;
  int ny_1 = __cellShape.ny() - 1;
  int nz_1 = __cellShape.nz() - 1;

  i = (i>nx_1) ? nx_1 : i;
  j = (j>ny_1) ? ny_1 : j;
  k = (k>nz_1) ? nz_1 : k;

  Index I(i,j,k);

  return I;
}

inline Index SpectralMesh::vertexIndex(const Vertex& V) const
{
  const real_t& x = V.x();
  const real_t& y = V.y();
  const real_t& z = V.z();

  const real_t& x0 = __verticesShape.a(0);
  const real_t& y0 = __verticesShape.a(1);
  const real_t& z0 = __verticesShape.a(2);

  int i = int((x - x0)/__verticesShape.hx()+0.5);
  int j = int((y - y0)/__verticesShape.hy()+0.5);
  int k = int((z - z0)/__verticesShape.hz()+0.5);

  bool foundCell = true;
  if (((i<0)||(i>static_cast<int>(__verticesShape.nx()-1)))
      &&((j<0)||(j>static_cast<int>(__verticesShape.ny()-1)))
      &&((k<0)||(k>static_cast<int>(__verticesShape.nz()-1)))) {
    foundCell = false;
  }

  Index I(i,j,k);

  return I;
}

inline SpectralMesh::const_iterator 
SpectralMesh::find(const double& x,
                   const double& y,
                   const double& z) const
{
  bool foundCell = inside(x,y,z);

  if (!foundCell)
    return SpectralMesh::const_iterator(*this,
                                        SpectralMesh::const_iterator::End);

  const real_t& x0 = __verticesShape.a(0);
  const real_t& y0 = __verticesShape.a(1);
  const real_t& z0 = __verticesShape.a(2);

  int i = int(std::floor((x - x0)/__verticesShape.hx()));
  int j = int(std::floor((y - y0)/__verticesShape.hy()));
  int k = int(std::floor((z - z0)/__verticesShape.hz()));

  {
    int nx_1 = __cellShape.nx() - 1;
    int ny_1 = __cellShape.ny() - 1;
    int nz_1 = __cellShape.nz() - 1;

    // this transformation is to operate if the (x,y,z) is a point on
    // one of the faces of the box.

    i = (i>nx_1) ? nx_1 : i;
    j = (j>ny_1) ? ny_1 : j;
    k = (k>nz_1) ? nz_1 : k;
  }

  return SpectralMesh::const_iterator(*this,
                                      __cellShape(i, j, k));
}

#endif // SPECTRAL_MESH_HPP

---