6 #if defined( __MINGW32__ )
46 void compute_linear_coords(
const int nvpe,
double* elemcoords,
double* naturals,
double* linearcoords );
70 int main(
int argc,
char* argv[] )
73 MPI_Init( &argc, &argv );
79 int degree = 2,
dim = 0;
96 infile = std::string( argv[1] );
99 for(
int i = 2; i < argc; ++i )
103 if( std::string( argv[i] ) ==
"-degree" )
105 degree = atoi( argv[++i] );
107 else if( std::string( argv[i] ) ==
"-interp" )
109 interp = atoi( argv[++i] );
111 else if( std::string( argv[i] ) ==
"-dim" )
113 dim = atoi( argv[++i] );
118 std::cout << argv[i] << std::endl;
119 std::cout <<
"usage: " << argv[0]
120 <<
" <mesh file> -degree <degree> -interp <0=least square, "
121 "1=interpolation> -dim <mesh dimension>"
130 std::cout <<
"Dimension of input mesh should be provided, positive and less than 3" << std::endl;
134 if( degree <= 0 || dim > 2 ||
dim <= 0 )
136 std::cout <<
"Input degree should be positive number;\n";
137 std::cout <<
"Input dimesion should be positive and less than 3;" << std::endl;
141 std::cout <<
"Testing on " << infile <<
" with dimension " <<
dim <<
"\n";
142 std::string opts = interp ?
"interpolation" :
"least square fitting";
143 std::cout <<
"High order reconstruction with degree " << degree <<
" " << opts << std::endl;
165 MPI_Comm_size( comm, &nprocs );
166 MPI_Comm_rank( comm, &
rank );
184 read_options =
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
185 "SHARED_ENTS;PARALLEL_GHOST=2.0.";
189 read_options =
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_"
190 "SHARED_ENTS;PARALLEL_GHOST=1.0.";
194 MB_SET_ERR( MB_FAILURE,
"unsupported dimension" );
201 read_options =
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
212 assert( !pc && degree &&
dim );
228 double center[3] = { 0, 0, 0 };
229 double R = 1, r = 0.3;
248 double mxdist = 0, errl1 = 0, errl2 = 0, errli = 0;
249 double* pnts_proj =
new double[elems.
size() * 3];
256 double w = 1.0 / (double)nvpe;
257 std::vector< double > naturalcoords2fit( nvpe, w );
258 double newcoords[3], linearcoords[3];
260 pnts_proj[3 * ( *ielem - *elems.
begin() )] = newcoords[0];
261 pnts_proj[3 * ( *ielem - *elems.
begin() ) + 1] = newcoords[1];
262 pnts_proj[3 * ( *ielem - *elems.
begin() ) + 2] = newcoords[2];
263 std::vector< double > coords( 3 * nvpe );
272 std::cout <<
"Errors using exact torus for degree " << degree <<
" fit : L1 = " << errl1 <<
", L2 = " << errl2
273 <<
", Linf = " << errli << std::endl;
274 std::cout <<
"Maximum projection lift is " << mxdist << std::endl;
282 MB_SET_ERR( MB_FAILURE,
"n must be at least 2" );
286 std::vector< EntityHandle > verts( n * n );
287 size_t istr = tris.size();
288 tris.resize( istr + 2 * ( n - 1 ) * ( n - 1 ) );
289 double istep = 1.0 / (double)( n - 1 );
291 for(
size_t j = 0; j < n; ++j )
293 for(
size_t i = 0; i < n; ++i )
295 double coord[3] = { i * istep, j * istep, 0 };
298 verts[j * n + i] = temp;
302 for(
size_t jj = 0; jj < n - 1; ++jj )
304 for(
size_t ii = 0; ii < n - 1; ++ii )
306 EntityHandle conn[3] = { verts[jj * n + ii], verts[( jj + 1 ) * n + ii + 1], verts[( jj + 1 ) * n + ii] };
308 conn[0] = verts[jj * n + ii];
309 conn[1] = verts[jj * n + ii + 1];
310 conn[2] = verts[( jj + 1 ) * n + ii + 1];
322 MB_SET_ERR( MB_FAILURE,
"n must be at least 2" );
326 std::vector< EntityHandle > verts( n * n );
327 size_t istr = quads.size();
328 quads.resize( istr + ( n - 1 ) * ( n - 1 ) );
329 double istep = 1.0 / (double)( n - 1 );
331 for(
size_t j = 0; j < n; ++j )
333 for(
size_t i = 0; i < n; ++i )
335 double coord[3] = { i * istep, j * istep, 0 };
340 for(
size_t jj = 0; jj < n - 1; ++jj )
342 for(
size_t ii = 0; ii < n - 1; ++ii )
344 EntityHandle conn[4] = { verts[jj * n + ii], verts[jj * n + ii + 1], verts[( jj + 1 ) * n + ii + 1],
345 verts[( jj + 1 ) * n + ii] };
357 for(
size_t n = 2; n <= 2; ++n )
361 std::vector< EntityHandle > tris;
366 for(
int degree = 2; degree <= 2; ++degree )
373 for(
size_t itri = 0; itri < tris.size(); ++itri )
376 double naturalcoords2fit[nvpe] = { 1.0 / (double)nvpe, 1.0 / (
double)nvpe, 1.0 / (double)nvpe },
379 std::vector< EntityHandle > conn;
381 double coords[3 * nvpe], linearcoords[3];
387 std::cout <<
"triangulated unit square n= " << n <<
" degree= " << degree <<
" interpolation:\n";
388 std::cout <<
"maximum projection lift is " << mxdist << std::endl;
396 for(
size_t itri = 0; itri < tris.size(); ++itri )
399 double naturalcoords2fit[nvpe] = { 1.0 / (double)nvpe, 1.0 / (
double)nvpe, 1.0 / (double)nvpe },
402 std::vector< EntityHandle > conn;
404 double coords[3 * nvpe], linearcoords[3];
410 std::cout <<
"unit square n= " << n <<
" degree= " << degree <<
" least square:\n";
411 std::cout <<
"maximum projection lift is " << mxdist << std::endl;
420 assert( elemcoords && linearcoords );
422 for(
int i = 0; i < 3; ++i )
426 for(
int j = 0; j < nvpe; ++j )
428 linearcoords[i] += naturals[j] * elemcoords[3 * j + i];
437 for(
size_t n = 2; n <= 8; ++n )
441 std::vector< EntityHandle > quads;
446 for(
int degree = 1; degree <= 6; ++degree )
453 for(
size_t iquad = 0; iquad < quads.size(); ++iquad )
456 double w = 1.0 / (double)nvpe;
457 double naturalcoords2fit[nvpe] = { w, w, w, w }, newcoords[3];
459 std::vector< EntityHandle > conn;
461 double coords[3 * nvpe], linearcoords[3];
467 std::cout <<
"quadrilateral unit square n= " << n <<
" degree= " << degree <<
" interpolation:\n";
468 std::cout <<
"maximum projection lift is " << mxdist << std::endl;
474 for(
size_t iquad = 0; iquad < quads.size(); ++iquad )
477 double w = 1.0 / (double)nvpe;
478 double naturalcoords2fit[nvpe] = { w, w, w, w }, newcoords[3];
480 std::vector< EntityHandle > conn;
482 double coords[3 * nvpe], linearcoords[3];
488 std::cout <<
"quadrilateral unit square n= " << n <<
" degree= " << degree <<
" least square:\n";
489 std::cout <<
"maximum projection lift is " << mxdist << std::endl;
500 std::string filenames[] = { TestDir +
"unittest/sphere_tris_5.vtk", TestDir +
"unittest/sphere_tris_20.vtk",
501 TestDir +
"unittest/sphere_quads_5.vtk", TestDir +
"unittest/sphere_quads_20.vtk" };
505 for(
int ifile = 0; ifile < nfiles; ++ifile )
519 for(
int degree = 1; degree <= maxdeg; ++degree )
523 double mxdist = 0, mxerr = 0;
530 double w = 1.0 / (double)nvpe;
531 std::vector< double > naturalcoords2fit( nvpe, w );
532 double newcoords[3], linearcoords[3];
534 std::vector< double > coords( 3 * nvpe );
541 std::cout << filenames[ifile] <<
": unit sphere"
542 <<
" degree= " << degree <<
" interpolation:\n";
543 std::cout <<
"maximum projection lift is " << mxdist <<
", maximum error is " << mxerr << std::endl;
554 double w = 1.0 / (double)nvpe;
555 std::vector< double > naturalcoords2fit( nvpe, w );
556 double newcoords[3], linearcoords[3];
558 std::vector< double > coords( 3 * nvpe );
565 std::cout << filenames[ifile] <<
": unit sphere"
566 <<
" degree= " << degree <<
" least square:\n";
567 std::cout <<
"maximum projection lift is " << mxdist <<
", maximum error is " << mxerr << std::endl;
578 std::string filenames[] = { TestDir +
"unittest/circle_3.vtk", TestDir +
"unittest/circle_4.vtk",
579 TestDir +
"unittest/circle_10.vtk", TestDir +
"unittest/circle_20.vtk" };
583 for(
int ifile = 0; ifile < nfiles; ++ifile )
598 for(
int degree = 1; degree <= maxdeg; ++degree )
602 double mxdist = 0, mxerr = 0;
609 double w = 1.0 / (double)nvpe;
610 std::vector< double > naturalcoords2fit( nvpe, w );
611 double newcoords[3], linearcoords[3];
613 std::vector< double > coords( 3 * nvpe );
620 std::cout << filenames[ifile] <<
": unit circle"
621 <<
" degree= " << degree <<
" interpolation:\n";
622 std::cout <<
"maximum projection lift is " << mxdist <<
", maximum error is " << mxerr << std::endl;
633 double w = 1.0 / (double)nvpe;
634 std::vector< double > naturalcoords2fit( nvpe, w );
635 double newcoords[3], linearcoords[3];
637 std::vector< double > coords( 3 * nvpe );
644 std::cout << filenames[ifile] <<
": unit circle"
645 <<
" degree= " << degree <<
" least square:\n";
646 std::cout <<
"maximum projection lift is " << mxdist <<
", maximum error is " << mxerr << std::endl;
664 double pnts[3] = { 0, 0, 0 }, cnt[3] = { 0, 0, 0 }, nrms[3] = { 0, 0, 0 };
665 double x = 0, y = 0, z = 0, d1 = 0, d2 = 0;
667 for(
int i = 0; i < (int)verts.
size(); i++ )
674 d1 = sqrt( x * x + y * y );
678 d2 = sqrt( ( x - cnt[0] ) * ( x - cnt[0] ) + ( y - cnt[1] ) * ( y - cnt[1] ) + z * z );
679 nrms[0] = ( x - cnt[0] ) / d2;
680 nrms[1] = ( y - cnt[1] ) / d2;
682 pnts[0] = cnt[0] + r * nrms[0];
683 pnts[1] = cnt[1] + r * nrms[1];
684 pnts[2] = cnt[2] + r * nrms[2];
703 double x = 0, y = 0, z = 0;
704 double error_single = 0;
706 for(
int i = 0; i < npnts; i++ )
708 x = pnts[3 * i] -
center[0];
709 y = pnts[3 * i + 1] -
center[1];
710 z = pnts[3 * i + 2] -
center[2];
711 error_single = fabs( sqrt( pow( sqrt( pow( x, 2 ) + pow( y, 2 ) ) -
R, 2 ) + pow( z, 2 ) ) - r );
712 error_l1 = error_l1 + error_single;
713 error_l2 = error_l2 + error_single * error_single;
715 if( error_li < error_single )
717 error_li = error_single;
721 error_l1 = error_l1 / (double)npnts;
722 error_l2 = sqrt( error_l2 / (
double)npnts );