writedatapng2.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: writedatapng.c 1246 2008-01-30 18:50:24Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "writepng.h"
#include "util/util.h"

static int  doinit=1;
static int  fieldon=0;
static int  pix[ND];   
static Mat  Afe2sq;    
static Vec  VIpol;     
int WritePNG2init(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);

  int         nmax;
  PetscTruth  flg;

  PetscReal slice_n[ND];
  nmax=ND;
  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-slice2_n",slice_n,(PetscInt*)&nmax,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    slice_n[0]=slice_n[1]=0;slice_n[2]=1;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -slice2_n not found, using default value: (%g,%g,%g)\n",
      slice_n[0],slice_n[1],slice_n[2]
    );
  }
  else if (nmax!=ND) {
    SETERRQ(PETSC_ERR_ARG_INCOMP,
      "Wrong number of values for option -slice2_n!\n");
  }

  PetscReal slice_p[ND];
  nmax=ND;
  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-slice2_p",slice_p,(PetscInt*)&nmax,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    slice_p[0]=slice_p[1]=slice_p[2]=PETSC_MAX;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -slice2_p not found, using default value: (%g,%g,%g)\n",
      slice_p[0],slice_p[1],slice_p[2]
    );

    PetscPrintf(PETSC_COMM_WORLD,"PNG output off\n");
    MagparFunctionLogReturn(0);
  }
  else if (nmax!=ND) {
    SETERRQ(PETSC_ERR_ARG_INCOMP,
      "Wrong number of values for option -slice2_p!\n");
  }

  int slice_g;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-slice2_g",(PetscInt*)&slice_g,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    slice_g=0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -slice2_g not found, using default value: %i\n",
      slice_g
    );
  }
  if (PetscAbsInt(slice_g)>gdata->n_prop) {
    SETERRQ(PETSC_ERR_ARG_INCOMP,
      "slice2_g > number property ids\n");
  }
  if (slice_g==0) {
    PetscPrintf(PETSC_COMM_WORLD,
      "slice2_g == 0: drawing magnetization of any volume\n"
    );
  }
  else if (slice_g<0) {
    PetscPrintf(PETSC_COMM_WORLD,
      "slice2_g < 0: drawing magnetization of any volume except pid %i\n",
      slice_g
    );
  }
  else {
    PetscPrintf(PETSC_COMM_WORLD,
      "slice2_g > 0: drawing only magnetization of volumes with pid %i\n",
      slice_g
    );
  }

  int res;
  ierr = PetscOptionsGetInt (PETSC_NULL,"-res",(PetscInt*)&res,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    res=50;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -res not found, using default value: %i\n",
      res
    );
  }

  PetscPrintf(PETSC_COMM_WORLD,"shift and tilt to avoid cutting 'between' elements\n");
#define DIST 1e-12
  slice_p[0]+=DIST;slice_n[0]+=DIST;
  slice_p[1]+=DIST;slice_n[1]+=DIST;
  slice_p[2]+=DIST;slice_n[2]+=DIST;
  my_dscal(ND,1.0+DIST,slice_n,1);
  my_dscal(ND,1.0+DIST,slice_p,1);

  PetscReal *vertxyz2;
  PetscReal axes2[ND*ND];
  ierr = PetscMalloc(ND*gdata->n_vert*sizeof(PetscReal),&vertxyz2);CHKERRQ(ierr);
  ierr = AxesRotation(slice_p,slice_n,axes2,0,gdata->n_vert,gdata->vertxyz,vertxyz2);CHKERRQ(ierr);

  /* find bounding box */
  slice_g--;
  PetscReal bbox[2*ND];
  pix[0]=res;
  ierr = bbox2(gdata,vertxyz2,1,&slice_g,bbox,pix);CHKERRQ(ierr);

  /* reset z-values to force slice cut through z=0 (vertxyz2 coordinates) */
  bbox[2]=0;bbox[ND+2]=D_EPS;

  ierr = CalAfe2sq(gdata,vertxyz2,bbox,pix,slice_g,&Afe2sq);CHKERRQ(ierr);
  ierr = PetscFree(vertxyz2);CHKERRQ(ierr);

  if (Afe2sq==PETSC_NULL) {
    MagparFunctionLogReturn(0);
  }
  fieldon=1;

  ierr = VecCreate(PETSC_COMM_WORLD,&VIpol);CHKERRQ(ierr);
  if (!rank) {
    ierr = VecSetSizes(VIpol,
      ND*pix[0]*pix[1],
      ND*pix[0]*pix[1]
    );CHKERRQ(ierr);
  }
  else {
    ierr = VecSetSizes(VIpol,0,ND*pix[0]*pix[1]);CHKERRQ(ierr);
  }
  ierr = VecSetFromOptions(VIpol);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}


int WritePNG2(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);

  if (doinit) {
    ierr = WritePNG2init(gdata);CHKERRQ(ierr);
    doinit=0;
  }
  if (!fieldon) {
    MagparFunctionInfoReturn(0);
  }

  ierr = MatMult(Afe2sq,gdata->M,VIpol);CHKERRQ(ierr);
  if (rank) {
    MagparFunctionInfoReturn(0);
  }

  png_bytepp image;
  ierr = PetscMalloc(pix[1]*sizeof(png_bytep),&image);CHKERRQ(ierr);
  ierr = PetscMalloc(3*pix[1]*pix[0]*sizeof(png_byte),&image[0]);CHKERRQ(ierr);

  for (int i=1;i<pix[1];i++) {
    image[i] = (png_bytep) (image[0]+3*pix[0]*i);
    assert(image[i]-image[i-1]==3*pix[0]);
  }

  PetscReal *a;
  VecGetArray(VIpol,&a);

  for (int k=0;k<ND;k++) {
    int mcomp;

    mcomp=k;

    for (int i=0;i<pix[0];i++) {
      for (int j=0;j<pix[1];j++) {
        png_byte red,green,blue;

        /* if the |M|<0.1 then we are outside any magnetic material and
           set the background color white */
        if (a[ND*(i*pix[1]+pix[1]-1-j)+0]*a[ND*(i*pix[1]+pix[1]-1-j)+0]+
            a[ND*(i*pix[1]+pix[1]-1-j)+1]*a[ND*(i*pix[1]+pix[1]-1-j)+1]+
            a[ND*(i*pix[1]+pix[1]-1-j)+2]*a[ND*(i*pix[1]+pix[1]-1-j)+2]
            < D_EPS
           ) {
          red=green=blue=255;
        }
        else {
          PetscReal t_val;
          switch(mcomp) {
            case 0:
              /* Mx */
              t_val=a[ND*(i*pix[1]+pix[1]-1-j)+0];
            break;
            case 1:
              /* My */
              t_val=a[ND*(i*pix[1]+pix[1]-1-j)+1];
            break;
            case 2:
              /* Mz */
              t_val=a[ND*(i*pix[1]+pix[1]-1-j)+2];
            break;
            default:
              SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"mcomp has illegal value\n");
            break;
          }

          /* t_val in [-1,+1] */
/*
          red=  255-(png_byte)((val-vmin)/(vmax-vmin)*255.0);
          green=255-(png_byte)PetscAbsInt((val-(vmax+vmin)/2.0)/(vmax-vmin)*2.0*255.0);
          blue=     (png_byte)((val-vmin)/(vmax-vmin)*255.0);
*/

          /* colormap:  -1 ... blue, 0 ... green, 1 ... red */
          /* have green in the middle */
          PetscReal    t_r,t_g,t_b;
          PetscReal    t_min=-1.0,t_max=+1.0,t_mid=0.5;

          t_b=(t_val < (t_max+t_min)*t_mid) ?
            t_val/(-(t_max+t_min)*t_mid+t_min)-
              (t_max+t_min)*t_mid/(-(t_max+t_min)*t_mid+t_min) :
            0.0;
          t_g=(t_val < (t_max+t_min)*t_mid) ?
            t_val/(+(t_max+t_min)*t_mid-t_min)-
              t_min/(+(t_max+t_min)*t_mid-t_min) :
            t_val/(-t_max+(t_max+t_min)*t_mid)-
              t_max/(-t_max+(t_max+t_min)*t_mid);
          t_r=(t_val > (t_max+t_min)*t_mid) ?
            t_val/(t_max-(t_max+t_min)*t_mid)-
              (t_max+t_min)*t_mid/(t_max-(t_max+t_min)*t_mid) :
            0.0;

          /* fix cases when |M|>1.0 */
          if (t_b>1.0) t_b=1.0; if (t_b<0.0) t_b=0.0;
          if (t_g>1.0) t_g=1.0; if (t_g<0.0) t_g=0.0;
          if (t_r>1.0) t_r=1.0; if (t_r<0.0) t_r=0.0;

          red=  (png_byte)(t_r*254.0);
          green=(png_byte)(t_g*254.0);
          blue= (png_byte)(t_b*254.0);
        }

/*      gcc: comparison is always 1 due to limited range of data type
        assert(red  >=0 && red   <= 255);
        assert(green>=0 && green <= 255);
        assert(blue >=0 && blue  <= 255);
*/

        image[j][3*i+0]=red;
        image[j][3*i+1]=green;
        image[j][3*i+2]=blue;
      }
    }

    char fmesh[256],desc[256];
    ierr = PetscSNPrintf(fmesh,255,"%s_%01i.%04i.%i.png",gdata->simname,2,gdata->inp,mcomp);CHKERRQ(ierr);
    ierr = PetscSNPrintf(desc,255,"t=%g ns, Etot=%g (J/m^3)",gdata->time*gdata->tscale*1e9,gdata->Etot*gdata->escale);CHKERRQ(ierr);
    WritePNGfile(fmesh,pix[0],pix[1],image,desc);
  }

  VecRestoreArray(VIpol,&a);

  ierr = PetscFree(image[0]);CHKERRQ(ierr);
  ierr = PetscFree(image);CHKERRQ(ierr);


  MagparFunctionInfoReturn(0);
}