/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4660 by sshaw, Fri Feb 7 03:08:49 2014 UTC revision 4675 by sshaw, Mon Feb 17 01:07:12 2014 UTC
# Line 490  void Rectangle::dump(const string& fileN Line 490  void Rectangle::dump(const string& fileN
490          fn+=".silo";          fn+=".silo";
491      }      }
492    
493      int driver=DB_HDF5;          int driver=DB_HDF5;
494      DBfile* dbfile = NULL;      DBfile* dbfile = NULL;
495      const char* blockDirFmt = "/block%04d";      const char* blockDirFmt = "/block%04d";
496    
# Line 1870  void Rectangle::interpolateNodesOnFaces( Line 1870  void Rectangle::interpolateNodesOnFaces(
1870  namespace  namespace
1871  {  {
1872      // Calculates a guassian blur colvolution matrix for 2D      // Calculates a guassian blur colvolution matrix for 2D
1873      // See wiki article on the subject          // See wiki article on the subject
1874      double* get2DGauss(unsigned radius, double sigma)      double* get2DGauss(unsigned radius, double sigma)
1875      {      {
1876          double* arr=new double[(radius*2+1)*(radius*2+1)];          double* arr=new double[(radius*2+1)*(radius*2+1)];
1877          double common=M_1_PI*0.5*1/(sigma*sigma);          double common=M_1_PI*0.5*1/(sigma*sigma);
1878      double total=0;          double total=0;
1879      int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
1880      for (int y=-r;y<=r;++y)          for (int y=-r;y<=r;++y)
1881      {          {
1882          for (int x=-r;x<=r;++x)              for (int x=-r;x<=r;++x)
1883          {                      {
1884              arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));                  arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));
1885  // cout << (x+y*(r*2+1)) << " " << arr[(x+r)+(y+r)*(r*2+1)] << endl;  // cout << (x+y*(r*2+1)) << " " << arr[(x+r)+(y+r)*(r*2+1)] << endl;
1886              total+=arr[(x+r)+(y+r)*(r*2+1)];                  total+=arr[(x+r)+(y+r)*(r*2+1)];
1887          }              }
1888      }          }
1889      double invtotal=1/total;          double invtotal=1/total;
1890  //cout << "Inv total is "    << invtotal << endl;  //cout << "Inv total is "         << invtotal << endl;
1891      for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)          for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)
1892      {          {
1893          arr[p]*=invtotal;              arr[p]*=invtotal;
1894  //cout << p << "->" << arr[p] << endl;        //cout << p << "->" << arr[p] << endl;
1895      }          }
1896      return arr;          return arr;
1897      }      }
1898        
1899      // applies conv to source to get a point.      // applies conv to source to get a point.
1900      // (xp, yp) are the coords in the source matrix not the destination matrix      // (xp, yp) are the coords in the source matrix not the destination matrix
1901      double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)      double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)
1902      {      {
1903          size_t bx=xp-radius, by=yp-radius;          size_t bx=xp-radius, by=yp-radius;
1904      size_t sbase=bx+by*width;          size_t sbase=bx+by*width;
1905      double result=0;          double result=0;
1906      for (int y=0;y<2*radius+1;++y)          for (int y=0;y<2*radius+1;++y)
1907      {              {
1908          for (int x=0;x<2*radius+1;++x)              for (int x=0;x<2*radius+1;++x)
1909          {              {
1910              result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];                  result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];
1911          }              }
1912      }          }
1913          return result;                return result;
1914      }      }
1915  }  }
1916    
# Line 1922  A parameter radius  gives the size of th Line 1922  A parameter radius  gives the size of th
1922  A point on the left hand edge for example, will still require `radius` extra points to the left  A point on the left hand edge for example, will still require `radius` extra points to the left
1923  in order to complete the stencil.  in order to complete the stencil.
1924    
1925  All local calculation is done on an array called `src`, with  All local calculation is done on an array called `src`, with
1926  dimensions = ext[0] * ext[1].  dimensions = ext[0] * ext[1].
1927  Where ext[i]= internal[i]+2*radius.  Where ext[i]= internal[i]+2*radius.
1928    
1929  Now for MPI there is overlap to deal with. We need to share both the overlapping  Now for MPI there is overlap to deal with. We need to share both the overlapping
1930  values themselves but also the external region.  values themselves but also the external region.
1931    
1932  In a hypothetical 1-D case:  In a hypothetical 1-D case:
# Line 1945  need to be known. Line 1945  need to be known.
1945    
1946  pp123(4)56   23(4)567pp  pp123(4)56   23(4)567pp
1947    
1948  Now in our case, we wout set all the values 23456 on the left rank and send them to the  Now in our case, we wout set all the values 23456 on the left rank and send them to the
1949  right hand rank.  right hand rank.
1950    
1951  So the edges _may_ need to be shared at a distance `inset` from all boundaries.  So the edges _may_ need to be shared at a distance `inset` from all boundaries.
1952  inset=2*radius+1      inset=2*radius+1
1953  This is to ensure that values at distance `radius` from the shared/overlapped element  This is to ensure that values at distance `radius` from the shared/overlapped element
1954  that ripley has.  that ripley has.
1955    
# Line 1965  escript::Data Rectangle::randomFill(long Line 1965  escript::Data Rectangle::randomFill(long
1965          throw RipleyException("Unsupported random filter");          throw RipleyException("Unsupported random filter");
1966      }      }
1967      boost::python::extract<string> ex(filter[0]);      boost::python::extract<string> ex(filter[0]);
1968      if (!ex.check() || (ex()!="gaussian"))      if (!ex.check() || (ex()!="gaussian"))
1969      {      {
1970          throw RipleyException("Unsupported random filter");          throw RipleyException("Unsupported random filter");
1971      }      }
# Line 1980  escript::Data Rectangle::randomFill(long Line 1980  escript::Data Rectangle::randomFill(long
1980      if (!ex2.check() || (sigma=ex2())<=0)      if (!ex2.check() || (sigma=ex2())<=0)
1981      {      {
1982          throw RipleyException("Sigma must be a postive floating point number.");          throw RipleyException("Sigma must be a postive floating point number.");
1983      }          }
1984        
1985      size_t internal[2];      size_t internal[2];
1986      internal[0]=m_NE[0]+1;  // number of points in the internal region      internal[0]=m_NE[0]+1;        // number of points in the internal region
1987      internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of      internal[1]=m_NE[1]+1;        // that is, the ones we need smoothed versions of
1988      size_t ext[2];      size_t ext[2];
1989      ext[0]=internal[0]+2*radius;    // includes points we need as input      ext[0]=internal[0]+2*radius;        // includes points we need as input
1990      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=internal[1]+2*radius;        // for smoothing
1991        
1992      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
1993      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
1994    
1995      if ((2*radius>=internal[0]) || (2*radius>=internal[1]))      if ((2*radius>=internal[0]) || (2*radius>=internal[1]))
1996      {      {
1997          throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");          throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");
1998      }      }
1999        
2000      
2001      double* src=new double[ext[0]*ext[1]];      double* src=new double[ext[0]*ext[1]];
2002      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]);        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]);
2003        
2004      
2005  #ifdef ESYS_MPI      #ifdef ESYS_MPI
2006      size_t inset=2*radius+1;      size_t inset=2*radius+1;
2007      size_t Eheight=ext[1]-2*inset;  // how high the E (shared) region is      size_t Eheight=ext[1]-2*inset;        // how high the E (shared) region is
2008      size_t Swidth=ext[0]-2*inset;      size_t Swidth=ext[0]-2*inset;
2009    
2010      double* SWin=new double[inset*inset];  memset(SWin, 0, inset*inset*sizeof(double));      double* SWin=new double[inset*inset];  memset(SWin, 0, inset*inset*sizeof(double));
# Line 2017  escript::Data Rectangle::randomFill(long Line 2017  escript::Data Rectangle::randomFill(long
2017      unsigned int base=ext[0]-inset+(ext[1]-inset)*ext[0];      unsigned int base=ext[0]-inset+(ext[1]-inset)*ext[0];
2018      for (unsigned int i=0;i<inset;++i)      for (unsigned int i=0;i<inset;++i)
2019      {      {
2020      memcpy(NEout+inset*i, src+base, inset*sizeof(double));          memcpy(NEout+inset*i, src+base, inset*sizeof(double));
2021      base+=ext[0];          base+=ext[0];
2022      }      }
2023      double* NWout=new double[inset*inset];  memset(NWout, 0, inset*inset*sizeof(double));      double* NWout=new double[inset*inset];  memset(NWout, 0, inset*inset*sizeof(double));
2024      base=(ext[1]-inset)*ext[0];      base=(ext[1]-inset)*ext[0];
2025      for (unsigned int i=0;i<inset;++i)      for (unsigned int i=0;i<inset;++i)
2026      {      {
2027      memcpy(NWout+inset*i, src+base, inset*sizeof(double));          memcpy(NWout+inset*i, src+base, inset*sizeof(double));
2028      base+=ext[0];          base+=ext[0];
2029      }      }
2030        
2031      double* SEout=new double[inset*inset];  memset(SEout, 0, inset*inset*sizeof(double));      double* SEout=new double[inset*inset];  memset(SEout, 0, inset*inset*sizeof(double));
2032      base=ext[0]-inset;      base=ext[0]-inset;
2033      for (int i=0;i<inset;++i)      for (int i=0;i<inset;++i)
2034      {      {
2035      memcpy(SEout+inset*i, src+base, inset*sizeof(double));          memcpy(SEout+inset*i, src+base, inset*sizeof(double));
2036      base+=ext[0];          base+=ext[0];
2037      }      }
2038      double* Nout=new double[inset*Swidth];  memset(Nout, 0, inset*Swidth*sizeof(double));      double* Nout=new double[inset*Swidth];  memset(Nout, 0, inset*Swidth*sizeof(double));
2039      base=inset+(ext[1]-inset)*ext[0];      base=inset+(ext[1]-inset)*ext[0];
2040      for (unsigned int i=0;i<inset;++i)      for (unsigned int i=0;i<inset;++i)
2041      {      {
2042      memcpy(Nout+Swidth*i, src+base, Swidth*sizeof(double));          memcpy(Nout+Swidth*i, src+base, Swidth*sizeof(double));
2043      base+=ext[0];          base+=ext[0];
2044      }      }
2045        
2046      double* Eout=new double[inset*Eheight];  memset(Eout, 0, inset*Eheight*sizeof(double));      double* Eout=new double[inset*Eheight];  memset(Eout, 0, inset*Eheight*sizeof(double));
2047      base=ext[0]-inset+inset*ext[0];      base=ext[0]-inset+inset*ext[0];
2048      for (unsigned int i=0;i<Eheight;++i)      for (unsigned int i=0;i<Eheight;++i)
2049      {      {
2050      memcpy(Eout+i*inset, src+base, inset*sizeof(double));          memcpy(Eout+i*inset, src+base, inset*sizeof(double));
2051      base+=ext[0];          base+=ext[0];
2052      }        }
2053    
2054      MPI_Request reqs[10];      MPI_Request reqs[10];
2055      MPI_Status stats[10];      MPI_Status stats[10];
2056      short rused=0;      short rused=0;
2057        
2058      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
2059      dim_t Y=m_mpiInfo->rank/m_NX[0];      dim_t Y=m_mpiInfo->rank/m_NX[0];
2060      dim_t row=m_NX[0];      dim_t row=m_NX[0];
2061      bool swused=false;      // These vars will be true if data needs to be copied out of them      bool swused=false;                // These vars will be true if data needs to be copied out of them
2062      bool seused=false;      bool seused=false;
2063      bool nwused=false;      bool nwused=false;
2064      bool sused=false;      bool sused=false;
2065      bool wused=false;          bool wused=false;
2066        
2067      // Tags:      // Tags:
2068      // 10 : EW transfer (middle)      // 10 : EW transfer (middle)
2069      // 8 : NS transfer (middle)      // 8 : NS transfer (middle)
2070      // 7 : NE corner -> to N, E and NE      // 7 : NE corner -> to N, E and NE
2071      // 11 : NW corner to SW corner (only used on the left hand edge      // 11 : NW corner to SW corner (only used on the left hand edge
2072      // 12 : SE corner to SW corner (only used on the bottom edge      // 12 : SE corner to SW corner (only used on the bottom edge
2073        
2074    
2075    
2076      int comserr=0;      int comserr=0;
2077      if (Y!=0)   // not on bottom row,      if (Y!=0)        // not on bottom row,
2078      {      {
2079      if (X!=0)   // not on the left hand edge          if (X!=0)        // not on the left hand edge
2080      {          {
2081          // recv bottomleft from SW              // recv bottomleft from SW
2082          comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, (X-1)+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));              comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, (X-1)+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));
2083          swused=true;              swused=true;
2084          comserr|=MPI_Irecv(Win, Eheight*inset, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));              comserr|=MPI_Irecv(Win, Eheight*inset, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));
2085          wused=true;              wused=true;
2086      }          }
2087      else    // on the left hand edge          else        // on the left hand edge
2088      {          {
2089          comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, (Y-1)*row, 11, m_mpiInfo->comm, reqs+(rused++));              comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, (Y-1)*row, 11, m_mpiInfo->comm, reqs+(rused++));
2090          swused=true;              swused=true;
2091      }          }
2092      comserr|=MPI_Irecv(Sin, Swidth*inset, MPI_DOUBLE, X+(Y-1)*row, 8, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Irecv(Sin, Swidth*inset, MPI_DOUBLE, X+(Y-1)*row, 8, m_mpiInfo->comm, reqs+(rused++));
2093      sused=true;          sused=true;
2094      comserr|=MPI_Irecv(SEin, inset*inset, MPI_DOUBLE, X+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Irecv(SEin, inset*inset, MPI_DOUBLE, X+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));
2095      seused=true;          seused=true;
2096    
2097          
     }  
     else        // on the bottom row  
     {  
     if (X!=0)  
     {  
         comserr|=MPI_Irecv(Win, Eheight*inset, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));  
         wused=true;  
         // Need to use tag 12 here because SW is coming from the East not South East  
         comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, X-1+Y*row, 12, m_mpiInfo->comm, reqs+(rused++));  
         swused=true;  
     }  
     if (X!=(row-1))  
     {  
         comserr|=MPI_Isend(SEout, inset*inset, MPI_DOUBLE, X+1+(Y)*row, 12, m_mpiInfo->comm, reqs+(rused++));    
     }  
     }  
       
     if (Y!=(m_NX[1]-1)) // not on the top row  
     {  
     comserr|=MPI_Isend(Nout, inset*Swidth, MPI_DOUBLE, X+(Y+1)*row, 8, m_mpiInfo->comm, reqs+(rused++));  
     comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));  
     if (X!=(row-1)) // not on right hand edge  
     {  
         comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+1+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));  
     }  
     if (X==0)   // left hand edge  
     {  
         comserr|=MPI_Isend(NWout, inset*inset, MPI_DOUBLE, (Y+1)*row,11, m_mpiInfo->comm, reqs+(rused++));        
     }    
2098      }      }
2099      if (X!=(row-1)) // not on right hand edge      else                // on the bottom row
2100      {      {
2101      comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+1+(Y)*row, 7, m_mpiInfo->comm, reqs+(rused++));          if (X!=0)
2102      comserr|=MPI_Isend(Eout, Eheight*inset, MPI_DOUBLE, X+1+(Y)*row, 10, m_mpiInfo->comm, reqs+(rused++));          {
2103                comserr|=MPI_Irecv(Win, Eheight*inset, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));
2104                wused=true;
2105                // Need to use tag 12 here because SW is coming from the East not South East
2106                comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, X-1+Y*row, 12, m_mpiInfo->comm, reqs+(rused++));
2107                swused=true;
2108            }
2109            if (X!=(row-1))
2110            {
2111                comserr|=MPI_Isend(SEout, inset*inset, MPI_DOUBLE, X+1+(Y)*row, 12, m_mpiInfo->comm, reqs+(rused++));
2112            }
2113        }
2114    
2115        if (Y!=(m_NX[1]-1))        // not on the top row
2116        {
2117             comserr|=MPI_Isend(Nout, inset*Swidth, MPI_DOUBLE, X+(Y+1)*row, 8, m_mpiInfo->comm, reqs+(rused++));
2118            comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));
2119            if (X!=(row-1))        // not on right hand edge
2120            {
2121                comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+1+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));
2122            }
2123            if (X==0)        // left hand edge
2124            {
2125                comserr|=MPI_Isend(NWout, inset*inset, MPI_DOUBLE, (Y+1)*row,11, m_mpiInfo->comm, reqs+(rused++));
2126            }
2127        }
2128        if (X!=(row-1))        // not on right hand edge
2129        {
2130            comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+1+(Y)*row, 7, m_mpiInfo->comm, reqs+(rused++));
2131            comserr|=MPI_Isend(Eout, Eheight*inset, MPI_DOUBLE, X+1+(Y)*row, 10, m_mpiInfo->comm, reqs+(rused++));
2132      }      }
2133      if (X!=0)      if (X!=0)
2134      {      {
2135      comserr|=MPI_Irecv(NWin, inset*inset, MPI_DOUBLE, (X-1)+Y*row, 7, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Irecv(NWin, inset*inset, MPI_DOUBLE, (X-1)+Y*row, 7, m_mpiInfo->comm, reqs+(rused++));
2136      nwused=true;          nwused=true;
2137          
2138          
2139      }      }
2140        
2141      if (!comserr)      if (!comserr)
2142      {          {
2143          comserr=MPI_Waitall(rused, reqs, stats);          comserr=MPI_Waitall(rused, reqs, stats);
2144      }      }
2145    
2146      if (comserr)      if (comserr)
2147      {      {
2148      // Yes this is throwing an exception as a result of an MPI error.          // Yes this is throwing an exception as a result of an MPI error.
2149      // and no we don't inform the other ranks that we are doing this.          // and no we don't inform the other ranks that we are doing this.
2150      // however, we have no reason to believe coms work at this point anyway          // however, we have no reason to believe coms work at this point anyway
2151          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");
2152      }      }
2153        
2154    
2155      // Need to copy the values back from the buffers      // Need to copy the values back from the buffers
2156      // Copy from SW      // Copy from SW
2157        
2158      if (swused)      if (swused)
2159      {      {
2160      base=0;          base=0;
2161      for (unsigned int i=0;i<inset;++i)          for (unsigned int i=0;i<inset;++i)
2162      {          {
2163          memcpy(src+base, SWin+i*inset, inset*sizeof(double));              memcpy(src+base, SWin+i*inset, inset*sizeof(double));
2164          base+=ext[0];              base+=ext[0];
2165      }          }
2166      }      }
2167      if (seused)      if (seused)
2168      {      {
2169          base=ext[0]-inset;          base=ext[0]-inset;
2170      for (unsigned int i=0;i<inset;++i)          for (unsigned int i=0;i<inset;++i)
2171      {          {
2172          memcpy(src+base, SEin+i*inset, inset*sizeof(double));              memcpy(src+base, SEin+i*inset, inset*sizeof(double));
2173          base+=ext[0];              base+=ext[0];
2174      }          }
2175      }      }
2176      if (nwused)      if (nwused)
2177      {      {
2178          base=(ext[1]-inset)*ext[0];          base=(ext[1]-inset)*ext[0];
2179      for (unsigned int i=0;i<inset;++i)          for (unsigned int i=0;i<inset;++i)
2180      {          {
2181          memcpy(src+base, NWin+i*inset, inset*sizeof(double));              memcpy(src+base, NWin+i*inset, inset*sizeof(double));
2182          base+=ext[0];              base+=ext[0];
2183      }          }
2184      }      }
2185      if (sused)      if (sused)
2186      {      {
2187         base=inset;         base=inset;
2188         for (unsigned int i=0;i<inset;++i)         for (unsigned int i=0;i<inset;++i)
2189         {         {
2190         memcpy(src+base, Sin+i*Swidth, Swidth*sizeof(double));             memcpy(src+base, Sin+i*Swidth, Swidth*sizeof(double));
2191         base+=ext[0];             base+=ext[0];
2192         }         }
2193      }      }
2194      if (wused)      if (wused)
2195      {      {
2196      base=inset*ext[0];          base=inset*ext[0];
2197      for (unsigned int i=0;i<Eheight;++i)          for (unsigned int i=0;i<Eheight;++i)
2198      {          {
2199          memcpy(src+base, Win+i*inset, inset*sizeof(double));              memcpy(src+base, Win+i*inset, inset*sizeof(double));
2200          base+=ext[0];              base+=ext[0];
2201      }          }
2202          
2203      }      }
2204        
2205      delete[] SWin;      delete[] SWin;
2206      delete[] SEin;      delete[] SEin;
2207      delete[] NWin;      delete[] NWin;
# Line 2213  escript::Data Rectangle::randomFill(long Line 2213  escript::Data Rectangle::randomFill(long
2213      delete[] SEout;      delete[] SEout;
2214      delete[] Nout;      delete[] Nout;
2215      delete[] Eout;      delete[] Eout;
2216  #endif      #endif
2217      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
2218      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
2219      // don't need to check for exwrite because we just made it      // don't need to check for exwrite because we just made it
2220      escript::DataVector& dv=resdat.getExpandedVectorReference();      escript::DataVector& dv=resdat.getExpandedVectorReference();
2221      double* convolution=get2DGauss(radius, sigma);      double* convolution=get2DGauss(radius, sigma);
2222      for (size_t y=0;y<(internal[1]);++y)          for (size_t y=0;y<(internal[1]);++y)
2223      {      {
2224          for (size_t x=0;x<(internal[0]);++x)          for (size_t x=0;x<(internal[0]);++x)
2225      {              {
2226          dv[x+y*(internal[0])]=Convolve2D(convolution, src, x+radius, y+radius, radius, ext[0]);              dv[x+y*(internal[0])]=Convolve2D(convolution, src, x+radius, y+radius, radius, ext[0]);
2227            
2228      }          }
2229      }      }
2230      delete[] convolution;      delete[] convolution;
2231      delete[] src;      delete[] src;
# Line 2264  int Rectangle::findNode(const double *co Line 2264  int Rectangle::findNode(const double *co
2264              double ydist = (y - (ey + dy)*m_dx[1]);              double ydist = (y - (ey + dy)*m_dx[1]);
2265              double total = xdist*xdist + ydist*ydist;              double total = xdist*xdist + ydist*ydist;
2266              if (total < minDist) {              if (total < minDist) {
2267                  closest = INDEX2(ex+dx-m_offset[0], ey+dy-m_offset[1], m_NE[0] + 1);                  closest = INDEX2(ex+dx-m_offset[0], ey+dy-m_offset[1], m_NE[0] + 1);
2268                  minDist = total;                  minDist = total;
2269              }              }
2270          }          }

Legend:
Removed from v.4660  
changed lines
  Added in v.4675

  ViewVC Help
Powered by ViewVC 1.1.26