CO2INJECTION: fchessboard3d.cpp Source File

00001 #include "fchessboard3d.h"
00002 #include <fstream>
00003 #include "exception.h"
00004 #include "numericmethods.h"
00005 
00006 
00007 #define MAX_LINE_SIZE 10000
00008 using std::ifstream;
00009 FChessBoard3D::FChessBoard3D(Point3D p0,Point3D p1,Point3D lp0,Point3D lp1,std::string fileName) 
00010   :m_p0(p0),m_p1(p1)
00011  {
00012    assert(lp0[0] <= lp1[0]);
00013    assert(lp0[1] <= lp1[1]);
00014    assert(lp0[2] <= lp1[2]);
00015    assert(p0[0] <= p1[0]);
00016    assert(p0[1] <= p1[1]);
00017    assert(p0[2] <= p1[2]);
00018    assert(lp0[0] >= p0[0]);
00019    assert(lp0[1] >= p0[1]);
00020    assert(lp0[2] >= p0[2]);
00021    assert(lp1[0] <= p1[0]);
00022    assert(lp1[1] <= p1[1]);
00023    assert(lp1[2] <= p1[2]);
00024    double ddLixo;
00025 
00026 
00027 
00028    
00029    ifstream file(fileName.c_str());
00030    if (file.fail())
00031      throw new Exception("FChessBoard3D: File \"%s\" not found",fileName.c_str());
00032 
00033    char str[MAX_LINE_SIZE];
00034 
00035    //Ignore the first line
00036    file.getline(str,MAX_LINE_SIZE);
00037    
00038    //Get the second  line containing the dimension of the chess board
00039    //in the format "size_X x size_Y x size_Z"
00040    file.getline(str,MAX_LINE_SIZE);
00041    if (file.fail() || sscanf(str,"%d x %d x %d",&m_dimX,&m_dimY,&m_dimZ) != 3)
00042      throw new Exception("FChessBoard3D: File \"%s\", Second line  (\"%s\")\n\tnot in format <number> x <number> x <number>",fileName.c_str(),str);
00043 
00044 
00045 
00046    m_dX = (p1[0]-p0[0])/m_dimX;
00047    m_dY = (p1[1]-p0[1])/m_dimY;
00048    m_dZ = (p1[2]-p0[2])/m_dimZ;
00049 
00050    Xi[0] = static_cast<unsigned>(floor(lp0[0]/m_dX));
00051    Xi[1] = static_cast<unsigned>(floor(lp0[1]/m_dY));
00052    Xi[2] = static_cast<unsigned>(floor(lp0[2]/m_dZ));
00053 
00054    Xe[0] = static_cast<unsigned>(ceil(lp1[0]/m_dX));
00055    Xe[1] = static_cast<unsigned>(ceil(lp1[1]/m_dY));
00056    Xe[2] = static_cast<unsigned>(ceil(lp1[2]/m_dZ));
00057 
00058    //Just to garantee that roundoff errors would not occur.
00059    //Actually it never ocurred, but just to make sure.
00060    if (Xe[0] > m_dimX)
00061      Xe[0]=m_dimX;
00062    if (Xe[1] > m_dimY)
00063      Xe[1]=m_dimY;
00064    if (Xe[2] > m_dimZ)
00065      Xe[2]=m_dimZ;
00066 
00067 
00068    //Allocate the table to store the chess board
00069    TableIndices<3> sizes(Xe[2]-Xi[2],Xe[1]-Xi[1],Xe[0]-Xi[0]);
00070    //TableIndices<3> sizes(Xe[0]-Xi[0]+1,Xe[1]-Xi[1]+1,Xe[2]-Xi[2]+1);
00071    m_matrix3D.reinit(sizes);
00072 
00073 
00074    m_lp0[0]=Xi[0]*m_dX;
00075    m_lp0[1]=Xi[1]*m_dY;
00076    m_lp0[2]=Xi[2]*m_dZ;
00077 
00078    m_lp1[0]=Xe[0]*m_dX;
00079    m_lp1[1]=Xe[1]*m_dY;
00080    m_lp1[2]=Xe[2]*m_dZ;
00081 
00082    
00083    /*printf("========================\n%d) Points %g,%g,%g - %g,%g,%g => %g %g %g - %g %g %g\n",NetMPI::rank(),p0[0],p0[1],p0[2],p1[0],p1[1],p1[2],lp0[0],
00084             lp0[1],lp0[2],lp1[0],lp1[1],lp1[2]);
00085    printf("========================\n%d) Indices %d,%d,%d - %d,%d,%d\n",NetMPI::rank(),Xi[0],Xi[1],Xi[2],Xe[0],Xe[1],Xe[2]);
00086    printf("========================\n%d) Minimal Geo Mesh %g,%g,%g - %g,%g,%g\n",NetMPI::rank(),m_lp0[0],m_lp0[1],m_lp0[2],m_lp1[0],m_lp1[1],m_lp1[2]);
00087    */
00088    
00089    
00090    for (unsigned z=0;z<m_dimZ;z++)
00091      for (unsigned y=0;y<m_dimY;y++)
00092        for (unsigned x=0;x<m_dimX;x++)
00093        {
00094          if (x >= Xi[0] && x < Xe[0] &&
00095              y >= Xi[1] && y < Xe[1] &&
00096              z >= Xi[2] && z < Xe[2] )
00097          {
00098            //printf("%d) Reading %d,%d,%d\n",NetMPI::rank(),x,y,z); 
00099            file >> m_matrix3D(z-Xi[2],y-Xi[1],x-Xi[0]);
00100          }
00101          else
00102            file >> ddLixo;
00103          
00104          if (file.fail())
00105            throw new Exception("FChessBoard3D: File \"%s\", dimension of the matrix does not match with the numbers supplied in the file\n",fileName.c_str());
00106        }
00107 }
00108 
00109 
00110  FChessBoard3D::FChessBoard3D(Point3D p0,Point3D p1,int dimX,int dimY,int dimZ, std::vector<double> &data) 
00111    :m_p0(p0),m_p1(p1),m_lp0(p0),m_lp1(p1)
00112 {
00113 
00114   m_dimZ = dimZ;
00115   m_dimX = dimX;
00116   m_dimY = dimY;
00117 
00118   Xi = 0;
00119   Xe[0]=m_dimX;
00120   Xe[1]=m_dimY;
00121   Xe[2]=m_dimZ;
00122   
00123   TableIndices<3> sizes(m_dimZ,m_dimY,m_dimX);
00124   std::vector<double>::iterator it = data.begin();
00125 
00126   
00127   assert(data.size() == (unsigned) dimX*dimY*dimZ);
00128    m_matrix3D.reinit(sizes);
00129    for (unsigned z=0;z<m_dimZ;z++)
00130      for (unsigned y=0;y<m_dimY;y++)
00131        for (unsigned x=0;x<m_dimX;x++)
00132        {
00133          m_matrix3D(z,y,x) = *(it++);
00134        }
00135    m_dX = (p1[0]-p0[0])/m_dimX;
00136    m_dY = (p1[1]-p0[1])/m_dimY;
00137    m_dZ = (p1[2]-p0[2])/m_dimZ;
00138 
00139 }
00140 
00141 
00142 
00143 double FChessBoard3D::operator() (const VecDouble &p, unsigned int component) const 
00144 {
00145   assert(component == 0);
00146   assert(NumericMethods::isInCube(p,m_lp0,m_lp1));
00147 
00148   unsigned x = (unsigned) floor((p(0)-m_lp0(0))/m_dX);
00149   unsigned y = (unsigned) floor((p(1)-m_lp0(1))/m_dY);
00150   unsigned z = (unsigned) floor((p(2)-m_lp0(2))/m_dZ);
00151 
00152   
00153   return m_matrix3D(z,y,x);
00154 }
00155 
00156 
00157  FChessBoard3D::~FChessBoard3D() 
00158 {
00159 
00160 }
00161 
00162 
00163 
00164 
00165 
00166 
00179  FChessBoard3D::FChessBoard3D(Point3D p0,Point3D p1,Point3D lp0,Point3D lp1,double My,double CVy,std::string fileName) 
00180   :m_p0(p0),m_p1(p1),m_lp0(lp0),m_lp1(lp1)
00181  {
00182    
00183    assert(lp0[0] <= lp1[0]);
00184    assert(lp0[1] <= lp1[1]);
00185    assert(lp0[2] <= lp1[2]);
00186    assert(p0[0] <= p1[0]);
00187    assert(p0[1] <= p1[1]);
00188    assert(p0[2] <= p1[2]);
00189    assert(lp0[0] >= p0[0]);
00190    assert(lp0[1] >= p0[1]);
00191    assert(lp0[2] >= p0[2]);
00192    assert(lp1[0] <= p1[0]);
00193    assert(lp1[1] <= p1[1]);
00194    assert(lp1[2] <= p1[2]);
00195    
00196 
00197 
00198    
00199    ifstream file(fileName.c_str());
00200    if (file.fail())
00201      throw new Exception("FChessBoard3D: File \"%s\" not found",fileName.c_str());
00202 
00203    char str[MAX_LINE_SIZE];
00204 
00205    //Ignore the first line
00206    file.getline(str,MAX_LINE_SIZE);
00207    
00208    //Get the second  line containing the dimension of the chess board
00209    //in the format "size_X x size_Y x size_Z"
00210    file.getline(str,MAX_LINE_SIZE);
00211    if (file.fail() || sscanf(str,"%d x %d x %d",&m_dimX,&m_dimY,&m_dimZ) != 3)
00212      throw new Exception("FChessBoard3D: File \"%s\", Second line  (\"%s\")\n\tnot in format <number> x <number> x <number>",fileName.c_str(),str);
00213 
00214 
00215 
00216    m_dX = (p1[0]-p0[0])/m_dimX;
00217    m_dY = (p1[1]-p0[1])/m_dimY;
00218    m_dZ = (p1[2]-p0[2])/m_dimZ;
00219 
00220    Xi[0] = static_cast<unsigned>(floor(lp0[0]/m_dX));
00221    Xi[1] = static_cast<unsigned>(floor(lp0[1]/m_dY));
00222    Xi[2] = static_cast<unsigned>(floor(lp0[2]/m_dZ));
00223 
00224    Xe[0] = static_cast<unsigned>(ceil(lp1[0]/m_dX));
00225    Xe[1] = static_cast<unsigned>(ceil(lp1[1]/m_dY));
00226    Xe[2] = static_cast<unsigned>(ceil(lp1[2]/m_dZ));
00227 
00228    //Just to garantee that roundoff errors would not occur.
00229    //Actually it never ocurred, but just to make sure.
00230    if (Xe[0] > m_dimX)
00231      Xe[0]=m_dimX;
00232    if (Xe[1] > m_dimY)
00233      Xe[1]=m_dimY;
00234    if (Xe[2] > m_dimZ)
00235      Xe[2]=m_dimZ;
00236 
00237 
00238    //Allocate the table to store the chess board
00239    TableIndices<3> sizes(Xe[2]-Xi[2],Xe[1]-Xi[1],Xe[0]-Xi[0]);
00240    //TableIndices<3> sizes(Xe[0]-Xi[0]+1,Xe[1]-Xi[1]+1,Xe[2]-Xi[2]+1);
00241    m_matrix3D.reinit(sizes);
00242 
00243 
00244    m_lp0[0]=Xi[0]*m_dX;
00245    m_lp0[1]=Xi[1]*m_dY;
00246    m_lp0[2]=Xi[2]*m_dZ;
00247 
00248    m_lp1[0]=Xe[0]*m_dX;
00249    m_lp1[1]=Xe[1]*m_dY;
00250    m_lp1[2]=Xe[2]*m_dZ;
00251 
00252    
00253    /*printf("========================\n%d) Points %g,%g,%g - %g,%g,%g => %g %g %g - %g %g %g\n",NetMPI::rank(),p0[0],p0[1],p0[2],p1[0],p1[1],p1[2],lp0[0],
00254             lp0[1],lp0[2],lp1[0],lp1[1],lp1[2]);
00255    printf("========================\n%d) Indices %d,%d,%d - %d,%d,%d\n",NetMPI::rank(),Xi[0],Xi[1],Xi[2],Xe[0],Xe[1],Xe[2]);
00256    printf("========================\n%d) Minimal Geo Mesh %g,%g,%g - %g,%g,%g\n",NetMPI::rank(),m_lp0[0],m_lp0[1],m_lp0[2],m_lp1[0],m_lp1[1],m_lp1[2]);
00257    */
00258    double dd;
00259    double dpGauss = sqrt(log(CVy*CVy + 1));
00260    double C0  = My/(exp(0.5*dpGauss*dpGauss));
00261    
00262    
00263    for (unsigned z=0;z<m_dimZ;z++)
00264      for (unsigned y=0;y<m_dimY;y++)
00265        for (unsigned x=0;x<m_dimX;x++)
00266        {
00267          file >> dd;
00268          if (x >= Xi[0] && x < Xe[0] &&
00269              y >= Xi[1] && y < Xe[1] &&
00270              z >= Xi[2] && z < Xe[2] )
00271          {
00272            m_matrix3D(z-Xi[2],y-Xi[1],x-Xi[0]) = C0*exp(dpGauss*dd);
00273          }
00274          if (file.fail())
00275            throw new Exception("FChessBoard3D: File \"%s\", dimension of the matrix does not match with the numbers supplied in the file\n",fileName.c_str());
00276        }
00277  }
00278 
00279 /*
00280    assert(lp0[0] <= lp1[0]);
00281    assert(lp0[1] <= lp1[1]);
00282    assert(lp0[2] <= lp1[2]);
00283    assert(p0[0] <= p1[0]);
00284    assert(p0[1] <= p1[1]);
00285    assert(p0[2] <= p1[2]);
00286    assert(lp0[0] >= p0[0]);
00287    assert(lp0[1] >= p0[1]);
00288    assert(lp0[2] >= p0[2]);
00289    assert(lp1[0] <= p1[0]);
00290    assert(lp1[1] <= p1[1]);
00291    assert(lp1[2] <= p1[2]);
00292 
00293 
00294    
00295    ifstream file(fileName.c_str());
00296    if (file.fail())
00297      throw new Exception("FChessBoard3D: File \"%s\" not found",fileName.c_str());
00298 
00299    char str[MAX_LINE_SIZE];
00300 
00301    //Ignore the first line
00302    file.getline(str,MAX_LINE_SIZE);
00303    
00304    //Get the second  line containing the dimension of the chess board
00305    //in the format "size_X x size_Y x size_Z"
00306    file.getline(str,MAX_LINE_SIZE);
00307    if (file.fail() || sscanf(str,"%d x %d x %d",&m_dimX,&m_dimY,&m_dimZ) != 3)
00308      throw new Exception("FChessBoard3D: File \"%s\", Second line  (\"%s\")\n\tnot in format <number> x <number> x <number>",fileName.c_str(),str);
00309 
00310 
00311 
00312    m_dX = (p1[0]-p0[0])/m_dimX;
00313    m_dY = (p1[1]-p0[1])/m_dimY;
00314    m_dZ = (p1[2]-p0[2])/m_dimZ;
00315 
00316    Xi[0] = round(lp0[0]/m_dX);
00317    Xi[1] = round(lp0[1]/m_dY);
00318    Xi[2] = round(lp0[2]/m_dZ);
00319 
00320    Xe[0] = ceil(lp1[0]/m_dX);
00321    Xe[1] = ceil(lp1[1]/m_dY);
00322    Xe[2] = ceil(lp1[2]/m_dZ);
00323 
00324    //Allocate the table to store the chess board
00325    TableIndices<3> sizes(Xe[2]-Xi[2],Xe[1]-Xi[1],Xe[0]-Xi[0]);
00326    m_matrix3D.reinit(sizes);
00327 
00328    m_lp0[0]=Xi[0]*m_dX;
00329    m_lp0[1]=Xi[1]*m_dY;
00330    m_lp0[2]=Xi[2]*m_dZ;
00331 
00332    m_lp1[0]=Xe[0]*m_dX;
00333    m_lp1[1]=Xe[1]*m_dY;
00334    m_lp1[2]=Xe[2]*m_dZ;
00335 
00336 
00337    
00338    double dd;
00339    double dpGauss = sqrt(log(CVy*CVy + 1));
00340    double C0  = My/(exp(0.5*dpGauss*dpGauss));
00341 
00342    for (unsigned z=0;z<m_dimZ;z++)
00343      for (unsigned y=0;y<m_dimY;y++)
00344        for (unsigned x=0;x<m_dimX;x++)
00345        {
00346          file >> dd;
00347          if (x >= Xi[0] && x <= Xe[0] &&
00348              y >= Xi[1] && y <= Xe[1] &&
00349              z >= Xi[2] && z <= Xe[2] )
00350          {
00351            m_matrix3D(z-Xi[2],y-Xi[1],x-Xi[0]) = C0*exp(dpGauss*dd);
00352          }
00353          if (file.fail())
00354            throw new Exception("FChessBoard3D: File \"%s\", dimension of the matrix does not match with the numbers supplied in the file\n",fileName.c_str());
00355        }
00356 */