Core: added Vectormath function to re-distribute a vector-field on ne…

…w dimensions.

The Geometry API now uses this to fit the spins and effective field arrays to an updated geometry.
Note: a shift can be applied in order to move (or re-center) the system, e.g. when the dimensions are increased.

core/include/engine/Vectormath.hpp

@@ -28,6 +28,7 @@ namespace Engine
         /////////////////////////////////////////////////////////////////
         //////// Translating across the lattice
 
+        // Note: translations must lie within bounds of n_cells
         inline int idx_from_translations(const intfield & n_cells, const int n_cell_atoms, const std::array<int, 3> & translations)
         {
             int Na = n_cells[0];
@@ -397,6 +398,12 @@ namespace Engine
         //      This requires to know the underlying geometry, as well as the boundary conditions.
         void directional_gradient(const vectorfield & vf, const Data::Geometry & geometry, const intfield & boundary_conditions, const Vector3 & direction, vectorfield & gradient);
 
+        /////////////////////////////////////////////////////////////////
+
+        // Re-distribute a given vector-field according to a new set of dimensions.
+        vectorfield change_dimensions(vectorfield & sf, int n_cell_atoms, intfield n_cells,
+            intfield dimensions_new, std::array<int,3> shift=std::array<int,3>{0,0,0});
+
         /////////////////////////////////////////////////////////////////
         //////// Vectormath-like operations
 

core/src/Spirit/Geometry.cpp

@@ -10,18 +10,20 @@
 
 void Helper_System_Set_Geometry(std::shared_ptr<Data::Spin_System> system, const Data::Geometry & new_geometry)
 {
-    *system->geometry = new_geometry;
+    // *system->geometry = new_geometry;
     auto ge = system->geometry;
 
     // Spins
     int nos_old = system->nos;
-    int nos = ge->nos;
+    int nos = new_geometry.nos;
     system->nos = nos;
-    // TODO: ordering of spins should be considered and date potentially extrapolated -> write a function for this
-    system->spins->resize(nos);
-    system->effective_field.resize(nos);
-    for (int i = nos_old; i<nos; ++i) (*system->spins)[i] = Vector3{ 0, 0, 1 };
-    for (int i = nos_old; i<nos; ++i) system->effective_field[i] = Vector3{ 0, 0, 1 };
+
+    // Move the vector-fields to the new geometry
+    *system->spins = Engine::Vectormath::change_dimensions(*system->spins, ge->n_cell_atoms, ge->n_cells, new_geometry.n_cells);
+    system->effective_field = Engine::Vectormath::change_dimensions(system->effective_field, ge->n_cell_atoms, ge->n_cells, new_geometry.n_cells);
+
+    // Update the system geometry
+    *system->geometry = new_geometry;
 
     // Parameters
     // TODO: properly re-generate pinning

core/src/engine/Vectormath.cpp

@@ -370,6 +370,39 @@ namespace Engine
 
         /////////////////////////////////////////////////////////////////
 
+        vectorfield change_dimensions(vectorfield & sf, int n_cell_atoms, intfield n_cells,
+            intfield dimensions_new, std::array<int,3> shift)
+        {
+            int N_old = n_cell_atoms*n_cells[0]*n_cells[1]*dimensions_new[2];
+            int N_new = n_cell_atoms*dimensions_new[0]*dimensions_new[1]*dimensions_new[2];
+            vectorfield newfield(N_new);
+
+            #pragma omp parallel for collapse(3)
+            for (int i=0; i<dimensions_new[0]; ++i)
+            {
+                for (int j=0; j<dimensions_new[1]; ++j)
+                {
+                    for (int k=0; k<dimensions_new[2]; ++k)
+                    {
+                        for (int iatom=0; iatom<n_cell_atoms; ++iatom)
+                        {
+                            int idx_old = iatom + idx_from_translations(n_cells, n_cell_atoms, {i,j,k});
+
+                            int idx_new = iatom + idx_from_translations(dimensions_new, n_cell_atoms, {i,j,k}, shift);
+
+                            if ( (i>=n_cells[0]) || (j>=n_cells[1]) || (k>=n_cells[2]))
+                                newfield[idx_new] = {0,0,1};
+                            else
+                                newfield[idx_new] = sf[idx_old];
+                        }
+                    }
+                }
+            }
+            return newfield;
+        }
+
+        /////////////////////////////////////////////////////////////////
+
         void fill(scalarfield & sf, scalar s)
         {
             #pragma omp parallel for

core/src/engine/Vectormath.cu

@@ -621,6 +621,38 @@ namespace Engine
 
 
 
+        /////////////////////////////////////////////////////////////////
+
+        vectorfield change_dimensions(vectorfield & sf, int n_cell_atoms, intfield n_cells,
+            intfield dimensions_new, std::array<int,3> shift)
+        {
+            int N_old = n_cell_atoms*n_cells[0]*n_cells[1]*dimensions_new[2];
+            int N_new = n_cell_atoms*dimensions_new[0]*dimensions_new[1]*dimensions_new[2];
+            vectorfield newfield(N_new);
+
+            for (int i=0; i<dimensions_new[0]; ++i)
+            {
+                for (int j=0; j<dimensions_new[1]; ++j)
+                {
+                    for (int k=0; k<dimensions_new[2]; ++k)
+                    {
+                        for (int iatom=0; iatom<n_cell_atoms; ++iatom)
+                        {
+                            int idx_old = iatom + idx_from_translations(n_cells, n_cell_atoms, {i,j,k});
+
+                            int idx_new = iatom + idx_from_translations(dimensions_new, n_cell_atoms, {i,j,k}, shift.data());
+
+                            if ( (i>=n_cells[0]) || (j>=n_cells[1]) || (k>=n_cells[2]))
+                                newfield[idx_new] = {0,0,1};
+                            else
+                                newfield[idx_new] = sf[idx_old];
+                        }
+                    }
+                }
+            }
+            return newfield;
+        }
+
         /////////////////////////////////////////////////////////////////