1 /* pdnmesh - a 2D finite element solver
2 Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net>
3 This program is free software; you can redistribute it and/or modify
4 it under the terms of the GNU General Public License as published by
5 the Free Software Foundation; either version 2 of the License, or
6 (at your option) any later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16 $Id: spmat.c,v 1.11 2005/03/27 22:39:59 sarod Exp $
17 */
18
19 /* implementation of sparse solvers
20 using iterative methods = conjugate gradients, arnoldi etc
21 */
22
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26
27
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <math.h>
31 #include "types.h"
32
33
34 /* a function to solve a matrix equation by conjugate gradient */
35 /* method Ax=y , A - symmetric positive definite */
36 /* iter = no. of iterations, iter > N, N = size of A */
37 static void
solve_conjugate_sparse(SPMat * A,MY_DOUBLE * x,MY_DOUBLE * y,int iter)38 solve_conjugate_sparse( SPMat *A, MY_DOUBLE *x, MY_DOUBLE *y, int iter) {
39 MY_DOUBLE *r_old, *p_old, *v;
40 MY_DOUBLE a,b,alph,beta;
41 MY_INT i,k;
42 MY_INT N;
43
44 N=SPM_dim(A);
45 /* important step: build linked lists for speed
46 before calling SPM_vec() */
47 SPM_build_lists(A);
48
49 #ifdef DEBUG
50 int j;
51 for ( i =0; i < N ; i++ ) {
52 for ( j = 0 ; j <N; j++ )
53 printf(MDF" ",SPM_e(A,i,j));
54 printf("\n");
55 }
56 #endif
57 /* allocate memory for arrays */
58 if ((r_old= (MY_DOUBLE*)calloc(N,sizeof(MY_DOUBLE)))==0){
59 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
60 exit(1);
61 }
62 if ((p_old= (MY_DOUBLE*)calloc(N,sizeof(MY_DOUBLE)))==0){
63 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
64 exit(1);
65 }
66 if ((v= (MY_DOUBLE*)calloc(N,sizeof(MY_DOUBLE)))==0){
67 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
68 exit(1);
69 }
70
71 /* intialization */
72 /* x_old = 0 */
73 for (i=0; i<N; i++) {
74 p_old[i]=r_old[i]=y[i];
75 x[i]=0;
76 }
77 /* iteration make iteration approx size of matrix */
78 for (k=0; k < iter; k++) {
79 /* calculate Ap_old =>v*/
80 SPM_vec(A, p_old, v);
81 a=b=0;
82 for (i=0; i<N; i++) {
83 a+=r_old[i]*r_old[i];
84 b+=p_old[i]*v[i];
85 }
86 alph=a/b;
87 b=0;
88 for (i=0; i<N; i++) {
89 x[i]=x[i]+alph*p_old[i];
90 r_old[i]=r_old[i]-alph*v[i];
91 b+=r_old[i]*r_old[i];
92 }
93 beta=b/a;
94 for (i=0; i<N; i++) {
95 p_old[i]=r_old[i]+beta*p_old[i];
96 }
97
98 #ifdef DEBUG
99 printf("conj grad iteration %d\n",k);
100 for (i=0; i<N; i++) {
101 printf(MDF" ",x[i]);
102 }
103 printf("\n");
104 #endif
105 }
106
107 free(r_old);
108 free(p_old);
109 free(v);
110 }
111
112 /* a function to build the global matrix or the matrix */
113 /* equation for the system - Poisson Equation */
114 static void
build_sparse(SPMat * P,MY_DOUBLE * q,MY_DOUBLE * x,MY_DOUBLE * y,MY_DOUBLE * phi,int NUN,int NE,int NN,MY_INT * renumber,MY_INT * re_renumber)115 build_sparse( SPMat *P, MY_DOUBLE *q, MY_DOUBLE *x, MY_DOUBLE *y, MY_DOUBLE *phi, int NUN, int NE, int NN, MY_INT *renumber, MY_INT *re_renumber)
116 {
117 /* infd - input file descriptor */
118 /* P[][] - global matrix to build NUNxNUN */
119 /* q[] - global array to build NUNx1 */
120 /* x[] - x coordinates NNx1 */
121 /* y[] - y coordinates NNx1 */
122 /* phi[] - potentials NNx1 */
123 /* NUN - no of unknown potential nodes */
124 /* NE - no of elements */
125 /* NN - no of nodes */
126
127 /* now the running variables for each triangular element */
128 int v[3]; /* to keep vertex number of each node */
129 MY_DOUBLE pl[3][3]; /* local version of P[][] above */
130 MY_DOUBLE ql[3]; /* local version of q[] above */
131 MY_DOUBLE rho[3]; /* rho values */ /* conductivity */
132 triangle *t;
133
134 int i,j;
135
136 MY_DOUBLE
137 epsilon, mu, /* permittivity, permeability*/
138 A, /* area of a triangulr element */
139 k, /* value of the determinanat */
140 yy0, yy1, yy2;
141
142 MY_DOUBLE b[3]; /*local arrays */
143 MY_DOUBLE c[3];
144
145
146
147 update_status_bar("Begin building...");
148 /* now read node data with unknown potentials */
149 for ( i = 0; i < NUN; i++) {
150 x[i]= (Mx(re_renumber[i],M))*g_xscale-g_xoff;
151 y[i]= (My(re_renumber[i],M))*g_yscale-g_yoff;
152 }
153
154 /* now read node data with known potentials , which comes last */
155 for ( i = NUN; i < NN; i++ ) {
156 x[i] = (Mx(re_renumber[i],M))*g_xscale-g_xoff;
157 y[i] = (My(re_renumber[i],M))*g_yscale-g_yoff;
158 phi[i] = (Mz(re_renumber[i],M));
159 }
160
161 /* now initialize P[][] and q[] */
162 /* for ( i = 0; i < NUN; i++) {
163 q[i] = 0;
164 }*/
165
166 #ifdef DEBUG
167 for ( i = 0; i< NN; i++)
168 printf("build: %d "MDF", "MDF": "MDF"\n", i, x[i], y[i], phi[i]);
169 #endif
170
171 /* now read element data, whilst building the global matrices */
172 DAG_traverse_list_reset(&dt);
173 t=DAG_traverse_prune_list(&dt);
174 while(t) {
175 if ( t &&t->status==LIVE ) {
176
177 v[0] = renumber[t->p0];
178 v[1] = renumber[t->p1];
179 v[2] = renumber[t->p2]; /* ccw direction */
180 if (bound.poly_array[t->boundary].rhoexp) {
181 rho[0] = ex(bound.poly_array[t->boundary].rhoexp,t->p0);
182 rho[1] = ex(bound.poly_array[t->boundary].rhoexp,t->p1);
183 rho[2] = ex(bound.poly_array[t->boundary].rhoexp,t->p2);
184 } else {
185 rho[0] = bound.poly_array[t->boundary].rho;
186 rho[1] = bound.poly_array[t->boundary].rho;
187 rho[2] = bound.poly_array[t->boundary].rho;
188 }
189
190 /* for Poisson equation, calculate the ration epsilon/mu
191 (eps/mu) del^2 phi = -rho
192 */
193 if(bound.poly_array[t->boundary].muexp) {
194 mu=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].muexp,t->p0)
195 +ex(bound.poly_array[t->boundary].muexp,t->p1)
196 +ex(bound.poly_array[t->boundary].muexp,t->p2));
197 } else {
198 mu=bound.poly_array[t->boundary].mu;
199 }
200 if(bound.poly_array[t->boundary].epsexp) {
201 epsilon=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].epsexp,t->p0)
202 +ex(bound.poly_array[t->boundary].epsexp,t->p1)
203 +ex(bound.poly_array[t->boundary].epsexp,t->p2));
204 } else {
205 epsilon=bound.poly_array[t->boundary].epsilon;
206 }
207 /* find the ratio */
208 epsilon /=mu;
209
210
211 /* now have read one triangle , add it to global */
212 /* first the coefficients */
213 /* some temporary values, neede several times */
214 yy0 = y[v[1]]-y[v[2]];
215 yy1 = y[v[2]]-y[v[0]];
216 yy2 = y[v[0]]-y[v[1]];
217
218 /* area of the triangle */
219 A = 0.5*( x[v[0]]*yy0 + x[v[1]]*yy1 + x[v[2]]*yy2 );
220 /* A = ABS(A); */ /* points are in ccw order to ensure this */
221 /* another constant */
222 k = 0.25/A;
223 /* now calculate b's and c's */
224 #ifdef DEBUG
225 printf("adding element with nodes %d %d and %d...\n\n", v[0], v[1], v[2]);
226 #endif
227 b[0]= yy0;
228 #ifdef DEBUG
229 printf("b0="MDF" ", b[0]);
230 #endif
231 b[1]= yy1;
232 #ifdef DEBUG
233 printf("b1="MDF" ", b[1]);
234 #endif
235 b[2]= yy2;
236 #ifdef DEBUG
237 printf("b2="MDF" \n", b[2]);
238 #endif
239 c[0] = ( x[v[2]]-x[v[1]] );
240 #ifdef DEBUG
241 printf("c0="MDF" ", c[0]);
242 #endif
243 c[1] = ( x[v[0]]-x[v[2]] );
244 #ifdef DEBUG
245 printf("c1="MDF" ", c[1]);
246 #endif
247 c[2] = ( x[v[1]]-x[v[0]] );
248 #ifdef DEBUG
249 printf("c2="MDF" \n", c[2]);
250 #endif
251
252 /* now build pl[] and ql */
253 for ( i = 0; i < 3; i++ )
254 for ( j = 0; j < 3; j++ )
255 pl[i][j] = k*epsilon*(b[i]*b[j]+c[i]*c[j]);
256
257 ql[0] = A*(2*rho[0]+rho[1]+rho[2])*ONE_OVER_TWELVE;
258 ql[1] = A*(rho[0]+2*rho[1]+rho[2])*ONE_OVER_TWELVE;
259 ql[2] = A*(rho[0]+rho[1]+2*rho[2])*ONE_OVER_TWELVE;
260 #ifdef DEBUG
261 printf("local p and q ..\n");
262 for ( i = 0; i <3; i++) {
263 printf("| ");
264 for ( j = 0; j < 3; j++)
265 printf(" "MDF" ", pl[i][j]);
266 printf("|| "MDF" | \n", ql[i]);
267 }
268 #endif
269
270
271 /* better to use multiplication than division */
272 /* now augment the global data set */
273 for ( i = 0; i < 3; i++)
274 if ( v[i] < NUN ) { /* we index from 0, not from 1 */
275 q[v[i]] = q[v[i]] + ql[i];
276 for ( j = 0; j < 3; j++ ) {
277 if ( v[j] < NUN )
278 {
279 if ( v[j] <= v[i] )
280 SPM_add(P,v[i],v[j],pl[i][j]);
281 SPM_eq(P,v[j],v[i],SPM_e(P,v[i],v[j]));
282 /* we build both halfs of the matrix:FIXME: this can be done later */
283 }
284 else
285 {
286 q[v[i]] = q[v[i]] - (pl[i][j])*phi[v[j]];
287 }
288 }
289 } /*else printf("rejecting %d\n", v[i]);*/
290 #ifdef DEBUG
291 printf("building global matrix ..\n");
292 #endif
293 }
294 t=DAG_traverse_prune_list(&dt);
295 }
296
297 #ifdef DEBUG
298 for ( i = 0; i < NUN; i++){
299 printf("| ");
300 for( j = 0; j <= i; j++)
301 printf(" "MDF" ", SPM_e(P,i,j));
302 printf("= "MDF" |\n", q[i]);
303 }
304 for ( i = 0; i < NN; i++)
305 printf("node %d potential "MDF"\n", i, phi[i]);
306 #endif
307 update_status_bar("Done building.");
308 }
309
310 /* after genarating the mesh, need to re-number points */
311 /* so that unknown points are first */
312 /* then we solve the problem -poisson equation*/
313 int
re_number_and_solve_poisson_sparse(void)314 re_number_and_solve_poisson_sparse(void)
315 {
316 unsigned MY_INT i,j;
317 MY_INT * renumber; /* re-numbering array */
318 MY_INT * re_renumber; /* re-renumbering array */
319 /* used to map potentials back to the original points */
320
321
322 unsigned MY_INT unknowns=0; /* unknown points */
323 unsigned MY_INT total = 0; /* total points */
324 MY_INT triangles = 0; /* known points */
325
326
327 triangle *t;
328
329 /* first the global node data, all total size */
330 MY_DOUBLE *x ; /* x coordinates */
331 MY_DOUBLE *y ; /* y coordinates */
332 MY_DOUBLE *phi ; /* potentials at each node */
333
334 /* the global variables */
335 SPMat P ; /* sparse NN x NN matrix, only the unknowns */
336 MY_DOUBLE *q ; /* unknowns sizesame here */
337
338 #ifdef DEBUG
339 printf("starting sparse solver\n");
340 #endif
341
342 /* allocate memory for arrays */
343 if ((renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT )))==0){
344 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
345 exit(1);
346 }
347
348 if((re_renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
349 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
350 exit(1);
351 }
352
353 /* need to solve for only points referred by triangles */
354 /* SO -- */
355 DAG_traverse_list_reset(&dt);
356 t=DAG_traverse_prune_list(&dt);
357 while ( t ) {
358 if ( t && t->status==LIVE ) {
359 /* mark the points in this triangle */
360 /* use re_renumber[] array for this */
361 re_renumber[t->p0] = 1;
362 re_renumber[t->p1] = 1;
363 re_renumber[t->p2] = 1;
364
365 triangles++;
366 }
367 t=DAG_traverse_prune_list(&dt);
368 }
369
370 /* now find unknowns, or INSIDE points */
371 j = 0;
372 for ( i = 0; i < M.count; i++ )
373 if ( Mval(i,M) == VAR && re_renumber[i] == 1) {
374 renumber[i] = j++;
375 #ifdef DEBUG
376 printf("renumber_and_solve: renumbering variable point %d as %d\n",i,renumber[i]);
377 #endif
378 unknowns++;
379 }
380 /* now renumber the known points */
381 total = unknowns;
382 for ( i = 0; i < M.count; i++ )
383 if ( Mval(i,M) == FX && re_renumber[i]==1 ) {
384 renumber[i] = j++;
385 #ifdef DEBUG
386 printf("renumber_and_solve: renumbering fixed point %d as %d\n",i,renumber[i]);
387 #endif
388 total++;
389 }
390 /* finally give numbers to points 0,1,2 */
391 /* we will not use these points */
392 renumber[0]=j++;
393 renumber[1]=j++;
394 renumber[2]=j++;
395
396 /* now construct re-renmbering array */
397 i = 0;
398 /* do an insertion sort */
399 while ( i < M.count ) {
400 for ( j = 0; j < M.count; j++ ) {
401 if ( renumber[j] ==(int) i )
402 break;
403 }
404 re_renumber[i] = j;
405 i++;
406 }
407
408 #ifdef DEBUG
409 printf("build: renumbering %d points with %d unknowns out of %d total\n",total,unknowns,M.count);
410 #endif
411 /*renumber_points(renumber, re_renumber, total, unknowns); */
412
413 #ifdef DEBUG
414 printf("i, renum, rerenum\n");
415 for ( i = 0; i < total; i++ )
416 printf("%d %d %d\n",i,renumber[i],re_renumber[i]);
417 for ( i = 0; i < total; i++ )
418 printf("%d "MDF" "MDF" "MDF"\n",i,Mx(i,M),My(i,M),Mz(i,M));
419 #endif
420
421 /* we free these two arrays only after solving the problem */
422 /* now solve the problem */
423 /* first allocate memory for the arrays */
424 if((x=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
425 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
426 exit(1);
427 }
428
429 if((y=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
430 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
431 exit(1);
432 }
433
434 if((phi=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
435 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
436 exit(1);
437 }
438
439 if((q=(MY_DOUBLE*)calloc(unknowns, sizeof(MY_DOUBLE)))==0){
440 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
441 exit(1);
442 }
443
444
445 /* now the matrix */
446 SPM_init(&P,unknowns);
447
448 /* Poisson equation */
449 /* now the memory allocation is over */
450 build_sparse( &P, q, x, y, phi, unknowns, triangles, total, renumber, re_renumber);
451 /* if user has chosen so, solve by conj. grad method*/
452 /* by default solve by cholevsky decomposition */
453 solve_conjugate_sparse(&P,phi,q,unknowns+10);
454
455 /* after solution, array phi[] will contain the new potentials */
456 /* now put them back to the point data */
457 /* we also note the maximum and minimum values for contour plotting */
458 g_maxpot = -1000;
459 g_minpot = 1000;
460 for ( i = 0; i < total; i++ ) {
461 Mz(re_renumber[i],M) = phi[i];
462 if ( phi[i] > g_maxpot )
463 g_maxpot = phi[i];
464 if ( phi[i] < g_minpot )
465 g_minpot = phi[i];
466 }
467
468
469
470 /* now we can free the renumber[] array */
471 free(renumber);
472 /* now free this array too */
473 free(re_renumber);
474
475 /* free everything else */
476
477 /* first the arrays */
478 free(x);
479 free(y);
480 free(phi);
481 free(q);
482
483 /* next the matrix */
484 SPM_destroy(&P);
485 /* now printout a summery */
486 #ifdef DEBUG
487 for ( i = 0; i < total; i++ )
488 printf("renumber_and_solve: %d "MDF" "MDF" "MDF"\n",i,Mx(i,M),My(i,M),Mz(i,M));
489 #endif
490
491 return(0);
492 }
493