Structured3DMesh.cpp

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: Structured3DMesh.cpp,v 1.16 2008/03/04 00:17:56 delpinux Exp $


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

#include <fstream>
#include <cmath>

#include <ConnectivityBuilder.hpp>

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

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

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

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

Structured3DMesh::
Structured3DMesh(const Structured3DMeshShape& SMShape,
                 ReferenceCounting<VerticesCorrespondance> correspondances)
  : Mesh(Mesh::cartesianHexahedraMesh,
         Mesh::volume,
         new VerticesSet(SMShape.numberOfVertices()),
         correspondances),
    __verticesShape(SMShape),
    __cellShape(__verticesShape.shape()-1),
    __cells(SMShape.numberOfCells()),
    __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 nx = __verticesShape.nx();
  const size_t ny = __verticesShape.ny();
  const size_t nz = __verticesShape.nz();

  const size_t nx_1 = __cellShape.nx();
  const size_t ny_1 = __cellShape.ny();
  const size_t nz_1 = __cellShape.nz();

  for (size_t i=0; i<nx; i++) {
    real_t wi0 = static_cast<real_t>(i) / static_cast<real_t>(nx_1);
    real_t wi1 = static_cast<real_t>(nx_1-i) / static_cast<real_t>(nx_1);;
    for (size_t j=0; j<ny; j++) {
      real_t wj0 = static_cast<real_t>(j) / static_cast<real_t>(ny_1);
      real_t wj1 = static_cast<real_t>(ny_1-j) / static_cast<real_t>(ny_1);
      for (size_t k=0; k<nz; k++) {
        real_t wk0 = static_cast<real_t>(k) / static_cast<real_t>(nz_1);
        real_t wk1 = static_cast<real_t>(nz_1-k)/static_cast<real_t>(nz_1);
        Vertex& v = vertex(i,j,k);
        v.x() = wi0 * __verticesShape.b(0) + wi1 * __verticesShape.a(0);
        v.y() = wj0 * __verticesShape.b(1) + wj1 * __verticesShape.a(1);
        v.z() = wk0 * __verticesShape.b(2) + wk1 * __verticesShape.a(2);
      }
    }
  }

  // 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);
}

---