writefemavs.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 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: writefemavs.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

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

#ifdef ADDONS
#define NVDATA  10
#else
#define NVDATA  4
#endif
#define NEDATA  11
#define NEDATA2 0

#define D_FLT 50
#define D_EID 8

int WriteFEMAVS(GridData *gdata)
{
  MagparFunctionLogBegin;

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

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

#if defined(ZLIB) && !defined(_WIN32) && !defined(__CYGWIN32__)
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"%i %i ??? ??? ???\n",
    gdata->n_vert,gdata->n_ele);CHKERRQ(ierr);
#else
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"%i %i %i %i %i\n",
    gdata->n_vert,gdata->n_ele,NVDATA,NEDATA,0);CHKERRQ(ierr);
#endif

  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);
  }

  char *wbufs,*wbufe;
  int nmax;
  /* add 1 so we allocate at least once char even if gdata->ln_ele==0 */
  nmax=gdata->ln_ele*(D_EID+6*D_FLT)+1;
  ierr = PetscMalloc(nmax*sizeof(char),&wbufs);CHKERRQ(ierr);
  wbufe=wbufs;

  for (int i=0; i<gdata->ln_ele; i++) {
    int nchars;
    /* write into buffer */
    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;

#if defined(ZLIB) && !defined(_WIN32) && !defined(__CYGWIN32__)
  ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);


  ierr = PetscSNPrintf(fmesh,255,"%s.%04i.%s",gdata->simname,gdata->inp,"fedat");CHKERRQ(ierr);
  ierr = PetscFOpen(PETSC_COMM_WORLD,fmesh,"w",&fd);CHKERRQ(ierr);
  if (!rank && !fd) {
    SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",fmesh);
  }
  PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",fmesh);
#endif

  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"%i",NVDATA);CHKERRQ(ierr);
  for (int i=0;i<NVDATA;i++) {
    ierr = PetscFPrintf(PETSC_COMM_WORLD,fd," 1");CHKERRQ(ierr);
  }
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"\n");CHKERRQ(ierr);

  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"id, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"vert_vol, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"vert_pid, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"proc, none\n");CHKERRQ(ierr);
#ifdef ADDONS
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"H1x, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"H1y, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"H1z, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"H2x, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"H2y, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"H2z, none\n");CHKERRQ(ierr);
#endif

  PetscReal *ta_vvol;
  VecGetArray(gdata->vertvol,&ta_vvol);

#ifdef ADDONS
  PetscReal *ta_vh1;
  PetscReal *ta_vh2;
  if (gdata->VH1!=PETSC_NULL) {
    VecGetArray(gdata->VH1,&ta_vh1);
  }
  if (gdata->VH2!=PETSC_NULL) {
    VecGetArray(gdata->VH2,&ta_vh2);
  }
#endif

  for (int i=0; i<gdata->ln_vert; i++) {
    int nchars;
    /* write into buffer */
    nchars=
    snprintf(
      wbufe,
      nmax-(wbufe-wbufs),
      "%i %i %g %i %i",
      gdata->vertl2g[i]+1,
      gdata->vertl2g[i]+1,
      ta_vvol[i],
      gdata->vertprop[i],
      rank
    );
#ifdef ADDONS
    if (gdata->VH1!=PETSC_NULL) {
      nchars+=
      snprintf(
        wbufe+nchars,
        nmax-(wbufe+nchars-wbufs),
        " %g %g %g",
        ta_vh1[ND*i+0]*gdata->hscale/MU0,
        ta_vh1[ND*i+1]*gdata->hscale/MU0,
        ta_vh1[ND*i+2]*gdata->hscale/MU0
      );
    }
    else {
      nchars+=
      snprintf(
        wbufe+nchars,
        nmax-(wbufe+nchars-wbufs),
        " 0 0 0"
      );
    }
    if (gdata->VH2!=PETSC_NULL) {
      nchars+=
      snprintf(
        wbufe+nchars,
        nmax-(wbufe+nchars-wbufs),
        " %g %g %g",
        ta_vh2[ND*i+0]*gdata->hscale/MU0,
        ta_vh2[ND*i+1]*gdata->hscale/MU0,
        ta_vh2[ND*i+2]*gdata->hscale/MU0
      );
    }
    else {
      nchars+=
      snprintf(
        wbufe+nchars,
        nmax-(wbufe+nchars-wbufs),
        " 0 0 0"
      );
    }
    assert(NVDATA==10);
#else
    assert(NVDATA==4);
#endif
    nchars+=
    snprintf(
      wbufe+nchars,
      nmax-(wbufe+nchars-wbufs),
      "\n"
    );

    /* 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->vertvol,&ta_vvol);

#ifdef ADDONS
  if (gdata->VH1!=PETSC_NULL) {
    VecRestoreArray(gdata->VH1,&ta_vh1);
  }
  if (gdata->VH2!=PETSC_NULL) {
    VecRestoreArray(gdata->VH2,&ta_vh2);
  }
#endif

  SynchronizedFastFPrintf(fd,wbufe-wbufs,nmax,wbufs,PETSC_FALSE);
  wbufe=wbufs;


  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"%i",NEDATA);CHKERRQ(ierr);
  for (int i=0;i<NEDATA;i++) {
    ierr = PetscFPrintf(PETSC_COMM_WORLD,fd," 1");CHKERRQ(ierr);
  }
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"\n");CHKERRQ(ierr);

  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"id, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"proc, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"vol, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"prop, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"a_x, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"a_y, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"a_z, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"K1, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"K2, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"Js, none\n");CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,"A, none\n");CHKERRQ(ierr);

  PetscReal *ta_vol;
  VecGetArray(gdata->elevol,&ta_vol);

  for (int i=0; i<gdata->ln_ele; i++) {
    int nchars;
    /* write into buffer */
    nchars=
    snprintf(
      wbufe,
      nmax-(wbufe-wbufs),
      "%i %i %i %g %i %g %g %g %g %g %g %g\n",
      gdata->elel2g[i]+1,
      gdata->elel2g[i]+1,
      rank,
      PetscRealPart(ta_vol[i]),
      gdata->eleprop[i]+1,
      sin(gdata->propdat[NP*gdata->eleprop[i]+0])*cos(gdata->propdat[NP*gdata->eleprop[i]+1]),
      sin(gdata->propdat[NP*gdata->eleprop[i]+0])*sin(gdata->propdat[NP*gdata->eleprop[i]+1]),
      cos(gdata->propdat[NP*gdata->eleprop[i]+0]),
      gdata->propdat[NP*gdata->eleprop[i]+2],
      gdata->propdat[NP*gdata->eleprop[i]+3],
      gdata->propdat[NP*gdata->eleprop[i]+4],
      gdata->propdat[NP*gdata->eleprop[i]+5]
    );
    /* 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->elevol,&ta_vol);

  SynchronizedFastFPrintf(fd,wbufe-wbufs,nmax,wbufs,PETSC_FALSE);
  wbufe=wbufs;

  ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);

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

  MagparFunctionLogReturn(0);
}