OctreeMeshBuilder.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

//  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: OctreeMeshBuilder.cpp,v 1.3 2007/08/21 10:32:28 delpinux Exp $

#include <OctreeMeshBuilder.hpp>
#include <Structured3DMesh.hpp>

#include <ErrorHandler.hpp>

void OctreeMeshBuilder::
buildMesh()
{
  if (__mesh->type() != Mesh::cartesianHexahedraMesh) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "cannot build octree mesh on "
                       +__mesh->typeName()
                       +": structured mesh is mandatory!",
                       ErrorHandler::normal);
  }

  const Structured3DMesh& mesh
    = dynamic_cast<const Structured3DMesh&>(*__mesh);
  const Domain& domain = *__domain;

  ffout(2) << "Starting Octree mesh builder\n";
  Vector<bool> domainVertex(mesh.numberOfVertices());

  for (size_t i=0; i<mesh.numberOfVertices(); ++i) {
    domainVertex[i] = domain.inside(mesh.vertex(i));
  }

  std::vector<size_t> inCells;
  std::vector<size_t> outCells;
  std::vector<size_t> boundaryCells;

  for(Structured3DMesh::const_iterator icell(mesh); not icell.end(); ++icell) {
    bool allIn  = true;
    bool allOut = true;
    const CartesianHexahedron& K = *icell;
    for (size_t i=0; i<CartesianHexahedron::NumberOfVertices; ++i) {
      size_t vertexNumber = mesh.vertexNumber(K(i));
      if (domainVertex[vertexNumber]) {
        allOut = false;
      } else {
        allIn = false;
      }
    }
    if (allIn) {
      inCells.push_back(icell.number());
    } else if (allOut) {
      outCells.push_back(icell.number());
    } else {
      boundaryCells.push_back(icell.number());
    }
  }

  typedef std::map<size_t, TinyVector<3,real_t> > VerticesList;
  VerticesList verticesList;

  typedef TinyVector<3, size_t> NodeIndex;
  typedef TinyVector<3, real_t> NodeCoord;
  typedef TinyVector<8, NodeIndex> CellDescriptor;
  typedef std::map<NodeIndex, NodeCoord> NodeList;
  typedef std::map<NodeIndex, size_t> NodeNumber;
  typedef std::vector<CellDescriptor> CellList;

  std::vector<NodeList> nodeListSet(__maximumLevel+1);
  std::vector<CellList> cellListSet(__maximumLevel+1); 

  { // Building level 0
    NodeList& nodeList = nodeListSet[0];
    CellList& cellList = cellListSet[0];

    for(size_t i=0; i<inCells.size(); ++i) {
      CellDescriptor cell;
      const CartesianHexahedron& H = mesh.cell(inCells[i]);
      for(size_t j=0; j<CartesianHexahedron::NumberOfVertices; ++j) {
        NodeIndex nodeIndex = mesh.vertexIndex(H(j))*(1<<__maximumLevel);

        cell[j] = nodeIndex;

        NodeCoord nodeCoord;
        nodeCoord[0] = nodeIndex[0]*mesh.shape().hx()/(1<<__maximumLevel)+mesh.shape().a()[0];
        nodeCoord[1] = nodeIndex[1]*mesh.shape().hy()/(1<<__maximumLevel)+mesh.shape().a()[1];
        nodeCoord[2] = nodeIndex[2]*mesh.shape().hz()/(1<<__maximumLevel)+mesh.shape().a()[2];

        nodeList[nodeIndex] = nodeCoord;
      }
      cellList.push_back(cell);
    }
  }

  // list of the cells to examine at the current level
  CellList cellToCut;
  cellToCut.reserve(boundaryCells.size());
  // and their vertices
  NodeList cellToCutNode;

  { // Convert data structures to more flexible one
    for(size_t i=0; i<boundaryCells.size(); ++i) {
      const CartesianHexahedron& H = mesh.cell(boundaryCells[i]);
      CellDescriptor cell;
      Index index0 = mesh.vertexIndex(H(0));
      for(size_t j=0; j<CartesianHexahedron::NumberOfVertices; ++j) {
        Index index = mesh.vertexIndex(H(j));
        NodeIndex nodeIndex = mesh.vertexIndex(H(j))*(1<<__maximumLevel);

        NodeCoord nodeCoord;
        nodeCoord[0] = nodeIndex[0]*mesh.shape().hx()/(1<<__maximumLevel)+mesh.shape().a()[0];
        nodeCoord[1] = nodeIndex[1]*mesh.shape().hy()/(1<<__maximumLevel)+mesh.shape().a()[1];
        nodeCoord[2] = nodeIndex[2]*mesh.shape().hz()/(1<<__maximumLevel)+mesh.shape().a()[2];

        cell[j] = nodeIndex;

        cellToCutNode[nodeIndex] = nodeCoord;
        cellToCut.push_back(nodeIndex);
      }
    }
  }

  TinyVector<8,TinyVector<3,size_t> > shift;
  shift[0] = TinyVector<3,size_t>(0,0,0);
  shift[1] = TinyVector<3,size_t>(1,0,0);
  shift[2] = TinyVector<3,size_t>(1,1,0);
  shift[3] = TinyVector<3,size_t>(0,1,0);
  shift[4] = TinyVector<3,size_t>(0,0,1);
  shift[5] = TinyVector<3,size_t>(1,0,1);
  shift[6] = TinyVector<3,size_t>(1,1,1);
  shift[7] = TinyVector<3,size_t>(0,1,1);


  const TinyVector<3,real_t> meshStep(mesh.shape().hx()/(1<<(__maximumLevel)),
                                      mesh.shape().hy()/(1<<(__maximumLevel)),
                                      mesh.shape().hz()/(1<<(__maximumLevel)));

  for (size_t level=1; level<=__maximumLevel; ++level) {
    NodeList& nodeList = nodeListSet[level];
    CellList& cellList = cellListSet[level];

    ffout(3) << "- Creating level " << level << '\n';
    CellList cuttedCell;
    NodeList cuttedCellNode;
    for (CellList::const_iterator icell = cellToCut.begin();
         icell != cellToCut.end(); ++icell) {
      const CellDescriptor& oldCell = *icell;
      const NodeIndex& baseNode = oldCell[0];
      for (size_t i=0; i<8; ++i) {
        CellDescriptor cell;
        for (size_t j=0; j<8; ++j) {
          const NodeIndex newNode = baseNode+(shift[i]+shift[j])*(1<<(__maximumLevel-level));
          cell[j] = newNode;
          NodeCoord nodeCoord;
          nodeCoord[0] = newNode[0]*meshStep[0] + mesh.shape().a()[0];
          nodeCoord[1] = newNode[1]*meshStep[1] + mesh.shape().a()[1];
          nodeCoord[2] = newNode[2]*meshStep[2] + mesh.shape().a()[2];

          cuttedCellNode[newNode] = nodeCoord;
        }
        cuttedCell.push_back(cell);
      }
    }

    // We don't need this list anymore
    cellToCut.clear();

    // we now check which vertices are inside the domain
    typedef std::map<NodeIndex, bool> NodeInside;
    NodeInside nodeInside;
    for (NodeList::const_iterator inode = cuttedCellNode.begin();
         inode != cuttedCellNode.end(); ++inode) {
      nodeInside[inode->first] = domain.inside(inode->second);
    }

    for(CellList::const_iterator icell = cuttedCell.begin();
        icell != cuttedCell.end(); ++icell) {
      bool allIn  = true;
      bool allOut = true;
      const CellDescriptor& cell = *icell;

      for (size_t i=0; i<CartesianHexahedron::NumberOfVertices; ++i) {
        if (nodeInside[cell[i]]) {
          allOut = false;
        } else {
          allIn = false;
        }
      }
      if (allIn) {
        cellList.push_back(cell);
        for (size_t i=0; i<CartesianHexahedron::NumberOfVertices; ++i) {
          nodeList[cell[i]] = cuttedCellNode[cell[i]];
        }
      } else if ((not allOut) and (level<__maximumLevel)) {
        cellToCut.push_back(cell);
      }
    }
  }

  std::vector<NodeNumber> nodeNumberSet(__maximumLevel+1);

  size_t numberOfNodes=0;
  { // computes new node numbers
    for(size_t level = 0; level<__maximumLevel+1; ++level) {
      NodeNumber& nodeNumber = nodeNumberSet[level];
      NodeList& nodeList = nodeListSet[level];
      for (NodeList::const_iterator i=nodeList.begin();
           i != nodeList.end(); ++i) {
        nodeNumber[i->first] = numberOfNodes;
        numberOfNodes++;
      }
    }
  }

  ReferenceCounting<VerticesSet> pVerticesSet;
  ReferenceCounting<VerticesCorrespondance> pVerticesCorrespondance;

  {// Keeping required vertices
    pVerticesSet = new VerticesSet(numberOfNodes);
    VerticesSet& verticesSet = *pVerticesSet;
    size_t k=0;
    for(size_t level = 0; level<__maximumLevel+1; ++level) {
      NodeList& nodeList = nodeListSet[level];
      for (NodeList::const_iterator i=nodeList.begin();
           i != nodeList.end(); ++i) {
        verticesSet[k] = i->second;
        k++;
      }
    }

    pVerticesCorrespondance
      = new VerticesCorrespondance(verticesList.size());
  }

  VerticesSet& verticesSet = *pVerticesSet;

  size_t numberOfCells = 0;
  for(size_t level = 0; level<__maximumLevel+1; ++level) {
    const CellList& cellList = cellListSet[level];
    numberOfCells += cellList.size();
  }

  ReferenceCounting<Vector<CartesianHexahedron> > pCellsSet
    = new Vector<CartesianHexahedron>(numberOfCells);
  Vector<CartesianHexahedron>& cellsSet = *pCellsSet;
  {// Adding cells to the list
    size_t k=0;
    for(size_t level = 0; level<__maximumLevel+1; ++level) {
      CellList& cellList = cellListSet[level];
      NodeNumber& nodeNumber = nodeNumberSet[level];
      for (size_t i=0; i<cellList.size();++i) {
        const CellDescriptor& cell = cellList[i];

        cellsSet[k] = CartesianHexahedron(verticesSet[nodeNumber[cell[0]]],
                                          verticesSet[nodeNumber[cell[1]]],
                                          verticesSet[nodeNumber[cell[2]]],
                                          verticesSet[nodeNumber[cell[3]]],
                                          verticesSet[nodeNumber[cell[4]]],
                                          verticesSet[nodeNumber[cell[5]]],
                                          verticesSet[nodeNumber[cell[6]]],
                                          verticesSet[nodeNumber[cell[7]]]);
        ++k;
      }
    }
  }

  for (size_t level=0; level<__maximumLevel+1; ++level) {
    ffout(4) << "  - Level " << level << '\n';
    ffout(4) << "    - Number of cells: " << cellListSet[level].size() << '\n';
    ffout(4) << "    - Number of nodes: " << nodeNumberSet[level].size() << '\n';
  }

  ffout(3) << "- Number of cells: " << cellsSet.size() << '\n';
  ffout(3) << "- Number of nodes: " << verticesSet.numberOfVertices() << '\n';

  __octreeMesh = new OctreeMesh(pVerticesSet,
                                pVerticesCorrespondance,
                                pCellsSet);

  ffout(2) << "Octree mesh builder finished\n";
}

---