Homec4science
Diffusion CADD_mesher (dev/incorporate-qhull) Browse

rCADDMESH/test/mesh_test_2d.cc
58a104d7bd95dev/incorporate-qhull

Last Change
Raw File

mesh_test_2d.cc
View Options

#include "common.hh"
#include "cadd_mesh.hh"
#include "block.hh"
#include <stdlib.h>
#include <time.h>
#include <vector>
#include <iostream>
#include <cmath>
#include <limits>
#include <sys/time.h>
timeval start_time;
inline void start() {
gettimeofday(&start_time, NULL);
}
inline Real report() {
timeval now;
gettimeofday(&now, NULL);
Real diff = (now.tv_sec - start_time.tv_sec);
diff += 1e-6 * (now.tv_usec - start_time.tv_usec);
return diff;
}
int main(int argc, char *argv[])
{
const Uint dim = twoD;
Real total_elapsed = 0;
Real tol = 1e-4;
std::vector<Real> points;
std::vector<Real> nodal_vals;
Real lim = 32;
for (Uint i = 0; i < lim+1 ; ++i) {
for (Uint j = 0; j < lim+1 ; ++j) {
points.push_back(i/lim);
points.push_back(j/lim);
nodal_vals.push_back((i+j)/lim);
}
}
Real min_max[2*dim] = {.45, 10.6,
.45, .55};
Block<dim> exclude(min_max, "excludes (1,1,1) node)");
Real min_max2[2*dim] = {.4, 10.6,
.4, .6};
Block<dim> atomic(min_max2, "excludes (1,1,1) node)");
start();
Mesh<dim> mesh(points, 0, &exclude);
Real constructor_elapsed = report();
total_elapsed += constructor_elapsed;
start();
mesh.delaunay(1000);
Real delaunay_elapsed = report();
total_elapsed += delaunay_elapsed;
start();
mesh.cleanup(atomic);
Real cleanup_elapsed = report();
total_elapsed += cleanup_elapsed;
mesh.dump_paraview("after_cleanup_2D");
std::cout << mesh << std::endl;
start();
mesh.initInterpolation();
Real init_interpolation_elapsed = report();
total_elapsed += init_interpolation_elapsed;
Real coord[dim] = {0.334737, 0.495676};//{.76, .231, .2};
Real interpol_hihi;
Real interpol_elapsed;
try {
start();
interpol_hihi = mesh.interpolate(coord, nodal_vals);
interpol_elapsed = report();
total_elapsed += interpol_elapsed;
} catch (std::string& scheisse) {
std::cout << scheisse <<std::endl;
throw(scheisse);
}
Real interpol_correct = 0;
for (Uint i = 0; i < dim ; ++i) {
interpol_correct += coord[i];
}
std::cout << " Interpolated value = " << interpol_hihi;
if (fabs(interpol_hihi - interpol_correct)/interpol_correct <= 2*std::numeric_limits<Real>::epsilon()) {
std::cout << ". This is correct" << std::endl;
} else {
std::cout << ". should be " << interpol_correct << ". Test Failed" <<std::endl;
}
std::cout << "profiling data" <<std::endl;
std::cout.width(5);
std::cout << " Constuctor :" << constructor_elapsed << "s (" << std::fixed
<<constructor_elapsed / total_elapsed*100 << "%)" << std::endl;
std::cout << " Delaunay :" << delaunay_elapsed << "s (" << std::fixed
<<delaunay_elapsed / total_elapsed*100 << "%)" << std::endl;
std::cout << " Cleanup :" << cleanup_elapsed << "s (" << std::fixed
<<cleanup_elapsed / total_elapsed*100 << "%)" << std::endl;
std::cout << " init_inter :" << init_interpolation_elapsed << "s (" << std::fixed
<<init_interpolation_elapsed / total_elapsed*100 << "%)" << std::endl;
std::cout << " interpolate :" << interpol_elapsed << "s (" << std::fixed
<<interpol_elapsed / total_elapsed*100 << "%)" << std::endl;
std::cout << "\nadditional random point test: " << std::endl;
srand(time(NULL));
Uint nb_tests = 100;
Uint nb_failed = 0;
for (Uint i = 0; i < 100 ; ++i) {
try {
interpol_correct = 0;
for (Uint j = 0; j < dim; ++j) {
coord[j] = rand()/static_cast<Real>(RAND_MAX);
interpol_correct += coord[j];
}
if(atomic.is_inside(coord, tol)) {
--i;
continue;
}
interpol_hihi = mesh.interpolate(coord, nodal_vals);
if (fabs( interpol_hihi - interpol_correct)/interpol_correct > 2*std::numeric_limits<Real>::epsilon()) {
nb_failed += 1;
std::cout << "got " << interpol_hihi << " should have gotten " << interpol_correct << ", diff = " << fabs(interpol_correct - interpol_hihi) << " ";
print_array(coord, 3, "failed point");
std::cout << std::endl;
}
} catch (std::string &error ) {
std::cout << error << std::endl;
throw(error);
}
}
std::cout << "failed on " << nb_failed << "/" << nb_tests << "tests. ";
if (nb_failed != 0) {
std::cout << "FAILED!" << std::endl;
} else {
std::cout << "This is good" << std::endl;
}
return 0;
}
c4science · Help