6 #if defined( __MINGW32__ )
63 MPI_Comm_size( comm, &procs );
67 read_options =
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
83 std::vector< EntityHandle > meshes;
95 int nvorig = vorig.
size();
96 double l1err = 0, l2err = 0, linferr = 0;
103 std::vector< EntityHandle > parentEntities;
105 assert( parentEntities.size() );
117 std::vector< double > cornercoords( 3 * nvpe ), currcoords( 3 );
121 std::vector< double > naturalcoords2fit( 3 );
132 linferr = std::max( linferr, err );
135 l1err /= verts.
size() - nvorig;
136 l2err = sqrt( l2err / ( verts.
size() - nvorig ) );
137 std::cout <<
"L1 error " << l1err <<
" L2 error " << l2err <<
" Linf error " << linferr << std::endl;
144 std::vector< int >& degs2fit,
161 MPI_Comm_size( comm, &procs );
165 read_options =
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;";
178 std::vector< EntityHandle > meshes;
179 std::vector< int > meshsizes;
181 std::vector< std::vector< double > > geoml1errs( 1 + degs2fit.size() ), geoml2errs( 1 + degs2fit.size() ),
182 geomlinferrs( 1 + degs2fit.size() );
186 for(
size_t i = 0; i < meshes.size(); ++i )
196 double currcoords[3], exactcoords[3];
205 meshsizes.push_back( elems.
size() );
207 std::vector< double > testpnts, testnaturalcoords;
209 testnaturalcoords.reserve( nvpe * elems.
size() *
nsamples );
214 std::vector< EntityHandle > conn;
217 std::vector< double > elemcoords( 3 * conn.size() );
225 double a = (double)rand() / RAND_MAX, b = (double)rand() / RAND_MAX, c = (double)rand() / RAND_MAX,
239 testpnts.push_back( a * elemcoords[0] + b * elemcoords[3] + c * elemcoords[6] );
240 testpnts.push_back( a * elemcoords[1] + b * elemcoords[4] + c * elemcoords[7] );
241 testpnts.push_back( a * elemcoords[2] + b * elemcoords[5] + c * elemcoords[8] );
242 testnaturalcoords.push_back( a );
243 testnaturalcoords.push_back( b );
244 testnaturalcoords.push_back( c );
248 double xi = (double)rand() / RAND_MAX, eta = (double)rand() / RAND_MAX, N[4];
251 N[0] = ( 1 - xi ) * ( 1 - eta ) / 4, N[1] = ( 1 + xi ) * ( 1 - eta ) / 4,
252 N[2] = ( 1 + xi ) * ( 1 + eta ) / 4, N[3] = ( 1 - xi ) * ( 1 + eta ) / 4;
253 testpnts.push_back( N[0] * elemcoords[0] + N[1] * elemcoords[3] + N[2] * elemcoords[6] +
254 N[3] * elemcoords[9] );
255 testpnts.push_back( N[0] * elemcoords[1] + N[1] * elemcoords[4] + N[2] * elemcoords[7] +
256 N[3] * elemcoords[10] );
257 testpnts.push_back( N[0] * elemcoords[2] + N[1] * elemcoords[5] + N[2] * elemcoords[8] +
258 N[3] * elemcoords[11] );
259 testnaturalcoords.push_back( N[0] );
260 testnaturalcoords.push_back( N[1] );
261 testnaturalcoords.push_back( N[2] );
262 testnaturalcoords.push_back( N[3] );
266 throw std::invalid_argument(
"NOT SUPPORTED ELEMENT TYPE" );
272 double l1err, l2err, linferr;
274 geoml1errs[0].push_back( l1err );
275 geoml2errs[0].push_back( l2err );
276 geomlinferrs[0].push_back( linferr );
280 for(
size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
290 &( testnaturalcoords[nvpe *
nsamples * index] ),
296 geoml1errs[ideg + 1].push_back( l1err );
297 geoml2errs[ideg + 1].push_back( l2err );
298 geomlinferrs[ideg + 1].push_back( linferr );
302 std::cout <<
"Mesh Size: ";
304 for(
size_t i = 0; i < meshsizes.size(); ++i )
306 std::cout << meshsizes[i] <<
" ";
309 std::cout << std::endl;
310 std::cout <<
"Degrees: 0 ";
312 for(
size_t ideg = 0; ideg < degs2fit.size(); ++ideg )
314 std::cout << degs2fit[ideg] <<
" ";
317 std::cout << std::endl;
318 std::cout <<
"L1-norm error: \n";
320 for(
size_t i = 0; i < geoml1errs.size(); ++i )
322 for(
size_t j = 0; j < geoml1errs[i].size(); ++j )
324 std::cout << geoml1errs[i][j] <<
" ";
327 std::cout << std::endl;
330 std::cout <<
"L2-norm error: \n";
332 for(
size_t i = 0; i < geoml2errs.size(); ++i )
334 for(
size_t j = 0; j < geoml2errs[i].size(); ++j )
336 std::cout << geoml2errs[i][j] <<
" ";
339 std::cout << std::endl;
342 std::cout <<
"Linf-norm error: \n";
344 for(
size_t i = 0; i < geomlinferrs.size(); ++i )
346 for(
size_t j = 0; j < geomlinferrs[i].size(); ++j )
348 std::cout << geomlinferrs[i][j] <<
" ";
351 std::cout << std::endl;
359 std::cout <<
"usage: ./uref_hirec_test <mesh file> -degree <degree> -interp <0=least square, "
360 "1=interpolation> -dim <mesh dimension> -geom <s=sphere,t=torus>"
362 std::cout <<
"Example: ./uref_hirec_test $MOAB_SOURCE_DIR/MeshFiles/unittest/sphere_tris_5.vtk "
363 "-degree 2 -interp 0 -dim 2 -geom s"
367 int main(
int argc,
char* argv[] )
370 MPI_Init( &argc, &argv );
376 std::string infile = TestDir +
"unittest/sphere_tris_5.vtk";
378 int degree = 2,
dim = 2,
geom = 0;
384 infile = std::string( argv[1] );
387 for(
int i = 2; i < argc; ++i )
391 if( std::string( argv[i] ) ==
"-degree" )
393 degree = atoi( argv[++i] );
395 else if( std::string( argv[i] ) ==
"-interp" )
397 interp = atoi( argv[++i] );
399 else if( std::string( argv[i] ) ==
"-dim" )
401 dim = atoi( argv[++i] );
404 else if( std::string( argv[i] ) ==
"-geom" )
406 geom = std::string( argv[++i] ) ==
"s" ? 0 : 1;
432 std::cout <<
"Dimension of input mesh should be provided, positive and less than 3" << std::endl;
436 std::cout <<
"Dimension of input mesh should be provided, positive and less than 3" << std::endl;
441 if( degree <= 0 || dim > 2 ||
dim <= 0 )
447 std::cout <<
"Input degree should be positive number;\n";
448 std::cout <<
"Input dimesion should be positive and less than 3;" << std::endl;
452 std::cout <<
"Input degree should be positive number;\n";
453 std::cout <<
"Input dimesion should be positive and less than 3;" << std::endl;
462 std::cout <<
"Testing on " << infile <<
" with dimension " <<
dim <<
"\n";
463 std::string opts = interp ?
"interpolation" :
"least square fitting";
464 std::cout <<
"High order reconstruction with degree " << degree <<
" " << opts << std::endl;
468 std::cout <<
"Testing on " << infile <<
" with dimension " <<
dim <<
"\n";
469 std::string opts = interp ?
"interpolation" :
"least square fitting";
470 std::cout <<
"High order reconstruction with degree " << degree <<
" " << opts << std::endl;
483 int level_degrees[3] = { 2, 2, 2 };
495 std::vector< int > degs2fit( 6 );
496 for(
int d = 1; d <= 6; ++d )