Fixing bug in matrix transpose

src/CMakeLists.txt

@@ -28,7 +28,7 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake/modules )
 # prepare configuration file
 SET(VERSION_MAJOR "2")
 SET(VERSION_MINOR "0")
-SET(VERSION_MICRO "1")
+SET(VERSION_MICRO "2")
 configure_file(${CMAKE_SOURCE_DIR}/config.h.in ${CMAKE_SOURCE_DIR}/config.h @ONLY)
 if(EXISTS "${CMAKE_SOURCE_DIR}/config.h")
     MESSAGE("[USER] Creating ${CMAKE_SOURCE_DIR}/config.h from ${CMAKE_SOURCE_DIR}/config.h.in")
@@ -87,27 +87,18 @@ include_directories(${CMAKE_CURRENT_SOURCE_DIR}
 
 # Add sources
 file(GLOB SOURCES "*.cpp")
-add_executable(edp ${SOURCES})
+list(REMOVE_ITEM SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/edp.cpp)
 add_library(edpsources STATIC ${SOURCES})
+add_executable(edp ${CMAKE_CURRENT_SOURCE_DIR}/edp.cpp)
 
-# Set C++14
-add_definitions(-std=c++14)
+# Set c++17 and use march=native
+add_definitions(-std=c++17 -march=native)
 
 # Link libraries
-if(UNIX AND NOT APPLE)
-    SET(CMAKE_EXE_LINKER_FLAGS "-Wl,-rpath=\$ORIGIN/lib")
-endif()
-if(APPLE)
-    SET(CMAKE_MACOSX_RPATH TRUE)
-    SET_TARGET_PROPERTIES(edp PROPERTIES INSTALL_RPATH "@executable_path/lib")
-    target_link_libraries(edp ${Boost_LIBRARIES} ${CAIRO_LIBRARIES})
-else()
-    target_link_libraries(edp ${Boost_LIBRARIES} ${CAIRO_LIBRARIES})
-endif()
-
-# add Wno-literal-suffix to suppress warning messages
-set(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS}")
-
+target_link_libraries(edp
+                      edpsources
+                      ${Boost_LIBRARIES}
+                      ${CAIRO_LIBRARIES})
 
 if(NOT DISABLE_TEST)
     message("[USER] Performing unit tests")

src/edp.cpp

@@ -174,7 +174,7 @@ int main(int argc, char *argv[]) {
             v[2] = boost::lexical_cast<double>(what[3]);
         } else if(boost::regex_match(v_str, what, re_scalar_2)) {
             v = ( sf.get_atom_position(boost::lexical_cast<unsigned int>(what[2])-1) -
-                  sf.get_atom_position(boost::lexical_cast<unsigned int>(what[1])-1)).normalized();
+                  sf.get_atom_position(boost::lexical_cast<unsigned int>(what[1])-1));
 
             // if the user has supplied a "xyz" directives, check which of these
             // vector components should be retained
@@ -195,6 +195,8 @@ int main(int argc, char *argv[]) {
                 }
             }
 
+            v = v.normalized();
+
         } else {
             std::runtime_error("Could not obtain a vector v");
         }

src/scalar_field.cpp

@@ -561,16 +561,13 @@ fpt ScalarField::get_value(unsigned int i, unsigned int j, unsigned int k) const
  *
  */
 Vec3 ScalarField::grid_to_realspace(fpt i, fpt j, fpt k) const {
+    // build direct coordinates
     fpt dx = (fpt)i / (fpt)grid_dimensions[0];
     fpt dy = (fpt)j / (fpt)grid_dimensions[1];
     fpt dz = (fpt)k / (fpt)grid_dimensions[2];
+    Vec3 d(dx, dy, dz);
 
-    Vec3 r;
-    r[0] = this->mat(0,0) * dx + this->mat(1,0) * dy + this->mat(2,0) * dz;
-    r[1] = this->mat(0,1) * dx + this->mat(1,1) * dy + this->mat(2,1) * dz;
-    r[2] = this->mat(0,2) * dx + this->mat(1,2) * dy + this->mat(2,2) * dz;
-
-    return r;
+    return this->mat.transpose() * d;
 }
 
 /*
@@ -598,13 +595,10 @@ Vec3 ScalarField::realspace_to_grid(fpt i, fpt j, fpt k) const {
  * Convert 3d realspace vector to direct position.
  *
  */
-Vec3 ScalarField::realspace_to_direct(fpt i, fpt j, fpt k) const {
-    Vec3 r;
-    r[0] = this->imat(0,0) * i + this->imat(0,1) * j + this->imat(0,2) * k;
-    r[1] = this->imat(1,0) * i + this->imat(1,1) * j + this->imat(1,2) * k;
-    r[2] = this->imat(2,0) * i + this->imat(2,1) * j + this->imat(2,2) * k;
+Vec3 ScalarField::realspace_to_direct(fpt x, fpt y, fpt z) const {
+    Vec3 r(x,y,z);
 
-    return r;
+    return this->imat.transpose() * r;
 }
 
 fpt ScalarField::get_max() const {
@@ -617,7 +611,7 @@ fpt ScalarField::get_min() const {
 
 Vec3 ScalarField::get_atom_position(unsigned int atid) const {
     if(atid < this->atom_pos.size()) {
-        return this->mat * this->atom_pos[atid];
+        return this->mat.transpose() * this->atom_pos[atid];
     } else {
         throw std::runtime_error("Requested atom id lies outside bounds");
     }

src/scalar_field.h

@@ -154,7 +154,7 @@ class ScalarField{
 
     Vec3 grid_to_realspace(fpt i, fpt j, fpt k) const;
 
-    Vec3 realspace_to_grid(fpt i, fpt j, fpt k) const;
+    Vec3 realspace_to_grid(fpt x, fpt y, fpt z) const;
 
     Vec3 realspace_to_direct(fpt i, fpt j, fpt k) const;