writedatadat.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: writedatadat.c 3051 2010-04-14 17:58:54Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "magpario.h"
#include "util/util.h"
#ifdef ADDONS
#include "addons/addons.h"
#endif

#define NVDATA 18
#define NEDATA  1
#define NEDATA2 0

static int  doinit=1;
static Mat  AIpolline; 
static Vec  VIpolline; 
int WriteDatInit(GridData *gdata)
{
  int         i,j,k,l;

  int         nmax;
  PetscTruth  flg;

  int         res;
  PetscReal   *vertxyz2;
  PetscReal   line_p[ND], line_v[ND];

  PetscReal   axes2[ND*ND];

  int         cnt,*sliceele;
  PetscReal   sgnx,sgny;
  PetscReal   x2min[ND]={0,0,PETSC_MAX};
  PetscReal   x2max[ND]={0,0,PETSC_MIN};
  PetscReal   xmin[ND],xmax[ND];
  PetscReal   line_length;
  PetscReal   *d,p[ND];
  PetscReal   rhs[ND];
  int         t_ele;

  PetscReal   matele;


  MagparFunctionLogBegin;

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

  ierr = PetscMalloc(ND*gdata->n_vert*sizeof(PetscReal),&vertxyz2);CHKERRQ(ierr);

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

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

  PetscPrintf(PETSC_COMM_WORLD,
    "sampling line: (x,y,z)=(%g,%g,%g)+n*(%g,%g,%g)\n",
    line_p[0],line_p[1],line_p[2],
    line_v[0],line_v[1],line_v[2]
  );

  /* calculate unit vector */
  PetscReal nrm;
  nrm=1.0/my_dnrm2(ND,line_v,1);
  my_dscal(ND,nrm,line_v,1);
  /* this gives incorrect line_v (!!!): my_dscal(ND,1.0/my_dnrm2(ND,line_v,1),line_v,1); */

  PetscPrintf(PETSC_COMM_WORLD,
    "normalized sampling line vector: (%g,%g,%g)\n",
    line_v[0],line_v[1],line_v[2]
  );

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

  /* line_v will be z-axes in axes2 coordinates */
  AxesRotation(line_p,line_v,axes2,0,gdata->n_vert,gdata->vertxyz,vertxyz2);

  ierr = PetscMalloc(gdata->ln_ele*sizeof(int),&sliceele);CHKERRQ(ierr);

  cnt=0;
  for (i=0;i<gdata->ln_ele;i++) {
    /* get x- and y-coordinate of first vertex */
    sgnx=vertxyz2[ND*gdata->elevert[NV*i+0]+0];
    sgny=vertxyz2[ND*gdata->elevert[NV*i+0]+1];

    for (j=1;j<NV;j++) {
      if (sgnx*vertxyz2[ND*gdata->elevert[NV*i+j]+0] <= 0.0) break;
    }
    if (j<NV) {
      for (j=1;j<NV;j++) {
        /* multiply with other y-coordinate of other vertices
           if (product < 0) then they have opposite signs and
           the element must be cut by the sampling line
        */
        if (sgny*vertxyz2[ND*gdata->elevert[NV*i+j]+1] <= 0.0) break;
      }
      /* if an element cut by the sampling line has been found */
      if (j<NV) {
        /* add its local index to the array */
        sliceele[cnt++]=i;
      }
    }
  }
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,
    "<%i>Found %i elements cut by sampling line\n",
    rank,cnt
  );
  PetscSynchronizedFlush(PETSC_COMM_WORLD);

  int sum;
  ierr=MPI_Allreduce((int*)&cnt,(int*)&sum,1,MPI_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRQ(ierr);
  if (sum==0) {
    AIpolline=PETSC_NULL;

    PetscPrintf(PETSC_COMM_WORLD,
      "Sampling line outside the FE mesh.\nNo *.d files will be stored!\n"
    );

    ierr = PetscFree(sliceele);CHKERRQ(ierr);
    ierr = PetscFree(vertxyz2);CHKERRQ(ierr);

    MagparFunctionLogReturn(0);
  }

  /* get number of pixels */
  ierr=MPI_Allreduce((int*)&cnt,(int*)&res,1,MPI_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRQ(ierr);
  res*=5;

  PetscPrintf(PETSC_COMM_WORLD,"sampling points: %i\n",res);

  /* find bbox of probing line in axes2 coordinates */
  for (i=0;i<cnt;i++) {
    for (j=0;j<NV;j++) {
      x2min[2]=PetscMin(x2min[2],vertxyz2[ND*gdata->elevert[NV*sliceele[i]+j]+2]);
      x2max[2]=PetscMax(x2max[2],vertxyz2[ND*gdata->elevert[NV*sliceele[i]+j]+2]);
    }
  }

  ierr = PetscGlobalMin(&x2min[2],&xmin[2],PETSC_COMM_WORLD);CHKERRQ(ierr);
  x2min[2]=xmin[2];
  ierr = PetscGlobalMax(&x2max[2],&xmax[2],PETSC_COMM_WORLD);CHKERRQ(ierr);
  x2max[2]=xmax[2];

  PetscPrintf(PETSC_COMM_WORLD,
    "sampling line ends in local coordinates (axes2):\n"
    "(%g,%g,%g) (%g,%g,%g)\n",
    x2min[0],x2min[1],x2min[2],
    x2max[0],x2max[1],x2max[2]
  );

  /* add a little space on both ends */
  x2min[2] -= 1.0*(x2max[2]-x2min[2])/(res);
  x2max[2] += 2.0*(x2max[2]-x2min[2])/(res);
  PetscPrintf(PETSC_COMM_WORLD,
    "enlarged slice corners in local slice plane coordinates (axes2):\n"
    "(%g,%g,%g) (%g,%g,%g)\n",
    x2min[0],x2min[1],x2min[2],
    x2max[0],x2max[1],x2max[2]
  );
  /* get length of slice edges */
  line_length=x2max[2]-x2min[2];

  /* find bbox of cut probing line in axes coordinates (rotate back) */
  AxesRotation(line_p,PETSC_NULL,axes2,1,1,x2min,xmin);
  AxesRotation(line_p,PETSC_NULL,axes2,1,1,x2max,xmax);

  PetscPrintf(PETSC_COMM_WORLD,
    "sampling line ends in world coordinates (axes):\n"
    "(%g,%g,%g) (%g,%g,%g)\n",
    xmin[0],xmin[1],xmin[2],
    xmax[0],xmax[1],xmax[2]
  );

  if (!rank) {
    PetscPrintf(PETSC_COMM_WORLD,"length of sampling line: %g\n",line_length);

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

    ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
      "#.%4s %14s %14s %14s %14s\n"
      "#.%4s %14s %14s %14s %14s\n",
      "..id","dist","x","y","z",
      "...-","-","-","-","-"
    );CHKERRQ(ierr);

    for (i=0;i<res;i++) {
      /* calculate xyz coordinates of pixel */
      my_dcopy(ND,xmin,1,p,1);
      my_daxpy(ND,i*line_length/res,axes2+ND*2,1,p,1);

      ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
        "  %4d %14e %14e %14e %14e\n",
        i,i*line_length/res,
        p[0],p[1],p[2]
      );
    }

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


  /* why "MatAssemblyEnd_SeqAIJMaximum nonzeros in any row is 5" ? */
  if (!rank) {
    ierr = MatCreateMPIAIJ(
      PETSC_COMM_WORLD,
      ND*res,ND*gdata->ln_vert,
      ND*res,ND*gdata->n_vert,
      NV+1,PETSC_NULL,NV+1,PETSC_NULL,
      &AIpolline
    );CHKERRQ(ierr);
  }
  else {
    ierr = MatCreateMPIAIJ(
      PETSC_COMM_WORLD,
      0,ND*gdata->ln_vert,
      ND*res,ND*gdata->n_vert,
      NV,PETSC_NULL,NV,PETSC_NULL,
      &AIpolline
    );CHKERRQ(ierr);
  }
  ierr = MatSetFromOptions(AIpolline);CHKERRQ(ierr);

  if (cnt>0) {
    t_ele=sliceele[0];
    d=gdata->vertxyz+ND*gdata->elevert[NV*t_ele+3];

    for (i=0;i<res;i++) {
      /* calculate xyz coordinates of pixel */

      my_dcopy(ND,xmin,1,p,1);
      my_daxpy(ND,i*line_length/res,axes2+ND*2,1,p,1);

      /* find the tetrahedron, in which it is located */
      for (k=-1;k<cnt;k++) {
        /* let's try with the element t_ele of our predecessor
           first (k==-1) (we are close for sure!)
        */
        if (k>=0) {
          t_ele=sliceele[k];
          d=gdata->vertxyz+ND*gdata->elevert[NV*t_ele+3];
        }

        /* calculate barycentric coordinates
           info from: 3Ddata.pdf
           Course in4008 - Data Visualization
           http://visualization.tudelft.nl/frits.html
        */

        /* calculate rhs=p-d */
        my_dcopy(ND,p,1,rhs,1);
        my_daxpy(ND,-1.0,d,1,rhs,1);

        /* skip element, if we are too far away */
        if (my_dnrm2(ND,rhs,1)>gdata->elenmax) continue;

        if (
          barycent(
            p,
            gdata->vertxyz+ND*gdata->elevert[NV*t_ele+0],
            gdata->vertxyz+ND*gdata->elevert[NV*t_ele+1],
            gdata->vertxyz+ND*gdata->elevert[NV*t_ele+2],
            gdata->vertxyz+ND*gdata->elevert[NV*t_ele+3],
            rhs
          )
        ) {
          for (l=0;l<ND;l++) {
            matele=(1-rhs[0]-rhs[1]-rhs[2]);
            ierr = MatSetValue(AIpolline,
              ND*i+l,
              ND*gdata->elevert[NV*t_ele+3]+l,
            matele,ADD_VALUES);CHKERRQ(ierr);

            matele=rhs[0];
            ierr = MatSetValue(AIpolline,
              ND*i+l,
              ND*gdata->elevert[NV*t_ele+0]+l,
            matele,ADD_VALUES);CHKERRQ(ierr);

            matele=rhs[1];
            ierr = MatSetValue(AIpolline,
              ND*i+l,
              ND*gdata->elevert[NV*t_ele+1]+l,
            matele,ADD_VALUES);CHKERRQ(ierr);

            matele=rhs[2];
            ierr = MatSetValue(AIpolline,
              ND*i+l,
              ND*gdata->elevert[NV*t_ele+2]+l,
            matele,ADD_VALUES);CHKERRQ(ierr);
          }
          break;
        }
      }
    }
  }

  /* FIXME: large number of mallocs during MatSetValue!? */
  PetscInfo(0,"AIpolline matrix assembly complete\n");
  ierr = MatAssemblyBegin(AIpolline,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(AIpolline,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = PrintMatInfoAll(AIpolline);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD,&VIpolline);CHKERRQ(ierr);
  if (!rank) {
    ierr = VecSetSizes(VIpolline,ND*res,ND*res);CHKERRQ(ierr);
  }
  else {
    ierr = VecSetSizes(VIpolline,0,ND*res);CHKERRQ(ierr);
  }
  ierr = VecSetFromOptions(VIpolline);CHKERRQ(ierr);

  ierr = PetscFree(sliceele);CHKERRQ(ierr);
  ierr = PetscFree(vertxyz2);CHKERRQ(ierr);


  MagparFunctionLogReturn(0);
}



int WriteDat(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 = WriteDatInit(gdata);CHKERRQ(ierr);
    doinit=0;
  }
  if (AIpolline==PETSC_NULL) {
    MagparFunctionInfoReturn(0);
  }

  char fmesh[256];
  FILE *fd;
  ierr = PetscSNPrintf(fmesh,255,"%s.%04i.d",gdata->simname,gdata->inp);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);

#ifdef ADDONS
/* Htot (useful for measuring "real" field in airbox, where Hexchani==0)
*/
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "#.%4s %14s %14s %14s\n"
    "#.%4s %14s %14s %14s\n",
    "..id","Htot_x","Htot_y","Htot_z",
    "....","(T)","(T)","(T)"
  );CHKERRQ(ierr);

  ierr = MatMult(AIpolline,gdata->VHtot,VIpolline);CHKERRQ(ierr);
  ierr = VecScale(VIpolline,gdata->hscale);CHKERRQ(ierr);

#else
/* magnetization */
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "#.%4s %14s %14s %14s\n"
    "#.%4s %14s %14s %14s %s\n",
    "..id","Mx","My","Mz",
    "...-","(1)","(1)","(1)",
    "(normalized magnetization, not proportional to Js)"
  );CHKERRQ(ierr);

  ierr = MatMult(AIpolline,gdata->M,VIpolline);CHKERRQ(ierr);
#endif

  if (!rank) {
    int i,j;
    VecGetSize(VIpolline,(PetscInt*)&i);
    VecGetLocalSize(VIpolline,(PetscInt*)&j);
    /* the whole vector has to be on the first processor! */
    assert(i==j);
    PetscReal *ta_dat;
    VecGetArray(VIpolline,&ta_dat);

    for (i=0; i<j/ND; i++) {
      ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
        "  %4d %14e %14e %14e\n",
        i,
        ta_dat[ND*i+0],
        ta_dat[ND*i+1],
        ta_dat[ND*i+2]
      );
    }
    VecRestoreArray(VIpolline,&ta_dat);
  }

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

  MagparFunctionInfoReturn(0);
}