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

#include "util.h"

int CalAsq2fe(GridData *gdata,PetscReal *vertxyz2,PetscReal *bbox,int *pix,
              int slice_id,Mat *pAIpol)
{
  int         i,j,k,l;
  Mat         AIpol;


  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 sq->fe 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]
  );

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

  for (l=0;l<gdata->ln_vert;l++) {
    PetscReal   *v, p[ND];

    if (slice_id>=0 && gdata->vertprop[l] != slice_id) continue;

    /* calculate which pixel the x2 and y2 coordinates of the
       vertex correspond to; z2 value neglected -
       thus all nodes are effectively projected into the x2-y2 plane!
    */
    v=vertxyz2+ND*gdata->vertl2g[l];
    i=(int) floor(pix[0]*(v[0]-bbox[0])/(bbox[ND+0]-bbox[0]));
    j=(int) floor(pix[1]*(v[1]-bbox[1])/(bbox[ND+1]-bbox[1]));

    PetscReal a,b;
    a=b=-2;
    if (i<0) {i=0;a=-1;}
    if (j<0) {j=0;b=-1;}
    if (i>=pix[0]) {i=pix[0]-1;a=+1;}
    if (j>=pix[1]) {j=pix[1]-1;b=+1;}

    /* 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]=0;

    /* calculate weight factors unless they were set above */
    if (a<-1.5) a=pix[0]*(v[0]-p[0])/(bbox[ND+0]-bbox[0])*2.0-1.0;
    if (b<-1.5) b=pix[1]*(v[1]-p[1])/(bbox[ND+1]-bbox[1])*2.0-1.0;

    for (k=0;k<ND;k++) {
      PetscReal   matele;

      matele=(1-a)*(1-b)/4.0;
      ierr = MatSetValue(AIpol,
        ND*gdata->vertl2g[l]+k,
        ND*((i+0)*pix[1]+(j+0))+k,
      matele,ADD_VALUES);CHKERRQ(ierr);

      matele=(1-a)*(1+b)/4.0;
      ierr = MatSetValue(AIpol,
        ND*gdata->vertl2g[l]+k,
        ND*((i+0)*pix[1]+(j+1))+k,
      matele,ADD_VALUES);CHKERRQ(ierr);

      matele=(1+a)*(1-b)/4.0;
      ierr = MatSetValue(AIpol,
        ND*gdata->vertl2g[l]+k,
        ND*((i+1)*pix[1]+(j+0))+k,
      matele,ADD_VALUES);CHKERRQ(ierr);

      matele=(1+a)*(1+b)/4.0;
      ierr = MatSetValue(AIpol,
        ND*gdata->vertl2g[l]+k,
        ND*((i+1)*pix[1]+(j+1))+k,
      matele,ADD_VALUES);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);


  /* TODO: we should not miss any nodes any more
     remove this section (down to the assert) if this holds true
  */
  PetscPrintf(PETSC_COMM_WORLD,"Checking for missing nodes\n");

  Vec       test1,test2;
  int       cnt=0;

  ierr = VecDuplicate(gdata->M,&test2);
  ierr = VecCreate(PETSC_COMM_WORLD,&test1);CHKERRQ(ierr);
  if (!rank) {
    ierr = VecSetSizes(test1,ND*pix[0]*pix[1],ND*pix[0]*pix[1]);CHKERRQ(ierr);
  }
  else{
    ierr = VecSetSizes(test1,0,ND*pix[0]*pix[1]);CHKERRQ(ierr);
  }
  ierr = VecSetFromOptions(test1);CHKERRQ(ierr);

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

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

  PetscReal *atest2;
  ierr = VecGetArray(test2,&atest2);CHKERRQ(ierr);

  for (l=0;l<gdata->ln_vert;l++) {
    PetscReal   *v;

    if (slice_id>=0 && gdata->vertprop[l] != slice_id) continue;
    if (atest2[ND*l]>0.5) continue;

    /* calculate which pixel the x2 and y2 coordinates of the
       vertex correspond to; z2 value neglected -
       thus all nodes are effectively projected into the x2-y2 plane!
    */
    v=vertxyz2+ND*gdata->vertl2g[l];

    i=(int) floor(pix[0]*(v[0]-bbox[0])/(bbox[ND+0]-bbox[0]));
    j=(int) floor(pix[1]*(v[1]-bbox[1])/(bbox[ND+1]-bbox[1]));

/*
    PetscPrintf(PETSC_COMM_WORLD,
      "missed node %i (%g,%g,%g), would need pixel %i %i, matele %g\n",
      gdata->vertl2g[l],v[0],v[1],v[2],i,j,atest2[ND*l]
    );
*/

    if (i<0) i=0;
    if (j<0) j=0;
    if (i>=pix[0]) i=pix[0]-1;
    if (j>=pix[1]) j=pix[1]-1;

    for (k=0;k<ND;k++) {
      ierr = MatSetValue(AIpol,
        ND*gdata->vertl2g[l]+k,
        ND*((i+1)*pix[1]+(j+1))+k,
      1.0,ADD_VALUES);CHKERRQ(ierr);
    }
    cnt++;
  }

  ierr = VecRestoreArray(test2,&atest2);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);

  PetscSynchronizedPrintf(PETSC_COMM_WORLD,
    "<%i>Added %i nodes\n",
    rank,cnt
  );
  PetscSynchronizedFlush(PETSC_COMM_WORLD);

  assert(cnt==0);


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

  *pAIpol=AIpol;

  MagparFunctionLogReturn(0);
}