STLreaderForSubgridParticleWallContact.h

Go to the documentation of this file.

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2010-2024 Thomas Henn, Mathias J. Krause, Christoph Gaul
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  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
 *  of the License, 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., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/
 
#ifndef STL_READER_H
#define STL_READER_H
 
#include <string>
#include <vector>
#include <iostream>
#include <sstream>
#include <set>
#include <limits>
 
#include "communication/loadBalancer.h"
#include "geometry/cuboidDecomposition.h"
#include "functors/analytical/indicator/indicatorF3D.h"
#include "functors/analytical/indicator/indicatorBaseF3D.h"
#include "utilities/vectorHelpers.h"
#include "octree.h"
#include "core/vector.h"
 
 
// All OpenLB code is contained in this namespace.
namespace olb {
 
template<typename T>
class Octree;
 
template<typename T, unsigned D>
struct Vertex {
  Vertex() : coords() {};
  Vertex( Vector<T,D> coords_ ) : coords(coords_) {};
  Vertex<T,D>& operator=(Vertex<T,D> const& rhs)
  {
    = rhs;
    return *this;
  };
  Vertex(Vertex<T,D> const& rhs)(rhs) {};
 
  Vector<T,D> coords;
 
  int getDim() const
  {
    return D;
  }
 
};
 
template<typename T>
struct STLtriangle {
  bool testRayIntersect(const Vector<T,3>& pt,const Vector<T,3>& dir, Vector<T,3>& q, T& alpha, const T& rad = T(), bool print = false);
  Vector<T,3> closestPtPointTriangle(const Vector<T,3>& pt) const;
 
  std::vector<Vertex<T,3> > point;
 
  Vector<T,3> normal;
 
  Vector<T,3> uBeta, uGamma;
  T d, kBeta, kGamma;
 
public:
  STLtriangle():point(3, Vertex<T,3>()), normal(T()), uBeta(T()), uGamma(T()), d(T()), kBeta(T()), kGamma(T()) {};
  STLtriangle(STLtriangle<T> const& tri):point(tri.point), normal(tri.normal), uBeta(tri.uBeta), uGamma(tri.uGamma), d(tri.d), kBeta(tri.kBeta), kGamma(tri.kGamma) {};
  STLtriangle<T>& operator=(STLtriangle<T> const& tri)
  {
    point = tri.point;
    normal = tri.normal;
    uBeta = tri.uBeta;
    uGamma = tri.uGamma;
    d = tri.d;
    kBeta = tri.kBeta;
    kGamma = tri.kGamma;
    return *this;
  };
 
  ~STLtriangle() {};
 
  void init();
  inline Vector<T,3>& getNormal()
  {
    return normal;
  }
  inline const Vector<T,3>& getNormal() const
  {
    return normal;
  }
  Vector<T,3> getCenter();
  std::vector<T> getE0();
  std::vector<T> getE1();
  bool isPointInside(const PhysR<T,3>& pt) const;
  bool getPointToEdgeDistances (const Vector<T,3>& input, Vector<T,3> & output, T sensitivity = 1.e-15);
  bool isEdgePoint (const Vector<T,3> & input,  Vector<T,3>& P1, Vector<T,3> & P2, T sensitivity = 1.e-15);
  bool isVortexPoint (const Vector<T,3> & input,Vector<T,3>& P, T sensitivity = 1.e-15);
};
 
template<typename T>
class STLmesh {
  T distPoints(Vertex<T,3>& p1, Vertex<T,3>& p2);
  const std::string _fName;
  std::vector<STLtriangle<T> > _triangles;
  Vector<T,3> _min, _max;
  T _maxDist2;
  mutable OstreamManager clout;
 
public:
  STLmesh(std::string, T stlSize = 1.);
 
  STLmesh(const std::vector<std::vector<T>> meshPoints, T stlSize = 1.);
 
  inline STLtriangle<T>& getTri(unsigned int i)
  {
    return _triangles[i];
  }
  inline std::vector<STLtriangle<T> >& getTriangles()
  {
    return _triangles;
  }
  inline unsigned int triangleSize() const
  {
    return _triangles.size();
  }
  inline Vector<T,3>& getMin()
  {
    return _min;
  };
  inline Vector<T,3>& getMax()
  {
    return _max;
  };
  inline float maxDist2() const
  {
    return _maxDist2;
  }
  void print(bool full = false);
  void write(std::string fName);
  bool testRayIntersect(const std::set<unsigned int>& tris, const Vector<T,3>& pt,const Vector<T,3>& dir, Vector<T,3>& q, T& alpha);
};
 
enum class SignMode { EXACT, CACHED };
 
 
enum class RayMode : int {
  FastRayZ  = 1,  
  Robust    = 2,  
  FastRayX  = 3,  
  FastRayY  = 4,  
  DoubleRay = 5   
};
 
 
template<typename T>
class STLreader : public IndicatorF3D<T> {
private:
  /*
   *  Old indicate function (slower, more stable)
   *  Define three rays (X-, Y-, Z-direction) for each leaf and count intersections
   *  with STL for each ray. Odd number of intersection means inside (Majority vote).
   */
  void indicate1();
  /*
   *  New indicate function (faster, less stable)
   *  Define ray in Z-direction for each Voxel in XY-layer. Indicate all nodes on the fly.
   */
  void indicate2();
  /*
  *  New indicate function (faster, less stable)
  *  Define ray in X-direction for each Voxel in YZ-layer. Indicate all nodes on the fly.
  */
  void indicate2_Xray();
  /*
  *  New indicate function (faster, less stable)
  *  Define ray in Y-direction for each Voxel in XZ-layer. Indicate all nodes on the fly.
  */
  void indicate2_Yray();
  /*
   *  Double ray approach: two times (X-, Y-, Z-direction) for each leaf.
   *  Could be use to deal with double layer triangles and face intersections.
   */
 
  void indicate3();
 
  template <typename F>
  void iterateOverCloseTriangles(const PhysR<T,3>& pt, F func, Octree<T>* leafNode = nullptr);
 
  Vector<T,3> evalNormalOnSurface(const PhysR<T,3>& pt, const Vector<T,3>& fallbackNormal);
 
  template <SignMode SIGNMODE>
  Vector<T,3> evalSurfaceNormal(const Vector<T,3>& origin);
 
  template <SignMode SIGNMODE>
  short evalSignForSignedDistance(const Vector<T,3>& pt, [[maybe_unused]] const T distance, Vector<T,3>  vecdist = {} , STLtriangle<T> stlT = STLtriangle<T>() );
 
template<SignMode SIGNMODE>
Vector<T,3> evalSurfaceNormalForPseudoNormal(const Vector<T,3>& origin, Vector<T,3> & outputPointOnSurface);
 
 
  short evalSignForSignedDistanceFromNormal(const Vector<T,3>& normal, const Vector<T,3>& distance);
  short evalSignForSignedDistanceFromPseudonormal(const Vector<T,3>& pseudonormal, const Vector<T,3>& distance);
  short evalSignForSignedDistanceFromWindingNumber(const Vector<T,3>& pt);
  short evalSignForSignedDistanceFromCache(const Vector<T,3>& pt);
 
  T _voxelSize;
  T _stlSize;
  T _overlap;
  Octree<T>* _tree;
  const std::string _fName;
  STLmesh<T> _mesh;
  bool _verbose;
  mutable OstreamManager clout;
 
public:
  STLreader(const std::string fName, T voxelSize, T stlSize=1, RayMode method = RayMode::FastRayZ,
            bool verbose = false, T overlap=0., T max=0.);
  STLreader(const std::vector<std::vector<T>> meshPoints, T voxelSize, T stlSize=1, RayMode method= RayMode::FastRayZ,
            bool verbose = false, T overlap=0., T max=0.);
 
  ~STLreader() override;
  bool operator() (bool output[], const T input[]) override;
 
  bool isInsideRootTree(const T input[]);
 
  Vector<T,3> closestPointInBoundingBox(const Vector<T,3>& input);
 
  bool distance(T& distance,const Vector<T,3>& origin, const Vector<T,3>& direction, int iC=-1) override;
 
  T signedDistance(const Vector<T,3>& input) override;
  T signedDistanceExact(const Vector<T,3>& input) override;
  template <SignMode SIGNMODE>
  T signedDistance(const Vector<T,3>& input);
 
  Vector<T,3> surfaceNormal(const Vector<T,3>& pos, const T meshSize=0) override;
  Vector<T,3> surfaceNormalExact(const Vector<T,3>& pos, const T meshSize=0) override;
  template <SignMode SIGNMODE>
  Vector<T,3> surfaceNormal(const Vector<T,3>& pos, const T meshSize=0);
 
  void print();
 
  void writeSTL(std::string stlName="");
 
  void writeOctree();
 
  void setNormalsOutside();
 
  void setBoundaryInsideNodes();
 
  inline Octree<T>* getTree() const
  {
    return _tree;
  };
 
 
  inline STLmesh<T>& getMesh()
  {
    return _mesh;
  };
 
};
 
}  // namespace olb
 
#endif