00001 /******************************************************************************
00002 
00003  Copyright 2008 Departamento de Realidad Virtual
00004  y Unidad de Cómputo Aplicado DGSGA, UNAM.
00005 
00006 
00007  This file is part of RBF++.
00008 
00009  RBF++ is free software: you can redistribute it and/or modify
00010  it under the terms of the GNU General Public License as published by
00011  the Free Software Foundation, either version 3 of the License, or
00012  (at your option) any later version.
00013 
00014  RBF++ is distributed in the hope that it will be useful,
00015  but WITHOUT ANY WARRANTY; without even the implied warranty of
00016  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00017  GNU General Public License for more details.
00018 
00019  You should have received a copy of the GNU General Public License
00020  along with RBF++. If not, see <http://www.gnu.org/licenses/>.
00021 
00022 
00023 *******************************************************************************/
00024 
00025 /* 
00026  Author: Daniel Cervantes Cabrera
00027  Supervision: Pedro Gonzalez Casanova
00028  Project: RBF++
00029  Institution: DGSCA, UNAM.
00030  Date: 3/05/08
00031  Description: Poisson equation with boundary conditions solved with
00032  unsymetric colocation method.
00033  For references check (for now only in spanish):
00034  http://www.dci.dgsca.unam.mx/pderbf/media/PoissonPDERBF.pdf
00035  */
00036 
00037 
00038 #include <time.h>
00039 #include <math.h>
00040 #include <stdlib.h>
00041 #include <stdio.h>
00042 #include <iostream>
00043 #include <RBF/RBF.h>
00044 #include <RBF/GrammMatrix.h>
00045 #include <LA/Matrix.h>
00046 #include <LA/Point2dT.h>
00047 #include <LA/Point3dT.h>
00048 
00049 
00050 using namespace RBF;
00051 
00052 
00053 inline float f(float x, float y)
00054 {
00055   if (0 <= x &&  x <= 1 && 0 <= y && y <= 1)
00056     return (2*(5.0/4.0 + cos(5.4*y))* pow((108.0*x-36.0),2)   )/ ( pow( ( 6.0+ 6.0*  (pow((3.0*x-1.0),2))   )  ,3)  )  -(108.0*(5.0/4.0+cos(5.4*y)))  / pow( (6.0+6.0*(pow((3.0*x-1.0),2))) ,2)-  (29.16*cos(5.4*y)) / (6.0+6.0* pow((3.0*x-1),2) );   
00057   else
00058     return 0;             
00059 }
00060 
00061 inline float g1(float y){
00062   return 5.0/48.0+(1.0/12.0)*cos(5.4*y);
00063 }
00064 inline float g2(float y){
00065    return  1.0/24.0+1.0/30.0*cos(5.4*y);
00066 }
00067 
00068 inline float g3(float x){
00069   return 2.25/(6+6*pow((3*x-1),2));
00070 }
00071 inline float g4(float x){
00072   return 1.884692876/(6+6*pow((3*x-1),2));
00073 }
00074 
00075 
00076 float eval_sol(float x, float y,RBF::TPS4<float> tps)
00077 {
00078   int N = points2df.size();
00079   
00080   float interp = 0.0;  
00081 
00082   for(int i=0; i<N; i++)
00083     interp += lambdaf[i]*tps(x,y,points2df[i].X(),points2df[i].Y());
00084     
00085   interp += lambdaf[N]*x*x+lambdaf[N+1]*y*y+
00086             lambdaf[N+2]*x*y+lambdaf[N+3]*x+
00087             lambdaf[N+4]*y+lambdaf[N+5];
00088 
00089   return interp;  
00090 } 
00091 
00092 
00093 void knots(int NI, int NF, float epsilon){
00094   srand(time(0));  
00095   int nf = NF/4;
00096   int N = NI+NF;
00097   
00098   //Nodos interiores
00099   for (int i=0; i< NI; i++){   
00100     points2df[i].X() = epsilon + (float(rand())/RAND_MAX)*(1.0-2.0*epsilon);
00101     points2df[i].Y() = epsilon + (float(rand())/RAND_MAX)*(1.0-2.0*epsilon);  
00102   }
00103   
00104   //Nodos frontera
00105   for (int i=NI, k=0; k<nf; i++,k++){
00106     points2df[i].X() = float(k)/float(nf);
00107     points2df[i].Y() = 0.0;
00108   }
00109    
00110   for (int i=NI+nf, k=0; k<nf; i++,k++){
00111     points2df[i].X() = 1.0;
00112     points2df[i].Y() = float(k)/float(nf);
00113   }
00114  
00115   for (int i=NI+2*nf, k=nf; k>0; i++,k--){
00116     points2df[i].X() = float(k)/float(nf);
00117     points2df[i].Y() = 1.0;
00118   }
00119 
00120   for (int i=NI+3*nf, k=nf; k>0; i++,k--){
00121     points2df[i].X() = 0.0;
00122     points2df[i].Y() = float(k)/float(nf);
00123   } 
00124 
00125   FILE * fpt = fopen("puntos2d.txt","w");
00126   
00127   int size = points2df.size();
00128 
00129   fprintf(fpt,"%d %d\n",NI,NF);
00130 
00131   for(int i=0; i < size; i++) 
00132     fprintf(fpt,"%3.6f %3.6f\n",points2df[i].X(),points2df[i].Y());
00133   
00134   fclose(fpt);  
00135 
00136 }
00137 
00138 
00139 void knots(int & NI, int & NF, char nombre[]){
00140   FILE * fpt = fopen(nombre,"r");
00141   
00142   fscanf(fpt,"%d %d",&NI,&NF);
00143   // std::cout << NI <<" " << NF << std::endl;
00144   
00145   points2df.resize(NI+NF);
00146   //y.resize(NI+NF);
00147 
00148   for(int i=0; i< NI+NF; i++){
00149     fscanf(fpt,"%f %f",&points2df[i].X(),&points2df[i].Y());
00150     // std::cout << points2df[i].X() <<" "<< points2df[i].Y() <<std::endl;
00151   }
00152   
00153   fclose(fpt);
00154 } 
00155 
00156 int main(int argc, char * argv[]){
00157 
00158   LA::Matrixf Wl,Pl,Wb,Pb,P;  
00159   TPS_2DX<float> tps_2dxf;
00160   TPS_2DY<float> tps_2dyf;
00161   TPS4<float> tps;
00162 
00163   
00164   if (argc <= 2){ 
00165       std::cout << "usar ./sol_poisson na nif " << std::endl
00166                 << "na = numero de nodos aleatorios" << std::endl
00167                 << "nif = numero de intervalos en cada lado de la frontera" << std::endl;
00168       exit(1);
00169     }
00170 
00171 
00172   //Se lee el n'umero de puntos aleatorios interiores y de frontera.
00173   int NI = atoi(argv[1]);  
00174   int NF = 4*atoi(argv[2]);  
00175   
00176   int N = NI + NF;
00177   
00178   
00179   
00180   //Se establece el tama~o en los vectores de nodos
00181   points2df.resize(N);
00182   
00183   // Se establece el tama~o de la matrices de la 
00184   // de la matriz global del sistema algebraico.
00185   
00186   /*
00187  | Wl Pl | 
00188  | Wb Pb |
00189  | Pbt 0 | 
00190  */
00191 
00192   
00193   Wl.resize(NI,N);
00194   Pl.resize(NI,6); 
00195 
00196   Wb.resize(NF,N);
00197   Pb.resize(NF,6);
00198 
00199   P.resize(6,N);
00200   
00201  
00202 
00203   //Se estable el tama~o de los vectores lambda y u.
00204   lambdaf.resize(N+6);
00205   uf.resize(N+6);
00206 
00207   //Se generan los nodos aleatorios entre 0 y 1
00208   //asi como los nodos frontera.
00209   knots(NI,NF,0.01);
00210     
00211 
00212   //Se crea la matriz Wl y Pl 
00213   for(int i=0; i<NI; i++)
00214     for(int j=0; j<N; j++)
00215       Wl(i,j) = tps_2dxf(points2df[i].X(),points2df[i].Y(),points2df[j].X(),points2df[j].Y())+
00216                 tps_2dyf(points2df[i].X(),points2df[i].Y(),points2df[j].X(),points2df[j].Y());  
00217 
00218   
00219   for(int i=0; i<NI; i++)
00220     for(int j=0; j<6; j++){
00221       Pl(i,0) = 2.0; 
00222       Pl(i,1) = 2.0; 
00223       Pl(i,2) = 0.0; 
00224       Pl(i,3) = 0.0;
00225       Pl(i,4) = 0.0; 
00226       Pl(i,5) = 0.0;
00227     }
00228 
00229 
00230   //Se crea la matriz Wb y Pb 
00231    for(int i=0; i<NF; i++)
00232     for(int j=0; j<N; j++)
00233       Wb(i,j) = tps(points2df[i+NI].X(),points2df[i+NI].Y(),points2df[j].X(),points2df[j].Y());   
00234   
00235 
00236   for(int i=0; i<NF; i++){      
00237     Pb(i,0) = points2df[NI+i].X()*points2df[i+NI].X();
00238     Pb(i,1) = points2df[i+NI].Y()*points2df[i+NI].Y();
00239     Pb(i,2) = points2df[i+NI].X()*points2df[i+NI].Y();
00240     Pb(i,3) = points2df[i+NI].X();
00241     Pb(i,4) = points2df[i+NI].Y();
00242     Pb(i,5) = 1.0;      
00243   } 
00244   
00245 
00246   //Se crea la matriz que impone condiciones de momentos.
00247   for(int i=0; i<N; i++){
00248     P(0,i) = points2df[i].X()*points2df[i].X();
00249     P(1,i) = points2df[i].Y()*points2df[i].Y();
00250     P(2,i) = points2df[i].X()*points2df[i].Y();
00251     P(3,i) = points2df[i].X();
00252     P(4,i) = points2df[i].Y();
00253     P(5,i) = 1.0;
00254   }
00255 
00256 
00257   //Se crea la matriz de Gramm que integra las matrices anteriores.
00258   GrammMatrix<float> G(Wl,Pl,Wb,Pb,P);
00259   
00260           
00261   //Se calcula el vector de soluciones exactas
00262   for (int i=0; i<NI; i++)
00263     uf[i] = f(points2df[i].X(),points2df[i].Y());
00264 
00265   // Condiciones en las fronteras.
00266   int nf = NF/4;
00267   
00268   for(int i=NI; i<NI+nf; i++)
00269     uf[i] = g3(points2df[i].X());
00270   
00271   for(int i=NI+nf; i<NI+2*nf; i++)
00272     uf[i] = g2(points2df[i].Y());
00273   
00274   for(int i=NI+2*nf; i<NI+3*nf; i++)
00275     uf[i] = g4(points2df[i].X());
00276   
00277   for(int i=NI+3*nf; i<NI+4*nf; i++)
00278     uf[i] = g1(points2df[i].Y());  
00279   
00280   // Momentos
00281   for (int i=N; i<N+6; i++)
00282     uf[i] = 0.0;
00283     
00284   //Solucion del sistema utilizando eliminación Gaussiana.
00285    lambdaf = G/uf;
00286 
00287 
00288    //Evaluación del la solución sobre una malla de 100x100.
00289    float xmi,ymi;
00290         
00291    FILE * fpt  = fopen("sol_aprox.txt","w");
00292    
00293          
00294    for(int i=0; i <= 100; i++)
00295      for(int j=0; j <= 100; j++){
00296        ymi = float(i)/100;
00297        xmi = float(j)/100;
00298        fprintf(fpt, "%4.6f %4.6f %4.6f\n", xmi,ymi,eval_sol(xmi,ymi,tps));        
00299      }                        
00300 
00301    fclose(fpt);
00302 }