00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
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
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
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
00144
00145 points2df.resize(NI+NF);
00146
00147
00148 for(int i=0; i< NI+NF; i++){
00149 fscanf(fpt,"%f %f",&points2df[i].X(),&points2df[i].Y());
00150
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
00173 int NI = atoi(argv[1]);
00174 int NF = 4*atoi(argv[2]);
00175
00176 int N = NI + NF;
00177
00178
00179
00180
00181 points2df.resize(N);
00182
00183
00184
00185
00186
00187
00188
00189
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
00204 lambdaf.resize(N+6);
00205 uf.resize(N+6);
00206
00207
00208
00209 knots(NI,NF,0.01);
00210
00211
00212
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
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
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
00258 GrammMatrix<float> G(Wl,Pl,Wb,Pb,P);
00259
00260
00261
00262 for (int i=0; i<NI; i++)
00263 uf[i] = f(points2df[i].X(),points2df[i].Y());
00264
00265
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
00281 for (int i=N; i<N+6; i++)
00282 uf[i] = 0.0;
00283
00284
00285 lambdaf = G/uf;
00286
00287
00288
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 }