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 }
Superior

Tenemos 260 visitantes y ningun miembro en Línea