![]() |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
#include <DGMSolver.hpp>
Static Public Member Functions | |
static unsigned int | nchoosek (unsigned int n, unsigned int k) |
compute combinational number, n choose k, maximum output is std::numeric_limits<unsigned int>::max(); | |
static unsigned int | compute_numcols_vander_multivar (unsigned int kvars, unsigned int degree) |
compute the number of columns for a multivariate vandermonde matrix, given certen degree | |
static void | gen_multivar_monomial_basis (const int kvars, const double *vars, const int degree, std::vector< double > &basis) |
compute the monomial basis of mutiple variables, up to input degree, lexicographically ordered | |
static void | gen_vander_multivar (const int mrows, const int kvars, const double *us, const int degree, std::vector< double > &V) |
compute multivariate vandermonde matrix, monomial basis ordered in the same way as gen_multivar_monomial_basis | |
static void | rescale_matrix (int mrows, int ncols, double *V, double *ts) |
static void | compute_qtransposeB (int mrows, int ncols, const double *Q, int bncols, double *bs) |
static void | qr_polyfit_safeguarded (const int mrows, const int ncols, double *V, double *D, int *rank) |
static void | backsolve (int mrows, int ncols, double *R, int bncols, double *bs, double *ws) |
static void | backsolve_polyfit_safeguarded (int dim, int degree, const bool interp, int mrows, int ncols, double *R, int bncols, double *bs, const double *ws, int *degree_out) |
static void | vec_dotprod (const int len, const double *a, const double *b, double *c) |
static void | vec_scalarprod (const int len, const double *a, const double c, double *b) |
static void | vec_crossprod (const double a[3], const double b[3], double(&c)[3]) |
static double | vec_innerprod (const int len, const double *a, const double *b) |
static double | vec_2norm (const int len, const double *a) |
static double | vec_normalize (const int len, const double *a, double *b) |
static double | vec_distance (const int len, const double *a, const double *b) |
static void | vec_projoff (const int len, const double *a, const double *b, double *c) |
static void | vec_linear_operation (const int len, const double mu, const double *a, const double psi, const double *b, double *c) |
static void | get_tri_natural_coords (const int dim, const double *cornercoords, const int npts, const double *currcoords, double *naturalcoords) |
Private Member Functions | |
DGMSolver () | |
~DGMSolver () |
Definition at line 8 of file DGMSolver.hpp.
moab::DGMSolver::DGMSolver | ( | ) | [inline, private] |
Definition at line 10 of file DGMSolver.hpp.
{};
moab::DGMSolver::~DGMSolver | ( | ) | [inline, private] |
Definition at line 11 of file DGMSolver.hpp.
{};
void moab::DGMSolver::backsolve | ( | int | mrows, |
int | ncols, | ||
double * | R, | ||
int | bncols, | ||
double * | bs, | ||
double * | ws | ||
) | [static] |
Definition at line 256 of file DGMSolver.cpp.
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
{
#if 0
std::cout.precision(16);
std::cout<<"Before backsolve "<= 0; j-- )
{
for( int i = j + 1; i < ncols; ++i )
bs[mrows * k + j] = bs[mrows * k + j] - R[mrows * i + j] * bs[mrows * k + i];
assert( R[mrows * j + j] != 0 );
bs[mrows * k + j] = bs[mrows * k + j] / R[mrows * j + j];
}
}
for( int k = 0; k < bncols; k++ )
{
for( int j = 0; j < ncols; ++j )
bs[mrows * k + j] = bs[mrows * k + j] / ws[j];
}
}
void moab::DGMSolver::backsolve_polyfit_safeguarded | ( | int | dim, |
int | degree, | ||
const bool | interp, | ||
int | mrows, | ||
int | ncols, | ||
double * | R, | ||
int | bncols, | ||
double * | bs, | ||
const double * | ws, | ||
int * | degree_out | ||
) | [static] |
Definition at line 300 of file DGMSolver.cpp.
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
{
if( ncols < 1 ) std::cout << "ERROR: Invalid input to safeguarded polyfit backsolve routine.\n";
#if 0
std::cout.precision(12);
std::cout<<"Before backsolve "< bs_bak( numcols );
if( deg >= 2 )
{
for( int i = 0; i < numcols; i++ )
{
assert( mrows * k + i < mrows * bncols );
bs_bak.at( i ) = bs[mrows * k + i];
}
}
while( deg >= 1 )
{
int cend = numcols - 1;
bool downgrade = false;
// The reconstruction can be applied only on edges (2-d) or faces (3-d)
assert( cend >= 0 );
assert( dim > 0 && dim < 3 );
for( int d = deg; d >= 0; d-- )
{
int cstart = 0;
if( dim == 1 )
{
cstart = d;
}
else if( dim == 2 )
{
cstart = ( ( d + 1 ) * d ) / 2;
// cstart = ((d*(d+1))>>1)-interp;
}
// Solve for bs
for( int j = cend; j >= cstart; j-- )
{
assert( mrows * k + j < mrows * bncols );
for( int i = j + 1; i < numcols; ++i )
{
assert( mrows * k + i < mrows * bncols );
assert( mrows * i + j < mrows * ncols ); // check R
bs[mrows * k + j] = bs[mrows * k + j] - R[mrows * i + j] * bs[mrows * k + i];
}
assert( mrows * j + j < mrows * ncols ); // check R
bs[mrows * k + j] = bs[mrows * k + j] / R[mrows * j + j];
}
// Checking for change in the coefficient
if( d >= 2 && d < deg )
{
double tol;
if( dim == 1 )
{
tol = 1e-06;
assert( mrows * cstart + cstart < mrows * ncols ); // check R
double tb = bs_bak.at( cstart ) / R[mrows * cstart + cstart];
assert( mrows * k + cstart < mrows * bncols );
if( fabs( bs[mrows * k + cstart] - tb ) > ( 1 + tol ) * fabs( tb ) )
{
downgrade = true;
break;
}
}
else if( dim == 2 )
{
tol = 0.05;
std::vector< double > tb( cend - cstart + 1 );
for( int j = 0; j <= ( cend - cstart ); j++ )
{
tb.at( j ) = bs_bak.at( cstart + j );
}
for( int j = cend; j >= cstart; j-- )
{
int jind = j - cstart;
for( int i = j + 1; i <= cend; ++i )
{
assert( mrows * i + j < mrows * ncols ); // check R
tb.at( jind ) = tb.at( jind ) - R[mrows * i + j] * tb.at( i - cstart );
}
assert( mrows * j + j < mrows * ncols ); // check R
tb.at( jind ) = tb.at( jind ) / R[mrows * j + j];
assert( mrows * k + j < mrows * bncols );
double err = fabs( bs[mrows * k + j] - tb.at( jind ) );
if( ( err > tol ) && ( err >= ( 1 + tol ) * fabs( tb.at( jind ) ) ) )
{
downgrade = true;
break;
}
}
if( downgrade ) break;
}
}
cend = cstart - 1;
}
if( !downgrade )
break;
else
{
deg = deg - 1;
if( dim == 1 )
numcols = deg + 1;
else if( dim == 2 )
numcols = ( deg + 2 ) * ( deg + 1 ) / 2;
for( int i = 0; i < numcols; i++ )
{
assert( mrows * k + i < mrows * bncols );
bs[mrows * k + i] = bs_bak.at( i );
}
}
}
assert( k < bncols );
degree_out[k] = deg;
for( int i = 0; i < numcols; i++ )
{
// assert(mrows*k+i < mrows*bncols);
// assert(i < ncols);
bs[mrows * k + i] = bs[mrows * k + i] / ws[i];
}
for( int i = numcols; i < mrows; i++ )
{
// assert(mrows*k+i < mrows*bncols);
bs[mrows * k + i] = 0;
}
}
}
unsigned int moab::DGMSolver::compute_numcols_vander_multivar | ( | unsigned int | kvars, |
unsigned int | degree | ||
) | [static] |
compute the number of columns for a multivariate vandermonde matrix, given certen degree
If degree = 0, out put is 1; If kvars = 1, degree = k, output is k+1; If kvars = 2, degree = k, output is (k+2)*(k+1)/2;
Definition at line 42 of file DGMSolver.cpp.
References nchoosek().
Referenced by gen_multivar_monomial_basis(), and gen_vander_multivar().
{
unsigned int mcols = 0;
for( unsigned int i = 0; i <= degree; ++i )
{
unsigned int temp = nchoosek( kvars - 1 + i, kvars - 1 );
if( !temp )
{
std::cout << "overflow to compute nchoosek n= " << kvars - 1 + i << " k= " << kvars - 1 << std::endl;
return 0;
}
mcols += temp;
}
return mcols;
}
void moab::DGMSolver::compute_qtransposeB | ( | int | mrows, |
int | ncols, | ||
const double * | Q, | ||
int | bncols, | ||
double * | bs | ||
) | [static] |
Definition at line 174 of file DGMSolver.cpp.
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
{
for( int k = 0; k < ncols; k++ )
{
for( int j = 0; j < bncols; j++ )
{
double t2 = 0;
for( int i = k; i < mrows; i++ )
t2 += Q[mrows * k + i] * bs[mrows * j + i];
t2 = t2 + t2;
for( int i = k; i < mrows; i++ )
bs[mrows * j + i] -= t2 * Q[mrows * k + i];
}
}
}
void moab::DGMSolver::gen_multivar_monomial_basis | ( | const int | kvars, |
const double * | vars, | ||
const int | degree, | ||
std::vector< double > & | basis | ||
) | [static] |
compute the monomial basis of mutiple variables, up to input degree, lexicographically ordered
if degree = 0, output basis = {1} If kvars = 1, vars = {u}, degree = k, basis = {1,u,...,u^k} If kvars = 2, vars = {u,v}, degree = k, basis = {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 = {u,v,w}, degree = k, basis = {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}
kvars | Integer, number of variables |
vars | Pointer to array of doubles, size = kvars, variable values |
degree | Integer, maximum degree |
basis | Reference to vector, user input container to hold output which is appended to existing data; users don't have to preallocate for basis, this function will allocate interally |
Definition at line 58 of file DGMSolver.cpp.
References compute_numcols_vander_multivar().
{
unsigned int len = compute_numcols_vander_multivar( kvars, degree );
basis.reserve( len - basis.capacity() + basis.size() );
size_t iend = basis.size();
#ifndef NDEBUG
size_t istr = basis.size();
#endif
basis.push_back( 1 );
++iend;
if( !degree )
{
return;
}
std::vector< size_t > varspos( kvars );
// degree 1
for( int ivar = 0; ivar < kvars; ++ivar )
{
basis.push_back( vars[ivar] );
varspos[ivar] = iend++;
}
// degree 2 to degree
for( int ideg = 2; ideg <= degree; ++ideg )
{
size_t preend = iend;
for( int ivar = 0; ivar < kvars; ++ivar )
{
size_t varpreend = iend;
for( size_t ilast = varspos[ivar]; ilast < preend; ++ilast )
{
basis.push_back( vars[ivar] * basis[ilast] );
++iend;
}
varspos[ivar] = varpreend;
}
}
assert( len == iend - istr );
}
void moab::DGMSolver::gen_vander_multivar | ( | const int | mrows, |
const int | kvars, | ||
const double * | us, | ||
const int | degree, | ||
std::vector< double > & | V | ||
) | [static] |
compute multivariate vandermonde matrix, monomial basis ordered in the same way as gen_multivar_monomial_basis
if degree = 0, V = {1,...,1}'; If kvars = 1, us = {u1;u2;..,;um}, degree = k, V = {1 u1 u1^2 ... u1^k;1 u2 u2^2 ... u2^k;...;1 um um^2 ... um^k}; *If kvars = 2, us = {u1 v1;u2 v2;...;um vm}, degree = k, V = {1 u1 v1 u1^2 u1v1 v1^2;...;1 um vm um^2 umvm vm^2};
mrows | Integer, number of points to evaluate Vandermonde matrix |
kvars | Integer, number of variables |
us | Pointer to array of doubles, size = mrow*kvars, variable values for all points. Stored in row-wise, like {u1 v1 u2 v2 ...} |
degree | Integer, maximum degree |
basis | Reference to vector, user input container to hold Vandermonde matrix which is appended to existing data; users don't have to preallocate for basis, this function will allocate interally; the Vandermonde matrix is stored in an array, columnwise, like {1 ... 1 u1 ...um u1^2 ... um^2 ...} |
Definition at line 100 of file DGMSolver.cpp.
References compute_numcols_vander_multivar().
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
{
unsigned int ncols = compute_numcols_vander_multivar( kvars, degree );
V.reserve( mrows * ncols - V.capacity() + V.size() );
size_t istr = V.size(), icol = 0;
// add ones, V is stored in an single array, elements placed in columnwise order
for( int irow = 0; irow < mrows; ++irow )
{
V.push_back( 1 );
}
++icol;
if( !degree )
{
return;
}
std::vector< size_t > varspos( kvars );
// degree 1
for( int ivar = 0; ivar < kvars; ++ivar )
{
for( int irow = 0; irow < mrows; ++irow )
{
V.push_back( us[irow * kvars + ivar] ); // us stored in row-wise
}
varspos[ivar] = icol++;
}
// from 2 to degree
for( int ideg = 2; ideg <= degree; ++ideg )
{
size_t preendcol = icol;
for( int ivar = 0; ivar < kvars; ++ivar )
{
size_t varpreend = icol;
for( size_t ilast = varspos[ivar]; ilast < preendcol; ++ilast )
{
for( int irow = 0; irow < mrows; ++irow )
{
V.push_back( us[irow * kvars + ivar] * V[istr + irow + ilast * mrows] );
}
++icol;
}
varspos[ivar] = varpreend;
}
}
assert( icol == ncols );
}
void moab::DGMSolver::get_tri_natural_coords | ( | const int | dim, |
const double * | cornercoords, | ||
const int | npts, | ||
const double * | currcoords, | ||
double * | naturalcoords | ||
) | [static] |
Definition at line 667 of file DGMSolver.cpp.
References dim.
{
assert( dim == 2 || dim == 3 );
double a = 0, b = 0, d = 0, tol = 1e-12;
for( int i = 0; i < dim; ++i )
{
a += ( cornercoords[dim + i] - cornercoords[i] ) * ( cornercoords[dim + i] - cornercoords[i] );
b += ( cornercoords[dim + i] - cornercoords[i] ) * ( cornercoords[2 * dim + i] - cornercoords[i] );
d += ( cornercoords[2 * dim + i] - cornercoords[i] ) * ( cornercoords[2 * dim + i] - cornercoords[i] );
}
double det = a * d - b * b;
assert( det > 0 );
for( int ipt = 0; ipt < npts; ++ipt )
{
double e = 0, f = 0;
for( int i = 0; i < dim; ++i )
{
e += ( cornercoords[dim + i] - cornercoords[i] ) * ( currcoords[ipt * dim + i] - cornercoords[i] );
f += ( cornercoords[2 * dim + i] - cornercoords[i] ) * ( currcoords[ipt * dim + i] - cornercoords[i] );
}
naturalcoords[ipt * 3 + 1] = ( e * d - b * f ) / det;
naturalcoords[ipt * 3 + 2] = ( a * f - b * e ) / det;
naturalcoords[ipt * 3] = 1 - naturalcoords[ipt * 3 + 1] - naturalcoords[ipt * 3 + 2];
if( naturalcoords[ipt * 3] < -tol || naturalcoords[ipt * 3 + 1] < -tol || naturalcoords[ipt * 3 + 2] < -tol )
{
std::cout << "Corners: \n";
std::cout << cornercoords[0] << "\t" << cornercoords[1] << "\t" << cornercoords[3] << std::endl;
std::cout << cornercoords[3] << "\t" << cornercoords[4] << "\t" << cornercoords[5] << std::endl;
std::cout << cornercoords[6] << "\t" << cornercoords[7] << "\t" << cornercoords[8] << std::endl;
std::cout << "Candidate: \n";
std::cout << currcoords[ipt * dim] << "\t" << currcoords[ipt * dim + 1] << "\t" << currcoords[ipt * dim + 2]
<< std::endl;
exit( 0 );
}
assert( fabs( naturalcoords[ipt * 3] + naturalcoords[ipt * 3 + 1] + naturalcoords[ipt * 3 + 2] - 1 ) < tol );
for( int i = 0; i < dim; ++i )
{
assert( fabs( naturalcoords[ipt * 3] * cornercoords[i] +
naturalcoords[ipt * 3 + 1] * cornercoords[dim + i] +
naturalcoords[ipt * 3 + 2] * cornercoords[2 * dim + i] - currcoords[ipt * dim + i] ) < tol );
}
}
}
unsigned int moab::DGMSolver::nchoosek | ( | unsigned int | n, |
unsigned int | k | ||
) | [static] |
compute combinational number, n choose k, maximum output is std::numeric_limits<unsigned int>::max();
Definition at line 19 of file DGMSolver.cpp.
Referenced by compute_numcols_vander_multivar().
{
if( k > n )
{
return 0;
}
unsigned long long ans = 1;
if( k > ( n >> 1 ) )
{
k = n - k;
}
for( unsigned int i = 1; i <= k; ++i )
{
ans *= n--;
ans /= i;
if( ans > std::numeric_limits< unsigned int >::max() )
{
return 0;
}
}
return ans;
}
void moab::DGMSolver::qr_polyfit_safeguarded | ( | const int | mrows, |
const int | ncols, | ||
double * | V, | ||
double * | D, | ||
int * | rank | ||
) | [static] |
Definition at line 191 of file DGMSolver.cpp.
References t.
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
{
double tol = 1e-8;
*rank = ncols;
double* v = new double[mrows];
for( int k = 0; k < ncols; k++ )
{
int nv = mrows - k;
for( int j = 0; j < nv; j++ )
v[j] = V[mrows * k + ( j + k )];
double t2 = 0;
for( int j = 0; j < nv; j++ )
t2 = t2 + v[j] * v[j];
double t = sqrt( t2 );
double vnrm = 0;
if( v[0] >= 0 )
{
vnrm = sqrt( 2 * ( t2 + v[0] * t ) );
v[0] = v[0] + t;
}
else
{
vnrm = sqrt( 2 * ( t2 - v[0] * t ) );
v[0] = v[0] - t;
}
if( vnrm > 0 )
{
for( int j = 0; j < nv; j++ )
v[j] = v[j] / vnrm;
}
for( int j = k; j < ncols; j++ )
{
t2 = 0;
for( int i = 0; i < nv; i++ )
t2 = t2 + v[i] * V[mrows * j + ( i + k )];
t2 = t2 + t2;
for( int i = 0; i < nv; i++ )
V[mrows * j + ( i + k )] = V[mrows * j + ( i + k )] - t2 * v[i];
}
D[k] = V[mrows * k + k];
for( int i = 0; i < nv; i++ )
V[mrows * k + ( i + k )] = v[i];
if( ( fabs( D[k] ) ) < tol && ( *rank == ncols ) )
{
*rank = k;
break;
}
}
delete[] v;
}
void moab::DGMSolver::rescale_matrix | ( | int | mrows, |
int | ncols, | ||
double * | V, | ||
double * | ts | ||
) | [static] |
Definition at line 150 of file DGMSolver.cpp.
References vec_2norm().
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
{
// This function rescales the input matrix using the norm of each column.
double* v = new double[mrows];
for( int i = 0; i < ncols; i++ )
{
for( int j = 0; j < mrows; j++ )
v[j] = V[mrows * i + j];
// Compute norm of the column vector
double w = vec_2norm( mrows, v );
if( fabs( w ) == 0 )
ts[i] = 1;
else
{
ts[i] = w;
for( int j = 0; j < mrows; j++ )
V[mrows * i + j] = V[mrows * i + j] / ts[i];
}
}
delete[] v;
}
double moab::DGMSolver::vec_2norm | ( | const int | len, |
const double * | a | ||
) | [static] |
Definition at line 544 of file DGMSolver.cpp.
References MB_SET_ERR_RET_VAL.
Referenced by rescale_matrix(), and vec_projoff().
{
if( !a )
{
MB_SET_ERR_RET_VAL( "NULL Pointer", 0.0 );
}
double w = 0, s = 0;
for( int k = 0; k < len; k++ )
w = std::max( w, fabs( a[k] ) );
if( w == 0 )
{
return 0;
}
else
{
for( int k = 0; k < len; k++ )
{
s += ( a[k] / w ) * ( a[k] / w );
}
s = w * sqrt( s );
}
return s;
}
void moab::DGMSolver::vec_crossprod | ( | const double | a[3], |
const double | b[3], | ||
double(&) | c[3] | ||
) | [static] |
Definition at line 523 of file DGMSolver.cpp.
Referenced by moab::HiReconstruction::average_vertex_normal(), and moab::HiReconstruction::polyfit3d_surf_get_coeff().
{
c[0] = a[1] * b[2] - a[2] * b[1];
c[1] = a[2] * b[0] - a[0] * b[2];
c[2] = a[0] * b[1] - a[1] * b[0];
}
double moab::DGMSolver::vec_distance | ( | const int | len, |
const double * | a, | ||
const double * | b | ||
) | [static] |
Definition at line 604 of file DGMSolver.cpp.
{
double res = 0;
for( int i = 0; i < len; ++i )
{
res += ( a[i] - b[i] ) * ( a[i] - b[i] );
}
return sqrt( res );
}
void moab::DGMSolver::vec_dotprod | ( | const int | len, |
const double * | a, | ||
const double * | b, | ||
double * | c | ||
) | [static] |
Definition at line 499 of file DGMSolver.cpp.
References MB_SET_ERR_RET.
{
if( !a || !b || !c )
{
MB_SET_ERR_RET( "NULL Pointer" );
}
for( int i = 0; i < len; ++i )
{
c[i] = a[i] * b[i];
}
}
double moab::DGMSolver::vec_innerprod | ( | const int | len, |
const double * | a, | ||
const double * | b | ||
) | [static] |
Definition at line 530 of file DGMSolver.cpp.
References MB_SET_ERR_RET_VAL.
Referenced by moab::HiReconstruction::compute_weights(), moab::HiReconstruction::polyfit3d_curve_get_coeff(), moab::HiReconstruction::polyfit3d_surf_get_coeff(), vec_projoff(), moab::HiReconstruction::walf3d_curve_vertex_eval(), and moab::HiReconstruction::walf3d_surf_vertex_eval().
{
if( !a || !b )
{
MB_SET_ERR_RET_VAL( "NULL Pointer", 0.0 );
}
double ans = 0;
for( int i = 0; i < len; ++i )
{
ans += a[i] * b[i];
}
return ans;
}
void moab::DGMSolver::vec_linear_operation | ( | const int | len, |
const double | mu, | ||
const double * | a, | ||
const double | psi, | ||
const double * | b, | ||
double * | c | ||
) | [static] |
Definition at line 650 of file DGMSolver.cpp.
References MB_SET_ERR_RET.
Referenced by moab::HiReconstruction::average_vertex_normal(), moab::HiReconstruction::average_vertex_tangent(), moab::HiReconstruction::polyfit3d_curve_get_coeff(), moab::HiReconstruction::polyfit3d_surf_get_coeff(), and moab::HiReconstruction::walf3d_curve_vertex_eval().
{
if( !a || !b || !c )
{
MB_SET_ERR_RET( "NULL Pointer" );
}
for( int i = 0; i < len; ++i )
{
c[i] = mu * a[i] + psi * b[i];
}
}
double moab::DGMSolver::vec_normalize | ( | const int | len, |
const double * | a, | ||
double * | b | ||
) | [static] |
Definition at line 569 of file DGMSolver.cpp.
References MB_SET_ERR_RET_VAL.
Referenced by moab::HiReconstruction::average_vertex_normal(), moab::HiReconstruction::average_vertex_tangent(), and moab::HiReconstruction::polyfit3d_surf_get_coeff().
{
if( !a || !b )
{
MB_SET_ERR_RET_VAL( "NULL Pointer", 0.0 );
}
double nrm = 0, mx = 0;
for( int i = 0; i < len; ++i )
{
mx = std::max( fabs( a[i] ), mx );
}
if( mx == 0 )
{
for( int i = 0; i < len; ++i )
{
b[i] = 0;
}
return 0;
}
for( int i = 0; i < len; ++i )
{
nrm += ( a[i] / mx ) * ( a[i] / mx );
}
nrm = mx * sqrt( nrm );
if( nrm == 0 )
{
return nrm;
}
for( int i = 0; i < len; ++i )
{
b[i] = a[i] / nrm;
}
return nrm;
}
void moab::DGMSolver::vec_projoff | ( | const int | len, |
const double * | a, | ||
const double * | b, | ||
double * | c | ||
) | [static] |
Definition at line 614 of file DGMSolver.cpp.
References MB_SET_ERR_RET, vec_2norm(), and vec_innerprod().
Referenced by moab::HiReconstruction::polyfit3d_surf_get_coeff().
{
if( !a || !b || !c )
{
MB_SET_ERR_RET( "NULL Pointer" );
}
// c = a-b/;
double bnrm = vec_2norm( len, b );
if( bnrm == 0 )
{
for( int i = 0; i < len; ++i )
{
c[i] = a[i];
}
return;
}
double innerp = vec_innerprod( len, a, b ) / bnrm;
if( innerp == 0 )
{
if( c != a )
{
for( int i = 0; i < len; ++i )
{
c[i] = a[i];
}
}
return;
}
for( int i = 0; i < len; ++i )
{
c[i] = a[i] - innerp * b[i] / bnrm;
}
}
void moab::DGMSolver::vec_scalarprod | ( | const int | len, |
const double * | a, | ||
const double | c, | ||
double * | b | ||
) | [static] |
Definition at line 511 of file DGMSolver.cpp.
References MB_SET_ERR_RET.
{
if( !a || !b )
{
MB_SET_ERR_RET( "NULL Pointer" );
}
for( int i = 0; i < len; ++i )
{
b[i] = c * a[i];
}
}