From 28b882ce020c50db81326f56c344c2aa0f4d1b84 Mon Sep 17 00:00:00 2001 From: Jiaming Guo Date: Mon, 13 Feb 2017 15:57:08 -0600 Subject: [PATCH] new flow header --- stim/biomodels/flow.h | 612 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ stim/biomodels/flow_dep.h | 405 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 645 insertions(+), 372 deletions(-) create mode 100644 stim/biomodels/flow_dep.h diff --git a/stim/biomodels/flow.h b/stim/biomodels/flow.h index 7154ad6..4990c02 100644 --- a/stim/biomodels/flow.h +++ b/stim/biomodels/flow.h @@ -1,405 +1,273 @@ -#pragma once -#include // Required for ofstream, etc. -#include // Required for setw -#include // Required for cout, cin, etc. -#include // Required for returning multiple values from a function - -using namespace std; - - -class flow -{ -public: - void backupToTxt(unsigned int nL, double **D, char filename[]); - tuple copySrcDesRadLen(char filename[]); - void copyToArray(int *src, int *dest, double *radii, double *len); - int getDangleNodes(int datarow, int numNodes, int *row, int *column, int *dangleNodes); - void inversion(double **a, int n, double **b); - -protected: - float determinant(double **a, int n); - int minor(double **src, double **dest, int row, int col, int order); -}; - -/* Function to find the dangle nodes in a network */ -// Created by Cherub P. Harder (8/10/2015), U of Houston -// Modified by Cherub P. Harder on 8/12/2015 -int flow::getDangleNodes(int datarow, int numNodes, int *column1, int *column2, int *dangleNodes) -{ - int count = datarow, diff1 = 0, diff2 = 0, numPress = 0, st = 0; - - // Find matching nodes within column2 - for( int i = 0; i < count; i++ ) - { - for( int y = i+1; y < datarow; y++ ) - { - if( column2[i] == column2[y] ) // Is there a match? - { - st = column2[i]; // Save the matching node -// cout << endl << column2[i] << " = " << column2[y] << endl; // Display the matching nodes - memmove(column2+i, column2+i+1, (datarow-(i+1)) * sizeof(column2[0])); // Move up the rows - // taking the places of the rows before them starting - // with where the matching node is located - column2[datarow-1] = st; // Place the matching node at the very end of the array-- - // this is for comparison purpose so that the other match- - // ing node will be moved as well and then omitted later. - diff1++; // Count the matching node - - // Display the updated array (with the matching node moved to the bottommost row) -/* cout << "Updated array:" << endl; - for( int k = 0; k < datarow; k++ ) - cout << column2[k] << endl; -*/ - // Decrement the counters - // NOTE: The counters need to be decremented because the rows had been moved up, so the same - // locations need to be read again because they contain different values now after the move. - i--; // Decrement i to read the node that took over the place - // of the matching node. Otherwise, it will be skipped. - y--; // Decrement y to read the following node for comparison - count--; // The maximum count need to be decremented so that the - // matching nodes that had been moved will not be read again. - // However, the maximum count (datarow) for finding a match - // will not be decremented because the remaining matching - // node that has not been moved yet needs to be moved and - // the only way to do that is to match it with its duplicate. - } - } - } - - // Store the nodes that have no duplicates - // NOTE: This will only save the nodes that have not been moved to the bottom. -// cout << "\ndangleNodes array:" << endl; - for( int j = 0; j < datarow-diff1; j++ ) - { - dangleNodes[numPress] = column2[j]; -// cout << dangleNodes[j] << endl; // DELETE!!! - numPress++; // Count the non-duplicated node - } - - // Find if the non-duplicated nodes have a match from column1 - count = datarow-diff1; // Reinitialize the counter - - for( int i = 0; i < count; i++ ) - { - for( int j = 0; j < datarow; j++ ) - { - if( dangleNodes[i] == column1[j] ) // Is there a match? - { - st = column1[j]; // Save the matching node -// cout << endl << dangleNodes[i] << " = " << column1[j] << endl; // Display the matching nodes - memmove(dangleNodes+i, dangleNodes+i+1, (datarow-diff1-(i+1)) * sizeof(dangleNodes[0])); - dangleNodes[count-1] = st; // Move the matching node to the bottom of the array - diff2++; // Count the matching node - - // Display the updated array -/* cout << "Updated dangleNodes array:" << endl; - for( int k = 0; k < count-1; k++ ) - { - cout << dangleNodes[k] << endl; +#ifndef STIM_FLOW_H +#define STIM_FLOW_H + +#include +#include + +//STIM include +#include +#include +#include + +#ifdef __CUDACC__ +#include +#include +#endif + +namespace stim { + template + struct triple { + A first; + B second; + C third; + }; + + template + struct bridge { + std::vector v; // vertices' indices + std::vector > V; // vertices' coordinates + T l; // length + T r; // radii + T deltaP; // pressure drop + T Q; // volume flow rate + }; + + template + class flow { + + private: + + // calculate the cofactor of elemen[row][col] + void get_minor(T** src, T** dest, int row, int col, int order) { + + // index of element to be copied + int rowCount = 0; + int colCount = 0; + + for (int i = 0; i < order; i++) { + if (i != row) { + colCount = 0; + for (int j = 0; j < order; j++) { + // when j is not the element + if (j != col) { + dest[rowCount][colCount] = src[i][j]; + colCount++; + } + } + rowCount++; } -*/ - // Decrement the counters - i--; - j--; - count--; - numPress--; // Decrement to the exact number of dangle nodes } } - } - return numPress; // Number of dangle nodes -} + // calculate the det() + T determinant(T** mat, int order) { + // degenate case when n = 1 + if (order == 1) + return mat[0][0]; -// Function to make a backup copy of the contents of a matrix to a .txt file -// Created by Cherub P. Harder (8/10/2015), U of Houston -void flow::backupToTxt(unsigned int nL, double **D, char filename[]) -{ - ofstream output_file(filename); - - for( unsigned int i = 0; i < nL; i++ ) - { - for( int j = 0; j < 4; j++ ) - { - if( j < 3 ) - output_file << D[i][j] << "\t"; - - else - output_file << D[i][j]; - } + T det = 0.0; // determinant value - output_file << "\n"; - } + // allocate the cofactor matrix + T** minor = (T**)malloc((order - 1) * sizeof(T*)); + for (int i = 0; i < order - 1; i++) + minor[i] = (T*)malloc((order - 1) * sizeof(T)); - output_file.close( ); -} + for (int i = 0; i < order; i++) { -// Function to make separate copies of the source nodes, destination nodes, radii, and lengths -// Created by Cherub P. Harder (8/10/2015), U of Houston -tuple flow::copySrcDesRadLen(char filename[]) -{ - int cnt = 0, numElements = 0, numNodes = 0; - float number = 0.0; - ofstream srcData("srcCol.txt"); // A .txt file to store the source nodes - ofstream desData("destCol.txt"); // A .txt file to store the destination nodes - ofstream radiiData("radii.txt"); // A .txt file to store the radii - ofstream lenData("lengths.txt"); // A .txt file to store the lengths - FILE *fp = fopen(filename, "r"); // Create a variable of type FILE* and open the file using - // the fopen function and assign the file to the variable - // Check if the file exists - if(fp == NULL) // Alternative: if(!fp) - { - printf("Error! File does not exist.\n"); - getchar( ); // Pause - exit(-1); // NOTE: Must include stdlib.h. - } - - // Store data to their respective .txt files - while(fscanf(fp, "%f", &number) == 1) - { - cnt++; // Increment counter - - // Store to srcCol.txt - if(cnt == 1) - srcData << number << endl; - - // Store to destCol.txt - if(cnt == 2) - desData << number << endl; - - // Save the current number of nodes - if(cnt < 3) - { - if(number > numNodes) - numNodes = (int)number; - } + // get minor of element(0, i) + get_minor(mat, minor, 0, i, order); - // Store to radii.txt - if(cnt == 3) - radiiData << number << endl; + // recursion + det += (i % 2 == 1 ? -1.0 : 1.0) * mat[0][i] * determinant(minor, order - 1); + } - // Store to lengths.txt - if(cnt == 4) - { - lenData << number << endl; + // release memory + for (int i = 0; i < order - 1; i++) + free(minor[i]); + free(minor); - numElements++; // Count the elements - cnt = 0; // Reset counter + return det; } - } - - srcData.close( ); - desData.close( ); - radiiData.close( ); - lenData.close( ); - - return make_tuple(numNodes, numElements); // Return two values -} - - -// Function to copy data for .txt files to their respective arrays -// Created by Cherub P. Harder (8/11/2015), U of Houston -void flow::copyToArray(int *src, int *dest, double *radii, double *len) -{ - int v = 0; - double tmp = 0, R = 0, L = 0; - - // Store source node values to the array src - ifstream readSrc("srcCol.txt"); - - while( readSrc >> tmp ) - { - src[v] = (int)tmp; - v++; - } - - readSrc.close( ); - // Store destination node values to the array dest - v = 0; // Reset counter - ifstream readDest("destCol.txt"); + public: + std::vector > V; // velocity + std::vector > Q; // volume flow rate + std::vector > deltaP; // pressure drop + std::vector pressure; // pressure + + flow() {} // default constructor + + // calculate the flow rate of 3D model(circle cross section) + void calculate_flow_rate(unsigned e, T r) { + stim::triple tmp_Q; + tmp_Q.first = V[e].first; // copy the vertices information + tmp_Q.second = V[e].second; + tmp_Q.third = V[e].third * stim::PI * pow(r, 2); // UNITS: uL/s + Q.push_back(tmp_Q); // push back the volume flow rate information for every edge + } - while( readDest >> tmp ) - { - dest[v] = (int)tmp; - v++; - } + // calculate the flow rate of 2D model(rectangular cross section) + void calculate_flow_rate(unsigned e, T r, T h) { + stim::triple tmp_Q; + tmp_Q.first = V[e].first; // copy the vertices information + tmp_Q.second = V[e].second; + tmp_Q.third = V[e].third * h * r; // UNITS: uL/s = mm^3/s + Q.push_back(tmp_Q); // push back the volume flow rate information for every edge + } - readDest.close( ); + // calculate the pressure drop of 3D model(circle cross section) + void calculate_deltaP(unsigned e, T u, T l, T r) { + stim::triple tmp_deltaP; + tmp_deltaP.first = V[e].first; // copy the vertices information + tmp_deltaP.second = V[e].second; + tmp_deltaP.third = (8 * u * l * Q[e].third) / (stim::PI * pow(r, 4)); // UNITS: g/mm/s^2 = Pa + deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge + } - // Store radius values to the array radii - v = 0; // Reset counter - ifstream readRad("radii.txt"); + // calculate the pressure drop of 2D model(rectangular cross section) + void calculate_deltaP(unsigned e, T u, T l, T r, T h) { + stim::triple tmp_deltaP; + tmp_deltaP.first = V[e].first; // copy the vertices information + tmp_deltaP.second = V[e].second; + tmp_deltaP.third = (12 * u * l * Q[e].third) / (stim::PI * h * pow(r, 3)); // UNITS: g/mm/s^2 = Pa + deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge + } - while( readRad >> tmp ) - { - radii[v] = tmp; - v++; - } + // better way to do this??? + // find the maximum and minimum pressure positions + void find_max_min_pressure(size_t n_e, size_t n_v, unsigned& max, unsigned& min) { + std::vector P(n_v, FLT_MAX); + // set one to reference + P[Q[0].first] = 0.0; + unsigned first = 0; + unsigned second = 0; + // calculate all the relative pressure in brute force manner + for (unsigned e = 0; e < n_e; e++) { + // assuming the obj file stores in a straight order, in other words, like swc file + first = Q[e].first; + second = Q[e].second; + if (P[first] != FLT_MAX) // if pressure at start vertex is known + P[second] = P[first] - deltaP[e].third; + else if (P[second] != FLT_MAX) // if pressure at end vertex is known + P[first] = P[second] + deltaP[e].third; + } - readRad.close( ); + // find the maximum and minimum pressure position + auto m1 = std::max_element(P.begin(), P.end()); // temporarily max number + auto m2 = std::min_element(P.begin(), P.end()); // temporarily min number - // Store length values to the array len - v = 0; // Reset counter - ifstream readLen("lengths.txt"); + max = std::distance(P.begin(), m1); + min = std::distance(P.begin(), m2); - while( readLen >> tmp ) - { - len[v] = tmp; - v++; - } + T tmp_m = *m2; + // Now set the lowest pressure port to reference pressure(0.0 Pa) + for (unsigned i = 0; i < n_v; i++) + P[i] -= tmp_m; - readLen.close( ); -} + for (unsigned i = 0; i < n_v; i++) + pressure.push_back(P[i]); + } + void inversion(T** A, int order, T** C) { + +#ifdef __CUDACC__ + + // convert from double pointer to single pointer, make it flat + T* Aflat = (T*)malloc(order * order * sizeof(T)); + for (unsigned i = 0; i < order; i++) + for (unsigned j = 0; j < order; j++) + Aflat[i * order + j] = A[i][j]; + T* Cflat = (T*)malloc(order * order * sizeof(T)); + + // create device pointer + T* d_Aflat; // flat original matrix + T* d_Cflat; // flat inverse matrix + T** d_A; // put the flat original matrix into another array of pointer + T** d_C; + int *d_P; + int *d_INFO; + + // allocate memory on device + HANDLE_ERROR(cudaMalloc((void**)&d_Aflat, order * order * sizeof(T))); + HANDLE_ERROR(cudaMalloc((void**)&d_Cflat, order * order * sizeof(T))); + HANDLE_ERROR(cudaMalloc((void**)&d_A, sizeof(T*))); + HANDLE_ERROR(cudaMalloc((void**)&d_C, sizeof(T*))); + HANDLE_ERROR(cudaMalloc((void**)&d_P, order * 1 * sizeof(int))); + HANDLE_ERROR(cudaMalloc((void**)&d_INFO, 1 * sizeof(int))); + + // copy matrix from host to device + HANDLE_ERROR(cudaMemcpy(d_Aflat, Aflat, order * order * sizeof(T), cudaMemcpyHostToDevice)); + + // copy matrix from device to device + HANDLE_ERROR(cudaMemcpy(d_A, &d_Aflat, sizeof(T*), cudaMemcpyHostToDevice)); + HANDLE_ERROR(cudaMemcpy(d_C, &d_Cflat, sizeof(T*), cudaMemcpyHostToDevice)); -// Function to find the inverse of a square matrix -void flow::inversion(double **a, int n, double **b) -{ - // Get 1 over the determinant of A - double det = (double)(1.0/determinant(a, n)); -// cerr << "\n1/det(C) = " << det << endl; // DELETE!!! - - // Memory allocation - double *tmp = new double[(n-1) * (n-1)]; - double **m = new double * [n-1]; - for( int i = 0; i < n-1; i++ ) - m[i] = tmp + ( i * (n-1) ); - - for( int j = 0; j < n; j++) - { - for( int i = 0; i < n; i++ ) - { - // Get the cofactor (matrix) of a(j,i) - minor(a, m, j, i, n); - b[i][j] = det * determinant( m, n-1 ); - if( (i+j)%2 == 1 ) - b[i][j] = -b[i][j]; - } - } - - // Release memory - // Delete [] minor[0]; - delete [] tmp; - delete [] m; -} + // calculate the inverse of matrix based on cuBLAS + cublasHandle_t handle; + CUBLAS_HANDLE_ERROR(cublasCreate_v2(&handle)); // create cuBLAS handle object + CUBLAS_HANDLE_ERROR(cublasSgetrfBatched(handle, order, d_A, order, d_P, d_INFO, 1)); -// Function to find the determinant of a matrix using recursion -// Contribution by Edward Popko -// Modified by Cherub P. Harder (7/15/2015), U of Houston -// Arguments: a(double **) - pointer to a pointer of an arbitrary square matrix -// n(int) - dimension of the square matrix -float flow::determinant(double **a, int n) -{ - int i, j, j1, j2; // General loop and matrix subscripts - double det = 0; // Initialize determinant - double **m = NULL; // Pointer to pointer to implement 2D square array - - // Display contents of matrix C (DELETE!!!) -/* std::cout << "\nThe updated matrix C:\n"; - for( int j = 0; j < n; ++j ) - { - std::cerr << "\t"; - - for( int k = 0; k < n; ++k ) - std::cerr << left << setw(15) << a[j][k]; - - std::cerr << endl; - } - - getchar(); // DELETE!!!*/ - - if(n < 1) { } // Error condition - should never get here - - else if (n == 1) // Should never get here - { - det = a[0][0]; - } - - else if(n == 2) // Basic 2x2 sub-matrix determinate definition - { // When n == 2, this ends the recursion series - det = a[0][0] * a[1][1] - a[1][0] * a[0][1]; - } - // Recursion continues, solve next sub-matrix - else // Solve the next minor by building a sub-matrix - { - det = 0; // Initialize determinant of sub-matrix - - for (j1 = 0; j1 < n; j1++) // For each column in sub-matrix get space for the - { // pointer list - m = (double **) malloc((n-1) * sizeof(double *)); - - for (i = 0; i < n-1; i++) - m[i] = (double *) malloc((n-1)* sizeof(double)); - // i[0][1][2][3] first malloc - // m -> + + + + space for 4 pointers - // | | | | j second malloc - // | | | +-> _ _ _ [0] pointers to - // | | +----> _ _ _ [1] and memory for - // | +-------> _ a _ [2] 4 doubles - // +----------> _ _ _ [3] - // - // a[1][2] - // Build sub-matrix with minor elements excluded - - for (i = 1; i < n; i++) + int INFO = 0; + HANDLE_ERROR(cudaMemcpy(&INFO, d_INFO, sizeof(int), cudaMemcpyDeviceToHost)); + if (INFO == order) { - j2 = 0 ; // Start at first sum-matrix column position - // Loop to copy source matrix less one column - for (j = 0; j < n; j++) - { - if (j == j1) continue; // Do not copy the minor column element - - m[i-1][j2] = a[i][j]; // Copy source element into new sub-matrix - // i-1 because new sub-matrix is one row - // (and column) smaller with excluded minors - j2++; // Move to next sub-matrix column position + std::cout << "Factorization Failed : Matrix is singular." << std::endl; + cudaDeviceReset(); + exit(1); + } + + CUBLAS_HANDLE_ERROR(cublasSgetriBatched(handle, order, (const T **)d_A, order, d_P, d_C, order, d_INFO, 1)); + + CUBLAS_HANDLE_ERROR(cublasDestroy_v2(handle)); + + // copy inverse matrix from device to device + HANDLE_ERROR(cudaMemcpy(&d_Cflat, d_C, sizeof(T*), cudaMemcpyDeviceToHost)); + + // copy inverse matrix from device to host + HANDLE_ERROR(cudaMemcpy(Cflat, d_Cflat, order * order * sizeof(T), cudaMemcpyDeviceToHost)); + + for(unsigned i = 0; i < order; i++) + memcpy(C[i], &Cflat[i*order], order * sizeof(T*)); + + // clear up + free(Aflat); + free(Cflat); + HANDLE_ERROR(cudaFree(d_Aflat)); + HANDLE_ERROR(cudaFree(d_Cflat)); + HANDLE_ERROR(cudaFree(d_A)); + HANDLE_ERROR(cudaFree(d_C)); + HANDLE_ERROR(cudaFree(d_P)); + HANDLE_ERROR(cudaFree(d_INFO)); + +#else + // get the determinant of a + double det = 1.0 / determinant(A, order); + + // memory allocation + T* tmp = (T*)malloc((order - 1)*(order - 1) * sizeof(T)); + T** minor = (T**)malloc((order - 1) * sizeof(T*)); + for (int i = 0; i < order - 1; i++) + minor[i] = tmp + (i * (order - 1)); + + for (int j = 0; j < order; j++) { + for (int i = 0; i < order; i++) { + // get the co-factor (matrix) of A(j,i) + get_minor(A, minor, j, i, order); + C[i][j] = det * determinant(minor, order - 1); + if ((i + j) % 2 == 1) + C[i][j] = -C[i][j]; } } - - det += (double)pow(-1.0, 1.0 + j1 + 1.0) * a[0][j1] * determinant(m, n-1); - // Sum x raised to y power - // recursively get determinant of next - // sub-matrix which is now one - // row & column smaller - - for (i = 0; i < n-1; i++) free(m[i]); // Free the storage allocated to - // this minor's set of pointers - free(m); // Free the storage for the original - // pointer to pointer + + // release memory + free(tmp); + free(minor); +#endif } - } - - return(det); + }; } - -// Function to calculate the cofactor of element (row, col) -int flow::minor(double **src, double **dest, int row, int col, int order) -{ - // Indicate which col and row is being copied to dest - int colCount=0,rowCount=0; - - for(int i = 0; i < order; i++) - { - if(i != row) - { - colCount = 0; - for(int j = 0; j < order; j++) - { - // When j is not the element - if( j != col ) - { - dest[rowCount][colCount] = src[i][j]; - colCount++; - } - } - - rowCount++; - } - } - - return 1; -} \ No newline at end of file +#endif \ No newline at end of file diff --git a/stim/biomodels/flow_dep.h b/stim/biomodels/flow_dep.h new file mode 100644 index 0000000..7154ad6 --- /dev/null +++ b/stim/biomodels/flow_dep.h @@ -0,0 +1,405 @@ +#pragma once +#include // Required for ofstream, etc. +#include // Required for setw +#include // Required for cout, cin, etc. +#include // Required for returning multiple values from a function + +using namespace std; + + +class flow +{ +public: + void backupToTxt(unsigned int nL, double **D, char filename[]); + tuple copySrcDesRadLen(char filename[]); + void copyToArray(int *src, int *dest, double *radii, double *len); + int getDangleNodes(int datarow, int numNodes, int *row, int *column, int *dangleNodes); + void inversion(double **a, int n, double **b); + +protected: + float determinant(double **a, int n); + int minor(double **src, double **dest, int row, int col, int order); +}; + +/* Function to find the dangle nodes in a network */ +// Created by Cherub P. Harder (8/10/2015), U of Houston +// Modified by Cherub P. Harder on 8/12/2015 +int flow::getDangleNodes(int datarow, int numNodes, int *column1, int *column2, int *dangleNodes) +{ + int count = datarow, diff1 = 0, diff2 = 0, numPress = 0, st = 0; + + // Find matching nodes within column2 + for( int i = 0; i < count; i++ ) + { + for( int y = i+1; y < datarow; y++ ) + { + if( column2[i] == column2[y] ) // Is there a match? + { + st = column2[i]; // Save the matching node +// cout << endl << column2[i] << " = " << column2[y] << endl; // Display the matching nodes + memmove(column2+i, column2+i+1, (datarow-(i+1)) * sizeof(column2[0])); // Move up the rows + // taking the places of the rows before them starting + // with where the matching node is located + column2[datarow-1] = st; // Place the matching node at the very end of the array-- + // this is for comparison purpose so that the other match- + // ing node will be moved as well and then omitted later. + diff1++; // Count the matching node + + // Display the updated array (with the matching node moved to the bottommost row) +/* cout << "Updated array:" << endl; + for( int k = 0; k < datarow; k++ ) + cout << column2[k] << endl; +*/ + // Decrement the counters + // NOTE: The counters need to be decremented because the rows had been moved up, so the same + // locations need to be read again because they contain different values now after the move. + i--; // Decrement i to read the node that took over the place + // of the matching node. Otherwise, it will be skipped. + y--; // Decrement y to read the following node for comparison + count--; // The maximum count need to be decremented so that the + // matching nodes that had been moved will not be read again. + // However, the maximum count (datarow) for finding a match + // will not be decremented because the remaining matching + // node that has not been moved yet needs to be moved and + // the only way to do that is to match it with its duplicate. + } + } + } + + // Store the nodes that have no duplicates + // NOTE: This will only save the nodes that have not been moved to the bottom. +// cout << "\ndangleNodes array:" << endl; + for( int j = 0; j < datarow-diff1; j++ ) + { + dangleNodes[numPress] = column2[j]; +// cout << dangleNodes[j] << endl; // DELETE!!! + numPress++; // Count the non-duplicated node + } + + // Find if the non-duplicated nodes have a match from column1 + count = datarow-diff1; // Reinitialize the counter + + for( int i = 0; i < count; i++ ) + { + for( int j = 0; j < datarow; j++ ) + { + if( dangleNodes[i] == column1[j] ) // Is there a match? + { + st = column1[j]; // Save the matching node +// cout << endl << dangleNodes[i] << " = " << column1[j] << endl; // Display the matching nodes + memmove(dangleNodes+i, dangleNodes+i+1, (datarow-diff1-(i+1)) * sizeof(dangleNodes[0])); + dangleNodes[count-1] = st; // Move the matching node to the bottom of the array + diff2++; // Count the matching node + + // Display the updated array +/* cout << "Updated dangleNodes array:" << endl; + for( int k = 0; k < count-1; k++ ) + { + cout << dangleNodes[k] << endl; + } +*/ + // Decrement the counters + i--; + j--; + count--; + numPress--; // Decrement to the exact number of dangle nodes + } + } + } + + return numPress; // Number of dangle nodes +} + + +// Function to make a backup copy of the contents of a matrix to a .txt file +// Created by Cherub P. Harder (8/10/2015), U of Houston +void flow::backupToTxt(unsigned int nL, double **D, char filename[]) +{ + ofstream output_file(filename); + + for( unsigned int i = 0; i < nL; i++ ) + { + for( int j = 0; j < 4; j++ ) + { + if( j < 3 ) + output_file << D[i][j] << "\t"; + + else + output_file << D[i][j]; + } + + output_file << "\n"; + } + + output_file.close( ); +} + + +// Function to make separate copies of the source nodes, destination nodes, radii, and lengths +// Created by Cherub P. Harder (8/10/2015), U of Houston +tuple flow::copySrcDesRadLen(char filename[]) +{ + int cnt = 0, numElements = 0, numNodes = 0; + float number = 0.0; + ofstream srcData("srcCol.txt"); // A .txt file to store the source nodes + ofstream desData("destCol.txt"); // A .txt file to store the destination nodes + ofstream radiiData("radii.txt"); // A .txt file to store the radii + ofstream lenData("lengths.txt"); // A .txt file to store the lengths + FILE *fp = fopen(filename, "r"); // Create a variable of type FILE* and open the file using + // the fopen function and assign the file to the variable + // Check if the file exists + if(fp == NULL) // Alternative: if(!fp) + { + printf("Error! File does not exist.\n"); + getchar( ); // Pause + exit(-1); // NOTE: Must include stdlib.h. + } + + // Store data to their respective .txt files + while(fscanf(fp, "%f", &number) == 1) + { + cnt++; // Increment counter + + // Store to srcCol.txt + if(cnt == 1) + srcData << number << endl; + + // Store to destCol.txt + if(cnt == 2) + desData << number << endl; + + // Save the current number of nodes + if(cnt < 3) + { + if(number > numNodes) + numNodes = (int)number; + } + + // Store to radii.txt + if(cnt == 3) + radiiData << number << endl; + + // Store to lengths.txt + if(cnt == 4) + { + lenData << number << endl; + + numElements++; // Count the elements + cnt = 0; // Reset counter + } + } + + srcData.close( ); + desData.close( ); + radiiData.close( ); + lenData.close( ); + + return make_tuple(numNodes, numElements); // Return two values +} + + +// Function to copy data for .txt files to their respective arrays +// Created by Cherub P. Harder (8/11/2015), U of Houston +void flow::copyToArray(int *src, int *dest, double *radii, double *len) +{ + int v = 0; + double tmp = 0, R = 0, L = 0; + + // Store source node values to the array src + ifstream readSrc("srcCol.txt"); + + while( readSrc >> tmp ) + { + src[v] = (int)tmp; + v++; + } + + readSrc.close( ); + + // Store destination node values to the array dest + v = 0; // Reset counter + ifstream readDest("destCol.txt"); + + while( readDest >> tmp ) + { + dest[v] = (int)tmp; + v++; + } + + readDest.close( ); + + // Store radius values to the array radii + v = 0; // Reset counter + ifstream readRad("radii.txt"); + + while( readRad >> tmp ) + { + radii[v] = tmp; + v++; + } + + readRad.close( ); + + // Store length values to the array len + v = 0; // Reset counter + ifstream readLen("lengths.txt"); + + while( readLen >> tmp ) + { + len[v] = tmp; + v++; + } + + readLen.close( ); +} + + +// Function to find the inverse of a square matrix +void flow::inversion(double **a, int n, double **b) +{ + // Get 1 over the determinant of A + double det = (double)(1.0/determinant(a, n)); +// cerr << "\n1/det(C) = " << det << endl; // DELETE!!! + + // Memory allocation + double *tmp = new double[(n-1) * (n-1)]; + double **m = new double * [n-1]; + for( int i = 0; i < n-1; i++ ) + m[i] = tmp + ( i * (n-1) ); + + for( int j = 0; j < n; j++) + { + for( int i = 0; i < n; i++ ) + { + // Get the cofactor (matrix) of a(j,i) + minor(a, m, j, i, n); + b[i][j] = det * determinant( m, n-1 ); + if( (i+j)%2 == 1 ) + b[i][j] = -b[i][j]; + } + } + + // Release memory + // Delete [] minor[0]; + delete [] tmp; + delete [] m; +} + + +// Function to find the determinant of a matrix using recursion +// Contribution by Edward Popko +// Modified by Cherub P. Harder (7/15/2015), U of Houston +// Arguments: a(double **) - pointer to a pointer of an arbitrary square matrix +// n(int) - dimension of the square matrix +float flow::determinant(double **a, int n) +{ + int i, j, j1, j2; // General loop and matrix subscripts + double det = 0; // Initialize determinant + double **m = NULL; // Pointer to pointer to implement 2D square array + + // Display contents of matrix C (DELETE!!!) +/* std::cout << "\nThe updated matrix C:\n"; + for( int j = 0; j < n; ++j ) + { + std::cerr << "\t"; + + for( int k = 0; k < n; ++k ) + std::cerr << left << setw(15) << a[j][k]; + + std::cerr << endl; + } + + getchar(); // DELETE!!!*/ + + if(n < 1) { } // Error condition - should never get here + + else if (n == 1) // Should never get here + { + det = a[0][0]; + } + + else if(n == 2) // Basic 2x2 sub-matrix determinate definition + { // When n == 2, this ends the recursion series + det = a[0][0] * a[1][1] - a[1][0] * a[0][1]; + } + // Recursion continues, solve next sub-matrix + else // Solve the next minor by building a sub-matrix + { + det = 0; // Initialize determinant of sub-matrix + + for (j1 = 0; j1 < n; j1++) // For each column in sub-matrix get space for the + { // pointer list + m = (double **) malloc((n-1) * sizeof(double *)); + + for (i = 0; i < n-1; i++) + m[i] = (double *) malloc((n-1)* sizeof(double)); + // i[0][1][2][3] first malloc + // m -> + + + + space for 4 pointers + // | | | | j second malloc + // | | | +-> _ _ _ [0] pointers to + // | | +----> _ _ _ [1] and memory for + // | +-------> _ a _ [2] 4 doubles + // +----------> _ _ _ [3] + // + // a[1][2] + // Build sub-matrix with minor elements excluded + + for (i = 1; i < n; i++) + { + j2 = 0 ; // Start at first sum-matrix column position + // Loop to copy source matrix less one column + for (j = 0; j < n; j++) + { + if (j == j1) continue; // Do not copy the minor column element + + m[i-1][j2] = a[i][j]; // Copy source element into new sub-matrix + // i-1 because new sub-matrix is one row + // (and column) smaller with excluded minors + j2++; // Move to next sub-matrix column position + } + } + + det += (double)pow(-1.0, 1.0 + j1 + 1.0) * a[0][j1] * determinant(m, n-1); + // Sum x raised to y power + // recursively get determinant of next + // sub-matrix which is now one + // row & column smaller + + for (i = 0; i < n-1; i++) free(m[i]); // Free the storage allocated to + // this minor's set of pointers + free(m); // Free the storage for the original + // pointer to pointer + } + } + + return(det); +} + + +// Function to calculate the cofactor of element (row, col) +int flow::minor(double **src, double **dest, int row, int col, int order) +{ + // Indicate which col and row is being copied to dest + int colCount=0,rowCount=0; + + for(int i = 0; i < order; i++) + { + if(i != row) + { + colCount = 0; + for(int j = 0; j < order; j++) + { + // When j is not the element + if( j != col ) + { + dest[rowCount][colCount] = src[i][j]; + colCount++; + } + } + + rowCount++; + } + } + + return 1; +} \ No newline at end of file -- libgit2 0.21.4