calAfe2fe.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: calAfe2fe.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 CalAfe2fe(GridData *gdata,PetscReal *shift,int *slice_id,Mat *pAIpol)
{
  MagparFunctionLogBegin;

  PetscPrintf(PETSC_COMM_WORLD,"Calculating interpolation matrix fe->fe\n");

  Mat AIpol;
  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ND*gdata->ln_vert,ND*gdata->ln_vert,
    ND*gdata->n_vert,ND*gdata->n_vert,
    NV,PETSC_NULL,NV,PETSC_NULL,
    &AIpol
  );CHKERRQ(ierr);
  ierr = MatSetFromOptions(AIpol);CHKERRQ(ierr);

  /* Set diagonal matrix elements to 1 to make sure we do not interpolate any vectors to 0.
     This also takes care of the boundary nodes, which "fall off the map".
     Though, if the shift is too large (node not within neighboring elements)
     it also gets interpolated to itself (i.e. not shifted)!
  */
  for (int i=0;i<gdata->ln_vert;i++) {
    for (int k=0;k<ND;k++) {
      ierr = MatSetValue(
        AIpol,
        ND*gdata->vertl2g[i]+k,
        ND*gdata->vertl2g[i]+k,
        1.0,
        INSERT_VALUES
      );CHKERRQ(ierr);
    }
  }

#define DIST 1e-12
  PetscReal tshift[ND];
  my_dcopy(ND,shift,1,tshift,1);
  my_dscal(ND,1+DIST,tshift,1);
  tshift[0]+=DIST;
  tshift[1]+=DIST;
  tshift[2]+=DIST;
  PetscPrintf(PETSC_COMM_WORLD,"Shifting nodes by (%g,%g,%g)\n",tshift[0],tshift[1],tshift[2]);

  PetscPrintf(PETSC_COMM_WORLD,"Calculating interpolation matrix...\n");
  for (int i=0;i<gdata->ln_ele;i++) {
    ierr = ProgressBar(i,gdata->ln_ele,10);

    int j;
    for (j=0;j<gdata->n_prop;j++) {
      if (gdata->eleprop[i]==slice_id[j]) break;
    }
    if (j==gdata->n_prop) continue;

    for (int k=0;k<NV;k++) {
      int v;
      v=gdata->elevert[NV*i+k];

      PetscReal p[ND],rhs[ND];

      /* copy Cartesian coordinates */
      for (int m=0;m<ND;m++) {
        p[m]=gdata->vertxyz[ND*v+m];
        p[m]+=tshift[m];
      }

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

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

          matele=rhs[0];
          ierr = MatSetValue(AIpol,
            ND*v+l,
            ND*gdata->elevert[NV*i+0]+l,
          matele,INSERT_VALUES);CHKERRQ(ierr);

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

          matele=rhs[2];
          ierr = MatSetValue(AIpol,
            ND*v+l,
            ND*gdata->elevert[NV*i+2]+l,
          matele,INSERT_VALUES);CHKERRQ(ierr);
        }
      }
    }
  }

  ierr = ProgressBar(1,1,10);

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

#if 0
  /* This is taken care of by setting the diagonal to 1 above. */

  /* Make sure all nodes are getting interpolated */
  /* Since all nodes are shifted only a small amount so they are within an element they already
   * belong to, there has to be a non-zero diagonal matrix element. */
  Vec Vdiag;
  ierr = VecDuplicate(gdata->M,&Vdiag);CHKERRQ(ierr);

  /* This could be rather slow for big matrices: */
  ierr = MatGetDiagonal(AIpol,Vdiag);CHKERRQ(ierr);

  PetscReal *d;
  ierr = VecGetArray(Vdiag,&d);CHKERRQ(ierr);
  for (int i=0;i<ND*gdata->ln_vert;i++) {
    if (d[gdata->vertl2g[i]]==0.0) d[gdata->vertl2g[i]]=1.0;
  }
  ierr = VecRestoreArray(Vdiag,&d);CHKERRQ(ierr);

  ierr = MatDiagonalSet(AIpol,Vdiag,INSERT_VALUES);CHKERRQ(ierr);
#endif

  ierr = PrintMatInfoAll(AIpol);CHKERRQ(ierr);

/* DEBUG
  PetscPrintf(PETSC_COMM_WORLD,"deb01: AIpol\n");
  ierr = MatView(AIpol,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
*/

  *pAIpol=AIpol;

  MagparFunctionLogReturn(0);
}