49 double out = sin( M_PI * x ) * cos( M_PI * y ) * sin( M_PI * z );
54 int main(
int argc,
char* argv[] )
56 int proc_id = 0,
size = 1;
58 MPI_Init( &argc, &argv );
59 MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
60 MPI_Comm_size( MPI_COMM_WORLD, &
size );
67 opts.
A = opts.
B = opts.
C = 1;
68 opts.
M = opts.
N = opts.
K = 1;
80 popts.
addOpt<
int >( string(
"blockSize,b" ), string(
"Block size of mesh (default=4)" ), &opts.
blockSize );
81 popts.
addOpt<
int >( string(
"xproc,M" ), string(
"Number of processors in x dir (default=1)" ), &opts.
M );
82 popts.
addOpt<
int >( string(
"yproc,N" ), string(
"Number of processors in y dir (default=1)" ), &opts.
N );
83 popts.
addOpt<
int >( string(
"zproc,K" ), string(
"Number of processors in z dir (default=1)" ), &opts.
K );
85 popts.
addOpt<
int >( string(
"xblocks,A" ), string(
"Number of blocks on a task in x dir (default=2)" ), &opts.
A );
86 popts.
addOpt<
int >( string(
"yblocks,B" ), string(
"Number of blocks on a task in y dir (default=2)" ), &opts.
B );
87 popts.
addOpt<
int >( string(
"zblocks,C" ), string(
"Number of blocks on a task in x dir (default=2)" ), &opts.
C );
89 popts.
addOpt<
double >( string(
"xsize,x" ), string(
"Total size in x direction (default=1.)" ), &opts.
xsize );
90 popts.
addOpt<
double >( string(
"ysize,y" ), string(
"Total size in y direction (default=1.)" ), &opts.
ysize );
91 popts.
addOpt<
double >( string(
"zsize,z" ), string(
"Total size in z direction (default=1.)" ), &opts.
zsize );
95 popts.
addOpt<
void >(
"quadratic,q",
"use hex 27 elements", &opts.
quadratic );
97 popts.
addOpt<
void >(
"keep_skins,k",
"keep skins with shared entities", &opts.
keep_skins );
99 popts.
addOpt<
void >(
"tetrahedrons,t",
"generate tetrahedrons", &opts.
tetra );
101 popts.
addOpt<
void >(
"faces_edges,f",
"create all faces and edges", &opts.
adjEnts );
103 popts.
addOpt<
int >( string(
"ghost_layers,g" ), string(
"Number of ghost layers (default=0)" ), &opts.
GL );
105 popts.
addOpt<
void >(
"parallel_merge,p",
"use parallel mesh merge, not vertex ID based merge", &opts.
parmerge );
109 double toler = 1.e-6;
110 popts.
addOpt<
double >( string(
"eps,e" ), string(
"tolerance for coupling, used in locating points" ), &toler );
112 bool writeMeshes =
false;
113 popts.
addOpt<
void >(
"print,p",
"write meshes", &writeMeshes );
127 double instance_time = MPI_Wtime();
128 double current = instance_time;
129 if( !proc_id ) std::cout <<
" instantiate first mesh " << instance_time -
start_time <<
"\n";
131 std::string interpTag(
"interp_tag" );
146 double fieldValue =
physField( vertPos[0], vertPos[1], vertPos[2] );
151 double setTag_time = MPI_Wtime();
152 if( !proc_id ) std::cout <<
" set tag " << setTag_time - current;
153 current = instance_time;
166 double instance_second = MPI_Wtime();
167 if( !proc_id ) std::cout <<
" instance second mesh" << instance_second - current <<
"\n";
168 current = instance_second;
173 rval =
mb->
write_file(
"mesh1.h5m", 0,
";;PARALLEL=WRITE_PART;CPUTIME;PARALLEL_COMM=0;", &fileset1, 1 );
MB_CHK_SET_ERR( rval,
"Can't write in parallel mesh 1" );
174 rval =
mb->
write_file(
"mesh2.h5m", 0,
";;PARALLEL=WRITE_PART;CPUTIME;PARALLEL_COMM=1;", &fileset2, 1 );
MB_CHK_SET_ERR( rval,
"Can't write in parallel mesh 1" );
175 double write_files = MPI_Wtime();
176 if( !proc_id ) std::cout <<
" write files " << write_files - current <<
"\n";
177 current = write_files;
183 double instancecoupler = MPI_Wtime();
184 if( !proc_id ) std::cout <<
" instance coupler " << instancecoupler - current <<
"\n";
185 current = instancecoupler;
190 std::vector< double > vpos;
191 int numPointsOfInterest = 0;
203 targ_verts =
subtract( targ_verts, tmp_verts );
205 numPointsOfInterest = (int)targ_verts.
size();
206 vpos.resize( 3 * targ_verts.
size() );
212 double locatetime = MPI_Wtime();
213 if( !proc_id ) std::cout <<
" locate points: " << locatetime - current <<
"\n";
214 current = locatetime;
217 std::vector< double > field( numPointsOfInterest );
223 for(
int i = 0; i < numPointsOfInterest; i++ )
225 double trval =
physField( vpos[3 * i], vpos[3 * i + 1], vpos[3 * i + 2] );
226 double err2 = fabs( trval - field[i] );
227 if( err2 > err_max ) err_max = err2;
230 double interpolateTime = MPI_Wtime();
231 if( !proc_id ) std::cout <<
" interpolate points: " << interpolateTime - current <<
"\n";
232 current = interpolateTime;
235 MPI_Allreduce( &err_max, &gerr, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
236 if( 0 == proc_id ) std::cout <<
"max err " << gerr <<
"\n";