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: wave.c,v 1.10 2005/04/24 05:49:17 sarod Exp $
17 */
18
19 #ifdef HAVE_CONFIG_H
20 #include <config.h>
21 #endif
22
23 /* formulation of inhomogeneous wave equation */
24 #include <stdlib.h>
25 #include <string.h>
26 #include <stdio.h>
27 #include <math.h>
28 #include "types.h"
29
30 /* cutoff frequency and propagation constant used
31 in waveguide problems */
32 MY_DOUBLE global_wave_k0, global_wave_beta;
33 /* if 1 use full matrix solver, else use half the matrix, for symmetric
34 positive definite eigen equation */
35 int use_full_eigensolver;
36
37 typedef struct hash_edge_{
38 MY_INT p0; /* point 0*/
39 MY_INT p1; /* point 1*/
40 MY_INT n; /* unique number */
41 } hash_edge;
42
43
44 static unsigned int
hashfn(const void * data)45 hashfn(const void *data)
46 {
47 unsigned int k,l;
48 hash_edge *e=(hash_edge*)data;
49
50 k=(int)e->p0;
51 l=(int)e->p1;
52 #ifdef DEBUG
53 printf("hash : key value=%d,%d\n",k,l);
54 #endif
55
56 k+=~(k<<15);
57 l^=(k>>10);
58 k+=(l<<3);
59 l^=(k>>6);
60 k+=~(l<<11);
61 l^=(k>>16);
62 k+=l;
63 #ifdef DEBUG
64 printf("hash : value=%d\n",k);
65 #endif
66 return(k);
67 }
68
69 /* 64 bit hash function */
70 /*long longhash1(long key)
71 {
72 key += ~(key << 32);
73 key ^= (key >>> 22);
74 key += ~(key << 13);
75 key ^= (key >>> 8);
76 key += (key << 3);
77 key ^= (key >>> 15);
78 key += ~(key << 27);
79 key ^= (key >>> 31);
80 return key;
81 } */
82 /* generic list functions */
83
84 /* comparison function */
85 /* compare two edges */
86 int
match(const void * key1,const void * key2)87 match(const void *key1, const void *key2)
88 {
89 hash_edge *e1,*e2;
90 e1=(hash_edge*)key1;
91 e2=(hash_edge*)key2;
92
93 #ifdef DEBUG
94 printf("comparing (%d,%d) with (%d,%d)\n",e1->p0,e1->p1,e2->p0,e2->p1);
95 #endif
96 if((e1->p0==e2->p0) && (e1->p1==e2->p1));
97 return(0);
98 /* else */
99 return(-1);
100 }
101
102 /* freeing function */
103 void
destroy(void * data)104 destroy(void *data)
105 {
106 hash_edge *t=(hash_edge *)data;
107 free(t);
108 }
109
110
111
112 /* if edge (p0-p1) is on a dirichlet boundary, return 1.
113 * else return 0 */
114 static MY_INT
is_this_edge_on_a_dirichlet_boundary(p0,p1)115 is_this_edge_on_a_dirichlet_boundary(p0,p1)
116 {
117 /* find the ratio using x=(x1+lambda. x2)/(1+lambda) */
118 /* if x is on (x1-x2), lambda >0 , equal to y */
119 /* try to do robust computations */
120 int i;
121 double xp0,xp1,xn0,xn1,yp0,yp1,yn0,yn1;
122 xp0=Mx(p0,M);
123 yp0=My(p0,M);
124 xp1=Mx(p1,M);
125 yp1=My(p1,M);
126
127 #ifdef DEBUG
128 printf("is_this_edge_on_a_boundary:consider points(%d,%d)\n",p0,p1);
129 #endif
130 for (i=0; i< nedges; i++) {
131 xn0=Mx(edge_array[i].p1,M);
132 yn0=My(edge_array[i].p1,M);
133 xn1=Mx(edge_array[i].p2,M);
134 yn1=My(edge_array[i].p2,M);
135 #ifdef DEBUG
136 printf("is_this_edge_on_a_boundary: edge %d\n",i);
137 #endif
138 if ( IS_ZERO((xn0-xp0)*(yp0-yn1)-(yn0-yp0)*(xp0-xn1))
139 /* ratio is equal magnitude, check sign */
140 && ( ((yn0-yp0)*(yp0-yn1)>0) /* if true, equal and +ve lambda. stop */
141 || (IS_ZERO((yn0-yp0)*(yp0-yn1))
142 && (xn0-xp0)*(xp0-xn1) >= 0))){
143 /* if we are here, point p0 is on line */
144 /* do the same check for point p1 */
145 #ifdef DEBUG
146 printf("is_this_edge_on_a_boundary: point %d is on line %d\n",p0,i);
147 #endif
148 if ( IS_ZERO((xn0-xp1)*(yp1-yn1)-(yn0-yp1)*(xp1-xn1))
149 /* ratio is equal magnitude, check sign */
150 && ( ((yn0-yp1)*(yp1-yn1)>0) /* if true, equal and +ve lambda. stop */
151 || (IS_ZERO((yn0-yp1)*(yp1-yn1))
152 && (xn0-xp1)*(xp1-xn1) >= 0))){
153 if ( edge_array[i].type == DR ) {
154 return(1);
155 } else {
156 return(0);
157 }
158 }
159 }
160
161 }
162
163 return(0); /* not found */
164
165 }
166
167
168
169
170 static MY_INT
get_edge_number(MY_INT p0,MY_INT p1,htbl H)171 get_edge_number(MY_INT p0,MY_INT p1, htbl H)
172 {
173 hash_edge *ee,*gg;
174
175 if ( (ee=(hash_edge*)malloc(sizeof(hash_edge)))==0) {
176 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
177 exit(1);
178 }
179 if ( p0 > p1 ) {
180 ee->p0=p1;
181 ee->p1=p0;
182 } else {
183 ee->p0=p0;
184 ee->p1=p1;
185 }
186 gg=(hash_edge*)htbl_lookup(&H, (void *)ee);
187 free(ee);
188 if (gg) {
189 return(gg->n);
190 } else {
191 printf("get_edge_number: something SERIOUSLY wrong! (%d,%d)\n",p0,p1);
192 return(-1);
193 }
194 }
195
196 /* solve for cutoff frequency given propagation constant */
197 static void
build_equations_from_hash_table_given_beta(htbl H,SPMat * S,SPMat * T,MY_INT * renumber,MY_INT Ntotal,MY_INT Nedges,MY_DOUBLE beta)198 build_equations_from_hash_table_given_beta(htbl H, SPMat *S, SPMat *T,
199 MY_INT *renumber, MY_INT Ntotal, MY_INT Nedges, MY_DOUBLE beta)
200 {
201 /*
202 S, T - Ntotal by Ntotal matrices
203 Nedges - number of edges, the rest points
204 renumber - map points to matrix location
205 edges are mapper using hash table H
206 */
207 triangle *t;
208 MY_INT i,j,m,n;
209 MY_DOUBLE x[3],y[3],a[3],b[3],c[3],x_bar,y_bar,x_var,y_var;
210 MY_DOUBLE A[3],B[3],C[3],D[3],L[3],Ar,k;
211 MY_DOUBLE S_tt[3][3],S_tz[3][3],S_zt[3][3],S_zz[3][3],T_tt[3][3],T_zz[3][3];
212 MY_DOUBLE I1[3][3],I2[3][3],I3[3][3],I4[3][3],I5[3][3];
213 MY_DOUBLE epsilon,sumI,mu;
214 MY_INT e[3],mult[3];
215
216 DAG_traverse_list_reset(&dt);
217 t=DAG_traverse_prune_list(&dt);
218 while(t) {
219 if (t && t->status==LIVE) {
220 /* epsilon= bound.poly_array[t->boundary].epsilon;
221 mu= bound.poly_array[t->boundary].mu; */
222 if(bound.poly_array[t->boundary].muexp) {
223 mu=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].muexp,t->p0)
224 +ex(bound.poly_array[t->boundary].muexp,t->p1)
225 +ex(bound.poly_array[t->boundary].muexp,t->p2));
226 } else {
227 mu=bound.poly_array[t->boundary].mu;
228 }
229 if(bound.poly_array[t->boundary].epsexp) {
230 epsilon=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].epsexp,t->p0)
231 +ex(bound.poly_array[t->boundary].epsexp,t->p1)
232 +ex(bound.poly_array[t->boundary].epsexp,t->p2));
233 } else {
234 epsilon=bound.poly_array[t->boundary].epsilon;
235 }
236
237 /* get three edges */
238 /* p0
239 /\
240 / \
241 edge 1 / \ edge 3
242 / \
243 /________\
244 p1 edge 2 p2
245 */
246
247 /* get coords */
248 /* convert to global coords */
249 x[0]=Mx(t->p0,M)*g_xscale-g_xoff;
250 y[0]=My(t->p0,M)*g_yscale-g_yoff;
251 x[1]=Mx(t->p1,M)*g_xscale-g_xoff;
252 y[1]=My(t->p1,M)*g_yscale-g_yoff;
253 x[2]=Mx(t->p2,M)*g_xscale-g_xoff;
254 y[2]=My(t->p2,M)*g_yscale-g_yoff;
255
256 x_bar=(x[1]+x[2]+x[0])*ONE_OVER_THREE;
257 y_bar=(y[1]+y[2]+y[0])*ONE_OVER_THREE;
258 x_var=(x[1]*x[1]+x[2]*x[2]+x[0]*x[0]);
259 y_var=(y[1]*y[1]+y[2]*y[2]+y[0]*y[0]);
260
261 a[0]=x[1]*y[2]-x[2]*y[1];
262 a[1]=x[2]*y[0]-x[0]*y[2];
263 a[2]=x[0]*y[1]-x[1]*y[0];
264 b[0]=y[1]-y[2];
265 b[1]=y[2]-y[0];
266 b[2]=y[0]-y[1];
267 c[0]=x[2]-x[1];
268 c[1]=x[0]-x[2];
269 c[2]=x[1]-x[0];
270
271 A[0]=a[0]*b[1]-a[1]*b[0];
272 A[1]=a[1]*b[2]-a[2]*b[1];
273 A[2]=a[2]*b[0]-a[0]*b[2];
274 B[0]=c[0]*b[1]-c[1]*b[0];
275 B[1]=c[1]*b[2]-c[2]*b[1];
276 B[2]=c[2]*b[0]-c[0]*b[2];
277 C[0]=a[0]*c[1]-a[1]*c[0];
278 C[1]=a[1]*c[2]-a[2]*c[1];
279 C[2]=a[2]*c[0]-a[0]*c[2];
280 D[0]=-B[0];
281 D[1]=-B[1];
282 D[2]=-B[2];
283 /* edge lengths */
284 L[0]=sqrt(ABS(c[2]*c[2]+b[2]*b[2]));
285 L[1]=sqrt(ABS(c[0]*c[0]+b[0]*b[0]));
286 L[2]=sqrt(ABS(c[1]*c[1]+b[1]*b[1]));
287 /* area */
288 Ar=0.5*(a[1]+a[2]+a[0]);
289
290 for (m=0; m<3; m++) {
291 for (n=0; n<3; n++) {
292 I1[m][n]=A[m]*A[n]+C[m]*C[n];
293 I2[m][n]=(C[m]*D[n]+C[n]*D[m])*x_bar;
294 I3[m][n]=(A[m]*B[n]+A[n]*B[m])*y_bar;
295 I4[m][n]=ONE_OVER_TWELVE*B[m]*B[n]*(y_var+9*y_bar*y_bar);
296 I5[m][n]=ONE_OVER_TWELVE*D[m]*D[n]*(x_var+9*x_bar*x_bar);
297 }
298 }
299
300
301
302 k=1/(16*Ar*Ar*Ar);
303 for (m=0; m<3; m++) {
304 for (n=0; n<3; n++) {
305 sumI=I1[m][n]+I2[m][n]+I3[m][n]+I4[m][n]+I5[m][n];
306 S_tt[m][n]=k*L[m]*L[n]*(D[m]*D[n]+sumI)/(mu);
307 T_tt[m][n]=k*epsilon*L[m]*L[n]*sumI;
308 }
309 }
310
311 k=beta*beta/(mu*8*Ar*Ar);
312 for (m=0; m<3; m++) {
313 for (n=0; n<3; n++) {
314 S_tz[m][n]=k*L[m]*(b[n]*(A[m]+B[m]*y_bar)+c[n]*(C[m]+D[m]*x_bar));
315 S_zt[n][m]=k*L[m]*(b[n]*(A[m]+B[m]*y_bar)+c[n]*(C[m]+D[m]*x_bar));
316 }
317 }
318
319 k=beta*beta/(mu*4*Ar);
320 for (m=0; m<3; m++) {
321 for (n=0; n<3; n++) {
322 S_zz[m][n]=k*(b[m]*b[n]+c[m]*c[n]);
323 }
324 }
325
326 k=beta*beta*epsilon*Ar*ONE_OVER_TWELVE;
327 for (m=0; m<3; m++) {
328 for (n=0; n<3; n++) {
329 T_zz[m][n]=k;
330 }
331 }
332 T_zz[0][0]=T_zz[1][1]=T_zz[2][2]=2*k; /* 1/6 */
333
334 /* hashes for the edges */
335 e[0]=get_edge_number(t->p0,t->p1,H);
336 e[1]=get_edge_number(t->p1,t->p2,H);
337 e[2]=get_edge_number(t->p2,t->p0,H);
338
339 /*
340 || Stt Stz || || Ttt 0||
341 || Szt Szz || || 0 Tzz||
342 S T
343 Sx = lambda Tx
344 */
345 /* consistencay with sign change */
346 mult[0]=mult[1]=mult[2]=1;
347 if (t->p0>t->p1) {
348 mult[0]=-1;
349 }
350 if (t->p1>t->p2) {
351 mult[1]=-1;
352 }
353 if (t->p2>t->p0) {
354 mult[2]=-1;
355 }
356 /* signs changes
357 Stt, and Ttt - both rows and columns
358 Stz rows
359 Szt columns
360 */
361 for (i=0; i<3; i++) {
362 for (j=0; j<3; j++) {
363 if ( mult[i]*mult[j] < 0 ) {
364 S_tt[i][j]=-S_tt[i][j]; T_tt[i][j]=-T_tt[i][j];
365 }
366 }
367 if ( mult[i] < 0 ) {
368 for (j=0; j<3; j++) {
369 S_tz[i][j]=-S_tz[i][j];
370 S_zt[j][i]=-S_zt[j][i];
371 }
372 }
373 }
374
375 /* assemble global matrix */
376 /* if any of the edges are Dirichlet, we do not include them in the
377 global matrix */
378 #ifdef DEBUG
379 printf("build_equations: triangle %d, (%d,%d,%d)->(%d,%d,%d)/%d, edges (%d,%d,%d)\n",t->n,t->p0,t->p1,t->p2,renumber[t->p0],renumber[t->p1],renumber[t->p2], Nedges,e[0],e[1],e[2]);
380 for (i=0;i<3;i++) {
381 printf("||");
382 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",S_tt[i][0],S_tt[i][1],S_tt[i][2],S_tz[i][0],S_tz[i][1],S_tz[i][2]);
383 printf("|| ||");
384 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",T_tt[i][0],T_tt[i][1],T_tt[i][2],0.0,0.0,0.0);
385 printf("||\n");
386 }
387 for (i=0;i<3;i++) {
388 printf("||");
389 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",S_zt[i][0],S_zt[i][1],S_zt[i][2],S_zz[i][0],S_zz[i][1],S_zz[i][2]);
390 printf("|| ||");
391 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",0.0,0.0,0.0,T_zz[i][0],T_zz[i][1],T_zz[i][2]);
392 printf("||\n");
393 }
394 #endif
395
396 /* assembling global matrix */
397 for (i=0; i<3; i++) {
398 if ( e[i] != -1 ) {
399 for (j=0;j<3;j++) {
400 if (e[j] != -1 ) {
401 /*T[e[i]][e[j]]+=(T_tt[i][j]); */
402 SPM_add(T,e[i],e[j],T_tt[i][j]);
403 /*S[e[i]][e[j]]+=(S_tt[i][j]); */
404 SPM_add(S,e[i],e[j],S_tt[i][j]);
405 }
406 }
407 }
408 }
409 i=renumber[t->p0];
410 if (i>=Nedges) { /* fixed points have i == 0 */
411 for (j=0; j<3; j++) {
412 if ( e[j]!=-1) { SPM_add(S,i,e[j],S_zt[0][j]); }
413 if ( e[j]!=-1) { SPM_add(S,e[j],i,S_tz[j][0]); }
414 }
415 }
416 i=renumber[t->p1];
417 if (i>=Nedges) { /* fixed points have i == 0 */
418 for (j=0; j<3; j++) {
419 if ( e[j]!=-1) { SPM_add(S,i,e[j],S_zt[1][j]); }
420 if ( e[j]!=-1) { SPM_add(S,e[j],i,S_tz[j][1]); }
421 }
422 }
423 i=renumber[t->p2];
424 if (i>=Nedges) { /* fixed points have i == 0 */
425 for (j=0; j<3; j++) {
426 if ( e[j]!=-1) { SPM_add(S,i,e[j],S_zt[2][j]); }
427 if ( e[j]!=-1) { SPM_add(S,e[j],i,S_tz[j][2]); }
428 }
429 }
430 /* user e array to store point numbers temporarily */
431 e[0]=renumber[t->p0];e[1]=renumber[t->p1];e[2]=renumber[t->p2];
432 for (i=0; i<3; i++) {
433 if ( e[i] >=Nedges ) {
434 for (j=0;j<3; j++) {
435 if ( e[j] >=Nedges ) {
436 /*S[e[i]][e[j]]+=S_zz[i][j]; */
437 SPM_add(S,e[i],e[j],S_zz[i][j]);
438 /*T[e[i]][e[j]]+=T_zz[i][j]; */
439 SPM_add(T,e[i],e[j],T_zz[i][j]);
440 }
441 }
442 }
443 }
444
445 }
446
447 t=DAG_traverse_prune_list(&dt);
448 }
449 }
450
451 /* solve for propagation constant given frequency */
452 static void
build_equations_from_hash_table_given_freq(htbl H,SPMat * S,SPMat * T,MY_INT * renumber,MY_INT Ntotal,MY_INT Nedges,MY_DOUBLE k0)453 build_equations_from_hash_table_given_freq(htbl H, SPMat *S, SPMat *T,
454 MY_INT *renumber, MY_INT Ntotal, MY_INT Nedges, MY_DOUBLE k0)
455 {
456 triangle *t;
457 MY_INT i,j,m,n;
458 MY_DOUBLE x[3],y[3],a[3],b[3],c[3],x_bar,y_bar,x_var,y_var;
459 MY_DOUBLE A[3],B[3],C[3],D[3],L[3],Ar,k;
460 MY_DOUBLE S_tt[3][3],T_tz[3][3],T_zt[3][3],T_zz[3][3],T_tt[3][3],S_zz[3][3];
461 MY_DOUBLE I1[3][3],I2[3][3],I3[3][3],I4[3][3],I5[3][3];
462 MY_DOUBLE epsilon,sumI,mu;
463 MY_INT e[3],mult[3];
464
465
466 DAG_traverse_list_reset(&dt);
467 t=DAG_traverse_prune_list(&dt);
468 while(t) {
469 if (t && t->status==LIVE) {
470 /* epsilon = bound.poly_array[t->boundary].epsilon;
471 mu= bound.poly_array[t->boundary].mu; */
472 if(bound.poly_array[t->boundary].muexp) {
473 mu=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].muexp,t->p0)
474 +ex(bound.poly_array[t->boundary].muexp,t->p1)
475 +ex(bound.poly_array[t->boundary].muexp,t->p2));
476 } else {
477 mu=bound.poly_array[t->boundary].mu;
478 }
479 if(bound.poly_array[t->boundary].epsexp) {
480 epsilon=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].epsexp,t->p0)
481 +ex(bound.poly_array[t->boundary].epsexp,t->p1)
482 +ex(bound.poly_array[t->boundary].epsexp,t->p2));
483 } else {
484 epsilon=bound.poly_array[t->boundary].epsilon;
485 }
486
487 /* get three edges */
488 /* p0
489 /\
490 / \
491 edge 1 / \ edge 3
492 / \
493 /________\
494 p1 edge 2 p2
495 */
496
497 /* get coords */
498 /* convert to global coords */
499 x[0]=Mx(t->p0,M)*g_xscale-g_xoff;
500 y[0]=My(t->p0,M)*g_yscale-g_yoff;
501 x[1]=Mx(t->p1,M)*g_xscale-g_xoff;
502 y[1]=My(t->p1,M)*g_yscale-g_yoff;
503 x[2]=Mx(t->p2,M)*g_xscale-g_xoff;
504 y[2]=My(t->p2,M)*g_yscale-g_yoff;
505
506 x_bar=(x[1]+x[2]+x[0])*ONE_OVER_THREE;
507 y_bar=(y[1]+y[2]+y[0])*ONE_OVER_THREE;
508 x_var=(x[1]*x[1]+x[2]*x[2]+x[0]*x[0]);
509 y_var=(y[1]*y[1]+y[2]*y[2]+y[0]*y[0]);
510
511 a[0]=x[1]*y[2]-x[2]*y[1];
512 a[1]=x[2]*y[0]-x[0]*y[2];
513 a[2]=x[0]*y[1]-x[1]*y[0];
514 b[0]=y[1]-y[2];
515 b[1]=y[2]-y[0];
516 b[2]=y[0]-y[1];
517 c[0]=x[2]-x[1];
518 c[1]=x[0]-x[2];
519 c[2]=x[1]-x[0];
520
521 A[0]=a[0]*b[1]-a[1]*b[0];
522 A[1]=a[1]*b[2]-a[2]*b[1];
523 A[2]=a[2]*b[0]-a[0]*b[2];
524 B[0]=c[0]*b[1]-c[1]*b[0];
525 B[1]=c[1]*b[2]-c[2]*b[1];
526 B[2]=c[2]*b[0]-c[0]*b[2];
527 C[0]=a[0]*c[1]-a[1]*c[0];
528 C[1]=a[1]*c[2]-a[2]*c[1];
529 C[2]=a[2]*c[0]-a[0]*c[2];
530 D[0]=-B[0];
531 D[1]=-B[1];
532 D[2]=-B[2];
533 /* edge lengths */
534 L[0]=sqrt(ABS(c[2]*c[2]+b[2]*b[2]));
535 L[1]=sqrt(ABS(c[0]*c[0]+b[0]*b[0]));
536 L[2]=sqrt(ABS(c[1]*c[1]+b[1]*b[1]));
537 /* area */
538 Ar=0.5*(a[1]+a[2]+a[0]);
539
540 for (m=0; m<3; m++) {
541 for (n=0; n<3; n++) {
542 I1[m][n]=A[m]*A[n]+C[m]*C[n];
543 I2[m][n]=(C[m]*D[n]+C[n]*D[m])*x_bar;
544 I3[m][n]=(A[m]*B[n]+A[n]*B[m])*y_bar;
545 I4[m][n]=ONE_OVER_TWELVE*B[m]*B[n]*(y_var+9*y_bar*y_bar);
546 I5[m][n]=ONE_OVER_TWELVE*D[m]*D[n]*(x_var+9*x_bar*x_bar);
547 }
548 }
549
550
551
552 k=1/(16*Ar*Ar*Ar);
553 for (m=0; m<3; m++) {
554 for (n=0; n<3; n++) {
555 sumI=I1[m][n]+I2[m][n]+I3[m][n]+I4[m][n]+I5[m][n];
556 S_tt[m][n]=k*L[m]*L[n]*(D[m]*D[n]-k0*k0*epsilon*sumI)/(mu);
557 T_tt[m][n]=k*epsilon*L[m]*L[n]*sumI;
558 }
559 }
560
561 k=1/(mu*8*Ar*Ar);
562 for (m=0; m<3; m++) {
563 for (n=0; n<3; n++) {
564 T_tz[m][n]=k*L[m]*(b[n]*(A[m]+B[m]*y_bar)+c[n]*(C[m]+D[m]*x_bar));
565 T_zt[n][m]=k*L[m]*(b[n]*(A[m]+B[m]*y_bar)+c[n]*(C[m]+D[m]*x_bar));
566 }
567 }
568
569 k=1/(mu*4*Ar);
570 for (m=0; m<3; m++) {
571 for (n=0; n<3; n++) {
572 S_zz[m][n]=k*(b[m]*b[n]+c[m]*c[n]);
573 }
574 }
575
576 k=k0*k0*epsilon*Ar*ONE_OVER_TWELVE;
577 for (m=0; m<3; m++) {
578 for (n=0; n<3; n++) {
579 T_zz[m][n]=k;
580 }
581 }
582 T_zz[0][0]=T_zz[1][1]=T_zz[2][2]=2*k; /* 1/6 */
583 for (m=0; m<3; m++) {
584 for (n=0; n<3; n++) {
585 T_zz[m][n]+=S_zz[m][n];
586 }
587 }
588
589
590
591
592 /* hashes for the edges */
593 e[0]=get_edge_number(t->p0,t->p1,H);
594 e[1]=get_edge_number(t->p1,t->p2,H);
595 e[2]=get_edge_number(t->p2,t->p0,H);
596
597 /* consistencay with sign change */
598 mult[0]=mult[1]=mult[2]=1;
599 if (t->p0>t->p1) {
600 mult[0]=-1;
601 }
602 if (t->p1>t->p2) {
603 mult[1]=-1;
604 }
605 if (t->p2>t->p0) {
606 mult[2]=-1;
607 }
608 for (i=0; i<3; i++) {
609 for (j=0; j<3; j++) {
610 if ( mult[i]*mult[j] < 0 ) {
611 S_tt[i][j]=-S_tt[i][j]; T_tt[i][j]=-T_tt[i][j];
612 }
613 }
614 if ( mult[i] < 0 ) {
615 for (j=0; j<3; j++) {
616 T_tz[i][j]=-T_tz[i][j];
617 T_zt[j][i]=-T_zt[j][i];
618 }
619 }
620 }
621 /* assemble global matrix */
622 /* if any of the edges are Dirichlet, we do not include them in the
623 global matrix */
624 #ifdef DEBUG
625 printf("build_equations: triangle %d, (%d,%d,%d)->(%d,%d,%d)/%d, edges (%d,%d,%d)\n",t->n,t->p0,t->p1,t->p2,renumber[t->p0],renumber[t->p1],renumber[t->p2], Nedges,e[0],e[1],e[2]);
626 for (i=0;i<3;i++) {
627 printf("||");
628 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",S_tt[i][0],S_tt[i][1],S_tt[i][2],0.0,0.0,0.0);
629 printf("|| ||");
630 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",T_tt[i][0],T_tt[i][1],T_tt[i][2],T_tz[i][2],T_tz[i][2],T_tz[i][2]);
631 printf("||\n");
632 }
633 for (i=0;i<3;i++) {
634 printf("||");
635 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",0.0,0.0,0.0,0.0,0.0,0.0);
636 printf("|| ||");
637 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",T_zt[i][0],T_zt[i][1],T_zt[i][2],T_zz[i][2],T_zz[i][2],T_zz[i][2]);
638 printf("||\n");
639 }
640 #endif
641
642 /* assembling global matrix */
643 for (i=0; i<3; i++) {
644 if ( e[i] != -1 ) {
645 for (j=0;j<3;j++) {
646 if (e[j] != -1 ) {
647 /*S[e[i]][e[j]]+=(S_tt[i][j]); */
648 SPM_add(S,e[i],e[j],S_tt[i][j]);
649 /*T[e[i]][e[j]]+=(T_tt[i][j]); */
650 SPM_add(T,e[i],e[j],T_tt[i][j]);
651 }
652 }
653 }
654 }
655 i=renumber[t->p0];
656 if (i>=Nedges) { /* fixed points have i == 0 */
657 for (j=0; j<3; j++) {
658 if ( e[j]!=-1) { SPM_add(T,i,e[j],T_tz[0][j]); }
659 if ( e[j]!=-1) { SPM_add(T,e[j],i,T_zt[j][0]); }
660 }
661 }
662 i=renumber[t->p1];
663 if (i>=Nedges) { /* fixed points have i == 0 */
664 for (j=0; j<3; j++) {
665 if ( e[j]!=-1) { SPM_add(T,i,e[j],T_tz[1][j]); }
666 if ( e[j]!=-1) { SPM_add(T,e[j],i,T_zt[j][1]); }
667 }
668 }
669 i=renumber[t->p2];
670 if (i>=Nedges) { /* fixed points have i == 0 */
671 for (j=0; j<3; j++) {
672 if ( e[j]!=-1) { SPM_add(T,i,e[j],T_tz[2][j]); }
673 if ( e[j]!=-1) { SPM_add(T,e[j],i,T_zt[j][2]); }
674 }
675 }
676 /* user e array to store point numbers temporarily */
677 e[0]=renumber[t->p0];e[1]=renumber[t->p1];e[2]=renumber[t->p2];
678 for (i=0; i<3; i++) {
679 if ( e[i] >=Nedges ) {
680 for (j=0;j<3; j++) {
681 if ( e[j] >=Nedges ) {
682 SPM_add(T,e[i],e[j],T_zz[i][j]);
683 }
684 }
685 }
686 }
687
688
689 }
690
691 t=DAG_traverse_prune_list(&dt);
692 }
693 }
694
695
696 /* solve for cutoff frequency */
697 static void
build_equations_from_hash_table(htbl H,SPMat * S,SPMat * T,MY_INT * renumber,MY_INT Ntotal,MY_INT Nedges)698 build_equations_from_hash_table(htbl H, SPMat *S, SPMat *T,
699 MY_INT *renumber, MY_INT Ntotal, MY_INT Nedges)
700 {
701 /*
702 S, T - Ntotal by Ntotal matrices
703 Nedges - number of edges, the rest points
704 renumber - map points to matrix location
705 edges are mapper using hash table H
706 */
707 triangle *t;
708 MY_INT i,j,m,n;
709 MY_DOUBLE x[3],y[3],a[3],b[3],c[3],x_bar,y_bar,x_var,y_var;
710 MY_DOUBLE A[3],B[3],C[3],D[3],L[3],Ar,k;
711 MY_DOUBLE S_tt[3][3],S_zz[3][3],T_tt[3][3],T_zz[3][3];
712 MY_DOUBLE I1[3][3],I2[3][3],I3[3][3],I4[3][3],I5[3][3];
713 MY_DOUBLE epsilon,sumI,mu;
714 MY_INT e[3],mult[3];
715
716 DAG_traverse_list_reset(&dt);
717 t=DAG_traverse_prune_list(&dt);
718 while(t) {
719 if (t && t->status==LIVE) {
720 /* epsilon= bound.poly_array[t->boundary].epsilon;
721 mu= bound.poly_array[t->boundary].mu; */
722 if(bound.poly_array[t->boundary].muexp) {
723 mu=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].muexp,t->p0)
724 +ex(bound.poly_array[t->boundary].muexp,t->p1)
725 +ex(bound.poly_array[t->boundary].muexp,t->p2));
726 } else {
727 mu=bound.poly_array[t->boundary].mu;
728 }
729 if(bound.poly_array[t->boundary].epsexp) {
730 epsilon=ONE_OVER_THREE*(ex(bound.poly_array[t->boundary].epsexp,t->p0)
731 +ex(bound.poly_array[t->boundary].epsexp,t->p1)
732 +ex(bound.poly_array[t->boundary].epsexp,t->p2));
733 } else {
734 epsilon=bound.poly_array[t->boundary].epsilon;
735 }
736
737 /* get three edges */
738 /* p0
739 /\
740 / \
741 edge 1 / \ edge 3
742 / \
743 /________\
744 p1 edge 2 p2
745 */
746
747 /* get coords */
748 /* convert to global coords */
749 x[0]=Mx(t->p0,M)*g_xscale-g_xoff;
750 y[0]=My(t->p0,M)*g_yscale-g_yoff;
751 x[1]=Mx(t->p1,M)*g_xscale-g_xoff;
752 y[1]=My(t->p1,M)*g_yscale-g_yoff;
753 x[2]=Mx(t->p2,M)*g_xscale-g_xoff;
754 y[2]=My(t->p2,M)*g_yscale-g_yoff;
755
756 x_bar=(x[1]+x[2]+x[0])*ONE_OVER_THREE;
757 y_bar=(y[1]+y[2]+y[0])*ONE_OVER_THREE;
758 x_var=(x[1]*x[1]+x[2]*x[2]+x[0]*x[0]);
759 y_var=(y[1]*y[1]+y[2]*y[2]+y[0]*y[0]);
760
761 a[0]=x[1]*y[2]-x[2]*y[1];
762 a[1]=x[2]*y[0]-x[0]*y[2];
763 a[2]=x[0]*y[1]-x[1]*y[0];
764 b[0]=y[1]-y[2];
765 b[1]=y[2]-y[0];
766 b[2]=y[0]-y[1];
767 c[0]=x[2]-x[1];
768 c[1]=x[0]-x[2];
769 c[2]=x[1]-x[0];
770
771 A[0]=a[0]*b[1]-a[1]*b[0];
772 A[1]=a[1]*b[2]-a[2]*b[1];
773 A[2]=a[2]*b[0]-a[0]*b[2];
774 B[0]=c[0]*b[1]-c[1]*b[0];
775 B[1]=c[1]*b[2]-c[2]*b[1];
776 B[2]=c[2]*b[0]-c[0]*b[2];
777 C[0]=a[0]*c[1]-a[1]*c[0];
778 C[1]=a[1]*c[2]-a[2]*c[1];
779 C[2]=a[2]*c[0]-a[0]*c[2];
780 D[0]=-B[0];
781 D[1]=-B[1];
782 D[2]=-B[2];
783 /* edge lengths */
784 L[0]=sqrt(ABS(c[2]*c[2]+b[2]*b[2]));
785 L[1]=sqrt(ABS(c[0]*c[0]+b[0]*b[0]));
786 L[2]=sqrt(ABS(c[1]*c[1]+b[1]*b[1]));
787 /* area */
788 Ar=0.5*(a[1]+a[2]+a[0]);
789
790 for (m=0; m<3; m++) {
791 for (n=0; n<3; n++) {
792 I1[m][n]=A[m]*A[n]+C[m]*C[n];
793 I2[m][n]=(C[m]*D[n]+C[n]*D[m])*x_bar;
794 I3[m][n]=(A[m]*B[n]+A[n]*B[m])*y_bar;
795 I4[m][n]=ONE_OVER_TWELVE*B[m]*B[n]*(y_var+9*y_bar*y_bar);
796 I5[m][n]=ONE_OVER_TWELVE*D[m]*D[n]*(x_var+9*x_bar*x_bar);
797 }
798 }
799
800
801
802 k=1/(16*Ar*Ar*Ar);
803 for (m=0; m<3; m++) {
804 for (n=0; n<3; n++) {
805 sumI=I1[m][n]+I2[m][n]+I3[m][n]+I4[m][n]+I5[m][n];
806 S_tt[m][n]=4*k*L[m]*L[n]*(D[m]*D[n])/(mu);
807 T_tt[m][n]=k*epsilon*L[m]*L[n]*sumI;
808 }
809 }
810
811 k=1/(mu*4*Ar);
812 for (m=0; m<3; m++) {
813 for (n=0; n<3; n++) {
814 S_zz[m][n]=k*(b[m]*b[n]+c[m]*c[n]);
815 }
816 }
817
818 k=epsilon*Ar*ONE_OVER_TWELVE;
819 for (m=0; m<3; m++) {
820 for (n=0; n<3; n++) {
821 T_zz[m][n]=k;
822 }
823 }
824 T_zz[0][0]=T_zz[1][1]=T_zz[2][2]=2*k; /* 1/6 */
825
826 /* hashes for the edges */
827 e[0]=get_edge_number(t->p0,t->p1,H);
828 e[1]=get_edge_number(t->p1,t->p2,H);
829 e[2]=get_edge_number(t->p2,t->p0,H);
830
831 /*
832 || Stt 0 || || Ttt 0||
833 || 0 Szz || || 0 Tzz||
834 S T
835 Sx = lambda Tx
836 */
837 /* consistencay with sign change */
838 mult[0]=mult[1]=mult[2]=1;
839 if (t->p0>t->p1) {
840 mult[0]=-1;
841 }
842 if (t->p1>t->p2) {
843 mult[1]=-1;
844 }
845 if (t->p2>t->p0) {
846 mult[2]=-1;
847 }
848 /* signs changes
849 Stt, and Ttt - both rows and columns
850 Stz rows
851 Szt columns
852 */
853 for (i=0; i<3; i++) {
854 for (j=0; j<3; j++) {
855 if ( mult[i]*mult[j] < 0 ) {
856 S_tt[i][j]=-S_tt[i][j]; T_tt[i][j]=-T_tt[i][j];
857 }
858 }
859 }
860
861 /* assemble global matrix */
862 /* if any of the edges are Dirichlet, we do not include them in the
863 global matrix */
864 #ifdef DEBUG
865 printf("build_equations: triangle %d, (%d,%d,%d)->(%d,%d,%d)/%d, edges (%d,%d,%d)\n",t->n,t->p0,t->p1,t->p2,renumber[t->p0],renumber[t->p1],renumber[t->p2], Nedges,e[0],e[1],e[2]);
866 for (i=0;i<3;i++) {
867 printf("||");
868 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",S_tt[i][0],S_tt[i][1],S_tt[i][2],0.0,0.0,0.0);
869 printf("|| ||");
870 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",T_tt[i][0],T_tt[i][1],T_tt[i][2],0.0,0.0,0.0);
871 printf("||\n");
872 }
873 for (i=0;i<3;i++) {
874 printf("||");
875 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",0.0,0.0,0.0,S_zz[i][0],S_zz[i][1],S_zz[i][2]);
876 printf("|| ||");
877 printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",0.0,0.0,0.0,T_zz[i][0],T_zz[i][1],T_zz[i][2]);
878 printf("||\n");
879 }
880 #endif
881
882 /* assembling global matrix */
883 for (i=0; i<3; i++) {
884 if ( e[i] != -1 ) {
885 for (j=0;j<3;j++) {
886 if (e[j] != -1 ) {
887 /*T[e[i]][e[j]]+=(T_tt[i][j]); */
888 SPM_add(T,e[i],e[j],T_tt[i][j]);
889 /*S[e[i]][e[j]]+=(S_tt[i][j]); */
890 SPM_add(S,e[i],e[j],S_tt[i][j]);
891 }
892 }
893 }
894 }
895
896 /* user e array to store point numbers temporarily */
897 e[0]=renumber[t->p0];e[1]=renumber[t->p1];e[2]=renumber[t->p2];
898 for (i=0; i<3; i++) {
899 if ( e[i] >=Nedges ) {
900 for (j=0;j<3; j++) {
901 if ( e[j] >=Nedges ) {
902 /*S[e[i]][e[j]]+=S_zz[i][j]; */
903 SPM_add(S,e[i],e[j],S_zz[i][j]);
904 /*T[e[i]][e[j]]+=T_zz[i][j]; */
905 SPM_add(T,e[i],e[j],T_zz[i][j]);
906 }
907 }
908 }
909 }
910
911 }
912
913 t=DAG_traverse_prune_list(&dt);
914 }
915 }
916
917 int
buildup_hash_table_and_equation(Mesh Ms,DAG_tree Dt)918 buildup_hash_table_and_equation(Mesh Ms, DAG_tree Dt) {
919
920 triangle *t;
921 htbl H;
922 hash_edge *ee,*eep;
923 #ifdef DEBUG
924 hash_edge *gg;
925 #endif
926 MY_INT j,Ni,temp;
927 unsigned int i;
928 SPMat S, T;
929 MY_DOUBLE **v;
930 MY_INT *renumber, *re_renumber;
931 MY_DOUBLE *max_value;
932 MY_INT NN;
933
934 MY_DOUBLE x[3],y[3],a[3],b[3],c[3],et[2],Ar;
935 MY_DOUBLE A[3],B[3],C[3],D[3],L[3];
936 MY_INT e1,e2,e3;
937
938 /* make degree of freedom=1 */
939 degree_of_freedom=requested_degree_of_freedom;
940 /* find size needed for hash table */
941 i=Ms.count;
942 #ifdef DEBUG
943 printf("buildup_hash: creating table size %d\n",i);
944 #endif
945 htbl_init(&H, i, sizeof(hash_edge), &hashfn, &match, &destroy);
946
947 /* create renumbering array for points */
948 /* keep point numbers also in edge number array.
949 we first go through the edges and number them.
950 next we run another triangle loop to number the points.
951 We alwasy remember last edge number. So everthyng after
952 that is a point. */
953 i=0;
954 if ( (ee=(hash_edge*)malloc(sizeof(hash_edge)))==0) {
955 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
956 exit(1);
957 }
958
959 /* insert data */
960 DAG_traverse_list_reset(&Dt);
961 t=DAG_traverse_prune_list(&Dt);
962 while (t) {
963 if ( t && t->status==LIVE ) {
964 ee->p0=t->p0;
965 ee->p1=t->p1;
966 /* always p0 <= p1 */
967 if ( ee->p0 > ee->p1 ) {
968 temp=ee->p0;
969 ee->p0=ee->p1;
970 ee->p1=temp;
971 }
972 #ifdef DEBUG
973 printf("**build: tring to insert (%d,%d) into table\n",ee->p0,ee->p1);
974 #endif
975
976 if ( (!htbl_insert(&H, (const void *)ee)) ) {
977 #ifdef DEBUG
978 printf("sucessfully inserted (%d,%d) into table\n",ee->p0,ee->p1);
979 #endif
980 /* check whether this edge is on a Dirichlet boundary */
981 eep=(hash_edge*)htbl_lookup(&H, (void *)ee);
982 if ( !is_this_edge_on_a_dirichlet_boundary(ee->p0,ee->p1) ) {
983 /* give a unique contiguous number */
984 eep->n=i++;
985 } else {
986 #ifdef DEBUG
987 printf("buildup_hash: edge on a Dirichlet boundary\n");
988 #endif
989 eep->n=-1;
990 }
991
992 } else {
993 #ifdef DEBUG
994 printf("failed inserting (%d,%d) into table ",ee->p0,ee->p1);
995 /* lookup */
996 if ( (gg=(hash_edge*)htbl_lookup(&H, (void *)ee)) ) {
997 printf("it exists as (%d,%d)->%d.\n",gg->p0,gg->p1,gg->n);
998 }
999 #endif
1000 }
1001 ee->p0=t->p1;
1002 ee->p1=t->p2;
1003 /* always p0 <= p1 */
1004 if ( ee->p0 > ee->p1 ) {
1005 temp=ee->p0;
1006 ee->p0=ee->p1;
1007 ee->p1=temp;
1008 }
1009 #ifdef DEBUG
1010 printf("**build: tring to insert (%d,%d) into table\n",ee->p0,ee->p1);
1011 #endif
1012
1013 if ( (!htbl_insert(&H, (const void *)ee)) ) {
1014 #ifdef DEBUG
1015 printf("sucessfully inserted (%d,%d) into table\n",ee->p0,ee->p1);
1016 #endif
1017 eep=(hash_edge*)htbl_lookup(&H, (void *)ee);
1018 if ( !is_this_edge_on_a_dirichlet_boundary(ee->p0,ee->p1) ) {
1019 /* give a unique contiguous number */
1020 eep->n=i++;
1021 } else {
1022 #ifdef DEBUG
1023 printf("buildup_hash: edge on a Dirichlet boundary\n");
1024 #endif
1025 eep->n=-1;
1026 }
1027 } else {
1028 #ifdef DEBUG
1029 printf("failed inserting (%d,%d) into table ",ee->p0,ee->p1);
1030 /* lookup */
1031 if ( (gg=(hash_edge*)htbl_lookup(&H, (void *)ee)) ) {
1032 printf("it exists as (%d,%d)->%d.\n",gg->p0,gg->p1,gg->n);
1033 }
1034 #endif
1035 }
1036 ee->p0=t->p2;
1037 ee->p1=t->p0;
1038 /* always p0 <= p1 */
1039 if ( ee->p0 > ee->p1 ) {
1040 temp=ee->p0;
1041 ee->p0=ee->p1;
1042 ee->p1=temp;
1043 }
1044 #ifdef DEBUG
1045 printf("**build: tring to insert (%d,%d) into table\n",ee->p0,ee->p1);
1046 #endif
1047
1048 if ( (!htbl_insert(&H, (const void *)ee)) ) {
1049 #ifdef DEBUG
1050 printf("sucessfully inserted (%d,%d) into table\n",ee->p0,ee->p1);
1051 #endif
1052 eep=(hash_edge*)htbl_lookup(&H, (void *)ee);
1053 if ( !is_this_edge_on_a_dirichlet_boundary(ee->p0,ee->p1) ) {
1054 /* give a unique contiguous number */
1055 eep->n=i++;
1056 } else {
1057 #ifdef DEBUG
1058 printf("buildup_hash: edge on a Dirichlet boundary\n");
1059 #endif
1060 eep->n=-1;
1061 }
1062
1063 } else {
1064 #ifdef DEBUG
1065 printf("failed inserting (%d,%d) into table ",ee->p0,ee->p1);
1066 /* lookup */
1067 if ( (gg=(hash_edge*)htbl_lookup(&H, (void *)ee)) ) {
1068 printf("it exists as (%d,%d)->%d.\n",gg->p0,gg->p1,gg->n);
1069 }
1070 #endif
1071 }
1072 }
1073
1074 t=DAG_traverse_prune_list(&Dt);
1075 }
1076 #ifdef DEBUG
1077 printf("buildup_hash %d of %d edges, hash table filled %lf\n",i,H.positions,H.size/(double)H.positions);
1078 #endif
1079 /* remember no of edges */
1080 Ni=i;
1081
1082 /******************/
1083 /* setup point to matrix mapping array */
1084 if ((renumber = (MY_INT *)calloc(Ms.count,sizeof(MY_INT )))==0){
1085 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1086 exit(1);
1087 }
1088 if((re_renumber = (MY_INT *)calloc(Ms.count,sizeof(MY_INT)))==0){
1089 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1090 exit(1);
1091 }
1092
1093 /* need to solve for only points referred by triangles */
1094 /* SO -- */
1095 DAG_traverse_list_reset(&Dt);
1096 t=DAG_traverse_prune_list(&Dt);
1097 while ( t ) {
1098 if ( t && t->status==LIVE ) {
1099 /* mark the points in this triangle */
1100 /* use re_renumber[] array for this */
1101 re_renumber[t->p0] = 1;
1102 re_renumber[t->p1] = 1;
1103 re_renumber[t->p2] = 1;
1104
1105 }
1106 t=DAG_traverse_prune_list(&Dt);
1107 }
1108
1109 /* now give numbers to variable points*/
1110 /* start numbering from last edge */
1111 j = Ni;
1112 for ( i = 0; i < Ms.count; i++ )
1113 if ( Mval(i,Ms) == VAR && re_renumber[i] == 1) {
1114 renumber[i] = j++;
1115 #ifdef DEBUG
1116 printf("buildup_hash: renumbering variable point %d as %d\n",i,renumber[i]);
1117 #endif
1118 }
1119 free(re_renumber);
1120 /******************/
1121 /* update total size of matrix */
1122 NN=j;
1123 /* Ni edges and NN-Ni points */
1124 /* set up matrices */
1125 SPM_init(&S,NN);
1126 SPM_init(&T,NN);
1127
1128 /* depending on equation type, we formulate a different problem */
1129 if (solve_equation == HELMHOLTZ_INHOMO ) {
1130 build_equations_from_hash_table(H, &S, &T, renumber, NN, Ni);
1131 } else if (solve_equation == HELMHOLTZ_FREQ) {
1132 build_equations_from_hash_table_given_freq(H, &S, &T, renumber, NN, Ni, global_wave_k0);
1133 } else if (solve_equation == HELMHOLTZ_BETA) {
1134 build_equations_from_hash_table_given_beta(H, &S, &T, renumber, NN, Ni, global_wave_beta);
1135 }
1136
1137 #ifdef DEBUG
1138 printf("====================\n");
1139 printf("S=\n");
1140 for (i=0;i<NN;i++) {
1141 printf("|");
1142 for (j=0;j<NN;j++) {
1143 if(SPM_e(&S,i,j)) {printf(" "MDF,SPM_e(&S,i,j));} else {printf(" 0");}
1144 }
1145 printf("|\n");
1146 }
1147 printf("====================\n");
1148 printf("T=\n");
1149 for (i=0;i<NN;i++) {
1150 printf("|");
1151 for (j=0;j<NN;j++) {
1152 if(SPM_e(&T,i,j)) {printf(" "MDF,SPM_e(&T,i,j));} else {printf(" 0");}
1153 }
1154 printf("|\n");
1155 }
1156 #endif
1157
1158 /* structure to store all eigenvectors - matrix */
1159 if((v=(MY_DOUBLE **)calloc(NN, sizeof(MY_DOUBLE*)))==0){
1160 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1161 exit(1);
1162 }
1163 for (i = 0; i <(unsigned)NN; i++)
1164 if((v[i] = ( MY_DOUBLE * ) calloc( (size_t) degree_of_freedom, sizeof(MY_DOUBLE)))==0){
1165 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1166 exit(1);
1167 }
1168
1169 /* now solve eigenvalue problem (size Ni by Ni )
1170 S. x=k^2 T.x
1171 or (S-lambda T)x= 0
1172 call
1173 solve_helmholtz(MY_DOUBLE **P,MY_DOUBLE **T, MY_DOUBLE **x, MY_INT NUN, MY_INT n_to_find)
1174 */
1175
1176 #if defined (USE_LAPACK) || defined (USE_MKL)
1177 if( use_arpack_solver ) {
1178 arpack_find_eigs(&S,&T,v,eig_array,degree_of_freedom,0);
1179 } else {
1180 if(use_full_eigensolver) {
1181 solve_generalized_eig_lapack(S,T,v,eig_array, NN,degree_of_freedom);
1182 }else {
1183 solve_helmholtz_lapack_sparse(&S,&T,v,eig_array,NN,degree_of_freedom);
1184 }
1185 }
1186 #endif
1187
1188 /* reallocate memory for z array and copy all eigenvectors */
1189 /* we need 2 positions for x,y components per degree of freedom */
1190 /* and one more for the z component */
1191 for ( i = 0; i < Ms.count; i++ ) {
1192 if (((Ms.Narray[i])->z=(MY_DOUBLE *)realloc((void*)(Ms.Narray[i])->z,(size_t)3*degree_of_freedom*sizeof(MY_DOUBLE)))==0) {
1193 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1194 exit(2);
1195 }
1196 }
1197
1198 /* array to remember max values of potentials */
1199 if ((max_value=(MY_DOUBLE *)malloc((size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0) {
1200 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1201 exit(2);
1202 }
1203 for (i=0; i< (unsigned)degree_of_freedom; i++) {
1204 max_value[i]=-10000; /* low value */
1205 }
1206
1207 /* update nodal values from the calculated edge values */
1208 DAG_traverse_list_reset(&Dt);
1209 t=DAG_traverse_prune_list(&Dt);
1210 while(t) {
1211 if (t && t->status==LIVE) {
1212 x[0]=Mx(t->p0,Ms)*g_xscale-g_xoff;
1213 y[0]=My(t->p0,Ms)*g_yscale-g_yoff;
1214 x[1]=Mx(t->p1,Ms)*g_xscale-g_xoff;
1215 y[1]=My(t->p1,Ms)*g_yscale-g_yoff;
1216 x[2]=Mx(t->p2,Ms)*g_xscale-g_xoff;
1217 y[2]=My(t->p2,Ms)*g_yscale-g_yoff;
1218
1219 a[0]=x[1]*y[2]-x[2]*y[1];
1220 a[1]=x[2]*y[0]-x[0]*y[2];
1221 a[2]=x[0]*y[1]-x[1]*y[0];
1222 b[0]=y[1]-y[2];
1223 b[1]=y[2]-y[0];
1224 b[2]=y[0]-y[1];
1225 c[0]=x[2]-x[1];
1226 c[1]=x[0]-x[2];
1227 c[2]=x[1]-x[0];
1228
1229 A[0]=a[0]*b[1]-a[1]*b[0];
1230 A[1]=a[1]*b[2]-a[2]*b[1];
1231 A[2]=a[2]*b[0]-a[0]*b[2];
1232 B[0]=c[0]*b[1]-c[1]*b[0];
1233 B[1]=c[1]*b[2]-c[2]*b[1];
1234 B[2]=c[2]*b[0]-c[0]*b[2];
1235 C[0]=a[0]*c[1]-a[1]*c[0];
1236 C[1]=a[1]*c[2]-a[2]*c[1];
1237 C[2]=a[2]*c[0]-a[0]*c[2];
1238 D[0]=-B[0];
1239 D[1]=-B[1];
1240 D[2]=-B[2];
1241 /* edge lengths */
1242 L[0]=sqrt(ABS(c[2]*c[2]+b[2]*b[2]));
1243 L[1]=sqrt(ABS(c[0]*c[0]+b[0]*b[0]));
1244 L[2]=sqrt(ABS(c[1]*c[1]+b[1]*b[1]));
1245 /* area */
1246 Ar=0.5*(a[1]+a[2]+a[0]);
1247
1248 #ifdef DEBUG
1249 printf("build_hash: adding tria "MIF" X:("MDF","MDF","MDF")Y:("MDF","MDF","MDF") ",t->n,x[0],x[1],x[2],y[0],y[1],y[2]);
1250 #endif
1251 /* normalize edges: make L -> L/(4A^2) */
1252 Ar=1/(4*Ar*Ar);
1253 L[0]*=Ar;
1254 L[1]*=Ar;
1255 L[2]*=Ar;
1256
1257 /* hashes for the edges */
1258 e1=get_edge_number(t->p0,t->p1,H);
1259 e2=get_edge_number(t->p1,t->p2,H);
1260 e3=get_edge_number(t->p2,t->p0,H);
1261
1262 for (i=0;i<(unsigned)degree_of_freedom;i++) {
1263 if (e1!=-1) {
1264 et[0]=v[e1][i];
1265 #ifdef DEBUG
1266 printf("<=="MDF","MDF"==>",et[0],v[e1][i]);
1267 #endif
1268 } else {
1269 et[0]=0.0;
1270 }
1271 if (e2!=-1) {
1272 et[1]=v[e2][i];
1273 #ifdef DEBUG
1274 printf("<=="MDF","MDF"==>",et[1],v[e2][i]);
1275 #endif
1276
1277 } else {
1278 et[1]=0.0;
1279 }
1280 if (e3!=-1) {
1281 et[2]=v[e3][i];
1282 #ifdef DEBUG
1283 printf("<=="MDF","MDF"==>",et[2],v[e3][i]);
1284 #endif
1285
1286 } else {
1287 et[2]=0.0;
1288 }
1289
1290 #ifdef DEBUG
1291 printf("edges %d(%d,%d), %d(%d,%d), %d(%d,%d),t1="MDF" t2="MDF" t3="MDF" ",e1,t->p0,t->p1,e2,t->p1,t->p2,e3,t->p2,t->p0,et[0],et[1],et[2]);
1292 #endif
1293 Mzz(t->p0,3*i,Ms)=0.0;
1294 Mzz(t->p0,3*i+1,Ms)=0.0;
1295 Mzz(t->p1,3*i,Ms)=0.0;
1296 Mzz(t->p1,3*i+1,Ms)=0.0;
1297 Mzz(t->p2,3*i,Ms)=0.0;
1298 Mzz(t->p2,3*i+1,Ms)=0.0;
1299 for (j=0; j<3; j++) {
1300 Mzz(t->p0,3*i,Ms)+=et[j]*L[j]*(A[j]+B[j]*y[0]);
1301 Mzz(t->p0,3*i+1,Ms)+=et[j]*L[j]*(C[j]+D[j]*x[0]);
1302 Mzz(t->p1,3*i,Ms)+=et[j]*L[j]*(A[j]+B[j]*y[1]);
1303 Mzz(t->p1,3*i+1,Ms)+=et[j]*L[j]*(C[j]+D[j]*x[1]);
1304 Mzz(t->p2,3*i,Ms)+=et[j]*L[j]*(A[j]+B[j]*y[2]);
1305 Mzz(t->p2,3*i+1,Ms)+=et[j]*L[j]*(C[j]+D[j]*x[2]);
1306 }
1307
1308 /* z component */
1309 Mzz(t->p0,3*i+2,Ms)=v[renumber[t->p0]][i];
1310 Mzz(t->p1,3*i+2,Ms)=v[renumber[t->p1]][i];
1311 Mzz(t->p2,3*i+2,Ms)=v[renumber[t->p2]][i];
1312
1313 /* find min, max value for plotting */
1314 if ( ABS(Mzz(t->p0,3*i,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p0,3*i,Ms));
1315 if ( ABS(Mzz(t->p0,3*i+1,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p0,3*i+1,Ms));
1316 if ( ABS(Mzz(t->p0,3*i+2,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p0,3*i+2,Ms));
1317 if ( ABS(Mzz(t->p1,3*i,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p1,3*i,Ms));
1318 if ( ABS(Mzz(t->p1,3*i+1,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p1,3*i+1,Ms));
1319 if ( ABS(Mzz(t->p1,3*i+2,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p1,3*i+2,Ms));
1320 if ( ABS(Mzz(t->p2,3*i,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p2,3*i,Ms));
1321 if ( ABS(Mzz(t->p2,3*i+1,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p2,3*i+1,Ms));
1322 if ( ABS(Mzz(t->p2,3*i+2,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p2,3*i+2,Ms));
1323
1324
1325
1326
1327 #ifdef DEBUG
1328 printf("Z:("MDF","MDF","MDF")("MDF","MDF","MDF")("MDF","MDF","MDF")\n",Mzz(t->p0,3*i,Ms),Mzz(t->p1,3*i,Ms),Mzz(t->p2,3*i,Ms),Mzz(t->p0,3*i+1,Ms),Mzz(t->p1,3*i+1,Ms),Mzz(t->p2,3*i+1,Ms),Mzz(t->p0,3*i+2,Ms),Mzz(t->p1,3*i+2,Ms),Mzz(t->p2,3*i+2,Ms));
1329 #endif
1330 }
1331
1332 }
1333 t=DAG_traverse_prune_list(&Dt);
1334 }
1335
1336
1337 /* normalize all potentials by dividing their max values */
1338 for (i=0; i<Ms.count; i++) {
1339 for (j=0; j<degree_of_freedom; j++) {
1340 if ( max_value[j]>TOL) {
1341 Mzz(i,3*j,Ms) /=max_value[j];
1342 Mzz(i,3*j+1,Ms) /=max_value[j];
1343 Mzz(i,3*j+2,Ms) /=max_value[j];
1344 }
1345 }
1346 }
1347 g_maxpot = 1.0;
1348 g_minpot = -1.0;
1349
1350
1351
1352 #ifdef DEBUG
1353 printf("potentials max=%lf min=%lf\n",g_maxpot,g_minpot);
1354 #endif
1355 /**********************************************************/
1356 /* destroy hash table */
1357 htbl_destroy(&H);
1358
1359
1360 for (j = 0; j < NN; j++) {
1361 free(v[j]);
1362 }
1363 SPM_destroy(&S);
1364 SPM_destroy(&T);
1365 free(v);
1366 free(renumber);
1367 free(max_value);
1368 free(ee);
1369 return(1);
1370
1371 }
1372