tps_poisson.cpp
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 .
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
00039 #include
00040 #include
00041 #include
00042 #include
00043 #include
00044 #include
00045 #include
00046 #include
00047 #include
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 tps)
00077 {
00078 int N = points2df.size();
00079
00080 float interp = 0.0;
00081
00082 for(int i=0; i0; 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() < tps_2dxf;
00160 TPS_2DY tps_2dyf;
00161 TPS4 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úmero 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 G(Wl,Pl,Wb,Pb,P);
00259
00260
00261 //Se calcula el vector de soluciones exactas
00262 for (int i=0; i