Click here to Skip to main content
15,902,112 members
Please Sign up or sign in to vote.
0.00/5 (No votes)
See more:
I am trying to do finite element matrix assembly using MPI in c++.First of all instead of dividing mesh into pieces, I copied entire mesh into all processors in order to avoid using hole data. I know this is not an optimized way for doing MPI but at the time being it works for me, because in the next step I am going to use CSR format for matrix assembly. Also I broadcasted element connectivity matrix to all processors. Let's say we have 2 or more processors, n node and N elements. In this case the global stiffness matrix would be n*n. I define the local stifness matrix n*n as well and each processor will fill its own part during runtime.After assembly I get strange number in global stiffness matrix which is -1.25549e+67 and I don't know why. I would be grateful if somebody could help me.

What I have tried:

<pre>
#include <assert.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <time.h>
#include <cstdlib> 
#include <stdio.h>
#include <mpi.h> 
#include <vector> 
#include<iterator>
#include<sstream>
#include<string>
#include<cmath>
#include<algorithm>
#include<iomanip>

#define nconec    3
#define dimension 2
const int n_node = 246;
const int n_elem = 434;
double NN_e  [nconec][nconec];
double NN  [n_node][n_node];
double Nx  [n_elem][nconec];
double Ny  [n_elem][nconec];
double x   [n_node];
double y   [n_node];
double elem_area[n_elem];
double RHS[n_node][1];


using namespace std;


int main(int argc, char *argv[])
{

	int num_procs, myrank;
	MPI_Init(&argc, &argv);
	MPI_Status status;
	MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
	int local_num = n_elem / num_procs;
	double gx1;
	double gy1;
	double gx2;
	double gy2;
	double gx3;
	double gy3;


	double a, b, c, d;
	int col = 4;

	double **coordinate=new double *[n_node];
	for (int i = 0; i < n_node; i++)
		coordinate[i] = new double[col];

	int **connectivity=new int *[n_elem];
	for (int j = 0; j < n_elem; j++)
		connectivity[j] = new int[col];

	double **global_NN = new double*[n_node];
	for (int i = 0; i < n_node; i++)
		global_NN[i] = new double[n_node];

	double **local_NN = new double*[n_node];
	for (int i = 0; i < n_node; i++)
		local_NN[i] = new double[n_node];


	if (myrank == 0)
	{

	ifstream file{ "connec.txt" };
	if (!file.is_open()) return -1;

	for (int i{}; i != n_elem; ++i)
	{
		for (int j{}; j != 4; ++j)
		{
			file >> connectivity[i][j];
		}
	}


	for (int i = 0; i < n_elem; i++)
	{
		for (int j = 0; j < 4; j++)
		{
			connectivity[i][j] = connectivity[i][j] - 1;
		}
	}
	file.close();

	file.open("mesh.txt");

	if (file.is_open())
	{
		while (!file.eof())
		{
			for (int i = 0; i < n_node; i++)
			{
				for (int j = 0; j < 4; j++)
				{
					file >> coordinate[i][j];
				}
			}
		}
	}

	for (int i = 0; i < n_node; i++)
	{
		for (int j = 0; j < 4; j++)
		{
			if (j == 0)
			{
				coordinate[i][j] = coordinate[i][j] - 1;
			}
		}
	}


	for (int i = 0; i < n_elem; i++)
	{
		for (int j = 0; j < 4; j++)
		{
			cout << connectivity[i][j] << "\t";
		}
		cout << endl;
	}


		for (int i = 0; i < n_node; i++)
		{
			for (int j = 0; j < 4; j++)
			{
				cout << coordinate[i][j] << "\t";
			}
			cout << endl;
		}
	}

	MPI_Bcast(&local_num, 1, MPI_INT, 0, MPI_COMM_WORLD);
	for (int i = 0; i < n_elem; i++)
	{
		MPI_Bcast((int **)&(connectivity[i][0]), 4, MPI_INT, 0, MPI_COMM_WORLD);
	}


	for (int i = 0; i < n_node; i++)
	{
		MPI_Bcast((double **)&(coordinate[i][0]), 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	}

	MPI_Barrier(MPI_COMM_WORLD);


	printf("process %d prainting connectivity:\n", myrank);

	for (int i = 0; i < n_node; i++)
	{
		for (int j = 0; j < col; j++)
			printf("%d \t", connectivity[i][j], myrank);
		printf("\n");
	}


	printf("process %d prainting coordinate:\n", myrank);

	for (int i = 0; i < n_node; i++)
	{
		for (int j = 0; j < col; j++)
			printf("%g \t", coordinate[i][j], myrank);
		printf("\n");
	}

	MPI_Barrier(MPI_COMM_WORLD);


	/***********************************************************/


	for (int i = 0; i < local_num; i++)
	{
		//local_NN[n_node][n_node] = {};
		int node_0 = connectivity[i][1];
		int node_1 = connectivity[i][2];
		int node_2 = connectivity[i][3];

		// element area will drive from detJ= x13*y23- y13*x23, Ae=0.5*|J|


		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];


		local_NN[node_0][node_0] = NN[node_0][node_0] + NN_e[0][0];
		local_NN[node_0][node_1] = NN[node_0][node_1] + NN_e[0][1];
		local_NN[node_0][node_2] = NN[node_0][node_2] + NN_e[0][2];
		local_NN[node_1][node_0] = NN[node_1][node_0] + NN_e[1][0];
		local_NN[node_1][node_1] = NN[node_1][node_1] + NN_e[1][1];
		local_NN[node_1][node_1] = NN[node_1][node_1] + NN_e[1][2];
		local_NN[node_2][node_0] = NN[node_2][node_0] + NN_e[2][0];
		local_NN[node_2][node_1] = NN[node_2][node_1] + NN_e[2][1];
		local_NN[node_2][node_2] = NN[node_2][node_2] + NN_e[2][2];
		 
	}

	/***********************************************************/

	//MPI_Reduce(&local_NN[0][0], &global_NN[0][0], n_node*n_node, MPI_DOUBLE,MPI_SUM, 0, MPI_COMM_WORLD);

	for (int i = 0; i < n_node; i++)
	{
		MPI_Reduce((double **)&(local_NN[i][0]),(double **)&(global_NN[i][0]), n_node, MPI_DOUBLE,MPI_SUM, 0, MPI_COMM_WORLD);
	}

	/****************************************/
	for (int i = 0; i < n_elem; i++)
	{
		delete[] connectivity[i];
	}
	delete[] connectivity;
	/***************************************/
	for (int i = 0; i < n_node; i++)
	{
		delete[] coordinate[i];
	}
	delete[] coordinate;
	/**************************************/
	for (int i = 0; i < n_node; i++)
	{
		delete[] local_NN [i];
	}
	delete[] local_NN;
	/**************************************/

	if (myrank == 0)
	{
		fstream myfile_NN;
		myfile_NN.open("NN.txt", fstream::out);

		myfile_NN << std::endl;

		for (int i = 0; i < n_node; i++)
		{
			for (int j = 0; j < n_node; j++)
			{

				myfile_NN << global_NN[i][j] << "\t";

			}

			myfile_NN << std::endl;
		}
	}

	/*****************************************/
	for (int i = 0; i < n_node; i++)
	{
		delete[] global_NN[i];
	}
	delete[] global_NN;
	/*****************************************/
	MPI_Finalize();
	return 0;
}
 Edit & Run
Posted
Updated 18-May-21 18:51pm

1 solution

Quote:
After assembly I get strange number in global stiffness matrix which is -1.25549e+67 and I don't know why.

Since you didn't provided the input files, it may get complicated to give you specific help.
Advice: change your code to print the input data and intermediate result after each step to follow what is going on.
Or use the dedicated tool the debugger.

Your code do not behave the way you expect, or you don't understand why !

There is an almost universal solution: Run your code on debugger step by step, inspect variables.
The debugger is here to show you what your code is doing and your task is to compare with what it should do.
There is no magic in the debugger, it don't know what your code is supposed to do, it don't find bugs, it just help you to by showing you what is going on. When the code don't do what is expected, you are close to a bug.
To see what your code is doing: Just set a breakpoint and see your code performing, the debugger allow you to execute lines 1 by 1 and to inspect variables as it execute.

Debugger - Wikipedia, the free encyclopedia[^]

Mastering Debugging in Visual Studio 2010 - A Beginner's Guide[^]
Basic Debugging with Visual Studio 2010 - YouTube[^]

1.11 — Debugging your program (stepping and breakpoints) | Learn C++[^]

The debugger is here to only show you what your code is doing and your task is to compare with what it should do.
 
Share this answer
 
Comments
Member 15181211 19-May-21 12:06pm    
Thnaks for your answer. How in MPI you can debug line by line?
Patrice T 19-May-21 13:33pm    
Before going to multi processing, make sure your code is correct with mono processing.
Member 15181211 19-May-21 16:25pm    
The code is working property on single processor. Actually I found the problem now and the code is working. Thanks
Patrice T 19-May-21 17:13pm    
Since you solved the problem, you can post a solution explaining what was the problem and your solution.

This content, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)



CodeProject, 20 Bay Street, 11th Floor Toronto, Ontario, Canada M5J 2N8 +1 (416) 849-8900