MeshOfHexahedra.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: MeshOfHexahedra.cpp,v 1.21 2008/03/04 00:17:58 delpinux Exp $

#include <MeshOfHexahedra.hpp>
#include <ConnectivityBuilder.hpp>
#include <set>
#include <queue>

#include <ConformTransformation.hpp>
#include <FacesBuilder.hpp>

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

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


void MeshOfHexahedra::
buildLocalizationTools()
{
  if (this->numberOfVertices()==0) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "the mesh contains no vertices\n",
                       ErrorHandler::normal);
  }

  {
    // Localization inside a hexahedra mesh requires the cell to cell
    // connectivity
    ConnectivityBuilder<MeshOfHexahedra> builder(*this);
    builder.generates(Connectivity<MeshOfHexahedra>::CellToCells);
  }

  if (__borderMesh != 0) {
    ConnectivityBuilder<MeshOfHexahedra> builder(*this);
    builder.borderMesh(__borderMesh);
    __borderMesh->setBackgroundMesh(this);
  }

  // getting bounding box size
  __a = (this->vertex(0));
  __b = __a;
  for (size_t i=1; i<this->numberOfVertices(); ++i) {
    Vertex& x = this->vertex(i);
    for (size_t k=0; k<3; ++k) {
      __a[k] = std::min(__a[k],x[k]);
      __b[k] = std::max(__b[k],x[k]);
    }
  }
  ffout(3) << "- Bounding box is " << __a << ',' << __b << '\n';

  __a -= TinyVector<3,real_t>(1,1,1);
  __b += TinyVector<3,real_t>(1,1,1);

  ffout(3) << "- Building the octree\n";
  __octree = new Octree<size_t, 3>(__a,__b);

  for (size_t i = 0; i < this->numberOfCells(); ++i) {

    ConformTransformationQ1Hexahedron T(this->cell(i));

    for (size_t k=0; k<QuadratureFormulaQ1Hexahedron::numberOfQuadraturePoints; ++k) {
      TinyVector<3,real_t> X;
      T.value(QuadratureFormulaQ1Hexahedron::instance()[k], X);
      __octree->add(i,X);
    }
  }
}


MeshOfHexahedra::const_iterator
MeshOfHexahedra::find(const double& x,
                      const double& y,
                      const double& z) const
{
  TinyVector<3,real_t> X(x,y,z);

  Octree<size_t, 3>::iterator i = (*__octree).fuzzySearch(X);

  size_t cellNumber = (*i).value();

  MeshOfHexahedra::const_iterator h0(*this, cellNumber);

  bool found=false;

  std::set<MeshOfHexahedra::const_iterator> visited;
  std::set<MeshOfHexahedra::const_iterator> toVisitSet;
  std::queue<MeshOfHexahedra::const_iterator> toVisitQueue;

  toVisitSet.insert(h0);
  toVisitQueue.push(h0);
  MeshOfHexahedra::const_iterator h(*this);

  TinyVector<3, real_t> Xhat;
  size_t cpt = 0;
  do {
    if (cpt > 125) {
      break; // test 3 layers arround the seeked point then consideres that point was missed
    }
    // treating next cell
    h = toVisitQueue.front();

    toVisitSet.erase(h);
    toVisitQueue.pop();
    visited.insert(h);

    const Hexahedron& H = *h;

    ConformTransformationQ1Hexahedron T(H);
    found = T.invertT(x,y,z, Xhat);

    if (not found) {
      for (size_t k=0; k<6; ++k) {
        MeshOfHexahedra::const_iterator icell(*this,const_iterator::End);
        icell = __connectivity.cells(*h)[k];
        if (not(icell.end())) {
          if (visited.find(icell) == visited.end()) {
            if (toVisitSet.find(icell) == toVisitSet.end()) {
              toVisitSet.insert(icell);
              toVisitQueue.push(icell);
            }
          }
        }
      }
    }
    cpt++;
  } while (not(found) and not(toVisitQueue.empty()));

  if (not found) {
    return MeshOfHexahedra::const_iterator(*this,const_iterator::End);
  }
  return  h;
}

---