Multidimensional basis and other improvements

.gitignore

@@ -18,6 +18,7 @@ install_manifest.txt
 compile_commands.json
 CTestTestfile.cmake
 /.vs/*
+/.vscode/*
 
 *.md
 LICENSE

CMakeLists.txt

@@ -2,6 +2,7 @@ cmake_minimum_required(VERSION 3.5)
 project(Schroedinger)
 
 # Force C++17 standard
+set(CMAKE_BUILD_TYPE             Debug)
 set(CMAKE_CXX_STANDARD              17)
 set(CMAKE_CXX_STANDARD_REQUIRED     ON)
 set(CMAKE_CXX_FLAGS -D_EMULATE_GLIBC=0)

src/Basis/Base.cpp

@@ -35,3 +35,24 @@ std::vector<ContinuousBase> Base::getContinuous() {
 std::vector<DiscreteBase> Base::getDiscrete() {
 	return this->discrete;
 }
+
+std::ostream& operator<<(std::ostream& stream, Base& base) {
+
+	// Print continuous dimension values (if present)
+	if (base.getContinuous().size() > 0)
+		for (int i = 0; i < base.getContinuous().size(); i++) {
+			for (int coord_counter = 0; coord_counter < base.getContinuous().at(i).getCoords().size(); coord_counter++)
+				stream << base.getContinuous().at(i).getCoords().at(coord_counter) << "; ";
+		}
+
+	stream << "\n\n";
+
+	// Print discrete dimension values (if present)
+	if (base.getDiscrete().size() > 0)
+		for (int i = 0; i < base.getContinuous().size(); i++) {
+			for (int coord_counter = 0; coord_counter < base.getDiscrete().at(i).getCoords().size(); coord_counter++)
+				stream << base.getDiscrete().at(i).getCoords().at(coord_counter) << "; ";
+		}
+
+    return stream;
+ }

src/Basis/Base.h

@@ -19,6 +19,7 @@ class Base
 	int getDim();
 	std::vector<ContinuousBase> getContinuous();
 	std::vector<DiscreteBase> getDiscrete();
+	friend std::ostream& operator<< (std::ostream& stream, Base& base);
 
 private:
 	std::vector< DiscreteBase > discrete;

src/Basis/BasisManager.cpp

@@ -21,8 +21,8 @@ void BasisManager::addBase(Base b) {
 	this-> bases.push_back(b);
 
 	// If it's there's just this one in the vector, then it's automatically selected
-		if (this->bases.size() == 1)
-		this->selectBase(b);
+	if (this->bases.size() == 1)
+	this->selectBase(b);
 }
 
 std::vector<Base> BasisManager::getBasisList() {
@@ -42,39 +42,39 @@ std::vector<Base> BasisManager::getBasisList(Source s) {
 	return this->bases;
 }
 
-Base BasisManager::Builder::build() {
+Base BasisManager::Builder::build(int dimension) {
 	//TODO: Eventually add controls...
-    int dimension = 0;
-
 	return Base(Base::basePreset::Custom, dimension, c_base, d_base);
 }
 
 Base BasisManager::Builder::build(Base::basePreset b, int dimension) {
     //TODO: Eventually add controls...
-        return Base(b, dimension, c_base, d_base);
-    }
+    return Base(b, dimension, c_base, d_base);
+}
 
 // --- Factory Methods --- //
 Base BasisManager::Builder::build(Base::basePreset b, int dimension, double mesh, int nbox) {
     //TODO: Eventually add controls...
 
     switch(b){
-
-        case Base::basePreset::Custom: throw std::invalid_argument("Custom basis not meaningful with parameters!");
-            break;
         case Base::basePreset::Cartesian:
             std::cout<< "Building Cartesian basis in " << dimension << "dimensions, with nbox = " << nbox << ", mesh = " << mesh << std::endl;
-            for(int i = 0; i < dimension; i++)
-            {
+            for(int i = 0; i < dimension; i++) 
                 this->addContinuous(mesh, nbox);
-            }
             break;
+
         // todo: find criterium criterium to define Lmax given only r. Harmonic oscillator r/l relation?
-        case Base::basePreset::Cylindrical: std::invalid_argument("Wrong parameters for Cylindrical basis");
+        case Base::basePreset::Cylindrical: 
+			throw std::invalid_argument("Wrong parameters for Cylindrical basis");
+            break;
+        case Base::basePreset::Spherical: 
+			throw std::invalid_argument("Wrong parameters for Spherical basis");
             break;
-        case Base::basePreset::Spherical: std::invalid_argument("Wrong parameters for Spherical basis");
+        case Base::basePreset::Custom: 
+			throw std::invalid_argument("Custom basis not meaningful with parameters!");
             break;
-        default: throw std::invalid_argument("Wrong basis type or initialization meaningless!");
+        default: 
+			throw std::invalid_argument("Wrong basis type or initialization meaningless!");
             break;
 
     }
@@ -83,9 +83,12 @@ Base BasisManager::Builder::build(Base::basePreset b, int dimension, double mesh
 }
 
 Base BasisManager::Builder::build(SphericalInitializer ini) {
-    if(ini.Lmax == 0) std::invalid_argument("Spherical basis with Lmax = 0 does not have sense");
-    if(ini.mesh <= 0 ) std::invalid_argument("mesh < 0 does not have sense");
-    if(ini.end <= ini.mesh ) std::invalid_argument("Spherical basis with r < mesh does not have sense");
+    if(ini.Lmax == 0) 
+		throw std::invalid_argument("Spherical basis with Lmax = 0 does not have sense");
+    if(ini.mesh <= 0 ) 
+		throw std::invalid_argument("mesh < 0 does not have sense");
+    if(ini.end <= ini.mesh ) 
+		throw std::invalid_argument("Spherical basis with r < mesh does not have sense");
 
     this->addContinuous(ini.start,ini.end,ini.mesh);
     this->addDiscrete(ini.Lmin,ini.Lmax,ini.Lstep);
@@ -94,16 +97,17 @@ Base BasisManager::Builder::build(SphericalInitializer ini) {
 }
 
 Base BasisManager::Builder::build(ContinuousInitializer ini) {
-	if(ini.mesh <= 0 ) std::invalid_argument("mesh < 0 does not have sense");
-	if(ini.end <= ini.mesh ) std::invalid_argument("xmax < mesh does not have sense");
 
-	std::cout<< "Building Cartesian basis between " << ini.start << ", " << ini.end << "  mesh = " << ini.mesh << std::endl;
+	if(ini.mesh <= 0 ) 
+		throw std::invalid_argument("mesh < 0 does not have sense");
+	if(ini.end <= ini.mesh ) 
+		throw std::invalid_argument("xmax < mesh does not have sense");
 
 	this->addContinuous(ini.start,ini.end,ini.mesh);
 
+	std::cout<< "Building Cartesian basis between " << ini.start << ", " << ini.end << "  mesh = " << ini.mesh << std::endl;
 	return Base(Base::basePreset::Spherical, 3, c_base, d_base);
 }
-
 // --- End Factory --- //
 
 

src/Basis/BasisManager.h

@@ -1,3 +1,6 @@
+#ifndef BASISMANAGER_H
+#define BASISMANAGER_H
+
 #include <Base.h>
 #include <Initializer.h>
 
@@ -18,7 +21,7 @@ class BasisManager
 		std::vector< DiscreteBase > d_base;
 		std::vector< ContinuousBase > c_base;
 	public:
-        Base build();
+        Base build(int dimension);
         Base build(ContinuousInitializer);
         Base build(SphericalInitializer);
         Base build(Base::basePreset, int dimension);
@@ -38,3 +41,5 @@ class BasisManager
 	static BasisManager* instance;
 	BasisManager() {}
 };
+
+#endif

src/Basis/ContinuousBase.cpp

@@ -20,10 +20,10 @@ ContinuousBase::ContinuousBase(double start, double end, double mesh)
 		std::invalid_argument("CountinousBase starting-end = 0");
 	}
 
-	this->start = start;
-	this->end = end;
-	this->mesh = mesh;
-	this->nbox = (unsigned int)((end - start) / mesh);
+	this->start  = start;
+	this->end    = end;
+	this->mesh   = mesh;
+	this->nbox   = (unsigned int)((end - start) / mesh);
 	this->coords = evaluate();
 }
 
@@ -33,10 +33,10 @@ ContinuousBase::ContinuousBase(double start, double end, unsigned int nbox)
 		std::invalid_argument("CountinousBase starting-end = 0");
 	}
 
-	this->start = start;
-	this->end = end;
-	this->mesh = (end - start) / nbox;
-	this->nbox = nbox;
+	this->start  = start;
+	this->end    = end;
+	this->mesh   = (end - start) / nbox;
+	this->nbox   = nbox;
 	this->coords = evaluate();
 }
 

src/Potential/Potential.cpp

@@ -1,46 +1,129 @@
 #include "Potential.h"
 
-Potential::Potential(std::vector<double> coord, std::string type, double k, double width, double height)
+Potential::Potential() {
+    throw std::invalid_argument("Can't initialize a potential without specifying parameters");
+}
+Potential::Potential(Base base, PotentialType type, double k, double width, double height, bool separable)
 {
-    this->x        = coord;
-    this->v        = std::vector<double>(x.size());
-    this->k        = k;
-    this->width    = width;
-    this->height   = height;
-    this->type     = type;
+    this->base      = base; 
+    this->k         = k;
+    this->width     = width;
+    this->height    = height;
+    this->type      = type;
+    this->separable = separable;
+    this->v         = std::vector<double>();
 
-    if(type.compare("box potential") == 0 || type.compare("box") == 0 || type.compare("0") == 0)
-        this->box_potential();
+    if (!this->separable) {
+        switch(type) {
+            case BOX_POTENTIAL:
+                this->box_potential();
+                break;
+            case HARMONIC_OSCILLATOR:
+                this->ho_potential();
+                break;
+            case FINITE_WELL_POTENTIAL:
+                this->finite_well_potential();
+                break;
+            default:
+                throw std::invalid_argument("Wrong potential type or initialization meaningless!");
+        }
+    }
+    // If the potential is separable, then n potentials (one for each base's representation)
+    else {
+        this->separated_potentials = std::vector<Potential>();
 
-    else if(type.compare("harmonic oscillator") == 0 || type.compare("ho") == 0 ||  type.compare("1") == 0)
-        this->ho_potential();
+        for (int i = 0; i < base.getContinuous().size(); i++) {
+            // Extract continuous representation from the main base
+            std::vector<ContinuousBase> c_base = std::vector<ContinuousBase>(1);
+            c_base.push_back(base.getContinuous().at(i));
 
-    else if(type.compare("finite well potential") == 0 || type.compare("well") == 0 || type.compare("2") == 0 )
-        this->finite_well_potential();
+            // Create new base object with the extracted representation
+            Base monodimensional = Base(Base::basePreset::Custom, 1, c_base, std::vector<DiscreteBase>());
 
-    else {
-        throw std::invalid_argument("Wrong potential type ("+type+") or initialization meaningless!");
+            // Create new potential with the new base
+            Potential separated_potential = Potential(monodimensional, this->type, this->k, this->width, this->height, false);
+
+            // Add the new Potential to the separated potenial vector associated to the main potential
+            this->separated_potentials.push_back(separated_potential);
+        }
+
+        for (int i = 0; i < base.getDiscrete().size(); i++) {
+            // Extract discrete representation from the main base
+            std::vector<DiscreteBase> d_base = std::vector<DiscreteBase>(1);
+            d_base.push_back(base.getDiscrete().at(i));
+
+            // Create new base object with the extracted representation
+            Base monodimensional = Base(Base::basePreset::Custom, 1, std::vector<ContinuousBase>(), d_base);
+
+            // Create new potential with the new base
+            Potential separated_potential = Potential(monodimensional, this->type, this->k, this->width, this->height, false);
+
+            // Add the new Potential to the separated potenial vector associated to the main potential
+            this->separated_potentials.push_back(separated_potential);
+        }
     }
 }
 
 void Potential::ho_potential()
 {
+    std::vector<double> x = this->getCoordsFromBase();
+    this->v.reserve(x.size());
     for(std::vector<int>::size_type i = 0; i < x.size(); i++)
-        this->v[i] = this->x[i] * this->x[i] * this->k;
+            this->v.push_back(x.at(i) * x.at(i) * this->k);
 }
 
 void Potential::box_potential()
 {
+    this->v = this->getCoordsFromBase();
     std::fill(this->v.begin(), this->v.end(), 0.0);
 }
 
 void Potential::finite_well_potential()
 {
-    for(std::vector<int>::size_type i = 0; i < x.size(); i++)
-        this->v[i] = (this->x[i] > -this->width/2.0 && this->x[i] < this->width/2.0) ? 0.0 : this->height;
+    std::vector<double> x = this->getCoordsFromBase();
+    this->v.reserve(x.size());
+    for(std::vector<int>::size_type i = 0; i < x.size(); i++) {
+        double value = (x.at(i) > -this->width/2.0 && x.at(i) < this->width/2.0) ? 0.0 : this->height;
+        this->v.push_back(value);
+    }
 }
 
 std::vector<double> Potential::getValues()
 {
+    this->printToFile();
     return this->v;
 }
+
+std::vector<double> Potential::getCoordsFromBase() 
+{
+    if (this->base.getContinuous().size() == 1)
+        return this->base.getContinuous().at(0).getCoords();
+
+    else if (this->base.getDiscrete().size() == 1) {
+        // tricky conversion taking each std::vector<int> value and returning a final std::vector<double> 
+        std::vector<int> original_coords = this->base.getDiscrete().at(0).getCoords();
+        return std::vector<double>(original_coords.begin(), original_coords.end());
+    }
+}
+
+bool Potential::isSeparated() {
+    return this->separable;
+}
+
+std::vector<Potential> Potential::getSeparatedPotentials() {
+    if (this->separable == true && this->separated_potentials.size() > 0)
+        return this->separated_potentials;
+}
+void Potential::printToFile() {
+  std::ofstream myfile ("potential.dat");
+  if (myfile.is_open())
+  {
+
+    for(int i = 0; i < this->v.size(); i ++){
+        myfile << i <<" " << this->v.at(i)<< std::endl ;
+    }
+    myfile.close();
+  }
+}
+
+