ScalarDegreeOfFreedomPositionsSet.cpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001-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: ScalarDegreeOfFreedomPositionsSet.cpp,v 1.3 2007/08/21 16:25:42 delpinux Exp $

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

#include <OctreeMesh.hpp>

#include <FiniteElementTraits.hpp>

#include <DegreeOfFreedomPositionsSet.hpp>
#include <ScalarDegreeOfFreedomPositionsSet.hpp>

#include <ScalarDiscretizationTypeSpectral.hpp>

#include <ConnectivityBuilder.hpp>

class ScalarDegreeOfFreedomPositionsSet::Builder
{
private:
  ScalarDegreeOfFreedomPositionsSet&
  __dofPositionsSet;            
  template <typename FiniteElementType,
            typename MeshType>
  size_t __computeFEMNbDOF(const MeshType& mesh)
  {
    size_t nbDOF = 0;
    size_t neededConnectivities = 0;

    // adds vertices dof
    if (FiniteElementType::numberOfVertexDegreesOfFreedom > 0) {
      nbDOF
        += FiniteElementType::numberOfVertexDegreesOfFreedom
        *  mesh.numberOfVertices();
    }
    
    // adds edges dof
    if (FiniteElementType::numberOfEdgeDegreesOfFreedom > 0) {
      if (not(mesh.hasEdges())) {
        const_cast<MeshType&>(mesh).buildEdges();
      }
      if (not(mesh.connectivity().hasCellToEdges())) {
        neededConnectivities += Connectivity<MeshType>::CellToEdges;
      }

      nbDOF
        += FiniteElementType::numberOfEdgeDegreesOfFreedom
        *  mesh.numberOfEdges();
    }

    // adds faces dof
    if (FiniteElementType::numberOfFaceDegreesOfFreedom > 0) {
      if (not(mesh.hasFaces())) {
        const_cast<MeshType&>(mesh).buildFaces();
      }
      if (not(mesh.connectivity().hasCellToFaces())) {
        neededConnectivities += Connectivity<MeshType>::CellToFaces;
      }
      nbDOF
        += FiniteElementType::numberOfFaceDegreesOfFreedom
        *  mesh.numberOfFaces();
    }

    // adds cells dof
    if (FiniteElementType::numberOfVolumeDegreesOfFreedom > 0) {
      nbDOF
        += FiniteElementType::numberOfVolumeDegreesOfFreedom
        *  mesh.numberOfCells();
    }

    if (neededConnectivities > 0) {
      ConnectivityBuilder<MeshType> builder(mesh);
      builder.generates(neededConnectivities);
    }

    ffout(3) << "- number of degrees of freedom positions: " << nbDOF << '\n';
    return nbDOF;
  }

  template<typename MeshType,
           typename FiniteElementType>
  void __buildFEM(const MeshType& mesh)
  {
    typedef typename MeshType::CellType CellType;

    __dofPositionsSet.__numberOfDOFPerCell = FiniteElementType::numberOfDegreesOfFreedom;
    __dofPositionsSet.__dofNumber.resize(mesh.numberOfCells()*FiniteElementType::numberOfDegreesOfFreedom);
    // fills the correspondance table with wrong values
    __dofPositionsSet.__dofNumber = std::numeric_limits<size_t>::max();

    // Computes the number of degrees of freedom
    const size_t nbDOF = __computeFEMNbDOF<FiniteElementType>(mesh);

    __dofPositionsSet.__positions.resize(nbDOF);

    size_t begin = 0;
    size_t end = 0;

    // vertices
    if (FiniteElementType::numberOfVertexDegreesOfFreedom > 0) {
      begin = end;
      end += mesh.numberOfVertices();
      for (size_t i=begin; i < end; ++i) {
        __dofPositionsSet.__positions[i] = mesh.vertex(i);
      }
      for (size_t i=0; i<mesh.numberOfCells(); ++i) {
        const CellType& cell = mesh.cell(i);
        for (size_t j=0; j<CellType::NumberOfVertices; ++j) {
          const size_t vertexNumber = mesh.vertexNumber(cell(j));
          __dofPositionsSet.__dofNumber[i*FiniteElementType::numberOfDegreesOfFreedom+j]
            = vertexNumber;
        }
      }
    }

    // edges
    if (FiniteElementType::numberOfEdgeDegreesOfFreedom > 0) {
      begin = end;
      end += mesh.numberOfEdges();
      for (size_t i=begin; i < end; ++i) {
        const Edge& edge = mesh.edge(i-begin);
        __dofPositionsSet.__positions[i] = 0.5*(edge(0)+edge(1));
      }

      // Take into account vertices degrees of freedom in numbering
      const size_t inputShift
        = FiniteElementType::numberOfVertexDegreesOfFreedom
        * CellType::NumberOfVertices;
      const size_t outputShift
        = mesh.numberOfVertices()
        * FiniteElementType::numberOfVertexDegreesOfFreedom;

      const Connectivity<MeshType>&
        connectivity = mesh.connectivity();

      for (size_t i=0; i<mesh.numberOfCells(); ++i) {
        const CellType& cell = mesh.cell(i);
        const typename Connectivity<MeshType>::CellToEdgesType&
          cellEdges = connectivity.edges(cell);
        for (size_t j=0; j<CellType::NumberOfEdges; ++j) {
          const Edge& edge = *(cellEdges[j]);
          const size_t edgeNumber = mesh.edgeNumber(edge);
          __dofPositionsSet.__dofNumber[i*FiniteElementType::numberOfDegreesOfFreedom+inputShift+j]
            = edgeNumber + outputShift;
        }
      }
    }

    // faces
    if (FiniteElementType::numberOfFaceDegreesOfFreedom > 0) {
      begin = end;
      end += mesh.numberOfFaces();
      for (size_t i=begin; i < end; ++i) {
        typedef typename MeshType::CellType::FaceType FaceType;

        const FaceType& face = mesh.face(i-begin);
        __dofPositionsSet.__positions[i] = 0;
        for (size_t j=0; j<FaceType::NumberOfVertices; ++j) {
          __dofPositionsSet.__positions[i] += face(j);
        }
        __dofPositionsSet.__positions[i] *= 1./FaceType::NumberOfVertices;
      }

      // Take into account vertices and edges degrees of freedom in
      // numbering
      const size_t inputShift
        = FiniteElementType::numberOfVertexDegreesOfFreedom
        * CellType::NumberOfVertices
        + FiniteElementType::numberOfEdgeDegreesOfFreedom
        * CellType::NumberOfEdges;
      const size_t outputShift
        = mesh.numberOfVertices()
        * FiniteElementType::numberOfVertexDegreesOfFreedom
        + (mesh.hasEdges() ? mesh.numberOfEdges() : 0)
        * FiniteElementType::numberOfEdgeDegreesOfFreedom;

      const Connectivity<MeshType>&
        connectivity = mesh.connectivity();

      for (size_t i=0; i<mesh.numberOfCells(); ++i) {
        const CellType& cell = mesh.cell(i);
        const typename Connectivity<MeshType>::CellToFacesType&
          cellFaces = connectivity.faces(cell);

        for (size_t j=0; j<CellType::NumberOfFaces; ++j) {
          const typename CellType::FaceType& face = *(cellFaces[j]);
          const size_t faceNumber = mesh.faceNumber(face);
          __dofPositionsSet.__dofNumber[i*FiniteElementType::numberOfDegreesOfFreedom+inputShift+j]
            = faceNumber + outputShift;
        }
      }
    }

    // cells
    if (FiniteElementType::numberOfVolumeDegreesOfFreedom > 0) {
      begin = end;
      end += mesh.numberOfCells();

      for (size_t i=0; i<mesh.numberOfCells(); ++i) {
        const CellType& cell = mesh.cell(i);
        TinyVector<3, real_t> position = 0;
        for (size_t l = 0; l<MeshType::CellType::NumberOfVertices; ++l) {
          position += cell(l);
        }
        position /= MeshType::CellType::NumberOfVertices;
        __dofPositionsSet.__positions[i+begin] = position;
      }

      // Take into account vertices, edges and face degrees of freedom
      // in numbering
      const size_t inputShift
        = FiniteElementType::numberOfVertexDegreesOfFreedom
        * CellType::NumberOfVertices
        + FiniteElementType::numberOfEdgeDegreesOfFreedom
        * CellType::NumberOfEdges
        + FiniteElementType::numberOfFaceDegreesOfFreedom
        * CellType::NumberOfFaces;
      const size_t outputShift
        = mesh.numberOfVertices()
        * FiniteElementType::numberOfVertexDegreesOfFreedom
        + (mesh.hasEdges() ? mesh.numberOfEdges() : 0)
        * FiniteElementType::numberOfEdgeDegreesOfFreedom
        + (mesh.hasFaces() ? mesh.numberOfFaces() : 0)
        * FiniteElementType::numberOfFaceDegreesOfFreedom;

      for (size_t i=0; i<mesh.numberOfCells(); ++i) {
        __dofPositionsSet.__dofNumber[i*FiniteElementType::numberOfDegreesOfFreedom+inputShift]
          = i + outputShift;
      }
    }
  }

  template<typename MeshType>
  void __buildLegendre(const MeshType& mesh,
                       const ScalarDiscretizationTypeSpectral& spectralDiscretizationType)
  {
    throw ErrorHandler(__FILE__,__LINE__,
                       "Cannot build spectral degrees of freefom on '"
                       +mesh.typeName()+"'",
                       ErrorHandler::unexpected);
  }

  template <typename MeshType>
  void __build(const ScalarDiscretizationTypeBase& discretizationType,
               const MeshType& mesh);
public:
  Builder(ScalarDegreeOfFreedomPositionsSet& dofPositionsSet)
    : __dofPositionsSet(dofPositionsSet)
  {
    ;
  }

  ~Builder()
  {
    ;
  }

  void build(const ScalarDiscretizationTypeBase& discretizationType,
             const Mesh& mesh);
};


template<>
void ScalarDegreeOfFreedomPositionsSet::Builder::
__buildLegendre<SpectralMesh>(const SpectralMesh& mesh,
                              const ScalarDiscretizationTypeSpectral& spectralDiscretizationType)
{
  // All the degrees of freedom are contained in 1 spectral cell...
  const TinyVector<3,size_t>& degrees = spectralDiscretizationType.degrees();
  __dofPositionsSet.__numberOfDOFPerCell = (degrees[0]+1)*(degrees[1]+1)*(degrees[2]+1);
  __dofPositionsSet.__dofNumber.resize(__dofPositionsSet.__numberOfDOFPerCell);

  for (size_t i=0; i<__dofPositionsSet.__dofNumber.size(); ++i) {
    __dofPositionsSet.__dofNumber[i]=i;
  }

  // legendre degrees of freedom do not really have a position...
  __dofPositionsSet.__positions.resize(__dofPositionsSet.__numberOfDOFPerCell);
  __dofPositionsSet.__positions = TinyVector<3,real_t>(0,0,0);
}

template<>
void ScalarDegreeOfFreedomPositionsSet::Builder::
__buildLegendre<OctreeMesh>(const OctreeMesh& mesh,
                            const ScalarDiscretizationTypeSpectral& spectralDiscretizationType)
{
  // All the degrees of freedom are contained in 1 spectral cell...
  const TinyVector<3,size_t>& degrees = spectralDiscretizationType.degrees();
  __dofPositionsSet.__numberOfDOFPerCell = (degrees[0]+1)*(degrees[1]+1)*(degrees[2]+1);
  __dofPositionsSet.__dofNumber.resize(__dofPositionsSet.__numberOfDOFPerCell);

  for (size_t i=0; i<__dofPositionsSet.__dofNumber.size(); ++i) {
    __dofPositionsSet.__dofNumber[i]=i;
  }

  // legendre degrees of freedom do not really have a position...
  __dofPositionsSet.__positions.resize(__dofPositionsSet.__numberOfDOFPerCell);
  __dofPositionsSet.__positions = TinyVector<3,real_t>(0,0,0);
}


template <typename MeshType>
void ScalarDegreeOfFreedomPositionsSet::Builder::
__build(const ScalarDiscretizationTypeBase& discretizationType,
        const MeshType& mesh)
{
  typedef typename MeshType::CellType CellType;
  switch(discretizationType.type()) {
  case ScalarDiscretizationTypeBase::lagrangianFEM0: {
    typedef
      typename FiniteElementTraits<CellType,
      ScalarDiscretizationTypeBase::lagrangianFEM0>::Type
      FiniteElementType;
    __buildFEM<MeshType, FiniteElementType>(mesh);
    break;
  }
  case ScalarDiscretizationTypeBase::lagrangianFEM1: {
    typedef
      typename FiniteElementTraits<CellType,
      ScalarDiscretizationTypeBase::lagrangianFEM1>::Type
      FiniteElementType;
    __buildFEM<MeshType, FiniteElementType>(mesh);
    break;
  }
  case ScalarDiscretizationTypeBase::lagrangianFEM2: {
    typedef
      typename FiniteElementTraits<CellType,
      ScalarDiscretizationTypeBase::lagrangianFEM2>::Type
      FiniteElementType;
    __buildFEM<MeshType, FiniteElementType>(mesh);
    break;
  }
  case ScalarDiscretizationTypeBase::spectralLegendre: {
    switch (discretizationType.type()) {
    case ScalarDiscretizationTypeBase::spectralLegendre: {
      const ScalarDiscretizationTypeSpectral& spectralDiscretizationType
        = dynamic_cast<const ScalarDiscretizationTypeSpectral&>(discretizationType);
      __buildLegendre<MeshType>(mesh, spectralDiscretizationType);
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "Unknown spectral discretization type: '"+
                         ScalarDiscretizationTypeBase::name(discretizationType)+"'",
                         ErrorHandler::unexpected);
    }
    }
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "not implemented",
                       ErrorHandler::unexpected);
  }
  }
}


void ScalarDegreeOfFreedomPositionsSet::Builder::
build(const ScalarDiscretizationTypeBase& discretizationType,
      const Mesh& mesh)
{
  ffout(2) << "Computing degrees of freedom positions...\n";
  switch(mesh.type()) {
  case Mesh::cartesianHexahedraMesh: {
    __build(discretizationType, static_cast<const Structured3DMesh&>(mesh));
    break;
  }
  case Mesh::octreeMesh: {
    __build(discretizationType, static_cast<const OctreeMesh&>(mesh));
    break;
  }
  case Mesh::hexahedraMesh: {
    __build(discretizationType, static_cast<const MeshOfHexahedra&>(mesh));
    break;
  }
  case Mesh::spectralMesh: {
    __build(discretizationType, static_cast<const SpectralMesh&>(mesh));
    break;
  }
  case Mesh::tetrahedraMesh: {
    __build(discretizationType, static_cast<const MeshOfTetrahedra&>(mesh));
    break;
  }
  case Mesh::trianglesMesh: {
    __build(discretizationType, static_cast<const MeshOfTriangles&>(mesh));
    break;
  }
  case Mesh::surfaceMeshTriangles:
  case Mesh::surfaceMeshQuadrangles: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "Cannot solve problems on surface meshes",
                       ErrorHandler::normal);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "type of mesh not supported",
                       ErrorHandler::unexpected);
  }
  }
  ffout(2) << "Computing degrees of freedom positions: done\n";
}

void ScalarDegreeOfFreedomPositionsSet::
__build(const ScalarDiscretizationTypeBase& discretizationType,
        const Mesh& mesh)
{
  Builder builder(*this);

  builder.build(discretizationType, mesh);
}

---