Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Untitled

By: a guest on May 22nd, 2012  |  syntax: C++  |  size: 5.74 KB  |  hits: 13  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1.  
  2. void GetPlane(double x1, double y1, double x2, double y2, double *coef)
  3. {
  4.     coef[0] = y2 - y1;
  5.     coef[1] = x1 - x2;
  6.     coef[2] = y1 * (x2 - x1) - x1 * (y2 - y1);
  7.  
  8.     return;
  9. }
  10.  
  11. double GetPlaneValue(double x, double y, double *coef)
  12. {
  13.     return coef[0] * x  +  coef[1] * y  +  coef[2];
  14. }
  15.  
  16. void GetBasisPlane(double x, double y, double x1, double y1, double x2, double y2, double *coef)
  17. {
  18.  
  19.     GetPlane(x1, y1, x2, y2, coef);
  20.  
  21.     double tmp = 1 / GetPlaneValue(x,y,coef);
  22.  
  23.     coef[0] *= tmp;
  24.     coef[1] *= tmp;
  25.     coef[2] *= tmp;
  26. }
  27.  
  28. void InitializeBasis( double *coef)
  29. {
  30.  
  31.  
  32.     GetBasisPlane(0, 0, 0, 1, -1, 0, coef);
  33.     GetBasisPlane(0, 0, 1, 1, 1, 0, coef + 6);
  34.     GetBasisPlane(0, 0, 0, -1, -1, -1, coef + 12);
  35.  
  36.     GetBasisPlane(0, 0, 0, 1, 1, 1, coef + 3);
  37.     GetBasisPlane(0, 0, 1, 0, 0, -1, coef + 9);
  38.     GetBasisPlane(0, 0, -1, 0, -1, -1, coef + 15);
  39.  
  40.     /*GetBasisPlane(0, 0, -x, 0, 0, y, coef);
  41.     GetBasisPlane(0, 0, x, 0, x, y, coef + 6);
  42.     GetBasisPlane(0, 0, 0, -y, -x, -y, coef + 12);
  43.  
  44.     GetBasisPlane(0, 0, 0, y, x, y, coef + 3);
  45.     GetBasisPlane(0, 0, x, 0, 0, -y, coef + 9);
  46.     GetBasisPlane(0, 0, -x, 0, -x, -y, coef + 15);*/
  47.  
  48.     for (int i=0; i<6 ;i++)
  49.         qDebug()<<coef[i*3]<<coef[1+i*3]<<coef[2+i*3];
  50. }
  51.  
  52.  
  53.  
  54.  
  55. double GetBasisFunctionValue(double x, double y, int x_center, int y_center, double *coef)
  56. {
  57.     int type = -1;
  58.   /*  if (fabs(x - x_center) > 1) return 0;
  59.     if (fabs(y - y_center) > 1) return 0;*/
  60.  
  61.     if (fabs(x - x_center) + fabs(y - y_center) < 1e-15) return 1;
  62.     if (x <= x_center)
  63.     {
  64.         if (y >= y_center) type = 0;
  65.         else if ( fabs(x_center - x) <= fabs(y_center - y) ) type = 4;
  66.             else type = 5;
  67.     }
  68.     else
  69.     {
  70.         if (y <= y_center) type = 3;
  71.         else if ( fabs(x - x_center) <= fabs(y - y_center) ) type = 1;
  72.             else type = 2;
  73.     }
  74.  
  75.  
  76.     double x_tmp = x - x_center;
  77.     double y_tmp = y - y_center;
  78.     double res = GetPlaneValue(x_tmp, y_tmp, coef + type * 3);
  79.     if (res > 1 || res < 0) res=0;
  80.    // qDebug()<<"!!!"<< x_tmp << y_tmp << type << "COEF"<< coef[type*3]<< coef[type*3 +1]<< coef[type*3 +2];
  81.     return res;
  82. }
  83. void PointSurround(int p_number, int m, int n, int *xi, int *yi)
  84. {
  85.     for (int i=0; i<6; i++)
  86.     {
  87.         xi[i] = yi[i] = -900;
  88.     }
  89.  
  90.     //  point number -> coordinates
  91.  
  92.     int x = p_number % m;
  93.     int y = n - (p_number / m) -1;
  94.  
  95.     xi[6] = x;
  96.     yi[6] = y;
  97.  
  98.     qDebug()<<"("<<x<<y<<")";
  99.     if (x != 0)
  100.     {
  101.         xi[0] = x - 1;
  102.         yi[0] = y;
  103.     }
  104.  
  105.     if (y != n-1)
  106.     {
  107.         xi[1] = x;
  108.         yi[1] = y + 1;
  109.     }
  110.  
  111.     if ((x != m-1) && (y != n-1))
  112.     {
  113.         xi[2] = x + 1;
  114.         yi[2] = y + 1;
  115.     }
  116.  
  117.     if (x != m-1)
  118.     {
  119.         xi[3] = x + 1;
  120.         yi[3] = y;
  121.     }
  122.  
  123.     if (y != 0)
  124.     {
  125.         xi[4] = x;
  126.         yi[4] = y - 1;
  127.     }
  128.  
  129.     if ((x != 0) && (y != 0))
  130.     {
  131.         xi[5] = x - 1;
  132.         yi[5] = y - 1;
  133.     }
  134.  
  135.     //*count = ind;
  136.     return;
  137. }
  138.  
  139.  
  140.  
  141.  
  142.  
  143.  
  144. void Plane3p(double* x, double *y, double *z, double *coef)
  145. {
  146.         double a[9];
  147.  
  148.         a[0] =  x[2];
  149.         a[1] =  y[2];
  150.         a[2] =  z[2];
  151.  
  152.         a[3] = x[0] - x[2];
  153.         a[4] = y[0] - y[2];
  154.         a[5] = z[0] - z[2];
  155.  
  156.         a[6] = x[1] - x[2];
  157.         a[7] = y[1] - y[2];
  158.         a[8] = z[1] - z[2];
  159.  
  160.         double tmp = 1 / (a[4] * a[6]  -  a[3] * a[7]);
  161.         coef[0] = tmp * (-a[5] * a[7]  +  a[4] * a[8] );
  162.         coef[1] = tmp * (a[5] * a[6]  -  a[3] * a[8] );
  163.         coef[2] = tmp * (a[2] * a[4] * a[6] - a[1] * a[5] * a[6] - a[2] * a[3] * a[7] + a[0] * a[5] * a[7] + a[1] * a[3] * a[8] - a[0] * a[4] * a[8]);
  164. }
  165.  
  166. double Integrate0(double a1, double b1, double c1, double a2, double b2, double c2, double x0, double y0)
  167. {
  168.         return
  169. (a1 * (-4 * (c2 + 3*c2*x0) + 2 * a2 * (1 + 4*x0 + 6*x0*x0) + b2*(-1 - 4*x0 + 4*y0 + 12 * x0 * y0)) + 2 * (2 * c1 * (b2 + 3 * c2 - 3 * b2 * y0) + b1 * (b2 + 2*c2 - 4*b2*y0 - 6 *c2 * y0 + 6* b2 * y0 * y0)) + a2 * (-4 * (c1 + 3 * c1 * x0) + b1 * (-1 + 4*y0 + 4*x0 * (-1 + 3*y0))))/24;
  170. }
  171.  
  172. double Integrate1(double a1, double b1, double c1, double a2, double b2, double c2, double x0, double y0)
  173. {
  174.         return
  175. (a1*(4*c2*(1 - 3* x0) + 2* a2* (1 - 4* x0 + 6* x0*x0) +
  176.       b2* (3 - 8* x0 - 4* y0 + 12* x0* y0)) +
  177.          a2*(4* c1* (1 - 3* x0) + b1* (3 - 4* y0 + 4* x0* (-2 + 3* y0))) +
  178.             2* (4* b2* c1 + 6* c1* c2 - 6* b2* c1* y0 +
  179.                   b1* (4* c2 - 6* c2* y0 + b2* (3 - 8* y0 + 6* y0*y0))))/24;
  180. }
  181.  
  182. void init_f(int m, int n, double *coef, double *f) // Initialize right side
  183. {
  184.  
  185.         int p_x[7]={-1};        // point surroundings x
  186.         int p_y[7]={-1};        // point surroundings y
  187.         double x[3], y[3], z[3];// point for plane
  188.         double fcoef[3];        // approximated plane coefficients
  189.         double x_val, y_val;
  190.         for (int i1=0; i1<n; i1++)
  191.         {
  192.                 for(int i2=0; i2<m; i2++)
  193.                 {
  194.                         PointSurround(m*i1+i2, m, n, p_x, p_y);
  195.                         for (int z=0; z<7; z++)
  196.                             qDebug()<<(m*i1+i2)<<"COORD:"<<p_x[z]<<p_y[z];
  197.  
  198.  
  199.                         qDebug("------");
  200.                         for (int j=0; j<6; j++)
  201.                         {
  202.                                 if ( (p_x[j] >= 0) && (p_x[ (j+1)%6 ] >= 0) )
  203.                                 {
  204.                                         //x_val = x[0]
  205.                                         x[0] = p_x[j] - p_x[6];
  206.                                         y[0] = p_y[j] - p_y[6];
  207.                                         z[0] = func( p_x[j], p_y[j] );
  208.  
  209.                                         x[1] = p_x[(j+1)%6] - p_x[6];
  210.                                         y[1] = p_y[(j+1)%6] - p_y[6];
  211.                                         z[1] = func( p_x[(j+1)%6], p_y[(j+1)%6] );
  212.  
  213.                                         x[2] = 0;
  214.                                         y[2] = 0;
  215.                                         z[2] = func(i1, i2);
  216.  
  217.                                         Plane3p(x, y, z, fcoef);
  218.  
  219.                                         double x0 = 0;
  220.                                         double y0 = 0;
  221.  
  222.                                         if (j % 2 ==0)
  223.                                         {
  224.                                                 if (j == 2)     x0 = -1;
  225.                                                 if (j == 4)     y0 = 1;
  226.  
  227.                                                 f[i2 + i1 * m] += 24 * Integrate0(fcoef[0], fcoef[1], fcoef[2], coef[j*3], coef[j*3 + 1], coef[j*3 + 2], x0, y0);
  228.                                         }
  229.                                         else
  230.                                         {
  231.                                                 if (j == 3)     y0 = 1;
  232.                                                 if (j == 5)
  233.                                                 {
  234.                                                         x0 = 1;
  235.                                                         y0 = 1;
  236.                                                 }
  237.  
  238.                                                 f[i2 + i1 * m] += 24 * Integrate1(fcoef[0], fcoef[1], fcoef[2], coef[j*3], coef[j*3 + 1], coef[j*3 + 2], x0, y0);
  239.                                         }
  240.  
  241.                                 }
  242.                         }
  243.                 }
  244.         }
  245.  
  246.         for (int i=0; i<m*n;i++)
  247.             qDebug()<<i<<f[i];
  248. }
create a new version of this paste RAW Paste Data