diff --git a/FacetData.h b/FacetData.h
index 0ee36f1ff83cec19bb5bd65c5b3fa02c13ff6c25..2f1a723159a1006a4c5a75222dfe91bf32bbd6d9 100644
--- a/FacetData.h
+++ b/FacetData.h
@@ -33,51 +33,44 @@ public:
     };
 };
 
-struct Facet : public RTPrimitive {
-    Facet() : RTPrimitive(), sh(0){ surf = nullptr; };
-    Facet(size_t nbIndex) : RTPrimitive(), sh(nbIndex) { surf = nullptr; };
-    ~Facet(){
+struct GeomPrimitive : public RTPrimitive {
+    GeomPrimitive() : RTPrimitive(), sh(0){ surf = nullptr; };
+    explicit GeomPrimitive(size_t nbIndex) : RTPrimitive(), sh(nbIndex){ surf = nullptr; };
+    ~GeomPrimitive() override{
         if (surf) {
             //delete surf;
             // don' t delete, origin is an unreferenced shared ptr
             surf = nullptr;
         }
     }
+    size_t globalId; //Global index (to identify when superstructures are present)
     FacetProperties sh;
     std::vector<size_t>      indices;          // Indices (Reference to geometry vertex)
     std::vector<Vector2d> vertices2;        // Vertices (2D plane space, UV coordinates)
     Surface* surf;
+};
+
+struct Facet : public GeomPrimitive {
+    Facet() : GeomPrimitive(){ };
+    Facet(size_t nbIndex) : GeomPrimitive(nbIndex) { };
+    virtual ~Facet(){ }
 
-    size_t globalId; //Global index (to identify when superstructures are present)
     //size_t iSCount{0};
 
-    void ComputeBB() { bb = sh.bb;};
+    void ComputeBB() override { bb = sh.bb;};
     bool Intersect(Ray &r) override;
 
 };
 
-struct TriangleFacet : public RTPrimitive {
-    TriangleFacet() : RTPrimitive(), sh(0){ surf = nullptr; };
-    TriangleFacet(size_t nbIndex) : RTPrimitive(), sh(nbIndex) { surf = nullptr; };
-    ~TriangleFacet(){
-        if (surf) {
-            //delete surf;
-            // don' t delete, origin is an unreferenced shared ptr
-            surf = nullptr;
-        }
-    }
-    FacetProperties sh;
-    std::vector<size_t>      indices;          // Indices (Reference to geometry vertex)
-    std::vector<Vector2d> vertices2;        // Vertices (2D plane space, UV coordinates)
+struct TriangleFacet : public GeomPrimitive {
+    TriangleFacet() : GeomPrimitive(3){ };
+    virtual ~TriangleFacet() { }
     std::vector<Vector3d>* vertices3;        // Vertices (2D plane space, UV coordinates)
     Vector2d* texCoord;
 
-    Surface* surf;
-
-    size_t globalId; //Global index (to identify when superstructures are present)
     //size_t iSCount{0};
 
-    void ComputeBB() { bb = sh.bb;};
+    void ComputeBB() override { bb = sh.bb;};
     bool Intersect(Ray &r) override;
 };
 
diff --git a/IntersectAABB_shared.cpp b/IntersectAABB_shared.cpp
index 6f097de437a5d121a0a2c84e9ebaac0f4278c633..668872707253ce40c4b9286c043a5591cd03784c 100644
--- a/IntersectAABB_shared.cpp
+++ b/IntersectAABB_shared.cpp
@@ -50,9 +50,9 @@ std::tuple<size_t,size_t,size_t> AABBNODE::FindBestCuttingPlane() {
 	size_t nbLeft, nbRight;
 
 	for (const auto& f : facets) {
-		if (f->sh.center.x > centerX) rightFromCenterZ++;
-		if (f->sh.center.y > centerY) rightFromCenterY++;
-		if (f->sh.center.z > centerZ) rightFromCenterX++;
+		if (f->prim->sh.center.x > centerX) rightFromCenterZ++;
+		if (f->prim->sh.center.y > centerY) rightFromCenterY++;
+		if (f->prim->sh.center.z > centerZ) rightFromCenterX++;
 	}
 
 	double deviationFromHalfHalf_X = fabs((double)rightFromCenterZ - (double)(facets.size()) / 2.0);
@@ -87,12 +87,12 @@ void AABBNODE::ComputeBB() {
 	bb.min=Vector3d(1e100,1e100,1e100);
 
 	for (const auto& f : facets) {
-		bb.min.x = std::min(f->sh.bb.min.x,bb.min.x);
-		bb.min.y = std::min(f->sh.bb.min.y, bb.min.y);
-		bb.min.z = std::min(f->sh.bb.min.z, bb.min.z);
-		bb.max.x = std::max(f->sh.bb.max.x, bb.max.x);
-		bb.max.y = std::max(f->sh.bb.max.y, bb.max.y);
-		bb.max.z = std::max(f->sh.bb.max.z, bb.max.z);
+		bb.min.x = std::min(f->prim->sh.bb.min.x,bb.min.x);
+		bb.min.y = std::min(f->prim->sh.bb.min.y, bb.min.y);
+		bb.min.z = std::min(f->prim->sh.bb.min.z, bb.min.z);
+		bb.max.x = std::max(f->prim->sh.bb.max.x, bb.max.x);
+		bb.max.y = std::max(f->prim->sh.bb.max.y, bb.max.y);
+		bb.max.z = std::max(f->prim->sh.bb.max.z, bb.max.z);
 	}
 
 }
@@ -120,7 +120,7 @@ AABBNODE *BuildAABBTree(const std::vector<SubprocessFacet *> &facets, const size
 		case 1: // yz
 			m = (newNode->bb.min.x + newNode->bb.max.x) / 2.0;
 			for (const auto& f : newNode->facets) {
-				if (f->sh.center.x > m) rList[nbr++] = f;
+				if (f->prim->sh.center.x > m) rList[nbr++] = f;
 				else                   lList[nbl++] = f;
 			}
 			break;
@@ -128,7 +128,7 @@ AABBNODE *BuildAABBTree(const std::vector<SubprocessFacet *> &facets, const size
 		case 2: // xz
 			m = (newNode->bb.min.y + newNode->bb.max.y) / 2.0;
 			for (const auto& f : newNode->facets) {
-				if (f->sh.center.y > m) rList[nbr++] = f;
+				if (f->prim->sh.center.y > m) rList[nbr++] = f;
 				else                   lList[nbl++] = f;
 			}
 			break;
@@ -136,7 +136,7 @@ AABBNODE *BuildAABBTree(const std::vector<SubprocessFacet *> &facets, const size
 		case 3: // xy
 			m = (newNode->bb.min.z + newNode->bb.max.z) / 2.0;
 			for (const auto& f : newNode->facets) {
-				if (f->sh.center.z > m) rList[nbr++] = f;
+				if (f->prim->sh.center.z > m) rList[nbr++] = f;
 				else                   lList[nbl++] = f;
 			}
 			break;
@@ -299,9 +299,9 @@ IntersectTree(MFSim::Particle &currentParticle, const AABBNODE &node, const Vect
 			if (f == lastHitBefore)
 				continue;
 
-			double det = Dot(f->sh.Nuv, rayDirOpposite);
+			double det = Dot(f->prim->sh.Nuv, rayDirOpposite);
 			// Eliminate "back facet"
-			if ((f->sh.is2sided) || (det > 0.0)) { //If 2-sided or if ray going opposite facet normal
+			if ((f->prim->sh.is2sided) || (det > 0.0)) { //If 2-sided or if ray going opposite facet normal
 
 				double u, v, d;
 				// Ray/rectangle instersection. Find (u,v,dist) and check 0<=u<=1, 0<=v<=1, dist>=0
@@ -309,27 +309,27 @@ IntersectTree(MFSim::Particle &currentParticle, const AABBNODE &node, const Vect
 				if (det != 0.0) {
 
 					double iDet = 1.0 / det;
-					Vector3d intZ = rayPos - f->sh.O;
+					Vector3d intZ = rayPos - f->prim->sh.O;
 
-					u = iDet * DET33(intZ.x, f->sh.V.x, rayDirOpposite.x,
-						intZ.y, f->sh.V.y, rayDirOpposite.y,
-						intZ.z, f->sh.V.z, rayDirOpposite.z);
+					u = iDet * DET33(intZ.x, f->prim->sh.V.x, rayDirOpposite.x,
+						intZ.y, f->prim->sh.V.y, rayDirOpposite.y,
+						intZ.z, f->prim->sh.V.z, rayDirOpposite.z);
 
 					if (u >= 0.0 && u <= 1.0) {
 
-						v = iDet * DET33(f->sh.U.x, intZ.x, rayDirOpposite.x,
-							f->sh.U.y, intZ.y, rayDirOpposite.y,
-							f->sh.U.z, intZ.z, rayDirOpposite.z);
+						v = iDet * DET33(f->prim->sh.U.x, intZ.x, rayDirOpposite.x,
+							f->prim->sh.U.y, intZ.y, rayDirOpposite.y,
+							f->prim->sh.U.z, intZ.z, rayDirOpposite.z);
 
 						if (v >= 0.0 && v <= 1.0) {
 
-							d = iDet * Dot(f->sh.Nuv, intZ);
+							d = iDet * Dot(f->prim->sh.Nuv, intZ);
 
 							if (d>0.0) {
 
 								// Now check intersection with the facet polygon (in the u,v space)
 								// This check could be avoided on rectangular facet.
-								if (IsInFacet(*f, u, v)) {
+								if (IsInFacet(*currentParticle.model->primitives[f->globalId], u, v)) {
 									bool hardHit;
 #if defined(MOLFLOW)
 									double time = currentParticle.particle.time + d / 100.0 / currentParticle.velocity;
@@ -385,7 +385,7 @@ IntersectTree(MFSim::Particle &currentParticle, const AABBNODE &node, const Vect
 	}
 }
 
-bool IsInFacet(const SubprocessFacet &f, const double &u, const double &v) {
+bool IsInFacet(const GeomPrimitive &f, const double &u, const double &v) {
 
 	/*
 
diff --git a/IntersectAABB_shared.h b/IntersectAABB_shared.h
index aa6d0a4c0c49e241fbdb75196bb864dad2c9e067..6f886a5d85de2aa01db4dcae67efd082660e4501 100644
--- a/IntersectAABB_shared.h
+++ b/IntersectAABB_shared.h
@@ -51,6 +51,6 @@ std::tuple<bool, SubprocessFacet *, double>
 Intersect(MFSim::Particle &currentParticle, const Vector3d &rayPos, const Vector3d &rayDir, const AABBNODE *bvh);
 /*bool Visible(Simulation *sHandle, Vector3d *c1, Vector3d *c2, SubprocessFacet *f1, SubprocessFacet *f2,
              CurrentParticleStatus &currentParticle);*/
-bool IsInFacet(const SubprocessFacet &f,const double &u,const double &v);
+bool IsInFacet(const GeomPrimitive &f, const double &u, const double &v);
 //Vector3d PolarToCartesian(const SubprocessFacet *const collidedFacet, const double& theta, const double& phi, const bool& reverse); //sets sHandle->currentParticle.direction
 //std::tuple<double, double> CartesianToPolar(const Vector3d& incidentDir, const Vector3d& normU, const Vector3d& normV, const Vector3d& normN);
\ No newline at end of file
diff --git a/RayTracing/BVH.cpp b/RayTracing/BVH.cpp
index 67835f02a8e3b24a431d5fbfa25682f58b3d8662..7d9d965da123451c597d0a62448808098e0fe5e5 100644
--- a/RayTracing/BVH.cpp
+++ b/RayTracing/BVH.cpp
@@ -361,7 +361,7 @@ BVHAccel::BVHAccel(std::vector<std::shared_ptr<Primitive>> p,
         if (primitives.size() > probabilities.size())
             return;
         for (size_t i = 0; i < primitives.size(); ++i)
-            primitiveInfo[i] = {i, primitives[i]->sh.bb, probabilities[primitives[i]->globalId]};
+            primitiveInfo[i] = {i, primitives[i]->sh.bb, probabilities[i]};
     } else {
         for (size_t i = 0; i < primitives.size(); ++i)
             primitiveInfo[i] = {i, primitives[i]->sh.bb};
diff --git a/RayTracing/BVH.h b/RayTracing/BVH.h
index 52aad4328b40b248b94b4df0878b9d52e9798b80..5f3ba38f030bb52197aeb2dcde9156d69d75c344 100644
--- a/RayTracing/BVH.h
+++ b/RayTracing/BVH.h
@@ -10,7 +10,7 @@
 #include <FacetData.h>
 #include "Primitive.h"
 
-using Primitive = Facet;
+using Primitive = GeomPrimitive;
 
 struct BVHBuildNode;
 
diff --git a/RayTracing/KDTree.h b/RayTracing/KDTree.h
index 7d4bad1212490661c3f67ca831f6ebd12f2c03b0..6e6518d2e2deb710077d44037eb7dccfbe62d9b3 100644
--- a/RayTracing/KDTree.h
+++ b/RayTracing/KDTree.h
@@ -10,7 +10,7 @@
 #include <FacetData.h>
 #include "Primitive.h"
 
-using Primitive = Facet;
+using Primitive = GeomPrimitive;
 
 // KdTreeAccel Forward Declarations
 struct KdAccelNode;
diff --git a/Vector.cpp b/Vector.cpp
index b137044ee34484dbe7fc40804367687085e90e71..70e872a0d408b63e915fe3647728878dfcb46b67 100644
--- a/Vector.cpp
+++ b/Vector.cpp
@@ -19,7 +19,8 @@ Full license text: https://www.gnu.org/licenses/old-licenses/gpl-2.0.en.html
 */
 #include "Vector.h"
 #include "Helper/MathTools.h" //PI
-#include <math.h> //sqrt
+#include <cmath> //sqrt
+#include <array>
 
 Vector2d::Vector2d() {}