Matrix assembly in fem using mpi

following my previous question regarding using MPI for finite element
http://www.cplusplus.com/forum/general/278132/ I have made a few changes in my code. 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. 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.




1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
		


#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;
}





Last edited on
@resabzr
The whole point of using distributed-memory computing is that you don't send everything to all processors.

Also, you have already been told (several times) on StackOverflow to use contiguous memory - not your (current) **double pointer arrangement. The latter will hopelessly slow down any MPI transfers. Either flatten your dynamic arrays to 1-d, or use the data buffer of std::vector.

If you absolutely must access elements in a 2-d array manner then, for contiguous memory you will have to allocate memory as @dutch does in this post:
http://www.cplusplus.com/forum/beginner/267829/#msg1152277
Note that this is essentially a flattened 1-d array "under the hood".

You shouldn't use MPI_Bcast node-by-node; that would be horribly slow. You also can't use whole-array MPI transfers unless your data is in contiguous memory, which it won't be with your (current) double-pointer arrangement. Similarly, I wouldn't guarantee that your MPI_Reduce will work (though I can't see what it would be used for in this code, anyway).

FWIW you have a data indexing error in two of the terms on line 249:
local_NN[node_1][node_1] = NN[node_1][node_1] + NN_e[1][2];
but there's no way of knowing whether this is the cause of your error, because there is no way of testing your code on an input file.

If you want to learn how to program the finite-element method then you can read
"Programming the Finite Element Method" by Smith, Griffiths and Margetts:
https://onlinelibrary.wiley.com/doi/epub/10.1002/9781119189237
and yeah, two of the authors are my (retired or current) work colleagues.
Last edited on
I do appreciate your answer. But I think you are mistaken. I have NOT been told in stuckoverflow several times about 1d contiguous array. anyway, I need the 2d array in advance. If you look at my previous code I defined 1d contiguous array, but I needed a 2d array, that's why I am using double pointer. About mpi reduce, imagine we have 2 matrix in 2 process,

A={{1,2},{0,0}};
B={{0,0},{1,2}};
I want mpi reduce do this for me,
A+B={{1,2},{1,2}}
Last edited on
@resabzr
No, you do not understand.

Instead of a 2-d array, say
A = [ [ 1, 2 ], [ 3, 4 ] ]
you can simply "flatten" this to
A = [ 1, 2, 3, 4 ]
and if you want the equivalent 2-d indices then, as you know the number of columns, then you can easily work them out. But the key advantage is that the latter form can be MPI-transferred in a SINGLE communication - you don't have to do it row by row. MPI communications are slow. Imagine what would happen if you were working in 3-d, not 2-d, with your multiple-pointer notation.

I can see that your MPI_Reduce will (in this code) do a sum. However, at the moment, all the processors seem to know everything, so all you will be doing is summing the same thing up lots of times. You haven't decomposed your domain.

I assume that the index problem that I pointed out didn't solve the problem? Does the problem still occur if you just use a single processor?



resabzr wrote:
But I think you are mistaken. I have NOT been told in stuckoverflow several times about 1d contiguous array

I had a look at the various answers you've had in Stack Overflow (to the same codes that you have posted here). Those answers include these gems:
all these raw poiters in modern C++ make my head spin...

The matrices are not contiguous

At the very least, you need to allocate your 2D array(s) in contiguous memory. 

Assuming your matrix is in contiguous memory, all you need is MPI_Reduce(..., MPI_SUM, ...)



BTW, the term is "halo data", not "hole data" (this thread).
Likewise, "contiguous", not "contagious" (one of your SO posts).
Last edited on
I got your point about single block of memory, I understand this
And I used this method http://www.cplusplus.com/forum/beginner/267829/#msg1152277.
I fixed that index problem but it did not make any change in the results. I wrote the same code in c++ and I'm getting the right answer, but here even with 1 process everything is different. Even when I print out the area of each element it gives me zero, or when I print out the gradient of shape functions it gives me non(ind). How it can be possible even on one process results change like this in comparison with c++ code without mpi implementation.
Supply:
(1) your latest code (below this post)
(2) the input file (preferably small)

Make sure that any hard-coded cell counts correspond to the numbers in the input file.
this is the final code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
#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 gx1;
double gy1;
double gx2;
double gy2;
double gx3;
double gy3;

double x[n_node];
double y[n_node];
double elem_area[n_elem];

double RHS[n_node][1];


using namespace std;

void assembly(int n, int **connectivity, double **coordinate);
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;

	// allocate
	int **connectivity = new int*[n_elem];     
	connectivity[0] = new int[n_elem*4];   
	for (int r = 1; r < n_elem; ++r) 
		connectivity[r] = connectivity[r - 1] + col;

	// allocate
	double **coordinate = new double*[n_node];  
	coordinate[0] = new double[n_node * 4];   
	for (int r = 1; r < n_node; ++r) 
		coordinate[r] = coordinate[r - 1] + col;


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


	double **global_NN = new double*[n_node];    
	global_NN[0] = new double[n_node * n_node];   
	for (int r = 1; r < n_node; ++r) 
		global_NN[r] = global_NN[r - 1] + 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);
	MPI_Bcast(&(connectivity[0]), n_elem * 4, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&(coordinate[0]), n_node * 4, MPI_DOUBLE, 0, 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");
	}



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

	assembly(local_num, connectivity, coordinate);
	
	/***********************************************************/
	if (myrank == 0)
	{
		ofstream fout;
		fout.open("elem_area.txt");
		fout << "variables = elem_area " << endl;

		fout << "zone I = " << n_elem << endl;
		for (int i = 0; i < n_elem; i++)
		{
			fout << i << "       " << elem_area[i] << ' ';
			fout << endl;
		}
		fout.close();

		fout.open("gradient_x.txt");
		fout << "variables = Nx " << endl;

		fout << "zone I = " << n_elem << " " << "zone J = " << 3 << endl;
		for (int i = 0; i < n_elem; i++)
		{

			fout << i << "    " << Nx[i][0] << "    " << Nx[i][1] << "    " << Nx[i][2] << ' ';
			fout << endl;
		}
		fout.close();

		fout.open("gradient_y.txt");
		fout << "variables = Ny " << endl;

		fout << "zone I = " << n_elem << " " << "zone J = " << 3 << endl;
		for (int i = 0; i < n_elem; i++)
		{

			fout << i << "    " << Ny[i][0] << "    " << Ny[i][1] << "    " << Ny[i][2] << ' ';
			fout << endl;
		}
		fout.close();

	}

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

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

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

	delete[] connectivity[0];
	delete[] connectivity;   


	delete[] coordinate[0];
	delete[] coordinate;  

	delete[] local_NN[0]; 
	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;
		}
	}


	delete[] global_NN[0];
	delete[] global_NN;    

	MPI_Finalize();
	return 0;
}


void assembly(int n, int **connectivity, double **coordinate)
{

	int num_procs, myrank;
	MPI_Status status;
	MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);


	double elem_area[n_elem];
	double a, b, c, d;

	double **local_NN = new double*[n_node];      // allocate row pointers
	local_NN[0] = new double[n_node * n_node];   // allocate data
	for (int r = 1; r < n_node; ++r) // set row pointers
		local_NN[r] = local_NN[r - 1] + n_node;

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

		// 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_2] = 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];

	}
}




Last edited on
and this the mesh and connectivity

MESH dimension 3 ElemType Triangle Nnode 3
Coordinates

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 1               2               0               0
    2             1.5               0               0
    3               2             0.5               0
    4      1.28472222     0.501860758               0
    5               1               0               0
    6               2               1               0
    7             1.5               1               0
    8     0.708333333      0.51116455               0
    9               1               1               0
   10             0.5               0               0
   11             0.5               1               0
   12               0               0               0
   13               0             0.5               0
   14               0               1               0
End Coordinates

Elements
1 12 10 13
2 1 3 2
3 9 11 8
4 8 11 13
5 9 8 4
6 4 8 5
7 9 4 7
8 7 4 3
9 5 2 4
10 6 7 3
11 14 13 11
12 10 5 8
13 4 2 3
14 10 8 13
End Elements
Last edited on
What exactly do your files connec.txt and mesh.txt contain. These are the files being read by your program.

Note, you should remove the lines
1
2
	while (!file.eof())
		{

and the corresponding end brace, or the code will attempt (and fail) to read your coordinates a second time.
Last edited on
the data I put here as coordinate and element actually are two separate file named connec.txt and mesh.txt. I just copy them here so you can see the data, I didn't know how put files here. I just did what you said but no difference.
Please just post those files. As is.
Last edited on
Just cut and paste it from notepad or whatever.
Add outpad tags around it (second item on the format menu).
OK, change the lines
1
2
	MPI_Bcast(&(connectivity[0]), n_elem * 4, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&(coordinate[0]), n_node * 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);

to
1
2
	MPI_Bcast(&(connectivity[0][0]), n_elem * 4, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&(coordinate[0][0]), n_node * 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);


Make sure that the number of elements corresponds to that in the input file.
Last edited on
not working, I did these changes but not working. I think the problem start when area of each element yield zero, when I print the element area it's all zero
Last edited on
The areas will come out 0 because you don't set arrays x[] or y[]. Did you mean something to do with coordinate[][]?


With two processors your code runs and produces (from your 14-element files):

C:\c++>"C:\Program Files\Microsoft MPI\bin"\mpiexec -n 2 temp 
0	11	9	12	
1	0	2	1	
2	8	10	7	
3	7	10	12	
4	8	7	3	
5	3	7	4	
6	8	3	6	
7	6	3	2	
8	4	1	3	
9	5	6	2	
10	13	12	10	
11	9	4	7	
12	3	1	2	
13	9	7	12	
0	2	0	0	
1	1.5	0	0	
2	2	0.5	0	
3	1.28472	0.501861	0	
4	1	0	0	
5	2	1	0	
6	1.5	1	0	
7	0.708333	0.511165	0	
8	1	1	0	
9	0.5	0	0	
10	0.5	1	0	
11	0	0	0	
12	0	0.5	0	
13	0	1	0	
process 0 prainting connectivity:
0 	11 	9 	12 	
1 	0 	2 	1 	
2 	8 	10 	7 	
3 	7 	10 	12 	
4 	8 	7 	3 	
5 	3 	7 	4 	
6 	8 	3 	6 	
7 	6 	3 	2 	
8 	4 	1 	3 	
9 	5 	6 	2 	
10 	13 	12 	10 	
11 	9 	4 	7 	
12 	3 	1 	2 	
13 	9 	7 	12 	
process 0 prainting coordinate:
0 	2 	0 	0 	
1 	1.5 	0 	0 	
2 	2 	0.5 	0 	
3 	1.28472 	0.501861 	0 	
4 	1 	0 	0 	
5 	2 	1 	0 	
6 	1.5 	1 	0 	
7 	0.708333 	0.511165 	0 	
8 	1 	1 	0 	
9 	0.5 	0 	0 	
10 	0.5 	1 	0 	
11 	0 	0 	0 	
12 	0 	0.5 	0 	
13 	0 	1 	0 	
process 1 prainting connectivity:
0 	11 	9 	12 	
1 	0 	2 	1 	
2 	8 	10 	7 	
3 	7 	10 	12 	
4 	8 	7 	3 	
5 	3 	7 	4 	
6 	8 	3 	6 	
7 	6 	3 	2 	
8 	4 	1 	3 	
9 	5 	6 	2 	
10 	13 	12 	10 	
11 	9 	4 	7 	
12 	3 	1 	2 	
13 	9 	7 	12 	
process 1 prainting coordinate:
0 	2 	0 	0 	
1 	1.5 	0 	0 	
2 	2 	0.5 	0 	
3 	1.28472 	0.501861 	0 	
4 	1 	0 	0 	
5 	2 	1 	0 	
6 	1.5 	1 	0 	
7 	0.708333 	0.511165 	0 	
8 	1 	1 	0 	
9 	0.5 	0 	0 	
10 	0.5 	1 	0 	
11 	0 	0 	0 	
12 	0 	0.5 	0 	
13 	0 	1 	0 	


What other output do you expect?

Note that some of your printf statements have the wrong number of parameters.
Last edited on
The different number is because I changed the mesh to post here. It was not 14 points, it is 246 as you can see from the link I've shared. The problem is not the Broadcasting part. It works properly. The problem is in matrix assembly I get strange number. like this -6.27744e+66. All the elements area are zero, all the shape function gradient are non(ind). and I don't know why. The same code without MPI is working properly
Last edited on
Your number of elements (hard-coded) must match that in the input file.

Your code now runs.

Areas will work out as 0 (or, possibly, garbage) because you have not put anything in arrays x[] and y[].
many thanks for your help. the problem was because of x and y.
@lastchance
could you please take a look at this, https://www.cplusplus.com/forum/general/278531/
Last edited on
Registered users can post here. Sign in or register to post.