MeshOfTriangles.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: MeshOfTriangles.cpp,v 1.1 2006/09/20 16:45:05 delpinux Exp $

#include <set>
#include <stack>

#include <Types.hpp>
#include <TinyVector.hpp>
#include <TinyMatrix.hpp>

#include <ConnectivityBuilder.hpp>
#include <ConformTransformation.hpp>

#include <P1Triangle3DFiniteElement.hpp>

#include <MeshOfTriangles.hpp>
#include <EdgesBuilder.hpp>

MeshOfTriangles::
MeshOfTriangles(ReferenceCounting<VerticesSet> vertices,
                ReferenceCounting<VerticesCorrespondance> correspondances,
                ReferenceCounting<Vector<Triangle> > triangles)
    : Mesh(Mesh::trianglesMesh,
           Mesh::plane,
           vertices,
           correspondances),
      __cells(triangles),
      __octree(0),
      __connectivity(*this)
{
  ;
}

void MeshOfTriangles::buildLocalizationTools()
{
  ffout(3) << "Generating localization in triangles mesh...\n";
  if (this->numberOfVertices()==0) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "the mesh contains no vertices\n",
                       ErrorHandler::normal);
  }

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

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

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

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

    TinyVector<3,real_t> X;
    T.value(P1Triangle3DFiniteElement::massCenter(), X);

    TinyVector<2,real_t> x;
    x[0]=X[0];
    x[1]=X[1];

    __octree->add(i,x);
  }
  ffout(3) << "Generating localization in triangles mesh: done\n";
}

MeshOfTriangles::const_iterator
MeshOfTriangles::find(const double& x,
                      const double& y,
                      const double& z) const
{
  if (__octree == 0) { // Builds the octree "on-demand"
    const_cast<MeshOfTriangles*>(this)->buildLocalizationTools();
  }

  // starts the search
  TinyVector<2,real_t> X;
  X[0] = x;
  X[1] = y;
  Octree<size_t, 2>::iterator i = __octree->fuzzySearch(X);

  size_t guessedCellNumber = (*i).value();
  const_iterator t0(*this, guessedCellNumber);

  bool found=false;

  std::set<MeshOfTriangles::const_iterator> visited;
  std::set<MeshOfTriangles::const_iterator> toVisitSet;
  std::stack<MeshOfTriangles::const_iterator> toVisitStack;

  toVisitSet.insert(t0);
  toVisitStack.push(t0);
  const_iterator c (*this);

  do {
    // treating next cell
    c = toVisitStack.top();

    toVisitSet.erase(c);
    toVisitStack.pop();
    visited.insert(c);

    const Triangle& t = *c;

    TinyVector<3, real_t> lambda;

    t.getBarycentricCoordinates(X, lambda);

    found = true;
    for (size_t i=0; i<3; ++i) {
#warning Should use a relevent epsilon here ...
      if (lambda[i] < -1E-6) {
        found = false;
        const Triangle* newT = __connectivity.cells(t)[i];
        if (newT != 0) {
          const_iterator t1(*this, this->cellNumber(*newT));
  
          if (visited.find(t1) == visited.end()) {
            if (toVisitSet.find(t1) == toVisitSet.end()) {
              toVisitSet.insert(t1);
              toVisitStack.push(t1);
            }
          }
        }
      }
    }
  } while (not(found) and not(toVisitStack.empty()));

  if (not(found)) {
    c = const_iterator(*this,const_iterator::End);
  }
  return  c;
}

---