1 #ifndef DGMSOLVER_HPP
2 #define DGMSOLVER_HPP
3 #include <vector>
4
5 namespace moab
6 {
7
8 class DGMSolver
9 {
10 DGMSolver(){};
11 ~DGMSolver(){};
12
13 public:
14 //! \brief compute combinational number, n choose k, maximum output is
15 //! std::numeric_limits<unsigned int>::max();
16 // If overflows, return 0
17 static unsigned int nchoosek( unsigned int n, unsigned int k );
18
19 //! \brief compute the number of columns for a multivariate vandermonde matrix, given certen
20 //! degree
21 /** If degree = 0, out put is 1;
22 *If kvars = 1, degree = k, output is k+1;
23 *If kvars = 2, degree = k, output is (k+2)*(k+1)/2;
24 */
25 static unsigned int compute_numcols_vander_multivar( unsigned int kvars, unsigned int degree );
26
27 //! \brief compute the monomial basis of mutiple variables, up to input degree,
28 //! lexicographically ordered
29 /** if degree = 0, output basis = {1}
30 *If kvars = 1, vars = {u}, degree = k, basis = {1,u,...,u^k}
31 *If kvars = 2, vars = {u,v}, degree = k, basis =
32 *{1,u,v,u^2,uv,v^2,u^3,u^2*v,uv^2,v^3,...,u^k,u^k-1*v,...,uv^k-1,v^k} If kvars = 3, vars =
33 *{u,v,w}, degree = k, basis =
34 *{1,u,v,w,u^2,uv,uw,v^2,v*w,w^2,...,u^k,u^k-1v,u^k-1w,...,v^k,v^k-1w,...,vw^k-1,w^k} \param
35 *kvars Integer, number of variables \param vars Pointer to array of doubles, size = kvars,
36 *variable values \param degree Integer, maximum degree \param basis Reference to vector, user
37 *input container to hold output which is appended to existing data; users don't have to
38 *preallocate for basis, this function will allocate interally
39 */
40 static void gen_multivar_monomial_basis( const int kvars,
41 const double* vars,
42 const int degree,
43 std::vector< double >& basis );
44
45 //! \brief compute multivariate vandermonde matrix, monomial basis ordered in the same way as
46 //! gen_multivar_monomial_basis
47 /** if degree = 0, V = {1,...,1}';
48 *If kvars = 1, us = {u1;u2;..,;um}, degree = k, V = {1 u1 u1^2 ... u1^k;1 u2 u2^2 ...
49 u2^k;...;1 um um^2 ... um^k}; *If kvars = 2, us = {u1 v1;u2 v2;...;um vm}, degree = k, V = {1
50 u1 v1 u1^2 u1v1 v1^2;...;1 um vm um^2 umvm vm^2};
51 * \param mrows Integer, number of points to evaluate Vandermonde matrix
52 * \param kvars Integer, number of variables
53 * \param us Pointer to array of doubles, size = mrow*kvars, variable values for all points.
54 Stored in row-wise, like {u1 v1 u2 v2 ...}
55 * \param degree Integer, maximum degree
56 * \param basis Reference to vector, user input container to hold Vandermonde matrix which is
57 appended to existing data; users don't have to preallocate for basis, this function will
58 allocate interally; the Vandermonde matrix is stored in an array, columnwise, like {1 ... 1
59 u1 ...um u1^2 ... um^2 ...}
60 */
61 static void gen_vander_multivar( const int mrows,
62 const int kvars,
63 const double* us,
64 const int degree,
65 std::vector< double >& V );
66
67 static void rescale_matrix( int mrows, int ncols, double* V, double* ts );
68
69 static void compute_qtransposeB( int mrows, int ncols, const double* Q, int bncols, double* bs );
70
71 static void qr_polyfit_safeguarded( const int mrows, const int ncols, double* V, double* D, int* rank );
72
73 static void backsolve( int mrows, int ncols, double* R, int bncols, double* bs, double* ws );
74
75 static void backsolve_polyfit_safeguarded( int dim,
76 int degree,
77 const bool interp,
78 int mrows,
79 int ncols,
80 double* R,
81 int bncols,
82 double* bs,
83 const double* ws,
84 int* degree_out );
85
86 static void vec_dotprod( const int len, const double* a, const double* b, double* c );
87
88 static void vec_scalarprod( const int len, const double* a, const double c, double* b );
89
90 static void vec_crossprod( const double a[3], const double b[3], double ( &c )[3] );
91
92 static double vec_innerprod( const int len, const double* a, const double* b );
93
94 static double vec_2norm( const int len, const double* a );
95
96 static double vec_normalize( const int len, const double* a, double* b );
97
98 static double vec_distance( const int len, const double* a, const double* b );
99
100 static void vec_projoff( const int len, const double* a, const double* b, double* c );
101
102 static void vec_linear_operation( const int len,
103 const double mu,
104 const double* a,
105 const double psi,
106 const double* b,
107 double* c );
108
109 static void get_tri_natural_coords( const int dim,
110 const double* cornercoords,
111 const int npts,
112 const double* currcoords,
113 double* naturalcoords );
114 };
115
116 } // namespace moab
117 #endif