Skip to content

Commit

Permalink
3d sanity test works adding mesh mtx to mtx read from file.
Browse files Browse the repository at this point in the history
  • Loading branch information
Brig Bagley committed Jun 7, 2016
1 parent 70b9226 commit 60971cf
Show file tree
Hide file tree
Showing 5 changed files with 76 additions and 211 deletions.
114 changes: 27 additions & 87 deletions src/FEMSolver.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,10 @@
#include "cutil.h"
#include "amg.h"
#include <time.h>
#include "cusp/elementwise.h"

FEMSolver::FEMSolver(
std::string fname, bool verbose) :
std::string fname, bool isTetMesh, bool verbose) :
verbose_(verbose), // output verbosity
filename_(fname), // mesh file name
maxLevels_(100), // the maximum number of levels
Expand All @@ -31,7 +32,15 @@ FEMSolver::FEMSolver(
blockSize_(256), // maximum size of a block
tetMesh_(NULL), // the tetmesh pointer
triMesh_(NULL) // the trimesh pointer
{}
{
if (isTetMesh) {
this->tetMesh_ = TetMesh::read((this->filename_ + ".node").c_str(),
(this->filename_ + ".ele").c_str(), verbose);
} else {
this->triMesh_ = TriMesh::read(this->filename_.c_str());
}
this->getMatrixFromMesh();
}

FEMSolver::~FEMSolver() {
if (this->tetMesh_ != NULL)
Expand All @@ -45,11 +54,12 @@ FEMSolver::~FEMSolver() {
* @data The set of options for the Eikonal algorithm.
* The defaults are used if nothing is provided.
*/
void FEMSolver::solveFEM(Matrix_ell_h* A_h,
void FEMSolver::solveFEM(
Vector_h_CG* x_h, Vector_h_CG* b_h) {
this->checkMatrixForValidContents(&this->A_h_);
clock_t starttime, endtime;
starttime = clock();
Matrix_ell_d_CG A_device(*A_h);
Matrix_ell_d_CG A_device(this->A_h_);
//copy back to the host
cudaThreadSynchronize();
//register configuration parameters
Expand Down Expand Up @@ -81,6 +91,10 @@ void FEMSolver::solveFEM(Matrix_ell_h* A_h,
printf("Computing time : %.10lf ms\n", duration);
}

size_t FEMSolver::getMatrixRows() {
return this->A_h_.num_rows;
}

bool FEMSolver::InitCUDA() {
int count = 0;
bool found = false;
Expand Down Expand Up @@ -123,7 +137,7 @@ void FEMSolver::checkMatrixForValidContents(Matrix_ell_h* A_h) {
}
}

void FEMSolver::getMatrixFromMesh(Matrix_ell_h* A_h) {
void FEMSolver::getMatrixFromMesh() {
if (this->triMesh_ == NULL && this->tetMesh_ == NULL)
exit(0);
Matrix_ell_d_CG A_device;
Expand Down Expand Up @@ -152,14 +166,14 @@ void FEMSolver::getMatrixFromMesh(Matrix_ell_h* A_h) {
fem3d.assemble(this->tetMesh_, A_device, RHS, true);
}
cudaThreadSynchronize();
*A_h = Matrix_ell_h(A_device);
this->A_h_ = Matrix_ell_h(A_device);
}

bool FEMSolver::compare_sparse_entry(SparseEntry_t a, SparseEntry_t b) {
return ((a.row_ != b.row_) ? (a.row_ < b.row_) : (a.col_ < b.col_));
}

int FEMSolver::readMatlabSparseMatrix(const std::string &filename, Matrix_ell_h *A_h) {
int FEMSolver::readMatlabSparseMatrix(const std::string &filename) {
//read in the description header
std::ifstream in(filename.c_str(), std::ios::binary);
if (!in.is_open()) {
Expand Down Expand Up @@ -323,11 +337,12 @@ int FEMSolver::readMatlabSparseMatrix(const std::string &filename, Matrix_ell_h
row_count = 0;
}
}
*A_h = Matrix_ell_h(A);
Matrix_ell_h original(this->A_h_);
cusp::add(A, original, this->A_h_);
return 0;
}

int FEMSolver::readMatlabNormalMatrix(const std::string &filename, std::vector<double> *A_h) {
int FEMSolver::readMatlabArray(const std::string &filename, Vector_h_CG* rhs) {
//read in the description header
std::ifstream in(filename.c_str());
if (!in.is_open()) {
Expand Down Expand Up @@ -413,15 +428,15 @@ int FEMSolver::readMatlabNormalMatrix(const std::string &filename, std::vector<d
double readInDouble;
unsigned int numValues = arrayData_length / byte_per_element;

A_h->clear();
rhs->clear();
for (int j = 0; j < numValues; ++j) {
in.read(buffer, byte_per_element);
memcpy(&readInDouble, buffer, sizeof(double));
A_h->push_back(readInDouble);
rhs->push_back(readInDouble);
}

in.close();
return numValues;
return 0;
}

int FEMSolver::writeMatlabArray(const std::string &filename, const Vector_h_CG &array) {
Expand Down Expand Up @@ -556,78 +571,3 @@ void FEMSolver::writeVTK(std::vector <float> values, std::string fname)
fclose(vtkfile);
}
}

void FEMSolver::printElementWithHeader(std::vector<double>& test, unsigned int index)
{
std::cout << "element #" << index << " = " << test[index] << std::endl;
}

void FEMSolver::printMatlabReadContents(std::vector<double>& test)
{
std::cout << "test result vector is size = " << test.size() << std::endl;
int incr = test.size() / 5;
for (int j = 0; j < test.size(); j += incr) {
printElementWithHeader(test, j);
}
if (test.size() > 0)
std::cout << "last element = " << test[test.size() - 1] << std::endl;
}

int FEMSolver::importRhsVectorFromFile(std::string filename,
Vector_h_CG& targetVector, bool verbose)
{
if (filename.empty()) {
std::string errMsg = "No matlab file provided for RHS (b) vector.";
errMsg += " Specify the file using argument at commandline using -b switch.";
std::cerr << errMsg << std::endl;
return -1;
}
std::vector<double> source;
if (FEMSolver::readMatlabNormalMatrix(filename, &source) < 0) {
std::cerr << "Failed to read matlab file for RHS (b)." << std::endl;
return -1;
}
// targetVector = static_cast<const vector<double> >(*source);
targetVector = source;
if (verbose) {
int sizeRead = targetVector.size();
std::cout << "Finished reading RHS (b) data file with ";
std::cout << sizeRead << " values." << std::endl;
}
return 0;
}

int FEMSolver::importStiffnessMatrixFromFile(std::string filename,
Matrix_ell_h* targetMatrix, bool verbose)
{
if (filename.empty()) {
std::string errMsg = "No matlab file provided for stiffness matrix (A).";
errMsg += " Specify the file using argument at commandline using -A switch.";
std::cerr << errMsg << std::endl;
return -1;
}
if (FEMSolver::readMatlabSparseMatrix(filename, targetMatrix) != 0) {
std::cerr << "Failed to read matlab file for stiffness matrix (A)." << std::endl;
return -1;
}
if (verbose) {
std::string msg = "Finished reading stiffness matrix.";
std::cout << msg << std::endl;
}
return 0;
}

void FEMSolver::debugPrintMatlabels(TetMesh* mesh)
{
std::cout << "Found " << mesh->matlabels.size() << " elements in matlabels." << std::endl;
unsigned int numZeros = 0;
for (std::vector<int>::iterator it = mesh->matlabels.begin(); it != mesh->matlabels.end(); ++it)
{
if ((*it) == 0)
numZeros++;
else
std::cout << (*it) << std::endl;
}
std::cout << numZeros << " zero values found." << std::endl;
}

33 changes: 13 additions & 20 deletions src/FEMSolver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@

/** The class that represents all of the available options for FEM */
class FEMSolver {
private:
class SparseEntry_t {
public:
int32_t row_;
Expand All @@ -21,31 +22,21 @@ class FEMSolver {
static_cast<float>(v)) {}
~SparseEntry_t() {}
};

public:
bool InitCUDA();
static bool compare_sparse_entry(SparseEntry_t a, SparseEntry_t b);
public:
FEMSolver(std::string fname = "../src/test/test_data/simple",
bool verbose = false);
bool isTetMesh = true, bool verbose = false);
virtual ~FEMSolver();
public:
void solveFEM(Matrix_ell_h* A_h, Vector_h_CG* x_h, Vector_h_CG* b_h);
void getMatrixFromMesh(Matrix_ell_h* A_h);
static int readMatlabSparseMatrix(const std::string &filename, Matrix_ell_h *A_h);
static int readMatlabNormalMatrix(const std::string &filename, std::vector<double> *A_h);
void solveFEM(Vector_h_CG* x_h, Vector_h_CG* b_h);
void getMatrixFromMesh();
int readMatlabSparseMatrix(const std::string &filename);
int readMatlabArray(const std::string &filename, Vector_h_CG* rhs);
int writeMatlabArray(const std::string &filename, const Vector_h_CG &array);
void checkMatrixForValidContents(Matrix_ell_h* A_h);
void writeVTK(std::vector <float> values, std::string fname);
void printElementWithHeader(std::vector<double>& test, unsigned int index);
void printMatlabReadContents(std::vector<double>& test);
int importRhsVectorFromFile(std::string filename,
Vector_h_CG& targetVector, bool verbose);
int importStiffnessMatrixFromFile(std::string filename,
Matrix_ell_h* targetMatrix, bool verbose);
void debugPrintMatlabels(TetMesh* mesh);
private:
bool InitCUDA();
static bool compare_sparse_entry(SparseEntry_t a, SparseEntry_t b);
public:
//data
size_t getMatrixRows();
//data members
bool verbose_; // output verbosity
std::string filename_; // mesh file name
int maxLevels_; // the maximum number of levels
Expand All @@ -71,6 +62,8 @@ class FEMSolver {
//The pointers to the meshes
TetMesh * tetMesh_;
TriMesh * triMesh_;
//The A matrix used by the solver
Matrix_ell_h A_h_;
};

#endif
65 changes: 18 additions & 47 deletions src/examples/example1.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,18 +15,15 @@

int main(int argc, char** argv)
{
//Verbose option
//Our main configuration object. We will set aspects where the
// default values are not what we desire.
FEMSolver cfg;
std::string Aname = "", bName;
cfg.filename_ = "../src/test/test_data/CubeMesh_size256step16";
//option
std::string Aname = "", bName, fname = "../src/test/test_data/CubeMesh_size256step16";
bool verbose = false;
for (int i = 0; i < argc; i++) {
if (strcmp(argv[i], "-v") == 0) {
cfg.verbose_ = true;
verbose = true;
} else if (strcmp(argv[i], "-i") == 0) {
if (i + 1 >= argc) break;
cfg.filename_ = std::string(argv[i + 1]);
fname = std::string(argv[i + 1]);
i++;
} else if (strcmp(argv[i], "-b") == 0) {
if (i + 1 >= argc) break;
Expand All @@ -38,51 +35,25 @@ int main(int argc, char** argv)
i++;
}
}
//assuming our device is zero...
cfg.device_ = 0;
// Make sure part_max_size is representative of harware limits by default
// when compiling the library, up to 64 registers were seen to be used on the
// device. We can set our max allocation based on that number
/*int max_registers_used = 64;
cfg.partitionMaxSize_ = getMaxThreads(max_registers_used, cfg.device_);
cfg.tolerance_= 1e-8;
cfg.smootherWeight_ = 0.7;
cfg.proOmega_ = 0.7;
cfg.maxIters_ = 10;
cfg.preInnerIters_ = 2;
cfg.postInnerIters_ = 3;
cfg.aggregatorType_ = 1;
cfg.solverType_ = PCG_SOLVER;
cfg.cycleType_ = V_CYCLE;
cfg.convergeType_ = ABSOLUTE_CONVERGENCE;*/
//Now we read in the mesh of choice
//TriMesh* meshPtr = TriMesh::read("mesh.ply"); //-----if we were reading a Triangle mesh
//read in the Tetmesh
cfg.tetMesh_ = TetMesh::read(
(cfg.filename_ + ".node").c_str(),
(cfg.filename_ + ".ele").c_str(), cfg.verbose_);
//The stiffness matrix A
Matrix_ell_h A_h;
if (Aname.empty()) {
//get the basic stiffness matrix (constant) by creating the mesh matrix
cfg.getMatrixFromMesh(&A_h);
} else {
//Our main configuration object. We will set aspects where the
// default values are not what we desire.
FEMSolver cfg(fname, true, verbose);
if (!Aname.empty()) {
//Import stiffness matrix (A)
if (cfg.importStiffnessMatrixFromFile(Aname, &A_h, cfg.verbose_) < 0)
return 0;
if (cfg.readMatlabSparseMatrix(Aname) != 0)
std::cerr << "Failed to read in A matrix: " << Aname << std::endl;
}
//intialize the b matrix to ones for now.
Vector_h_CG b_h(A_h.num_rows, 1.0);
Vector_h_CG b_h(cfg.getMatrixRows(), 1.0);
if (!bName.empty()) {
//Import right-hand-side single-column array (b)
if (cfg.importRhsVectorFromFile(bName, b_h, cfg.verbose_) < 0)
return 0;
if (cfg.readMatlabArray(bName, &b_h) != 0)
std::cerr << "Failed to read in b array: " << bName << std::endl;
}
//The answer vector.
Vector_h_CG x_h(A_h.num_rows, 0.0); //intial X vector
Vector_h_CG x_h(cfg.getMatrixRows(), 0.0); //intial X vector
//The final call to the solver
cfg.checkMatrixForValidContents(&A_h);
cfg.solveFEM(&A_h, &x_h, &b_h);
cfg.solveFEM(&x_h, &b_h);
//At this point, you can do what you need with the matrices.
if (cfg.writeMatlabArray("output.mat", x_h)) {
std::cerr << "failed to write matlab file." << std::endl;
Expand All @@ -95,8 +66,8 @@ int main(int argc, char** argv)
int pos = cfg.filename_.find_last_of("/");
if (pos == std::string::npos)
pos = cfg.filename_.find_last_of("\\");
cfg.filename_ = cfg.filename_.substr(pos + 1,
std::string outname = cfg.filename_.substr(pos + 1,
cfg.filename_.size() - 1);
cfg.writeVTK(vals, cfg.filename_);
cfg.writeVTK(vals, outname);
return 0;
}
Loading

0 comments on commit 60971cf

Please sign in to comment.