From d8ae1c14a5c84dde6f8abdd4dd6955199cefa59a Mon Sep 17 00:00:00 2001
From: Markus Frank <Markus.Frank@cern.ch>
Date: Mon, 30 Aug 2021 13:28:11 +0200
Subject: [PATCH] Fix problems with exporting boolean shapes to CAD

---
 DDCAD/CMakeLists.txt               |   2 +-
 DDCAD/include/DDCAD/ASSIMPReader.h |   2 +
 DDCAD/include/DDCAD/ASSIMPWriter.h |  10 +-
 DDCAD/include/DDCAD/OutputWriter.h |  14 +-
 DDCAD/src/ASSIMPReader.cpp         |  67 +++++--
 DDCAD/src/ASSIMPWriter.cpp         | 297 +++++++++++++++++++++--------
 DDCAD/src/plugins/CADPlugins.cpp   |  23 ++-
 DDCore/include/XML/UnicodeValues.h |   1 +
 examples/DDCAD/CMakeLists.txt      |  22 ++-
 9 files changed, 323 insertions(+), 115 deletions(-)

diff --git a/DDCAD/CMakeLists.txt b/DDCAD/CMakeLists.txt
index 04d390539..a85572aff 100644
--- a/DDCAD/CMakeLists.txt
+++ b/DDCAD/CMakeLists.txt
@@ -95,7 +95,7 @@ target_include_directories(DDCAD
   $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include/>
   $<INSTALL_INTERFACE:include>
 )
-target_link_libraries(DDCAD PUBLIC DD4hep::DDCore assimp::assimp )
+target_link_libraries(DDCAD PUBLIC DD4hep::DDCore assimp::assimp ROOT::RCsg)
 
 dd4hep_add_plugin(DDCADPlugins
   SOURCES src/plugins/*.cpp
diff --git a/DDCAD/include/DDCAD/ASSIMPReader.h b/DDCAD/include/DDCAD/ASSIMPReader.h
index 80fdf5ada..73f10b816 100644
--- a/DDCAD/include/DDCAD/ASSIMPReader.h
+++ b/DDCAD/include/DDCAD/ASSIMPReader.h
@@ -34,6 +34,8 @@ namespace dd4hep {
      *  \ingroup DD4HEP_DDCAD
      */
     class ASSIMPReader : public InputReader   {
+    public:
+      long flags = 0;
     public:
       using InputReader::InputReader;
 
diff --git a/DDCAD/include/DDCAD/ASSIMPWriter.h b/DDCAD/include/DDCAD/ASSIMPWriter.h
index 4d08da066..7fb23190b 100644
--- a/DDCAD/include/DDCAD/ASSIMPWriter.h
+++ b/DDCAD/include/DDCAD/ASSIMPWriter.h
@@ -34,16 +34,18 @@ namespace dd4hep {
      *  \ingroup DD4HEP_DDCAD
      */
     class ASSIMPWriter : public OutputWriter   {
+    public:
+      long flags { 1 };
     public:
       using OutputWriter::OutputWriter;
       /// Default destructor
       virtual ~ASSIMPWriter() = default;
       /// Write output file
       virtual int write(const std::string& output_file,
-			const std::string& output_type,
-			const VolumePlacements& places,
-			bool   recursive,
-			double unit_scale = 1.0)  const  override;      
+                        const std::string& output_type,
+                        const VolumePlacements& places,
+                        bool   recursive,
+                        double unit_scale = 1.0)  const  override;      
     };
     
   }        /* End namespace cad                      */
diff --git a/DDCAD/include/DDCAD/OutputWriter.h b/DDCAD/include/DDCAD/OutputWriter.h
index 2f98ca215..8640646c1 100644
--- a/DDCAD/include/DDCAD/OutputWriter.h
+++ b/DDCAD/include/DDCAD/OutputWriter.h
@@ -41,14 +41,12 @@ namespace dd4hep {
     class OutputWriter   {
     public:
       enum output_flags   {
-	EXPORT_POINT_CLOUDS = 1 << 1,
-	LAST = 0
+                           EXPORT_POINT_CLOUDS = 1 << 1,
+                           LAST = 0
       };
     public:
       /// Reference to the detector object
       Detector& detector;
-      /// Output configuration flags
-      unsigned long flags  { 0 };
 
     public:
       typedef std::vector<PlacedVolume> VolumePlacements;
@@ -58,10 +56,10 @@ namespace dd4hep {
       virtual ~OutputWriter();
       /// Write output file
       virtual int write(const std::string& output_file,
-			const std::string& output_type,
-			const VolumePlacements& places,
-			bool   recursive,
-			double unit_scale = 1.0)  const = 0;      
+                        const std::string& output_type,
+                        const VolumePlacements& places,
+                        bool   recursive,
+                        double unit_scale = 1.0)  const = 0;      
     };
     
   }        /* End namespace cad                      */
diff --git a/DDCAD/src/ASSIMPReader.cpp b/DDCAD/src/ASSIMPReader.cpp
index 7d0c7e25b..42dd39746 100644
--- a/DDCAD/src/ASSIMPReader.cpp
+++ b/DDCAD/src/ASSIMPReader.cpp
@@ -23,6 +23,9 @@
 #include "assimp/scene.h"
 #include "assimp/postprocess.h"
 
+/// C/C++ include files
+#include <sstream>
+
 using namespace dd4hep;
 using namespace dd4hep::cad;
 
@@ -33,12 +36,13 @@ ASSIMPReader::readShapes(const std::string& source, double unit_length)  const
   typedef TessellatedSolid::Vertex Vertex;
   std::vector<std::unique_ptr<TGeoTessellated> > result;
   std::unique_ptr<Assimp::Importer> importer = std::make_unique<Assimp::Importer>();
-  int flags = aiProcess_Triangulate|aiProcess_JoinIdenticalVertices|aiProcess_CalcTangentSpace;
-  auto scene = importer->ReadFile( source.c_str(), flags);
+  int aiflags = aiProcess_Triangulate|aiProcess_JoinIdenticalVertices|aiProcess_CalcTangentSpace;
+  auto scene = importer->ReadFile( source.c_str(), aiflags);
   if ( !scene )  {
     except("ASSIMPReader","+++ FileNotFound: %s",source.c_str());
   }
   double unit = unit_length;
+  bool   dump_facets = ((flags>>8)&0x1) == 1;
   for (unsigned int index = 0; index < scene->mNumMeshes; index++)   {
     aiMesh* mesh = scene->mMeshes[index];
     if ( mesh->mNumFaces > 0 )   {
@@ -52,6 +56,14 @@ ASSIMPReader::readShapes(const std::string& source, double unit_length)  const
         Vertex b(v[idx[1]].x*unit, v[idx[1]].y*unit, v[idx[1]].z*unit); 
         Vertex c(v[idx[2]].x*unit, v[idx[2]].y*unit, v[idx[2]].z*unit); 
         shape->AddFacet(a,b,c);
+        if ( dump_facets )   {
+          size_t which = shape->GetNfacets()-1;
+          const auto& facet = shape->GetFacet(which);
+          std::stringstream str;
+          str << facet;
+          printout(ALWAYS,"ASSIMPReader","++ Facet %4ld : %s",
+                   which, str.str().c_str());
+        }
       }
       if ( shape->GetNfacets() > 2 )   {
         result.emplace_back(std::unique_ptr<TGeoTessellated>(shape.ptr()));
@@ -72,9 +84,10 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length)  const
   typedef TessellatedSolid::Vertex Vertex;
   std::vector<std::unique_ptr<TGeoVolume> > result;
   std::unique_ptr<Assimp::Importer> importer = std::make_unique<Assimp::Importer>();
-  int flags = aiProcess_Triangulate|aiProcess_JoinIdenticalVertices|aiProcess_CalcTangentSpace;
-  auto scene = importer->ReadFile( source.c_str(), flags);
-  char text[128];
+  bool dump_facets = ((flags>>8)&0x1) == 1;
+  int aiflags = aiProcess_Triangulate|aiProcess_JoinIdenticalVertices|aiProcess_CalcTangentSpace;
+  auto scene = importer->ReadFile( source.c_str(), aiflags);
+  char text[1048];
   
   if ( !scene )  {
     except("ASSIMPReader","+++ FileNotFound: %s",source.c_str());
@@ -83,9 +96,19 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length)  const
   for (unsigned int index = 0; index < scene->mNumMeshes; index++)   {
     aiMesh* mesh = scene->mMeshes[index];
     if ( mesh->mNumFaces > 0 )   {
-      auto name = mesh->mName.C_Str();
-      TessellatedSolid shape(name, mesh->mNumFaces);
+      std::string name = mesh->mName.C_Str();
+      std::vector<Vertex> vertices;
       const aiVector3D* v = mesh->mVertices;
+      for(unsigned int i=0; i < mesh->mNumVertices; i++)  {
+        vertices.emplace_back(Vertex(v[i].x*unit, v[i].y*unit, v[i].z*unit));
+      }
+      TessellatedSolid shape(name,vertices);
+      for(unsigned int i=0; i < mesh->mNumFaces; i++)  {
+        const unsigned int* idx  = mesh->mFaces[i].mIndices;
+        shape->AddFacet(idx[0], idx[1], idx[2]);
+      }
+#if 0
+      TessellatedSolid shape(name, mesh->mNumFaces);
       for(unsigned int i=0; i < mesh->mNumFaces; i++)  {
         const aiFace&     face = mesh->mFaces[i];
         const unsigned int* idx = face.mIndices;
@@ -94,24 +117,28 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length)  const
         Vertex c(v[idx[2]].x*unit, v[idx[2]].y*unit, v[idx[2]].z*unit); 
         shape->AddFacet(a,b,c);
       }
+#endif
       if ( shape->GetNfacets() > 2 )   {
-        std::string nam;
+        std::string mat_name;
         Material mat;
         VisAttr  vis;
         if ( scene->HasMaterials() )   {
           aiMaterial* ai_mat = scene->mMaterials[mesh->mMaterialIndex];
-          nam = ai_mat->GetName().C_Str();
-          mat = detector.material(nam);
+          mat_name = ai_mat->GetName().C_Str();
+          mat = detector.material(mat_name);
         }
         if ( !mat.isValid() )   {
           printout(ERROR, "ASSIMPReader",
                    "+++ %s: No material named '%s' FOUND. Will use Air. [Missing material]",
-                   text, nam.c_str());
+                   text, mat_name.c_str());
           mat = detector.air();
         }
-        ::snprintf(text,sizeof(text),"tessellated_%ld", result.size());
-        text[sizeof(text)-1] = 0;
-        Volume vol(text, Solid(shape.ptr()), mat);
+        if ( name.empty() )  {
+          ::snprintf(text,sizeof(text),"tessellated_%ld", result.size());
+          text[sizeof(text)-1] = 0;
+          name = text;
+        }
+        Volume vol(name, Solid(shape.ptr()), mat);
         if ( mesh->HasVertexColors(0) )   {
           const aiColor4D* col = mesh->mColors[0];
           if ( col )   {
@@ -125,7 +152,7 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length)  const
               }
             }
             if ( !vis.isValid() )   {
-              ::snprintf(text,sizeof(text),"vis_tessellated_%p", (void*)vol.ptr());
+              ::snprintf(text,sizeof(text),"vis_%s_%p", name.c_str(), (void*)vol.ptr());
               text[sizeof(text)-1] = 0;
               vis = VisAttr(text);
               vis.setColor(col->a,col->r,col->g,col->b);
@@ -137,6 +164,16 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length)  const
         printout(INFO,"ASSIMPReader",
                  "+++ %-17s Material: %-16s  Viualization: %s",
                  vol.name(), mat.name(), vis.isValid() ? vis.name() : "NONE");
+        shape->CloseShape(true,true,true);
+        if ( dump_facets )   {
+          for( size_t i=0, n=shape->GetNfacets(); i < n; ++i )   {
+            const auto& facet = shape->GetFacet(i);
+            std::stringstream str;
+            str << facet;
+            printout(ALWAYS,"ASSIMPReader","++ Facet %4ld : %s",
+                     i, str.str().c_str());
+          }
+        }
         result.emplace_back(std::unique_ptr<TGeoVolume>(vol.ptr()));
         continue;
       }
diff --git a/DDCAD/src/ASSIMPWriter.cpp b/DDCAD/src/ASSIMPWriter.cpp
index 85ac143be..972239eb1 100644
--- a/DDCAD/src/ASSIMPWriter.cpp
+++ b/DDCAD/src/ASSIMPWriter.cpp
@@ -14,6 +14,7 @@
 /// Framework include files
 #include <DD4hep/Shapes.h>
 #include <DD4hep/Objects.h>
+#include <DD4hep/Detector.h>
 #include <DD4hep/Printout.h>
 #include <DDCAD/ASSIMPWriter.h>
 
@@ -26,6 +27,9 @@
 #include <TBuffer3DTypes.h>
 #include <TGeoBoolNode.h>
 #include <TGeoMatrix.h>
+#include <TClass.h>
+#include <CsgOps.h>
+#include <GL/glu.h>
 
 /// C/C++ include files
 #include <set>
@@ -34,6 +38,8 @@ using namespace std;
 using namespace dd4hep;
 using namespace dd4hep::cad;
 
+using _Vertex = Tessellated::Vertex_t;
+
 namespace  {
 
   void _collect(std::vector<std::pair<PlacedVolume,TGeoHMatrix*> >& cont,
@@ -41,10 +47,8 @@ namespace  {
   {
     Volume v = pv.volume();
     for(Int_t i=0; i<v->GetNdaughters(); ++i)  {
-      PlacedVolume p = v->GetNode(i);
-      Solid sol = p.volume().solid();
-      // TessellatedSolid sol = p.volume().solid();
-      // if ( sol.isValid() ) cont.push_back(p);
+      PlacedVolume p   = v->GetNode(i);
+      Solid        sol = p.volume().solid();
       unique_ptr<TGeoHMatrix> mother(new TGeoHMatrix(to_global));
       mother->Multiply(p->GetMatrix());
 
@@ -56,27 +60,174 @@ namespace  {
         mother.release();
     }
   }
-  struct vertex{ double x,y,z;
-    vertex() = default;
-    vertex(const vertex& copy) = default;
-    vertex(double vx, double vy, double vz) : x(vx), y(vy), z(vz) {}
-    vertex(const Tessellated::Vertex_t& v) : x(v.x()), y(v.y()), z(v.z()) {}
-    vertex& operator=(const vertex& v) = default;
-    vertex& operator=(const Tessellated::Vertex_t& v)
-    { x = v.x(); y = v.y(); z = v.z(); return *this;  }
-    operator Tessellated::Vertex_t () const
-    { return Tessellated::Vertex_t(x,y,z); }
+  bool equals(_Vertex const &lhs, _Vertex const &rhs)  {
+    constexpr double kTolerance = 1.e-32;
+    return TMath::Abs(lhs[0] - rhs[0]) < kTolerance &&
+                                         TMath::Abs(lhs[1] - rhs[1]) < kTolerance &&
+                                                                       TMath::Abs(lhs[2] - rhs[2]) < kTolerance;
+  }
+
+  struct TessellateShape   {
+  public:
+    TessellateShape() = default;
+    virtual ~TessellateShape() = default;
+    RootCsg::TBaseMesh* make_mesh(TGeoShape* sh)  const;
+    RootCsg::TBaseMesh* collect_composite(TGeoCompositeShape* sh)    const;
+    unique_ptr<TGeoTessellated> build_mesh(int id, const std::string& name, TGeoShape* shape);
+    unique_ptr<TGeoTessellated> tessellate_primitive(const std::string& name, Solid solid);
+    unique_ptr<TGeoTessellated> close_tessellated(int id, TGeoShape* shape, int nskip, unique_ptr<TGeoTessellated>&& tes);
   };
-#if 0
-  vertex operator-(const vertex& v1, const vertex& v2)
-  {  return vertex(v1.x-v2.x, v1.y-v2.y, v1.z-v2.z);   }
-  vertex operator+(const vertex& v1, const vertex& v2)
-  {  return vertex(v1.x+v2.x, v1.y+v2.y, v1.z+v2.z);   }
+
+  unique_ptr<TGeoTessellated> TessellateShape::close_tessellated(int id, TGeoShape* shape, int nskip, unique_ptr<TGeoTessellated>&& tes)   {
+    string nam = shape->GetName(), typ = "["+string(shape->IsA()->GetName())+"]";
+    nam = nam.substr(0, nam.find("_0x"));
+    tes->CloseShape(true, true, true);
+    if ( nskip > 0 )   {
+      printout(ALWAYS,"ASSIMPWriter","+++ %3d %-48s %-24s Skipped %3ld/%-4d degenerate facets %4d vertices closed:%s defined:%s",
+               id, nam.c_str(), typ.c_str(), nskip, tes->GetNfacets(),  tes->GetNvertices(),
+               yes_no(tes->IsClosedBody()), yes_no(tes->IsDefined()));
+    }
+    else   {
+      printout(ALWAYS,"ASSIMPWriter","+++ %3d %-48s %-24s Tessellated %4d facets %4d vertices closed:%s defined:%s",
+               id, nam.c_str(), typ.c_str(),
+               tes->GetNfacets(), tes->GetNvertices(),
+               yes_no(tes->IsClosedBody()),
+               yes_no(tes->IsDefined()));
+    }
+    cout << flush;
+    return move(tes);
+  }
+
+  unique_ptr<TGeoTessellated> TessellateShape::build_mesh(int id, const std::string& name, TGeoShape* shape)      {
+    auto mesh = unique_ptr<RootCsg::TBaseMesh>(this->make_mesh(shape));
+    vector<_Vertex> vertices;
+    vertices.reserve(mesh->NumberOfVertices());
+    size_t nskip = 0;
+    map<size_t,size_t> vtx_index_replacements;
+    for( size_t ipoint = 0, npoints = mesh->NumberOfVertices(); ipoint < npoints; ++ipoint )   {
+      long found = -1;
+      const double* v = mesh->GetVertex(ipoint);
+      _Vertex vtx(v[0], v[1], v[2]);
+      for(size_t i=0; i < vertices.size(); ++i)   {
+        if ( equals(vertices[i],vtx) )  {
+          vtx_index_replacements[ipoint] = found = i;
+          break;
+        }
+      }
+      if ( found < 0 )   {
+        vtx_index_replacements[ipoint] = vertices.size();
+        vertices.emplace_back(v[0], v[1], v[2]);
+      }
+    }
+    size_t vtx_len = vertices.size();
+    unique_ptr<TGeoTessellated> tes = make_unique<TGeoTessellated>(name.c_str(), vertices);
+    for( size_t ipoly = 0, npols = mesh->NumberOfPolys(); ipoly < npols; ++ipoly)    {
+      size_t npoints = mesh->SizeOfPoly(ipoly);
+      if ( npoints >= 3 )   {
+        printout(DEBUG,"ASSIMPWriter","+++ Got polygon with %ld points",npoints);
+        ///
+        /// 3-vertex polygons automatically translate to GL_TRIANGLES
+        /// See Kronos documentation to glBegin / glEnd from the glu library:
+        /// https://www.khronos.org/registry/OpenGL-Refpages/gl2.1/xhtml/glBegin.xml
+        ///
+        /// Otherwise:
+#if 1
+        /// Apparently this is the correct choice:
+        ///
+        /// Interprete as FAN:  GL_TRIANGLE_FAN
+        /// One triangle is defined for each vertex presented after the first two vertices.
+        /// Vertices 1 , n + 1 , and n + 2 define triangle n.
+        /// N - 2 triangles are drawn.
+        size_t v0  = mesh->GetVertexIndex(ipoly, 0);
+        size_t vv0 = vtx_index_replacements[v0];
+        for( size_t ipoint = 0; ipoint < npoints-2; ++ipoint )   {
+          size_t v1 = mesh->GetVertexIndex(ipoly, ipoint+1);
+          size_t v2 = mesh->GetVertexIndex(ipoly, ipoint+2);
+          size_t vv1 = vtx_index_replacements[v1];
+          size_t vv2 = vtx_index_replacements[v2];
+          if ( vv0 > vtx_len || vv1 > vtx_len || vv2 > vtx_len )  {
+            ++nskip;
+            continue;
+          }
+          if ( vv0 == vv1 || vv0 == vv2 || vv1 == vv2 )   {
+            ++nskip;
+            continue;
+          }
+          _Vertex w[3] = {vertices[vv0],vertices[vv1],vertices[vv2]};
+          if ( TGeoFacet::CompactFacet(w, 3) < 3 )   {
+            ++nskip;
+            continue;
+          }
+          bool degenerated = true;
+          TGeoFacet f(&vertices, 3, vv0, vv1, vv2);
+          f.ComputeNormal(degenerated);
+          if ( degenerated )    {
+            ++nskip;
+            continue;
+          }
+          tes->AddFacet(vv0, vv1, vv2);
+        }
+#else
+        /// Interprete as STRIP: GL_TRIANGLE_STRIP
+        /// One triangle is defined for each vertex presented after the first two vertices.
+        /// For odd n, vertices n, n + 1 , and n + 2 define triangle n.
+        /// For even n, vertices n + 1 , n, and n + 2 define triangle n.
+        /// N - 2 triangles are drawn.
+        for( size_t ipoint = 0; ipoint < npoints-2; ++ipoint )   {
+          vtx_t v0(mesh->GetVertex(mesh->GetVertexIndex(ipoly, ipoint)));
+          vtx_t v1(mesh->GetVertex(mesh->GetVertexIndex(ipoly, ipoint+1)));
+          vtx_t v2(mesh->GetVertex(mesh->GetVertexIndex(ipoly, ipoint+2)));
+          ((ipoint%2) == 0) ? tes->AddFacet(v1, v0, v2) : tes->AddFacet(v0, v1, v2);
+        }
 #endif
+      }
+    }
+    return close_tessellated(id, shape, nskip, move(tes));
+  }
+  
+  RootCsg::TBaseMesh* TessellateShape::make_mesh(TGeoShape* sh)   const   {
+    if (TGeoCompositeShape *shape = dynamic_cast<TGeoCompositeShape *>(sh))   {
+      return collect_composite(shape);
+    }
+    UInt_t  flags = TBuffer3D::kCore|TBuffer3D::kBoundingBox|TBuffer3D::kRawSizes|TBuffer3D::kRaw|TBuffer3D::kShapeSpecific;
+    const TBuffer3D& buffer = sh->GetBuffer3D(flags, kFALSE);
+    return RootCsg::ConvertToMesh(buffer);
+  }
+  
+  RootCsg::TBaseMesh* TessellateShape::collect_composite(TGeoCompositeShape* sh)  const  {
+    TGeoBoolNode* node  = sh->GetBoolNode();
+    TGeoShape*    left  = node->GetLeftShape();
+    TGeoShape*    right = node->GetRightShape();
+    TGeoHMatrix*  glmat = (TGeoHMatrix*)TGeoShape::GetTransform();
+    UInt_t        oper  = node->GetBooleanOperator();
+    TGeoHMatrix   copy(*glmat); // keep a copy
+
+    // Do not wonder about this logic.
+    // GetBuffer3D (->make_mesh) uses static variable fgTransform of TGeoShape!
+    glmat->Multiply(node->GetLeftMatrix());
+    auto left_mesh  = unique_ptr<RootCsg::TBaseMesh>(make_mesh(left));
+    *glmat = &copy;
 
-  unique_ptr<TGeoTessellated> _tessellate_primitive(const std::string& name, Solid solid)   {
-    typedef vertex vtx_t;
-    const TBuffer3D& buf3D = solid->GetBuffer3D(TBuffer3D::kAll, false);
+    glmat->Multiply(node->GetRightMatrix());
+    auto right_mesh = unique_ptr<RootCsg::TBaseMesh>(make_mesh(right));
+    *glmat = &copy;
+
+    switch (oper) {
+    case TGeoBoolNode::kGeoUnion:
+      return RootCsg::BuildUnion(left_mesh.get(), right_mesh.get());
+    case TGeoBoolNode::kGeoIntersection:
+      return RootCsg::BuildIntersection(left_mesh.get(), right_mesh.get());
+    case TGeoBoolNode::kGeoSubtraction:
+      return RootCsg::BuildDifference(left_mesh.get(), right_mesh.get());
+    default:
+      Error("BuildComposite", "Wrong boolean operation code %d\n", oper);
+      return 0;
+    }
+  }
+
+  unique_ptr<TGeoTessellated> TessellateShape::tessellate_primitive(const std::string& name, Solid solid)   {
+    using  vtx_t = _Vertex;
+    const  TBuffer3D& buf3D = solid->GetBuffer3D(TBuffer3D::kAll, false);
     struct pol_t { int c, n; int segs[1]; } *pol = nullptr;
     struct seg_t { int c, _1, _2;         };
     const  seg_t* segs = (seg_t*)buf3D.fSegs;
@@ -125,42 +276,6 @@ namespace  {
     }
     return tes;
   }
-
-  struct TessellateComposite   {
-    TBuffer3D buffer    { TBuffer3DTypes::kComposite };
-    std::vector<unique_ptr<TGeoTessellated> > meshes;
-  public:
-    TessellateComposite() = default;
-    virtual ~TessellateComposite() = default;
-    void collect_mesh(TGeoShape* s)  {
-      if (TGeoCompositeShape *shape = dynamic_cast<TGeoCompositeShape *>(s))
-        collect_composite(shape);
-      else
-        meshes.emplace_back(_tessellate_primitive(s->GetName(),s));
-    }
-    void collect_composite(TGeoCompositeShape* sh)   {
-      TGeoBoolNode* node  = sh->GetBoolNode();
-      TGeoShape*    left  = node->GetLeftShape();
-      TGeoShape*    right = node->GetRightShape();
-      TGeoHMatrix  *glmat = (TGeoHMatrix*)sh->GetTransform();
-      TGeoHMatrix   mat;
-      mat = glmat; // keep a copy
-
-      glmat->Multiply(node->GetLeftMatrix());
-      collect_mesh(left);
-      *glmat = &mat;
-
-      glmat->Multiply(node->GetRightMatrix());
-      collect_mesh(right);
-      *glmat = &mat;
-    }
-    unique_ptr<TGeoTessellated> build_composite(const std::string& name, TGeoCompositeShape* shape)     {
-      int num_facet = 0;
-      this->collect_composite(shape);
-      unique_ptr<TGeoTessellated> tes = make_unique<TGeoTessellated>(name.c_str(), num_facet);
-      return tes;
-    }
-  };
 }
 
 /// Write output file
@@ -171,9 +286,11 @@ int ASSIMPWriter::write(const std::string& file_name,
                         double unit_scale)  const
 {
   std::vector<std::pair<PlacedVolume,TGeoHMatrix*> > placements;
+  int  build_mode  = ((flags&0x1) != 0) ? 1 : 0;
+  bool dump_facets = ((flags&0x2) != 0);
   vector<Material>        materials;
   TGeoHMatrix             toGlobal;
-    
+
   for( auto pv : places )
     _collect(placements, recursive, toGlobal, pv);
 
@@ -192,9 +309,12 @@ int ASSIMPWriter::write(const std::string& file_name,
   root->mNumChildren  = 0;
   root->mChildren     = new aiNode* [num_mesh];
   root->mMeshes       = 0;//new unsigned int[root->mNumMeshes];
+  auto* geo_transform = TGeoShape::GetTransform();
 
+  TGeoHMatrix identity;
   for( size_t imesh=0; imesh < num_mesh; ++imesh )   {
-    unique_ptr<TGeoHMatrix> trafo(placements[imesh].second);
+    unique_ptr<TGeoHMatrix>     trafo(placements[imesh].second);
+    unique_ptr<TGeoTessellated> shape_holder;
     PlacedVolume     pv  = placements[imesh].first;
     Volume           v   = pv.volume();
     Solid            s   = v.solid();
@@ -202,16 +322,23 @@ int ASSIMPWriter::write(const std::string& file_name,
     TessellatedSolid tes = s;
     aiString         node_name(v.name());
 
-    unique_ptr<TGeoTessellated> buf;
+    identity.Clear();
+    TGeoShape::SetTransform(&identity);
+
+    /// If the shape is already tessellated, no need to create another one!
     if ( !tes.isValid() )   {
-      if ( auto* shape=dynamic_cast<TGeoCompositeShape*>(s.ptr()) )   {
-        TessellateComposite helper;
-        buf = helper.build_composite(v.name(), shape);
+      TessellateShape helper;
+      auto* shape=dynamic_cast<TGeoCompositeShape*>(s.ptr());
+      if ( build_mode || shape )   {  // Always use this method!
+        auto* paintVol = detector.manager().GetPaintVolume();
+        detector.manager().SetPaintVolume(v.ptr());
+        shape_holder = helper.build_mesh(imesh, v.name(), s.ptr());
+        detector.manager().SetPaintVolume(paintVol);
       }
       else   {
-        buf = _tessellate_primitive(v.name(), s);
+        shape_holder = helper.tessellate_primitive(v.name(), s);
       }
-      tes = buf.get();
+      tes = shape_holder.get();
     }
     if ( tes->GetNfacets() == 0 )   {
       continue;
@@ -258,7 +385,7 @@ int ASSIMPWriter::write(const std::string& file_name,
       face.mNumIndices = 0;
       face.mIndices = nullptr;
     }
-    vertex vtx, tmp, norm;
+    _Vertex vtx, tmp, norm;
     for( long j=0, nvx=0, n=tes->GetNfacets(); j < n; ++j )  {
       bool degenerated  = false;
       const auto& facet = tes->GetFacet(j);
@@ -268,28 +395,36 @@ int ASSIMPWriter::write(const std::string& file_name,
         double  u     = unit_scale;
 
         face.mIndices = new unsigned int[facet.GetNvert()];
-        trafo->LocalToMaster(&tmp.x, &norm.x);
+        trafo->LocalToMaster(tmp.fVec, norm.fVec);
         face.mNumIndices = 0;
         for( long k=0; k < facet.GetNvert(); ++k )  {
           tmp = facet.GetVertex(k);
-          trafo->LocalToMaster(&tmp.x, &vtx.x);
+          trafo->LocalToMaster(tmp.fVec, vtx.fVec);
           face.mIndices[face.mNumIndices] = nvx;
-          mesh->mNormals[nvx]  = aiVector3D(norm.x, norm.y, norm.z);
-          mesh->mVertices[nvx] = aiVector3D(vtx.x*u,vtx.y*u,vtx.z*u);
+          mesh->mNormals[nvx]  = aiVector3D(norm.x(), norm.y(), norm.z());
+          mesh->mVertices[nvx] = aiVector3D(vtx.x()*u,vtx.y()*u,vtx.z()*u);
           ++mesh->mNumVertices;
           ++face.mNumIndices;
           ++nvx;
         }
         ++mesh->mNumFaces;
+        if ( dump_facets )   {
+          stringstream str;
+          const auto* id = face.mIndices;
+          const auto* vv = mesh->mVertices;
+          TGeoFacet fac(_Vertex(vv[id[0]].x,vv[id[0]].y,vv[id[0]].z),
+                        _Vertex(vv[id[1]].x,vv[id[1]].y,vv[id[1]].z),
+                        _Vertex(vv[id[2]].x,vv[id[2]].y,vv[id[2]].z));
+          str << fac;
+          printout(ALWAYS,"ASSIMPWriter","++ Facet %4ld : %s", j, str.str().c_str());
+        }
+      }
+      else   {
+        printout(ALWAYS,"ASSIMPWriter",
+                 "+++ Found degenerate facet while writing [Should not happen]");
       }
     }
-    if ( imesh == 122585 )
-      break;
-    else if ( imesh == 122586 )
-      mesh->mNumFaces = 0;
-    else if ( imesh == 122587 )
-      mesh->mNumFaces = 0;
-
+    
     /// Check if we have here a valid mesh
     if ( 0 == mesh->mNumFaces || 0 == mesh->mNumVertices )    {
       delete [] mesh->mVertices;
@@ -312,13 +447,13 @@ int ASSIMPWriter::write(const std::string& file_name,
     ++root->mNumChildren;
     ++scene.mNumMeshes;
   }
-  printout(ALWAYS,"ASSIMPWriter","+++ Analysed %ld of %ld meshes.",
+  TGeoShape::SetTransform(geo_transform);
+  printout(ALWAYS,"ASSIMPWriter","+++ Analysed %ld out of %ld meshes.",
            scene.mNumMeshes, placements.size());
   if ( scene.mNumMeshes > 0 )   {
     Assimp::Exporter exporter;
     Assimp::ExportProperties *props = new Assimp::ExportProperties;
-    props->SetPropertyBool(AI_CONFIG_EXPORT_POINT_CLOUDS,
-                           flags&EXPORT_POINT_CLOUDS ? true : false);
+    props->SetPropertyBool(AI_CONFIG_EXPORT_POINT_CLOUDS, flags&EXPORT_POINT_CLOUDS ? true : false);
     exporter.Export(&scene, file_type.c_str(), file_name.c_str(), 0, props);
     return 1;
   }
diff --git a/DDCAD/src/plugins/CADPlugins.cpp b/DDCAD/src/plugins/CADPlugins.cpp
index 7fdc25c0f..3d4820d7c 100644
--- a/DDCAD/src/plugins/CADPlugins.cpp
+++ b/DDCAD/src/plugins/CADPlugins.cpp
@@ -60,9 +60,13 @@ DECLARE_DD4HEP_CONSTRUCTOR(DD4hep_read_CAD_volumes,read_CAD_Volume)
 
 static Handle<TObject> create_CAD_Shape(Detector& dsc, xml_h e)   {
   xml_elt_t elt(e);
+  cad::ASSIMPReader rdr(dsc);
   string fname = elt.attr<string>(_U(ref));
-  double unit  = elt.hasAttr(_U(unit)) ? elt.attr<double>(_U(unit)) : dd4hep::cm;
-  auto shapes = cad::ASSIMPReader(dsc).readShapes(fname, unit);
+  long   flags = elt.hasAttr(_U(flags)) ? elt.attr<long>(_U(flags))  : 0;
+  double unit  = elt.hasAttr(_U(unit))  ? elt.attr<double>(_U(unit)) : dd4hep::cm;
+
+  if ( flags ) rdr.flags = flags;
+  auto shapes = rdr.readShapes(fname, unit);
   if ( shapes.empty() )   {
     except("CAD_Shape","+++ CAD file: %s does not contain any "
            "understandable tessellated shapes.", fname.c_str());
@@ -149,13 +153,19 @@ DECLARE_XML_VOLUME(CAD_Assembly__volume_constructor,create_CAD_Assembly)
  *       <physvolid name="slice" value="10"/>
  *     </volume>
  *
+ *     If flags: (flags>>8)&1 == 1 (257): dump facets
+ *
  *   </XXX>
  */
 static Handle<TObject> create_CAD_Volume(Detector& dsc, xml_h e)   {
   xml_elt_t elt(e);
   string fname = elt.attr<string>(_U(ref));
   double unit  = elt.attr<double>(_U(unit));
-  auto volumes = cad::ASSIMPReader(dsc).readVolumes(fname, unit);
+  long   flags = elt.hasAttr(_U(flags)) ? elt.attr<long>(_U(flags))  : 0;
+  cad::ASSIMPReader rdr(dsc);
+
+  if ( flags ) rdr.flags = flags;
+  auto volumes = rdr.readVolumes(fname, unit);
   if ( volumes.empty() )   {
     except("CAD_Volume","+++ CAD file: %s does not contain any "
            "understandable tessellated volumes.", fname.c_str());
@@ -257,9 +267,10 @@ DECLARE_XML_VOLUME(CAD_MultiVolume__volume_constructor,create_CAD_Volume)
  *
  */
 static long CAD_export(Detector& description, int argc, char** argv)   {
-  bool recursive = false, help = false;
+  bool   recursive = false, help = false;
   string volume, detector, fname, ftype;
   double scale = 1.0;
+  int    flags = 0;
   
   for(int i = 0; i < argc && argv[i]; ++i)  {
     if (      0 == ::strncmp( "-output",argv[i],4) )    fname     = argv[++i];
@@ -274,6 +285,8 @@ static long CAD_export(Detector& description, int argc, char** argv)   {
     else if ( 0 == ::strncmp("--recursive",argv[i],5) ) recursive = true;
     else if ( 0 == ::strncmp( "-scale",argv[i],4) )     scale     = ::atof(argv[++i]);
     else if ( 0 == ::strncmp("--scale",argv[i],5) )     scale     = ::atof(argv[++i]);
+    else if ( 0 == ::strncmp( "-flags",argv[i],4) )     flags     = ::atol(argv[++i]);
+    else if ( 0 == ::strncmp("--flags",argv[i],5) )     flags     = ::atol(argv[++i]);
     else if ( 0 == ::strncmp( "-help",argv[i],2) )      help      = true;
     else if ( 0 == ::strncmp("--help",argv[i],3) )      help      = true;
   }
@@ -290,6 +303,7 @@ static long CAD_export(Detector& description, int argc, char** argv)   {
       "     -detector <string> Path to the detector element to be exported.     \n"
       "     -help              Print this help output.                          \n"
       "     -scale    <number> Unit scale before writing output data.           \n"
+      "     -flags    <number> Flagsging helper to pass args -- Experts only.   \n"
       "     Arguments given: " << arguments(argc,argv) << endl << flush;
     ::exit(EINVAL);
   }
@@ -316,6 +330,7 @@ static long CAD_export(Detector& description, int argc, char** argv)   {
     }
   }
   cad::ASSIMPWriter wr(description);
+  if ( flags ) wr.flags = flags;
   std::vector<PlacedVolume> places {pv};
   auto num_mesh = wr.write(fname, ftype, places, recursive, scale);
   if ( num_mesh < 0 )   {
diff --git a/DDCore/include/XML/UnicodeValues.h b/DDCore/include/XML/UnicodeValues.h
index f45140618..dd01285c0 100644
--- a/DDCore/include/XML/UnicodeValues.h
+++ b/DDCore/include/XML/UnicodeValues.h
@@ -165,6 +165,7 @@ UNICODE (finish);
 UNICODE (first);
 UNICODE (firstposition);
 UNICODE (firstrotation);
+UNICODE (flags);
 UNICODE (formula);
 UNICODE (fraction);
 UNICODE (fractions);
diff --git a/examples/DDCAD/CMakeLists.txt b/examples/DDCAD/CMakeLists.txt
index 29f85ec64..0a728f585 100644
--- a/examples/DDCAD/CMakeLists.txt
+++ b/examples/DDCAD/CMakeLists.txt
@@ -65,7 +65,7 @@ dd4hep_add_test_reg( DDCAD_export_sid_vertex
   	                  -input ${CLICSiDEx_INSTALL}/compact/SiD_Vertex.xml
 			  -plugin DD4hep_CAD_export -output clicsid_vertex.collada
 			  -type collada -recursive -detector /world -recursive -scale 1.0 
-  REGEX_PASS "Analysed 383 of 383 meshes"
+  REGEX_PASS "Analysed 383 out of 383 meshes"
   REGEX_FAIL "Exception"
   REGEX_FAIL "FAILED"
 )
@@ -86,7 +86,7 @@ dd4hep_add_test_reg( DDCAD_export_cal_endcaps
   EXEC_ARGS  geoPluginRun -input ${ClientTestsEx_INSTALL}/compact/CaloEndcapReflection.xml
   -plugin DD4hep_CAD_export -output endcap_reflection.collada
   -type collada -recursive -detector /world -recursive -scale 1.0
-  REGEX_PASS "Analysed 96 of 96 meshes"
+  REGEX_PASS "Analysed 96 out of 96 meshes"
   REGEX_FAIL "Exception"
   REGEX_FAIL "FAILED"
 )
@@ -100,6 +100,24 @@ dd4hep_add_test_reg( DDCAD_import_cal_endcaps
   REGEX_FAIL "Exception"
 )
 #
+#  Test CAD export of FCC machine part
+dd4hep_add_test_reg( DDCAD_export_FCC_machine
+  COMMAND    "${CMAKE_INSTALL_PREFIX}/bin/run_test_DDCAD.sh"
+  EXEC_ARGS  geoPluginRun -input ${ClientTestsEx_INSTALL}/compact/FCCmachine/FCCee_DectMaster.xml
+  -plugin DD4hep_CAD_export -output Machine.collada
+  -type collada -recursive -detector /world -recursive -scale 1.0
+  REGEX_PASS "Analysed 49 out of 49 meshes"
+  REGEX_FAIL "Exception"
+  REGEX_FAIL "FAILED"
+)
 #
+#  Test CAD import of FCC machine part
+dd4hep_add_test_reg( DDCAD_import_FCC_machine
+  COMMAND    "${CMAKE_INSTALL_PREFIX}/bin/run_test_DDCAD.sh"
+  EXEC_ARGS  geoDisplay -input ${DDCADEx_INSTALL}/compact/Import_FCC_Machine.xml -load -destroy
+  DEPENDS    DDCAD_export_FCC_machine
+  REGEX_PASS "Read 49 meshes"
+  REGEX_FAIL "Exception"
+)
 #
 
-- 
GitLab