Improved Geometry's dimensionality info.

Fixed windows & linux bug with gl where a glBindVertexArray would error.
Fixed windows bug with qhull where switching between debug and release in VS would somehow link/build against the wrong qhull libs.
Renamed Widgets' update functions to updateData, so as to separate them from the QT base class function.

core/CMakeLists.txt

@@ -38,7 +38,7 @@ include(ExternalProject)
 ExternalProject_add(qhullstatic_r
 	SOURCE_DIR    "${CMAKE_SOURCE_DIR}/thirdparty/qhull"
 	# GIT_REPOSITORY https://github.com/qhull/qhull
-	CMAKE_ARGS "-DCMAKE_INSTALL_PREFIX=${PROJECT_BINARY_DIR}/thirdparty-install;-Dqhull_TARGETS_INSTALL=qhullstatic_r"
+	CMAKE_ARGS "-DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE};-DCMAKE_INSTALL_PREFIX=${PROJECT_BINARY_DIR}/thirdparty-install;-Dqhull_TARGETS_INSTALL=qhullstatic_r"
 	PREFIX qhull-prefix
 )
 

core/include/data/Geometry.hpp

@@ -66,7 +66,8 @@ namespace Data
 
     const std::vector<tetrahedron_t>& triangulation();
 
-    bool is2D() const;
+    int dimensionality;
+	int calculateDimensionality() const;
 
   private:
     std::vector<tetrahedron_t> _triangulation;

core/include/interface/Interface_Geometry.h

@@ -22,8 +22,8 @@ DLLEXPORT void Geometry_Get_N_Cells(State *state, int * n_cells, int idx_image=-
 // Get translation vectors ta, tb, tc
 DLLEXPORT void Geometry_Get_Translation_Vectors(State *state, float * ta, float * tb, float * tc, int idx_image=-1, int idx_chain=-1);
 
-// Find out if the system is a true 2D system
-DLLEXPORT bool Geometry_Is_2D(State * state, int idx_image=-1, int idx_chain=-1);
+// Retrieve dimensionality of the system (0, 1, 2, 3)
+DLLEXPORT int Geometry_Get_Dimensionality(State * state, int idx_image=-1, int idx_chain=-1);
 
 // Get the 3D Delaunay triangulation. Returns the number of tetrahedrons and
 // sets *indices_ptr to point to a list of index 4-tuples.

core/src/data/Geometry.cpp

@@ -1,4 +1,10 @@
 #include <Geometry.hpp>
+#include <Neighbours.hpp>
+#include <Exception.hpp>
+
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
 #include <glm/vec3.hpp>
 #include <glm/glm.hpp>
 #include "Qhull.h"
@@ -33,6 +39,9 @@ namespace Data
 		{
 			this->center[dim] = (this->bounds_min[dim] + this->bounds_max[dim]) / 2.0;
 		}
+
+		// Calculate dimensionality
+		this->dimensionality = calculateDimensionality();
 	}
 
     std::vector<tetrahedron_t> compute_delaunay_triangulation(const std::vector<vector_t> &points) {
@@ -57,61 +66,188 @@ namespace Data
         return tetrahedra;
     }
 
-  const std::vector<tetrahedron_t>& Geometry::triangulation() {
-    if (is2D()) {
-      _triangulation.clear();
-      return _triangulation;
-    }
-    if (_triangulation.size() == 0) {
-      bool is_simple_regular_geometry = true;
-      if (is_simple_regular_geometry) {
-        _triangulation.clear();
-        int cell_indices[] = {
-          0, 1, 5, 3,
-          1, 3, 2, 5,
-          3, 2, 5, 6,
-          7, 6, 5, 3,
-          4, 7, 5, 3,
-          0, 4, 3, 5
-        };
-        int x_offset = 1;
-        int y_offset = n_cells[0];
-        int z_offset = n_cells[0]*n_cells[1];
-        int offsets[] = {
-          0, x_offset, x_offset+y_offset, y_offset,
-          z_offset, x_offset+z_offset, x_offset+y_offset+z_offset, y_offset+z_offset
-        };
+	const std::vector<tetrahedron_t>& Geometry::triangulation()
+	{
+		if (dimensionality == 2)
+		{
+			_triangulation.clear();
+			return _triangulation;
+		}
+		if (_triangulation.size() == 0)
+		{
+			bool is_simple_regular_geometry = true;
+			if (is_simple_regular_geometry)
+			{
+				_triangulation.clear();
+				int cell_indices[] = {
+					0, 1, 5, 3,
+					1, 3, 2, 5,
+					3, 2, 5, 6,
+					7, 6, 5, 3,
+					4, 7, 5, 3,
+					0, 4, 3, 5
+					};
+				int x_offset = 1;
+				int y_offset = n_cells[0];
+				int z_offset = n_cells[0]*n_cells[1];
+				int offsets[] = {
+					0, x_offset, x_offset+y_offset, y_offset,
+					z_offset, x_offset+z_offset, x_offset+y_offset+z_offset, y_offset+z_offset
+					};
 
-        for (int ix = 0; ix < n_cells[0]-1; ix++) {
-          for (int iy = 0; iy < n_cells[1]-1; iy++) {
-            for (int iz = 0; iz < n_cells[2]-1; iz++) {
-              int base_index = ix*x_offset+iy*y_offset+iz*z_offset;
-              for (int j = 0; j < 6; j++) {
-                tetrahedron_t tetrahedron;
-                for (int k = 0; k < 4; k++) {
-                  int index = base_index + offsets[cell_indices[j*4+k]];
-                  tetrahedron[k] = index;
-                }
-                _triangulation.push_back(tetrahedron);
-              }
-            }
-          }
-        }
-      } else {
-        std::vector<vector_t> points;
-        points.resize(spin_pos.size()/3);
-        for (std::vector<vector_t>::size_type i = 0; i < points.size(); i++) {
-          points[i].x = spin_pos[i];
-          points[i].y = spin_pos[points.size()+i];
-          points[i].z = spin_pos[points.size()*2+i];
-        }
-        _triangulation = compute_delaunay_triangulation(points);
-      }
-    }
-    return _triangulation;
-  }
-  bool Geometry::is2D() const {
-    return (n_spins_basic_domain == 1) && (n_cells[0] == 1 || n_cells[1] == 1 || n_cells[2] == 1);
-  }
+				for (int ix = 0; ix < n_cells[0]-1; ix++)
+				{
+					for (int iy = 0; iy < n_cells[1]-1; iy++)
+					{
+						for (int iz = 0; iz < n_cells[2]-1; iz++)
+						{
+							int base_index = ix*x_offset+iy*y_offset+iz*z_offset;
+							for (int j = 0; j < 6; j++)
+							{
+								tetrahedron_t tetrahedron;
+								for (int k = 0; k < 4; k++)
+								{
+									int index = base_index + offsets[cell_indices[j*4+k]];
+									tetrahedron[k] = index;
+								}
+								_triangulation.push_back(tetrahedron);
+							}
+						}
+					}
+				}
+			}
+			else 
+			{
+				std::vector<vector_t> points;
+				points.resize(spin_pos.size()/3);
+				for (std::vector<vector_t>::size_type i = 0; i < points.size(); i++)
+				{
+					points[i].x = spin_pos[i];
+					points[i].y = spin_pos[points.size()+i];
+					points[i].z = spin_pos[points.size()*2+i];
+				}
+				_triangulation = compute_delaunay_triangulation(points);
+			}
+		}
+	return _triangulation;
+	}
+
+	int Geometry::calculateDimensionality() const
+	{
+		int dims_basis = 0, dims_translations = 0;
+		Eigen::Vector3d test_vec_basis, test_vec_translations;
+
+		// ----- Find dimensionality of the basis -----
+		if (n_spins_basic_domain == 1) dims_basis = 0;
+		else if (n_spins_basic_domain == 2) dims_basis = 1;
+		else if (n_spins_basic_domain == 3) dims_basis = 2;
+		else
+		{
+			// Get basis atoms relative to the first atom
+			Eigen::Vector3d v0 = { basis_atoms[0][0], basis_atoms[1][0], basis_atoms[2][0] };
+			std::vector<Eigen::Vector3d> b_vectors(n_spins_basic_domain-1);
+			for (int i = 1; i < n_spins_basic_domain; ++i)
+			{
+				b_vectors[i-1] = Eigen::Vector3d{ basis_atoms[0][i], basis_atoms[1][i], basis_atoms[2][i] } - v0;
+			}
+			// Calculate basis dimensionality
+			// test vec is along line
+			test_vec_basis = b_vectors[0];
+			//		is it 1D?
+			int n_parallel = 0;
+			for (unsigned int i = 1; i < b_vectors.size(); ++i)
+			{
+				if (std::abs(b_vectors[i].dot(test_vec_basis) - 1.0) < 1e-9) ++n_parallel;
+				// else n_parallel will give us the last parallel vector
+				// also the if-statement for dims_basis=1 wont be met
+				else break;
+			}
+			if (n_parallel == b_vectors.size() - 1)
+			{
+				dims_basis = 1;
+			}
+			else
+			{
+				// test vec is normal to plane
+				test_vec_basis = b_vectors[0].cross(b_vectors[n_parallel+1]);
+				//		is it 2D?
+				int n_in_plane = 0;
+				for (unsigned int i = 2; i < b_vectors.size(); ++i)
+				{
+					if (std::abs(b_vectors[i].dot(test_vec_basis)) < 1e-9) ++n_in_plane;
+				}
+				if (n_in_plane == b_vectors.size() - 2)
+				{
+					dims_basis = 2;
+				}
+				else return 3;
+			}
+		}
+
+
+		// ----- Find dimensionality of the translations -----
+		std::vector<Eigen::Vector3d> t_vectors(3);
+		for (int i = 0; i < 3; ++i)
+		{
+			t_vectors[i] = { translation_vectors[0][i], translation_vectors[1][i], translation_vectors[2][i] };
+		}
+		//		The following are zero if the corresponding pair is parallel
+		double t01, t02, t12;
+		t01 = std::abs(t_vectors[0].dot(t_vectors[1]) - 1.0);
+		t02 = std::abs(t_vectors[0].dot(t_vectors[2]) - 1.0);
+		t12 = std::abs(t_vectors[1].dot(t_vectors[2]) - 1.0);
+		//		Check if pairs are linearly independent
+		int n_independent_pairs = 0;
+		if (t01>1e-9 && n_cells[0] > 1 && n_cells[1] > 1) ++n_independent_pairs;
+		if (t02>1e-9 && n_cells[0] > 1 && n_cells[2] > 1) ++n_independent_pairs;
+		if (t12>1e-9 && n_cells[1] > 1 && n_cells[2] > 1) ++n_independent_pairs;
+		//		Calculate translations dimensionality
+		if (n_cells[0] == 1 && n_cells[1] == 1 && n_cells[2] == 1) dims_translations = 0;
+		else if (n_independent_pairs == 0)
+		{
+			dims_translations = 1;
+			// test vec is along the line
+			for (int i=0; i<3; ++i) if (n_cells[i] > 1) test_vec_translations = t_vectors[i];
+		}
+		else if (n_independent_pairs < 3)
+		{
+			dims_translations = 2;
+			// test vec is normal to plane
+			int n = 0;
+			std::vector<Eigen::Vector3d> plane(2);
+			for (int i = 0; i < 3; ++i)
+			{
+				if (n_cells[i] > 1) plane[n] = t_vectors[i];
+				++n;
+			}
+			test_vec_translations = plane[0].cross(plane[1]);
+		}
+		else return 3;
+
+
+		// ----- Calculate dimensionality of system -----
+		test_vec_basis.normalize();
+		test_vec_translations.normalize();
+		//		If one dimensionality is zero, only the other counts
+		if (dims_basis == 0) return dims_translations;
+		else if (dims_translations == 0) return dims_basis;
+		//		If both are linear or both are planar, the test vectors should be parallel if the geometry is 1D or 2D
+		else if (dims_basis == dims_translations)
+		{
+			if (std::abs(test_vec_basis.dot(test_vec_translations) - 1.0) < 1e-9) return dims_basis;
+			else if (dims_basis == 1) return 2;
+			else if (dims_basis == 2) return 3;
+		}
+		//		If one is linear (1D), and the other planar (2D) then the test vectors should be orthogonal if the geometry is 2D
+		else if ( (dims_basis == 1 && dims_translations == 2) || (dims_basis == 2 && dims_translations == 1) )
+		{
+			if (std::abs(test_vec_basis.dot(test_vec_translations)) < 1e-9) return 2;
+			else return 3;
+		}
+
+		// We should never get here
+		throw Utility::Exception::Unknown_Exception;
+		return 0;
+	}
 }
 

core/src/interface/Interface_Geometry.cpp

@@ -91,14 +91,14 @@ void Geometry_Get_Translation_Vectors(State *state, float * ta, float * tb, floa
     }
 }
 
-bool Geometry_Is_2D(State * state, int idx_image, int idx_chain)
+int Geometry_Get_Dimensionality(State * state, int idx_image, int idx_chain)
 {
 	std::shared_ptr<Data::Spin_System> image;
 	std::shared_ptr<Data::Spin_System_Chain> chain;
 	from_indices(state, idx_image, idx_chain, image, chain);
 
-  auto g = image->geometry;
-  return g->is2D();
+	auto g = image->geometry;
+	return g->dimensionality;
 }
 
 

core/src/utility/IO_Configparser.cpp

@@ -272,6 +272,7 @@ namespace Utility
 
 			// Return geometry
 			auto geometry = std::unique_ptr<Data::Geometry>(new Data::Geometry(basis, translation_vectors, n_cells, no_spins_basic_domain, basis_atoms, spin_pos));
+			Log(Log_Level::Parameter, Log_Sender::IO, "Geometry dims: " + std::to_string(geometry->dimensionality) + "D"); 
 			Log(Log_Level::Info, Log_Sender::IO, "Geometry: built");
 			return geometry;
 		}// end Geometry from Config
@@ -623,7 +624,7 @@ namespace Utility
 			bool anisotropy_from_file = false;
 			std::vector<int> anisotropy_index(geometry.nos);				// [nos]
 			std::vector<double> anisotropy_magnitude(geometry.nos, 0.0);	// [nos]
-			std::vector<std::vector<double>> anisotropy_normal(geometry.nos, { 0.0, 0.0, 1.0 });	// [nos][3]
+			std::vector<std::vector<double>> anisotropy_normal(geometry.nos, K_normal);	// [nos][3]
 
 			// ------------ Two Spin Interactions ------------
 			int n_pairs = 0;

core/src/utility/Vectormath.cpp

@@ -108,9 +108,9 @@ namespace Utility
 
 			// Check for erronous input placing two spins on the same location
 			std::vector<double> sp(3);
-			for (int i = 0; i < basis_atoms[0].size(); ++i)
+			for (unsigned int i = 0; i < basis_atoms[0].size(); ++i)
 			{
-				for (int j = 0; j < basis_atoms[0].size(); ++j)
+				for (unsigned int j = 0; j < basis_atoms[0].size(); ++j)
 				{
 					for (int k1 = -2; k1 < 3; ++k1)
 					{

gl/src/SphereSpinRenderer.cpp

@@ -67,10 +67,11 @@ void SphereSpinRenderer::optionsHaveChanged(const std::vector<int>& changedOptio
 void SphereSpinRenderer::updateSpins(const std::vector<glm::vec3>& positions,
                                      const std::vector<glm::vec3>& directions) {
   CHECK_GL_ERROR;
-  glBindVertexArray(_vao1);
+  // Binding _vao1 here seems optional, binding it throws a gl error on some systems... idk why
+  glBindVertexArray(0);
+  //glBindVertexArray(_vao1);
   glBindBuffer(GL_ARRAY_BUFFER, _instanceDirectionVbo);
   glBufferData(GL_ARRAY_BUFFER, sizeof(glm::vec3) * directions.size(), directions.data(), GL_STREAM_DRAW);
-  glBufferData(GL_ARRAY_BUFFER, sizeof(glm::vec3) * directions.size(), directions.data(), GL_STREAM_DRAW);
   _numInstances = directions.size();
   CHECK_GL_ERROR;
 }

ui-qt/include/ControlWidget.hpp

@@ -25,7 +25,7 @@ class ControlWidget : public QWidget, private Ui::ControlWidget
 
 public:
 	ControlWidget(std::shared_ptr<State> state, SpinWidget *spinWidget);
-	void update();
+	void updateData();
 
 public slots:
 	void play_pause();

ui-qt/include/DebugWidget.hpp

@@ -16,7 +16,7 @@ class DebugWidget : public QWidget, private Ui::DebugWidget
 
 public:
 	DebugWidget(std::shared_ptr<State> state);
-	void update();
+	void updateData();
 
 	void LoadFromLog();
 	void UpdateFromLog();

ui-qt/include/PlotWidget.hpp

@@ -15,7 +15,7 @@ class PlotWidget : public QtCharts::QChartView	// We need a proper 2D plotting s
 
 public:
 	PlotWidget(std::shared_ptr<State> state, bool plot_interpolated=false);
-	void update();
+	void updateData();
 	void set_interpolating(bool b);
 
 private:

ui-qt/include/SettingsWidget.hpp

@@ -17,7 +17,7 @@ class SettingsWidget : public QWidget, private Ui::SettingsWidget
 
 public:
 	SettingsWidget(std::shared_ptr<State> state, SpinWidget *spinWidget);
-	void update();
+	void updateData();
 	void SelectTab(int index);
 
 	std::shared_ptr<State> state;