SpectralMesh.cpp

Go to the documentation of this file.

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

//  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.cpp,v 1.8 2008/03/10 21:58:10 delpinux Exp $


// This class allow to Manipulate 3D Scalar Variable.

#include <fstream>
#include <cmath>

#include <ConnectivityBuilder.hpp>

#include <CartesianHexahedron.hpp>
#include <SpectralMesh.hpp>
#include <Scene.hpp>

#include <FacesBuilder.hpp>
#include <EdgesBuilder.hpp>

#include <GaussLobattoManager.hpp>

#include <SpectralConformTransformation.hpp>

void SpectralMesh::
buildEdges()
{
  EdgesBuilder<SpectralMesh> edgesBuilder(*this);
  __edgesSet = edgesBuilder.edgesSet();
}

void SpectralMesh::
buildFaces()
{
  FacesBuilder<SpectralMesh> facesBuilder(*this);
  __facesSet = facesBuilder.facesSet();
}

SpectralMesh::
SpectralMesh(const Structured3DMeshShape& degrees,
             ReferenceCounting<VerticesCorrespondance> correspondances)
  : Mesh(Mesh::spectralMesh,
         Mesh::volume,
         new VerticesSet((degrees.shape().nx()+2)*(degrees.shape().ny()+2)*(degrees.shape().nz()+2)),
         correspondances),
    __degrees(TinyVector<3,size_t>(degrees.shape().nx(),degrees.shape().ny(),degrees.shape().nz())),
    __verticesShape(TinyVector<3,size_t>(degrees.shape().nx()+2,degrees.shape().ny()+2,degrees.shape().nz()+2), degrees.a(), degrees.b()),
    __cellShape(TinyVector<3,size_t>(degrees.shape().nx()+1,degrees.shape().ny()+1,degrees.shape().nz()+1)),
    __cells((degrees.shape().nx()+1)*(degrees.shape().ny()+1)*(degrees.shape().nz()+1)),
    __facesSet(0),
    __connectivity(*this)
{
  // w[ijk][01] variable respresent the coefficient needed to 
  // compute vertices position between a and b.
  // for example wj1 is the weight computed for the jth point in the
  // 'y' direction for v1
  // const size_t degreeX = degrees.shape().nx();
  //   const size_t degreeY = degrees.shape().ny();
  //   const size_t degreeZ = degrees.shape().nz();
  
  const size_t degreeX = degrees.nx();
  const size_t degreeY = degrees.ny();
  const size_t degreeZ = degrees.nz();
  
  const size_t nx_1 = __cellShape.nx();
  const size_t ny_1 = __cellShape.ny();
  const size_t nz_1 = __cellShape.nz();
  Interval intervalX(__verticesShape.a()[0],__verticesShape.b()[0]) ;
  Interval intervalY(__verticesShape.a()[1],__verticesShape.b()[1]) ;
  Interval intervalZ(__verticesShape.a()[2],__verticesShape.b()[2]) ;
  const GaussLobatto& gaussLobattoX = GaussLobattoManager::instance().get(degreeX+1); 
  const GaussLobatto& gaussLobattoY = GaussLobattoManager::instance().get(degreeY+1);
  const GaussLobatto& gaussLobattoZ = GaussLobattoManager::instance().get(degreeZ+1);

  SpectralConformTransformation transformationX(intervalX);
  SpectralConformTransformation transformationY(intervalY);
  SpectralConformTransformation transformationZ(intervalZ);
  for (size_t i=0; i<__verticesShape.nx(); i++) {
    for (size_t j=0; j<__verticesShape.ny(); j++) {
      for (size_t k=0; k<__verticesShape.nz(); k++) {
        Vertex& v = vertex(i,j,k);

        v.x() = transformationX(gaussLobattoX(i));
        v.y() = transformationY(gaussLobattoY(j));
        v.z() = transformationZ(gaussLobattoZ(k));
      }
    }
  }

  // Generating Mesh's cells.
  // As Cells generation needs Vertices's references
  // This must be done after.

  size_t n=0;
  for (size_t i=0; i<nx_1; i++) {
    for (size_t j=0; j<ny_1; j++) {
      for (size_t k=0; k<nz_1; k++) {
        __cells[n] = CartesianHexahedron(vertex(  i,  j,  k),
                                         vertex(i+1,  j,  k),
                                         vertex(i+1,j+1,  k),
                                         vertex(  i,j+1,  k),
                                         vertex(  i,  j,k+1),
                                         vertex(i+1,  j,k+1),
                                         vertex(i+1,j+1,k+1),
                                         vertex(  i,j+1,k+1));
        n++;
      }
    }
  }

  Vector<Quadrangle>* pQuadrangles
    = new Vector<Quadrangle>(2*ny_1*nz_1+2*nx_1*nz_1+2*nx_1*ny_1);
  Vector<Quadrangle>& quadrangles = *pQuadrangles;

  size_t currentCell = 0;
  for (size_t j=0; j<ny_1; j++) {
    for (size_t k=0; k<nz_1; k++) {
      const Vertex& V0 = vertex(0,  j,  k);
      const Vertex& V1 = vertex(0,  j,k+1);
      const Vertex& V2 = vertex(0,j+1,k+1);
      const Vertex& V3 = vertex(0,j+1,  k);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 0);
      Q.setMother(&(cell(0,j,k)),4);

      currentCell++;
    }
  }

  for (size_t j=0; j<ny_1; j++) {
    for (size_t k=0; k<nz_1; k++) {
      const Vertex& V0 = vertex(nx_1,  j,  k);
      const Vertex& V1 = vertex(nx_1,j+1,  k);
      const Vertex& V2 = vertex(nx_1,j+1,k+1);
      const Vertex& V3 = vertex(nx_1,  j,k+1);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 1);
      Q.setMother(&(cell(nx_1-1,j,k)),2);

      currentCell++;
    }
  }

  for (size_t i=0; i<nx_1; i++) {
    for (size_t k=0; k<nz_1; k++) {
      const Vertex& V0 = vertex(  i,0,  k);
      const Vertex& V1 = vertex(i+1,0,  k);
      const Vertex& V2 = vertex(i+1,0,k+1);
      const Vertex& V3 = vertex(  i,0,k+1);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 2);
      Q.setMother(&(cell(i,0,k)),1);

      currentCell++;
    }
  }

  for (size_t i=0; i<nx_1; i++) {
    for (size_t k=0; k<nz_1; k++) {
      const Vertex& V0 = vertex(  i,ny_1,  k);
      const Vertex& V1 = vertex(  i,ny_1,k+1);
      const Vertex& V2 = vertex(i+1,ny_1,k+1);
      const Vertex& V3 = vertex(i+1,ny_1,  k);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 3);
      Q.setMother(&(cell(i,ny_1-1,k)),3);

      currentCell++;
    }
  }

  for (size_t i=0; i<nx_1; i++) {
    for (size_t j=0; j<ny_1; j++) {
      const Vertex& V0 = vertex(  i,  j,0);
      const Vertex& V1 = vertex(  i,j+1,0);
      const Vertex& V2 = vertex(i+1,j+1,0);
      const Vertex& V3 = vertex(i+1,  j,0);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 4);
      Q.setMother(&(cell(i,j,0)),0);

      currentCell++;
    }
  }


  for (size_t i=0; i<nx_1; i++) {
    for (size_t j=0; j<ny_1; j++) {
      const Vertex& V0 = vertex(  i,  j,nz_1);
      const Vertex& V1 = vertex(i+1,  j,nz_1);
      const Vertex& V2 = vertex(i+1,j+1,nz_1);
      const Vertex& V3 = vertex(  i,j+1,nz_1);
      Quadrangle& Q = quadrangles[currentCell];

      Q = Quadrangle(V0, V1, V2, V3, 5);
      Q.setMother(&(cell(i,j,nz_1-1)),5);
      currentCell++;
    }
  }
  __borderMesh = new SurfaceMeshOfQuadrangles(__verticesSet,
                                              __verticesCorrespondance,
                                              pQuadrangles);
  __borderMesh->setBackgroundMesh(this);
}

---