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

#include "init.h"
#include "util/util.h"

PetscReal tetgrad(PetscReal *x1,PetscReal *x2,PetscReal *x3,PetscReal *x4,PetscReal D_etaj[][ND])
{
  PetscReal t_elevol;

  t_elevol=tetvol(x1,x2,x3,x4);

  /* calculate gradients */
  D_etaj[0][0]=-(x3[1]*x2[2]-x4[1]*x2[2]-x2[1]*x3[2]+x4[1]*x3[2]+x2[1]*x4[2]-x3[1]*x4[2]);
  D_etaj[0][1]=+(x3[0]*x2[2]-x4[0]*x2[2]-x2[0]*x3[2]+x4[0]*x3[2]+x2[0]*x4[2]-x3[0]*x4[2]);
  D_etaj[0][2]=-(x3[0]*x2[1]-x4[0]*x2[1]-x2[0]*x3[1]+x4[0]*x3[1]+x2[0]*x4[1]-x3[0]*x4[1]);

  D_etaj[1][0]=+(x4[1]*x3[2]-x1[1]*x3[2]-x3[1]*x4[2]+x1[1]*x4[2]+x3[1]*x1[2]-x4[1]*x1[2]);
  D_etaj[1][1]=-(x4[0]*x3[2]-x1[0]*x3[2]-x3[0]*x4[2]+x1[0]*x4[2]+x3[0]*x1[2]-x4[0]*x1[2]);
  D_etaj[1][2]=+(x4[0]*x3[1]-x1[0]*x3[1]-x3[0]*x4[1]+x1[0]*x4[1]+x3[0]*x1[1]-x4[0]*x1[1]);

  D_etaj[2][0]=-(x1[1]*x4[2]-x2[1]*x4[2]-x4[1]*x1[2]+x2[1]*x1[2]+x4[1]*x2[2]-x1[1]*x2[2]);
  D_etaj[2][1]=+(x1[0]*x4[2]-x2[0]*x4[2]-x4[0]*x1[2]+x2[0]*x1[2]+x4[0]*x2[2]-x1[0]*x2[2]);
  D_etaj[2][2]=-(x1[0]*x4[1]-x2[0]*x4[1]-x4[0]*x1[1]+x2[0]*x1[1]+x4[0]*x2[1]-x1[0]*x2[1]);

  D_etaj[3][0]=+(x2[1]*x1[2]-x3[1]*x1[2]-x1[1]*x2[2]+x3[1]*x2[2]+x1[1]*x3[2]-x2[1]*x3[2]);
  D_etaj[3][1]=-(x2[0]*x1[2]-x3[0]*x1[2]-x1[0]*x2[2]+x3[0]*x2[2]+x1[0]*x3[2]-x2[0]*x3[2]);
  D_etaj[3][2]=+(x2[0]*x1[1]-x3[0]*x1[1]-x1[0]*x2[1]+x3[0]*x2[1]+x1[0]*x3[1]-x2[0]*x3[1]);

  for (int i=0;i<NV;i++) {
    for (int j=0;j<ND;j++) {
      D_etaj[i][j] /= 6.0*t_elevol;
    }
  }

  /* sum of gradients has to be zero in each element, proof!? */
  for (int j=0;j<ND;j++) {
    if (PetscAbsReal(D_etaj[0][j]+D_etaj[1][j]+D_etaj[2][j]+D_etaj[3][j])>=D_EPS*1000) {
      PetscPrintf(PETSC_COMM_WORLD,
        "etaj: %g %g %g %g sum= %g\n",
        D_etaj[0][j],D_etaj[1][j],D_etaj[2][j],D_etaj[3][j],
        D_etaj[0][j]+D_etaj[1][j]+D_etaj[2][j]+D_etaj[3][j]
      );
    }
  }

  return(0);
}


PetscReal tetqual(PetscReal *x1,PetscReal *x2,PetscReal *x3,PetscReal *x4)
{
  PetscReal result;

  /* calculate ratio of the volume of the current element and the volume
     of a regular tetrahedron whose edge length is the average over all
     six edges of the current tetrahedron
     formulas: http://www.linmpi.mpg.de/~daly/CSDS/t4h/tetra.htm

     cf. VO Schrefl, Numerische Methoden, p. 44
  */

  result=tetvol(x1,x2,x3,x4)/(sqrt(2.0)/12.0*
    pow(
      (
        sqrt((x1[0]-x2[0])*(x1[0]-x2[0])+(x1[1]-x2[1])*(x1[1]-x2[1])+(x1[2]-x2[2])*(x1[2]-x2[2]))+
        sqrt((x1[0]-x3[0])*(x1[0]-x3[0])+(x1[1]-x3[1])*(x1[1]-x3[1])+(x1[2]-x3[2])*(x1[2]-x3[2]))+
        sqrt((x1[0]-x4[0])*(x1[0]-x4[0])+(x1[1]-x4[1])*(x1[1]-x4[1])+(x1[2]-x4[2])*(x1[2]-x4[2]))+
        sqrt((x2[0]-x3[0])*(x2[0]-x3[0])+(x2[1]-x3[1])*(x2[1]-x3[1])+(x2[2]-x3[2])*(x2[2]-x3[2]))+
        sqrt((x2[0]-x4[0])*(x2[0]-x4[0])+(x2[1]-x4[1])*(x2[1]-x4[1])+(x2[2]-x4[2])*(x2[2]-x4[2]))+
        sqrt((x3[0]-x4[0])*(x3[0]-x4[0])+(x3[1]-x4[1])*(x3[1]-x4[1])+(x3[2]-x4[2])*(x3[2]-x4[2]))
      )/6.0,
      3.0
    )
  );

  return(result);
}


PetscReal edgestat(int nele,int *ele,PetscReal *xyz,PetscReal *min,PetscReal *max,PetscReal *avg)
{
  PetscReal   *x1,*x2,*x3,*x4;
  PetscReal   edge[ND];
  PetscReal   elen;
  PetscReal   elenavg=0;
  PetscReal   elenmin=PETSC_MAX;
  PetscReal   elenmax=PETSC_MIN;

  MagparFunctionInfoBegin;

  for (int i=0; i<nele; i++) {

    /* Get Vertex coordinates */
    int *t_v;
    t_v=ele+i*NV;

    x1=xyz+t_v[0]*ND;
    x2=xyz+t_v[1]*ND;
    x3=xyz+t_v[2]*ND;
    x4=xyz+t_v[3]*ND;

    /* Calculate edge lengths */
    my_dcopy(ND,x1,1,edge,1);
    my_daxpy(ND,-1.0,x2,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);
    elenavg += elen;

    my_dcopy(ND,x1,1,edge,1);
    my_daxpy(ND,-1.0,x3,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);
    elenavg += elen;

    my_dcopy(ND,x1,1,edge,1);
    my_daxpy(ND,-1.0,x4,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);
    elenavg += elen;

    my_dcopy(ND,x2,1,edge,1);
    my_daxpy(ND,-1.0,x3,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);
    elenavg += elen;

    my_dcopy(ND,x2,1,edge,1);
    my_daxpy(ND,-1.0,x4,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);
    elenavg += elen;

    my_dcopy(ND,x3,1,edge,1);
    my_daxpy(ND,-1.0,x4,1,edge,1);
    elen=my_dnrm2(ND,edge,1);
    elenmax=PetscMax(elenmax,elen);
    elenmin=PetscMin(elenmin,elen);
    elenavg += elen;
  }

  ierr = PetscGlobalMin(&elenmin,min,PETSC_COMM_WORLD);CHKERRQ(ierr);
  ierr = PetscGlobalMax(&elenmax,max,PETSC_COMM_WORLD);CHKERRQ(ierr);
  ierr = PetscGlobalSum(&elenavg,avg,PETSC_COMM_WORLD);CHKERRQ(ierr);

  MagparFunctionInfoReturn(0);
}


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

  PetscReal *vprop;
  ierr = PetscMalloc(gdata->n_prop*sizeof(PetscReal),&vprop);CHKERRQ(ierr);
  for (int i=0;i<gdata->n_prop;i++) {
    vprop[i]=0.0;
  }

  ierr = VecCreate(PETSC_COMM_WORLD,&gdata->elevol);CHKERRQ(ierr);
  ierr = VecSetSizes(gdata->elevol,gdata->ln_ele,gdata->n_ele);CHKERRQ(ierr);
  ierr = VecSetFromOptions(gdata->elevol);CHKERRQ(ierr);

  Vec elequal;
  ierr = VecDuplicate(gdata->elevol,&elequal);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD,&gdata->vertvol);CHKERRQ(ierr);
  ierr = VecSetSizes(gdata->vertvol,gdata->ln_vert,gdata->n_vert);CHKERRQ(ierr);
  ierr = VecSetFromOptions(gdata->vertvol);CHKERRQ(ierr);

  ierr = VecDuplicate(gdata->M,&gdata->VMs3);CHKERRQ(ierr);

  PetscReal Javg=0.0;
  PetscReal Vmag=0.0;
  for (int i=0; i<gdata->ln_ele; i++) {
    /* get vertex coordinates */
    int *t_v;
    t_v=gdata->elevert+i*NV;

    PetscReal *x1,*x2,*x3,*x4;
    x1=gdata->vertxyz+t_v[0]*ND;
    x2=gdata->vertxyz+t_v[1]*ND;
    x3=gdata->vertxyz+t_v[2]*ND;
    x4=gdata->vertxyz+t_v[3]*ND;

    /* get element volume */
    PetscReal T;
    T=tetvol(x1,x2,x3,x4);

    /* fix elements with negative volume by switching node numbers */
    if (T<0.0) {
/*
      PetscPrintf(PETSC_COMM_WORLD,
        "element (local: %i, global: %i) (%i,%i,%i,%i) has negative volume - fixed!\n",
        i,gdata->elel2g[i],t_v[0],t_v[1],t_v[2],t_v[3]
      );
*/
      int j;
      j=t_v[0];
      t_v[0]=t_v[1];
      t_v[1]=j;
      T=-T;

      x1=gdata->vertxyz+t_v[0]*ND;
      x2=gdata->vertxyz+t_v[1]*ND;
    }

    if (gdata->propdat[NP*gdata->eleprop[i]+4]>0.0) {
      Vmag += T;
    }

    ierr = VecSetValue(
      gdata->elevol,
      gdata->elel2g[i],
      T,
      INSERT_VALUES
    );CHKERRQ(ierr);

    for (int j=0;j<NV;j++) {
      ierr = VecSetValue(
        gdata->vertvol,
        t_v[j],
        T/NV,
        ADD_VALUES
      );CHKERRQ(ierr);

      PetscReal matele;
      matele=T/NV*gdata->propdat[NP*gdata->eleprop[i]+4];
      if (PetscAbsReal(matele)>D_EPS) {
        for (int k=0;k<ND;k++) {
          ierr = VecSetValue(
            gdata->VMs3,
            gdata->elevert[NV*i+j]*ND+k,
            matele,
            ADD_VALUES
          );CHKERRQ(ierr);
        }
      }
    }
    Javg += T*gdata->propdat[NP*gdata->eleprop[i]+4]*gdata->hscale;
    vprop[gdata->eleprop[i]] += T;

    PetscReal t_qual;
    t_qual=tetqual(x1,x2,x3,x4);
    ierr = VecSetValue(
      elequal,
      gdata->elel2g[i],
      t_qual,
      INSERT_VALUES
    );CHKERRQ(ierr);
  }

  ierr = VecAssemblyBegin(gdata->elevol);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(gdata->elevol);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(gdata->vertvol);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(gdata->vertvol);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(gdata->VMs3);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(gdata->VMs3);CHKERRQ(ierr);

  for (int i=0;i<gdata->n_prop;i++) {
    PetscReal temp;
    ierr = PetscGlobalSum(&vprop[i],&temp,PETSC_COMM_WORLD);CHKERRQ(ierr);
    vprop[i]=temp;
  }

  /* save *.felog --------------------------------------- */

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

  /* mesh stats --------------------------------------- */

  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "simname: %s\n"
    "n_vert:     %7d\n"
    "n_ele:      %7d\n"
    "n_vert_bnd: %7d\n"
    "n_bnd_fac:  %7d\n"
    "n_prop:     %7d\n",
    gdata->simname,
    gdata->n_vert,
    gdata->n_ele,
    gdata->n_vert_bnd,
    gdata->n_bnd_fac,
    gdata->n_prop
  );

  /* bounding box --------------------------------------- */
  PetscReal   bbox[2*ND];
  int         pid[1]={-1};

  ierr = calcbbox(gdata->ln_ele,gdata->elevert,PETSC_NULL,gdata->vertxyz,1,pid,bbox);CHKERRQ(ierr);

  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "\nbounding box: (%g,%g,%g) (%g,%g,%g)\n",
    bbox[0],bbox[1],bbox[2],
    bbox[ND+0],bbox[ND+1],bbox[ND+2]
  );CHKERRQ(ierr);

  /* elevol --------------------------------------- */

  PetscReal Vtot,Vmin,Vmax;
  int t_idxmin,t_idxmax;
  ierr = VecMin(gdata->elevol,(PetscInt*)&t_idxmin,&Vmin);CHKERRQ(ierr);
  ierr = VecMax(gdata->elevol,(PetscInt*)&t_idxmax,&Vmax);CHKERRQ(ierr);
  ierr = VecSum(gdata->elevol,&Vtot);CHKERRQ(ierr);
  gdata->totvol=Vtot;

  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "\nelevol:\n"
    "  id_min: %i\n"
    "  id_max: %i\n"
    "  Vmax:   %g\n"
    "  Vmin:   %g\n"
    "  ratio:  %g\n"
    "  Vavg:   %g\n"
    "  Vtot:   %g\n",
    t_idxmin,t_idxmax,Vmax,Vmin,Vmax/Vmin,Vtot/gdata->n_ele,Vtot
  );

  /* vertvol --------------------------------------- */

  ierr = VecMin(gdata->vertvol,(PetscInt*)&t_idxmin,&Vmin);CHKERRQ(ierr);
  ierr = VecMax(gdata->vertvol,(PetscInt*)&t_idxmax,&Vmax);CHKERRQ(ierr);
  ierr = VecSum(gdata->vertvol,&Vtot);CHKERRQ(ierr);

  /* check that elevoltot == vertvoltot */
  assert(PetscAbsReal(gdata->totvol-Vtot) < Vtot*1e-7);

  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "\nvertvol:\n"
    "  id_min: %i\n"
    "  id_max: %i\n"
    "  Vmax:   %g\n"
    "  Vmin:   %g\n"
    "  ratio:  %g\n"
    "  Vavg:   %g\n"
    "  Vtot:   %g\n",
    t_idxmin,t_idxmax,Vmax,Vmin,Vmax/Vmin,Vtot/gdata->n_ele,Vtot
  );

  /* edgestat --------------------------------------- */

  PetscReal elenmax,elenmin,elenavg;
  edgestat(gdata->ln_ele,gdata->elevert,gdata->vertxyz,&elenmin,&elenmax,&elenavg);
  elenavg /= 6*gdata->n_ele;
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "\n"
    "edge_len_min: %g\n"
    "edge_len_max: %g\n"
    "edge_len_avg: %g\n",
    elenmin,elenmax,elenavg
  );
  gdata->elenmax=elenmax;

  /* elequal --------------------------------------- */

  /* TODO: abusing Vmin, Vmax, Vtot */
  ierr = VecMin(elequal,(PetscInt*)&t_idxmin,&Vmin);CHKERRQ(ierr);
  ierr = VecMax(elequal,(PetscInt*)&t_idxmax,&Vmax);CHKERRQ(ierr);
  ierr = VecSum(elequal,&Vtot);CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "\nelequal:\n"
    "  min:    %g\n"
    "  max:    %g\n"
    "  avg:    %g\n",
    Vmin,Vmax,Vtot/gdata->n_ele
  );

  ierr = VecDestroy(elequal);CHKERRQ(ierr);

  /* vol by prop --------------------------------------- */

  Vtot=0.0;
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "\nvolume by property id\n\n"
    "  pid | vol\n"
    "-----------\n"
  );
  for (int i=0;i<gdata->n_prop;i++) {
    ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
      "%5d   %g\n",
      i+1,vprop[i]
    );
    Vtot += vprop[i];
  }
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "-----------\n"
    "  sum = %g\n",
    Vtot
  );
#ifdef PETSC_USE_SINGLE
  assert(PetscAbsReal(gdata->totvol-Vtot) < Vtot*1e-3);
#else
  assert(PetscAbsReal(gdata->totvol-Vtot) < Vtot*1e-7);
#endif

  /* magnetic Vtot ------------------------------ */

  ierr = PetscGlobalSum(&Vmag,&Vtot,PETSC_COMM_WORLD);CHKERRQ(ierr);
  gdata->totvol=Vtot;

  /* Javg --------------------------------------- */

  PetscReal sum;
  ierr = PetscGlobalSum(&Javg,&sum,PETSC_COMM_WORLD);CHKERRQ(ierr);
  Javg=sum/Vtot;
  ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
    "\naverage magnetization: %g T\n",
    Javg
  );

  /* done --------------------------------------- */

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

  MagparFunctionLogReturn(0);
}