MPI for finite element stiffness matrix assembly

Question

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 18-May-21 17:57pm

Member 15181211

Updated 18-May-21 18:51pm

Add a Solution

1 solution

Add a Solution

Add your solution here

Treat my content as plain text, not as HTML

Preview 0

…

Existing Members

Sign in to your account

...or Join us

Download, Vote, Comment, Publish.

Your Email
Password
Forgot your password?

Your Email
This email is in use. Do you need your password?
Optional Password

I have read and agree to the Terms of Service and Privacy Policy
Please subscribe me to the CodeProject newsletters

When answering a question please:

Read the question carefully.
Understand that English isn't everyone's first language so be lenient of bad spelling and grammar.
If a question is poorly phrased then either ask for clarification, ignore it, or edit the question and fix the problem. Insults are not welcome.
Don't tell someone to read the manual. Chances are they have and don't get it. Provide an answer or move on to the next question.

Let's work to help developers, not make them feel stupid.

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

Patrice T · Answer 1 · 2021-05-18T18:51:00

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.