278 lines
11 KiB
C++
278 lines
11 KiB
C++
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <iostream>
|
|
#include <exception>
|
|
#include <stdexcept>
|
|
#include <fstream>
|
|
|
|
#include "shared_ptr.h"
|
|
#include "common/UnitTest.h"
|
|
#include "common/Utilities.h"
|
|
#include "common/MPI_Helpers.h"
|
|
#include "IO/MeshDatabase.h"
|
|
#include "IO/Reader.h"
|
|
#include "IO/Writer.h"
|
|
|
|
|
|
inline bool approx_equal( const Point& A, const Point& B )
|
|
{
|
|
double tol = 1e-8*sqrt(A.x*A.x+A.y*A.y+A.z*A.z);
|
|
return fabs(A.x-B.x)<=tol && fabs(A.y-B.y)<=tol && fabs(A.z-B.z)<=tol;
|
|
}
|
|
inline bool approx_equal( const double& A, const double& B )
|
|
{
|
|
return fabs(A-B)<=1e-8*fabs(A+B);
|
|
}
|
|
|
|
|
|
inline double distance( const Point& p )
|
|
{
|
|
return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
|
|
}
|
|
|
|
|
|
|
|
int main(int argc, char **argv)
|
|
{
|
|
int rank,nprocs;
|
|
MPI_Init(&argc,&argv);
|
|
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
|
|
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
|
|
Utilities::setAbortBehavior(true,true,false);
|
|
Utilities::setErrorHandlers();
|
|
UnitTest ut;
|
|
|
|
// Create some points
|
|
const int N_points = 8;
|
|
const int N_tri = 12;
|
|
double x[8] = { 0, 1, 0, 1, 0, 1, 0, 1 };
|
|
double y[8] = { 0, 0, 1, 1, 0, 0, 1, 1 };
|
|
double z[8] = { 0, 0, 0, 0, 1, 1, 1, 1 };
|
|
int tri[N_tri][3] = {
|
|
{0,1,3}, {0,3,2}, {4,5,7}, {4,7,6}, // z faces
|
|
{0,1,4}, {1,4,5}, {2,3,6}, {3,6,7}, // y faces
|
|
{0,2,4}, {2,4,6}, {1,3,5}, {3,5,7} // x faces
|
|
};
|
|
|
|
// Create the meshes
|
|
std::shared_ptr<IO::PointList> set1( new IO::PointList(N_points) );
|
|
for (int i=0; i<N_points; i++) {
|
|
set1->points[i].x = x[i];
|
|
set1->points[i].y = y[i];
|
|
set1->points[i].z = z[i];
|
|
}
|
|
std::shared_ptr<IO::TriMesh> trimesh( new IO::TriMesh(N_tri,set1) );
|
|
for (int i=0; i<N_tri; i++) {
|
|
trimesh->A[i] = tri[i][0];
|
|
trimesh->B[i] = tri[i][1];
|
|
trimesh->C[i] = tri[i][2];
|
|
}
|
|
std::shared_ptr<IO::TriList> trilist( new IO::TriList(*trimesh) );
|
|
for (int i=0; i<N_tri; i++) {
|
|
Point A(x[tri[i][0]],y[tri[i][0]],z[tri[i][0]]);
|
|
Point B(x[tri[i][1]],y[tri[i][1]],z[tri[i][1]]);
|
|
Point C(x[tri[i][2]],y[tri[i][2]],z[tri[i][2]]);
|
|
if ( !approx_equal(trilist->A[i],A) || !approx_equal(trilist->B[i],B) || !approx_equal(trilist->C[i],C) )
|
|
{
|
|
printf("Failed to create trilist\n");
|
|
return -1;
|
|
}
|
|
}
|
|
RankInfoStruct rank_data( rank, nprocs, 1, 1 );
|
|
std::shared_ptr<IO::DomainMesh> domain( new IO::DomainMesh(rank_data,6,7,8,1.0,1.0,1.0) );
|
|
|
|
// Create the variables
|
|
std::shared_ptr<IO::Variable> dist_set1( new IO::Variable() );
|
|
std::shared_ptr<IO::Variable> dist_list( new IO::Variable() );
|
|
std::shared_ptr<IO::Variable> domain_var( new IO::Variable() );
|
|
dist_set1->dim = 1;
|
|
dist_list->dim = 1;
|
|
domain_var->dim = 1;
|
|
dist_set1->name = "Distance";
|
|
dist_list->name = "Distance";
|
|
domain_var->name = "Distance";
|
|
dist_set1->type = IO::NodeVariable;
|
|
dist_list->type = IO::NodeVariable;
|
|
domain_var->type = IO::VolumeVariable;
|
|
dist_set1->data.resize( N_points );
|
|
for (int i=0; i<N_points; i++)
|
|
dist_set1->data(i) = distance(set1->points[i]);
|
|
dist_list->data.resize( 3, N_tri );
|
|
for (int i=0; i<N_tri; i++) {
|
|
dist_list->data(0,i) = distance(trilist->A[i]);
|
|
dist_list->data(1,i) = distance(trilist->B[i]);
|
|
dist_list->data(2,i) = distance(trilist->C[i]);
|
|
}
|
|
domain_var->data.resize(domain->nx,domain->ny,domain->nz);
|
|
for (int i=0; i<domain->nx; i++) {
|
|
for (int j=0; j<domain->ny; j++) {
|
|
for (int k=0; k<domain->nz; k++) {
|
|
domain_var->data(i,j,k) = distance(Point(i,j,k));
|
|
}
|
|
}
|
|
}
|
|
|
|
// Create the MeshDataStruct
|
|
std::vector<IO::MeshDataStruct> meshData(4);
|
|
meshData[0].meshName = "pointmesh";
|
|
meshData[0].mesh = set1;
|
|
meshData[0].vars.push_back(dist_set1);
|
|
meshData[1].meshName = "trimesh";
|
|
meshData[1].mesh = trimesh;
|
|
meshData[1].vars.push_back(dist_set1);
|
|
meshData[2].meshName = "trilist";
|
|
meshData[2].mesh = trilist;
|
|
meshData[2].vars.push_back(dist_list);
|
|
meshData[3].meshName = "domain";
|
|
meshData[3].mesh = domain;
|
|
meshData[3].vars.push_back(domain_var);
|
|
|
|
// Write the data
|
|
std::vector<int> format(2,0);
|
|
format[0] = 2;
|
|
format[1] = 1;
|
|
IO::writeData( 0, meshData, format[0] );
|
|
IO::writeData( 3, meshData, format[1] );
|
|
MPI_Barrier(MPI_COMM_WORLD);
|
|
|
|
// Get a list of the timesteps
|
|
std::vector<std::string> timesteps = IO::readTimesteps("summary.LBM");
|
|
if ( timesteps.size()==2 )
|
|
ut.passes("Corrent number of timesteps");
|
|
else
|
|
ut.failure("Incorrent number of timesteps");
|
|
|
|
// Check the mesh lists
|
|
for (size_t i=0; i<timesteps.size(); i++) {
|
|
// Load the list of meshes and check its size
|
|
std::vector<IO::MeshDatabase> list = IO::getMeshList(".",timesteps[i]);
|
|
if ( list.size()==meshData.size() )
|
|
ut.passes("Corrent number of meshes found");
|
|
else
|
|
ut.failure("Incorrent number of meshes found");
|
|
// Check the number of domains for each mesh
|
|
bool pass = true;
|
|
for (size_t j=0; j<list.size(); j++)
|
|
pass = pass && (int)list[j].domains.size()==nprocs;
|
|
if ( pass ) {
|
|
ut.passes("Corrent number of domains for mesh");
|
|
} else {
|
|
ut.failure("Incorrent number of domains for mesh");
|
|
continue;
|
|
}
|
|
// For each domain, load the mesh and check its data
|
|
for (size_t j=0; j<list.size(); j++) {
|
|
for (size_t k=0; k<list[j].domains.size(); k++) {
|
|
std::shared_ptr<IO::Mesh> mesh = IO::getMesh(".",timesteps[i],list[j],k);
|
|
if ( mesh.get()==NULL ) {
|
|
printf("Failed to load %s\n",list[j].name.c_str());
|
|
pass = false;
|
|
break;
|
|
}
|
|
if ( list[j].name=="pointmesh" ) {
|
|
// Check the pointmesh
|
|
std::shared_ptr<IO::PointList> pmesh = IO::getPointList(mesh);
|
|
if ( pmesh.get()==NULL ) {
|
|
pass = false;
|
|
break;
|
|
}
|
|
if ( pmesh->points.size()!=N_points ) {
|
|
pass = false;
|
|
break;
|
|
}
|
|
}
|
|
if ( list[j].name=="trimesh" || list[j].name=="trilist" ) {
|
|
// Check the trimesh/trilist
|
|
std::shared_ptr<IO::TriMesh> mesh1 = IO::getTriMesh(mesh);
|
|
std::shared_ptr<IO::TriList> mesh2 = IO::getTriList(mesh);
|
|
if ( mesh1.get()==NULL || mesh2.get()==NULL ) {
|
|
pass = false;
|
|
break;
|
|
}
|
|
if ( mesh1->A.size()!=N_tri || mesh1->B.size()!=N_tri || mesh1->C.size()!=N_tri ||
|
|
mesh2->A.size()!=N_tri || mesh2->B.size()!=N_tri || mesh2->C.size()!=N_tri )
|
|
{
|
|
pass = false;
|
|
break;
|
|
}
|
|
const std::vector<Point>& P1 = mesh1->vertices->points;
|
|
const std::vector<int>& A1 = mesh1->A;
|
|
const std::vector<int>& B1 = mesh1->B;
|
|
const std::vector<int>& C1 = mesh1->C;
|
|
const std::vector<Point>& A2 = mesh2->A;
|
|
const std::vector<Point>& B2 = mesh2->B;
|
|
const std::vector<Point>& C2 = mesh2->C;
|
|
const std::vector<Point>& A = trilist->A;
|
|
const std::vector<Point>& B = trilist->B;
|
|
const std::vector<Point>& C = trilist->C;
|
|
for (size_t i=0; i<N_tri; i++) {
|
|
if ( !approx_equal(P1[A1[i]],A[i]) || !approx_equal(P1[B1[i]],B[i]) || !approx_equal(P1[C1[i]],C[i]) )
|
|
pass = false;
|
|
if ( !approx_equal(A2[i],A[i]) || !approx_equal(B2[i],B[i]) || !approx_equal(C2[i],C[i]) )
|
|
pass = false;
|
|
}
|
|
}
|
|
if ( list[j].name=="domain" && format[i]>2 ) {
|
|
// Check the domain mesh
|
|
const IO::DomainMesh& mesh1 = *std::dynamic_pointer_cast<IO::DomainMesh>(mesh);
|
|
if ( mesh1.nprocx!=domain->nprocx || mesh1.nprocy!=domain->nprocy || mesh1.nprocz!=domain->nprocz )
|
|
pass = false;
|
|
if ( mesh1.nx!=domain->nx || mesh1.ny!=domain->ny || mesh1.nz!=domain->nz )
|
|
pass = false;
|
|
if ( mesh1.Lx!=domain->Lx || mesh1.Ly!=domain->Ly || mesh1.Lz!=domain->Lz )
|
|
pass = false;
|
|
}
|
|
}
|
|
}
|
|
if ( pass ) {
|
|
ut.passes("Meshes loaded correctly");
|
|
} else {
|
|
ut.failure("Meshes did not load correctly");
|
|
continue;
|
|
}
|
|
// For each domain, load the variables and check their data
|
|
if ( format[i]==1 )
|
|
continue; // Format 1 does not support variables
|
|
for (size_t j=0; j<list.size(); j++) {
|
|
const IO::MeshDataStruct* mesh0 = NULL;
|
|
for (size_t k=0; k<meshData.size(); k++) {
|
|
if ( meshData[k].meshName == list[j].name ) {
|
|
mesh0 = &meshData[k];
|
|
break;
|
|
}
|
|
}
|
|
for (size_t v=0; v<list[j].variables.size(); v++) {
|
|
for (size_t k=0; k<list[j].domains.size(); k++) {
|
|
std::shared_ptr<const IO::Variable> variable =
|
|
IO::getVariable(".",timesteps[i],list[j],k,list[j].variables[v].name);
|
|
const IO::Variable& var1 = *mesh0->vars[v];
|
|
const IO::Variable& var2 = *variable;
|
|
pass = var1.name == var2.name;
|
|
pass = pass && var1.dim == var2.dim;
|
|
pass = pass && var1.type == var2.type;
|
|
pass = pass && var1.data.length() == var2.data.length();
|
|
if ( pass ) {
|
|
for (size_t m=0; m<var1.data.length(); m++)
|
|
pass = pass && approx_equal(var1.data(m),var2.data(m));
|
|
}
|
|
if ( pass ) {
|
|
ut.passes("Variables match");
|
|
} else {
|
|
ut.failure("Variables did not match");
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Finished
|
|
ut.report();
|
|
int N_errors = ut.NumFailGlobal();
|
|
MPI_Barrier(MPI_COMM_WORLD);
|
|
MPI_Finalize();
|
|
return N_errors;
|
|
}
|
|
|
|
|