SPModel: plugins/isostaticFlexure/flex2d.c Source File
VPAC - Computational Software Development
Main | SPModel | StGermain FrameWork |
Main Page | Alphabetical List | Class List | Directories | File List | Class Members | File Members

flex2d.c

Go to the documentation of this file.
00001 #include <mpi.h>
00002 #include <StGermain/StGermain.h>
00003 #include <Cascade/cascade.h>
00004 
00005 #include "SPModel/SPModel.h"
00006 
00007 #include <stdio.h>
00008 #include <stdlib.h>
00009 #include <string.h>
00010 #include <assert.h>
00011 #include <math.h>
00012 
00013 void rlft3(double ***data, double **speq, long nn1, long nn2,
00014        long nn3, int isign)
00015 {
00016   void fourn(double data[], long nn[], int ndim, int isign);
00017   void nrerror(char error_text[]);
00018   long i1,i2,i3,j1,j2,j3,nn[4],ii3;
00019   double theta,wi,wpi,wpr,wr,wtemp;
00020   double c1,c2,h1r,h1i,h2r,h2i;
00021 
00022   if (1+&data[nn1][nn2][nn3]-&data[1][1][1] != nn1*nn2*nn3)
00023     printf("rlft3: problem with dimensions or contiguity of data array\n");
00024   c1=0.5;
00025   c2 = -0.5*isign;
00026   theta=isign*(6.28318530717959/nn3);
00027   wtemp=sin(0.5*theta);
00028   wpr = -2.0*wtemp*wtemp;
00029   wpi=sin(theta);
00030   nn[1]=nn1;
00031   nn[2]=nn2;
00032   nn[3]=nn3 >> 1;
00033   if (isign == 1) {
00034     fourn(&data[1][1][1]-1,nn,3,isign);
00035     for (i1=1;i1<=nn1;i1++)
00036       for (i2=1,j2=0;i2<=nn2;i2++) {
00037     speq[i1][++j2]=data[i1][i2][1];
00038     speq[i1][++j2]=data[i1][i2][2];
00039       }
00040   }
00041   for (i1=1;i1<=nn1;i1++) {
00042     j1=(i1 != 1 ? nn1-i1+2 : 1);
00043     wr=1.0;
00044     wi=0.0;
00045     for (ii3=1,i3=1;i3<=(nn3>>2)+1;i3++,ii3+=2) {
00046       for (i2=1;i2<=nn2;i2++) {
00047     if (i3 == 1) {
00048       j2=(i2 != 1 ? ((nn2-i2)<<1)+3 : 1);
00049       h1r=c1*(data[i1][i2][1]+speq[j1][j2]);
00050       h1i=c1*(data[i1][i2][2]-speq[j1][j2+1]);
00051       h2i=c2*(data[i1][i2][1]-speq[j1][j2]);
00052       h2r= -c2*(data[i1][i2][2]+speq[j1][j2+1]);
00053       data[i1][i2][1]=h1r+h2r;
00054       data[i1][i2][2]=h1i+h2i;
00055       speq[j1][j2]=h1r-h2r;
00056       speq[j1][j2+1]=h2i-h1i;
00057     } else {
00058       j2=(i2 != 1 ? nn2-i2+2 : 1);
00059       j3=nn3+3-(i3<<1);
00060       h1r=c1*(data[i1][i2][ii3]+data[j1][j2][j3]);
00061       h1i=c1*(data[i1][i2][ii3+1]-data[j1][j2][j3+1]);
00062       h2i=c2*(data[i1][i2][ii3]-data[j1][j2][j3]);
00063       h2r= -c2*(data[i1][i2][ii3+1]+data[j1][j2][j3+1]);
00064       data[i1][i2][ii3]=h1r+wr*h2r-wi*h2i;
00065       data[i1][i2][ii3+1]=h1i+wr*h2i+wi*h2r;
00066       data[j1][j2][j3]=h1r-wr*h2r+wi*h2i;
00067       data[j1][j2][j3+1]= -h1i+wr*h2i+wi*h2r;
00068     }
00069       }
00070       wr=(wtemp=wr)*wpr-wi*wpi+wr;
00071       wi=wi*wpr+wtemp*wpi+wi;
00072     }
00073   }
00074   if (isign == -1)
00075     fourn(&data[1][1][1]-1,nn,3,isign);
00076 }
00077 
00078 
00079 #define nSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
00080 
00081 void fourn(double data[], long nn[], int ndim, int isign)
00082 {
00083   int idim;
00084   long i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
00085   long ibit,k1,k2,n,nprev,nrem,ntot;
00086   double tempi,tempr;
00087   double theta,wi,wpi,wpr,wr,wtemp;
00088 
00089   for (ntot=1,idim=1;idim<=ndim;idim++)
00090     ntot *= nn[idim];
00091   nprev=1;
00092   for (idim=ndim;idim>=1;idim--) {
00093     n=nn[idim];
00094     nrem=ntot/(n*nprev);
00095     ip1=nprev << 1;
00096     ip2=ip1*n;
00097     ip3=ip2*nrem;
00098     i2rev=1;
00099     for (i2=1;i2<=ip2;i2+=ip1) {
00100       if (i2 < i2rev) {
00101     for (i1=i2;i1<=i2+ip1-2;i1+=2) {
00102       for (i3=i1;i3<=ip3;i3+=ip2) {
00103         i3rev=i2rev+i3-i2;
00104         nSWAP(data[i3],data[i3rev]);
00105         nSWAP(data[i3+1],data[i3rev+1]);
00106       }
00107     }
00108       }
00109       ibit=ip2 >> 1;
00110       while (ibit >= ip1 && i2rev > ibit) {
00111     i2rev -= ibit;
00112     ibit >>= 1;
00113       }
00114       i2rev += ibit;
00115     }
00116     ifp1=ip1;
00117     while (ifp1 < ip2) {
00118       ifp2=ifp1 << 1;
00119       theta=isign*6.28318530717959/(ifp2/ip1);
00120       wtemp=sin(0.5*theta);
00121       wpr = -2.0*wtemp*wtemp;
00122       wpi=sin(theta);
00123       wr=1.0;
00124       wi=0.0;
00125       for (i3=1;i3<=ifp1;i3+=ip1) {
00126     for (i1=i3;i1<=i3+ip1-2;i1+=2) {
00127       for (i2=i1;i2<=ip3;i2+=ifp2) {
00128         k1=i2;
00129         k2=k1+ifp1;
00130         tempr=(double)wr*data[k2]-(double)wi*data[k2+1];
00131         tempi=(double)wr*data[k2+1]+(double)wi*data[k2];
00132         data[k2]=data[k1]-tempr;
00133         data[k2+1]=data[k1+1]-tempi;
00134         data[k1] += tempr;
00135         data[k1+1] += tempi;
00136       }
00137     }
00138     wr=(wtemp=wr)*wpr-wi*wpi+wr;
00139     wi=wi*wpr+wtemp*wpi+wi;
00140       }
00141       ifp1=ifp2;
00142     }
00143     nprev *= n;
00144   }
00145 }
00146 #undef nSWAP
00147 
00148 
00149 void cosfilt(double l3,long nn3,double l2,long nn2,double **qtot,
00150          double dflex,double drho,double grav,double **w,double ***load,double **speq) {
00151 
00152 
00153   /*solves the 2D elastic flexural equation on a rectangular region using a cosine transform*/
00154   /*the solution has zero slope on the boundary, perpendicular to the boundary of the domain */
00155 
00156   double fac,f2,f3,f2i2,f3i3,fac1,fac2;
00157   long i,j,i1=1,i2,i3,nn1=1,isign;
00158 
00159   /*the cosine transform is obtained in the brute force fashion by extendnig the load*/
00160   /*periodically, enlarging the data set by a factor of 4; waste of time, but what the heck it's fast enough anyway*/
00161 
00162   /*printf ( "Lx %lf, Ly %lf, rnx %li, rny %li, grav %lf, flex %lf, rho %lf\n", l3, l2, nn3, nn2, grav, dflex, drho );*/
00163   
00164   /*this is moved to global level  
00165   load=d3tensor(1,nn1,1,2*nn2,1,2*nn3);
00166   speq=dmatrix(1,nn1,1,2*2*nn2);
00167   */
00168 /*  printf( "printing the input\n" );
00169   for (i=1;i<=nn3;i++){
00170     for (j=1;j<=nn2;j++){
00171         printf( "loadg[%d][%d] = %lf\n", i, j, qtot[i][j] );
00172     }
00173   }*/
00174   
00175   /*copy load into right locations*/
00176   /*lower left*/
00177   for (i=1;i<=nn3;i++)
00178     for (j=1;j<=nn2;j++) load[1][j][i]=qtot[i][j];
00179   
00180   /*lower right*/
00181   for (i=1;i<=nn3;i++)
00182     for (j=1;j<=nn2;j++) load[1][j][nn3+i]=qtot[nn3-i+1][j];
00183 
00184   /*upper left*/
00185   for (i=1;i<=nn3;i++)
00186     for (j=1;j<=nn2;j++) load[1][nn2+j][i]=qtot[i][nn2-j+1];
00187 
00188   /*upper right*/
00189   for (i=1;i<=nn3;i++)
00190     for (j=1;j<=nn2;j++) load[1][nn2+j][nn3+i]=qtot[nn3-i+1][nn2-j+1];
00191 
00192 
00193   
00194   /*fft of load*/
00195   isign=1;
00196   rlft3(load,speq,nn1,2*nn2,2*nn3,isign);
00197   
00198 
00199   fac1=nn1*2*nn2*2*nn3/2.0;
00200   fac2=drho*grav;
00201   f2=2.0*PI/(2*l2);
00202   f2*=f2;
00203   f3=2.0*PI/(2*l3);
00204   f3*=f3;
00205   
00206   for (i2=1;i2<=2*nn2/2;i2++)     /*y*/
00207     for (i3=1;i3<=2*nn3/2;i3++) {/*x*/
00208       f2i2=f2*(i2-1)*(i2-1);
00209       f3i3=f3*(i3-1)*(i3-1);
00210       fac=dflex*(f2i2*f2i2+f3i3*f3i3+2.0*f2i2*f3i3)+fac2;
00211       fac=1.0/fac/fac1;
00212       load[i1][i2][2*i3-1]*=fac;
00213       load[i1][i2][2*i3  ]*=fac;
00214     }/*j*/
00215 
00216   for (i2=2*nn2;i2>=2*nn2/2+1;i2--)   /*y*/
00217     for (i3=1;i3<=2*nn3/2;i3++) {      /*x*/
00218       f2i2=f2*(2*nn2+1-i2)*(2*nn2+1-i2);
00219       f3i3=f3*(i3-1)*(i3-1);
00220       fac=dflex*(f2i2*f2i2+f3i3*f3i3+2.0*f2i2*f3i3)+fac2;
00221       fac=1.0/fac/fac1;
00222       load[i1][i2][2*i3-1]*=fac;
00223       load[i1][i2][2*i3  ]*=fac;
00224     }/*j*/
00225 
00226   for (i2=1;i2<=2*nn2/2;i2++)       /*y*/
00227     for (i3=2*nn3/2+1;i3<=2*nn3/2+1;i3++) {/*x*/
00228       f2i2=f2*(i2-1)*(i2-1);
00229       f3i3=f3*(i3-1)*(i3-1);
00230       fac=dflex*(f2i2*f2i2+f3i3*f3i3+2.0*f2i2*f3i3)+fac2;
00231       fac=1.0/fac/fac1;
00232       speq[i1][2*i2-1]*=fac;
00233       speq[i1][2*i2  ]*=fac;
00234     }/*j*/
00235 
00236   for (i2=2*nn2;i2>=2*nn2/2+1;i2--)   /*y*/
00237     for (i3=2*nn3/2+1;i3<=2*nn3/2+1;i3++) {      /*x*/
00238       f2i2=f2*(2*nn2+1-i2)*(2*nn2+1-i2);
00239       f3i3=f3*(i3-1)*(i3-1);
00240       fac=dflex*(f2i2*f2i2+f3i3*f3i3+2.0*f2i2*f3i3)+fac2;
00241       fac=1.0/fac/fac1;
00242       speq[i1][2*i2-1]*=fac;
00243       speq[i1][2*i2  ]*=fac;
00244     }/*j*/
00245 
00246 
00247   /*reverse fft of load*/
00248   isign=-1;
00249   rlft3(load,speq,nn1,2*nn2,2*nn3,isign);
00250 
00251   /*copy deflection into right locations*/
00252   for (i=1;i<=nn3;i++)
00253     for (j=1;j<=nn2;j++) w[i][j]=-load[1][j][i];
00254   
00255 /*  printf( "printing the output\n" );
00256   for (i=1;i<=nn3;i++)
00257     for (j=1;j<=nn2;j++)
00258         printf( "w[%d][%d] = %lf\n", i, j, w[i][j] );*/
00259   
00260 }/*cosfilt*/
00261 
00262 
00263 void filt(double l3, long nn3, double l2,long nn2,double **qtot,
00264       double dflex,double drho,double grav,double**w) {
00265 
00266   double fac,f2,f3,f2i2,f3i3,fac1,fac2,pi=3.1415926;
00267   double  ***load,**speq;
00268   long i,j,i1=1,i2,i3,nn1=1,isign;
00269 
00270   load=d3tensor(1,nn1,1,nn2,1,nn3);
00271   speq=dmatrix(1,nn1,1,2*nn2);
00272 
00273   /*copy load into right locations*/
00274   for (i=1;i<=nn3;i++)
00275     for (j=1;j<=nn2;j++) load[1][j][i]=qtot[i][j];
00276 
00277 
00278   /*fft of load*/
00279   isign=1;
00280   rlft3(load,speq,nn1,nn2,nn3,isign);
00281 
00282 
00283   fac1=nn1*nn2*nn3/2.0;
00284   fac2=drho*grav;
00285   f2=2.0*pi/l2;
00286   f2*=f2;
00287   f3=2.0*pi/l3;
00288   f3*=f3;
00289 
00290   for (i2=1;i2<=nn2/2;i2++)     /*y*/
00291     for (i3=1;i3<=nn3/2;i3++) {/*x*/
00292       f2i2=f2*(i2-1)*(i2-1);
00293       f3i3=f3*(i3-1)*(i3-1);
00294       fac=dflex*(f2i2*f2i2+f3i3*f3i3+2.0*f2i2*f3i3)+fac2;
00295       fac=1.0/fac/fac1;
00296       load[i1][i2][2*i3-1]*=fac;
00297       load[i1][i2][2*i3  ]*=fac;
00298     }/*j*/
00299 
00300   for (i2=nn2;i2>=nn2/2+1;i2--)   /*y*/
00301     for (i3=1;i3<=nn3/2;i3++) {      /*x*/
00302       f2i2=f2*(nn2+1-i2)*(nn2+1-i2);
00303       f3i3=f3*(i3-1)*(i3-1);
00304       fac=dflex*(f2i2*f2i2+f3i3*f3i3+2.0*f2i2*f3i3)+fac2;
00305       fac=1.0/fac/fac1;
00306       load[i1][i2][2*i3-1]*=fac;
00307       load[i1][i2][2*i3  ]*=fac;
00308     }/*j*/
00309 
00310   for (i2=1;i2<=nn2/2;i2++)       /*y*/
00311     for (i3=nn3/2+1;i3<=nn3/2+1;i3++) {/*x*/
00312       f2i2=f2*(i2-1)*(i2-1);
00313       f3i3=f3*(i3-1)*(i3-1);
00314       fac=dflex*(f2i2*f2i2+f3i3*f3i3+2.0*f2i2*f3i3)+fac2;
00315       fac=1.0/fac/fac1;
00316       speq[i1][2*i2-1]*=fac;
00317       speq[i1][2*i2  ]*=fac;
00318     }/*j*/
00319 
00320   for (i2=nn2;i2>=nn2/2+1;i2--)   /*y*/
00321     for (i3=nn3/2+1;i3<=nn3/2+1;i3++) {      /*x*/
00322       f2i2=f2*(nn2+1-i2)*(nn2+1-i2);
00323       f3i3=f3*(i3-1)*(i3-1);
00324       fac=dflex*(f2i2*f2i2+f3i3*f3i3+2.0*f2i2*f3i3)+fac2;
00325       fac=1.0/fac/fac1;
00326       speq[i1][2*i2-1]*=fac;
00327       speq[i1][2*i2  ]*=fac;
00328     }/*j*/
00329 
00330 
00331   /*reverse fft of load*/
00332   isign=-1;
00333   rlft3(load,speq,nn1,nn2,nn3,isign);
00334 
00335   /*copy deflection into right locations*/
00336   for (i=1;i<=nn3;i++)
00337     for (j=1;j<=nn2;j++) w[i][j]=-load[1][j][i];
00338   
00339   free_d3tensor(load,1,nn1,1,nn2,1,nn3);
00340   free_dmatrix(speq,1,nn1,1,2*nn2);
00341 }