#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(); More... | |
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 More... | |
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 More... | |
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 More... | |
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.
|
inlineprivate |
Definition at line 10 of file DGMSolver.hpp.
10 {};
|
inlineprivate |
Definition at line 11 of file DGMSolver.hpp.
11 {};
|
static |
Definition at line 256 of file DGMSolver.cpp.
257 {
258 #if 0
259 std::cout.precision(16);
260 std::cout<<"Before backsolve "<<std::endl;
261 std::cout<<" V = "<<std::endl;
262 for (int k=0; k< ncols; k++){
263 for (int j=0; j<mrows; ++j){
264 std::cout<<R[mrows*k+j]<<std::endl;
265 }
266 std::cout<<std::endl;
267 }
268 std::cout<<std::endl;
269
270 //std::cout<<"#pnts = "<<mrows<<std::endl;
271 std::cout<<"bs = "<<std::endl;
272 std::cout<<std::endl;
273 for (int k=0; k< bncols; k++){
274 for (int j=0; j<mrows; ++j){
275 std::cout<<" "<<bs[mrows*k+j]<<std::endl;
276 }
277 }
278 std::cout<<std::endl;
279 #endif
280
281 for( int k = 0; k < bncols; k++ )
282 {
283 for( int j = ncols - 1; j >= 0; j-- )
284 {
285 for( int i = j + 1; i < ncols; ++i )
286 bs[mrows * k + j] = bs[mrows * k + j] - R[mrows * i + j] * bs[mrows * k + i];
287
288 assert( R[mrows * j + j] != 0 );
289 bs[mrows * k + j] = bs[mrows * k + j] / R[mrows * j + j];
290 }
291 }
292
293 for( int k = 0; k < bncols; k++ )
294 {
295 for( int j = 0; j < ncols; ++j )
296 bs[mrows * k + j] = bs[mrows * k + j] / ws[j];
297 }
298 }
References moab::R.
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
|
static |
Definition at line 300 of file DGMSolver.cpp.
310 {
311 if( ncols < 1 ) std::cout << "ERROR: Invalid input to safeguarded polyfit backsolve routine.\n";
312 #if 0
313 std::cout.precision(12);
314 std::cout<<"Before backsolve "<<std::endl;
315 std::cout<<" V = "<<std::endl;
316 for (int k=0; k< ncols; k++){
317 for (int j=0; j<mrows; ++j){
318 std::cout<<R[mrows*k+j]<<std::endl;
319 }
320 std::cout<<std::endl;
321 }
322 std::cout<<std::endl;
323
324 //std::cout<<"#pnts = "<<mrows<<std::endl;
325 std::cout<<"bs = "<<std::endl;
326 std::cout<<std::endl;
327 for (int k=0; k< bncols; k++){
328 for (int j=0; j<mrows; ++j){
329 std::cout<<" "<<bs[mrows*k+j]<<std::endl;
330 }
331 }
332 std::cout<<std::endl;
333
334 //std::cout<<" ] "<<std::endl;
335
336 std::cout<<"Input ws = [ ";
337 for (int k=0; k< ncols; k++){
338 std::cout<<ws[k]<<", ";
339 }
340 std::cout<<" ] "<<std::endl;
341
342 std::cout << "R: " << R << "size: [" << mrows << "," << ncols << "]" << std::endl;
343 std::cout << "bs: " << bs << "size: [" << mrows << "," << bncols << "]" << std::endl;
344 std::cout << "ws: " << ws << "size: [" << ncols << "," << 1 << "]" << std::endl;
345 std::cout << "degree_out: " << degree_out << std::endl;
346 #endif
347
348 int deg, numcols = 0;
349
350 for( int k = 0; k < bncols; k++ )
351 {
352 deg = degree;
353 /* I think we should consider interp = true/false -Xinglin*/
354 if( dim == 1 )
355 numcols = deg + 1 - interp;
356 else if( dim == 2 )
357 numcols = ( deg + 2 ) * ( deg + 1 ) / 2 - interp;
358
359 assert( numcols <= ncols );
360
361 std::vector< double > bs_bak( numcols );
362
363 if( deg >= 2 )
364 {
365 for( int i = 0; i < numcols; i++ )
366 {
367 assert( mrows * k + i < mrows * bncols );
368 bs_bak.at( i ) = bs[mrows * k + i];
369 }
370 }
371
372 while( deg >= 1 )
373 {
374 int cend = numcols - 1;
375 bool downgrade = false;
376
377 // The reconstruction can be applied only on edges (2-d) or faces (3-d)
378 assert( cend >= 0 );
379 assert( dim > 0 && dim < 3 );
380
381 for( int d = deg; d >= 0; d-- )
382 {
383 int cstart = 0;
384 if( dim == 1 )
385 {
386 cstart = d;
387 }
388 else if( dim == 2 )
389 {
390 cstart = ( ( d + 1 ) * d ) / 2;
391 // cstart = ((d*(d+1))>>1)-interp;
392 }
393
394 // Solve for bs
395 for( int j = cend; j >= cstart; j-- )
396 {
397 assert( mrows * k + j < mrows * bncols );
398 for( int i = j + 1; i < numcols; ++i )
399 {
400 assert( mrows * k + i < mrows * bncols );
401 assert( mrows * i + j < mrows * ncols ); // check R
402 bs[mrows * k + j] = bs[mrows * k + j] - R[mrows * i + j] * bs[mrows * k + i];
403 }
404 assert( mrows * j + j < mrows * ncols ); // check R
405 bs[mrows * k + j] = bs[mrows * k + j] / R[mrows * j + j];
406 }
407
408 // Checking for change in the coefficient
409 if( d >= 2 && d < deg )
410 {
411 double tol;
412
413 if( dim == 1 )
414 {
415 tol = 1e-06;
416 assert( mrows * cstart + cstart < mrows * ncols ); // check R
417 double tb = bs_bak.at( cstart ) / R[mrows * cstart + cstart];
418 assert( mrows * k + cstart < mrows * bncols );
419 if( fabs( bs[mrows * k + cstart] - tb ) > ( 1 + tol ) * fabs( tb ) )
420 {
421 downgrade = true;
422 break;
423 }
424 }
425
426 else if( dim == 2 )
427 {
428 tol = 0.05;
429
430 std::vector< double > tb( cend - cstart + 1 );
431 for( int j = 0; j <= ( cend - cstart ); j++ )
432 {
433 tb.at( j ) = bs_bak.at( cstart + j );
434 }
435
436 for( int j = cend; j >= cstart; j-- )
437 {
438 int jind = j - cstart;
439
440 for( int i = j + 1; i <= cend; ++i )
441 {
442 assert( mrows * i + j < mrows * ncols ); // check R
443 tb.at( jind ) = tb.at( jind ) - R[mrows * i + j] * tb.at( i - cstart );
444 }
445 assert( mrows * j + j < mrows * ncols ); // check R
446 tb.at( jind ) = tb.at( jind ) / R[mrows * j + j];
447 assert( mrows * k + j < mrows * bncols );
448 double err = fabs( bs[mrows * k + j] - tb.at( jind ) );
449
450 if( ( err > tol ) && ( err >= ( 1 + tol ) * fabs( tb.at( jind ) ) ) )
451 {
452 downgrade = true;
453 break;
454 }
455 }
456
457 if( downgrade ) break;
458 }
459 }
460
461 cend = cstart - 1;
462 }
463
464 if( !downgrade )
465 break;
466 else
467 {
468 deg = deg - 1;
469 if( dim == 1 )
470 numcols = deg + 1;
471 else if( dim == 2 )
472 numcols = ( deg + 2 ) * ( deg + 1 ) / 2;
473
474 for( int i = 0; i < numcols; i++ )
475 {
476 assert( mrows * k + i < mrows * bncols );
477 bs[mrows * k + i] = bs_bak.at( i );
478 }
479 }
480 }
481 assert( k < bncols );
482 degree_out[k] = deg;
483
484 for( int i = 0; i < numcols; i++ )
485 {
486 // assert(mrows*k+i < mrows*bncols);
487 // assert(i < ncols);
488 bs[mrows * k + i] = bs[mrows * k + i] / ws[i];
489 }
490
491 for( int i = numcols; i < mrows; i++ )
492 {
493 // assert(mrows*k+i < mrows*bncols);
494 bs[mrows * k + i] = 0;
495 }
496 }
497 }
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
|
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.
43 {
44 unsigned int mcols = 0;
45 for( unsigned int i = 0; i <= degree; ++i )
46 {
47 unsigned int temp = nchoosek( kvars - 1 + i, kvars - 1 );
48 if( !temp )
49 {
50 std::cout << "overflow to compute nchoosek n= " << kvars - 1 + i << " k= " << kvars - 1 << std::endl;
51 return 0;
52 }
53 mcols += temp;
54 }
55 return mcols;
56 }
References nchoosek().
Referenced by gen_multivar_monomial_basis(), and gen_vander_multivar().
|
static |
Definition at line 174 of file DGMSolver.cpp.
175 {
176 for( int k = 0; k < ncols; k++ )
177 {
178 for( int j = 0; j < bncols; j++ )
179 {
180 double t2 = 0;
181 for( int i = k; i < mrows; i++ )
182 t2 += Q[mrows * k + i] * bs[mrows * j + i];
183 t2 = t2 + t2;
184
185 for( int i = k; i < mrows; i++ )
186 bs[mrows * j + i] -= t2 * Q[mrows * k + i];
187 }
188 }
189 }
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
|
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.
62 {
63 unsigned int len = compute_numcols_vander_multivar( kvars, degree );
64 basis.reserve( len - basis.capacity() + basis.size() );
65 size_t iend = basis.size();
66 #ifndef NDEBUG
67 size_t istr = basis.size();
68 #endif
69 basis.push_back( 1 );
70 ++iend;
71 if( !degree )
72 {
73 return;
74 }
75 std::vector< size_t > varspos( kvars );
76 // degree 1
77 for( int ivar = 0; ivar < kvars; ++ivar )
78 {
79 basis.push_back( vars[ivar] );
80 varspos[ivar] = iend++;
81 }
82 // degree 2 to degree
83 for( int ideg = 2; ideg <= degree; ++ideg )
84 {
85 size_t preend = iend;
86 for( int ivar = 0; ivar < kvars; ++ivar )
87 {
88 size_t varpreend = iend;
89 for( size_t ilast = varspos[ivar]; ilast < preend; ++ilast )
90 {
91 basis.push_back( vars[ivar] * basis[ilast] );
92 ++iend;
93 }
94 varspos[ivar] = varpreend;
95 }
96 }
97 assert( len == iend - istr );
98 }
References compute_numcols_vander_multivar().
|
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.
105 {
106 unsigned int ncols = compute_numcols_vander_multivar( kvars, degree );
107 V.reserve( mrows * ncols - V.capacity() + V.size() );
108 size_t istr = V.size(), icol = 0;
109 // add ones, V is stored in an single array, elements placed in columnwise order
110 for( int irow = 0; irow < mrows; ++irow )
111 {
112 V.push_back( 1 );
113 }
114 ++icol;
115 if( !degree )
116 {
117 return;
118 }
119 std::vector< size_t > varspos( kvars );
120 // degree 1
121 for( int ivar = 0; ivar < kvars; ++ivar )
122 {
123 for( int irow = 0; irow < mrows; ++irow )
124 {
125 V.push_back( us[irow * kvars + ivar] ); // us stored in row-wise
126 }
127 varspos[ivar] = icol++;
128 }
129 // from 2 to degree
130 for( int ideg = 2; ideg <= degree; ++ideg )
131 {
132 size_t preendcol = icol;
133 for( int ivar = 0; ivar < kvars; ++ivar )
134 {
135 size_t varpreend = icol;
136 for( size_t ilast = varspos[ivar]; ilast < preendcol; ++ilast )
137 {
138 for( int irow = 0; irow < mrows; ++irow )
139 {
140 V.push_back( us[irow * kvars + ivar] * V[istr + irow + ilast * mrows] );
141 }
142 ++icol;
143 }
144 varspos[ivar] = varpreend;
145 }
146 }
147 assert( icol == ncols );
148 }
References compute_numcols_vander_multivar().
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
|
static |
Definition at line 667 of file DGMSolver.cpp.
672 {
673 assert( dim == 2 || dim == 3 );
674 double a = 0, b = 0, d = 0, tol = 1e-12;
675 for( int i = 0; i < dim; ++i )
676 {
677 a += ( cornercoords[dim + i] - cornercoords[i] ) * ( cornercoords[dim + i] - cornercoords[i] );
678 b += ( cornercoords[dim + i] - cornercoords[i] ) * ( cornercoords[2 * dim + i] - cornercoords[i] );
679 d += ( cornercoords[2 * dim + i] - cornercoords[i] ) * ( cornercoords[2 * dim + i] - cornercoords[i] );
680 }
681 double det = a * d - b * b;
682 assert( det > 0 );
683 for( int ipt = 0; ipt < npts; ++ipt )
684 {
685 double e = 0, f = 0;
686 for( int i = 0; i < dim; ++i )
687 {
688 e += ( cornercoords[dim + i] - cornercoords[i] ) * ( currcoords[ipt * dim + i] - cornercoords[i] );
689 f += ( cornercoords[2 * dim + i] - cornercoords[i] ) * ( currcoords[ipt * dim + i] - cornercoords[i] );
690 }
691 naturalcoords[ipt * 3 + 1] = ( e * d - b * f ) / det;
692 naturalcoords[ipt * 3 + 2] = ( a * f - b * e ) / det;
693 naturalcoords[ipt * 3] = 1 - naturalcoords[ipt * 3 + 1] - naturalcoords[ipt * 3 + 2];
694 if( naturalcoords[ipt * 3] < -tol || naturalcoords[ipt * 3 + 1] < -tol || naturalcoords[ipt * 3 + 2] < -tol )
695 {
696 std::cout << "Corners: \n";
697 std::cout << cornercoords[0] << "\t" << cornercoords[1] << "\t" << cornercoords[3] << std::endl;
698 std::cout << cornercoords[3] << "\t" << cornercoords[4] << "\t" << cornercoords[5] << std::endl;
699 std::cout << cornercoords[6] << "\t" << cornercoords[7] << "\t" << cornercoords[8] << std::endl;
700 std::cout << "Candidate: \n";
701 std::cout << currcoords[ipt * dim] << "\t" << currcoords[ipt * dim + 1] << "\t" << currcoords[ipt * dim + 2]
702 << std::endl;
703 exit( 0 );
704 }
705 assert( fabs( naturalcoords[ipt * 3] + naturalcoords[ipt * 3 + 1] + naturalcoords[ipt * 3 + 2] - 1 ) < tol );
706 for( int i = 0; i < dim; ++i )
707 {
708 assert( fabs( naturalcoords[ipt * 3] * cornercoords[i] +
709 naturalcoords[ipt * 3 + 1] * cornercoords[dim + i] +
710 naturalcoords[ipt * 3 + 2] * cornercoords[2 * dim + i] - currcoords[ipt * dim + i] ) < tol );
711 }
712 }
713 }
References dim.
|
static |
compute combinational number, n choose k, maximum output is std::numeric_limits<unsigned int>::max();
Definition at line 19 of file DGMSolver.cpp.
20 {
21 if( k > n )
22 {
23 return 0;
24 }
25 unsigned long long ans = 1;
26 if( k > ( n >> 1 ) )
27 {
28 k = n - k;
29 }
30 for( unsigned int i = 1; i <= k; ++i )
31 {
32 ans *= n--;
33 ans /= i;
34 if( ans > std::numeric_limits< unsigned int >::max() )
35 {
36 return 0;
37 }
38 }
39 return ans;
40 }
Referenced by compute_numcols_vander_multivar().
|
static |
Definition at line 191 of file DGMSolver.cpp.
192 {
193 double tol = 1e-8;
194 *rank = ncols;
195 double* v = new double[mrows];
196
197 for( int k = 0; k < ncols; k++ )
198 {
199 int nv = mrows - k;
200
201 for( int j = 0; j < nv; j++ )
202 v[j] = V[mrows * k + ( j + k )];
203
204 double t2 = 0;
205
206 for( int j = 0; j < nv; j++ )
207 t2 = t2 + v[j] * v[j];
208
209 double t = sqrt( t2 );
210 double vnrm = 0;
211
212 if( v[0] >= 0 )
213 {
214 vnrm = sqrt( 2 * ( t2 + v[0] * t ) );
215 v[0] = v[0] + t;
216 }
217 else
218 {
219 vnrm = sqrt( 2 * ( t2 - v[0] * t ) );
220 v[0] = v[0] - t;
221 }
222
223 if( vnrm > 0 )
224 {
225 for( int j = 0; j < nv; j++ )
226 v[j] = v[j] / vnrm;
227 }
228
229 for( int j = k; j < ncols; j++ )
230 {
231 t2 = 0;
232 for( int i = 0; i < nv; i++ )
233 t2 = t2 + v[i] * V[mrows * j + ( i + k )];
234
235 t2 = t2 + t2;
236
237 for( int i = 0; i < nv; i++ )
238 V[mrows * j + ( i + k )] = V[mrows * j + ( i + k )] - t2 * v[i];
239 }
240
241 D[k] = V[mrows * k + k];
242
243 for( int i = 0; i < nv; i++ )
244 V[mrows * k + ( i + k )] = v[i];
245
246 if( ( fabs( D[k] ) ) < tol && ( *rank == ncols ) )
247 {
248 *rank = k;
249 break;
250 }
251 }
252
253 delete[] v;
254 }
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
|
static |
Definition at line 150 of file DGMSolver.cpp.
151 {
152 // This function rescales the input matrix using the norm of each column.
153 double* v = new double[mrows];
154 for( int i = 0; i < ncols; i++ )
155 {
156 for( int j = 0; j < mrows; j++ )
157 v[j] = V[mrows * i + j];
158
159 // Compute norm of the column vector
160 double w = vec_2norm( mrows, v );
161
162 if( fabs( w ) == 0 )
163 ts[i] = 1;
164 else
165 {
166 ts[i] = w;
167 for( int j = 0; j < mrows; j++ )
168 V[mrows * i + j] = V[mrows * i + j] / ts[i];
169 }
170 }
171 delete[] v;
172 }
References vec_2norm().
Referenced by moab::HiReconstruction::eval_vander_bivar_cmf(), and moab::HiReconstruction::eval_vander_univar_cmf().
|
static |
Definition at line 544 of file DGMSolver.cpp.
545 {
546 if( !a )
547 {
548 MB_SET_ERR_RET_VAL( "NULL Pointer", 0.0 );
549 }
550 double w = 0, s = 0;
551 for( int k = 0; k < len; k++ )
552 w = std::max( w, fabs( a[k] ) );
553
554 if( w == 0 )
555 {
556 return 0;
557 }
558 else
559 {
560 for( int k = 0; k < len; k++ )
561 {
562 s += ( a[k] / w ) * ( a[k] / w );
563 }
564 s = w * sqrt( s );
565 }
566 return s;
567 }
References MB_SET_ERR_RET_VAL.
Referenced by rescale_matrix(), and vec_projoff().
|
static |
Definition at line 523 of file DGMSolver.cpp.
524 { 525 c[0] = a[1] * b[2] - a[2] * b[1]; 526 c[1] = a[2] * b[0] - a[0] * b[2]; 527 c[2] = a[0] * b[1] - a[1] * b[0]; 528 }
Referenced by moab::HiReconstruction::average_vertex_normal(), and moab::HiReconstruction::polyfit3d_surf_get_coeff().
|
static |
Definition at line 604 of file DGMSolver.cpp.
605 {
606 double res = 0;
607 for( int i = 0; i < len; ++i )
608 {
609 res += ( a[i] - b[i] ) * ( a[i] - b[i] );
610 }
611 return sqrt( res );
612 }
|
static |
Definition at line 499 of file DGMSolver.cpp.
500 {
501 if( !a || !b || !c )
502 {
503 MB_SET_ERR_RET( "NULL Pointer" );
504 }
505 for( int i = 0; i < len; ++i )
506 {
507 c[i] = a[i] * b[i];
508 }
509 }
References MB_SET_ERR_RET.
|
static |
Definition at line 530 of file DGMSolver.cpp.
531 {
532 if( !a || !b )
533 {
534 MB_SET_ERR_RET_VAL( "NULL Pointer", 0.0 );
535 }
536 double ans = 0;
537 for( int i = 0; i < len; ++i )
538 {
539 ans += a[i] * b[i];
540 }
541 return ans;
542 }
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().
|
static |
Definition at line 650 of file DGMSolver.cpp.
656 {
657 if( !a || !b || !c )
658 {
659 MB_SET_ERR_RET( "NULL Pointer" );
660 }
661 for( int i = 0; i < len; ++i )
662 {
663 c[i] = mu * a[i] + psi * b[i];
664 }
665 }
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().
|
static |
Definition at line 569 of file DGMSolver.cpp.
570 {
571 if( !a || !b )
572 {
573 MB_SET_ERR_RET_VAL( "NULL Pointer", 0.0 );
574 }
575 double nrm = 0, mx = 0;
576 for( int i = 0; i < len; ++i )
577 {
578 mx = std::max( fabs( a[i] ), mx );
579 }
580 if( mx == 0 )
581 {
582 for( int i = 0; i < len; ++i )
583 {
584 b[i] = 0;
585 }
586 return 0;
587 }
588 for( int i = 0; i < len; ++i )
589 {
590 nrm += ( a[i] / mx ) * ( a[i] / mx );
591 }
592 nrm = mx * sqrt( nrm );
593 if( nrm == 0 )
594 {
595 return nrm;
596 }
597 for( int i = 0; i < len; ++i )
598 {
599 b[i] = a[i] / nrm;
600 }
601 return nrm;
602 }
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().
|
static |
Definition at line 614 of file DGMSolver.cpp.
615 {
616 if( !a || !b || !c )
617 {
618 MB_SET_ERR_RET( "NULL Pointer" );
619 }
620 // c = a-<a,b>b/<b,b>;
621 double bnrm = vec_2norm( len, b );
622 if( bnrm == 0 )
623 {
624 for( int i = 0; i < len; ++i )
625 {
626 c[i] = a[i];
627 }
628 return;
629 }
630 double innerp = vec_innerprod( len, a, b ) / bnrm;
631
632 if( innerp == 0 )
633 {
634 if( c != a )
635 {
636 for( int i = 0; i < len; ++i )
637 {
638 c[i] = a[i];
639 }
640 }
641 return;
642 }
643
644 for( int i = 0; i < len; ++i )
645 {
646 c[i] = a[i] - innerp * b[i] / bnrm;
647 }
648 }
References MB_SET_ERR_RET, vec_2norm(), and vec_innerprod().
Referenced by moab::HiReconstruction::polyfit3d_surf_get_coeff().
|
static |
Definition at line 511 of file DGMSolver.cpp.
512 {
513 if( !a || !b )
514 {
515 MB_SET_ERR_RET( "NULL Pointer" );
516 }
517 for( int i = 0; i < len; ++i )
518 {
519 b[i] = c * a[i];
520 }
521 }
References MB_SET_ERR_RET.