DegreeOfFreedomSetBuilder.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: DegreeOfFreedomSetBuilder.cpp,v 1.17 2007/08/21 16:23:26 delpinux Exp $

#include <DegreeOfFreedomSetBuilder.hpp>
#include <Mesh.hpp>

#include <DegreeOfFreedomSetManager.hpp>

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

#include <OctreeMesh.hpp>

#include <FiniteElementTraits.hpp>

#include <ErrorHandler.hpp>

#include <ScalarDiscretizationTypeSpectral.hpp>


template <size_t ItemDOFNumber>
struct DegreeOfFreedomSetBuilder::ItemTraits
{
  typedef TinyVector<ItemDOFNumber,size_t> DOFList;
};

template <>
struct DegreeOfFreedomSetBuilder::ItemTraits<0>
{
  // This Traits is built to avoid the TinyVector<0,size_t> which is undefined
  // This typedef only acts for correct syntaxe
  typedef TinyVector<1,size_t> DOFList;
};

template <typename MeshType,
          typename FiniteElementType>
size_t DegreeOfFreedomSetBuilder::
__buildFEMCorrespondance(const MeshType& mesh,
                         DegreeOfFreedomSet::Correspondance& correspondance)
{
  size_t numberOfUsedDOF = 0;
  size_t dofPositionNumber = 0;

  const VerticesCorrespondance& verticesCorrespondance = *mesh.verticesCorrespondance();

  // vertices dof management
  if (FiniteElementType::numberOfVertexDegreesOfFreedom > 0) {
    for (size_t i=0; i<mesh.numberOfVertices(); ++i) {
      correspondance[i] = verticesCorrespondance[i];
      numberOfUsedDOF = std::max(verticesCorrespondance[i],numberOfUsedDOF);
      dofPositionNumber++;
    }
  }

  // edge dof management
  if (FiniteElementType::numberOfEdgeDegreesOfFreedom > 0) {
    ASSERT(mesh.hasEdges());

    typedef TinyVector<2,size_t> VerticesPair;

    typedef
      typename ItemTraits<FiniteElementType::numberOfEdgeDegreesOfFreedom>::DOFList
      EdgeDOFNumbers;

    typedef std::map<VerticesPair, EdgeDOFNumbers> SortedEdgeSet;
    SortedEdgeSet sortedEdgeSet;

    for (size_t i=0; i<mesh.numberOfEdges(); ++i) {
      const Edge& edge = mesh.edge(i);
      const size_t v0 = verticesCorrespondance[mesh.vertexNumber(edge(0))];
      const size_t v1 = verticesCorrespondance[mesh.vertexNumber(edge(1))];

      TinyVector<2,size_t> verticesPair;
      verticesPair[0] = std::min(v0,v1);
      verticesPair[1] = std::max(v0,v1);

      typename SortedEdgeSet::iterator iedge = sortedEdgeSet.find(verticesPair);
      if (iedge == sortedEdgeSet.end()) {
        EdgeDOFNumbers dofNumbers;
        for (size_t j=0; j<FiniteElementType::numberOfEdgeDegreesOfFreedom; ++j) {
          numberOfUsedDOF++;
          dofNumbers[j]=numberOfUsedDOF;
        }
        
        iedge = sortedEdgeSet.insert(iedge,std::make_pair(verticesPair, dofNumbers));
      }
      for (size_t j=0; j<FiniteElementType::numberOfEdgeDegreesOfFreedom; ++j) {
        correspondance[dofPositionNumber] = iedge->second[j];
        dofPositionNumber++;
      }
    }
  }

  // faces dof management
  if (FiniteElementType::numberOfFaceDegreesOfFreedom > 0) {
    ASSERT(mesh.hasFaces());

    typedef TinyVector<MeshType::FaceType::NumberOfVertices,size_t> VerticesList;

    typedef
      typename ItemTraits<FiniteElementType::numberOfFaceDegreesOfFreedom>::DOFList
      FaceDOFNumbers;

    typedef std::map<VerticesList, FaceDOFNumbers> SortedFaceSet;
    SortedFaceSet sortedFaceSet;

    for (size_t i=0; i<mesh.numberOfFaces(); ++i) {
      const typename MeshType::FaceType& face = mesh.face(i);
      std::set<size_t> verticesSet;

      for (size_t j=0; j<MeshType::FaceType::NumberOfVertices; ++j) {
        verticesSet.insert(verticesCorrespondance[mesh.vertexNumber(face(j))]);
      }
      VerticesList verticesList;
      {
        size_t j=0;
        for (std::set<size_t>::const_iterator k = verticesSet.begin();
             k != verticesSet.end(); ++k, ++j) {
          verticesList[j] = *k;
        }
      }

      typename SortedFaceSet::iterator iface = sortedFaceSet.find(verticesList);
      if (iface == sortedFaceSet.end()) {
        FaceDOFNumbers dofNumbers;
        for (size_t j=0; j<FiniteElementType::numberOfFaceDegreesOfFreedom; ++j) {
          numberOfUsedDOF++;
          dofNumbers[j]=numberOfUsedDOF;
        }
        
        iface = sortedFaceSet.insert(iface,std::make_pair(verticesList, dofNumbers));
      }
      for (size_t j=0; j<FiniteElementType::numberOfFaceDegreesOfFreedom; ++j) {
        correspondance[dofPositionNumber] = iface->second[j];
        dofPositionNumber++;
      }
    }
  }

  // cell dof management
  if (FiniteElementType::numberOfVolumeDegreesOfFreedom > 0) {
    for (size_t i=0;i<mesh.numberOfCells(); ++i) {
      for (size_t j=0;j<FiniteElementType::numberOfVolumeDegreesOfFreedom; ++j) {
        correspondance[dofPositionNumber] = numberOfUsedDOF;
        dofPositionNumber++;
        numberOfUsedDOF++;
      }
    }
  }
  return numberOfUsedDOF+1;
}

template <typename MeshType>
void DegreeOfFreedomSetBuilder::
__buildCorrespondance(const MeshType& mesh,
                      const DiscretizationType& discretization,
                      DegreeOfFreedomSet::Correspondance& correspondance)
{
  if (not mesh.isPeriodic()) {
    for (size_t i=0; i<correspondance.size(); ++i) {
      correspondance[i]=i;
    }
    return;
  }

  typedef typename MeshType::CellType CellType;

  switch (discretization[0].type()) {
  case ScalarDiscretizationTypeBase::lagrangianFEM1: {
    typedef
      typename FiniteElementTraits<CellType,
      ScalarDiscretizationTypeBase::lagrangianFEM1>::Type
      FiniteElementType;

    this->__buildFEMCorrespondance<MeshType,FiniteElementType>(mesh,correspondance);
    break;
  }
  case ScalarDiscretizationTypeBase::lagrangianFEM2: {
    typedef
      typename FiniteElementTraits<CellType,
      ScalarDiscretizationTypeBase::lagrangianFEM2>::Type
      FiniteElementType;

    this->__buildFEMCorrespondance<MeshType,FiniteElementType>(mesh,correspondance);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "not implemented",
                       ErrorHandler::unexpected);
  }
  }

  {
#warning this only works for Pk-Pk elements
    size_t offset = 0;
    for (size_t i=1; i<__dofPositionSet->number(); ++i) {
      offset += (*__dofPositionSet)[i-1].number();
      for (size_t j=0; j<(*__dofPositionSet)[i].number(); ++j) {
        correspondance[j+offset] = correspondance[j]+offset;
      }
    }
  }
}

template <typename MeshType,
          typename FiniteElementType>
void DegreeOfFreedomSetBuilder::
__buildFEMFictitious(const size_t& numberOfVariables,
                     const MeshType& mesh,
                     const Domain& domain)
{
  typedef DegreeOfFreedomSet::Correspondance Correspondance;

  size_t numberOfCorrespondances = 0;
  for (size_t i=0; i<__dofPositionSet->number(); ++i) {
    numberOfCorrespondances += (*__dofPositionSet)[i].number();
  }

  ReferenceCounting<Correspondance> pCorrespondance
    = new Correspondance(numberOfCorrespondances);

  Correspondance& correspondance = *pCorrespondance;

  // ONLY VALID FOR Pk-Pk elements
  this->__buildFEMCorrespondance<MeshType,FiniteElementType>(mesh,correspondance);

  {
#warning this only works for Pk-Pk elements
    size_t offset = 0;
    for (size_t i=1; i<__dofPositionSet->number(); ++i) {
      offset += (*__dofPositionSet)[i-1].number();
      for (size_t j=0; j<(*__dofPositionSet)[i].number(); ++j) {
        correspondance[j+offset] = correspondance[j]+offset;
      }
    }
  }

  Vector<bool> nonFictitiousDOF(numberOfCorrespondances);
  nonFictitiousDOF = false;

  {
    size_t offset = 0;
    for(size_t i=0; i<__dofPositionSet->number(); ++i) {
      const ScalarDegreeOfFreedomPositionsSet& scalarDOFPositions = (*__dofPositionSet)[i];
      if (i>0) offset += (*__dofPositionSet)[i-1].number();
      for (size_t j=0; j<scalarDOFPositions.number(); ++j) {
        nonFictitiousDOF[offset+j] = domain.inside(scalarDOFPositions.vertex(j));
      }
    }
  }

  typedef typename MeshType::CellType CellType;

  size_t offset = 0;
  for(size_t i=0; i<(*__dofPositionSet).number(); ++i) {
    if (i>0) offset += (*__dofPositionSet)[i-1].number();
    const ScalarDegreeOfFreedomPositionsSet& scalarDOFPositions = (*__dofPositionSet)[i];
    for (size_t j=0; j<mesh.numberOfCells();++j) {
      size_t numberIn = 0;
      const CellType& C = mesh.cell(j);
      const size_t cellNumer = mesh.cellNumber(C);
      // Computation of the caracteristic function
      for (size_t k=0; k<FiniteElementType::numberOfDegreesOfFreedom; ++k) {
        numberIn += (nonFictitiousDOF(offset+scalarDOFPositions(cellNumer,k)))?1:0;
      }

      // reduction of degrees of freedom
      if (numberIn != 0) {
        for (size_t k=0; k<FiniteElementType::numberOfDegreesOfFreedom; ++k) {
          const size_t n = scalarDOFPositions(cellNumer,k);
          nonFictitiousDOF[offset+n] = true;
        }
      }
    }
  }

  // numbering correspondant vertices remapping
  std::map<size_t, size_t> newVerticesCorrespondance;
  for (size_t i=0; i<numberOfCorrespondances; ++i) {
    if (nonFictitiousDOF[i]) {
      newVerticesCorrespondance[correspondance[i]] = 0;
    }
  }

  size_t numberOfDLVertices = 0;
  for (std::map<size_t, size_t>::iterator i=newVerticesCorrespondance.begin();
       i != newVerticesCorrespondance.end(); ++i) {
    i->second = numberOfDLVertices++;
  }

  ReferenceCounting<Correspondance> pFDMCorrespondance
    = new Correspondance(numberOfCorrespondances);

  DegreeOfFreedomSet::Correspondance& fdmCorrespondance = *pFDMCorrespondance;
  fdmCorrespondance = -1;

  for(size_t i=0; i<correspondance.size(); ++i) {
    if (nonFictitiousDOF[i]) {
      fdmCorrespondance[i] = newVerticesCorrespondance[correspondance[i]];
    }
  }
  
  __degreeOfFreedomSet
    = new DegreeOfFreedomSet(__dofPositionSet,
                             pFDMCorrespondance);
}

template <typename MeshType>
void DegreeOfFreedomSetBuilder::
__buildFictitious(const size_t& numberOfVariables,
                  const ScalarDiscretizationTypeBase& discretization,
                  const MeshType& mesh,
                  const Domain& domain)
{
  typedef typename MeshType::CellType CellType;

  switch (discretization.type()) {
  case ScalarDiscretizationTypeBase::lagrangianFEM1: {
    typedef
      typename FiniteElementTraits<CellType,
      ScalarDiscretizationTypeBase::lagrangianFEM1>::Type
      FiniteElementType;

    this->__buildFEMFictitious<MeshType,FiniteElementType>(numberOfVariables,
                                                           mesh,
                                                           domain);
    break;
  }
  case ScalarDiscretizationTypeBase::lagrangianFEM2: {
    typedef
      typename FiniteElementTraits<CellType,
      ScalarDiscretizationTypeBase::lagrangianFEM2>::Type
      FiniteElementType;

    this->__buildFEMFictitious<MeshType,FiniteElementType>(numberOfVariables,
                                                           mesh,
                                                           domain);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "not implemented",
                       ErrorHandler::unexpected);
  }
  }
}

DegreeOfFreedomSetBuilder::
DegreeOfFreedomSetBuilder(const DiscretizationType& discretizationType,
                          const Mesh& mesh)
  : __dofPositionSet(new DegreeOfFreedomPositionsSet)
{
  for (size_t i=0; i<discretizationType.number(); ++i) {
    __dofPositionSet->add(DegreeOfFreedomSetManager::instance().getDOFPositionsSet(mesh,
                                                                                   discretizationType[i]));
  }

  typedef DegreeOfFreedomSet::Correspondance Correspondance;

  size_t numberOfCorrespondances = 0;
  for (size_t i=0; i<__dofPositionSet->number(); ++i) {
    numberOfCorrespondances += (*__dofPositionSet)[i].number();
  }

  ReferenceCounting<Correspondance> pCorrespondance
    = new Correspondance(numberOfCorrespondances);
  Correspondance& correspondance = *pCorrespondance;

  switch(mesh.type()) {
  case Mesh::cartesianHexahedraMesh: {
    this->__buildCorrespondance(static_cast<const Structured3DMesh&>(mesh),
                                discretizationType,
                                correspondance);
    break;
  }
  case Mesh::hexahedraMesh: {
    this->__buildCorrespondance(static_cast<const MeshOfHexahedra&>(mesh),
                                discretizationType,
                                correspondance);
    break;
  }
  case Mesh::tetrahedraMesh: {
    this->__buildCorrespondance(static_cast<const MeshOfTetrahedra&>(mesh),
                                discretizationType,
                                correspondance);
    break;
  }
  case Mesh::spectralMesh: {
    this->__buildCorrespondance(static_cast<const SpectralMesh&>(mesh),
                                discretizationType,
                                correspondance);
    break;
  }
  case Mesh::octreeMesh: {
    this->__buildCorrespondance(static_cast<const OctreeMesh&>(mesh),
                                discretizationType,
                                correspondance);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "not implemented yet",
                       ErrorHandler::unexpected);
  }
  }

  __degreeOfFreedomSet
    = new DegreeOfFreedomSet(__dofPositionSet,
                             pCorrespondance);

}

DegreeOfFreedomSetBuilder::
DegreeOfFreedomSetBuilder(const DiscretizationType& discretizationType,
                          const Mesh& mesh,
                          const Domain& domain)
  : __dofPositionSet(new DegreeOfFreedomPositionsSet)
{
  for (size_t i=0; i<discretizationType.number(); ++i) {
    __dofPositionSet->add(DegreeOfFreedomSetManager::instance().getDOFPositionsSet(mesh,
                                                                                   discretizationType[i]));
  }

  switch(mesh.type()) {
  case Mesh::cartesianHexahedraMesh: {
    this->__buildFictitious(discretizationType.number(),
                            discretizationType[0],
                            static_cast<const Structured3DMesh&>(mesh),
                            domain);
    break;
  }
  case Mesh::tetrahedraMesh: {
    this->__buildFictitious(discretizationType.number(),
                            discretizationType[0],
                            static_cast<const Structured3DMesh&>(mesh),
                            domain);
    break;
  }
  case Mesh::hexahedraMesh: {
    this->__buildFictitious(discretizationType.number(),
                            discretizationType[0],
                            static_cast<const Structured3DMesh&>(mesh),
                            domain);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "not implemented yet",
                       ErrorHandler::unexpected);
  }
  }
}


DegreeOfFreedomSetBuilder::
~DegreeOfFreedomSetBuilder()
{
  ;
}

---