fork download
copy 
  1. #include <type_traits>
  2. #include <initializer_list>
  3.  
  4. #include <vector>
  5. #include <iostream>
  6.  
  7. // ------------------------------- MESHDIMENSION----------------------- //
  8. struct MeshDimension
  9. {
  10.     MeshDimension(const int numCells, const double min, const double max):
  11.         numCells_(numCells),
  12.         minVal_(min),
  13.         maxVal_(max)
  14.     {}
  15.     const int numCells_;
  16.     const double minVal_, maxVal_;
  17. };
  18. #include <cassert>
  19. namespace detail {
  20.     enum class enabler {};
  21. }
  22. template <bool Condition>
  23. using EnableIf =
  24.     typename std::enable_if<Condition, detail::enabler>::type;
  25. template<typename... Conds>
  26. struct and_ : std::true_type {};
  27. template<typename Cond, typename... Conds>
  28. struct and_<Cond, Conds...>
  29.         : std::conditional<Cond::value, and_<Conds...>,
  30.         std::false_type>::type {};
  31. template<typename Target, typename... Ts>
  32. using areT = and_<std::is_same<Ts, Target>...>;
  33.  
  34. // -------------------------------------- MESH.H -------------------------//
  35. template<size_t meshDim>
  36. class Mesh
  37. {
  38. public:
  39.     // Constructor
  40.     template<class... DimArgs>
  41.     Mesh(DimArgs... dimensions) {
  42.         // Check the dimensions are consistent
  43.         checkDims<DimArgs...>();
  44.         // Push values to the dimSize vector
  45.         pushDimensions<DimArgs...>(dimensions...);
  46.         // Store positions (templated to call the correct function)
  47.         prepareVectors<DimArgs...>();
  48.     }
  49.  
  50.     const size_t numCells() const;
  51.     // Always defined (can't have a 0D mesh)
  52.     const size_t xCells() const { return dimSize_[0]; }
  53.     // Get size if 2D or greater, return 0 otherwise
  54.     template<size_t meshD = meshDim, EnableIf<meshD>=2>...>
  55.     const size_t yCells() const { return dimSize_[1]; }
  56.     template<size_t meshD = meshDim, EnableIf<meshD<2>...>
  57.     const size_t yCells() const { return 0; }
  58.     template<size_t meshD = meshDim, EnableIf<meshD>=3>...>
  59.     const size_t zCells() const { return dimSize_[2]; }
  60.     template<size_t meshD = meshDim, EnableIf<meshD<3>...>
  61.     const size_t zCells() const { return 0; }
  62.  
  63.     const double x(const int &i) const { return getPosition(0, i); }
  64.     template<size_t meshD = meshDim, EnableIf<meshD>=2>...>
  65.     const double y(const int &i) const { return getPosition(1, i); }
  66.     template<size_t meshD = meshDim, EnableIf<meshD>=3>...>
  67.     const double z(const int &i) const { return getPosition(2, i); }
  68.  
  69.     template<size_t rhsMD>
  70.     bool operator==(const Mesh<rhsMD> &rhs) const;
  71.  
  72. private:
  73.     bool constSpacing_ = true;
  74.     std::vector<double> dimMin_;
  75.     std::vector<double> dimMax_;
  76.     std::vector<size_t> dimSize_;
  77.     std::vector<double> position_[meshDim];
  78.     void placeCentres(const int d, const bool constSpacing);
  79.     const double getPosition(const size_t d, const int &i) const;
  80.  
  81.     template<class... numArgs>
  82.     void checkDims() const;
  83.     void checkBounds(std::vector<int> idxs) const;
  84.  
  85.     template<typename... T, EnableIf<areT<MeshDimension, T...>::value>...>
  86.     void pushDimensions(T... dims);
  87.     template<typename... T, EnableIf<areT<int, T...>::value>...>
  88.     void pushDimensions(T... dims);
  89.  
  90.     template<typename... T, EnableIf<areT<MeshDimension, T...>::value>...>
  91.     void prepareVectors();
  92.     template<typename... T, EnableIf<areT<int, T...>::value>...>
  93.     void prepareVectors();
  94.  
  95.     enum class ScalingType {
  96.         Pivot,
  97.         Hyperbolic,
  98.         Exponential
  99.     };
  100. };
  101.  
  102. // ------------------------- MESH.HXX ---------------------------------//
  103. #include <math.h>
  104.  
  105. template<size_t meshDim>
  106. template<size_t rhsMD>
  107. bool Mesh<meshDim>::operator==(const Mesh<rhsMD> &rhs) const
  108. {
  109.     if (meshDim != rhsMD) {
  110.         return false;
  111.     }
  112.     if ((dimSize_ != rhs.dimSize_)
  113.             || (dimMin_ != rhs.dimMin_)
  114.             || (dimMax_ != rhs.dimMax_)) {
  115.         return false;
  116.     }
  117.     if (constSpacing_ != rhs.constSpacing_) {
  118.         return false;
  119.     }
  120.     return true;
  121. }
  122.  
  123.  
  124. // Assert that number of arguments is the same as number of dimensions.
  125. template<size_t meshDim>
  126. template<class... numArgs>
  127. void Mesh<meshDim>::checkDims() const {
  128.     static_assert(sizeof...(numArgs)==meshDim, "Wrong number of dimensions");
  129. }
  130.  
  131. // pushDimensions function, for int and MeshDimension type arguments
  132. template<size_t meshDim>
  133. template<typename... T, EnableIf<areT<MeshDimension, T...>::value>...>
  134. void Mesh<meshDim>::pushDimensions(T... dims) {
  135.     std::cout << "Constructed with MeshDimension structs" << std::endl;
  136.     std::vector<MeshDimension> list{dims...};
  137.     for (unsigned int i=0; i<meshDim; i++) {
  138.         dimSize_.push_back(list[i].numCells_);
  139.         dimMin_.push_back(list[i].minVal_);
  140.         dimMax_.push_back(list[i].maxVal_);
  141.     }
  142. }
  143. template<size_t meshDim>
  144. template<typename... T, EnableIf<areT<int, T...>::value>...>
  145. void Mesh<meshDim>::pushDimensions(T... dims) {
  146.     std::cout << "Constructed with ints. Must also add min/max values."
  147.               << std::endl;
  148.     std::vector<int> list{dims...};
  149.     for (unsigned int i=0; i<meshDim; i++) {
  150.         dimSize_.push_back(list[i]);
  151.     }
  152. }
  153.  
  154. // prepareVectors function, no-op if ints passed
  155. template<size_t meshDim>
  156. template<typename... T, EnableIf<areT<MeshDimension, T...>::value>...>
  157. void Mesh<meshDim>::prepareVectors() {
  158.     std::cout << "meshDim structs" << std::endl;
  159.     for (unsigned int d=0; d<meshDim; d++) {
  160.         placeCentres(d, constSpacing_);
  161.     }
  162. }
  163. template<size_t meshDim>
  164. template<typename... T, EnableIf<areT<int, T...>::value>...>
  165. void Mesh<meshDim>::prepareVectors() {
  166.     std::cout << "no-op: ints" << std::endl;
  167. }
  168.  
  169. template<size_t meshDim>
  170. void Mesh<meshDim>::placeCentres(const int d, const bool constSpacing) {
  171.     position_[d].resize(dimSize_[d]);
  172.     const double range = dimMax_[d] - dimMin_[d];
  173.     if (constSpacing) {
  174.         // Constant spacing
  175.         const double dx = range / dimSize_[d];
  176.         for (unsigned int i=0; i<dimSize_[d]; i++) {
  177.             position_[d][i] = (dimMin_[d] + (dx*(i+0.5)));
  178.         }
  179.     } else {
  180.         // Set the choice of method
  181.         constexpr ScalingType scalingType = ScalingType::Pivot;
  182.         switch (scalingType)
  183.         {
  184.             case ScalingType::Exponential:
  185.             {
  186.                 // Exponential spacing
  187.                 constexpr double b = 3;
  188.                 constexpr double dom = exp(b)-1;
  189.                 for (unsigned int i=0; i<dimSize_[d]; i++) {
  190.                     const double u = static_cast<double>(i+0.5) / dimSize_[d];
  191.                     const double x = ((exp(b*u) - 1)/dom);
  192.                     position_[d][i] = (x*range) + dimMin_[d];
  193.                 }
  194.             }
  195.             break;
  196.  
  197.             case ScalingType::Hyperbolic:
  198.             {
  199.                 // Hyperbolic tangent spacing
  200.                 constexpr double b = 3;
  201.                 constexpr double c = 0.5; // 1-c is the location of the 'attractor'
  202.                 constexpr double A = tanh(b*-1*c);
  203.                 constexpr double B = tanh(b*(1-c));
  204.                 for (unsigned int i=0; i<dimSize_[d]; i++) {
  205.                     const double u = static_cast<double>(i+0.5) / dimSize_[d];
  206.                     const double x = (tanh(b*(u-c)) - A) / (B-A);
  207.                     position_[d][i] = (x*range) + dimMin_[d];
  208.                 }
  209.             }
  210.             break;
  211.  
  212.             case ScalingType::Pivot:
  213.             {
  214.                 // Pivoted spacing - requires an even number of cells
  215.                 assert(dimSize_[d]%2==0);
  216.                 constexpr double p=0.35;
  217.                 constexpr double b = 2;
  218.                 constexpr double dom = exp(b)-1;
  219.                 const double mid = (range * p) + dimMin_[d];
  220.                 const size_t dim2 = dimSize_[d]/2;
  221.                 for (unsigned int i=0; i<dim2; i++) {
  222.                     const double u = static_cast<double>(i+0.5) / dim2;
  223.                     const double x = ((exp(b*u) - 1)/dom);
  224.                     const double step_lower = x*range*(p < 0.5 ? p : 1-p);
  225.                     const double step_upper = x*range*(p < 0.5 ? 1-p : p);
  226.                     // Upper half
  227.                     position_[d][i+dim2] = mid + step_upper;
  228.                     // Lower half
  229.                     position_[d][dim2-i-1] = mid - step_lower;
  230.                 }
  231.             }
  232.             break;
  233.  
  234.         }
  235. //        for (unsigned int i=0; i<dimSize_[d]; i++) {
  236. //            std::cout << position_[d][i] << " 1" << std::endl;
  237. //        }
  238.     }
  239. }
  240.  
  241. template<size_t meshDim>
  242. const double Mesh<meshDim>::getPosition(const size_t d, const int &i) const
  243. {
  244.     // Bounds checking if NDEBUG not defined
  245.     assert(i<dimSize_[d]);
  246.     assert(i>=0);
  247.     return position_[d][i];
  248. }
  249.  
  250. template<size_t meshDim>
  251. void Mesh<meshDim>::checkBounds(std::vector<int> idxs) const
  252. {
  253.     for (unsigned int d=0; d<meshDim; d++) {
  254.         assert(idxs[d] < dimSize_[d]);
  255.     }
  256. }
  257.  
  258. template<size_t meshDim>
  259. const size_t Mesh<meshDim>::numCells() const
  260. {
  261.     int numcells = 1;
  262.     for (unsigned int d=0; d<meshDim; d++) {
  263.         numcells *= dimSize_[d];
  264.     }
  265.     return numcells;
  266. }
  267.  
  268.  
  269. // ----------------------------------- FIELD.H ----------------------------//
  270. template <typename T, size_t fD, size_t mD>
  271. class Field
  272. {
  273. public:
  274.     // Constructor
  275.     Field(const Mesh<mD>& mesh, const std::string &fileName):
  276.     fileName_(fileName),
  277.     mesh_(mesh)
  278.     {
  279.         for (unsigned int d=0; d<fD; d++) {
  280.             field_[d].reserve(mesh.numCells());
  281.         }
  282.     }
  283.  
  284.     // Copy constructor
  285.     Field(const Field<T,fD,mD>& refToCopy):
  286.         field_(refToCopy.field_),
  287.         fileName_(refToCopy.fileName_),
  288.         mesh_(refToCopy.mesh_)
  289.     {}
  290.  
  291.     Field(const Field<T,fD,mD>& refToCopy, const std::string& name):
  292.         fileName_(name),
  293.         mesh_(refToCopy.mesh_)
  294.     {
  295.         for (size_t d=0; d<fD; d++) {
  296.             field_[d] = refToCopy.field_[d];
  297.             for (size_t i=0; i<mesh_.numCells(); i++) {
  298.                 std::cout << "field_[" << d << "][" << i << "] = " << field_[d][i] << std::endl;
  299.             }
  300.         }
  301.     }
  302.     // Set all values
  303.     void setZero();
  304.     void setFixed(const T &val);
  305.  
  306.     // Get the value at some point - now used only for 1D case.
  307.     template<size_t fd=fD, EnableIf<fd==1>..., typename... PosArgs, EnableIf<areT<size_t, PosArgs...>::value>...>
  308.     T& operator()(const PosArgs... indices) { return field_[0][singleIndex(0, indices...)]; }
  309.  
  310.     const Mesh<mD> &mesh() const { return mesh_; }
  311.  
  312.     // Return entire (const) vector<T>&
  313.     const std::vector<T> &x() const { return field_[0]; }
  314.     // Non-const
  315.     std::vector<T> &x() { return field_[0]; }
  316.  
  317.     // Only specific values
  318.     template<typename... PosArgs, EnableIf<areT<size_t, PosArgs...>::value>...>
  319.     T &x(const size_t xIdx, const PosArgs... indices);
  320. private:
  321.     std::vector<T> field_[fD];
  322.     std::string fileName_;
  323.     const Mesh<mD> &mesh_;
  324. public:
  325.     // Index lookup - component dependent.
  326.     template<typename... PosArgs, EnableIf<areT<size_t, PosArgs...>::value>...>
  327.     const size_t singleIndex(const size_t component, PosArgs... indices) const;
  328. };
  329. // #include "Field.hxx" // The functions implemented in this aren't used
  330. // Copying the singleIndex function for example
  331. template<typename T, size_t D, size_t M>
  332. template<typename... PosArgs, EnableIf<areT<size_t, PosArgs...>::value>...>
  333. const size_t Field<T,D,M>::singleIndex(const size_t component, PosArgs... indices) const
  334. {
  335.     //assert(component >= 0); // Commented out since MVCE doesn't #include <cassert>
  336.     //assert(component < M);
  337.     const size_t c0 = component;
  338.     const size_t c1 = (component+1)%3;
  339.     const size_t c2 = (component+2)%3;
  340.     const std::vector<size_t> &C = {mesh_.xCells(), mesh_.yCells(), mesh_.zCells()};
  341.     std::vector<size_t> idx{indices...};
  342.     idx.resize(3);
  343.     return (idx[c0] + (idx[c1]*C[c0]) + (idx[c2]*C[c1]*C[c0]));
  344. }
  345.  
  346. // ----------------------------- FIELD.cpp ------------------------------//
  347.  
  348. template<typename T, size_t fD, size_t mD>
  349. void Field<T,fD,mD>::setZero()
  350. {
  351.     setFixed(T(0));
  352. }
  353.  
  354. // This might need to be explicit.
  355. // Currently, 'a' can be converted for Field<double, D>
  356. template<typename T, size_t fD, size_t mD>
  357. void Field<T,fD,mD>::setFixed(const T &val)
  358. {
  359.     for (size_t d=0; d<fD; d++) {
  360.         std::fill(field_[d].begin(), field_[d].end(), val);
  361.     }
  362. }
  363.  
  364.  
  365. // ------------------------------- MAIN ----------------------------------//
  366. using vectorField = Field<double, 2, 2>;
  367. using scalarField = Field<double, 1, 2>;
  368.  
  369. int main()
  370. {
  371.     MeshDimension xDim(20, 0, 1);
  372.     MeshDimension yDim(30, 0, 1);
  373.     Mesh<2> mesh(xDim, yDim);
  374.  
  375.     vectorField U(mesh, "U");
  376.     U.setZero();
  377.     vectorField Phi(U, "copy-U");
  378.     std::cout << "Phi.x()[5] = " << Phi.x()[5] << std::endl;
  379. }
  380.  
  381.  
Runtime error #stdin #stdout 0s 3236KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Constructed with MeshDimension structs
meshDim structs
https://ideone.com/pzWCpU
language:
C++14 (gcc 8.3)
created:
10 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!