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: solve.c,v 1.44 2005/04/24 05:32:24 sarod Exp $
17 */
18
19
20
21 #ifdef HAVE_CONFIG_H
22 #include <config.h>
23 #endif
24
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <math.h> /* for sqrt() */
28 #include <string.h> /* for memcpy */
29 #include "types.h"
30
31
32 MY_DOUBLE g_maxpot,g_minpot;
33
34 /* type of equation POISSON, HELMHOLTZ */
35 eq_type solve_equation;
36 int use_arpack_solver=0;
37
38 /* a function to solve a matrix equation by symmetric LU */
39 /* decomposition or cholevsky decomposition */
40 static void
solve_cholevsky(MY_DOUBLE ** a,MY_DOUBLE * x,MY_DOUBLE * y,int N)41 solve_cholevsky( MY_DOUBLE **a, MY_DOUBLE *x, MY_DOUBLE *y, int N)
42 {
43 /* solves a system in the form [a][x] = [y] order N where [a] is symmetric*/
44 /* first cholevsky decomposition */
45 int i,j,k;
46 MY_DOUBLE temp;
47
48 #ifdef DEBUG
49 for ( i = 0; i < N ; i++ ) {
50 for ( j = 0 ; j <= i; j++ )
51 printf(MDF" ",a[i][j]);
52 printf("\n");
53 }
54 #endif
55
56 update_status_bar("Begin solving...");
57 for ( j = 0; j < N; j++) {
58 /* first diagonal element */
59 temp = 0;
60 for ( k = 0; k < j; k++)
61 temp += a[j][k]*a[j][k];
62 a[j][j] = sqrtf((float)( a[j][j] - temp) );
63 #ifdef DEBUG
64 printf("solve_chol: col %d diag "MDF"\n",j,a[j][j]);
65 #endif
66 /* check for singularities */
67 if ( a[j][j] <= 0.0f ) {
68 a[j][j] = TOL; /* a small value */
69 fprintf(stderr,"%s: %d: warning: solve_cholevsky: singular matrix\n",__FILE__,__LINE__);
70 }
71 /* now handle the i th column */
72 for( i = j+1; i < N; i++) {
73 temp = 0;
74 for ( k = 0; k < j; k++)
75 temp += a[i][k]*a[j][k];
76 a[i][j] = ( a[i][j] - temp )/a[j][j];
77 #ifdef DEBUG
78 printf("solve_chol: row %d elem "MDF"\n",i,a[i][j]);
79 #endif
80 }
81 }
82
83 #ifdef DEBUG
84 for ( i =0; i < N ; i++ ) {
85 for ( j = 0 ; j <= i; j++ )
86 printf(MDF" ",a[i][j]);
87 printf("\n");
88 }
89 #endif
90
91 /* forward elimination */
92 for ( i = 0; i < N; i++) {
93 temp = 0;
94 for ( k = 0; k < i; k++)
95 temp += a[i][k]*y[k];
96 /* since we have eliminated singularities, no need to check a[i][i]*/
97 y[i] = ( y[i] - temp )/a[i][i];
98 }
99
100 /* backward substitution */
101 for ( i = N-1; i >= 0; i--) {
102 temp = 0;
103 for ( j = i+1; j< N; j++)
104 temp += x[j]*a[j][i];
105 /* since we have eliminated singularities, no need to check a[i][i]*/
106 x[i] = ( y[i] - temp )/a[i][i];
107 }
108
109 update_status_bar("Done solving.");
110 }
111
112
113 /* a function to build the global matrix or the matrix */
114 /* equation for the system - Poisson Equation */
115 static void
build(MY_DOUBLE ** 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 build( MY_DOUBLE **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)
117 {
118 /* infd - input file descriptor */
119 /* P[][] - global matrix to build NUNxNUN */
120 /* q[] - global array to build NUNx1 */
121 /* x[] - x coordinates NNx1 */
122 /* y[] - y coordinates NNx1 */
123 /* phi[] - potentials NNx1 */
124 /* NUN - no of unknown potential nodes */
125 /* NE - no of elements */
126 /* NN - no of nodes */
127
128 /* now the running variables for each triangular element */
129 int v[3]; /* to keep vertex number of each node */
130 MY_DOUBLE pl[3][3]; /* local version of P[][] above */
131 MY_DOUBLE ql[3]; /* local version of q[] above */
132 MY_DOUBLE rho[3]; /* rho values */ /* conductivity */
133 triangle *t;
134
135 int i,j;
136
137 MY_DOUBLE
138 epsilon, mu, /* permittivity, permeability */
139 A, /* area of a triangulr element */
140 k, /* value of the determinanat */
141 yy0, yy1, yy2;
142
143 MY_DOUBLE b[3]; /*local arrays */
144 MY_DOUBLE c[3];
145
146
147
148 update_status_bar("Begin building...");
149 /* now read node data with unknown potentials */
150 for ( i = 0; i < NUN; i++) {
151 x[i]= (Mx(re_renumber[i],M))*g_xscale-g_xoff;
152 y[i]= (My(re_renumber[i],M))*g_yscale-g_yoff;
153 }
154
155 /* now read node data with known potentials , which comes last */
156 for ( i = NUN; i < NN; i++ ) {
157 x[i] = (Mx(re_renumber[i],M))*g_xscale-g_xoff;
158 y[i] = (My(re_renumber[i],M))*g_yscale-g_yoff;
159 phi[i] = (Mz(re_renumber[i],M));
160 }
161
162
163 #ifdef DEBUG
164 for ( i = 0; i< NN; i++)
165 printf("build: %d "MDF", "MDF": "MDF"\n", i, x[i], y[i], phi[i]);
166 #endif
167
168 /* now read element data, whilst building the global matrices */
169 DAG_traverse_list_reset(&dt);
170 t=DAG_traverse_prune_list(&dt);
171 while(t) {
172 if ( t &&t->status==LIVE ) {
173 #ifdef DEBUG
174 printf("%s, %s, %s\n",
175 bound.poly_array[t->boundary].mustr,
176 bound.poly_array[t->boundary].epsstr,
177 bound.poly_array[t->boundary].rhostr);
178 #endif
179
180 v[0] = renumber[t->p0];
181 v[1] = renumber[t->p1];
182 v[2] = renumber[t->p2]; /* ccw direction */
183 if (bound.poly_array[t->boundary].rhoexp) {
184 rho[0] = ex(bound.poly_array[t->boundary].rhoexp,t->p0);
185 rho[1] = ex(bound.poly_array[t->boundary].rhoexp,t->p1);
186 rho[2] = ex(bound.poly_array[t->boundary].rhoexp,t->p2);
187 } else {
188 rho[0] = bound.poly_array[t->boundary].rho;
189 rho[1] = bound.poly_array[t->boundary].rho;
190 rho[2] = bound.poly_array[t->boundary].rho;
191 }
192 /* for Poisson equation, calculate the ration epsilon/mu
193 (eps/mu) del^2 phi = -rho
194 */
195 if(bound.poly_array[t->boundary].muexp) {
196 mu=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].muexp,t->p0)
197 +ex(bound.poly_array[t->boundary].muexp,t->p1)
198 +ex(bound.poly_array[t->boundary].muexp,t->p2));
199 } else {
200 mu=bound.poly_array[t->boundary].mu;
201 }
202 if(bound.poly_array[t->boundary].epsexp) {
203 epsilon=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].epsexp,t->p0)
204 +ex(bound.poly_array[t->boundary].epsexp,t->p1)
205 +ex(bound.poly_array[t->boundary].epsexp,t->p2));
206 } else {
207 epsilon=bound.poly_array[t->boundary].epsilon;
208 }
209 /* find the ratio */
210 epsilon /=mu;
211
212 /* now have read one triangle , add it to global */
213 /* first the coefficients */
214 /* some temporary values, need several times */
215 yy0 = y[v[1]]-y[v[2]];
216 yy1 = y[v[2]]-y[v[0]];
217 yy2 = y[v[0]]-y[v[1]];
218
219 /* area of the triangle */
220 A = 0.5*( x[v[0]]*yy0 + x[v[1]]*yy1 + x[v[2]]*yy2 );
221 /* A = ABS(A); */ /* points are in ccw order to ensure this */
222 /* another constant */
223 k = 0.25/A;
224 /* now calculate b's and c's */
225 #ifdef DEBUG
226 printf("adding element with nodes %d %d and %d...\n\n", v[0], v[1], v[2]);
227 #endif
228 b[0]= yy0;
229 #ifdef DEBUG
230 printf("b0="MDF" ", b[0]);
231 #endif
232 b[1]= yy1;
233 #ifdef DEBUG
234 printf("b1="MDF" ", b[1]);
235 #endif
236 b[2]= yy2;
237 #ifdef DEBUG
238 printf("b2="MDF" \n", b[2]);
239 #endif
240 c[0] = ( x[v[2]]-x[v[1]] );
241 #ifdef DEBUG
242 printf("c0="MDF" ", c[0]);
243 #endif
244 c[1] = ( x[v[0]]-x[v[2]] );
245 #ifdef DEBUG
246 printf("c1="MDF" ", c[1]);
247 #endif
248 c[2] = ( x[v[1]]-x[v[0]] );
249 #ifdef DEBUG
250 printf("c2="MDF" \n", c[2]);
251 #endif
252
253 /* now build pl[] and ql */
254 for ( i = 0; i < 3; i++ )
255 for ( j = 0; j < 3; j++ )
256 pl[i][j] = k*epsilon*(b[i]*b[j]+c[i]*c[j]);
257
258 ql[0] = A*(2*rho[0]+rho[1]+rho[2])*ONE_OVER_TWELVE;
259 ql[1] = A*(rho[0]+2*rho[1]+rho[2])*ONE_OVER_TWELVE;
260 ql[2] = A*(rho[0]+rho[1]+2*rho[2])*ONE_OVER_TWELVE;
261 #ifdef DEBUG
262 printf("local p and q ..\n");
263 for ( i = 0; i <3; i++) {
264 printf("| ");
265 for ( j = 0; j < 3; j++)
266 printf(" "MDF" ", pl[i][j]);
267 printf("|| "MDF" | \n", ql[i]);
268 }
269 #endif
270
271
272 /* better to use multiplication than division */
273 /* now augment the global data set */
274 for ( i = 0; i < 3; i++)
275 if ( v[i] < NUN ) { /* we index from 0, not from 1 */
276 q[v[i]] = q[v[i]] + ql[i];
277 for ( j = 0; j < 3; j++ ) {
278 if ( v[j] < NUN )
279 {
280 if ( v[j] <= v[i] )
281 P[v[i]][v[j]] = P[v[i]][v[j]] + pl[i][j];
282 /* we only build half since we have half only */
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" ", 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 /* a function to build the global matrix or the matrix */
311 /* equation for the system - Helmholtz Equation */
312 /* FIXME - need to solve only homogeneous case. dont need q */
313 static void
build_helmholtz(MY_DOUBLE ** P,MY_DOUBLE ** T,MY_DOUBLE * x,MY_DOUBLE * y,int NUN,int NE,int NN,MY_INT * renumber,MY_INT * re_renumber)314 build_helmholtz( MY_DOUBLE **P, MY_DOUBLE **T, MY_DOUBLE *x, MY_DOUBLE *y, int NUN, int NE, int NN, MY_INT *renumber, MY_INT *re_renumber)
315 {
316 /* infd - input file descriptor */
317 /* P[][] - global matrix to build NUNxNUN */
318 /* T[][] - global matrix to build NUNxNUN */
319 /* x[] - x coordinates NNx1 */
320 /* y[] - y coordinates NNx1 */
321 /* NUN - no of unknown potential nodes */
322 /* NE - no of elements */
323 /* NN - no of nodes */
324
325 /* now the running variables for each triangular element */
326 int v[3]; /* to keep vertex number of each node */
327 MY_DOUBLE pl[3][3]; /* local version of P[][] above */
328 MY_DOUBLE tl[3][3]; /* local version of T[][] above */
329 /*MY_DOUBLE ql[3];*/ /* local version of q[] above */
330 /*long double rho[3];*/ /* rho values */ /* conductivity */
331 triangle *t;
332
333 int i,j;
334
335 MY_DOUBLE
336 epsilon, mu, /* permittivity and permeability */
337 A, /* area of a triangulr element */
338 k, /* value of the determinanat */
339 yy0, yy1, yy2;
340
341 MY_DOUBLE b[3]; /*local arrays */
342 MY_DOUBLE c[3];
343
344
345
346 update_status_bar("Start building Eigenproblem.");
347 /* now read node data with unknown potentials */
348 for ( i = 0; i < NUN; i++) {
349 x[i]= (Mx(re_renumber[i],M))*g_xscale-g_xoff;
350 y[i]= (My(re_renumber[i],M))*g_yscale-g_yoff;
351 }
352
353 /* now read node data with known potentials , which comes last */
354 for ( i = NUN; i < NN; i++ ) {
355 x[i] = (Mx(re_renumber[i],M))*g_xscale-g_xoff;
356 y[i] = (My(re_renumber[i],M))*g_yscale-g_yoff;
357 }
358
359 #ifdef DEBUG
360 for ( i = 0; i< NN; i++)
361 printf("build: %d->%d "MDF", "MDF"\n", i,re_renumber[i], x[i], y[i]);
362 #endif
363
364 /* now read element data, whilst building the global matrices */
365 DAG_traverse_list_reset(&dt);
366 t=DAG_traverse_prune_list(&dt);
367 while(t) {
368 if ( t &&t->status==LIVE ) {
369
370 v[0] = renumber[t->p0];
371 v[1] = renumber[t->p1];
372 v[2] = renumber[t->p2]; /* ccw direction */
373 /* rho[0] = bound.poly_array[t->boundary].rho;
374 rho[1] = bound.poly_array[t->boundary].rho;
375 rho[2] = bound.poly_array[t->boundary].rho; */
376 /* for Homogeneous Helmholtz equation, calculate the ratio epsilon/mu
377 */
378 if(bound.poly_array[t->boundary].muexp) {
379 mu=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].muexp,t->p0)
380 +ex(bound.poly_array[t->boundary].muexp,t->p1)
381 +ex(bound.poly_array[t->boundary].muexp,t->p2));
382 } else {
383 mu=bound.poly_array[t->boundary].mu;
384 }
385 if(bound.poly_array[t->boundary].epsexp) {
386 epsilon=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].epsexp,t->p0)
387 +ex(bound.poly_array[t->boundary].epsexp,t->p1)
388 +ex(bound.poly_array[t->boundary].epsexp,t->p2));
389 } else {
390 epsilon=bound.poly_array[t->boundary].epsilon;
391 }
392 /* find the ratio */
393 epsilon /=mu;
394
395
396 /* now have read one triangle , add it to global */
397 /* first the coefficients */
398 /* some temporary values, neede several times */
399 yy0 = y[v[1]]-y[v[2]];
400 yy1 = y[v[2]]-y[v[0]];
401 yy2 = y[v[0]]-y[v[1]];
402
403 /* area of the triangle */
404 A = 0.5*( x[v[0]]*yy0 + x[v[1]]*yy1 + x[v[2]]*yy2 );
405 /*A = ABS(A); *//* we keep points (0,1,2) in ccw order, A > 0 */
406 /* another constant */
407 k = 0.25/(A*A);
408 /* now calculate b's and c's */
409 #ifdef DEBUG
410 printf("adding element with nodes %d %d and %d area "MDF"...?("MDF","MDF")("MDF","MDF")("MDF","MDF")\n\n", v[0], v[1], v[2],A,x[v[0]],y[v[0]],x[v[1]],y[v[1]],x[v[2]],y[v[2]]);
411 #endif
412 b[0]= yy0;
413 #ifdef DEBUG
414 printf("b0="MDF" ", b[0]);
415 #endif
416 b[1]= yy1;
417 #ifdef DEBUG
418 printf("b1="MDF" ", b[1]);
419 #endif
420 b[2]= yy2;
421 #ifdef DEBUG
422 printf("b2="MDF" \n", b[2]);
423 #endif
424 c[0] = ( x[v[2]]-x[v[1]] );
425 #ifdef DEBUG
426 printf("c0="MDF" ", c[0]);
427 #endif
428 c[1] = ( x[v[0]]-x[v[2]] );
429 #ifdef DEBUG
430 printf("c1="MDF" ", c[1]);
431 #endif
432 c[2] = ( x[v[1]]-x[v[0]] );
433 #ifdef DEBUG
434 printf("c2="MDF" \n", c[2]);
435 #endif
436
437 /* now build pl[] and ql */
438 for ( i = 0; i < 3; i++ )
439 for ( j = 0; j < 3; j++ )
440 pl[i][j] = k*epsilon*(b[i]*b[j]+c[i]*c[j]);
441
442 /* build tl */
443 tl[0][0] = tl[1][1] = tl[2][2] = A*ONE_OVER_SIX;
444 tl[0][1] = tl[0][2] = tl[1][0] = tl[1][2] = tl[2][0] = tl[2][1] = A*ONE_OVER_TWELVE;
445 /* ql[0] = A*(2*rho[0]+rho[1]+rho[2])/12.0;
446 ql[1] = A*(rho[0]+2*rho[1]+rho[2])/12.0;
447 ql[2] = A*(rho[0]+rho[1]+2*rho[2])/12.0; */
448 #ifdef DEBUG
449 printf("local p and q and t..\n");
450 for ( i = 0; i <3; i++) {
451 printf("| ");
452 for ( j = 0; j < 3; j++)
453 printf(" "MDF" ", pl[i][j]);
454 printf("||");
455 for ( j = 0; j < 3; j++)
456 printf(" "MDF" ", tl[i][j]);
457 printf("|\n");
458 }
459 #endif
460
461
462 /* better to use multiplication than division */
463 /* now augment the global data set */
464 for ( i = 0; i < 3; i++)
465 if ( v[i] < NUN ) { /* we index from 0, not from 1 */
466 /* q[v[i]] = q[v[i]] + ql[i]; */
467 for ( j = 0; j < 3; j++ ) {
468 if ( v[j] < NUN )
469 {
470 if ( v[j] <= v[i] ) { /* matricex P, T are symmetric */
471 P[v[i]][v[j]] = P[v[i]][v[j]] + pl[i][j];
472 T[v[i]][v[j]] = T[v[i]][v[j]] + tl[i][j];
473 }
474 /* we only build half since we have half only */
475 }/* else {
476 q[v[i]] = q[v[i]] - (pl[i][j])*phi[v[j]];
477 } */
478 }
479 } /*else printf("rejecting %d\n", v[i]);*/
480 #ifdef DEBUG
481 printf("building global matrix ..\n");
482 printf("global matrices P and T\n");
483 for ( i = 0; i < NUN; i++){
484 printf("| ");
485 for( j = 0; j <= i; j++)
486 printf(" "MDF" ", P[i][j]);
487 printf("|\n");
488 }
489 for ( i = 0; i < NUN; i++){
490 printf("| ");
491 for( j = 0; j <= i; j++)
492 printf(" "MDF" ", T[i][j]);
493 printf("|\n");
494 }
495
496 #endif
497 }
498 t=DAG_traverse_prune_list(&dt);
499 }
500
501 #ifdef DEBUG
502 printf("global matrices P and T\n");
503 for ( i = 0; i < NUN; i++){
504 printf("| ");
505 for( j = 0; j <= i; j++)
506 printf(" "MDF" ", P[i][j]);
507 printf("|\n");
508 }
509 for ( i = 0; i < NUN; i++){
510 printf("| ");
511 for( j = 0; j <= i; j++)
512 printf(" "MDF" ", T[i][j]);
513 printf("|\n");
514 }
515 #endif
516
517 update_status_bar("Done building Eigenproblem.");
518 }
519
520 /* a function to build the global matrix or the matrix */
521 /* equation for the system - Helmholtz Equation */
522 static void
build_helmholtz_sparse(SPMat * P,SPMat * T,MY_DOUBLE * x,MY_DOUBLE * y,int NUN,int NE,int NN,MY_INT * renumber,MY_INT * re_renumber)523 build_helmholtz_sparse( SPMat *P, SPMat *T, MY_DOUBLE *x, MY_DOUBLE *y, int NUN, int NE, int NN, MY_INT *renumber, MY_INT *re_renumber)
524 {
525 /* infd - input file descriptor */
526 /* P[][] - global matrix to build NUNxNUN */
527 /* T[][] - global matrix to build NUNxNUN */
528 /* x[] - x coordinates NNx1 */
529 /* y[] - y coordinates NNx1 */
530 /* NUN - no of unknown potential nodes */
531 /* NE - no of elements */
532 /* NN - no of nodes */
533
534 /* now the running variables for each triangular element */
535 int v[3]; /* to keep vertex number of each node */
536 MY_DOUBLE pl[3][3]; /* local version of P[][] above */
537 MY_DOUBLE tl[3][3]; /* local version of T[][] above */
538 /*MY_DOUBLE ql[3];*/ /* local version of q[] above */
539 /*long double rho[3];*/ /* rho values */ /* conductivity */
540 triangle *t;
541
542 int i,j;
543
544 MY_DOUBLE
545 epsilon, mu, /* permittivity and permeability */
546 A, /* area of a triangulr element */
547 k, /* value of the determinanat */
548 yy0, yy1, yy2;
549
550 MY_DOUBLE b[3]; /*local arrays */
551 MY_DOUBLE c[3];
552
553
554
555 update_status_bar("Start building Eigenproblem.");
556 /* now read node data with unknown potentials */
557 for ( i = 0; i < NUN; i++) {
558 x[i]= (Mx(re_renumber[i],M))*g_xscale-g_xoff;
559 y[i]= (My(re_renumber[i],M))*g_yscale-g_yoff;
560 }
561
562 /* now read node data with known potentials , which comes last */
563 for ( i = NUN; i < NN; i++ ) {
564 x[i] = (Mx(re_renumber[i],M))*g_xscale-g_xoff;
565 y[i] = (My(re_renumber[i],M))*g_yscale-g_yoff;
566 }
567
568 #ifdef DEBUG
569 for ( i = 0; i< NN; i++)
570 printf("build: %d->%d "MDF", "MDF"\n", i,re_renumber[i], x[i], y[i]);
571 #endif
572
573 /* now read element data, whilst building the global matrices */
574 DAG_traverse_list_reset(&dt);
575 t=DAG_traverse_prune_list(&dt);
576 while(t) {
577 if ( t &&t->status==LIVE ) {
578
579 v[0] = renumber[t->p0];
580 v[1] = renumber[t->p1];
581 v[2] = renumber[t->p2]; /* ccw direction */
582 /* rho[0] = bound.poly_array[t->boundary].rho;
583 rho[1] = bound.poly_array[t->boundary].rho;
584 rho[2] = bound.poly_array[t->boundary].rho; */
585 /* for Homogeneous Helmholtz equation, calculate the ratio epsilon/mu
586 */
587 if(bound.poly_array[t->boundary].muexp) {
588 mu=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].muexp,t->p0)
589 +ex(bound.poly_array[t->boundary].muexp,t->p1)
590 +ex(bound.poly_array[t->boundary].muexp,t->p2));
591 } else {
592 mu=bound.poly_array[t->boundary].mu;
593 }
594 if(bound.poly_array[t->boundary].epsexp) {
595 epsilon=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].epsexp,t->p0)
596 +ex(bound.poly_array[t->boundary].epsexp,t->p1)
597 +ex(bound.poly_array[t->boundary].epsexp,t->p2));
598 } else {
599 epsilon=bound.poly_array[t->boundary].epsilon;
600 }
601 /* find the ratio */
602 epsilon /=mu;
603
604
605 /* now have read one triangle , add it to global */
606 /* first the coefficients */
607 /* some temporary values, neede several times */
608 yy0 = y[v[1]]-y[v[2]];
609 yy1 = y[v[2]]-y[v[0]];
610 yy2 = y[v[0]]-y[v[1]];
611
612 /* area of the triangle */
613 A = 0.5*( x[v[0]]*yy0 + x[v[1]]*yy1 + x[v[2]]*yy2 );
614 /*A = ABS(A); *//* we keep points (0,1,2) in ccw order, A > 0 */
615 /* another constant */
616 k = 0.25/(A*A);
617 /* now calculate b's and c's */
618 #ifdef DEBUG
619 printf("adding element with nodes %d %d and %d area "MDF"...?("MDF","MDF")("MDF","MDF")("MDF","MDF")\n\n", v[0], v[1], v[2],A,x[v[0]],y[v[0]],x[v[1]],y[v[1]],x[v[2]],y[v[2]]);
620 #endif
621 b[0]= yy0;
622 #ifdef DEBUG
623 printf("b0="MDF" ", b[0]);
624 #endif
625 b[1]= yy1;
626 #ifdef DEBUG
627 printf("b1="MDF" ", b[1]);
628 #endif
629 b[2]= yy2;
630 #ifdef DEBUG
631 printf("b2="MDF" \n", b[2]);
632 #endif
633 c[0] = ( x[v[2]]-x[v[1]] );
634 #ifdef DEBUG
635 printf("c0="MDF" ", c[0]);
636 #endif
637 c[1] = ( x[v[0]]-x[v[2]] );
638 #ifdef DEBUG
639 printf("c1="MDF" ", c[1]);
640 #endif
641 c[2] = ( x[v[1]]-x[v[0]] );
642 #ifdef DEBUG
643 printf("c2="MDF" \n", c[2]);
644 #endif
645
646 /* now build pl[] and ql */
647 for ( i = 0; i < 3; i++ )
648 for ( j = 0; j < 3; j++ )
649 pl[i][j] = k*epsilon*(b[i]*b[j]+c[i]*c[j]);
650
651 /* build tl */
652 tl[0][0] = tl[1][1] = tl[2][2] = A*ONE_OVER_SIX;
653 tl[0][1] = tl[0][2] = tl[1][0] = tl[1][2] = tl[2][0] = tl[2][1] = A*ONE_OVER_TWELVE;
654 /* ql[0] = A*(2*rho[0]+rho[1]+rho[2])/12.0;
655 ql[1] = A*(rho[0]+2*rho[1]+rho[2])/12.0;
656 ql[2] = A*(rho[0]+rho[1]+2*rho[2])/12.0; */
657 #ifdef DEBUG
658 printf("local p and q and t..\n");
659 for ( i = 0; i <3; i++) {
660 printf("| ");
661 for ( j = 0; j < 3; j++)
662 printf(" "MDF" ", pl[i][j]);
663 printf("||");
664 for ( j = 0; j < 3; j++)
665 printf(" "MDF" ", tl[i][j]);
666 printf("|\n");
667 }
668 #endif
669
670
671 /* better to use multiplication than division */
672 /* now augment the global data set */
673 for ( i = 0; i < 3; i++)
674 if ( v[i] < NUN ) { /* we index from 0, not from 1 */
675 /* q[v[i]] = q[v[i]] + ql[i]; */
676 for ( j = 0; j < 3; j++ ) {
677 if ( v[j] < NUN )
678 {
679 if ( v[j] <= v[i] ) { /* matricex P, T are symmetric */
680 /*P[v[i]][v[j]] = P[v[i]][v[j]] + pl[i][j]; */
681 SPM_add(P,v[i],v[j],pl[i][j]);
682 /*T[v[i]][v[j]] = T[v[i]][v[j]] + tl[i][j]; */
683 SPM_add(T,v[i],v[j],tl[i][j]);
684 }
685 /* we only build half since we have half only */
686 }
687 }
688 } /*else printf("rejecting %d\n", v[i]);*/
689 #ifdef DEBUG
690 printf("building global matrix ..\n");
691 printf("global matrices P and T\n");
692 for ( i = 0; i < NUN; i++){
693 printf("| ");
694 for( j = 0; j <= i; j++)
695 printf(" "MDF" ", SPM_e(P,i,j));
696 printf("|\n");
697 }
698 for ( i = 0; i < NUN; i++){
699 printf("| ");
700 for( j = 0; j <= i; j++)
701 printf(" "MDF" ", SPM_e(T,i,j));
702 printf("|\n");
703 }
704
705 #endif
706 }
707 t=DAG_traverse_prune_list(&dt);
708 }
709
710 #ifdef DEBUG
711 printf("global matrices P and T\n");
712 for ( i = 0; i < NUN; i++){
713 printf("| ");
714 for( j = 0; j <= i; j++)
715 printf(" "MDF" ", SPM_e(P,i,j));
716 printf("|\n");
717 }
718 for ( i = 0; i < NUN; i++){
719 printf("| ");
720 for( j = 0; j <= i; j++)
721 printf(" "MDF" ", SPM_e(T,i,j));
722 printf("|\n");
723 }
724 #endif
725
726 update_status_bar("Done building Eigenproblem.");
727 }
728
729 #if defined (USE_LAPACK) || defined (USE_MKL)
730 /* nothing */
731 #else
732 static void
solve_helmholtz(MY_DOUBLE ** P,MY_DOUBLE ** T,MY_DOUBLE ** x,MY_DOUBLE * ev,MY_INT NUN,MY_INT n_to_find)733 solve_helmholtz(MY_DOUBLE **P,MY_DOUBLE **T, MY_DOUBLE **x, MY_DOUBLE *ev, MY_INT NUN, MY_INT n_to_find)
734 {
735 /* solves the symmetrix, generalized eigenvalue problem
736 * (P-lambda. T)x= 0
737 * P, T = NUN by NUN matrices, only lower triangle
738 * x = eigenvectors to be found size (m x n_to_find)
739 * ev = eigenvalues to be found, n_to_find by 1 array
740 */
741 MY_INT i,j,k;
742 MY_DOUBLE **L,**Pt;
743 MY_DOUBLE temp;
744 /* first find the cholesky factorization of T */
745 /* T=L L' */
746 if((L= (MY_DOUBLE**) calloc(NUN, sizeof(MY_DOUBLE*)))==0){
747 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
748 exit(1);
749 }
750 for (i = 0; i <NUN; ++i) {
751 if((L[i] = ( MY_DOUBLE * ) calloc( (size_t)i+1, sizeof(MY_DOUBLE)))==0){
752 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
753 exit(1);
754 }
755 memcpy((void*)L[i],(void*)T[i],(size_t)(i+1)*sizeof(MY_DOUBLE));
756 }
757
758 for ( j = 0; j < NUN; j++) {
759 /* first diagonal element */
760 temp = 0;
761 for ( k = 0; k < j; k++)
762 temp += L[j][k]*L[j][k];
763 L[j][j] = sqrtf((float)( L[j][j] - temp) );
764 #ifdef DEBUG
765 printf("solve_helmholtz: col %d diag "MDF"\n",j,L[j][j]);
766 #endif
767 /* check for singularities */
768 if ( L[j][j] <= 0.0f ) {
769 L[j][j] = TOL; /* a small value */
770 fprintf(stderr,"warning: solve_cholevsky: singular matrix\n");
771 }
772 /* now handle the i th column */
773 for( i = j+1; i < NUN; i++) {
774 temp = 0;
775 for ( k = 0; k < j; k++)
776 temp += L[i][k]*L[j][k];
777 L[i][j] = ( L[i][j] - temp )/L[j][j];
778 #ifdef DEBUG
779 printf("solve_helm: row %d elem "MDF"\n",i,L[i][j]);
780 #endif
781 }
782 }
783
784 #ifdef DEBUG
785 for ( i =0; i < NUN ; i++ ) {
786 for ( j = 0 ; j <= i; j++ )
787 printf(""MDF" ",L[i][j]);
788 printf("\n");
789 }
790 #endif
791
792 /* now we solve the system
793 * (L^{-1} P L^{-T}-lambda I)y =0
794 * and Ly =x
795 * to get x
796 */
797 /* first: create a full matrix for P */
798 if((Pt= (MY_DOUBLE**) calloc(NUN, sizeof(MY_DOUBLE*)))==0){
799 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
800 exit(1);
801 }
802 for (i = 0; i <NUN; ++i) {
803 if((Pt[i] = ( MY_DOUBLE * ) calloc( (size_t)NUN, sizeof(MY_DOUBLE)))==0){
804 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
805 exit(1);
806 }
807 memcpy((void*)Pt[i],(void*)P[i],(size_t)(i+1)*sizeof(MY_DOUBLE));
808 /* fill up the rest of Pt */
809 for (j=0;j<i;j++) {
810 Pt[j][i]=P[i][j];
811 }
812 }
813
814
815 /* calculation of L^{-1}P */
816 /* column by column */
817 /* forward elimination */
818 for (k=0;k<NUN;k++) {
819 for (i=0; i<NUN; i++) {
820 temp=0;
821 for (j=0; j<i;j++) {
822 temp+=L[i][j]*Pt[j][k];
823 }
824 Pt[i][k]=(Pt[i][k]-temp)/L[i][i];
825 }
826 }
827
828 /* calculation of P L^{-T} */
829 /* row by row */
830 for (k=0;k<NUN;k++) {
831 for (i=0;i<NUN;i++) {
832 temp=0;
833 for (j=0; j<i;j++) {
834 temp+=L[i][j]*Pt[k][j];
835 }
836 Pt[k][i]=(Pt[k][i]-temp)/L[i][i];
837 }
838 }
839 #ifdef DEBUG
840 for ( i =0; i < NUN ; i++ ) {
841 for ( j = 0 ; j <NUN; j++ )
842 printf(""MDF" ",Pt[i][j]);
843 printf("\n");
844 }
845 #endif
846
847 /* find the eigenvector */
848 /*find_eigenvectors(Pt, x, ev, NUN-2, NUN); */
849 /*x - m x n_to_find */
850 find_eigenvectors(Pt, x, ev, NUN, n_to_find);
851
852 /* now we have L^{T} x = eigenvector
853 * hence solve to get true eigenvector */
854 for (k=0;k<n_to_find;k++) {
855 for (i=NUN-1;i>=0;i--) {
856 temp=0;
857 for (j=NUN-1; j>i;j--) {
858 temp+=L[j][i]*x[i][k];
859 }
860 x[i][k]=(x[i][k]-temp)/L[i][i];
861
862 #ifdef DEBUG
863 printf("x[%d][%d]="MDF" ",i,k,x[i][k]);
864 #endif
865 }
866 }
867 #ifdef DEBUG
868 printf("\n");
869 #endif
870 for (i=0; i<NUN; i++) {
871 free(L[i]);
872 free(Pt[i]);
873 }
874 free(L);
875 free(Pt);
876 }
877
878 static void
solve_helmholtz_sparse(SPMat * P,SPMat * T,MY_DOUBLE ** x,MY_DOUBLE * ev,MY_INT NUN,MY_INT n_to_find)879 solve_helmholtz_sparse(SPMat *P,SPMat *T, MY_DOUBLE **x, MY_DOUBLE *ev, MY_INT NUN, MY_INT n_to_find)
880 {
881 /* solves the symmetrix, generalized eigenvalue problem
882 * (P-lambda. T)x= 0
883 * P, T = NUN by NUN matrices, only lower triangle
884 * x = eigenvectors to be found size (m x n_to_find)
885 * ev = eigenvalues to be found, n_to_find by 1 array
886 */
887 MY_INT i,j,k;
888 MY_DOUBLE **L,**Pt;
889 MY_DOUBLE temp;
890 /* first find the cholesky factorization of T */
891 /* T=L L' */
892 if((L= (MY_DOUBLE**) calloc(NUN, sizeof(MY_DOUBLE*)))==0){
893 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
894 exit(1);
895 }
896 for (i = 0; i <NUN; ++i) {
897 if((L[i] = ( MY_DOUBLE * ) calloc( (size_t)i+1, sizeof(MY_DOUBLE)))==0){
898 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
899 exit(1);
900 }
901 //memcpy((void*)L[i],(void*)T[i],(size_t)(i+1)*sizeof(MY_DOUBLE));
902 for (j=0; j<=i; j++)
903 L[i][j]=SPM_e(T,i,j);
904 }
905
906 for ( j = 0; j < NUN; j++) {
907 /* first diagonal element */
908 temp = 0;
909 for ( k = 0; k < j; k++)
910 temp += L[j][k]*L[j][k];
911 L[j][j] = sqrtf((float)( L[j][j] - temp) );
912 #ifdef DEBUG
913 printf("solve_helmholtz: col %d diag "MDF"\n",j,L[j][j]);
914 #endif
915 /* check for singularities */
916 if ( L[j][j] <= 0.0f ) {
917 L[j][j] = TOL; /* a small value */
918 fprintf(stderr,"warning: solve_cholevsky: singular matrix\n");
919 }
920 /* now handle the i th column */
921 for( i = j+1; i < NUN; i++) {
922 temp = 0;
923 for ( k = 0; k < j; k++)
924 temp += L[i][k]*L[j][k];
925 L[i][j] = ( L[i][j] - temp )/L[j][j];
926 #ifdef DEBUG
927 printf("solve_helm: row %d elem "MDF"\n",i,L[i][j]);
928 #endif
929 }
930 }
931
932 #ifdef DEBUG
933 for ( i =0; i < NUN ; i++ ) {
934 for ( j = 0 ; j <= i; j++ )
935 printf(""MDF" ",L[i][j]);
936 printf("\n");
937 }
938 #endif
939
940 /* now we solve the system
941 * (L^{-1} P L^{-T}-lambda I)y =0
942 * and Ly =x
943 * to get x
944 */
945 /* first: create a full matrix for P */
946 if((Pt= (MY_DOUBLE**) calloc(NUN, sizeof(MY_DOUBLE*)))==0){
947 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
948 exit(1);
949 }
950 for (i = 0; i <NUN; ++i) {
951 if((Pt[i] = ( MY_DOUBLE * ) calloc( (size_t)NUN, sizeof(MY_DOUBLE)))==0){
952 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
953 exit(1);
954 }
955 //memcpy((void*)Pt[i],(void*)P[i],(size_t)(i+1)*sizeof(MY_DOUBLE));
956 for(j=0;j<=i;j++)
957 Pt[j][i]=Pt[i][j]=SPM_e(P,i,j);
958 }
959
960
961 /* calculation of L^{-1}P */
962 /* column by column */
963 /* forward elimination */
964 for (k=0;k<NUN;k++) {
965 for (i=0; i<NUN; i++) {
966 temp=0;
967 for (j=0; j<i;j++) {
968 temp+=L[i][j]*Pt[j][k];
969 }
970 Pt[i][k]=(Pt[i][k]-temp)/L[i][i];
971 }
972 }
973
974 /* calculation of P L^{-T} */
975 /* row by row */
976 for (k=0;k<NUN;k++) {
977 for (i=0;i<NUN;i++) {
978 temp=0;
979 for (j=0; j<i;j++) {
980 temp+=L[i][j]*Pt[k][j];
981 }
982 Pt[k][i]=(Pt[k][i]-temp)/L[i][i];
983 }
984 }
985 #ifdef DEBUG
986 for ( i =0; i < NUN ; i++ ) {
987 for ( j = 0 ; j <NUN; j++ )
988 printf(""MDF" ",Pt[i][j]);
989 printf("\n");
990 }
991 #endif
992
993 /* find the eigenvector */
994 /*find_eigenvectors(Pt, x, ev, NUN-2, NUN); */
995 /*x - m x n_to_find */
996 find_eigenvectors(Pt, x, ev, NUN, n_to_find);
997
998 /* now we have L^{T} x = eigenvector
999 * hence solve to get true eigenvector */
1000 for (k=0;k<n_to_find;k++) {
1001 for (i=NUN-1;i>=0;i--) {
1002 temp=0;
1003 for (j=NUN-1; j>i;j--) {
1004 temp+=L[j][i]*x[i][k];
1005 }
1006 x[i][k]=(x[i][k]-temp)/L[i][i];
1007
1008 #ifdef DEBUG
1009 printf("x[%d][%d]="MDF" ",i,k,x[i][k]);
1010 #endif
1011 }
1012 }
1013 #ifdef DEBUG
1014 printf("\n");
1015 #endif
1016 for (i=0; i<NUN; i++) {
1017 free(L[i]);
1018 free(Pt[i]);
1019 }
1020 free(L);
1021 free(Pt);
1022 }
1023
1024 #endif /* !USE_LAPACK || USE_MKL */
1025
1026 /* after genarating the mesh, need to re-number points */
1027 /* so that unknown points are first */
1028 /* then we solve the problem -poisson equation*/
1029 int
re_number_and_solve_poisson(void)1030 re_number_and_solve_poisson(void)
1031 {
1032 unsigned MY_INT i,j;
1033 MY_INT * renumber; /* re-numbering array */
1034 MY_INT * re_renumber; /* re-renumbering array */
1035 /* used to map potentials back to the original points */
1036
1037
1038 unsigned MY_INT unknowns=0; /* unknown points */
1039 unsigned MY_INT total = 0; /* total points */
1040 MY_INT triangles = 0; /* known points */
1041
1042
1043 triangle *t;
1044
1045 /* first the global node data, all total size */
1046 MY_DOUBLE *x ; /* x coordinates */
1047 MY_DOUBLE *y ; /* y coordinates */
1048 MY_DOUBLE *phi ; /* potentials at each node */
1049
1050 /* the global variables */
1051 MY_DOUBLE **P ; /* size unknownsxunknowns note we do not need a NN x NN matrix, only the unknowns
1052 */
1053 MY_DOUBLE *q ; /* unknowns sizesame here */
1054
1055
1056 /* allocate memory for arrays */
1057 if ((renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT )))==0){
1058 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1059 exit(1);
1060 }
1061
1062 if((re_renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
1063 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1064 exit(1);
1065 }
1066
1067 /* need to solve for only points referred by triangles */
1068 /* SO -- */
1069 DAG_traverse_list_reset(&dt);
1070 t=DAG_traverse_prune_list(&dt);
1071 while ( t ) {
1072 if ( t && t->status==LIVE ) {
1073 /* mark the points in this triangle */
1074 /* use re_renumber[] array for this */
1075 re_renumber[t->p0] = 1;
1076 re_renumber[t->p1] = 1;
1077 re_renumber[t->p2] = 1;
1078
1079 triangles++;
1080 }
1081 t=DAG_traverse_prune_list(&dt);
1082 }
1083
1084 /* now find unknowns, or INSIDE points */
1085 j = 0;
1086 for ( i = 0; i < M.count; i++ )
1087 if ( Mval(i,M) == VAR && re_renumber[i] == 1) {
1088 renumber[i] = j++;
1089 #ifdef DEBUG
1090 printf("renumber_and_solve: renumbering variable point %d as %d\n",i,renumber[i]);
1091 #endif
1092 unknowns++;
1093 }
1094 /* now renumber the known points */
1095 total = unknowns;
1096 for ( i = 0; i < M.count; i++ )
1097 if ( Mval(i,M) == FX && re_renumber[i]==1 ) {
1098 renumber[i] = j++;
1099 #ifdef DEBUG
1100 printf("renumber_and_solve: renumbering fixed point %d as %d\n",i,renumber[i]);
1101 #endif
1102 total++;
1103 }
1104 /* finally give numbers to points 0,1,2 */
1105 /* we will not use these points */
1106 renumber[0]=j++;
1107 renumber[1]=j++;
1108 renumber[2]=j++;
1109
1110 /* now construct re-renmbering array */
1111 i = 0;
1112 /* do an insertion sort */
1113 while ( i < M.count ) {
1114 for ( j = 0; j < M.count; j++ ) {
1115 if ( i==(unsigned)renumber[j] )
1116 break;
1117 }
1118 re_renumber[i] = j;
1119 i++;
1120 }
1121
1122 #ifdef DEBUG
1123 printf("build: renumbering %d points with %d unknowns out of %d total\n",total,unknowns,M.count);
1124 #endif
1125 /*renumber_points(renumber, re_renumber, total, unknowns); */
1126
1127 #ifdef DEBUG
1128 printf("i, renum, rerenum\n");
1129 for ( i = 0; i < total; i++ )
1130 printf("%d %d %d\n",i,renumber[i],re_renumber[i]);
1131 for ( i = 0; i < total; i++ )
1132 printf("%d "MDF" "MDF" "MDF"\n",i,Mx(i,M),My(i,M),Mz(i,M));
1133 #endif
1134
1135
1136 /* we free these two arrays only after solving the problem */
1137 /* now solve the problem */
1138 /* first allocate memory for the arrays */
1139 if((x=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
1140 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1141 exit(1);
1142 }
1143
1144 if((y=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
1145 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1146 exit(1);
1147 }
1148
1149 if((phi=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
1150 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1151 exit(1);
1152 }
1153
1154 if((q=(MY_DOUBLE*)calloc(unknowns, sizeof(MY_DOUBLE)))==0){
1155 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1156 exit(1);
1157 }
1158
1159
1160 /* now the matrix */
1161 /* STEP 1: SET UP THE ROWS. */
1162 if((P= (MY_DOUBLE**) calloc(unknowns , sizeof(MY_DOUBLE*)))==0){
1163 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1164 exit(1);
1165 }
1166 /* STEP 2: SET UP THE COLUMNS. */
1167 for (i = 0; i < unknowns; ++i)
1168 /* variable length columns */
1169 if((P[i] = ( MY_DOUBLE * ) calloc( (size_t)i+1, sizeof(MY_DOUBLE)))==0){
1170 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1171 exit(1);
1172 }
1173
1174 /* Poisson equation */
1175 /* now the memory allocation is over */
1176 build( P, q, x, y, phi, unknowns, triangles, total, renumber, re_renumber);
1177 /* if user has chosen so, solve by conj. grad method*/
1178 /* by default solve by cholevsky decomposition */
1179 solve_cholevsky(P,phi, q, unknowns);
1180
1181 /* after solution, array phi[] will contain the new potentials */
1182 /* now put them back to the point data */
1183 /* we also note the maximum and minimum values for contour plotting */
1184 g_maxpot = -1000;
1185 g_minpot = 1000;
1186 for ( i = 0; i < total; i++ ) {
1187 Mz(re_renumber[i],M) = phi[i];
1188 if ( phi[i] > g_maxpot )
1189 g_maxpot = phi[i];
1190 if ( phi[i] < g_minpot )
1191 g_minpot = phi[i];
1192 }
1193
1194
1195
1196 /* now we can free the renumber[] array */
1197 free(renumber);
1198 /* now free this array too */
1199 free(re_renumber);
1200
1201 /* free everything else */
1202
1203 /* first the arrays */
1204 free(x);
1205 free(y);
1206 free(phi);
1207 free(q);
1208
1209 /* next the matrix */
1210 for ( i = 0; i < unknowns ; i++)
1211 free(P[i]);
1212 free(P);
1213
1214 /* now printout a summery */
1215 #ifdef DEBUG
1216 for ( i = 0; i < total; i++ )
1217 printf("renumber_and_solve: %d "MDF" "MDF" "MDF"\n",i,Mx(i,M),My(i,M),Mz(i,M));
1218 #endif
1219
1220
1221 return(0);
1222 }
1223
1224 /* solve helmholtz equation */
1225 int
re_number_and_solve_helmholtz(void)1226 re_number_and_solve_helmholtz(void)
1227 {
1228 unsigned MY_INT i,j;
1229 MY_INT k;
1230 MY_INT * renumber; /* re-numbering array */
1231 MY_INT * re_renumber; /* re-renumbering array */
1232 /* used to map potentials back to the original points */
1233
1234
1235 unsigned MY_INT unknowns=0; /* unknown points */
1236 unsigned MY_INT total = 0; /* total points */
1237 MY_INT triangles = 0; /* known points */
1238
1239
1240 triangle *t;
1241
1242 /* first the global node data, all total size */
1243 MY_DOUBLE *x ; /* x coordinates */
1244 MY_DOUBLE *y ; /* y coordinates */
1245 MY_DOUBLE **v; /* eigenvector matrix */
1246
1247 MY_DOUBLE *max_value_array;
1248
1249 /* the global variables */
1250 MY_DOUBLE **P ; /* size unknowns x unknowns note we do not need a NN x NN matrix, only the unknowns
1251 */
1252 MY_DOUBLE **T=0; /* used in Helmholtz equation */
1253
1254
1255 /* allocate memory for arrays */
1256 if ((renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT )))==0){
1257 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1258 exit(1);
1259 }
1260
1261 if((re_renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
1262 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1263 exit(1);
1264 }
1265
1266 /* need to solve for only points referred by triangles */
1267 /* SO -- */
1268 DAG_traverse_list_reset(&dt);
1269 t=DAG_traverse_prune_list(&dt);
1270 while ( t ) {
1271 if ( t && t->status==LIVE ) {
1272 /* mark the points in this triangle */
1273 /* use re_renumber[] array for this */
1274 re_renumber[t->p0] = 1;
1275 re_renumber[t->p1] = 1;
1276 re_renumber[t->p2] = 1;
1277
1278 triangles++;
1279 }
1280 t=DAG_traverse_prune_list(&dt);
1281 }
1282
1283 /* now find unknowns, or INSIDE points */
1284 j = 0;
1285 for ( i = 0; i < M.count; i++ )
1286 if ( Mval(i,M) == VAR && re_renumber[i] == 1) {
1287 renumber[i] = j++;
1288 #ifdef DEBUG
1289 printf("renumber_and_solve: renumbering variable point %d as %d\n",i,renumber[i]);
1290 #endif
1291 unknowns++;
1292 }
1293 /* now renumber the known points */
1294 total = unknowns;
1295 for ( i = 0; i < M.count; i++ )
1296 if ( Mval(i,M) == FX && re_renumber[i]==1 ) {
1297 renumber[i] = j++;
1298 #ifdef DEBUG
1299 printf("renumber_and_solve: renumbering fixed point %d as %d\n",i,renumber[i]);
1300 #endif
1301 total++;
1302 }
1303 /* finally give numbers to points 0,1,2 */
1304 /* we will not use these points */
1305 renumber[0]=j++;
1306 renumber[1]=j++;
1307 renumber[2]=j++;
1308
1309 /* now construct re-renmbering array */
1310 i = 0;
1311 /* do an insertion sort */
1312 while ( i < M.count ) {
1313 for ( j = 0; j < M.count; j++ ) {
1314 if ( renumber[j] == (int)i )
1315 break;
1316 }
1317 re_renumber[i] = j;
1318 i++;
1319 }
1320
1321 #ifdef DEBUG
1322 printf("build: renumbering %d points with %d unknowns out of %d total\n",total,unknowns,M.count);
1323 #endif
1324 /*renumber_points(renumber, re_renumber, total, unknowns); */
1325
1326 #ifdef DEBUG
1327 printf("i, renum, rerenum\n");
1328 for ( i = 0; i < total; i++ )
1329 printf("%d %d %d\n",i,renumber[i],re_renumber[i]);
1330 for ( i = 0; i < total; i++ )
1331 printf("%d "MDF" "MDF" "MDF"\n",i,Mx(i,M),My(i,M),Mz(i,M));
1332 #endif
1333
1334
1335 /* we free these two arrays only after solving the problem */
1336 /* now solve the problem */
1337 /* first allocate memory for the arrays */
1338 if((x=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
1339 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1340 exit(1);
1341 }
1342
1343 if((y=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
1344 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1345 exit(1);
1346 }
1347
1348 /* structure to store all eigenvectors - matrix */
1349 if((v=(MY_DOUBLE **)calloc(unknowns, sizeof(MY_DOUBLE*)))==0){
1350 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1351 exit(1);
1352 }
1353 for (i = 0; i <unknowns ; i++)
1354 if((v[i] = ( MY_DOUBLE * ) calloc( (size_t) degree_of_freedom, sizeof(MY_DOUBLE)))==0){
1355 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1356 exit(1);
1357 }
1358
1359
1360 /* now the matrix */
1361 /* STEP 1: SET UP THE ROWS. */
1362 if((P= (MY_DOUBLE**) calloc(unknowns , sizeof(MY_DOUBLE*)))==0){
1363 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1364 exit(1);
1365 }
1366 /* STEP 2: SET UP THE COLUMNS. */
1367 for (i = 0; i < unknowns; ++i)
1368 /* variable length columns */
1369 if((P[i] = ( MY_DOUBLE * ) calloc( (size_t)i+1, sizeof(MY_DOUBLE)))==0){
1370 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1371 exit(1);
1372 }
1373
1374 if((T= (MY_DOUBLE**) calloc(unknowns , sizeof(MY_DOUBLE*)))==0){
1375 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1376 exit(1);
1377 }
1378 for (i = 0; i < unknowns; ++i)
1379 if((T[i] = ( MY_DOUBLE * ) calloc( (size_t)i+1, sizeof(MY_DOUBLE)))==0){
1380 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1381 exit(1);
1382 }
1383 #ifdef DEBUG
1384 printf("solving problem with %d unknowns, %d triangles\n",unknowns,triangles);
1385 #endif
1386
1387 /*build_helmholtz( MY_DOUBLE **P, MY_DOUBLE **T, MY_DOUBLE *x, MY_DOUBLE *y, int NUN, int NE, int NN) */
1388 build_helmholtz( P, T, x, y, unknowns, triangles, total, renumber, re_renumber);
1389
1390 /*solve_helmholtz(MY_DOUBLE **P,MY_DOUBLE **T, MY_DOUBLE **x, MY_INT NUN, MY_INT n_to_find) */
1391 #if defined (USE_LAPACK) || defined (USE_MKL)
1392 solve_helmholtz_lapack(P,T,v,eig_array,unknowns,degree_of_freedom);
1393 #else
1394 solve_helmholtz(P,T,v,eig_array,unknowns,degree_of_freedom);
1395 #endif /* ! USE_LAPACK || USE_MKL */
1396
1397 /* all eigenvectors are normalized, hence normlization has to be done individually */
1398 /* we copy back elements of eigenvectors to z arrays of each point */
1399 /* in order to do that, we need to change size of those z arrays */
1400 /* find ranges for potentials */
1401 g_maxpot = -1.0;
1402 g_minpot = 1.0;
1403 /* normalize each eigenvector such that the absolute max
1404 value is 1 */
1405 if((max_value_array=( MY_DOUBLE * ) calloc( (size_t) degree_of_freedom, sizeof(MY_DOUBLE)))==0){
1406 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1407 exit(1);
1408 }
1409
1410 for(i=0; i<unknowns;i++) {
1411 for (k=0;k< degree_of_freedom;k++) {
1412 if ( ABS(v[i][k])> max_value_array[k]) max_value_array[k]=ABS(v[i][k]);
1413 }
1414 }
1415 /* find reciprocals of max values for division */
1416 for(k=0;k<degree_of_freedom;k++) {
1417 if(!IS_ZERO(max_value_array[k]))
1418 max_value_array[k]=1/max_value_array[k];
1419 }
1420
1421 /* after finding the max values, divide each column
1422 of v by its max abs value */
1423 for(i=0; i<unknowns;i++) {
1424 for (k=0;k< degree_of_freedom;k++) {
1425 v[i][k] *=max_value_array[k];
1426 }
1427 }
1428
1429 /* reallocate memory for z array and copy all eigenvectors */
1430 for ( i = 0; i < unknowns; i++ ) {
1431 if (((M.Narray[re_renumber[i]])->z=(MY_DOUBLE *)realloc((void*)(M.Narray[re_renumber[i]])->z,(size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0) {
1432 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1433 exit(2);
1434 }
1435 memcpy((void*)(M.Narray[re_renumber[i]])->z,(void*)v[i],(size_t)degree_of_freedom*sizeof(MY_DOUBLE));
1436
1437 }
1438 /* now consider fixed points: fill their z arrays with zeros */
1439 for ( i = unknowns; i < total; i++ ) {
1440 if (((M.Narray[re_renumber[i]])->z=(MY_DOUBLE *)realloc((void*)(M.Narray[re_renumber[i]])->z,(size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0) {
1441 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1442 exit(2);
1443 }
1444 /* fill zeros */
1445 memset((void*)(M.Narray[re_renumber[i]])->z,0,(size_t)degree_of_freedom*sizeof(MY_DOUBLE));
1446 }
1447
1448
1449
1450 /* now we can free the renumber[] array */
1451 free(renumber);
1452 /* now free this array too */
1453 free(re_renumber);
1454
1455 /* free everything else */
1456
1457 /* first the arrays */
1458 free(x);
1459 free(y);
1460
1461 /* next the matrix */
1462 for ( i = 0; i < unknowns ; i++) {
1463 free(P[i]);
1464 free(v[i]);
1465 free(T[i]);
1466 }
1467 free(P);
1468 free(v);
1469 free(max_value_array);
1470 free(T);
1471
1472 /* now printout a summery */
1473 #ifdef DEBUG
1474 for ( i = 0; i < total; i++ )
1475 printf("renumber_and_solve: %d "MDF" "MDF" "MDF"\n",i,Mx(i,M),My(i,M),Mz(i,M));
1476 #endif
1477
1478 return(0);
1479 }
1480
1481 /* solve helmholtz equation using sparce matrices */
1482 int
re_number_and_solve_helmholtz_sparse(Mesh Ms,DAG_tree Dt)1483 re_number_and_solve_helmholtz_sparse(Mesh Ms, DAG_tree Dt)
1484 {
1485 unsigned MY_INT i,j;
1486 MY_INT k;
1487 MY_INT * renumber; /* re-numbering array */
1488 MY_INT * re_renumber; /* re-renumbering array */
1489 /* used to map potentials back to the original points */
1490
1491
1492 unsigned MY_INT unknowns=0; /* unknown points */
1493 unsigned MY_INT total = 0; /* total points */
1494 MY_INT triangles = 0; /* known points */
1495
1496
1497 triangle *t;
1498
1499 /* first the global node data, all total size */
1500 MY_DOUBLE *x ; /* x coordinates */
1501 MY_DOUBLE *y ; /* y coordinates */
1502 MY_DOUBLE **v; /* eigenvector matrix */
1503
1504 MY_DOUBLE *max_value_array;
1505 /* the global variables */
1506 SPMat P ; /* size unknowns x unknowns note we do not need a NN x NN matrix, only the unknowns
1507 */
1508 SPMat T; /* used in Helmholtz equation */
1509
1510
1511 /* allocate memory for arrays */
1512 if ((renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT )))==0){
1513 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1514 exit(1);
1515 }
1516
1517 if((re_renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
1518 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1519 exit(1);
1520 }
1521
1522 /* need to solve for only points referred by triangles */
1523 /* SO -- */
1524 DAG_traverse_list_reset(&Dt);
1525 t=DAG_traverse_prune_list(&Dt);
1526 while ( t ) {
1527 if ( t && t->status==LIVE ) {
1528 /* mark the points in this triangle */
1529 /* use re_renumber[] array for this */
1530 re_renumber[t->p0] = 1;
1531 re_renumber[t->p1] = 1;
1532 re_renumber[t->p2] = 1;
1533
1534 triangles++;
1535 }
1536 t=DAG_traverse_prune_list(&Dt);
1537 }
1538
1539 /* now find unknowns, or INSIDE points */
1540 j = 0;
1541 for ( i = 0; i < Ms.count; i++ )
1542 if ( Mval(i,Ms) == VAR && re_renumber[i] == 1) {
1543 renumber[i] = j++;
1544 #ifdef DEBUG
1545 printf("renumber_and_solve: renumbering variable point %d as %d\n",i,renumber[i]);
1546 #endif
1547 unknowns++;
1548 }
1549 /* now renumber the known points */
1550 total = unknowns;
1551 for ( i = 0; i < Ms.count; i++ )
1552 if ( Mval(i,Ms) == FX && re_renumber[i]==1 ) {
1553 renumber[i] = j++;
1554 #ifdef DEBUG
1555 printf("renumber_and_solve: renumbering fixed point %d as %d\n",i,renumber[i]);
1556 #endif
1557 total++;
1558 }
1559 /* finally give numbers to points 0,1,2 */
1560 /* we will not use these points */
1561 renumber[0]=j++;
1562 renumber[1]=j++;
1563 renumber[2]=j++;
1564
1565 /* now construct re-renmbering array */
1566 i = 0;
1567 /* do an insertion sort */
1568 while ( i < Ms.count ) {
1569 for ( j = 0; j < Ms.count; j++ ) {
1570 if ( renumber[j] == (int)i )
1571 break;
1572 }
1573 re_renumber[i] = j;
1574 i++;
1575 }
1576
1577 #ifdef DEBUG
1578 printf("build: renumbering %d points with %d unknowns out of %d total\n",total,unknowns,Ms.count);
1579 #endif
1580 /*renumber_points(renumber, re_renumber, total, unknowns); */
1581
1582 #ifdef DEBUG
1583 printf("i, renum, rerenum\n");
1584 for ( i = 0; i < total; i++ )
1585 printf("%d %d %d\n",i,renumber[i],re_renumber[i]);
1586 for ( i = 0; i < total; i++ )
1587 printf("%d "MDF" "MDF" "MDF"\n",i,Mx(i,Ms),My(i,Ms),Mz(i,Ms));
1588 #endif
1589
1590
1591 /* we free these two arrays only after solving the problem */
1592 /* now solve the problem */
1593 /* first allocate memory for the arrays */
1594 if((x=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
1595 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1596 exit(1);
1597 }
1598
1599 if((y=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
1600 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1601 exit(1);
1602 }
1603
1604 /* structure to store all eigenvectors - matrix */
1605 if((v=(MY_DOUBLE **)calloc(unknowns, sizeof(MY_DOUBLE*)))==0){
1606 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1607 exit(1);
1608 }
1609 for (i = 0; i <unknowns ; i++)
1610 if((v[i] = ( MY_DOUBLE * ) calloc( (size_t) degree_of_freedom, sizeof(MY_DOUBLE)))==0){
1611 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1612 exit(1);
1613 }
1614
1615
1616
1617 /* now the matrices */
1618 SPM_init(&P,unknowns);
1619 SPM_init(&T,unknowns);
1620 #ifdef DEBUG
1621 printf("solving problem with %d unknowns, %d triangles\n",unknowns,triangles);
1622 #endif
1623
1624 /*build_helmholtz( MY_DOUBLE **P, MY_DOUBLE **T, MY_DOUBLE *x, MY_DOUBLE *y, int NUN, int NE, int NN) */
1625 build_helmholtz_sparse(&P,&T,x,y, unknowns, triangles, total, renumber, re_renumber);
1626
1627 /*solve_helmholtz(MY_DOUBLE **P,MY_DOUBLE **T, MY_DOUBLE **x, MY_INT NUN, MY_INT n_to_find) */
1628 #if defined (USE_LAPACK) || defined (USE_MKL)
1629 if(!use_arpack_solver)
1630 solve_helmholtz_lapack_sparse(&P,&T,v,eig_array,unknowns,degree_of_freedom);
1631 else
1632 arpack_find_eigs(&P,&T,v,eig_array,degree_of_freedom,1);
1633 #else
1634 solve_helmholtz_sparse(&P,&T,v,eig_array,unknowns,degree_of_freedom);
1635 #endif /* !USE_LAPACK || USE_MKL */
1636 /* all eigenvectors are normalized, hence normlization has to be done individually */
1637 /* we copy back elements of eigenvectors to z arrays of each point */
1638 /* in order to do that, we need to change size of those z arrays */
1639 /* find ranges for potentials */
1640 g_maxpot = -1.0;
1641 g_minpot = 1.0;
1642 /* normalize each eigenvector such that the absolute max
1643 value is 1 */
1644 if((max_value_array=( MY_DOUBLE * ) calloc( (size_t) degree_of_freedom, sizeof(MY_DOUBLE)))==0){
1645 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1646 exit(1);
1647 }
1648
1649 for(i=0; i<unknowns;i++) {
1650 for (k=0;k< degree_of_freedom;k++) {
1651 if ( ABS(v[i][k])> max_value_array[k]) max_value_array[k]=ABS(v[i][k]);
1652 }
1653 }
1654 /* find reciprocals of max values for division */
1655 for(k=0;k<degree_of_freedom;k++) {
1656 if(!IS_ZERO(max_value_array[k]))
1657 max_value_array[k]=1/max_value_array[k];
1658 }
1659
1660 /* reallocate memory for z array and copy all eigenvectors */
1661 for ( i = 0; i < unknowns; i++ ) {
1662 if (((Ms.Narray[re_renumber[i]])->z=(MY_DOUBLE *)realloc((void*)(Ms.Narray[re_renumber[i]])->z,(size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0) {
1663 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1664 exit(2);
1665 }
1666 memcpy((void*)(Ms.Narray[re_renumber[i]])->z,(void*)v[i],(size_t)degree_of_freedom*sizeof(MY_DOUBLE));
1667
1668 /* find current limits */
1669 for (k=0;k< degree_of_freedom;k++) {
1670 if ( v[i][k]> g_maxpot ) g_maxpot=v[i][k];
1671 if ( v[i][k]< g_minpot ) g_minpot=v[i][k];
1672 }
1673 }
1674 /* now consider fixed points: fill their z arrays with zeros */
1675 for ( i = unknowns; i < total; i++ ) {
1676 if (((Ms.Narray[re_renumber[i]])->z=(MY_DOUBLE *)realloc((void*)(Ms.Narray[re_renumber[i]])->z,(size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0) {
1677 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1678 exit(2);
1679 }
1680 /* fill zeros */
1681 memset((void*)(Ms.Narray[re_renumber[i]])->z,0,(size_t)degree_of_freedom*sizeof(MY_DOUBLE));
1682 }
1683
1684
1685
1686 /* now we can free the renumber[] array */
1687 free(renumber);
1688 /* now free this array too */
1689 free(re_renumber);
1690
1691 /* free everything else */
1692
1693 /* first the arrays */
1694 free(x);
1695 free(y);
1696
1697 /* next the matrix */
1698 for ( i = 0; i < unknowns ; i++) {
1699 free(v[i]);
1700 }
1701 free(v);
1702 free(max_value_array);
1703 SPM_destroy(&P);
1704 SPM_destroy(&T);
1705
1706 /* now printout a summery */
1707 #ifdef DEBUG
1708 for ( i = 0; i < total; i++ )
1709 printf("renumber_and_solve: %d "MDF" "MDF" "MDF"\n",i,Mx(i,Ms),My(i,Ms),Mz(i,Ms));
1710 #endif
1711
1712 return(0);
1713 }
1714