Page MenuHomec4science
Files F92358068

exercice6.cpp
No OneTemporary
Actions

File Metadata

Created
Tue, Nov 19, 16:19

exercice6.cpp
View Options

#include <iostream>
#include <string>
#include <cmath>
#include <vector>
#include <fstream>
#include "ConfigFile.tpp"
#include "debug.hpp"
#include "subdiag.h"
#include "discret.h"
#include "findiff.h"
using namespace std;
// constants, defined in main
class Symulation
{
// symulation and output setup
ofstream output;
bool trivial;
// system variables
interval rmap;
std::size_t N; // total N
subdiag<long double> A;
vector<long double> B;
const double eps0;
const double V0;
const double b;
const double R;
const double err_b;
const double a0;
public:
Symulation(ConfigFile& config)
: output(config.get<std::string>("output"), ios::out),
trivial(config.get<bool>("trivial")),
rmap(config.get<interval>("interval")),
N(rmap.N()),
A(N+1),
B(N+1, 0),
eps0(config.get<double>("eps0")),
V0(config.get<double>("V0")),
b(rmap.border(0)),
R(rmap.border(1)),
err_b(10e-12 * b),
a0(1e4)
{
// set output precision
output << std::setprecision(15);
// take values from config
double p = config.get<double>("trapezium");
// initialize interval and A matrix
for (std::size_t k = 0; k < N; ++k)
{
// compose A matrix
findiff::compose_matrix(A, k, rmap,
[&](double r) { // this is a lambda, Armando
return epsilonr(r) * r;
}, p);
// compose b vector, initial conditions
findiff::compose_rhs(B, k, rmap,
[&](double r) {
return r * rho(r) / eps0;
}, p);
}
A[N][N-1] = 0;
A[N][N] = 1;
B[N] = V0; // initial potenzial
npdebug("Value of R: ", R)
npdebug("Value of b: ", b)
/*for (std::size_t k = 0; k < N+1; ++k)
cout << "N = " << k << ", r = " << rmap[k] << endl;
cout << endl;*/
/*
// printing A
for (std::size_t k = 0; k < N+1; ++k)
{
for (std::size_t j = 0; j < N+1; ++j)
cout << A[k][j] << " ";
cout << endl;
}
cout << endl;
// printing B
for (std::size_t k = 0; k < N+1; ++k)
cout << B[k] << " ";
cout << endl;
cout << endl;
*/
}
~Symulation()
{
output.close();
}
inline bool discontinuity(double r) const
{
return (abs(r-b)<=err_b && r<=b) || abs(r-b) <= 5e-14;
}
//funzione che definisce la permeabilità elettrica relativa del dielettrico in funzione di r
inline double epsilonr(const double& r) const
{
if (trivial || r<=(b-err_b) || discontinuity(r))
return 1;
else
return 8.0 - 6.0 * (r-b)/(R-b);
}
//funzione che definisce la densità di carica in funzione di r
inline double rho(const double& r) const
{
if (trivial)
return eps0;
else if (r>b)
return 0;
else
return eps0*a0*(1.0 - pow(r/b,2));
}
// Helper functions
template<class T>
inline void write_arg(const T& value)
{
output << value << " ";
//cout << value << " ";
}
template<class ...Args>
inline void write(const Args& ...args)
{
(write_arg(args), ...);
output << endl;
//cout << endl;
}
// symulation execution
int exec()
{
// solve system
std::vector<long double> phi = A.solve(B);
// save data
//write("r", "b", "phi");
for (std::size_t i = 0; i < N+1; ++i)
write(rmap[i], B[i], phi[i]);
return 0;
}
};
int main(int argc, char* argv[])
{
string inputPath("configuration.in"); // Fichier d'input par defaut
if(argc>1) // Fichier d'input specifie par l'utilisateur ("./Exercice5 config_perso.in")
inputPath = argv[1];
ConfigFile configFile(inputPath); // Les parametres sont lus et stockes dans une "map" de strings.
for(int i(2); i<argc; ++i) // Input complementaires ("./Exercice5 config_perso.in input_scan=[valeur]")
configFile.process(argv[i]);
Symulation s(configFile);
return s.exec();
}

Event Timeline

c4science · Help