sandbox/ariviere/tools/liner-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
    
    void lineX(char * fname,scalar s,double yp, double zp, int maxlevel){
      FILE *fpver =fopen (fname,"a"); 
      int nn = (1<<maxlevel);
      double ** field = matrix_new (1, nn, sizeof(double));
      double stp = L0/(double)nn;
      for (int i = 0; i < nn; i++){
          double xp = stp*i + X0 + stp/2.;
          Point point = locate (xp, yp,zp);
          field[0][i] = point.level >= 0 ? interpolate(s,xp,yp,zp) : nodata;
        }
      if (pid() == 0){ // master
    #if _MPI
        MPI_Reduce (MPI_IN_PLACE, field[0], nn, MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
        for (int i = 0; i < nn; i++) {
    			fprintf (fpver, "%g ", field[0][i]);
    		}
    		fputc ('\n', fpver);
        fflush (fpver);
      }
    #if _MPI
        else // slave
        MPI_Reduce (field[0], NULL, nn, MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
      matrix_free (field);
    }
    
    void lineY(char * fname,scalar s,double xp, double zp, int maxlevel){
      FILE *fpver =fopen (fname,"a"); 
      int nn = (1<<maxlevel);
      double ** field = matrix_new (1, nn, sizeof(double));
      double stp = L0/(double)nn;
      for (int i = 0; i < nn; i++)
        {
          double yp = stp*i + Y0 + stp/2.;
          Point point = locate (xp, yp,zp);
          field[0][i] = point.level >= 0 ? interpolate(s,xp,yp,zp) : nodata;
        }
      if (pid() == 0){ // master
    #if _MPI
        MPI_Reduce (MPI_IN_PLACE, field[0], nn, MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
        for (int i = 0; i < nn; i++) {
    			fprintf (fpver, "%g ", field[0][i]);
        }
    		fputc ('\n', fpver);
        fflush (fpver);
      }
    #if _MPI
      else // slave
        MPI_Reduce (field[0], NULL, nn, MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
      matrix_free (field);
    }
    
    void lineZ(char * fname,scalar s,double xp, double yp, int maxlevel){
      FILE *fpver =fopen (fname,"a"); 
      int nn = (1<<maxlevel);
      double ** field = matrix_new (1, nn, sizeof(double));
      double stp = L0/(double)nn;
      for (int i = 0; i < nn; i++)
        {
          double zp = stp*i + Z0 + stp/2.;
          Point point = locate (xp, yp,zp);
          field[0][i] = point.level >= 0 ? interpolate(s,xp,yp,zp) : nodata;
        }
      if (pid() == 0){ // master
    #if _MPI
        MPI_Reduce (MPI_IN_PLACE, field[0], nn, MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
        for (int i = 0; i < nn; i++) {
    			fprintf (fpver, "%g ", field[0][i]);
        }
        fputc ('\n', fpver);
        fflush (fpver);
      }
    #if _MPI
      else // slave
        MPI_Reduce (field[0], NULL, nn, MPI_DOUBLE, MPI_MIN, 0,
    		MPI_COMM_WORLD);
    #endif
      matrix_free (field);
    }