Implementation of the TGeoTessellated navigation functionality

geom/gdml/src/TGDMLParse.cxx

@@ -114,6 +114,9 @@ When most solids or volumes are added to the geometry they
 #include "TGeoSystemOfUnits.h"
 #include "TGeant4SystemOfUnits.h"
 
+#include "TGeoTypedefs.h"
+#include "Tessellated/TGeoMeshLoading.h"
+
 #include <cstdlib>
 #include <string>
 #include <sstream>
@@ -4311,6 +4314,7 @@ XMLNodePointer_t TGDMLParse::Tessellated(TXMLEngine *gdml, XMLNodePointer_t node
       local_name = TString::Format("%s_%s", name.Data(), fCurrentFile);
 
    auto tsl = new TGeoTessellated(NameShort(name));
+   Tessellated::MeshBuilder tslbuilder;
    TGeoTranslation *pos = nullptr;
    Tessellated::Vertex_t vertices[4];
 
@@ -4326,7 +4330,7 @@ XMLNodePointer_t TGDMLParse::Tessellated(TXMLEngine *gdml, XMLNodePointer_t node
          vertices[2] += vertices[0] + vertices[1];
          vertices[1] += vertices[0];
       }
-      tsl->AddFacet(vertices[0], vertices[1], vertices[2]);
+      tslbuilder.AddFacet(vertices[0], vertices[1], vertices[2]);
    };
 
    auto AddQuadrangularFacet = [&](bool relative) {
@@ -4335,7 +4339,7 @@ XMLNodePointer_t TGDMLParse::Tessellated(TXMLEngine *gdml, XMLNodePointer_t node
          vertices[2] += vertices[0] + vertices[1];
          vertices[1] += vertices[0];
       }
-      tsl->AddFacet(vertices[0], vertices[1], vertices[2], vertices[3]);
+      tslbuilder.AddFacet(vertices[0], vertices[1], vertices[2], vertices[3]);
    };
 
    // Get facet attributes
@@ -4436,7 +4440,11 @@ XMLNodePointer_t TGDMLParse::Tessellated(TXMLEngine *gdml, XMLNodePointer_t node
       }
       child = gdml->GetNext(child);
    }
-   tsl->CloseShape(false);
+   auto mesh = tslbuilder.CreateMesh();
+   if (!mesh->CheckClosure(false, false)) {
+      Fatal("Tessellated", "Mesh is not closed");
+   }
+   tsl->SetMesh(std::move(mesh));
 
    fsolmap[local_name.Data()] = tsl;
 

geom/gdml/src/TGDMLWrite.cxx

@@ -160,6 +160,7 @@ See that function for details.
 #include "TGeoOpticalSurface.h"
 #include "TMap.h"
 #include "TObjString.h"
+#include "Tessellated/TGeoTriangleMesh.h"
 #include <cstdlib>
 #include <string>
 #include <map>
@@ -1755,13 +1756,14 @@ XMLNodePointer_t TGDMLWrite::CreateTessellatedN(TGeoTessellated *geoShape)
 {
    // add all vertices to the define section
    TString genname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
-   for (int i = 0; i < geoShape->GetNvertices(); ++i) {
-      auto vertex = geoShape->GetVertex(i);
-      TString posName = TString::Format("%s_%d", genname.Data(), i);
+   const Tessellated::TGeoTriangleMesh *mesh = geoShape->GetTriangleMesh();
+   for (size_t i = 0; i < mesh->GetNumberOfVertices(); ++i) {
+      auto vertex = mesh->Point(i);
+      TString posName = TString::Format("%s_%zu", genname.Data(), i);
       Xyz nodPos;
-      nodPos.x = vertex[0];
-      nodPos.y = vertex[1];
-      nodPos.z = vertex[2];
+      nodPos.x = vertex.X();
+      nodPos.y = vertex.Y();
+      nodPos.z = vertex.Z();
       auto childN = CreatePositionN(posName.Data(), nodPos, "position", fDefault_lunit);
       fGdmlE->AddChild(fDefineNode, childN); // adding node to <define> node
    }
@@ -1770,17 +1772,14 @@ XMLNodePointer_t TGDMLWrite::CreateTessellatedN(TGeoTessellated *geoShape)
    fGdmlE->NewAttr(mainN, nullptr, "lunit", fDefault_lunit);
 
    XMLNodePointer_t childN;
-   for (Int_t it = 0; it < geoShape->GetNfacets(); it++) {
+   for (size_t it = 0; it < mesh->GetNumberOfTriangles(); it++) {
       // add section child node
-      auto facet = geoShape->GetFacet(it);
-      bool triangular = facet.GetNvert() == 3;
-      TString ntype = (triangular) ? "triangular" : "quadrangular";
+      auto facet = mesh->TriangleAt(it);
+      TString ntype = "triangular";
       childN = fGdmlE->NewChild(nullptr, nullptr, ntype.Data(), nullptr);
-      fGdmlE->NewAttr(childN, nullptr, "vertex1", TString::Format("%s_%d", genname.Data(), facet[0]));
-      fGdmlE->NewAttr(childN, nullptr, "vertex2", TString::Format("%s_%d", genname.Data(), facet[1]));
-      fGdmlE->NewAttr(childN, nullptr, "vertex3", TString::Format("%s_%d", genname.Data(), facet[2]));
-      if (!triangular)
-         fGdmlE->NewAttr(childN, nullptr, "vertex4", TString::Format("%s_%d", genname.Data(), facet[3]));
+      fGdmlE->NewAttr(childN, nullptr, "vertex1", TString::Format("%s_%d", genname.Data(), facet.Index(0)));
+      fGdmlE->NewAttr(childN, nullptr, "vertex2", TString::Format("%s_%d", genname.Data(), facet.Index(1)));
+      fGdmlE->NewAttr(childN, nullptr, "vertex3", TString::Format("%s_%d", genname.Data(), facet.Index(2)));
       fGdmlE->NewAttr(childN, nullptr, "type", "ABSOLUTE");
       fGdmlE->AddChild(mainN, childN);
    }

geom/geom/CMakeLists.txt

@@ -67,6 +67,13 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Geom
     TVirtualGeoPainter.h
     TVirtualGeoTrack.h
     TVirtualMagField.h
+    Tessellated/TBVH.h
+    Tessellated/TGeoMeshLoading.h
+    Tessellated/TOctant.h
+    Tessellated/TOctree.h
+    Tessellated/TGeoTriangle.h
+    Tessellated/TGeoTriangleMesh.h
+    Tessellated/TPartitioningI.h
   SOURCES
     src/TGDMLMatrix.cxx
     src/TGeoArb8.cxx
@@ -121,6 +128,12 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Geom
     src/TVirtualGeoPainter.cxx
     src/TVirtualGeoTrack.cxx
     src/TVirtualMagField.cxx
+    src/Tessellated/TBVH.cxx
+    src/Tessellated/TGeoMeshLoading.cxx
+    src/Tessellated/TOctant.cxx
+    src/Tessellated/TOctree.cxx
+    src/Tessellated/TGeoTriangle.cxx
+    src/Tessellated/TGeoTriangleMesh.cxx
   DEPENDENCIES
     Thread
     RIO

geom/geom/inc/LinkDef1.h

@@ -76,8 +76,13 @@
 #pragma link C++ class TGeoPolygon + ;
 #pragma link C++ class TGeoXtru + ;
 #pragma link C++ class ROOT::Geom::Vertex_t + ;
-#pragma link C++ class TGeoFacet + ;
-#pragma link C++ class TGeoTessellated + ;
+#pragma link C++ class TGeoTessellated - ;
+#pragma link C++ class Tessellated::TGeoTriangleMesh - ;
+#pragma link C++ class Tessellated::TGeoTriangle + ;
+#pragma link C++ class Tessellated::TOctree + ;
+#pragma link C++ class Tessellated::TOctant + ;
+#pragma link C++ class Tessellated::TPartitioningI + ;
+#pragma link C++ class Tessellated::TBVH + ;
 #pragma link C++ class TGeoShapeAssembly + ;
 #pragma link C++ class TGeoScaledShape + ;
 #pragma link C++ class TGeoVolume - ;

geom/geom/inc/TGeoTessellated.h

@@ -1,5 +1,4 @@
-// @(#)root/geom:$Id$
-// Author: Andrei Gheata   20/12/19
+// @(#)root/geom:$Id$// Author: Andrei Gheata   20/12/19
 
 /*************************************************************************
  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
@@ -12,125 +11,90 @@
 #ifndef ROOT_TGeoTessellated
 #define ROOT_TGeoTessellated
 
-#include <map>
-#include "TGeoVector3.h"
-#include "TGeoTypedefs.h"
-#include "TGeoBBox.h"
+#include <memory> // for unique_ptr
+#include <vector> // for vector
 
-class TGeoFacet {
-public:
-   using Vertex_t = Tessellated::Vertex_t;
-   using VertexVec_t = Tessellated::VertexVec_t;
-
-private:
-   int fIvert[4] = {0, 0, 0, 0}; // Vertex indices in the array
-   int fNvert = 0;               // number of vertices (can be 3 or 4)
+#include "Rtypes.h"     // for THashConsistencyHolder, ClassDefOverride
+#include "RtypesCore.h" // for Double_t, Bool_t, Int_t, UInt_t, kTRUE
+#include "TGeoBBox.h"   // for TGeoBBox
+#include "TStopwatch.h" // for TStopwatch
 
-private:
-   void SetVertices(int nvert = 0, int i0 = -1, int i1 = -1, int i2 = -1, int i3 = -1)
-   {
-      fNvert = nvert;
-      fIvert[0] = i0;
-      fIvert[1] = i1;
-      fIvert[2] = i2;
-      fIvert[3] = i3;
-   }
+#include "Tessellated/TPartitioningI.h"   // for TPartitioningI
+#include "Tessellated/TGeoTriangleMesh.h" // for TGeoTriangleMesh
 
-public:
-   TGeoFacet() {}
-   TGeoFacet(int i0, int i1, int i2) { SetVertices(3, i0, i1, i2); }
-   TGeoFacet(int i0, int i1, int i2, int i3) { SetVertices(4, i0, i1, i2, i3); }
-
-   int operator[](int ivert) const { return fIvert[ivert]; }
-   static int CompactFacet(Vertex_t *vert, int nvertices);
-   int GetNvert() const { return fNvert; }
-
-   void Flip()
-   {
-      int iv = fIvert[0];
-      fIvert[0] = fIvert[2];
-      fIvert[2] = iv;
-   }
-   bool IsNeighbour(const TGeoFacet &other, bool &flip) const;
-};
+class TBuffer3D;
+class TBuffer;
+class TClass;
+class TGeoMatrix;
+class TGeoShape;
+class TMemberInspector;
 
 class TGeoTessellated : public TGeoBBox {
 
-public:
-   using Vertex_t = Tessellated::Vertex_t;
+   using TPartitioningI = Tessellated::TPartitioningI;
+   using TGeoTriangleMesh = Tessellated::TGeoTriangleMesh;
+   using TGeoTriangle = Tessellated::TGeoTriangle;
 
 private:
-   int fNfacets = 0;                      // Number of facets
-   int fNvert = 0;                        // Number of vertices
-   int fNseg = 0;                         // Number of segments
-   bool fDefined = false;                 //! Shape fully defined
-   bool fClosedBody = false;              // The faces are making a closed body
-   std::vector<Vertex_t> fVertices;       // List of vertices
-   std::vector<TGeoFacet> fFacets;        // List of facets
-   std::multimap<long, int> fVerticesMap; //! Temporary map used to deduplicate vertices
-
-   TGeoTessellated(const TGeoTessellated &) = delete;
-   TGeoTessellated &operator=(const TGeoTessellated &) = delete;
+   std::unique_ptr<TGeoTriangleMesh> fMesh{nullptr};             ///<  triangle mesh
+   std::unique_ptr<TPartitioningI> fPartitioningStruct{nullptr}; ///<  partitioning structure
+   std::vector<UInt_t> fUsedTriangles{};                         ///<! vector of indices of valid triangles
+   mutable TStopwatch fTimer{}; ///<! timer to help determine timeconsuming TGeoTessellated instances
+   Bool_t fPrintTime{kFALSE};
 
-public:
-   // constructors
-   TGeoTessellated() {}
-   TGeoTessellated(const char *name, int nfacets = 0);
-   TGeoTessellated(const char *name, const std::vector<Vertex_t> &vertices);
-   // destructor
-   ~TGeoTessellated() override {}
-
-   void ComputeBBox() override;
-   void CloseShape(bool check = true, bool fixFlipped = true, bool verbose = true);
-
-   bool AddFacet(const Vertex_t &pt0, const Vertex_t &pt1, const Vertex_t &pt2);
-   bool AddFacet(const Vertex_t &pt0, const Vertex_t &pt1, const Vertex_t &pt2, const Vertex_t &pt3);
-   bool AddFacet(int i1, int i2, int i3);
-   bool AddFacet(int i1, int i2, int i3, int i4);
-   int AddVertex(const Vertex_t &vert);
-
-   bool FacetCheck(int ifacet) const;
-   Vertex_t FacetComputeNormal(int ifacet, bool &degenerated) const;
-
-   int GetNfacets() const { return fFacets.size(); }
-   int GetNsegments() const { return fNseg; }
-   int GetNvertices() const { return fNvert; }
-   bool IsClosedBody() const { return fClosedBody; }
-   Bool_t IsConvex() const final { return kFALSE; }
-   bool IsDefined() const { return fDefined; }
-
-   const TGeoFacet &GetFacet(int i) const { return fFacets[i]; }
-   const Vertex_t &GetVertex(int i) const { return fVertices[i]; }
+private:
+   void FillBuffer3DWithPoints(TBuffer3D &b) const;
+   void FillBuffer3DWithSegmentsAndPols(TBuffer3D &b, const std::vector<UInt_t> &indices) const;
+
+protected:
+   virtual void FillBuffer3D(TBuffer3D &b, Int_t reqSections, Bool_t localFrame) const override;
 
+public:
+   TGeoTessellated();
+   TGeoTessellated(const char *);
+   virtual ~TGeoTessellated() override;
+
+   virtual Bool_t Contains(const Double_t *pointa) const override;
+   virtual Double_t DistFromInside(const Double_t *pointa, const Double_t *dira, Int_t iact, Double_t step,
+                                   Double_t *safe) const override;
+   virtual Double_t DistFromOutside(const Double_t *pointa, const Double_t *dira, Int_t iact, Double_t step,
+                                    Double_t *safe) const override;
+   virtual Double_t Safety(const Double_t *pointa, bool inside) const override;
+
+   virtual void ComputeBBox() override;
+   virtual void ComputeNormal(const Double_t *pointa, const Double_t *dira, Double_t *norm) const override;
+   virtual void GetBoundingCylinder(Double_t *param) const override;
+   virtual void InspectShape() const override;
+   virtual Int_t GetByteCount() const override;
+   virtual Double_t Capacity() const override;
+   virtual TGeoShape *GetMakeRuntimeShape(TGeoShape *shape, TGeoMatrix *matrix) const override;
+   virtual Bool_t IsCylType() const override;
+   virtual Bool_t IsValidBox() const override;
+   virtual Bool_t GetPointsOnSegments(Int_t intl, Double_t *list) const override;
+   virtual TBuffer3D *MakeBuffer3D() const override;
+   virtual const TBuffer3D &GetBuffer3D(Int_t reqSections, Bool_t localFrame) const override;
    int DistancetoPrimitive(int, int) override { return 99999; }
-   const TBuffer3D &GetBuffer3D(int reqSections, Bool_t localFrame) const override;
-   void GetMeshNumbers(int &nvert, int &nsegs, int &npols) const override;
-   int GetNmeshVertices() const override { return fNvert; }
-   void InspectShape() const override {}
-   TBuffer3D *MakeBuffer3D() const override;
-   void Print(Option_t *option = "") const override;
    void SavePrimitive(std::ostream &, Option_t *) override {}
-   void SetPoints(double *points) const override;
-   void SetPoints(Float_t *points) const override;
-   void SetSegsAndPols(TBuffer3D &buff) const override;
    void Sizeof3D() const override {}
 
-   /// Resize and center the shape in a box of size maxsize
-   void ResizeCenter(double maxsize);
+   void SetPoints(Double_t *points) const override;
+   void SetPoints(Float_t *points) const override;
+   void GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const override;
+   Int_t GetNmeshVertices() const override { return GetTriangleMesh()->Points().size(); }
+
+   virtual void SetMesh(std::unique_ptr<TGeoTriangleMesh> mesh);
+   TGeoTriangleMesh const *GetTriangleMesh() const { return fMesh.get(); }
+   const std::vector<UInt_t> &GetUsedTriangleIndices() const { return fUsedTriangles; }
 
-   /// Flip all facets
-   void FlipFacets()
+   void SetPartitioningStruct(std::unique_ptr<TPartitioningI> &partitioningStruct)
    {
-      for (auto facet : fFacets)
-         facet.Flip();
+      fPartitioningStruct.reset(partitioningStruct.release());
    }
+   TPartitioningI const *GetPartitioningStruct() const { return fPartitioningStruct.get(); }
 
-   bool CheckClosure(bool fixFlipped = true, bool verbose = true);
-
-   /// Reader from .obj format
-   static TGeoTessellated *ImportFromObjFormat(const char *objfile, bool check = false, bool verbose = false);
-
-   ClassDefOverride(TGeoTessellated, 1) // tessellated shape class
+   void ResizeCenter(Double_t maxsize);
+   void PrintTime(Bool_t flag) { fPrintTime = flag; }
+   ClassDefOverride(TGeoTessellated, 1)
 };
 
-#endif
+#endif /*ROOT_TGeoTessellated*/

geom/geom/inc/TGeoTypedefs.h

@@ -14,6 +14,7 @@
 
 /// Typedefs used by the geometry group
 #include <vector>
+#include "TGeoVector3.h"
 
 namespace Tessellated {