diff --git a/physnum/rap6/include/ConfigFile.h b/physnum/rap6/include/ConfigFile.h new file mode 100644 index 0000000..a7fa42b --- /dev/null +++ b/physnum/rap6/include/ConfigFile.h @@ -0,0 +1,37 @@ +// Classe facilitant la lecture de fichiers de configuration. +// Contributeurs : K. Steiner, J. Dominski, N. Ohana +// Utilisation : Envoyer au constructeur le nom d'un fichier contenant +// les parametres sous la forme [param=valeur] sur chaque ligne, puis +// appeler get<type>("param") pour acceder a un parametre. + +#ifndef CONFIGFILE_H +#define CONFIGFILE_H + +#include <map> +#include <iostream> +#include <fstream> +#include <stdio.h> +#include <stdlib.h> +#include <sstream> + +class ConfigFile{ + + public: + ConfigFile(const std::string& filename); + ~ConfigFile(); + + template<typename T> T get(const std::string& key) const; + + void process(const std::string& lineread); + + std::string toString() const; + + void printOut(const std::string& path) const; + + private: + std::string trim(const std::string& str); + + std::map<std::string, std::string> configMap; +}; + +#endif diff --git a/physnum/rap6/include/ConfigFile.tpp b/physnum/rap6/include/ConfigFile.tpp new file mode 100644 index 0000000..f92b761 --- /dev/null +++ b/physnum/rap6/include/ConfigFile.tpp @@ -0,0 +1,98 @@ +#include <iostream> +#include <fstream> +#include <stdio.h> +#include <stdlib.h> +#include <sstream> +#include <iomanip> +#include "ConfigFile.h" + +ConfigFile::ConfigFile(const std::string& fileName){ + std::ifstream file; + file.open(fileName.c_str()); + if (!file){ + std::cerr << "[ConfigFile] Impossible d'ouvrir le fichier " << fileName << std::endl; + } + else{ + std::string lineread; + while(getline(file, lineread)){ + process(lineread); + } + file.close(); + } +} + +ConfigFile::~ConfigFile(){} + +void ConfigFile::printOut(const std::string& path) const { + std::ofstream outputFile(path.c_str()); + if (outputFile.is_open()) + { + outputFile << toString() << std::endl; + } + outputFile.close(); +} + +std::string ConfigFile::toString() const { + std::string strToReturn; + + for (std::map<std::string,std::string>::const_iterator iter = configMap.begin(); iter != configMap.end(); ++iter) { + strToReturn.append(iter->first); + strToReturn.append("="); + strToReturn.append(iter->second); + strToReturn.append("\n"); + } + return strToReturn; +} + +void ConfigFile::process(const std::string& lineread) { + size_t commentPosition=trim(lineread).find('%',0); + if(commentPosition!=0 && trim(lineread).length()>0){ // End of line is counted as a character on some architectures + size_t equalPosition = lineread.find('=',1); + if(equalPosition==std::string::npos){ + std::cerr << "[ConfigFile] Ligne sans '=' : \"" << trim(lineread) << "\"" << std::endl; + }else{ + std::string key = trim(lineread.substr(0,equalPosition)); + std::string value = trim(lineread.substr(equalPosition+1,lineread.length())); + std::map<std::string, std::string>::const_iterator val = configMap.find(key); + if (val != configMap.end()) { + configMap.erase(key); + } + configMap.insert( std::pair<std::string, std::string>(key,value) ); + } + } +} + +template<typename T> T ConfigFile::get(const std::string& key) const{ + std::map<std::string, std::string>::const_iterator val = configMap.find(key); + T out; + if ( val != configMap.end() ) { + std::istringstream iss(val->second); + iss >> std::setprecision(12); + iss >> out; + std::cout << std::setprecision(12) << "\t" << key << "=" << out << std::endl; + }else{ + std::cerr << "[ConfigFile] Le parametre suivant est manquant : " << key << std::endl; + } + return out; +} + +template<> bool ConfigFile::get<bool>(const std::string& key) const{ + std::istringstream iss(configMap.find(key)->second); + bool result(false); + iss >> result; + if (iss.fail()){ + iss.clear(); + iss >> std::boolalpha >> result; + } + std::cout << "\t" << key << "=" << result << std::endl; + return result; +} + +std::string ConfigFile::trim(const std::string& str) +{ // Remove tabs and spaces at the beginning and end of a string + size_t first = str.find_first_not_of(" \t"); + if(first==std::string::npos) + return ""; + size_t last = str.find_last_not_of(" \t\r"); + return str.substr(first, (last-first+1)); +} diff --git a/physnum/rap6/include/debug.hpp b/physnum/rap6/include/debug.hpp new file mode 100644 index 0000000..5ca4baa --- /dev/null +++ b/physnum/rap6/include/debug.hpp @@ -0,0 +1,53 @@ +#pragma once + +#ifndef NPDEBUG +#define NPDEBUG + +#ifdef DEBUG +#include <iostream> +#include <sstream> + +#define __FILENAME__ (\ + __builtin_strrchr(__FILE__, '/') ? \ + __builtin_strrchr(__FILE__, '/') + 1 : __FILE__) + +#define npdebug(...); { \ + std::cerr << "[" << __FILENAME__ \ + << ":" << __LINE__ \ + << ", " << __func__ \ + << "] " ; \ + _np::debug(__VA_ARGS__); \ +} + +#define npinspect(msg, expr, ...); _np::inspect(msg, expr, __VA_ARGS__); + +namespace _np { + template<typename... Args> + inline void debug(const Args&... args) { + (std::cerr << ... << args) << std::endl; + } + + template<typename Msg, typename Expr, typename... Args> + inline Expr& inspect(const Msg& msg, Expr& expr, const Args&... args) { + npdebug(msg, expr, args...); + return expr; + } +} + +#else + +#define npdebug(...); {} +#define npinspect(...); {} + +namespace _np { + template<typename... Args> + void debug(Args&... args) {} + + template<typename Msg, typename Expr, typename... Args> + inline Expr& inspect(const Msg& msg, Expr& expr, const Args&... args) { + return expr; + } +} + +#endif // defined DEBUG +#endif // not defined NPDEBUG diff --git a/physnum/rap6/include/discret.h b/physnum/rap6/include/discret.h new file mode 100644 index 0000000..b68c9cf --- /dev/null +++ b/physnum/rap6/include/discret.h @@ -0,0 +1,17 @@ +#ifndef __DISCRET_H__ +#define __DISCRET_H__ + +struct discret +{ + std::size_t N; + double L; +}; + +<template class ...Points> +class interval +{ + std::tuplePoints + std::vector< +}; + +#endif diff --git a/physnum/rap6/include/field.h b/physnum/rap6/include/field.h new file mode 100644 index 0000000..53a719a --- /dev/null +++ b/physnum/rap6/include/field.h @@ -0,0 +1,55 @@ +#pragma once + +#ifndef __FIELD_H__ +#define __FIELD_H__ + +#include <ostream> +#include <vector> + +/* + * Wrapper containing information for a discretized point + * value = Temperature + * flag = 0 := can be modified + * + * Simpler usage: for an instance f, + * + * - *f <-> f.value + * - if (f) <-> if (f.flag) + */ +struct field +{ + field(double _value = 0, bool _flag = false) : value(_value), flag(_flag) {} + + double value; + bool flag; + + inline operator bool() const + { + return flag; + } + + inline double& operator*() + { + return value; + } + + inline const double& operator*() const + { + return value; + } +}; + +// grid mapping definition +typedef std::vector<std::vector<field>> grid_t; + +/* + * allows to print scalar field discretized grid + * + * Warning: it applies transposition, such that + * x direction = j + * y direction = i + * So that in matlab the matrix can be directly applied as it is + */ +std::ostream& operator<<(std::ostream& os, const grid_t& grid); + +#endif diff --git a/physnum/rap6/include/findiff.h b/physnum/rap6/include/findiff.h new file mode 100644 index 0000000..d0552f0 --- /dev/null +++ b/physnum/rap6/include/findiff.h @@ -0,0 +1,22 @@ +#ifndef __FINDIFF_H__ +#define __FINDIFF_H__ + +#include "symbasic/dvector.h" + +/*template<int i> +inline int Lambda(double k, double h) +{ + return (k + i +}*/ + +typedef DVector<double> rmap; +typedef Matrix; // TODO + +template<class H> +inline Matrix& compose(Matrix& matrix, int k, const H& h) +{ + matrix[k][k] += 1.0 / h[k]; + matrix[k+1][k] = matrix[k][k+1] = - 1.0 / h[k]; +} + +#endif diff --git a/physnum/rap6/include/geom.h b/physnum/rap6/include/geom.h new file mode 100644 index 0000000..9bcd600 --- /dev/null +++ b/physnum/rap6/include/geom.h @@ -0,0 +1,86 @@ +#ifndef __GEOM_H__ +#define __GEOM_H__ + +#include "vec2ext.h" +#include <istream> +#include <ostream> + +struct geom +{ + // char identifier for configuration files + virtual char chr() const = 0; + + virtual bool inside(double x, double y) const = 0; + + virtual ~geom() {} +}; + +/* + * rectangle struct + * + * Simply a wrapper which contains + * the top left vertex and the right bottom vertex + * in order to define a rectangular domain + */ +struct rect : public geom +{ + // tl = top left + // br = bottom right + Vec2D tl, br; + + rect(const Vec2D& _tl = {0, 0}, const Vec2D& _br = {0, 0}) + : tl(_tl), br(_br) {} + + rect(const rect& r) = default; + + virtual char chr() const override + { + return 'R'; + } + + virtual bool inside(double x, double y) const override; +}; + +struct circle : public geom +{ + Vec2D center; + double radius; + + circle(const Vec2D& _center = {0, 0}, const double& _radius = 1) + : center(_center), radius(_radius) {} + + circle(const circle&) = default; + + virtual char chr() const override + { + return 'C'; + } + + virtual bool inside(double x, double y) const override; +}; + +/* + * Allows to load a rectangles using the form {(x, y), (a, b)}, + * where + * + * (x, y) = top-left + * (a, b) = bottom-right + */ +std::istream& operator>>(std::istream& is, rect &r); + +// allows to print a rectangle in the format {(tl[0], tl[1]), (br[0], br[1])} +std::ostream& operator<<(std::ostream& os, const rect &r); + + +/* + * Allows to load a circle using the form {(x, y), r}, + * where + * + * (x, y) = center + * r = radius + */ +std::istream& operator>>(std::istream& is, circle &r); + +std::ostream& operator<<(std::ostream& os, const circle &r); + +#endif diff --git a/physnum/rap6/include/scheme.h b/physnum/rap6/include/scheme.h new file mode 100644 index 0000000..89fa283 --- /dev/null +++ b/physnum/rap6/include/scheme.h @@ -0,0 +1,25 @@ +#ifndef __SCHEME_H__ +#define __SCHEME_H__ + + +// runge kutta order 4 +// x = position +// v = speed +// t = time +// a = acceleration function, owned by a Base class +// returns the delta of position dx +// +template<class T, class Base> +T runge_kutta4(const T &args, double dt, double t, T (Base::*f)(const T&, double) const, Base *base) +{ +#define ACC(x, t) (base->*f)(x, t) + + T k1 = dt * ACC(args, t); + T k2 = dt * ACC(args + k1 / 2, t + dt / 2); + T k3 = dt * ACC(args + k2 / 2, t + dt / 2); + T k4 = dt * ACC(args + k3, t + dt); + + return args + (k1 + 2 * k2 + 2 * k3 + k4) / 6; +} + +#endif diff --git a/physnum/rap6/include/source.h b/physnum/rap6/include/source.h new file mode 100644 index 0000000..5d19e3e --- /dev/null +++ b/physnum/rap6/include/source.h @@ -0,0 +1,42 @@ +#ifndef __SOURCE_H__ +#define __SOURCE_H__ + +#include "geom.h" + +struct source +{ + double temperature; + geom * bounds; + bool fixed; + + // warning, passes bounds pointer ownership + source(double _T = 0, geom * _bounds = NULL, bool _fixed = false) + : temperature(_T), bounds(_bounds), fixed(_fixed) + { + } + + source(const source&) = delete; + source& operator=(const source& s) = delete; + + source(source&& s) + : bounds(s.bounds), temperature(s.temperature), fixed(s.fixed) + { + s.bounds = NULL; + } + + ~source() + { + if (bounds != NULL) + delete bounds; + } + + bool inside(double x, double y) const + { + return bounds->inside(x, y); + } +}; + +std::istream& operator>>(std::istream& is, source& s); +std::ostream& operator<<(std::ostream& is, const source& s); + +#endif diff --git a/physnum/rap6/include/symbasic/dvector.h b/physnum/rap6/include/symbasic/dvector.h new file mode 100644 index 0000000..99dadee --- /dev/null +++ b/physnum/rap6/include/symbasic/dvector.h @@ -0,0 +1,80 @@ +#ifndef __DVECTOR_H__ +#define __DVECTOR_H__ + +#include <vector> +#include <initializer_list> + +#include "vector.h" + +template<class T> +class DVector : public Vector<T> +{ + +public : + + DVector(size_t dim = 0); + DVector(const std::vector<T>&); + DVector(const std::initializer_list<T> &init); + + virtual const size_t size() const override; + + virtual T& operator[](size_t i) override; + virtual const T& operator[](size_t i) const override; + + DVector& operator<<(const T&); + + /* Plus / minus */ + DVector<T>& operator+=(const Vector<T>&); + DVector<T>& operator-=(const Vector<T>&); + + /* Scalar multiplication / division */ + DVector<T>& operator*=(double k); + DVector<T>& operator/=(double k); + + /* Reverse operator */ + DVector<T>& operator~(); + + /* Unit vector */ + DVector<T>& unit(); + +private: + + std::vector<T> components; +}; + +template<class T> +const DVector<T> operator+(DVector<T> , const Vector<T>& ); + +template<class T> +const DVector<T> operator-(DVector<T> , const Vector<T>& ); + +template<class T> +const DVector<T> operator*(DVector<T> , double); + +template<class T> +const DVector<T> operator/(DVector<T> , double); + +/*template<class T> +const DVector<T> operator+(const Vector<T>&, DVector<T>); + +template<class T> +const DVector<T> operator-(const Vector<T>&, DVector<T>);*/ + +template<class T> +const DVector<T> operator*(double, DVector<T>); + +template<class T> +const DVector<T> operator/(double, DVector<T>); + + +template<class T> +const DVector<T> operator-(DVector<T>); + +inline const DVector<double> operator^(const DVector<double>& v, const Vector<double>&); + +/* + * Include definitions file + */ +#include "dvector.tpp" + +#endif diff --git a/physnum/rap6/include/symbasic/dvector.tpp b/physnum/rap6/include/symbasic/dvector.tpp new file mode 100644 index 0000000..3aa7923 --- /dev/null +++ b/physnum/rap6/include/symbasic/dvector.tpp @@ -0,0 +1,147 @@ +#ifndef __DVECTOR_TPP__ +#define __DVECTOR_TPP__ + +template<class T> +DVector<T>::DVector(size_t dim) + : components(dim) +{ + +} + +template<class T> +DVector<T>::DVector(const std::vector<T>& init) + : components(init) +{ + +} + +template<class T> +DVector<T>::DVector(const std::initializer_list<T> &init) + : components(init.begin(), init.end()) +{ + +} + +template<class T> +const size_t DVector<T>::size() const +{ + return components.size(); +} + +template<class T> +T& DVector<T>::operator[](size_t i) +{ + return components[i]; +} + +template<class T> +const T& DVector<T>::operator[](size_t i) const +{ + return components[i]; +} + +template<class T> +DVector<T>& DVector<T>::operator<<(const T &value) +{ + components.push_back(value); +} + +template<class T> +DVector<T>& DVector<T>::operator+=(const Vector<T> &w) +{ + return static_cast<DVector<T>&>(this->add(w)); +} + +template<class T> +DVector<T>& DVector<T>::operator-=(const Vector<T> &w) +{ + return static_cast<DVector<T>&>(this->sub(w)); +} + +template<class T> +DVector<T>& DVector<T>::operator*=(double k) +{ + return static_cast<DVector<T>&>(this->mult(k)); +} + +template<class T> +DVector<T>& DVector<T>::operator/=(double k) +{ + return static_cast<DVector<T>&>(this->divide(k)); +} + +template<class T> +DVector<T>& DVector<T>::operator~() +{ + return static_cast<DVector<T>&>(this->reverse()); +} + +template<class T> +DVector<T>& DVector<T>::unit() +{ + return static_cast<DVector<T>&>(this->_unit()); +} + +template<class T> +const DVector<T> operator-(DVector<T> v) +{ + return ~v; +} + +template<class T> +const DVector<T> operator+(DVector<T> v, const Vector<T>& w) +{ + return v += w; +} + +template<class T> +const DVector<T> operator-(DVector<T> v, const Vector<T>& w) +{ + return v -= w; +} + +template<class T> +const DVector<T> operator*(DVector<T> v, double k) +{ + return v *= k; +} + +template<class T> +const DVector<T> operator/(DVector<T> v, double k) +{ + return v /= k; +} + +template<class T> +const DVector<T> operator*(double k, DVector<T> v) +{ + return v *= k; +} + +template<class T> +const DVector<T> operator/(double k, DVector<T> v) +{ + return v /= k; +} + +inline const DVector<double> operator^(const DVector<double>& v, const Vector<double>& w) +{ + DVector<double> u; + + size_t n = (v.size() <= w.size()) ? v.size() : w.size(); + + if (n > 2) + { + u << v[1] * w[2] - v[2] * w[1]; + u << v[2] * w[0] - v[0] * w[2]; + u << v[0] * w[1] - v[1] * w[0]; + + } else if (n == 2) + u << v[0] * w[1] - v[1] * w[0]; + else if (n == 1) + u << v[0] * w[0]; + + return u; +} + +#endif diff --git a/physnum/rap6/include/symbasic/svector.h b/physnum/rap6/include/symbasic/svector.h new file mode 100644 index 0000000..c53ae0f --- /dev/null +++ b/physnum/rap6/include/symbasic/svector.h @@ -0,0 +1,78 @@ +#ifndef __STATIC_VECTOR_H__ +#define __STATIC_VECTOR_H__ + +#include <array> +#include "vector.h" + +template<class T, std::size_t N> +class SVector : public Vector<T> +{ + +public: + + /* Constructors */ + + SVector(const T &value = T()); + SVector(const std::initializer_list<T> &init); + SVector(const SVector<T, N>&); + + SVector& operator=(const SVector<T, N>&); + + virtual const size_t size() const override; + + virtual T& operator[](size_t i) override; + virtual const T& operator[](size_t i) const override; + + /* Plus / minus */ + SVector<T, N>& operator+=(const Vector<T>&); + SVector<T, N>& operator-=(const Vector<T>&); + + /* Scalar multiplication / division */ + SVector<T, N>& operator*=(double k); + SVector<T, N>& operator/=(double k); + + /* Reverse operator */ + SVector<T, N>& operator~(); + + /* Unit vector */ + SVector<T, N>& unit(); + +private: + + std::array<T, N> components; +}; + +template<class T, std::size_t N> +const SVector<T, N> operator+(SVector<T, N> , const Vector<T>& ); + +template<class T, std::size_t N> +const SVector<T, N> operator-(SVector<T, N> , const Vector<T>& ); + +template<class T, std::size_t N> +const SVector<T, N> operator*(SVector<T, N> , double); + +template<class T, std::size_t N> +const SVector<T, N> operator/(SVector<T, N> , double); + +template<class T, std::size_t N> +const SVector<T, N> operator*(double, SVector<T, N>); + +template<class T, std::size_t N> +const SVector<T, N> operator/(double, SVector<T, N>); + +template<class T, std::size_t N> +const SVector<T, N> operator-(SVector<T, N>); + +/* + * Prototipe for cross product + */ +inline const SVector<double, 3> operator^(const SVector<double, 3>& v, const SVector<double, 3>& w); +inline double operator^(const SVector<double, 2>& v, const SVector<double, 2>& w); +inline double operator^(const SVector<double, 1>& v, const SVector<double, 1>& w); + +/* + * Include definitions file + */ +#include "svector.tpp" + +#endif // __SVECTOR_H__ diff --git a/physnum/rap6/include/symbasic/svector.tpp b/physnum/rap6/include/symbasic/svector.tpp new file mode 100644 index 0000000..e054049 --- /dev/null +++ b/physnum/rap6/include/symbasic/svector.tpp @@ -0,0 +1,162 @@ +#ifndef __SVECTOR_TPP__ +#define __SVECTOR_TPP__ + +template<class T, std::size_t N> +SVector<T, N>::SVector(const T &value) +{ + components.fill(value); +} + +template<class T, std::size_t N> +SVector<T, N>::SVector(const std::initializer_list<T> &init) +{ + size_t n = (init.size() < N) ? init.size() : N; + + /* copy init list without include <algorithm> */ + size_t index = 0; + + /* copy the initializer list until its size */ + for (auto c : init) + { + if (index >= n) + break; + + components[index++] = c; + } + + /* fill the remaining components with zero */ + while(index < N) + components[index++] = 0; +} + +template<class T, std::size_t N> +SVector<T, N>::SVector(const SVector<T, N> &cpy) +{ + this->components = cpy.components; +} + +template<class T, std::size_t N> +SVector<T, N>& SVector<T, N>::operator=(const SVector<T, N> &v) +{ + this->components = v.components; +} + +template<class T, std::size_t N> +const size_t SVector<T, N>::size() const +{ + return N; +} + +/* Access operators */ +template<class T, std::size_t N> +T& SVector<T, N>::operator[](size_t i) +{ + return components[i]; +} + +template<class T, std::size_t N> +const T& SVector<T, N>::operator[](size_t i) const +{ + return components[i]; +} + +template<class T, std::size_t N> +SVector<T, N>& SVector<T, N>::operator+=(const Vector<T> &w) +{ + return static_cast<SVector<T, N>&>(this->add(w)); +} + +template<class T, std::size_t N> +SVector<T, N>& SVector<T, N>::operator-=(const Vector<T> &w) +{ + return static_cast<SVector<T, N>&>(this->sub(w)); +} + +template<class T, std::size_t N> +SVector<T, N>& SVector<T, N>::operator*=(double k) +{ + return static_cast<SVector<T, N>&>(this->mult(k)); +} + +template<class T, std::size_t N> +SVector<T, N>& SVector<T, N>::operator/=(double k) +{ + return static_cast<SVector<T, N>&>(this->divide(k)); +} + +template<class T, std::size_t N> +SVector<T, N>& SVector<T, N>::operator~() +{ + return static_cast<SVector<T, N>&>(this->reverse()); +} + +template<class T, std::size_t N> +SVector<T, N>& SVector<T, N>::unit() +{ + return static_cast<SVector<T, N>&>(this->_unit()); +} + +template<class T, std::size_t N> +const SVector<T, N> operator+(SVector<T, N> v, const Vector<T>& w) +{ + return v += w; +} + +template<class T, std::size_t N> +const SVector<T, N> operator-(SVector<T, N> v, const Vector<T>& w) +{ + return v -= w; +} + +template<class T, std::size_t N> +const SVector<T, N> operator*(SVector<T, N> v, double k) +{ + return v *= k; +} + +template<class T, std::size_t N> +const SVector<T, N> operator/(SVector<T, N> v, double k) +{ + return v /= k; +} + +template<class T, std::size_t N> +const SVector<T, N> operator*(double k, SVector<T, N> v) +{ + return v *= k; +} + +template<class T, std::size_t N> +const SVector<T, N> operator/(double k, SVector<T, N> v) +{ + return v /= k; +} + +template<class T, std::size_t N> +const SVector<T, N> operator-(SVector<T, N> v) +{ + return ~v; +} + +inline const SVector<double, 3> operator^(const SVector<double, 3>& v, const SVector<double, 3>& w) +{ + SVector<double, 3> u = {0, 0, 0}; + + u[0] = v[1] * w[2] - v[2] * w[1]; + u[1] = v[2] * w[0] - v[0] * w[2]; + u[2] = v[0] * w[1] - v[1] * w[0]; + + return u; +} + +inline double operator^(const SVector<double, 2>& v, const SVector<double, 2>& w) +{ + return v[0] * w[1] - v[1] * w[0]; +} + +inline double operator^(const SVector<double, 1>& v, const SVector<double, 1>& w) +{ + return v[0] * w[0]; +} + +#endif diff --git a/physnum/rap6/include/symbasic/vector.h b/physnum/rap6/include/symbasic/vector.h new file mode 100644 index 0000000..f8574dc --- /dev/null +++ b/physnum/rap6/include/symbasic/vector.h @@ -0,0 +1,133 @@ +#ifndef __VECTOR_H__ +#define __VECTOR_H__ + +#ifdef _USE_VECTOR_OSTREAM +#include <iosfwd> +#endif + +#include <cstddef> + +template <class T> +class Vector { + +public: + + typedef std::size_t size_t; + + class iterator + { + size_t i; + Vector<T> *v; + + public: + + iterator(Vector<T> *v, size_t i = 0) + : v(v), i(i) + { + } + + iterator& operator++() + { + ++i; + } + + iterator& operator--() + { + --i; + } + + bool operator!=(const iterator& it) + { + return this->i != it.i; + } + + bool operator==(const iterator& it) + { + return !(*this != it); + } + + T& operator*() + { + return (*v)[i]; + } + }; + + typedef const iterator const_iterator; + + /* For auto implementation */ + /* it allows to loop through components + * + * for (auto comp : v) + * { + * // loop content + * } + */ + + iterator begin(); + iterator end(); + + const_iterator begin() const; + const_iterator end() const; + + //Virtual accesseurs + + virtual const size_t size() const = 0; + + /* Accède à la i-ème valeur */ + virtual T& operator[](size_t i) = 0; + virtual const T& operator[](size_t i) const = 0; + + /* Return the module */ + double module() const; + double sq_module() const; // squared module + + //Operateurs + + /* Operateurs de comparaison */ + bool operator==(const Vector<T>&) const; + bool operator!=(const Vector<T>&) const; + + bool operator>=(const Vector<T>&) const; + bool operator<=(const Vector<T>&) const; + + bool operator>=(double) const; + bool operator<=(double) const; + + bool operator>(const Vector<T>&) const; + bool operator<(const Vector<T>&) const; + + bool operator>(double) const; + bool operator<(double) const; + + /* Dot product */ + double operator*(const Vector<T>&) const; + + /* For printing a Vector on a ostream */ +#ifdef _USE_VECTOR_OSTREAM + void affiche(std::ostream&) const; +#endif + +protected: // not stable operations + + /* Plus / minus */ + Vector<T>& add(const Vector<T>&); + Vector<T>& sub(const Vector<T>&); + + /* Scalar multiplication / division */ + Vector<T>& mult(double k); + Vector<T>& divide(double k); + + Vector<T>& reverse(); + Vector<T>& _unit() const; +}; + +#ifdef _USE_VECTOR_OSTREAM +// forward declaration of ostream + +template<class T> +std::ostream& operator<<(std::ostream&, const Vector<T>&); +#endif + +#include "vector.tpp" + +#endif diff --git a/physnum/rap6/include/symbasic/vector.tpp b/physnum/rap6/include/symbasic/vector.tpp new file mode 100644 index 0000000..d622940 --- /dev/null +++ b/physnum/rap6/include/symbasic/vector.tpp @@ -0,0 +1,213 @@ +#ifndef __VECTOR_TPP__ +#define __VECTOR_TPP__ + +/* for auto loop */ + +template<class T> +typename Vector<T>::iterator Vector<T>::begin() +{ + return iterator(this); +} + +template<class T> +typename Vector<T>::iterator Vector<T>::end() +{ + return iterator(this, this->size() - 1); +} + +template<class T> +typename Vector<T>::const_iterator Vector<T>::begin() const +{ + return iterator(this); +} + +template<class T> +typename Vector<T>::const_iterator Vector<T>::end() const +{ + return iterator(this, this->size() - 1); +} + +template<class T> +bool Vector<T>::operator==(const Vector<T>& v) const { + + if(this->size() != v.size()) + return false; + + //Retourne déjà faux si les deux vecteurs n'ont pas la même dimension + for(size_t i(0); i < v.size(); ++i){ + + if ((*this)[i] != v[i]) + return false; + //Lors que le programme rencontre deux coordonnées différenentes respectivement pour la même dimension des deux vecteurs + //la fonction retourne false + } + return true; +} + +template<class T> +bool Vector<T>::operator>=(const Vector<T>& v) const +{ + return sq_module() >= v.sq_module(); +} + +template<class T> +bool Vector<T>::operator<=(const Vector<T>& v) const +{ + return sq_module() <= v.sq_module(); +} + +template<class T> +bool Vector<T>::operator>=(double m) const +{ + return sq_module() >= m * m; +} + +template<class T> +bool Vector<T>::operator<=(double m) const +{ + return sq_module() <= m * m; +} + +template<class T> +bool Vector<T>::operator>(const Vector<T>& v) const +{ + return sq_module() > v.sq_module(); +} + +template<class T> +bool Vector<T>::operator<(const Vector<T>& v) const +{ + return sq_module() < v.sq_module(); +} + +template<class T> +bool Vector<T>::operator>(double m) const +{ + return sq_module() > m * m; +} + +template<class T> +bool Vector<T>::operator<(double m) const +{ + return sq_module() < m * m; +} + +template<class T> +bool Vector<T>::operator!=(const Vector<T>& w) const +{ + return !(*this == w); +} + +//Operateurs somme et produit + +template<class T> +Vector<T>& Vector<T>::add(const Vector<T>& w) +{ + size_t n = (this->size() <= w.size()) ? this->size() : w.size(); + + for (size_t i= 0; i < n; ++i) + (*this)[i] += w[i]; + + return *this; +} + + +template<class T> +Vector<T>& Vector<T>::sub(const Vector<T>& w) +{ + size_t n = (this->size() <= w.size()) ? this->size() : w.size(); + + for (size_t i = 0; i < n; ++i) + (*this)[i] -= w[i]; + + return *this; +} + +template<class T> +Vector<T>& Vector<T>::mult(double k) +{ + for (size_t i = 0; i < this->size(); ++i) + (*this)[i] *= k; + + return *this; +} + +template<class T> +Vector<T>& Vector<T>::divide(double k) +{ + for (size_t i = 0; i < this->size(); ++i) + (*this)[i] /= k; + + return *this; +} + +template<class T> +Vector<T>& Vector<T>::_unit() const +{ + return (*this) /= this->module(); +} + +/* Dot product */ +template<class T> +double Vector<T>::operator*(const Vector<T>& w) const +{ + double x = 0; + size_t n = (this->size() <= w.size()) ? this->size() : w.size(); + + for (size_t i = 0; i < n; ++i) + x += (*this)[i] * w[i]; + + return x; +} + +template<class T> +Vector<T>& Vector<T>::reverse() +{ + return (*this)*=-1.0; +} + +template<class T> +double Vector<T>::sq_module() const +{ + return (*this) * (*this); +} + +template<class T> +double Vector<T>::module() const +{ + return sqrt(this->sq_module()); +} + +/* Implementation of ostream operator overloading */ + +/* Le prototype je le mets dans un header séparé parce que + * il ne fait pas partie de l'implementation d'un vecteur + * Une classe l'utilisant n'as pas besoin de l'imprimer + */ + +#ifdef _USE_VECTOR_OSTREAM +#include <ostream> + +template<class T> +void Vector<T>::affiche(std::ostream& os) const +{ + os << "("; + + for (size_t i(0); i < (size() - 1); ++i) + { + os << (*this)[i] << ", "; + } + + os << (*this)[size()-1] << ")"; +} + +// Flux de sortie permettant d'afficher les coordonnées du vecteur +template<class T> +std::ostream& operator<<(std::ostream& os, const Vector<T>& v) +{ + v.affiche(os); + return os; +} +#endif + +#endif //__VECTOR_TPP__ diff --git a/physnum/rap6/include/vec2ext.h b/physnum/rap6/include/vec2ext.h new file mode 100644 index 0000000..21149c3 --- /dev/null +++ b/physnum/rap6/include/vec2ext.h @@ -0,0 +1,49 @@ +#pragma once + +#ifndef __VEC2EXT_H__ +#define __VEC2EXT_H__ + +#define _USE_VECTOR_OSTREAM +#include "symbasic/svector.h" + +#include <istream> + +// two dimensional vector definition +typedef SVector<double, 2> Vec2D; + +/* + * Allows to load any kind of bi-dimensional vectors + * using the form (x, y) + * + * Example: (2, 3) + * + * for SVector declaration, take a look to "symbasic/svector.h" + */ +template<class T> +std::istream& operator>>(std::istream& is, SVector<T, 2>& v) +{ + char delim; + is >> std::skipws; + is >> delim; + + if (delim != '(') + throw '('; + + is >> v[0]; + + is >> delim; + + if (delim != ',') + throw ','; + + is >> v[1]; + + is >> delim; + + if (delim != ')') + throw ')'; + + return is; +} + +#endif diff --git a/physnum/rap6/src/exercice6.cpp b/physnum/rap6/src/exercice6.cpp new file mode 100644 index 0000000..a1a2713 --- /dev/null +++ b/physnum/rap6/src/exercice6.cpp @@ -0,0 +1,273 @@ +#include <iostream> +#include <string> +#include <cmath> +#include <vector> +#include <fstream> +#include "ConfigFile.tpp" +#include "debug.hpp" + +#include "geom.h" +#include "field.h" +#include "source.h" + +// output mode definition +#define POINT 0 +#define GRID 1 +#define ALL_GRID 2 + +using namespace std; + +// constants, defined in main + +// output mode +int mode; + +// in case of mode=0, specify what is point to observe +SVector<int, 2> outpoint; + +// main output name +string output; + +// geometrie +double L; + +// Border temperature +//double Tc; +//double Tf; +double Tb; + +// conduction constant +double kappa; + +// final time and iterations +double tfin; +int max_iter; + +// falloff check error (eps) +double eps; + +// Discretization +int N; +double dt; +double h; +double alpha; + +// output file streams +ofstream output_T; + +double error(const double&, const double&); + +std::pair<grid_t, double> diffusion(const grid_t& grid) +{ + grid_t out(grid); + + double r = 0; + + for (int i = 1; i < N; ++i) + { + for (int j = 1; j < N; ++j) + { + // evolution step + /* + * Deduce if it can be modified, else + * simply copy the previous value + */ + *out[i][j] = (grid[i][j]) ? + *grid[i][j] : + *grid[i][j] * (1 - 4 * alpha) + alpha * (*grid[i+1][j] + *grid[i-1][j] + *grid[i][j+1] + *grid[i][j-1]); + + // error computation + double err = error(*grid[i][j], *out[i][j]); + + // take the maximum + if (err > r) + r = err; + } + } + + return std::pair<grid_t, double>(out, r); +} + +// computes the error +double error(const double& old, const double& nu) +{ + return abs(nu - old) / dt; +} + +/* + * Main loop function, + * finally the core of the symulation + */ +grid_t loop(grid_t&& init) +{ + std::pair<grid_t, double> p(init, 0); + int iter = 0; + + /* + * In case of hard debugging issues, + * print all the flag grid + * + * for(int i(0);i<N+1;++i) + { + for(int j(0);j<N+1;++j) + { + cout << bool(p.first[i][j]) << " "; + } + cout << endl; + }*/ + + if (mode == POINT || mode == GRID) + { + output_T.open((output + "_T.out").c_str()); // Temperature au temps final + output_T.precision(15); + } + + do { + p = diffusion(p.first); + + // print a single point, mode = 0 + if (mode == POINT) + output_T << iter*dt << " " << *p.first[outpoint[0]][outpoint[1]] << endl; + else if (mode == ALL_GRID) + { + std::stringstream ss; + ss << output << ".iter=" << (iter+1) << "_T.out"; + + output_T.open(ss.str().c_str()); // Temperature au temps final + output_T.precision(15); + + output_T << p.first; + output_T.close(); + } + + } while (p.second > eps && iter++ < max_iter); + + // print all the grid, mode = 1 + if (mode == GRID) + output_T << p.first; + + // close output file + if (mode == POINT || mode == GRID) + output_T.close(); + + return p.first; +} + +template<class T> +bool source_check(grid_t& grid, int i, int j, T&& arg) +{ + bool out; + + // if it finds the domain, it sets the temperature + // to the corresponding one + if (out = arg.inside(i * h, j * h)) + grid[i][j] = field(arg.temperature, arg.fixed); + + return out; +} + +template<class ...Args> +grid_t initial(Args&& ...args) +{ + grid_t grid(N+1, std::vector<field>(N+1)); + + // Border initialization + for (std::size_t i = 0; i <= N; ++i) + grid[0][i] = + grid[N][i] = + grid[i][0] = + grid[i][N] = field(Tb, true); + + /* + * Internal initialization + * Check whether a particular domain should not be + * affected by the simulation + */ + for (std::size_t i = 1; i < N; ++i) + { + for (std::size_t j = 1; j < N; ++j) + { + /* + * Expanding parameters pack arguments: + * generating array of bool. + * + * If any value of that array is true, then + * the point belongs to one of the sources. + */ + bool flags[] = { source_check(grid, i, j, std::forward<Args>(args))... }; + + bool somewhere = false; + + // Verify booleans + for (bool flag : flags) + somewhere = somewhere || flag; + + /* + * A point belongs to the border if and only if + * none of the surfaces is container of that point + */ + if (!somewhere) + grid[i][j] = field(Tb, false); + } + } + + return grid; +} + +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]); + + // mode + mode = configFile.get<int>("mode"); + outpoint = configFile.get<SVector<int, 2>>("outpoint"); + + // Geometrie: + L = configFile.get<double>("L"); + + source chaud(configFile.get<source>("chaud")); + source froid(configFile.get<source>("froid")); + + // Temperatures: + Tb = configFile.get<double>("Tb"); + kappa = configFile.get<double>("kappa"); + + // Duree de la simulation: + tfin = configFile.get<double>("tfin"); + eps = configFile.get<double>("eps"); // Condition d'arret si etat stationnaire + + // Discretisation: + N = configFile.get<int>("N"); // Nombre d'intervalles dans chaque dimension + dt = configFile.get<double>("dt"); + h = L/N; + alpha = kappa * dt / h / h; + + // iteration + max_iter = tfin / dt; + + // Fichiers de sortie: + output = configFile.get<string>("output"); + /* + * Start the symulation loop + * Return the resulting grid + * + * initial() return type is taken as rvalue by the loop function, + * which means that no copy constructor of grid_t is called + * + * initial() function take an arbitrary number of sources + * as parameter pack (evaluated at compile time) + * + * Warning: compile with standard c++17 + */ + grid_t T = loop( initial(chaud, froid) ); + + return 0; +} + diff --git a/physnum/rap6/src/field.cpp b/physnum/rap6/src/field.cpp new file mode 100644 index 0000000..ae70594 --- /dev/null +++ b/physnum/rap6/src/field.cpp @@ -0,0 +1,30 @@ +#include "field.h" + +/* + * allows to print scalar field discretized grid + * + * Warning: it applies transposition, such that + * x direction = j + * y direction = i + * So that in mathlab the matrix can be directly applied as it is + */ +std::ostream& operator<<(std::ostream& os, const grid_t& grid) +{ + /* + * Notice: print mapping order is placed such that + * x = j * h + * y = i * h + */ + size_t N = grid.size(); + + for (int j = 0; j < N; ++j) + { + for (int i = 0; i < N; ++i) + os << *grid[i][j] << " "; + + os << std::endl; + } + + return os; +} + diff --git a/physnum/rap6/src/geom.cpp b/physnum/rap6/src/geom.cpp new file mode 100644 index 0000000..2744713 --- /dev/null +++ b/physnum/rap6/src/geom.cpp @@ -0,0 +1,108 @@ +#include "geom.h" + +#include "debug.hpp" + +bool rect::inside(double x, double y) const +{ + //npdebug("Checking: (", x, ", ", y, ")", " in ", tl, br) + return tl[0] < x && x < br[0] && + br[1] < y && y < tl[1]; +} + +bool circle::inside(double x, double y) const +{ + //npdebug("Checking: (", x, ", ", y, ")", " in ", tl, br) + double rx = (x - center[0]); + double ry = (y - center[1]); + + return (rx * rx + ry * ry) < radius * radius; +} + +/* + * Allows to load a rectangles using the form {(x, y), (a, b)}, + * where + * + * (x, y) = top-left + * (a, b) = bottom-right + */ +std::istream& operator>>(std::istream& is, rect &r) +{ + is >> std::skipws; + char delim; + + // begin + is >> delim; + + if (delim != '{') + throw '{'; + + is >> r.tl; + + npdebug("") + + is >> delim; + + if (delim != ',') + throw ','; + + is >> r.br; + + npdebug("") + + // end + is >> delim; + + if (delim != '}') + throw '}'; + + npdebug("") + + is >> std::noskipws; + + return is; +} + +// allows to print a rectangle in the format {(tl[0], tl[1]), (br[0], br[1])} +std::ostream& operator<<(std::ostream& os, const rect &r) +{ + os << "{" << r.tl << ", " << r.br << "}"; + return os; +} + +std::istream& operator>>(std::istream& is, circle &r) +{ + is >> std::skipws; + char delim; + + // begin + is >> delim; + + if (delim != '{') + throw '{'; + + is >> r.center; + + is >> delim; + + if (delim != ',') + throw ','; + + is >> r.radius; + + // end + is >> delim; + + if (delim != '}') + throw '}'; + + is >> std::noskipws; + + return is; +} + +// allows to print a rectangle in the format {(tl[0], tl[1]), (br[0], br[1])} +std::ostream& operator<<(std::ostream& os, const circle &r) +{ + os << "{" << r.center << ", " << r.radius << "}"; + return os; +} diff --git a/physnum/rap6/src/source.cpp b/physnum/rap6/src/source.cpp new file mode 100644 index 0000000..5d0f6ca --- /dev/null +++ b/physnum/rap6/src/source.cpp @@ -0,0 +1,95 @@ +#include "source.h" + +#include "debug.hpp" + +std::istream& operator>>(std::istream& is, source &s) +{ + is >> std::skipws; + char delim; + + // begin + is >> delim; + + char type; + + if (delim != '{') + throw '{'; + + is >> type; + + is >> delim; + + if (delim != ',') + throw ','; + + is >> s.temperature; + + is >> delim; + + if (delim != ',') + throw ','; + + switch(type) + { + case 'R': + is >> *dynamic_cast<rect*>(s.bounds = new rect()); + break; + case 'C': + is >> *dynamic_cast<circle*>(s.bounds = new circle()); + break; + default: + throw type; + } + + // be sure + is >> std::skipws; + + is >> delim; + + if (delim != ',') + throw ','; + + int fixed_buf; + + // 0 or 1 + is >> fixed_buf; + + s.fixed = static_cast<int>(fixed_buf); + + // end + is >> delim; + + if (delim != '}') + throw '}'; + + is >> std::noskipws; + + return is; +} + +// allows to print a rectangle in the format {(tl[0], tl[1]), (br[0], br[1])} +std::ostream& operator<<(std::ostream& os, const source &s) +{ + char type = s.bounds->chr(); + + os << "{" << type << ", " << s.temperature << ", "; + + switch (type) + { + case 'R': + os << *dynamic_cast<rect*>(s.bounds); + break; + + case 'C': + os << *dynamic_cast<circle*>(s.bounds); + break; + + default: + throw "Cannot cast bounds for source"; + } + + os << ", " << static_cast<int>(s.fixed) << "}"; + + return os; +} +