area.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: area.c 2992 2010-03-10 21:57:41Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "util.h"

int AreaCal(Mat AIpol, PetscReal *p_area)
{
  PetscReal area;

  int ngrows,ngcols;
  int nlrows,nlcols;

  Vec       test1,test2;
  PetscReal vmin,vmax;
  int       imin,imax;

  MagparFunctionLogBegin;

  ierr = MatGetSize(AIpol,&ngrows,&ngcols);CHKERRQ(ierr);
  ierr = MatGetLocalSize(AIpol,&nlrows,&nlcols);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD,&test1);CHKERRQ(ierr);
  ierr = VecCreate(PETSC_COMM_WORLD,&test2);CHKERRQ(ierr);
  ierr = VecSetSizes(test1,nlcols,ngcols);CHKERRQ(ierr);
  ierr = VecSetSizes(test2,nlrows,ngrows);CHKERRQ(ierr);
  ierr = VecSetFromOptions(test1);CHKERRQ(ierr);
  ierr = VecSetFromOptions(test2);CHKERRQ(ierr);

  ierr = VecSet(test1,1.0);CHKERRQ(ierr);

  ierr = MatMult(AIpol,test1,test2);CHKERRQ(ierr);

  /* check result of interpolation */
  ierr = VecMin(test2,&imin,&vmin);CHKERRQ(ierr);
  ierr = VecMax(test2,&imax,&vmax);CHKERRQ(ierr);

  if (vmin<0.0 || vmax>1.0+D_EPS) {
    PetscPrintf(PETSC_COMM_WORLD,
      "need to rescale matrix: %i: %e  %i: %e\n",
      imin,vmin,imax,vmax
    );

    /* TODO: rescale AIpol (not sure why this is necessary)
       seems like some pixels belong to two elements (just on interface!?)
    */
    PetscReal *atest2;
    ierr = VecGetArray(test2,&atest2);CHKERRQ(ierr);
    for (int i=0;i<nlrows;i++) {
      if (atest2[i]>1.0) atest2[i]=1.0/atest2[i];
    }
    ierr = VecRestoreArray(test2,&atest2);CHKERRQ(ierr);

    ierr = MatDiagonalScale(AIpol,test2,PETSC_NULL);CHKERRQ(ierr);

    ierr = MatMult(AIpol,test1,test2);CHKERRQ(ierr);
  }

  /* double check result of interpolation */
  ierr = VecMin(test2,&imin,&vmin);CHKERRQ(ierr);
  ierr = VecMax(test2,&imax,&vmax);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,
    "min/max: %i: %e  %i: %e\n",
    imin,vmin,imax,vmax
  );

  assert(vmin>=0.0);
  assert(vmax<=1.0+D_EPS);

  /* area calculation */
  ierr = VecSum(test2,&area);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"area = %g/%i=%g%%\n",area,ngrows,area/ngrows*100);

  *p_area=area;

  ierr = VecDestroy(test1);
  ierr = VecDestroy(test2);

  MagparFunctionLogReturn(0);
}


/* FIXME
 * save elevertall in griddata.h during initialization
 * save/use connectivity matrices:
 *   elements sharing faces
 *   elements sharing nodes
 *   nodes connected by edges
 */

int PidSharedArea(GridData *gdata,int pid1,int pid2,PetscReal *sharedarea_area,PetscReal *sharedarea_point,PetscReal *sharedarea_norm)
{
  MagparFunctionLogBegin;

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

  if (pid1==pid2)
    SETERRQ2(PETSC_ERR_ARG_INCOMP,"Pids have to be different to calculate shared area (%i==%i)!\n",pid1+1,pid2+1);

  if (pid1>=gdata->n_prop)
    SETERRQ2(PETSC_ERR_ARG_INCOMP,"pid1 too large (%i>%i)!\n",pid1+1,gdata->n_prop);
  if (pid2>=gdata->n_prop)
    SETERRQ2(PETSC_ERR_ARG_INCOMP,"pid2 too large (%i>%i)!\n",pid1+1,gdata->n_prop);
  if (pid1<=-1 && pid2<=-1)
    SETERRQ2(PETSC_ERR_ARG_INCOMP,"pids too small (%i,%i)<1!\n",pid1+1,pid2+1);

  PetscPrintf(PETSC_COMM_WORLD,"Search shared area for pids %i, %i\n",pid1+1,pid2+1);

  /* get ownership range */
  int low,high;
  ierr = VecGetOwnershipRange(gdata->elevol,&low,&high);CHKERRQ(ierr);

  int *elevertall;

#ifdef HAVE_ELEVERTALL
  elevertall=gdata->elevertall;
#else
  if (size>1) {
    ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&elevertall);CHKERRQ(ierr);
    ierr = PetscMemcpy(elevertall+NV*low,gdata->elevert,NV*gdata->ln_ele*sizeof(int));CHKERRQ(ierr);

    int *p;
    p=elevertall;

    int recvcount[size];
    recvcount[rank]=NV*gdata->ln_ele;
    for (int i=0;i<size;i++) {
      ierr = MPI_Bcast(recvcount+i,1,MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
      /* we could also use MPI_Send/Recv to get everything just to the first processor */
      ierr = MPI_Bcast(p,recvcount[i],MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
      p+=recvcount[i];
    }
    assert(p==elevertall+NV*gdata->n_ele);
  }
  else {
    elevertall=gdata->elevert;
  }
#endif

  int *elepropall;
  if (size>1) {
    ierr = PetscMalloc(gdata->n_ele*sizeof(int),&elepropall);CHKERRQ(ierr);
    ierr = PetscMemcpy(elepropall+low,gdata->eleprop,gdata->ln_ele*sizeof(int));CHKERRQ(ierr);

    int *p;
    p=elepropall;

    int recvcount[size];
    recvcount[rank]=gdata->ln_ele;
    for (int i=0;i<size;i++) {
      ierr = MPI_Bcast(recvcount+i,1,MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);

      ierr = MPI_Bcast(p,recvcount[i],MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
      p+=recvcount[i];
    }
    assert(p==elepropall+gdata->n_ele);
  }
  else {
    elepropall=gdata->eleprop;
  }

  int *elevertcoupled;
  ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&elevertcoupled);CHKERRQ(ierr);
  ierr = PetscMemcpy(elevertcoupled,elevertall,NV*gdata->n_ele*sizeof(int));CHKERRQ(ierr);

#ifndef HAVE_ELEVERTALL
  if (size>1) {
    ierr = PetscFree(elevertall);CHKERRQ(ierr);
  }
#endif

  PetscReal sarea_area;
  sarea_area=0;

  PetscReal sarea_point[ND];
  PetscReal sarea_norm[ND];
  for (int i=0;i<ND;i++) {
    sarea_point[i]=sarea_norm[i]=0.0;
  }

  /* only first processor does everything */
  if (!rank) {

#ifdef ADDONS
    /* map all coupled vertices to common vertex (lowest id) */
    if (gdata->vertcoupl!=PETSC_NULL) {
      for (int i=0;i<NV*gdata->n_ele;i++) {
        if (gdata->vertcoupl[elevertcoupled[i]]>=0 &&
            gdata->vertcoupl[elevertcoupled[i]]<elevertcoupled[i]) {
          elevertcoupled[i]=gdata->vertcoupl[elevertcoupled[i]];
        }
      }
    }
#endif

    PetscPrintf(PETSC_COMM_WORLD,"Calculating mesh connectivity\n");

    int *ia,*ja;
    ierr = Mesh2Dual(gdata->n_ele,gdata->n_vert,elevertcoupled,&ia,&ja);CHKERRQ(ierr);

    for (int i=1;i<gdata->n_ele+1;i++) {
      if (ia[i]-ia[i-1]>NF) {
        SETERRQ2(PETSC_ERR_MEMC,"Found element with %i neighbors (>NF)!\n",ia[i]-ia[i-1],NF);
      }
    }

    /* check for array overrun */
    assert(ia[gdata->n_ele]<NF*gdata->n_ele);

    PetscPrintf(PETSC_COMM_WORLD,"Searching shared faces\n");
    for (int i=0; i<gdata->n_ele; i++) {
      if (elepropall[i]!=pid1) continue;

      if (pid2>=0) {
        for (int j2=ia[i];j2<ia[i+1];j2++) {
          int j;
          j=ja[j2];

          if (elepropall[j]!=pid2) continue;

          int s[NV-1];
          int k;
          k=TetSharedVertices(elevertcoupled+NV*i,elevertcoupled+NV*j,s);
          assert(k==NV-1);

          PetscReal a;
          PetscReal n[ND];
          a=triarea(
            gdata->vertxyz+ND*s[0],
            gdata->vertxyz+ND*s[1],
            gdata->vertxyz+ND*s[2],
            n
          );

          sarea_area+=a;

          for (int p=0;p<ND;p++) {
            if (n[2]>=0.0) {
              sarea_norm[p]+=n[p]*a;
            }
            else {
              sarea_norm[p]-=n[p]*a;
            }
            for (int q=0;q<NV-1;q++) {
              sarea_point[p]+=gdata->vertxyz[ND*s[q]+p]*a/(NV-1.0);
            }
          }
        }
      }
      else if (pid2==-1) {
        /* total surface of vol with pid1: complicated to implement (see facnb.c) */
        SETERRQ1(PETSC_ERR_ARG_INCOMP,"Program branch for pid2==%i not implemented!\n",pid2);
      }
      else if (pid2==-2) {
        /* free surface of vol with pid1: complicated to implement (see facnb.c) */
        SETERRQ1(PETSC_ERR_ARG_INCOMP,"Program branch for pid2==%i not implemented!\n",pid2);
      }
    }

    ierr = PetscFree(ia);CHKERRQ(ierr);
    ierr = PetscFree(ja);CHKERRQ(ierr);

  }

  ierr = MPI_Bcast(&sarea_area,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  ierr = MPI_Bcast(sarea_point,ND,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  ierr = MPI_Bcast(sarea_norm,ND,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);

  if (size>1) {
    ierr = PetscFree(elepropall);CHKERRQ(ierr);
  }

  *sharedarea_area=sarea_area;

  if (sharedarea_point!=PETSC_NULL) {
    my_dscal(ND,1.0/sarea_area,sarea_point,1);
    my_dcopy(ND,sarea_point,1,sharedarea_point,1);
  }

  if (sharedarea_norm!=PETSC_NULL) {
    my_dscal(ND,1.0/sarea_area,sarea_norm,1);
    my_dcopy(ND,sarea_norm,1,sharedarea_norm,1);
  }

  MagparFunctionLogReturn(0);
}