hstep_file.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 Werner Scholz
    Copyright (C) 2008 Stefan Tibus

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: hstep_file.c 1000 2008-07-31 16:37:20Z tibus $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "field.h"

#define HSTEPCOL 2

/*
input file structure:
-----------
n
x_0  y_0
x_1  y_1
x_2  y_2
x_3  y_3
...
x_n  y_n
-----------
*/


int read_hstep_file(FILE *fd,PetscReal **indata,const int col)
{
  int            nrows;
  PetscReal      *xydata;

  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  if (col<2) {
    SETERRQ1(PETSC_ERR_ARG_CORRUPT,"col=%i < 2!\n",col);
  }

  /* only the first processor reads the file */
  if (!rank) {
    char msg[256];

    /* get number of lines */
    fscanf(fd,"%i\n",&nrows);
    PetscPrintf(PETSC_COMM_WORLD,"Reading column %i, %i lines\n",col,nrows);

    /* account for additional first and last entry */
    nrows += 2;

    ierr = PetscMalloc(HSTEPCOL*nrows*sizeof(PetscReal),&xydata);CHKERRQ(ierr);

    /* get xydata from file (leave room for first ) */
    for (int i=1; i<nrows-1; i++) {
      int k;
      /* read x from first column */
#ifdef PETSC_USE_SINGLE
      k=fscanf(fd,"%f",xydata+i*HSTEPCOL+0);
#else
      k=fscanf(fd,"%lf",xydata+i*HSTEPCOL+0);
#endif
      if (k!=1) {
        SETERRQ2(PETSC_ERR_FILE_UNEXPECTED,
          "xydata in corrupt in row %i col %i\n",
          i,1
        );
      }

      /* read y from column col */
      for (int j=1; j<col; j++) {
#ifdef PETSC_USE_SINGLE
        k=fscanf(fd,"%f",xydata+i*HSTEPCOL+1);
#else
        k=fscanf(fd,"%lf",xydata+i*HSTEPCOL+1);
#endif

        if (k!=1) {
          SETERRQ2(PETSC_ERR_FILE_UNEXPECTED,
            "xydata in corrupt in row %i col %i\n",
            i,j
          );
        }
      }
      /* read remainder of this line */
      fgets(msg,256,fd);
    }

    /* set first xydata set */
    xydata[0]=PETSC_MIN;
    xydata[1]=0.0;

    /* set last xydata set */
    xydata[HSTEPCOL*(nrows-1)+0]=PETSC_MAX;
    xydata[HSTEPCOL*(nrows-1)+1]=0.0;

    PetscPrintf(PETSC_COMM_WORLD,
      "%14s %14s\n",
      "step","scale"
    );

    for (int i=0; i<nrows; i++) {
      PetscPrintf(PETSC_COMM_WORLD,
        "%14e %14e\n",
        xydata[HSTEPCOL*i+0],
        xydata[HSTEPCOL*i+1]
      );
    }
  }

  /* broadcast xydata to all processors */
  ierr = MPI_Bcast(&nrows,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  if (rank) {
    ierr = PetscMalloc(HSTEPCOL*nrows*sizeof(PetscReal),&xydata);CHKERRQ(ierr);
  }
  ierr = MPI_Bcast(xydata,nrows*HSTEPCOL,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  *indata=xydata;

  MagparFunctionLogReturn(0);
}

PetscReal get_scalf_hstep(PetscReal xydata[],PetscReal xval)
{
  int       tint;
  PetscReal result;

  tint=0;
  while (xval > xydata[HSTEPCOL*(tint+1)]) tint++;

  /* calculate scaling factor */
  result=xydata[HSTEPCOL*tint+1];

  return(result);
}