calAfe2sq.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2009 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: calAfe2sq.c 2741 2009-08-26 18:00:35Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "util.h"

int findeleslice(int ln_ele,int *elevert,PetscReal *vertxyz2,PetscReal zcut,int *nele,int *sliceele)
{

  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 cnt=0;
  for (int i=0;i<ln_ele;i++) {
    PetscReal sgn;
    /* get z-coordinate of first vertex */
    sgn=vertxyz2[ND*elevert[NV*i+0]+2]-zcut;

    for (int j=1;j<NV;j++) {
      /* multiply with other z-coordinate of other vertices
         if (product < 0) then they have opposite signs and
         the element must be cut by the slice plane
      */
      if (sgn*(vertxyz2[ND*elevert[NV*i+j]+2]-zcut) <= 0.0) {
        /* add its local index to the array */
        sliceele[cnt++]=i;
        break;
      }
    }
  }

  PetscSynchronizedPrintf(PETSC_COMM_WORLD,
    "<%i>Found %i elements cut by slice plane\n",
    rank,cnt
  );
  PetscSynchronizedFlush(PETSC_COMM_WORLD);

  *nele=cnt;

  MagparFunctionLogReturn(0);
}

int CalAfe2sq(GridData *gdata,PetscReal *vertxyz2,PetscReal *bbox,int *pix,int slice_id,Mat *pAIpol)
{
  MagparFunctionLogBegin;

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

  PetscPrintf(PETSC_COMM_WORLD,
    "Calculating interpolation matrix fe->sq for volume with pid %i\n"
    "in bbox (%g,%g,%g) (%g,%g,%g) with %i x %i x %i pixels\n",
    slice_id+1,
    bbox[0],bbox[1],bbox[2],bbox[ND+0],bbox[ND+1],bbox[ND+2],
    pix[0],pix[1],pix[2]
  );

  /* find elements which cut the slice plane
     (// to x-y plane: z=bbox[ND+2]-bbox[2])
  */
  int *sliceele;
  ierr = PetscMalloc(gdata->ln_ele*sizeof(int),&sliceele);CHKERRQ(ierr);
  int cnt;
  ierr = findeleslice(gdata->ln_ele,gdata->elevert,vertxyz2,(bbox[2]+bbox[ND+2])/2.0,&cnt,sliceele);CHKERRQ(ierr);

  int sum;
  ierr=MPI_Allreduce((int*)&cnt,(int*)&sum,1,MPI_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,
    "Found total %i elements cut by slice plane.\n",
    sum
  );

  if (sum==0) {
    PetscPrintf(PETSC_COMM_WORLD,"Slice plane outside the FE mesh.\n");

    ierr = PetscFree(sliceele);CHKERRQ(ierr);
    *pAIpol=PETSC_NULL;

    MagparFunctionLogReturn(0);
  }

  Mat AIpol;
  if (!rank) {
    ierr = MatCreateMPIAIJ(
      PETSC_COMM_WORLD,
      ND*pix[0]*pix[1],ND*gdata->ln_vert,
      ND*pix[0]*pix[1],ND*gdata->n_vert,
      2*NV,PETSC_NULL,2*NV,PETSC_NULL,
      &AIpol
    );CHKERRQ(ierr);
  }
  else {
    ierr = MatCreateMPIAIJ(
      PETSC_COMM_WORLD,
      0,ND*gdata->ln_vert,
      ND*pix[0]*pix[1],ND*gdata->n_vert,
      NV,PETSC_NULL,NV,PETSC_NULL,
      &AIpol
    );CHKERRQ(ierr);
  }
  ierr = MatSetFromOptions(AIpol);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,"Calculating interpolation matrix...\n");
  /* loop over all elements, which are cut by the slice plane */
  for (int k=0; k<cnt; k++) {
    PetscReal elebox[2*ND];
    int       pixbox[2*ND];
    int       t_ele;

    ierr = ProgressBar(k,cnt,10);

    /* get current element id */
    t_ele=sliceele[k];

    /* check for requested property id
     * slice_id >= 0:  only if pid==slice_id
     * slice_id == -1: any pid
     * slice_id <= -2: only if pid!=slice_id
     */
    if (slice_id>=0 && gdata->eleprop[t_ele]!=slice_id) continue;
    if (slice_id<=-2 && gdata->eleprop[t_ele]==-slice_id-1) continue;

    /* find bbox around element */
    for (int i=0;i<ND;i++) {
      elebox[i]=PETSC_MAX;
      elebox[ND+i]=PETSC_MIN;
    }
    for (int l=0; l<NV; l++) {
      for (int i=0;i<ND;i++) {
        elebox[i]=   PetscMin(elebox[i],   vertxyz2[ND*gdata->elevert[NV*t_ele+l]+i]);
        elebox[ND+i]=PetscMax(elebox[ND+i],vertxyz2[ND*gdata->elevert[NV*t_ele+l]+i]);
      }
    }
    /* check that bbox is valid */
    for (int i=0;i<ND;i++) {
      assert(elebox[ND+i]>elebox[i]);
    }

    /* find pixels in the elebox */
    for (int i=0;i<ND;i++) {
      pixbox[i]=   (int) floor(pix[i]*(elebox[   i]-bbox[i])/(bbox[ND+i]-bbox[i]));
      pixbox[ND+i]=(int)  ceil(pix[i]*(elebox[ND+i]-bbox[i])/(bbox[ND+i]-bbox[i]));
    }

    /* loop over pixels in pixbox */
    for (int i=pixbox[0];i<=pixbox[ND+0];i++) {
      /* check if we are outside the requested rectangle */
      if (i<0 || i>=pix[0]) continue;

      for (int j=pixbox[1];j<=pixbox[ND+1];j++) {
        /* check if we are outside the requested rectangle */
        if (j<0 || j>=pix[1]) continue;

        PetscReal p[ND],rhs[ND];

        /* calculate xyz coordinates of pixel
           NB: keep integers (i,pix[0]) separated by float
           otherwise division of integers always gives 0!
        */
        p[0]=bbox[0]+i*(bbox[ND+0]-bbox[0])/pix[0];
        p[1]=bbox[1]+j*(bbox[ND+1]-bbox[1])/pix[1];
        p[2]=(bbox[2]+bbox[ND+2])/2.0;

        if (
          barycent(
            p,
            vertxyz2+ND*gdata->elevert[NV*t_ele+0],
            vertxyz2+ND*gdata->elevert[NV*t_ele+1],
            vertxyz2+ND*gdata->elevert[NV*t_ele+2],
            vertxyz2+ND*gdata->elevert[NV*t_ele+3],
            rhs
          )
        ) {
          for (int l=0;l<ND;l++) {
            PetscReal   matele;

            matele=(1-rhs[0]-rhs[1]-rhs[2]);
            ierr = MatSetValue(AIpol,
              ND*(i*pix[1]+j)+l,
              ND*gdata->elevert[NV*t_ele+3]+l,
            matele,ADD_VALUES);CHKERRQ(ierr);

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

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

            matele=rhs[2];
            ierr = MatSetValue(AIpol,
              ND*(i*pix[1]+j)+l,
              ND*gdata->elevert[NV*t_ele+2]+l,
            matele,ADD_VALUES);CHKERRQ(ierr);
          }
        }
      }
    }
  }
  ierr = ProgressBar(1,1,10);

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

  PetscInfo(0,"AIpol matrix assembly complete\n");
  ierr = MatAssemblyBegin(AIpol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(AIpol,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = PrintMatInfoAll(AIpol);CHKERRQ(ierr);

  /* check interpolation matrix */
  PetscReal area;
  ierr = AreaCal(AIpol,&area);CHKERRQ(ierr);

  *pAIpol=AIpol;

  MagparFunctionLogReturn(0);
}