1 /*
2 * This source code is part of
3 *
4 * HelFEM
5 * -
6 * Finite element methods for electronic structure calculations on small systems
7 *
8 * Written by Susi Lehtola, 2018-
9 * Copyright (c) 2018- Susi Lehtola
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 */
16 #include <helfem.h>
17
get_grid(double rmax,int num_el,int igrid,double zexp)18 arma::vec helfem::utils::get_grid(double rmax, int num_el, int igrid, double zexp) {
19 // Boundary values
20 arma::vec bval;
21
22 // Get boundary values
23 switch (igrid) {
24 // linear grid
25 case (1):
26 if (helfem::verbose)
27 printf("Using linear grid\n");
28 bval = arma::linspace<arma::vec>(0, rmax, num_el + 1);
29 break;
30
31 // quadratic grid (Schweizer et al 1999)
32 case (2):
33 if (helfem::verbose)
34 printf("Using quadratic grid\n");
35 bval.zeros(num_el + 1);
36 for (int i = 0; i <= num_el; i++)
37 bval(i) = i * i * rmax / (num_el * num_el);
38 break;
39
40 case (3):
41 // generalized polynomial grid, monotonic decrease till zexp~3, after that fails to work
42 if (helfem::verbose)
43 printf("Using generalized polynomial grid, zexp = %e\n", zexp);
44 bval.zeros(num_el + 1);
45 for (int i = 0; i <= num_el; i++)
46 bval(i) = rmax * std::pow(i * 1.0 / num_el, zexp);
47 break;
48
49 // generalized exponential grid, monotonic decrease till zexp~2, after that fails to work
50 case (4):
51 if (helfem::verbose)
52 printf("Using generalized exponential grid, zexp = %e\n", zexp);
53 bval = arma::exp(arma::pow(arma::linspace<arma::vec>(
54 0, std::pow(log(rmax + 1), 1.0 / zexp), num_el + 1),
55 zexp)) -
56 arma::ones<arma::vec>(num_el + 1);
57 break;
58
59 default:
60 throw std::logic_error("Invalid choice for grid\n");
61 }
62
63 // Make sure start and end points are numerically exact
64 bval(0) = 0.0;
65 bval(bval.n_elem - 1) = rmax;
66
67 return bval;
68 }
69