/* pdnmesh - a 2D finite element solver Copyright (C) 2001-2005 Sarod Yatawatta This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA $Id: wave.c,v 1.7 2005/02/22 03:26:02 sarod Exp $ */ /* formulation of inhomogeneous wave equation */ #include #include #include #include #include "types.h" /* cutoff frequency and propagation constant used in waveguide problems */ MY_DOUBLE global_wave_k0, global_wave_beta; /* if 1 use full matrix solver, else use half the matrix, for symmetric positive definite eigen equation */ int use_full_eigensolver; typedef struct hash_edge_{ MY_INT p0; /* point 0*/ MY_INT p1; /* point 1*/ MY_INT n; /* unique number */ } hash_edge; static unsigned int hashfn(const void *data) { unsigned int k,l; hash_edge *e=(hash_edge*)data; k=(int)e->p0; l=(int)e->p1; #ifdef DEBUG printf("hash : key value=%d,%d\n",k,l); #endif k+=~(k<<15); l^=(k>>10); k+=(l<<3); l^=(k>>6); k+=~(l<<11); l^=(k>>16); k+=l; #ifdef DEBUG printf("hash : value=%d\n",k); #endif return(k); } /* 64 bit hash function */ /*long longhash1(long key) { key += ~(key << 32); key ^= (key >>> 22); key += ~(key << 13); key ^= (key >>> 8); key += (key << 3); key ^= (key >>> 15); key += ~(key << 27); key ^= (key >>> 31); return key; } */ /* generic list functions */ /* comparison function */ /* compare two edges */ int match(const void *key1, const void *key2) { hash_edge *e1,*e2; e1=(hash_edge*)key1; e2=(hash_edge*)key2; #ifdef DEBUG printf("comparing (%d,%d) with (%d,%d)\n",e1->p0,e1->p1,e2->p0,e2->p1); #endif if((e1->p0==e2->p0) && (e1->p1==e2->p1)); return(0); /* else */ return(-1); } /* freeing function */ void destroy(void *data) { hash_edge *t=(hash_edge *)data; free(t); } /* if edge (p0-p1) is on a dirichlet boundary, return 1. * else return 0 */ static MY_INT is_this_edge_on_a_dirichlet_boundary(p0,p1) { /* find the ratio using x=(x1+lambda. x2)/(1+lambda) */ /* if x is on (x1-x2), lambda >0 , equal to y */ /* try to do robust computations */ int i; double xp0,xp1,xn0,xn1,yp0,yp1,yn0,yn1; xp0=Mx(p0,M); yp0=My(p0,M); xp1=Mx(p1,M); yp1=My(p1,M); #ifdef DEBUG printf("is_this_edge_on_a_boundary:consider points(%d,%d)\n",p0,p1); #endif for (i=0; i< nedges; i++) { xn0=Mx(edge_array[i].p1,M); yn0=My(edge_array[i].p1,M); xn1=Mx(edge_array[i].p2,M); yn1=My(edge_array[i].p2,M); #ifdef DEBUG printf("is_this_edge_on_a_boundary: edge %d\n",i); #endif if ( IS_ZERO((xn0-xp0)*(yp0-yn1)-(yn0-yp0)*(xp0-xn1)) /* ratio is equal magnitude, check sign */ && ( ((yn0-yp0)*(yp0-yn1)>0) /* if true, equal and +ve lambda. stop */ || (IS_ZERO((yn0-yp0)*(yp0-yn1)) && (xn0-xp0)*(xp0-xn1) >= 0))){ /* if we are here, point p0 is on line */ /* do the same check for point p1 */ #ifdef DEBUG printf("is_this_edge_on_a_boundary: point %d is on line %d\n",p0,i); #endif if ( IS_ZERO((xn0-xp1)*(yp1-yn1)-(yn0-yp1)*(xp1-xn1)) /* ratio is equal magnitude, check sign */ && ( ((yn0-yp1)*(yp1-yn1)>0) /* if true, equal and +ve lambda. stop */ || (IS_ZERO((yn0-yp1)*(yp1-yn1)) && (xn0-xp1)*(xp1-xn1) >= 0))){ if ( edge_array[i].type == DR ) { return(1); } else { return(0); } } } } return(0); /* not found */ } static MY_INT get_edge_number(MY_INT p0,MY_INT p1, htbl H) { hash_edge *ee,*gg; if ( (ee=(hash_edge*)malloc(sizeof(hash_edge)))==0) { fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__); exit(1); } if ( p0 > p1 ) { ee->p0=p1; ee->p1=p0; } else { ee->p0=p0; ee->p1=p1; } gg=(hash_edge*)htbl_lookup(&H, (void *)ee); free(ee); if (gg) { return(gg->n); } else { printf("get_edge_number: something SERIOUSLY wrong! (%d,%d)\n",p0,p1); return(-1); } } /* solve for cutoff frequency given propagation constant */ 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) { /* S, T - Ntotal by Ntotal matrices Nedges - number of edges, the rest points renumber - map points to matrix location edges are mapper using hash table H */ triangle *t; MY_INT i,j,m,n; MY_DOUBLE x[3],y[3],a[3],b[3],c[3],x_bar,y_bar,x_var,y_var; MY_DOUBLE A[3],B[3],C[3],D[3],L[3],Ar,k; 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]; MY_DOUBLE I1[3][3],I2[3][3],I3[3][3],I4[3][3],I5[3][3]; MY_DOUBLE epsilon,sumI,mu; MY_INT e[3],mult[3]; DAG_traverse_list_reset(&dt); t=DAG_traverse_prune_list(&dt); while(t) { if (t && t->status==LIVE) { epsilon= bound.poly_array[t->boundary].epsilon; mu= bound.poly_array[t->boundary].mu; /* get three edges */ /* p0 /\ / \ edge 1 / \ edge 3 / \ /________\ p1 edge 2 p2 */ /* get coords */ /* convert to global coords */ x[0]=Mx(t->p0,M)*g_xscale-g_xoff; y[0]=My(t->p0,M)*g_yscale-g_yoff; x[1]=Mx(t->p1,M)*g_xscale-g_xoff; y[1]=My(t->p1,M)*g_yscale-g_yoff; x[2]=Mx(t->p2,M)*g_xscale-g_xoff; y[2]=My(t->p2,M)*g_yscale-g_yoff; x_bar=(x[1]+x[2]+x[0])*ONE_OVER_THREE; y_bar=(y[1]+y[2]+y[0])*ONE_OVER_THREE; x_var=(x[1]*x[1]+x[2]*x[2]+x[0]*x[0]); y_var=(y[1]*y[1]+y[2]*y[2]+y[0]*y[0]); a[0]=x[1]*y[2]-x[2]*y[1]; a[1]=x[2]*y[0]-x[0]*y[2]; a[2]=x[0]*y[1]-x[1]*y[0]; b[0]=y[1]-y[2]; b[1]=y[2]-y[0]; b[2]=y[0]-y[1]; c[0]=x[2]-x[1]; c[1]=x[0]-x[2]; c[2]=x[1]-x[0]; A[0]=a[0]*b[1]-a[1]*b[0]; A[1]=a[1]*b[2]-a[2]*b[1]; A[2]=a[2]*b[0]-a[0]*b[2]; B[0]=c[0]*b[1]-c[1]*b[0]; B[1]=c[1]*b[2]-c[2]*b[1]; B[2]=c[2]*b[0]-c[0]*b[2]; C[0]=a[0]*c[1]-a[1]*c[0]; C[1]=a[1]*c[2]-a[2]*c[1]; C[2]=a[2]*c[0]-a[0]*c[2]; D[0]=-B[0]; D[1]=-B[1]; D[2]=-B[2]; /* edge lengths */ L[0]=sqrt(ABS(c[2]*c[2]+b[2]*b[2])); L[1]=sqrt(ABS(c[0]*c[0]+b[0]*b[0])); L[2]=sqrt(ABS(c[1]*c[1]+b[1]*b[1])); /* area */ Ar=0.5*(a[1]+a[2]+a[0]); for (m=0; m<3; m++) { for (n=0; n<3; n++) { I1[m][n]=A[m]*A[n]+C[m]*C[n]; I2[m][n]=(C[m]*D[n]+C[n]*D[m])*x_bar; I3[m][n]=(A[m]*B[n]+A[n]*B[m])*y_bar; I4[m][n]=ONE_OVER_TWELVE*B[m]*B[n]*(y_var+9*y_bar*y_bar); I5[m][n]=ONE_OVER_TWELVE*D[m]*D[n]*(x_var+9*x_bar*x_bar); } } k=1/(16*Ar*Ar*Ar); for (m=0; m<3; m++) { for (n=0; n<3; n++) { sumI=I1[m][n]+I2[m][n]+I3[m][n]+I4[m][n]+I5[m][n]; S_tt[m][n]=k*L[m]*L[n]*(D[m]*D[n]+sumI)/(mu); T_tt[m][n]=k*epsilon*L[m]*L[n]*sumI; } } k=beta*beta/(mu*8*Ar*Ar); for (m=0; m<3; m++) { for (n=0; n<3; n++) { S_tz[m][n]=k*L[m]*(b[n]*(A[m]+B[m]*y_bar)+c[n]*(C[m]+D[m]*x_bar)); S_zt[n][m]=k*L[m]*(b[n]*(A[m]+B[m]*y_bar)+c[n]*(C[m]+D[m]*x_bar)); } } k=beta*beta/(mu*4*Ar); for (m=0; m<3; m++) { for (n=0; n<3; n++) { S_zz[m][n]=k*(b[m]*b[n]+c[m]*c[n]); } } k=beta*beta*epsilon*Ar*ONE_OVER_TWELVE; for (m=0; m<3; m++) { for (n=0; n<3; n++) { T_zz[m][n]=k; } } T_zz[0][0]=T_zz[1][1]=T_zz[2][2]=2*k; /* 1/6 */ /* hashes for the edges */ e[0]=get_edge_number(t->p0,t->p1,H); e[1]=get_edge_number(t->p1,t->p2,H); e[2]=get_edge_number(t->p2,t->p0,H); /* || Stt Stz || || Ttt 0|| || Szt Szz || || 0 Tzz|| S T Sx = lambda Tx */ /* consistencay with sign change */ mult[0]=mult[1]=mult[2]=1; if (t->p0>t->p1) { mult[0]=-1; } if (t->p1>t->p2) { mult[1]=-1; } if (t->p2>t->p0) { mult[2]=-1; } /* signs changes Stt, and Ttt - both rows and columns Stz rows Szt columns */ for (i=0; i<3; i++) { for (j=0; j<3; j++) { if ( mult[i]*mult[j] < 0 ) { S_tt[i][j]=-S_tt[i][j]; T_tt[i][j]=-T_tt[i][j]; } } if ( mult[i] < 0 ) { for (j=0; j<3; j++) { S_tz[i][j]=-S_tz[i][j]; S_zt[j][i]=-S_zt[j][i]; } } } /* assemble global matrix */ /* if any of the edges are Dirichlet, we do not include them in the global matrix */ #ifdef DEBUG 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]); for (i=0;i<3;i++) { printf("||"); 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]); printf("|| ||"); 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); printf("||\n"); } for (i=0;i<3;i++) { printf("||"); 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]); printf("|| ||"); 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]); printf("||\n"); } #endif /* assembling global matrix */ for (i=0; i<3; i++) { if ( e[i] != -1 ) { for (j=0;j<3;j++) { if (e[j] != -1 ) { /*T[e[i]][e[j]]+=(T_tt[i][j]); */ SPM_add(T,e[i],e[j],T_tt[i][j]); /*S[e[i]][e[j]]+=(S_tt[i][j]); */ SPM_add(S,e[i],e[j],S_tt[i][j]); } } } } i=renumber[t->p0]; if (i>=Nedges) { /* fixed points have i == 0 */ for (j=0; j<3; j++) { if ( e[j]!=-1) { SPM_add(S,i,e[j],S_zt[0][j]); } if ( e[j]!=-1) { SPM_add(S,e[j],i,S_tz[j][0]); } } } i=renumber[t->p1]; if (i>=Nedges) { /* fixed points have i == 0 */ for (j=0; j<3; j++) { if ( e[j]!=-1) { SPM_add(S,i,e[j],S_zt[1][j]); } if ( e[j]!=-1) { SPM_add(S,e[j],i,S_tz[j][1]); } } } i=renumber[t->p2]; if (i>=Nedges) { /* fixed points have i == 0 */ for (j=0; j<3; j++) { if ( e[j]!=-1) { SPM_add(S,i,e[j],S_zt[2][j]); } if ( e[j]!=-1) { SPM_add(S,e[j],i,S_tz[j][2]); } } } /* user e array to store point numbers temporarily */ e[0]=renumber[t->p0];e[1]=renumber[t->p1];e[2]=renumber[t->p2]; for (i=0; i<3; i++) { if ( e[i] >=Nedges ) { for (j=0;j<3; j++) { if ( e[j] >=Nedges ) { /*S[e[i]][e[j]]+=S_zz[i][j]; */ SPM_add(S,e[i],e[j],S_zz[i][j]); /*T[e[i]][e[j]]+=T_zz[i][j]; */ SPM_add(T,e[i],e[j],T_zz[i][j]); } } } } } t=DAG_traverse_prune_list(&dt); } } /* solve for propagation constant given frequency */ 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) { triangle *t; MY_INT i,j,m,n; MY_DOUBLE x[3],y[3],a[3],b[3],c[3],x_bar,y_bar,x_var,y_var; MY_DOUBLE A[3],B[3],C[3],D[3],L[3],Ar,k; 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]; MY_DOUBLE I1[3][3],I2[3][3],I3[3][3],I4[3][3],I5[3][3]; MY_DOUBLE epsilon,sumI,mu; MY_INT e[3],mult[3]; DAG_traverse_list_reset(&dt); t=DAG_traverse_prune_list(&dt); while(t) { if (t && t->status==LIVE) { epsilon = bound.poly_array[t->boundary].epsilon; mu= bound.poly_array[t->boundary].mu; /* get three edges */ /* p0 /\ / \ edge 1 / \ edge 3 / \ /________\ p1 edge 2 p2 */ /* get coords */ /* convert to global coords */ x[0]=Mx(t->p0,M)*g_xscale-g_xoff; y[0]=My(t->p0,M)*g_yscale-g_yoff; x[1]=Mx(t->p1,M)*g_xscale-g_xoff; y[1]=My(t->p1,M)*g_yscale-g_yoff; x[2]=Mx(t->p2,M)*g_xscale-g_xoff; y[2]=My(t->p2,M)*g_yscale-g_yoff; x_bar=(x[1]+x[2]+x[0])*ONE_OVER_THREE; y_bar=(y[1]+y[2]+y[0])*ONE_OVER_THREE; x_var=(x[1]*x[1]+x[2]*x[2]+x[0]*x[0]); y_var=(y[1]*y[1]+y[2]*y[2]+y[0]*y[0]); a[0]=x[1]*y[2]-x[2]*y[1]; a[1]=x[2]*y[0]-x[0]*y[2]; a[2]=x[0]*y[1]-x[1]*y[0]; b[0]=y[1]-y[2]; b[1]=y[2]-y[0]; b[2]=y[0]-y[1]; c[0]=x[2]-x[1]; c[1]=x[0]-x[2]; c[2]=x[1]-x[0]; A[0]=a[0]*b[1]-a[1]*b[0]; A[1]=a[1]*b[2]-a[2]*b[1]; A[2]=a[2]*b[0]-a[0]*b[2]; B[0]=c[0]*b[1]-c[1]*b[0]; B[1]=c[1]*b[2]-c[2]*b[1]; B[2]=c[2]*b[0]-c[0]*b[2]; C[0]=a[0]*c[1]-a[1]*c[0]; C[1]=a[1]*c[2]-a[2]*c[1]; C[2]=a[2]*c[0]-a[0]*c[2]; D[0]=-B[0]; D[1]=-B[1]; D[2]=-B[2]; /* edge lengths */ L[0]=sqrt(ABS(c[2]*c[2]+b[2]*b[2])); L[1]=sqrt(ABS(c[0]*c[0]+b[0]*b[0])); L[2]=sqrt(ABS(c[1]*c[1]+b[1]*b[1])); /* area */ Ar=0.5*(a[1]+a[2]+a[0]); for (m=0; m<3; m++) { for (n=0; n<3; n++) { I1[m][n]=A[m]*A[n]+C[m]*C[n]; I2[m][n]=(C[m]*D[n]+C[n]*D[m])*x_bar; I3[m][n]=(A[m]*B[n]+A[n]*B[m])*y_bar; I4[m][n]=ONE_OVER_TWELVE*B[m]*B[n]*(y_var+9*y_bar*y_bar); I5[m][n]=ONE_OVER_TWELVE*D[m]*D[n]*(x_var+9*x_bar*x_bar); } } k=1/(16*Ar*Ar*Ar); for (m=0; m<3; m++) { for (n=0; n<3; n++) { sumI=I1[m][n]+I2[m][n]+I3[m][n]+I4[m][n]+I5[m][n]; S_tt[m][n]=k*L[m]*L[n]*(D[m]*D[n]-k0*k0*epsilon*sumI)/(mu); T_tt[m][n]=k*epsilon*L[m]*L[n]*sumI; } } k=1/(mu*8*Ar*Ar); for (m=0; m<3; m++) { for (n=0; n<3; n++) { T_tz[m][n]=k*L[m]*(b[n]*(A[m]+B[m]*y_bar)+c[n]*(C[m]+D[m]*x_bar)); T_zt[n][m]=k*L[m]*(b[n]*(A[m]+B[m]*y_bar)+c[n]*(C[m]+D[m]*x_bar)); } } k=1/(mu*4*Ar); for (m=0; m<3; m++) { for (n=0; n<3; n++) { S_zz[m][n]=k*(b[m]*b[n]+c[m]*c[n]); } } k=k0*k0*epsilon*Ar*ONE_OVER_TWELVE; for (m=0; m<3; m++) { for (n=0; n<3; n++) { T_zz[m][n]=k; } } T_zz[0][0]=T_zz[1][1]=T_zz[2][2]=2*k; /* 1/6 */ for (m=0; m<3; m++) { for (n=0; n<3; n++) { T_zz[m][n]+=S_zz[m][n]; } } /* hashes for the edges */ e[0]=get_edge_number(t->p0,t->p1,H); e[1]=get_edge_number(t->p1,t->p2,H); e[2]=get_edge_number(t->p2,t->p0,H); /* consistencay with sign change */ mult[0]=mult[1]=mult[2]=1; if (t->p0>t->p1) { mult[0]=-1; } if (t->p1>t->p2) { mult[1]=-1; } if (t->p2>t->p0) { mult[2]=-1; } for (i=0; i<3; i++) { for (j=0; j<3; j++) { if ( mult[i]*mult[j] < 0 ) { S_tt[i][j]=-S_tt[i][j]; T_tt[i][j]=-T_tt[i][j]; } } if ( mult[i] < 0 ) { for (j=0; j<3; j++) { T_tz[i][j]=-T_tz[i][j]; T_zt[j][i]=-T_zt[j][i]; } } } /* assemble global matrix */ /* if any of the edges are Dirichlet, we do not include them in the global matrix */ #ifdef DEBUG 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]); for (i=0;i<3;i++) { printf("||"); 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); printf("|| ||"); 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]); printf("||\n"); } for (i=0;i<3;i++) { printf("||"); printf(MDF" "MDF" "MDF" "MDF" "MDF" "MDF" ",0.0,0.0,0.0,0.0,0.0,0.0); printf("|| ||"); 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]); printf("||\n"); } #endif /* assembling global matrix */ for (i=0; i<3; i++) { if ( e[i] != -1 ) { for (j=0;j<3;j++) { if (e[j] != -1 ) { /*S[e[i]][e[j]]+=(S_tt[i][j]); */ SPM_add(S,e[i],e[j],S_tt[i][j]); /*T[e[i]][e[j]]+=(T_tt[i][j]); */ SPM_add(T,e[i],e[j],T_tt[i][j]); } } } } i=renumber[t->p0]; if (i>=Nedges) { /* fixed points have i == 0 */ for (j=0; j<3; j++) { if ( e[j]!=-1) { SPM_add(T,i,e[j],T_tz[0][j]); } if ( e[j]!=-1) { SPM_add(T,e[j],i,T_zt[j][0]); } } } i=renumber[t->p1]; if (i>=Nedges) { /* fixed points have i == 0 */ for (j=0; j<3; j++) { if ( e[j]!=-1) { SPM_add(T,i,e[j],T_tz[1][j]); } if ( e[j]!=-1) { SPM_add(T,e[j],i,T_zt[j][1]); } } } i=renumber[t->p2]; if (i>=Nedges) { /* fixed points have i == 0 */ for (j=0; j<3; j++) { if ( e[j]!=-1) { SPM_add(T,i,e[j],T_tz[2][j]); } if ( e[j]!=-1) { SPM_add(T,e[j],i,T_zt[j][2]); } } } /* user e array to store point numbers temporarily */ e[0]=renumber[t->p0];e[1]=renumber[t->p1];e[2]=renumber[t->p2]; for (i=0; i<3; i++) { if ( e[i] >=Nedges ) { for (j=0;j<3; j++) { if ( e[j] >=Nedges ) { SPM_add(T,e[i],e[j],T_zz[i][j]); } } } } } t=DAG_traverse_prune_list(&dt); } } /* solve for cutoff frequency */ static void build_equations_from_hash_table(htbl H, SPMat *S, SPMat *T, MY_INT *renumber, MY_INT Ntotal, MY_INT Nedges) { /* S, T - Ntotal by Ntotal matrices Nedges - number of edges, the rest points renumber - map points to matrix location edges are mapper using hash table H */ triangle *t; MY_INT i,j,m,n; MY_DOUBLE x[3],y[3],a[3],b[3],c[3],x_bar,y_bar,x_var,y_var; MY_DOUBLE A[3],B[3],C[3],D[3],L[3],Ar,k; MY_DOUBLE S_tt[3][3],S_zz[3][3],T_tt[3][3],T_zz[3][3]; MY_DOUBLE I1[3][3],I2[3][3],I3[3][3],I4[3][3],I5[3][3]; MY_DOUBLE epsilon,sumI,mu; MY_INT e[3],mult[3]; DAG_traverse_list_reset(&dt); t=DAG_traverse_prune_list(&dt); while(t) { if (t && t->status==LIVE) { epsilon= bound.poly_array[t->boundary].epsilon; mu= bound.poly_array[t->boundary].mu; /* get three edges */ /* p0 /\ / \ edge 1 / \ edge 3 / \ /________\ p1 edge 2 p2 */ /* get coords */ /* convert to global coords */ x[0]=Mx(t->p0,M)*g_xscale-g_xoff; y[0]=My(t->p0,M)*g_yscale-g_yoff; x[1]=Mx(t->p1,M)*g_xscale-g_xoff; y[1]=My(t->p1,M)*g_yscale-g_yoff; x[2]=Mx(t->p2,M)*g_xscale-g_xoff; y[2]=My(t->p2,M)*g_yscale-g_yoff; x_bar=(x[1]+x[2]+x[0])*ONE_OVER_THREE; y_bar=(y[1]+y[2]+y[0])*ONE_OVER_THREE; x_var=(x[1]*x[1]+x[2]*x[2]+x[0]*x[0]); y_var=(y[1]*y[1]+y[2]*y[2]+y[0]*y[0]); a[0]=x[1]*y[2]-x[2]*y[1]; a[1]=x[2]*y[0]-x[0]*y[2]; a[2]=x[0]*y[1]-x[1]*y[0]; b[0]=y[1]-y[2]; b[1]=y[2]-y[0]; b[2]=y[0]-y[1]; c[0]=x[2]-x[1]; c[1]=x[0]-x[2]; c[2]=x[1]-x[0]; A[0]=a[0]*b[1]-a[1]*b[0]; A[1]=a[1]*b[2]-a[2]*b[1]; A[2]=a[2]*b[0]-a[0]*b[2]; B[0]=c[0]*b[1]-c[1]*b[0]; B[1]=c[1]*b[2]-c[2]*b[1]; B[2]=c[2]*b[0]-c[0]*b[2]; C[0]=a[0]*c[1]-a[1]*c[0]; C[1]=a[1]*c[2]-a[2]*c[1]; C[2]=a[2]*c[0]-a[0]*c[2]; D[0]=-B[0]; D[1]=-B[1]; D[2]=-B[2]; /* edge lengths */ L[0]=sqrt(ABS(c[2]*c[2]+b[2]*b[2])); L[1]=sqrt(ABS(c[0]*c[0]+b[0]*b[0])); L[2]=sqrt(ABS(c[1]*c[1]+b[1]*b[1])); /* area */ Ar=0.5*(a[1]+a[2]+a[0]); for (m=0; m<3; m++) { for (n=0; n<3; n++) { I1[m][n]=A[m]*A[n]+C[m]*C[n]; I2[m][n]=(C[m]*D[n]+C[n]*D[m])*x_bar; I3[m][n]=(A[m]*B[n]+A[n]*B[m])*y_bar; I4[m][n]=ONE_OVER_TWELVE*B[m]*B[n]*(y_var+9*y_bar*y_bar); I5[m][n]=ONE_OVER_TWELVE*D[m]*D[n]*(x_var+9*x_bar*x_bar); } } k=1/(16*Ar*Ar*Ar); for (m=0; m<3; m++) { for (n=0; n<3; n++) { sumI=I1[m][n]+I2[m][n]+I3[m][n]+I4[m][n]+I5[m][n]; S_tt[m][n]=4*k*L[m]*L[n]*(D[m]*D[n])/(mu); T_tt[m][n]=k*epsilon*L[m]*L[n]*sumI; } } k=1/(mu*4*Ar); for (m=0; m<3; m++) { for (n=0; n<3; n++) { S_zz[m][n]=k*(b[m]*b[n]+c[m]*c[n]); } } k=epsilon*Ar*ONE_OVER_TWELVE; for (m=0; m<3; m++) { for (n=0; n<3; n++) { T_zz[m][n]=k; } } T_zz[0][0]=T_zz[1][1]=T_zz[2][2]=2*k; /* 1/6 */ /* hashes for the edges */ e[0]=get_edge_number(t->p0,t->p1,H); e[1]=get_edge_number(t->p1,t->p2,H); e[2]=get_edge_number(t->p2,t->p0,H); /* || Stt 0 || || Ttt 0|| || 0 Szz || || 0 Tzz|| S T Sx = lambda Tx */ /* consistencay with sign change */ mult[0]=mult[1]=mult[2]=1; if (t->p0>t->p1) { mult[0]=-1; } if (t->p1>t->p2) { mult[1]=-1; } if (t->p2>t->p0) { mult[2]=-1; } /* signs changes Stt, and Ttt - both rows and columns Stz rows Szt columns */ for (i=0; i<3; i++) { for (j=0; j<3; j++) { if ( mult[i]*mult[j] < 0 ) { S_tt[i][j]=-S_tt[i][j]; T_tt[i][j]=-T_tt[i][j]; } } } /* assemble global matrix */ /* if any of the edges are Dirichlet, we do not include them in the global matrix */ #ifdef DEBUG 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]); for (i=0;i<3;i++) { printf("||"); 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); printf("|| ||"); 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); printf("||\n"); } for (i=0;i<3;i++) { printf("||"); 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]); printf("|| ||"); 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]); printf("||\n"); } #endif /* assembling global matrix */ for (i=0; i<3; i++) { if ( e[i] != -1 ) { for (j=0;j<3;j++) { if (e[j] != -1 ) { /*T[e[i]][e[j]]+=(T_tt[i][j]); */ SPM_add(T,e[i],e[j],T_tt[i][j]); /*S[e[i]][e[j]]+=(S_tt[i][j]); */ SPM_add(S,e[i],e[j],S_tt[i][j]); } } } } /* user e array to store point numbers temporarily */ e[0]=renumber[t->p0];e[1]=renumber[t->p1];e[2]=renumber[t->p2]; for (i=0; i<3; i++) { if ( e[i] >=Nedges ) { for (j=0;j<3; j++) { if ( e[j] >=Nedges ) { /*S[e[i]][e[j]]+=S_zz[i][j]; */ SPM_add(S,e[i],e[j],S_zz[i][j]); /*T[e[i]][e[j]]+=T_zz[i][j]; */ SPM_add(T,e[i],e[j],T_zz[i][j]); } } } } } t=DAG_traverse_prune_list(&dt); } } int buildup_hash_table_and_equation(Mesh Ms, DAG_tree Dt) { triangle *t; htbl H; hash_edge *ee,*eep; #ifdef DEBUG hash_edge *gg; #endif MY_INT j,Ni,temp; unsigned int i; SPMat S, T; MY_DOUBLE **v; MY_INT *renumber, *re_renumber; MY_DOUBLE *max_value; MY_INT NN; MY_DOUBLE x[3],y[3],a[3],b[3],c[3],et[2],Ar; MY_DOUBLE A[3],B[3],C[3],D[3],L[3]; MY_INT e1,e2,e3; /* make degree of freedom=1 */ degree_of_freedom=requested_degree_of_freedom; /* find size needed for hash table */ i=Ms.count; #ifdef DEBUG printf("buildup_hash: creating table size %d\n",i); #endif htbl_init(&H, i, sizeof(hash_edge), &hashfn, &match, &destroy); /* create renumbering array for points */ /* keep point numbers also in edge number array. we first go through the edges and number them. next we run another triangle loop to number the points. We alwasy remember last edge number. So everthyng after that is a point. */ i=0; if ( (ee=(hash_edge*)malloc(sizeof(hash_edge)))==0) { fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__); exit(1); } /* insert data */ DAG_traverse_list_reset(&Dt); t=DAG_traverse_prune_list(&Dt); while (t) { if ( t && t->status==LIVE ) { ee->p0=t->p0; ee->p1=t->p1; /* always p0 <= p1 */ if ( ee->p0 > ee->p1 ) { temp=ee->p0; ee->p0=ee->p1; ee->p1=temp; } #ifdef DEBUG printf("**build: tring to insert (%d,%d) into table\n",ee->p0,ee->p1); #endif if ( (!htbl_insert(&H, (const void *)ee)) ) { #ifdef DEBUG printf("sucessfully inserted (%d,%d) into table\n",ee->p0,ee->p1); #endif /* check whether this edge is on a Dirichlet boundary */ eep=(hash_edge*)htbl_lookup(&H, (void *)ee); if ( !is_this_edge_on_a_dirichlet_boundary(ee->p0,ee->p1) ) { /* give a unique contiguous number */ eep->n=i++; } else { #ifdef DEBUG printf("buildup_hash: edge on a Dirichlet boundary\n"); #endif eep->n=-1; } } else { #ifdef DEBUG printf("failed inserting (%d,%d) into table ",ee->p0,ee->p1); /* lookup */ if ( (gg=(hash_edge*)htbl_lookup(&H, (void *)ee)) ) { printf("it exists as (%d,%d)->%d.\n",gg->p0,gg->p1,gg->n); } #endif } ee->p0=t->p1; ee->p1=t->p2; /* always p0 <= p1 */ if ( ee->p0 > ee->p1 ) { temp=ee->p0; ee->p0=ee->p1; ee->p1=temp; } #ifdef DEBUG printf("**build: tring to insert (%d,%d) into table\n",ee->p0,ee->p1); #endif if ( (!htbl_insert(&H, (const void *)ee)) ) { #ifdef DEBUG printf("sucessfully inserted (%d,%d) into table\n",ee->p0,ee->p1); #endif eep=(hash_edge*)htbl_lookup(&H, (void *)ee); if ( !is_this_edge_on_a_dirichlet_boundary(ee->p0,ee->p1) ) { /* give a unique contiguous number */ eep->n=i++; } else { #ifdef DEBUG printf("buildup_hash: edge on a Dirichlet boundary\n"); #endif eep->n=-1; } } else { #ifdef DEBUG printf("failed inserting (%d,%d) into table ",ee->p0,ee->p1); /* lookup */ if ( (gg=(hash_edge*)htbl_lookup(&H, (void *)ee)) ) { printf("it exists as (%d,%d)->%d.\n",gg->p0,gg->p1,gg->n); } #endif } ee->p0=t->p2; ee->p1=t->p0; /* always p0 <= p1 */ if ( ee->p0 > ee->p1 ) { temp=ee->p0; ee->p0=ee->p1; ee->p1=temp; } #ifdef DEBUG printf("**build: tring to insert (%d,%d) into table\n",ee->p0,ee->p1); #endif if ( (!htbl_insert(&H, (const void *)ee)) ) { #ifdef DEBUG printf("sucessfully inserted (%d,%d) into table\n",ee->p0,ee->p1); #endif eep=(hash_edge*)htbl_lookup(&H, (void *)ee); if ( !is_this_edge_on_a_dirichlet_boundary(ee->p0,ee->p1) ) { /* give a unique contiguous number */ eep->n=i++; } else { #ifdef DEBUG printf("buildup_hash: edge on a Dirichlet boundary\n"); #endif eep->n=-1; } } else { #ifdef DEBUG printf("failed inserting (%d,%d) into table ",ee->p0,ee->p1); /* lookup */ if ( (gg=(hash_edge*)htbl_lookup(&H, (void *)ee)) ) { printf("it exists as (%d,%d)->%d.\n",gg->p0,gg->p1,gg->n); } #endif } } t=DAG_traverse_prune_list(&Dt); } #ifdef DEBUG printf("buildup_hash %d of %d edges, hash table filled %lf\n",i,H.positions,H.size/(double)H.positions); #endif /* remember no of edges */ Ni=i; /******************/ /* setup point to matrix mapping array */ if ((renumber = (MY_INT *)calloc(Ms.count,sizeof(MY_INT )))==0){ fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__); exit(1); } if((re_renumber = (MY_INT *)calloc(Ms.count,sizeof(MY_INT)))==0){ fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__); exit(1); } /* need to solve for only points referred by triangles */ /* SO -- */ DAG_traverse_list_reset(&Dt); t=DAG_traverse_prune_list(&Dt); while ( t ) { if ( t && t->status==LIVE ) { /* mark the points in this triangle */ /* use re_renumber[] array for this */ re_renumber[t->p0] = 1; re_renumber[t->p1] = 1; re_renumber[t->p2] = 1; } t=DAG_traverse_prune_list(&Dt); } /* now give numbers to variable points*/ /* start numbering from last edge */ j = Ni; for ( i = 0; i < Ms.count; i++ ) if ( Mval(i,Ms) == VAR && re_renumber[i] == 1) { renumber[i] = j++; #ifdef DEBUG printf("buildup_hash: renumbering variable point %d as %d\n",i,renumber[i]); #endif } free(re_renumber); /******************/ /* update total size of matrix */ NN=j; /* Ni edges and NN-Ni points */ /* set up matrices */ SPM_init(&S,NN); SPM_init(&T,NN); /* depending on equation type, we formulate a different problem */ if (solve_equation == HELMHOLTZ_INHOMO ) { build_equations_from_hash_table(H, &S, &T, renumber, NN, Ni); } else if (solve_equation == HELMHOLTZ_FREQ) { build_equations_from_hash_table_given_freq(H, &S, &T, renumber, NN, Ni, global_wave_k0); } else if (solve_equation == HELMHOLTZ_BETA) { build_equations_from_hash_table_given_beta(H, &S, &T, renumber, NN, Ni, global_wave_beta); } #ifdef DEBUG printf("====================\n"); printf("S=\n"); for (i=0;iz=(MY_DOUBLE *)realloc((void*)(Ms.Narray[i])->z,(size_t)3*degree_of_freedom*sizeof(MY_DOUBLE)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(2); } } /* array to remember max values of potentials */ if ((max_value=(MY_DOUBLE *)malloc((size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(2); } for (i=0; i< (unsigned)degree_of_freedom; i++) { max_value[i]=-10000; /* low value */ } /* update nodal values from the calculated edge values */ DAG_traverse_list_reset(&Dt); t=DAG_traverse_prune_list(&Dt); while(t) { if (t && t->status==LIVE) { x[0]=Mx(t->p0,Ms)*g_xscale-g_xoff; y[0]=My(t->p0,Ms)*g_yscale-g_yoff; x[1]=Mx(t->p1,Ms)*g_xscale-g_xoff; y[1]=My(t->p1,Ms)*g_yscale-g_yoff; x[2]=Mx(t->p2,Ms)*g_xscale-g_xoff; y[2]=My(t->p2,Ms)*g_yscale-g_yoff; a[0]=x[1]*y[2]-x[2]*y[1]; a[1]=x[2]*y[0]-x[0]*y[2]; a[2]=x[0]*y[1]-x[1]*y[0]; b[0]=y[1]-y[2]; b[1]=y[2]-y[0]; b[2]=y[0]-y[1]; c[0]=x[2]-x[1]; c[1]=x[0]-x[2]; c[2]=x[1]-x[0]; A[0]=a[0]*b[1]-a[1]*b[0]; A[1]=a[1]*b[2]-a[2]*b[1]; A[2]=a[2]*b[0]-a[0]*b[2]; B[0]=c[0]*b[1]-c[1]*b[0]; B[1]=c[1]*b[2]-c[2]*b[1]; B[2]=c[2]*b[0]-c[0]*b[2]; C[0]=a[0]*c[1]-a[1]*c[0]; C[1]=a[1]*c[2]-a[2]*c[1]; C[2]=a[2]*c[0]-a[0]*c[2]; D[0]=-B[0]; D[1]=-B[1]; D[2]=-B[2]; /* edge lengths */ L[0]=sqrt(ABS(c[2]*c[2]+b[2]*b[2])); L[1]=sqrt(ABS(c[0]*c[0]+b[0]*b[0])); L[2]=sqrt(ABS(c[1]*c[1]+b[1]*b[1])); /* area */ Ar=0.5*(a[1]+a[2]+a[0]); #ifdef DEBUG 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]); #endif /* normalize edges: make L -> L/(4A^2) */ Ar=1/(4*Ar*Ar); L[0]*=Ar; L[1]*=Ar; L[2]*=Ar; /* hashes for the edges */ e1=get_edge_number(t->p0,t->p1,H); e2=get_edge_number(t->p1,t->p2,H); e3=get_edge_number(t->p2,t->p0,H); for (i=0;i<(unsigned)degree_of_freedom;i++) { if (e1!=-1) { et[0]=v[e1][i]; #ifdef DEBUG printf("<=="MDF","MDF"==>",et[0],v[e1][i]); #endif } else { et[0]=0.0; } if (e2!=-1) { et[1]=v[e2][i]; #ifdef DEBUG printf("<=="MDF","MDF"==>",et[1],v[e2][i]); #endif } else { et[1]=0.0; } if (e3!=-1) { et[2]=v[e3][i]; #ifdef DEBUG printf("<=="MDF","MDF"==>",et[2],v[e3][i]); #endif } else { et[2]=0.0; } #ifdef DEBUG 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]); #endif Mzz(t->p0,3*i,Ms)=0.0; Mzz(t->p0,3*i+1,Ms)=0.0; Mzz(t->p1,3*i,Ms)=0.0; Mzz(t->p1,3*i+1,Ms)=0.0; Mzz(t->p2,3*i,Ms)=0.0; Mzz(t->p2,3*i+1,Ms)=0.0; for (j=0; j<3; j++) { Mzz(t->p0,3*i,Ms)+=et[j]*L[j]*(A[j]+B[j]*y[0]); Mzz(t->p0,3*i+1,Ms)+=et[j]*L[j]*(C[j]+D[j]*x[0]); Mzz(t->p1,3*i,Ms)+=et[j]*L[j]*(A[j]+B[j]*y[1]); Mzz(t->p1,3*i+1,Ms)+=et[j]*L[j]*(C[j]+D[j]*x[1]); Mzz(t->p2,3*i,Ms)+=et[j]*L[j]*(A[j]+B[j]*y[2]); Mzz(t->p2,3*i+1,Ms)+=et[j]*L[j]*(C[j]+D[j]*x[2]); } /* z component */ Mzz(t->p0,3*i+2,Ms)=v[renumber[t->p0]][i]; Mzz(t->p1,3*i+2,Ms)=v[renumber[t->p1]][i]; Mzz(t->p2,3*i+2,Ms)=v[renumber[t->p2]][i]; /* find min, max value for plotting */ if ( ABS(Mzz(t->p0,3*i,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p0,3*i,Ms)); if ( ABS(Mzz(t->p0,3*i+1,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p0,3*i+1,Ms)); if ( ABS(Mzz(t->p0,3*i+2,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p0,3*i+2,Ms)); if ( ABS(Mzz(t->p1,3*i,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p1,3*i,Ms)); if ( ABS(Mzz(t->p1,3*i+1,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p1,3*i+1,Ms)); if ( ABS(Mzz(t->p1,3*i+2,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p1,3*i+2,Ms)); if ( ABS(Mzz(t->p2,3*i,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p2,3*i,Ms)); if ( ABS(Mzz(t->p2,3*i+1,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p2,3*i+1,Ms)); if ( ABS(Mzz(t->p2,3*i+2,Ms)) > max_value[i] ) max_value[i]=ABS(Mzz(t->p2,3*i+2,Ms)); #ifdef DEBUG 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)); #endif } } t=DAG_traverse_prune_list(&Dt); } /* normalize all potentials by dividing their max values */ for (i=0; iTOL) { Mzz(i,3*j,Ms) /=max_value[j]; Mzz(i,3*j+1,Ms) /=max_value[j]; Mzz(i,3*j+2,Ms) /=max_value[j]; } } } g_maxpot = 1.0; g_minpot = -1.0; #ifdef DEBUG printf("potentials max=%lf min=%lf\n",g_maxpot,g_minpot); #endif /**********************************************************/ /* destroy hash table */ htbl_destroy(&H); for (j = 0; j < NN; j++) { free(v[j]); } SPM_destroy(&S); SPM_destroy(&T); free(v); free(renumber); free(max_value); free(ee); return(1); }