hello everyone.
I am trying to assemble the finite element stiffness matrix in CSR format using MPI. I drived the node connectivity and three csr 1d arrays(row_ptr,col_ind and noz_zero value) for each chunk of elements in the function attached. But I'm not getting the right answer. Beside I get the error " process exited without calling finalize". I wonder whether I'm using mpi in a right way or not.
What I have tried:
void csr_assembly(int n,double *x,double *y,int **connectivity, double **coordinate,double *NN_vallocal)
{
int num_procs, myrank;
MPI_Status status;
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
int node_0, node_1, node_2;
int counter[n_node];
int node_connec[n_node][n_node] = { 0 };
for (int i = 0; i < n_node; i++)
{
counter[i] = 0;
}
for (int i = 0; i < n; i++)
{
node_0 = connectivity[i][1];
node_1 = connectivity[i][2];
node_2 = connectivity[i][3];
counter[node_0] = counter[node_0] + 1;
node_connec[node_0][counter[node_0]] = node_0 + 1;
counter[node_0] = counter[node_0] + 1;
node_connec[node_0][counter[node_0]] = node_1 + 1;
counter[node_0] = counter[node_0] + 1;
node_connec[node_0][counter[node_0]] = node_2 + 1;
counter[node_1] = counter[node_1] + 1;
node_connec[node_1][counter[node_1]] = node_1 + 1;
counter[node_1] = counter[node_1] + 1;
node_connec[node_1][counter[node_1]] = node_2 + 1;
counter[node_1] = counter[node_1] + 1;
node_connec[node_1][counter[node_1]] = node_0 + 1;
counter[node_2] = counter[node_2] + 1;
node_connec[node_2][counter[node_2]] = node_2 + 1;
counter[node_2] = counter[node_2] + 1;
node_connec[node_2][counter[node_2]] = node_1 + 1;
counter[node_2] = counter[node_2] + 1;
node_connec[node_2][counter[node_2]] = node_0 + 1;
}
int r = sizeof(node_connec[0]) / sizeof(node_connec[0][0]);
int cc = sizeof(node_connec) / sizeof((node_connec[0]));
sortRowWise(node_connec, r, cc);
int nnz[n_node];
for (int i = 0; i < n_node; i++)
{
int counter = 0;
for (int j = 0; j < n_node; j++)
{
if (node_connec[i][j] != 0)
{
counter = counter + 1;
}
nnz[i] = counter;
}
}
vector<vector<int>> vec(n_node);
for (int i = 0; i < n_node; i++)
{
int col = { nnz[i] };
vec[i] = vector<int>(col);
for (int j = 0; j < col; j++)
{
for (int z = 0; z < n_node; z++)
{
if (node_connec[i][z] != 0)
{
vec[i][j] = node_connec[i][z];
j = j + 1;
}
}
}
}
for (int i = 0; i < n_node; i++)
{
int prev = 0;
for (int j = 0; j < vec[i].size(); j++)
{
vec[i].erase(std::unique(vec[i].begin(), vec[i].end()), vec[i].end());
}
}
for (int i = 0; i < n_node; i++)
{
for (int j = 0; j < vec[i].size(); j++)
{
vec[i][j] = vec[i][j] - 1;
}
}
fstream network;
network.open("node_connec.txt", fstream::out);
network << std::endl;
for (int i = 0; i < n_node; i++)
{
for (int j = 0; j < vec[i].size(); j++)
{
network << vec[i][j] << "\t";
}
network << std::endl;
}
network.close();
int row_ptr[n_node + 1];
row_ptr[0] = 0;
for (int i = 1; i < n_node + 1; i++)
{
row_ptr[i] = row_ptr[i - 1] + vec[i - 1].size();
}
for (int i = 0; i < n_node + 1; i++)
{
cout << "i= " << i << " " << row_ptr[i] << endl;
}
for (int i = 0; i < n_node; i++)
{
n_nz += vec[i].size();
}
cout << "The number of non-zero element is= " << n_nz << endl;
int *col_ind = VecToArr(vec);
for (int i = 0; i < n_node; i++)
{
std::cout << col_ind[i] << "\n";
}
for (int i = 0; i < n; i++)
{
int node_0 = connectivity[i][1];
int node_1 = connectivity[i][2];
int node_2 = connectivity[i][3];
double a, b, c, d;
elem_area[i] = abs(0.5*(((x[node_0] - x[node_2])*(y[node_1] - y[node_2])) - ((y[node_0] - y[node_2])*(x[node_1] - x[node_2]))));
a = x[node_0] - x[node_1];
b = y[node_0] - y[node_1];
c = x[node_0] - x[node_2];
d = y[node_0] - y[node_2];
gx1 = (d - b) / (a*d - c * b);
gy1 = (c - d) / (b*c - a * d);
a = x[node_1] - x[node_2];
b = y[node_1] - y[node_2];
c = x[node_1] - x[node_0];
d = y[node_1] - y[node_0];
gx2 = (d - b) / (a*d - c * b);
gy2 = (c - d) / (b*c - a * d);
a = x[node_2] - x[node_0];
b = y[node_2] - y[node_0];
c = x[node_2] - x[node_1];
d = y[node_2] - y[node_1];
gx3 = (d - b) / (a*d - c * b);
gy3 = (c - d) / (b*c - a * d);
Nx[i][0] = gx1;
Ny[i][0] = gy1;
Nx[i][1] = gx2;
Ny[i][1] = gy2;
Nx[i][2] = gx3;
Ny[i][2] = gy3;
NN_e[0][0] = (1.0 / 6) * elem_area[i];
NN_e[0][1] = (1.0 / 12) * elem_area[i];
NN_e[0][2] = (1.0 / 12) * elem_area[i];
NN_e[1][0] = (1.0 / 12) * elem_area[i];
NN_e[1][1] = (1.0 / 6) * elem_area[i];
NN_e[1][2] = (1.0 / 12) * elem_area[i];
NN_e[2][0] = (1.0 / 12) * elem_area[i];
NN_e[2][1] = (1.0 / 12) * elem_area[i];
NN_e[2][2] = (1.0 / 6) * elem_area[i];
for (int n_sh = 0; n_sh < 3; n_sh++)
{
int row = connectivity[i][n_sh + 1];
int col_begin = row_ptr[row];
int col_end = row_ptr[row + 1] - 1;
for (int n_v = 0; n_v < 3; n_v++)
{
int col = connectivity[i][n_v + 1];
for (int j = col_begin; j <= col_end; j++)
{
if (col_ind[j] == col)
{
NN_vallocal[j] += NN_e[n_sh][n_v];
}
}
}
}
}
}
int* VecToArr(std::vector<std::vector<int>> vec)
{
std::size_t totalsize = 0;
for (int i = 0; i < vec.size(); i++)
{
totalsize += vec[i].size();
}
int* newarr = new int[totalsize];
int* walkarr = newarr;
for (int i = 0; i < vec.size(); i++)
{
std::copy(vec[i].begin(), vec[i].end(), walkarr);
walkarr += vec[i].size();
}
return newarr;
}
void sortRowWise(int m[][n_node], int r, int c)
{
for (int i = 0; i < r; i++)
{
for (int j = 0; j < c; j++)
{
for (int k = 0; k < c - j - 1; k++)
{
if (m[i][k] > m[i][k + 1])
{
swap(m[i][k], m[i][k + 1]);
}
}
}
}
fstream rowwisesorted_nodeconnec;
rowwisesorted_nodeconnec.open("node_connecm1.txt", fstream::out);
rowwisesorted_nodeconnec << std::endl;
for (int i = 0; i < n_node; i++)
{
for (int j = 0; j < n_node; j++)
{
rowwisesorted_nodeconnec << m[i][j] << "\t";
}
rowwisesorted_nodeconnec << std::endl;
}
rowwisesorted_nodeconnec.close();
}
Submit an article or tip