sandbox/ariviere/tools/slicer-interp.h

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    
    void sliceXZ(char * fname,scalar s,double yp, int maxlevel){
      FILE *fpver =fopen (fname,"w"); 
      int nn = (1<<maxlevel);
      double ** field = matrix_new (nn, nn, sizeof(double));
      double stp = L0/(double)nn;
      for (int i = 0; i < nn; i++){
          double xp = stp*i + X0 + stp/2.;
          for (int j = 0; j < nn; j++) {
        	  double zp = stp*j + Z0 + stp/2.;
        	  Point point = locate (xp, yp,zp);
        	  field[i][j] = point.level >= 0 ? interpolate(s,xp,yp,zp) : nodata;
        	}
        }
      if (pid() == 0){ // master
    #if _MPI
        MPI_Reduce (MPI_IN_PLACE, field[0], sq(nn), MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
        for (int i = 0; i < nn; i++) {
          for (int j = 0; j < nn; j++) {
    	      fprintf (fpver, "%g\t", field[i][j]);
          }
          fputc ('\n', fpver);
        }
        fflush (fpver);
      }
    #if _MPI
      else // slave
        MPI_Reduce (field[0], NULL, nn*nn, MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
      matrix_free (field);
    }
    
    void sliceYZ(char * fname,scalar s,double xp, int maxlevel){
      FILE *fpver =fopen (fname,"w"); 
      int nn = (1<<maxlevel);
      double ** field = matrix_new (nn, nn, sizeof(double));
      double stp = L0/(double)nn;
      for (int i = 0; i < nn; i++)
        {
          double yp = stp*i + Y0 + stp/2.;
          for (int j = 0; j < nn; j++) 
        	{
        	  double zp = stp*j + Z0 + stp/2.;
        	  Point point = locate (xp, yp,zp);
        	  field[i][j] = point.level >= 0 ? interpolate(s,xp,yp,zp) : nodata;
        	}
        }
      if (pid() == 0){ // master
    #if _MPI
        MPI_Reduce (MPI_IN_PLACE, field[0], sq(nn), MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
        for (int i = 0; i < nn; i++) {
          for (int j = 0; j < nn; j++) {
    	fprintf (fpver, "%g\t", field[i][j]);
          }
          fputc ('\n', fpver);
        }
        fflush (fpver);
      }
    #if _MPI
      else // slave
        MPI_Reduce (field[0], NULL, nn*nn, MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
      matrix_free (field);
    }
    
    void sliceXY(char * fname,scalar s,double zp, int maxlevel){
      FILE *fpver =fopen (fname,"w"); 
      int nn = (1<<maxlevel);
      double ** field = matrix_new (nn, nn, sizeof(double));
      double stp = L0/(double)nn;
      for (int i = 0; i < nn; i++)
        {
          double xp = stp*i + X0 + stp/2.;
          for (int j = 0; j < nn; j++) 
    	{
    	  double yp = stp*j + Y0 + stp/2.;
    	  Point point = locate (xp, yp,zp);
    	  field[i][j] = point.level >= 0 ? interpolate(s,xp,yp,zp) : nodata;
    	}
        }
      if (pid() == 0){ // master
    #if _MPI
        MPI_Reduce (MPI_IN_PLACE, field[0], sq(nn), MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
        for (int i = 0; i < nn; i++) {
          for (int j = 0; j < nn; j++) {
    	fprintf (fpver, "%g\t", field[i][j]);
          }
          fputc ('\n', fpver);
        }
        fflush (fpver);
      }
    #if _MPI
      else // slave
        MPI_Reduce (field[0], NULL, nn*nn, MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
      matrix_free (field);
    }