writedataavs.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2009 Werner Scholz

    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: writedataavs.c 2681 2009-07-31 04:30:53Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "magpario.h"
#include "util/util.h"

#include <errno.h>

#define NEDATA  0
#define NEDATA2 0

#if defined(ZLIB) && !defined(_WIN32) && !defined(__CYGWIN32__)
EXTERN_C_BEGIN
#include "zlib.h"
EXTERN_C_END

#define FPRINTF2(a,b)   gzprintf(a,b)
#define FPRINTF3(a,b,c) gzprintf(a,b,c)
#else
#define FPRINTF2(a,b)   PetscFPrintf(PETSC_COMM_WORLD,a,b)
#define FPRINTF3(a,b,c) PetscFPrintf(PETSC_COMM_WORLD,a,b,c)
#endif

#define D_FLT 50
#define D_EID 8

/* Writes the grid in AVS format */

int WriteAVS(GridData *gdata)
{
  MagparFunctionInfoBegin;

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

#ifdef EXCH
#define nvdata 15
#else
#define nvdata 12
#endif

  int nmax;
  /* add 1 so we allocate at least once char even if gdata->ln_ele==0 */
  nmax=PetscMax(gdata->ln_ele,gdata->ln_vert)*(D_EID+nvdata*D_FLT)+1;

  char *wbufs,*wbufe;
  ierr = PetscMalloc(nmax*sizeof(char),&wbufs);CHKERRQ(ierr);
  wbufe=wbufs;

  FILE *fd;
  fd=NULL;

  char fmesh[256];
#if defined(ZLIB) && !defined(_WIN32) && !defined(__CYGWIN32__)
  ierr = PetscSNPrintf(fmesh,255,"%s.%04i.gz",gdata->simname,gdata->inp);CHKERRQ(ierr);
  if (!rank) {
    fd = (FILE*) gzopen(fmesh, "wb");
    if (fd == NULL) {
      SETERRQ2(PETSC_ERR_FILE_OPEN,"Cannot open %s error %i\n",fmesh,errno);
    }
  }
#else
  ierr = PetscSNPrintf(fmesh,255,"%s.%04i.inp",gdata->simname,gdata->inp);CHKERRQ(ierr);
  ierr = PetscFOpen(PETSC_COMM_WORLD,fmesh,"w",&fd);CHKERRQ(ierr);
  if (!rank && !fd) {
    SETERRQ2(PETSC_ERR_FILE_OPEN,"Cannot open %s error %i\n",fmesh,ierr);
  }
#endif

  PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",fmesh);

#if !defined(ZLIB) || defined(_WIN32) || defined(__CYGWIN32__)
    ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"%i %i %i %i %i\n",
      gdata->n_vert,gdata->n_ele,nvdata,NEDATA,NEDATA2);CHKERRQ(ierr);

    for (int i=0; i<gdata->n_vert; i++) {
      ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
        "%i %g %g %g\n",
        i+1,
        PetscRealPart(gdata->vertxyz[ND*i+0]),
        PetscRealPart(gdata->vertxyz[ND*i+1]),
        PetscRealPart(gdata->vertxyz[ND*i+2])
      );CHKERRQ(ierr);
    }

    for (int i=0; i<gdata->ln_ele; i++) {
      int nchars;
      nchars=
      snprintf(
        wbufe,
        nmax-(wbufe-wbufs),
        "%i %i tet %i %i %i %i\n",
        gdata->elel2g[i]+1,
        gdata->eleprop[i]+1,
        gdata->elevert[NV*i+0]+1,
        gdata->elevert[NV*i+1]+1,
        gdata->elevert[NV*i+2]+1,
        gdata->elevert[NV*i+3]+1
      );
      /* check for problems of snprintf */
      assert(nchars>0);

      /* update pointer to the end of the array */
      wbufe += nchars;

      /* check if we write beyond allocated area */
      assert(wbufe-wbufs < nmax);
    }
    SynchronizedFastFPrintf(fd,wbufe-wbufs,nmax,wbufs,PETSC_FALSE);
    wbufe=wbufs;
#endif

  if (!rank) {
    /* must not use CHKERRQ(ierr); here in order not to confuse gzprintf */
    ierr = FPRINTF3(fd,"%i",nvdata);
    for (int i=0;i<nvdata;i++) {
      ierr = FPRINTF2(fd," 1");
    }
    ierr = FPRINTF2(fd,"\n");

    ierr = FPRINTF2(fd,"M_x, none\n");
    ierr = FPRINTF2(fd,"M_y, none\n");
    ierr = FPRINTF2(fd,"M_z, none\n");
    ierr = FPRINTF2(fd,"Hdemag_x, none\n");
    ierr = FPRINTF2(fd,"Hdemag_y, none\n");
    ierr = FPRINTF2(fd,"Hdemag_z, none\n");
#ifdef EXCH
    ierr = FPRINTF2(fd,"Hani_x, none\n");
    ierr = FPRINTF2(fd,"Hani_y, none\n");
    ierr = FPRINTF2(fd,"Hani_z, none\n");
    ierr = FPRINTF2(fd,"Hexch_x, none\n");
    ierr = FPRINTF2(fd,"Hexch_y, none\n");
    ierr = FPRINTF2(fd,"Hexch_z, none\n");
#else
    ierr = FPRINTF2(fd,"Hexchani_x, none\n");
    ierr = FPRINTF2(fd,"Hexchani_y, none\n");
    ierr = FPRINTF2(fd,"Hexchani_z, none\n");
#endif
    ierr = FPRINTF2(fd,"Hext_x, none\n");
    ierr = FPRINTF2(fd,"Hext_y, none\n");
    ierr = FPRINTF2(fd,"Hext_z, none\n");
  }

  PetscReal *ta_M;
  VecGetArray(gdata->M,&ta_M);

  PetscReal *ta_VHdem;
  if (gdata->VHdem!=PETSC_NULL) {
    VecGetArray(gdata->VHdem,&ta_VHdem);
  }
  else {
    ta_VHdem=PETSC_NULL;
  }

#ifdef EXCH
  PetscReal *ta_VHexchonly;
  ierr = VecAXPY(gdata->VHexchani,-1.0,gdata->VHexchonly);CHKERRQ(ierr);
  VecGetArray(gdata->VHexchonly,&ta_VHexchonly);
#endif

  PetscReal *ta_VHexchani;
  VecGetArray(gdata->VHexchani,&ta_VHexchani);

  PetscReal *ta_VHext;
  VecGetArray(gdata->VHext,&ta_VHext);

  for (int i=0; i<gdata->ln_vert; i++) {
    int nchars;
    /* write into buffer */
    nchars=
    snprintf(
      wbufe,
      nmax-(wbufe-wbufs),
      "%i %g %g %g %g %g %g %g %g %g "
#ifdef EXCH
      "%g %g %g "
#endif
      "%g %g %g\n",
      gdata->vertl2g[i]+1,
      ta_M[ND*i+0],
      ta_M[ND*i+1],
      ta_M[ND*i+2],
      gdata->VHdem!=PETSC_NULL ? ta_VHdem[ND*i+0]*gdata->hscale : 0.0,
      gdata->VHdem!=PETSC_NULL ? ta_VHdem[ND*i+1]*gdata->hscale : 0.0,
      gdata->VHdem!=PETSC_NULL ? ta_VHdem[ND*i+2]*gdata->hscale : 0.0,
      ta_VHexchani[ND*i+0]*gdata->hscale,
      ta_VHexchani[ND*i+1]*gdata->hscale,
      ta_VHexchani[ND*i+2]*gdata->hscale,
#ifdef EXCH
      ta_VHexchonly[ND*i+0]*gdata->hscale,
      ta_VHexchonly[ND*i+1]*gdata->hscale,
      ta_VHexchonly[ND*i+2]*gdata->hscale,
#endif
      ta_VHext[ND*i+0]*gdata->hscale,
      ta_VHext[ND*i+1]*gdata->hscale,
      ta_VHext[ND*i+2]*gdata->hscale
    );

    /* check for problems of snprintf */
    assert(nchars>0);

    /* update pointer to the end of the array */
    wbufe += nchars;

    /* check if we write beyond allocated area */
    assert(wbufe-wbufs < nmax);
  }

  VecRestoreArray(gdata->M,&ta_M);
  if (gdata->VHdem!=PETSC_NULL) {
    VecRestoreArray(gdata->VHdem,&ta_VHdem);
  }
  VecRestoreArray(gdata->VHexchani,&ta_VHexchani);
#ifdef EXCH
  VecRestoreArray(gdata->VHexchonly,&ta_VHexchonly);
#endif
  VecRestoreArray(gdata->VHext,&ta_VHext);

#if defined(ZLIB) && !defined(_WIN32) && !defined(__CYGWIN32__)
  SynchronizedFastFPrintf(fd,wbufe-wbufs,nmax,wbufs,PETSC_TRUE);
#else
  SynchronizedFastFPrintf(fd,wbufe-wbufs,nmax,wbufs,PETSC_FALSE);
#endif

  wbufe=wbufs;

#if defined(ZLIB) && !defined(_WIN32) && !defined(__CYGWIN32__)
  if (!rank) {
    gzclose(fd);
  }
#else
  ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
#endif

  ierr = PetscFree(wbufs);CHKERRQ(ierr);


  MagparFunctionInfoReturn(0);
}