Commits

Rio Yokota committed 2bf9469

3-D version by Huda.

Comments (0)

Files changed (12)

 #include "petiga.h"
-#include "fmm.h"
+#include "fmm2D.h"
 
 EXTERN_C_BEGIN
 extern PetscErrorCode PCCreate_FMM(PC);
   PetscInt N[2] = {16,16}, nN = 2; 
   PetscInt p[2] = { 1, 1}, np = 2;
   PetscInt C[2] = {-1,-1}, nC = 2;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson2D Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson Options","IGA");CHKERRQ(ierr);
   ierr = PetscOptionsIntArray("-N","number of elements",     __FILE__,N,&nN,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsIntArray("-p","polynomial order",       __FILE__,p,&np,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsIntArray("-C","global continuity order",__FILE__,C,&nC,PETSC_NULL);CHKERRQ(ierr);
+#include "petiga.h"
+#include "fmm3D.h"
+
+EXTERN_C_BEGIN
+extern PetscErrorCode PCCreate_FMM(PC);
+EXTERN_C_END
+
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt dim = p->dim;
+  PetscInt nen = p->nen;
+
+  PetscReal *N0 = p->shape[0];
+  PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
+
+  PetscInt a,b;
+  for (a=0; a<nen; a++) {
+    for (b=0; b<nen; b++)
+      K[a*nen+b] = DOT(dim,N1[a],N1[b]);
+    F[a] = N0[a] * 1.0;
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  PetscInt dim = 3;
+
+  PetscInt N[3] = {16,16,16}, nN = 3; 
+  PetscInt p[3] = { 1, 1, 1}, np = 3;
+  PetscInt C[3] = {-1,-1,-1}, nC = 3;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsInt     ("-dim","dimension",            __FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray("-N","number of elements",     __FILE__,N,&nN,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray("-p","polynomial order",       __FILE__,p,&np,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray("-C","global continuity order",__FILE__,C,&nC,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+  if (nN == 1) N[2] = N[1] = N[0]; else if (nN == 2) N[2] = N[0];
+  if (np == 1) p[2] = p[1] = p[0]; else if (np == 2) p[2] = p[0];
+  if (nC == 1) C[2] = C[1] = C[0]; else if (nC == 2) C[2] = C[0];
+  if (C[0] == -1) C[0] = p[0]-1;
+  if (C[1] == -1) C[1] = p[1]-1;
+  if (C[2] == -1) C[2] = p[2]-1;
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+
+  PetscInt i;
+  for (i=0; i<dim; i++) {
+    IGAAxis axis;
+    ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
+    ierr = IGAAxisSetDegree(axis,p[i]);CHKERRQ(ierr);
+    ierr = IGAAxisInitUniform(axis,N[i],0.0,1.0,C[i]);CHKERRQ(ierr);
+  }
+  IGABoundary bnd;
+  PetscInt dir,side;
+  PetscScalar value = 0.0;
+  for (dir=0; dir<dim; dir++) {
+    for (side=0; side<2; side++) {
+      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+      ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
+    }
+  }
+
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  Mat A;
+  Vec x,b;
+  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
+  ierr = IGASetUserSystem(iga,System,PETSC_NULL);CHKERRQ(ierr);
+  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+  
+  KSP ksp;
+  ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
+
+#if 1
+  PC  pc;
+  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
+  ierr = PCRegister(PCFMM,PCCreate_FMM);CHKERRQ(ierr);
+  ierr = PCSetType(pc,PCFMM);CHKERRQ(ierr);
+  ierr = PCFMMSetIGA(pc,iga);CHKERRQ(ierr);
+#endif
+
+  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+  
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = MatDestroy(&A);CHKERRQ(ierr);
+  ierr = VecDestroy(&x);CHKERRQ(ierr);
+  ierr = VecDestroy(&b);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga);CHKERRQ(ierr);
+
+  PetscBool flag = PETSC_FALSE;
+  PetscReal secs = -1;
+  ierr = PetscOptionsHasName(0,"-sleep",&flag);CHKERRQ(ierr);
+  ierr = PetscOptionsGetReal(0,"-sleep",&secs,0);CHKERRQ(ierr);
+  if (flag) {ierr = PetscSleep(secs);CHKERRQ(ierr);}
+
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}

fmm/fmm.c

-#include <assert.h>
-#include "petiga.h"
-#include "fmm.h"
-#include "nbody.h"
-#if PETSC_VERSION_LE(3,2,0)
-#include <private/pcimpl.h>
-#else
-#include <petsc-private/pcimpl.h>
-#endif
-
-#define DENSE 0
-
-#undef  __FUNCT__
-#define __FUNCT__ "PCFMMSetIGA"
-PetscErrorCode PCFMMSetIGA(PC pc, IGA iga)
-{
-  PC_FMM *fmm = (PC_FMM*)pc->data;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,2);
-  if (iga) fmm->iga = iga;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "BEM_MatMult"
-static PetscErrorCode BEM_MatMult(Mat G, Vec x, Vec y)
-{
-  PC_FMM *fmm;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  ierr = MatShellGetContext(G,(void**)&fmm);CHKERRQ(ierr);
-  PetscInt i, nb = fmm->nb;
-  PetscReal *xl, *yl;
-  PetscReal *xm = fmm->xm;
-  PetscReal *ym = fmm->ym;
-  PetscReal *re = fmm->re;
-  PetscReal eps = fmm->eps;
-  VecGetSize(x,&i);
-  ierr = VecGetArray(x,&xl);CHKERRQ(ierr);
-  ierr = VecGetArray(y,&yl);CHKERRQ(ierr);
-  for (i=0; i<nb; i++) {
-    yl[i] = 0;
-  }
-  nbodyG(nb,xm,ym,yl,nb,xm,ym,xl,eps);
-  for (i=0; i<nb; i++) {
-    yl[i] += (1 - log(re[i] / 2) + log(eps)) * xl[i];
-  }
-  ierr = VecRestoreArray(x,&xl);CHKERRQ(ierr);
-  ierr = VecRestoreArray(y,&yl);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "PCSetUp_FMM"
-static PetscErrorCode PCSetUp_FMM(PC pc)
-{
-  PC_FMM *fmm = (PC_FMM*)pc->data;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  if (!fmm->iga) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCFMMSetIGA() must be called first");
-  ierr = IGACreateNodeDM(fmm->iga,1,&(fmm->da));CHKERRQ(ierr);
-  DMDALocalInfo info;
-  ierr = DMDAGetLocalInfo(fmm->da,&info);CHKERRQ(ierr);
-  fmm->nb = 0;
-  if(info.xs == 0){
-    fmm->nb += info.ym;
-    if(info.ys == 0)               fmm->nb -= 1;
-    if(info.my == info.ys+info.ym) fmm->nb -= 1;
-  }
-  if(info.mx == info.xs+info.xm){
-    fmm->nb += info.ym;
-    if(info.ys == 0)               fmm->nb -= 1;
-    if(info.my == info.ys+info.ym) fmm->nb -= 1;
-  }
-  if(info.ys == 0              ) fmm->nb += info.xm;
-  if(info.my == info.ys+info.ym) fmm->nb += info.xm;
-  fmm->ni = info.xm*info.ym - fmm->nb;
-
-  ierr = PetscMalloc(fmm->nb*sizeof(PetscInt),&(fmm->ip1));CHKERRQ(ierr);
-  ierr = PetscMalloc3(fmm->ni,PetscReal,&(fmm->xi),
-		      fmm->ni,PetscReal,&(fmm->yi),
-		      fmm->ni,PetscReal,&(fmm->ri));CHKERRQ(ierr);
-  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->xb),
-		      fmm->nb,PetscReal,&(fmm->yb),
-		      fmm->nb,PetscReal,&(fmm->ub));CHKERRQ(ierr);
-  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->xm),
-		      fmm->nb,PetscReal,&(fmm->ym),
-		      fmm->nb,PetscReal,&(fmm->um));CHKERRQ(ierr);
-  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->dxe),
-		      fmm->nb,PetscReal,&(fmm->dye),
-		      fmm->nb,PetscReal,&(fmm->re));CHKERRQ(ierr);
-  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->rhs),
-		      fmm->nb,PetscReal,&(fmm->un),
-		      fmm->ni,PetscReal,&(fmm->ui));CHKERRQ(ierr);
-  ierr = VecCreateMPI(PETSC_COMM_WORLD,fmm->nb,PETSC_DETERMINE,&(fmm->u));CHKERRQ(ierr);
-  ierr = VecDuplicate(fmm->u,&(fmm->f));CHKERRQ(ierr);
-  ierr = MatCreate(PETSC_COMM_WORLD,&fmm->G);CHKERRQ(ierr);
-  ierr = MatSetSizes(fmm->G,fmm->nb,fmm->nb,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
-#if DENSE
-  ierr = MatSetType(fmm->G,MATDENSE);CHKERRQ(ierr);
-#else
-  ierr = MatSetType(fmm->G,MATSHELL);CHKERRQ(ierr);
-  ierr = MatShellSetOperation(fmm->G,MATOP_MULT,(void(*)(void))BEM_MatMult);CHKERRQ(ierr);
-  ierr = MatShellSetContext(fmm->G,fmm);CHKERRQ(ierr);
-#endif
-  ierr = MatSetFromOptions(fmm->G);CHKERRQ(ierr);
-  ierr = MatSetUp(fmm->G);CHKERRQ(ierr);
-
-  PetscReal eps = .2 / (1<<((int)(log2(info.mx-1))));
-  FMM_Init(eps*eps, 0);
-  fmm->eps = eps;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "PCApply_FMM"
-static PetscErrorCode PCApply_FMM(PC pc, Vec x, Vec y)
-{
-  PC_FMM *fmm = (PC_FMM*)pc->data;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-
-  DM da = fmm->da;
-
-  PetscScalar **X,**Y;
-  ierr = DMDAVecGetArray(da,x,&X);CHKERRQ(ierr);
-  ierr = DMDAVecGetArray(da,y,&Y);CHKERRQ(ierr);
-  DMDALocalInfo info;
-  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
-
-  Vec u = fmm->u;
-  Vec f = fmm->f;
-  Mat G = fmm->G;
-  IGA iga = fmm->iga;
-  PetscInt ni = fmm->ni;
-  PetscInt nb = fmm->nb;
-  PetscReal *xi  = fmm->xi;
-  PetscReal *yi  = fmm->yi;
-  PetscReal *ri  = fmm->ri;
-  //PetscReal *ub  = fmm->ub;
-  PetscReal *xm  = fmm->xm;
-  PetscReal *ym  = fmm->ym;
-  PetscReal *um  = fmm->um;
-  PetscReal *dxe = fmm->dxe;
-  PetscReal *dye = fmm->dye;
-  PetscReal *re  = fmm->re;
-  PetscReal *rhs = fmm->rhs;
-  PetscReal *un  = fmm->un;
-  PetscReal *ui  = fmm->ui;
-  PetscReal  eps = fmm->eps;
-
-  PetscInt i, j, ic = 0;
-  for (i=info.xs; i<info.xs+info.xm; i++) { // Loop over interior points
-    for (j=info.ys; j<info.ys+info.ym; j++) {
-      if (i != 0 && i != info.mx-1 && j != 0 && j != info.my-1) {
-	xi[ic] = (PetscReal)i / ( (PetscReal)(info.mx-1) );
-	yi[ic] = (PetscReal)j / ( (PetscReal)(info.my-1) );
-	ri[ic] = X[j][i] / 2 / M_PI;
-	ic++;
-      }
-    }
-  }
-  assert( ni == ic );
-
-  ic = 0;
-  j = 0;
-  PetscInt q;
-  PetscInt *start = iga->elem_start;
-  PetscInt *width = iga->elem_width;
-  PetscInt *sizes = iga->elem_sizes;
-  IGABasis *BD    = iga->basis;
-  PetscReal val_dirichlet = 0;
-  PetscReal h     =  BD[0]->detJ[0]*2;
-  for (j=0; j<2; j++){ // loop bottom and top
-    if (start[1]+j*width[1] == j*sizes[1]) { // domain boundary?
-      for (i=start[0]; i<start[0]+width[0]; i++) {
-	for (q=0; q<1; q++) {
-	  xm[ic]  =  (((PetscReal)i)+0.5)*h; //BD[0]->point[i*BD[0]->nqp+q];
-	  ym[ic]  =  j*1;
-	  um[ic]  = -val_dirichlet*0.5/M_PI;
-	  dxe[ic] =  BD[0]->detJ[i]*2 *(j == 0 ? 1 : -1);
-	  dye[ic] =  0;
-	  re[ic]  =  fabs(dxe[ic]);
-	  rhs[ic] =  um[ic]/re[ic]*M_PI;
-	  //printf("%d%d %f %f %f %f %f %f %f\n",start[0],start[1],xm[ic],ym[ic],um[ic],dxe[ic],dye[ic],re[ic],rhs[ic]);
-	  ic += 1;
-	}
-      }
-    }
-  }
-  for (j=0; j<2; j++){ // loop left and right
-    if (start[0]+j*width[0] == j*sizes[0]) { // domain boundary?
-      for (i=start[1]; i<start[1]+width[1]; i++) {
-	for (q=0; q<1; q++) {
-	  xm[ic]  =  j*1;
-	  ym[ic]  =  (((PetscReal)i)+0.5)*h; //BD[1]->point[i*BD[1]->nqp+q];
-	  um[ic]  = -val_dirichlet*0.5/M_PI;
-	  dxe[ic] =  0;
-	  dye[ic] =  BD[1]->detJ[i]*2 *(j == 0 ? -1 : 1);
-	  re[ic]  =  fabs(dye[ic]);
-	  rhs[ic] =  um[ic]/re[ic]*M_PI;
-	  //printf("%d%d %f %f %f %f %f %f %f\n",start[0],start[1],xm[ic],ym[ic],um[ic],dxe[ic],dye[ic],re[ic],rhs[ic]);
-	  ic += 1;
-	}
-      }
-    }
-  }
-
-  nb = ic; // this is ok, but we should improve memory allocation
-  assert( nb == ic );
-
-  //FMM(nb,xm,ym,rhs,ni,xi,yi,ri);
-  nbodyG(nb,xm,ym,rhs,ni,xi,yi,ri,eps);
-  //nbodyGn(nb,xm,ym,rhs,nb,xm,ym,um,dxe,dye,re,eps);
-
-  for (i=0; i<nb; i++) {
-    VecSetValue(f,i,rhs[i],INSERT_VALUES);
-#if DENSE
-    for (j=0; j<i; j++) {
-      PetscReal dx = xm[i] - xm[j];
-      PetscReal dy = ym[i] - ym[j];
-      PetscReal r = sqrt(dx * dx + dy * dy + eps * eps);
-      PetscReal value = -log(r);
-      MatSetValue(G,i,j,value,INSERT_VALUES);
-      MatSetValue(G,j,i,value,INSERT_VALUES);
-    }
-    PetscReal value = 1-log(re[i]/2);
-    MatSetValue(G,i,i,value,INSERT_VALUES);
-#endif
-  }
-  ierr = VecAssemblyBegin(f);CHKERRQ(ierr);
-  ierr = VecAssemblyEnd(f);CHKERRQ(ierr);
-#if DENSE
-  ierr = MatAssemblyBegin(G,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-  ierr = MatAssemblyEnd(G,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-#endif
-
-  KSP ksp;
-  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
-  ierr = KSPSetOperators(ksp,G,G,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
-  ierr = KSPSetOptionsPrefix(ksp,"pc_fmm_");CHKERRQ(ierr);
-  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
-
-  ierr = KSPSolve(ksp,f,u);CHKERRQ(ierr);
-  ierr = VecGetArray(u,&un);CHKERRQ(ierr);
-
-  for (i=0; i<nb; i++) {
-    un[i] = -un[i];
-  }
-
-  for (i=0; i<ni; i++) {
-    ui[i] = 0;
-  }
-  //FMM(ni,xi,yi,ui,ni,xi,yi,ri);
-  //FMM(ni,xi,yi,ui,nb,xm,ym,un);
-  nbodyG(ni,xi,yi,ui,ni,xi,yi,ri,eps);
-  nbodyG(ni,xi,yi,ui,nb,xm,ym,un,eps);
-  //nbodyGn(ni,xi,yi,ui,nb,xm,ym,um,dxe,dye,re,eps);
-
-  ic = 0;
-  for (i=info.xs; i<info.xs+info.xm; i++) { // Loop over interior points
-    for (j=info.ys; j<info.ys+info.ym; j++) {
-      if (i != 0 && i != info.mx-1 && j != 0 && j != info.my-1) {
-        Y[j][i] = ui[ic];
-        ic++;
-      } else {
-        Y[j][i] = 0;
-      }
-    }
-  }
-
-#if 0
-  ic = 0;
-  j = 0;
-  for (i=info.xs; i<info.xs+info.xm-1; i++,ic++) { // Bottom points
-    Y[j][i] = ub[ic];
-  }
-  i = info.xs+info.xm-1;
-  for (j=info.ys; j<info.ys+info.ym-1; j++,ic++) { // Right points
-    Y[j][i] = ub[ic];
-  }
-  j = info.ys+info.ym-1;
-  for (i=info.xs+info.xm-1; i>info.xs; i--,ic++) { // Top points
-    Y[j][i] = ub[ic];
-  }
-  i = 0;
-  for (j=info.ys+info.ym-1; j>info.ys; j--,ic++) { // Left points
-    Y[j][i] = ub[ic];
-  }
-#endif
-
-  ierr = DMDAVecRestoreArray(da,x,&X);CHKERRQ(ierr);
-  ierr = DMDAVecRestoreArray(da,y,&Y);CHKERRQ(ierr);
-
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "PCView_FMM"
-static PetscErrorCode PCView_FMM(PC pc, PetscViewer viewer)
-{
-  PC_FMM *fmm = (PC_FMM*)pc->data;
-  PetscBool      isascii;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
-  if (!isascii) PetscFunctionReturn(0);
-  if (!fmm->da) PetscFunctionReturn(0);
-  ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
-  ierr = PetscViewerASCIIPrintf(viewer,"DAView for dof grid:\n");CHKERRQ(ierr);
-  ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
-  ierr = DMView(fmm->da,viewer);CHKERRQ(ierr);
-  ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
-  ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "PCReset_FMM"
-static PetscErrorCode PCReset_FMM(PC pc)
-{
-  PC_FMM *fmm = (PC_FMM*)pc->data;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  ierr = DMDestroy(&(fmm->da));CHKERRQ(ierr);
-  ierr = VecDestroy(&(fmm->u));CHKERRQ(ierr);
-  ierr = VecDestroy(&(fmm->f));CHKERRQ(ierr);
-  ierr = MatDestroy(&(fmm->G));CHKERRQ(ierr);
-  ierr = PetscFree3(fmm->xi,fmm->yi,fmm->ri);CHKERRQ(ierr);
-  ierr = PetscFree3(fmm->xb,fmm->yb,fmm->ub);CHKERRQ(ierr);
-  ierr = PetscFree3(fmm->xm,fmm->ym,fmm->um);CHKERRQ(ierr);
-  ierr = PetscFree3(fmm->dxe,fmm->dye,fmm->re);CHKERRQ(ierr);
-  ierr = PetscFree3(fmm->rhs,fmm->un,fmm->ui);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "PCDestroy_FMM"
-static PetscErrorCode PCDestroy_FMM(PC pc)
-{
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  ierr = PCReset_FMM(pc);CHKERRQ(ierr);
-  ierr = PetscFree(pc->data);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-EXTERN_C_BEGIN
-#undef  __FUNCT__
-#define __FUNCT__ "PCCreate_FMM"
-PetscErrorCode PCCreate_FMM(PC pc)
-{
-  PC_FMM *fmm = 0;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  ierr = PetscNewLog(pc,PC_FMM,&fmm);CHKERRQ(ierr);
-  pc->data = (void*)fmm;
-
-  pc->ops->setup               = PCSetUp_FMM;
-  pc->ops->reset               = PCReset_FMM;
-  pc->ops->destroy             = PCDestroy_FMM;
-  pc->ops->setfromoptions      = 0;/* PCSetFromOptions_FMM; */
-  pc->ops->view                = PCView_FMM;
-  pc->ops->apply               = PCApply_FMM;
-  pc->ops->applytranspose      = 0;/* PCApplyTranspose_FMM; */
-
-  PetscFunctionReturn(0);
-}
-EXTERN_C_END

fmm/fmm.h

-#ifndef __pcfmm_h
-#define __pcfmm_h
-
-typedef struct {
-  DM  da;
-  IGA iga;
-  Vec u,f;
-  Mat G;
-  PetscInt ni,nb,*ip1;
-  PetscReal *xi,*yi,*ri,*xb,*yb,*ub,*xm,*ym,*um,*dxe,*dye,*re,*rhs,*un,*ui,eps;
-} PC_FMM;
-
-#define PCFMM "fmm"
-
-PETSC_EXTERN PetscErrorCode PCFMMSetIGA(PC pc,IGA iga);
-
-extern "C" void FMM_Init(double eps2, int verbose);
-extern "C" void FMM(int ni, double * xi, double * yi, double * vi,
-  	            int nj, double * xj, double * yj, double * vj);
-#endif
+#include <assert.h>
+#include "petiga.h"
+#include "fmm2D.h"
+#include "nbody2D.h"
+#if PETSC_VERSION_LE(3,2,0)
+#include <private/pcimpl.h>
+#else
+#include <petsc-private/pcimpl.h>
+#endif
+
+#define DENSE 0
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCFMMSetIGA"
+PetscErrorCode PCFMMSetIGA(PC pc, IGA iga)
+{
+  PC_FMM *fmm = (PC_FMM*)pc->data;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,2);
+  if (iga) fmm->iga = iga;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "BEM_MatMult"
+static PetscErrorCode BEM_MatMult(Mat G, Vec x, Vec y)
+{
+  PC_FMM *fmm;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = MatShellGetContext(G,(void**)&fmm);CHKERRQ(ierr);
+  PetscInt i, nb = fmm->nb;
+  PetscReal *xl, *yl;
+  PetscReal *xm = fmm->xm;
+  PetscReal *ym = fmm->ym;
+  PetscReal *re = fmm->re;
+  PetscReal eps = fmm->eps;
+  VecGetSize(x,&i);
+  ierr = VecGetArray(x,&xl);CHKERRQ(ierr);
+  ierr = VecGetArray(y,&yl);CHKERRQ(ierr);
+  for (i=0; i<nb; i++) {
+    yl[i] = 0;
+  }
+
+  nbodyG(nb,xm,ym,yl,nb,xm,ym,xl,eps);
+
+  for (i=0; i<nb; i++) {
+    yl[i] += (1 - log(re[i] / 2) + log(eps)) * xl[i];
+  }
+
+  ierr = VecRestoreArray(x,&xl);CHKERRQ(ierr);
+  ierr = VecRestoreArray(y,&yl);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCSetUp_FMM"
+static PetscErrorCode PCSetUp_FMM(PC pc)
+{
+  PC_FMM *fmm = (PC_FMM*)pc->data;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  if (!fmm->iga) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCFMMSetIGA() must be called first");
+  ierr = IGACreateNodeDM(fmm->iga,1,&(fmm->da));CHKERRQ(ierr);
+  DMDALocalInfo info;
+  ierr = DMDAGetLocalInfo(fmm->da,&info);CHKERRQ(ierr);
+  fmm->nb = 0;
+  if(info.xs == 0){
+    fmm->nb += info.ym;
+    if(info.ys == 0)               fmm->nb -= 1;
+    if(info.my == info.ys+info.ym) fmm->nb -= 1;
+  }
+  if(info.mx == info.xs+info.xm){
+    fmm->nb += info.ym;
+    if(info.ys == 0)               fmm->nb -= 1;
+    if(info.my == info.ys+info.ym) fmm->nb -= 1;
+  }
+  if(info.ys == 0              ) fmm->nb += info.xm;
+  if(info.my == info.ys+info.ym) fmm->nb += info.xm;
+  fmm->ni = info.xm*info.ym - fmm->nb;
+
+  ierr = PetscMalloc(fmm->nb*sizeof(PetscInt),&(fmm->ip1));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->ni,PetscReal,&(fmm->xi),
+		      fmm->ni,PetscReal,&(fmm->yi),
+		      fmm->ni,PetscReal,&(fmm->ri));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->xb),
+		      fmm->nb,PetscReal,&(fmm->yb),
+		      fmm->nb,PetscReal,&(fmm->ub));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->xm),
+		      fmm->nb,PetscReal,&(fmm->ym),
+		      fmm->nb,PetscReal,&(fmm->um));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->dxe),
+		      fmm->nb,PetscReal,&(fmm->dye),
+		      fmm->nb,PetscReal,&(fmm->re));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->rhs),
+		      fmm->nb,PetscReal,&(fmm->un),
+		      fmm->ni,PetscReal,&(fmm->ui));CHKERRQ(ierr);
+  ierr = VecCreateMPI(PETSC_COMM_WORLD,fmm->nb,PETSC_DETERMINE,&(fmm->u));CHKERRQ(ierr);
+  ierr = VecDuplicate(fmm->u,&(fmm->f));CHKERRQ(ierr);
+  ierr = MatCreate(PETSC_COMM_WORLD,&fmm->G);CHKERRQ(ierr);
+  ierr = MatSetSizes(fmm->G,fmm->nb,fmm->nb,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
+#if DENSE
+  ierr = MatSetType(fmm->G,MATDENSE);CHKERRQ(ierr);
+#else
+  ierr = MatSetType(fmm->G,MATSHELL);CHKERRQ(ierr);
+  ierr = MatShellSetOperation(fmm->G,MATOP_MULT,(void(*)(void))BEM_MatMult);CHKERRQ(ierr);
+  ierr = MatShellSetContext(fmm->G,fmm);CHKERRQ(ierr);
+#endif
+  ierr = MatSetFromOptions(fmm->G);CHKERRQ(ierr);
+  ierr = MatSetUp(fmm->G);CHKERRQ(ierr);
+
+  PetscReal eps = .2 / (1<<((int)(log2(info.mx-1))));
+//  FMM_Init(eps*eps, 0);
+  fmm->eps = eps;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCApply_FMM"
+static PetscErrorCode PCApply_FMM(PC pc, Vec x, Vec y)
+{
+  PC_FMM *fmm = (PC_FMM*)pc->data;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+
+  DM da = fmm->da;
+
+  PetscScalar **X,**Y;
+  ierr = DMDAVecGetArray(da,x,&X);CHKERRQ(ierr);
+  ierr = DMDAVecGetArray(da,y,&Y);CHKERRQ(ierr);
+  DMDALocalInfo info;
+  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
+
+  Vec u = fmm->u;
+  Vec f = fmm->f;
+  Mat G = fmm->G;
+  IGA iga = fmm->iga;
+  PetscInt ni = fmm->ni;
+  PetscInt nb = fmm->nb;
+  PetscReal *xi  = fmm->xi;
+  PetscReal *yi  = fmm->yi;
+  PetscReal *ri  = fmm->ri;
+  //PetscReal *ub  = fmm->ub;
+  PetscReal *xm  = fmm->xm;
+  PetscReal *ym  = fmm->ym;
+  PetscReal *um  = fmm->um;
+  PetscReal *dxe = fmm->dxe;
+  PetscReal *dye = fmm->dye;
+  PetscReal *re  = fmm->re;
+  PetscReal *rhs = fmm->rhs;
+  PetscReal *un  = fmm->un;
+  PetscReal *ui  = fmm->ui;
+  PetscReal  eps = fmm->eps;
+
+  PetscInt i, j, ic = 0;
+  for (i=info.xs; i<info.xs+info.xm; i++) { // Loop over interior points
+    for (j=info.ys; j<info.ys+info.ym; j++) {
+      if (i != 0 && i != info.mx-1 && j != 0 && j != info.my-1) {
+	xi[ic] = (PetscReal)i / ( (PetscReal)(info.mx-1) );
+	yi[ic] = (PetscReal)j / ( (PetscReal)(info.my-1) );
+	ri[ic] = X[j][i] / 2 / M_PI;
+	ic++;
+      }
+    }
+  }
+  assert( ni == ic );
+
+  ic = 0;
+  j = 0;
+  PetscInt q;
+  PetscInt *start = iga->elem_start;
+  PetscInt *width = iga->elem_width;
+  PetscInt *sizes = iga->elem_sizes;
+  IGABasis *BD    = iga->basis;
+  PetscReal val_dirichlet = 0;
+  PetscReal h     =  BD[0]->detJ[0]*2;
+  for (j=0; j<2; j++){ // loop bottom and top
+    if (start[1]+j*width[1] == j*sizes[1]) { // domain boundary?
+      for (i=start[0]; i<start[0]+width[0]; i++) {
+	for (q=0; q<1; q++) {
+	  xm[ic]  =  (((PetscReal)i)+0.5)*h; //BD[0]->point[i*BD[0]->nqp+q];
+	  ym[ic]  =  j*1;
+	  um[ic]  = -val_dirichlet*0.5/M_PI;
+	  dxe[ic] =  BD[0]->detJ[i]*2 *(j == 0 ? 1 : -1);
+	  dye[ic] =  0;
+	  re[ic]  =  fabs(dxe[ic]);
+	  rhs[ic] =  um[ic]/re[ic]*M_PI;
+//	  printf("%d%d %f %f %f %f %f %f %f\n",start[0],start[1],xm[ic],ym[ic],um[ic],dxe[ic],dye[ic],re[ic],rhs[ic]);
+	  ic += 1;
+	}
+      }
+    }
+  }
+
+  for (j=0; j<2; j++){ // loop left and right
+    if (start[0]+j*width[0] == j*sizes[0]) { // domain boundary?
+      for (i=start[1]; i<start[1]+width[1]; i++) {
+	for (q=0; q<1; q++) {
+	  xm[ic]  =  j*1;
+	  ym[ic]  =  (((PetscReal)i)+0.5)*h; //BD[1]->point[i*BD[1]->nqp+q];
+	  um[ic]  = -val_dirichlet*0.5/M_PI;
+	  dxe[ic] =  0;
+	  dye[ic] =  BD[1]->detJ[i]*2 *(j == 0 ? -1 : 1);
+	  re[ic]  =  fabs(dye[ic]);
+	  rhs[ic] =  um[ic]/re[ic]*M_PI;
+//	  printf("%d%d %f %f %f %f %f %f %f\n",start[0],start[1],xm[ic],ym[ic],um[ic],dxe[ic],dye[ic],re[ic],rhs[ic]);
+	  ic += 1;
+	}
+      }
+    }
+  }
+
+  nb = ic; // this is ok, but we should improve memory allocation
+  assert( nb == ic );
+  
+  //FMM(nb,xm,ym,rhs,ni,xi,yi,ri);
+  nbodyG(nb,xm,ym,rhs,ni,xi,yi,ri,eps);
+  //nbodyGn(nb,xm,ym,rhs,nb,xm,ym,um,dxe,dye,re,eps);
+
+  for (i=0; i<nb; i++) {
+    VecSetValue(f,i,rhs[i],INSERT_VALUES);
+#if DENSE
+    for (j=0; j<i; j++) {
+      PetscReal dx = xm[i] - xm[j];
+      PetscReal dy = ym[i] - ym[j];
+      PetscReal r = sqrt(dx * dx + dy * dy + eps * eps);
+      PetscReal value = -log(r);
+      MatSetValue(G,i,j,value,INSERT_VALUES);
+      MatSetValue(G,j,i,value,INSERT_VALUES);
+    }
+    PetscReal value = 1-log(re[i]/2);
+    MatSetValue(G,i,i,value,INSERT_VALUES);
+#endif
+  }
+  ierr = VecAssemblyBegin(f);CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(f);CHKERRQ(ierr);
+#if DENSE
+  ierr = MatAssemblyBegin(G,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(G,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+#endif
+
+  KSP ksp;
+  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
+  ierr = KSPSetOperators(ksp,G,G,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
+  ierr = KSPSetOptionsPrefix(ksp,"pc_fmm_");CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+
+  ierr = KSPSolve(ksp,f,u);CHKERRQ(ierr);
+  ierr = VecGetArray(u,&un);CHKERRQ(ierr);
+
+  for (i=0; i<nb; i++) {
+    un[i] = -un[i];
+  }
+
+  for (i=0; i<ni; i++) {
+    ui[i] = 0;
+  }
+  //FMM(ni,xi,yi,ui,ni,xi,yi,ri);
+  //FMM(ni,xi,yi,ui,nb,xm,ym,un);
+  nbodyG(ni,xi,yi,ui,ni,xi,yi,ri,eps);
+  nbodyG(ni,xi,yi,ui,nb,xm,ym,un,eps);
+  //nbodyGn(ni,xi,yi,ui,nb,xm,ym,um,dxe,dye,re,eps);
+
+  ic = 0;
+  for (i=info.xs; i<info.xs+info.xm; i++) { // Loop over interior points
+    for (j=info.ys; j<info.ys+info.ym; j++) {
+      if (i != 0 && i != info.mx-1 && j != 0 && j != info.my-1) {
+        Y[j][i] = ui[ic];
+        ic++;
+      } else {
+        Y[j][i] = 0;
+      }
+    }
+  }
+
+#if 0
+  ic = 0;
+  j = 0;
+  for (i=info.xs; i<info.xs+info.xm-1; i++,ic++) { // Bottom points
+    Y[j][i] = ub[ic];
+  }
+  i = info.xs+info.xm-1;
+  for (j=info.ys; j<info.ys+info.ym-1; j++,ic++) { // Right points
+    Y[j][i] = ub[ic];
+  }
+  j = info.ys+info.ym-1;
+  for (i=info.xs+info.xm-1; i>info.xs; i--,ic++) { // Top points
+    Y[j][i] = ub[ic];
+  }
+  i = 0;
+  for (j=info.ys+info.ym-1; j>info.ys; j--,ic++) { // Left points
+    Y[j][i] = ub[ic];
+  }
+#endif
+
+  ierr = DMDAVecRestoreArray(da,x,&X);CHKERRQ(ierr);
+  ierr = DMDAVecRestoreArray(da,y,&Y);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCView_FMM"
+static PetscErrorCode PCView_FMM(PC pc, PetscViewer viewer)
+{
+  PC_FMM *fmm = (PC_FMM*)pc->data;
+  PetscBool      isascii;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
+  if (!isascii) PetscFunctionReturn(0);
+  if (!fmm->da) PetscFunctionReturn(0);
+  ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
+  ierr = PetscViewerASCIIPrintf(viewer,"DAView for dof grid:\n");CHKERRQ(ierr);
+  ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
+  ierr = DMView(fmm->da,viewer);CHKERRQ(ierr);
+  ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
+  ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCReset_FMM"
+static PetscErrorCode PCReset_FMM(PC pc)
+{
+  PC_FMM *fmm = (PC_FMM*)pc->data;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = DMDestroy(&(fmm->da));CHKERRQ(ierr);
+  ierr = VecDestroy(&(fmm->u));CHKERRQ(ierr);
+  ierr = VecDestroy(&(fmm->f));CHKERRQ(ierr);
+  ierr = MatDestroy(&(fmm->G));CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->xi,fmm->yi,fmm->ri);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->xb,fmm->yb,fmm->ub);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->xm,fmm->ym,fmm->um);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->dxe,fmm->dye,fmm->re);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->rhs,fmm->un,fmm->ui);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCDestroy_FMM"
+static PetscErrorCode PCDestroy_FMM(PC pc)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = PCReset_FMM(pc);CHKERRQ(ierr);
+  ierr = PetscFree(pc->data);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+EXTERN_C_BEGIN
+#undef  __FUNCT__
+#define __FUNCT__ "PCCreate_FMM"
+PetscErrorCode PCCreate_FMM(PC pc)
+{
+  PC_FMM *fmm = 0;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = PetscNewLog(pc,PC_FMM,&fmm);CHKERRQ(ierr);
+  pc->data = (void*)fmm;
+
+  pc->ops->setup               = PCSetUp_FMM;
+  pc->ops->reset               = PCReset_FMM;
+  pc->ops->destroy             = PCDestroy_FMM;
+  pc->ops->setfromoptions      = 0;/* PCSetFromOptions_FMM; */
+  pc->ops->view                = PCView_FMM;
+  pc->ops->apply               = PCApply_FMM;
+  pc->ops->applytranspose      = 0;/* PCApplyTranspose_FMM; */
+
+  PetscFunctionReturn(0);
+}
+EXTERN_C_END
+#ifndef __pcfmm2d_h
+#define __pcfmm2d_h
+
+typedef struct {
+  DM  da;
+  IGA iga;
+  Vec u,f;
+  Mat G;
+  PetscInt ni,nb,*ip1;
+  PetscReal *xi,*yi,*ri,*xb,*yb,*ub,*xm,*ym,*um,*dxe,*dye,*re,*rhs,*un,*ui,eps;
+} PC_FMM;
+
+#define PCFMM "fmm"
+
+PETSC_EXTERN PetscErrorCode PCFMMSetIGA(PC pc,IGA iga);
+
+#ifdef __cplusplus
+extern "C" void FMM_Init(double eps2, int verbose);
+extern "C" void FMM(int ni, double * xi, double * yi, double * vi,
+  	            int nj, double * xj, double * yj, double * vj);
+#endif
+#endif
+#include <assert.h>
+#include "petiga.h"
+#include "fmm3D.h"
+#include "nbody3D.h"
+#if PETSC_VERSION_LE(3,2,0)
+#include <private/pcimpl.h>
+#else
+#include <petsc-private/pcimpl.h>
+#endif
+
+#define DENSE 1
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCFMMSetIGA"
+PetscErrorCode PCFMMSetIGA(PC pc, IGA iga)
+{
+  PC_FMM *fmm = (PC_FMM*)pc->data;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,3);
+  if (iga) fmm->iga = iga;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "BEM_MatMult"
+static PetscErrorCode BEM_MatMult(Mat G, Vec x, Vec y)
+{
+  PC_FMM *fmm;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = MatShellGetContext(G,(void**)&fmm);CHKERRQ(ierr);
+  PetscInt i, nb = fmm->nb;
+  PetscReal *xl, *yl;
+  PetscReal *xm = fmm->xm;
+  PetscReal *ym = fmm->ym;
+  PetscReal *zm = fmm->zm;
+//  PetscReal *re = fmm->re;
+  PetscReal eps = fmm->eps;
+  VecGetSize(x,&i);
+  ierr = VecGetArray(x,&xl);CHKERRQ(ierr);
+  ierr = VecGetArray(y,&yl);CHKERRQ(ierr);
+  for (i=0; i<nb; i++) {
+    yl[i] = 0;
+  }
+//  nbodyG(nb,xm,ym,zm,yl,nb,xm,ym,zm,xl,eps);
+  for (i=0; i<nb; i++) {
+//    yl[i] += (1 - 1/(re[i] / 2) + (1/eps)) * xl[i];
+  }
+  ierr = VecRestoreArray(x,&xl);CHKERRQ(ierr);
+  ierr = VecRestoreArray(y,&yl);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCSetUp_FMM"
+static PetscErrorCode PCSetUp_FMM(PC pc)
+{
+  PC_FMM *fmm = (PC_FMM*)pc->data;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  if (!fmm->iga) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCFMMSetIGA() must be called first");
+  ierr = IGACreateNodeDM(fmm->iga,1,&(fmm->da));CHKERRQ(ierr);
+  DMDALocalInfo info;
+  ierr = DMDAGetLocalInfo(fmm->da,&info);CHKERRQ(ierr);
+  fmm->nb = 0;
+  if(info.xs == 0){
+    fmm->nb += info.ym*info.zm;
+    if(info.ys == 0)               fmm->nb -= info.zm;
+    if(info.zs == 0)               fmm->nb -= info.ym-2;
+    if(info.my == info.ys+info.ym) fmm->nb -= info.zm;
+    if(info.mz == info.zs+info.zm) fmm->nb -= info.ym-2;
+  }
+
+  if(info.mx == info.xs+info.xm){
+    fmm->nb += info.ym*info.zm;
+    if(info.ys == 0)               fmm->nb -= info.zm;
+    if(info.zs == 0)               fmm->nb -= info.ym-2;
+    if(info.my == info.ys+info.ym) fmm->nb -= info.zm;
+    if(info.mz == info.zs+info.zm) fmm->nb -= info.ym-2;
+  }
+
+  if(info.ys == 0){
+    fmm->nb += info.xm*info.zm;
+    if(info.xs == 0)               fmm->nb -= 2;
+    if(info.zs == 0)               fmm->nb -= info.xm-2;
+    if(info.mx == info.xs+info.xm) fmm->nb -= 2;
+    if(info.mz == info.zs+info.zm) fmm->nb -= info.xm-2;
+  }
+
+  if(info.my == info.ys+info.ym){
+    fmm->nb += info.xm*info.zm;
+    if(info.xs == 0)               fmm->nb -= 2;
+    if(info.zs == 0)               fmm->nb -= info.xm-2;
+    if(info.mx == info.xs+info.xm) fmm->nb -= 2;
+    if(info.mz == info.zs+info.zm) fmm->nb -= info.xm-2;
+  }
+
+  if(info.zs == 0              ) fmm->nb += info.xm*info.ym;
+  if(info.mz == info.zs+info.zm) fmm->nb += info.xm*info.ym;
+
+  fmm->ni = info.xm*info.ym*info.zm - fmm->nb;
+  fmm->nb = 12*(info.xm-1)*(info.ym-1);
+
+  ierr = PetscMalloc(fmm->ni*sizeof(PetscReal),&(fmm->ri));CHKERRQ(ierr);
+  ierr = PetscMalloc(fmm->nb*sizeof(PetscReal),&(fmm->jac));CHKERRQ(ierr);
+  ierr = PetscMalloc(fmm->nb*sizeof(PetscReal),&(fmm->BCV));CHKERRQ(ierr);
+  ierr = PetscMalloc(fmm->nb*sizeof(PetscInt),&(fmm->BCT));CHKERRQ(ierr);
+
+  ierr = PetscMalloc(16*sizeof(PetscReal),&(fmm->tg));CHKERRQ(ierr);
+  ierr = PetscMalloc(16*sizeof(PetscReal),&(fmm->vg));CHKERRQ(ierr);
+  ierr = PetscMalloc(fmm->nb*sizeof(PetscReal),&(fmm->phi));CHKERRQ(ierr);
+  ierr = PetscMalloc(fmm->nb*sizeof(PetscReal),&(fmm->dphi));CHKERRQ(ierr);
+
+  ierr = PetscMalloc3(fmm->ni,PetscReal,&(fmm->xi),
+		      fmm->ni,PetscReal,&(fmm->yi),
+		      fmm->ni,PetscReal,&(fmm->zi));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->xb0),
+		      fmm->nb,PetscReal,&(fmm->yb0),
+		      fmm->nb,PetscReal,&(fmm->zb0));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->xb1),
+                      fmm->nb,PetscReal,&(fmm->yb1),
+                      fmm->nb,PetscReal,&(fmm->zb1));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->xb2),
+                      fmm->nb,PetscReal,&(fmm->yb2),
+                      fmm->nb,PetscReal,&(fmm->zb2));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->xm),
+		      fmm->nb,PetscReal,&(fmm->ym),
+		      fmm->nb,PetscReal,&(fmm->zm));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->dxe),
+		      fmm->nb,PetscReal,&(fmm->dye),
+		      fmm->nb,PetscReal,&(fmm->dze));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->rhs),
+		      fmm->nb,PetscReal,&(fmm->un),
+		      fmm->ni,PetscReal,&(fmm->ui));CHKERRQ(ierr);
+  
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->PCv00),
+                      fmm->nb,PetscReal,&(fmm->PCv01),
+                      fmm->nb,PetscReal,&(fmm->PCv02));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->PCv10),
+                      fmm->nb,PetscReal,&(fmm->PCv11),
+                      fmm->nb,PetscReal,&(fmm->PCv12));CHKERRQ(ierr);
+  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->PCv20),
+                      fmm->nb,PetscReal,&(fmm->PCv21),
+                      fmm->nb,PetscReal,&(fmm->PCv22));CHKERRQ(ierr);
+
+  ierr = VecCreateMPI(PETSC_COMM_WORLD,fmm->nb,PETSC_DETERMINE,&(fmm->u));CHKERRQ(ierr);
+  ierr = VecDuplicate(fmm->u,&(fmm->f));CHKERRQ(ierr);
+  ierr = MatCreate(PETSC_COMM_WORLD,&fmm->G);CHKERRQ(ierr);
+  ierr = MatSetSizes(fmm->G,fmm->nb,fmm->nb,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
+
+#if DENSE
+  ierr = MatSetType(fmm->G,MATDENSE);CHKERRQ(ierr);
+#else
+  ierr = MatSetType(fmm->G,MATSHELL);CHKERRQ(ierr);
+  ierr = MatShellSetOperation(fmm->G,MATOP_MULT,(void(*)(void))BEM_MatMult);CHKERRQ(ierr);
+  ierr = MatShellSetContext(fmm->G,fmm);CHKERRQ(ierr);
+#endif
+  ierr = MatSetFromOptions(fmm->G);CHKERRQ(ierr);
+  ierr = MatSetUp(fmm->G);CHKERRQ(ierr);
+
+//TODO: check eps 
+  PetscReal eps = 0.0000001;//.2 / (1<<((int)(log2(info.mx-1))));
+//  FMM_Init(eps*eps, 0);
+  fmm->eps = eps;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCApply_FMM"
+static PetscErrorCode PCApply_FMM(PC pc, Vec x, Vec y)
+{
+  PC_FMM *fmm = (PC_FMM*)pc->data;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+
+  DM da = fmm->da;
+
+  PetscScalar ***X,***Y;
+  ierr = DMDAVecGetArray(da,x,&X);CHKERRQ(ierr);
+  ierr = DMDAVecGetArray(da,y,&Y);CHKERRQ(ierr);
+  DMDALocalInfo info;
+  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
+
+  Vec u = fmm->u;
+  Vec f = fmm->f;
+  Mat G = fmm->G;
+  IGA iga = fmm->iga;
+  PetscInt ni = fmm->ni;
+  PetscInt nb = fmm->nb;
+  PetscReal *xi  = fmm->xi;
+  PetscReal *yi  = fmm->yi;
+  PetscReal *zi  = fmm->zi;
+  PetscReal *ri  = fmm->ri;
+  PetscReal *xm  = fmm->xm;
+  PetscReal *ym  = fmm->ym;
+  PetscReal *zm  = fmm->zm;
+  PetscReal *dxe = fmm->dxe;
+  PetscReal *dye = fmm->dye;
+  PetscReal *dze = fmm->dze;
+  PetscReal *rhs = fmm->rhs;
+  PetscReal *un  = fmm->un;
+  PetscReal *ui  = fmm->ui;
+  PetscReal  eps = fmm->eps;
+
+  PetscInt i, j, k, ic = 0;
+  for (i=info.xs; i<info.xs+info.xm; i++) { // Loop over interior points
+    for (j=info.ys; j<info.ys+info.ym; j++) {
+      for (k=info.zs; k<info.zs+info.zm; k++) {
+        if (i != 0 && i != info.mx-1 && j != 0 && j != info.my-1 && k != 0 && k != info.mz-1) {
+	  xi[ic] = (PetscReal)i / ( (PetscReal)(info.mx-1) );
+	  yi[ic] = (PetscReal)j / ( (PetscReal)(info.my-1) );
+          zi[ic] = (PetscReal)k / ( (PetscReal)(info.mz-1) );
+	  ri[ic] = X[k][j][i];
+	  ic++;
+        }
+      }
+    }
+  }
+
+  assert( ni == ic );
+
+  ic = 0;
+  PetscInt ie, m;
+  PetscInt *start = iga->elem_start;
+  PetscInt *width = iga->elem_width;
+  PetscInt *sizes = iga->elem_sizes;
+  PetscReal *BCV = fmm->BCV;
+  PetscInt *BCT = fmm->BCT;
+
+  PetscReal x21, x31, y21, y31, z21, z31;
+  PetscReal xk, yk, zk, D1, D2, del, ut;
+  PetscInt  N0 = (info.mx-1);
+  PetscReal aa, bb, cc, dd, s1, s2, s3 = 1.0/sqrt(3.0), dl;
+  PetscReal *jac = fmm->jac;
+  PetscReal *tg = fmm->tg, *vg = fmm->vg;
+  PetscReal *phi = fmm->phi, *dphi = fmm->dphi;
+  PetscReal *xb0 = fmm->xb0, *xb1 = fmm->xb1, *xb2 = fmm->xb2;
+  PetscReal *yb0 = fmm->yb0, *yb1 = fmm->yb1, *yb2 = fmm->yb2;
+  PetscReal *zb0 = fmm->zb0, *zb1 = fmm->zb1, *zb2 = fmm->zb2;
+
+  PetscReal *PCv00 = fmm->PCv00, *PCv01 = fmm->PCv01, *PCv02 = fmm->PCv02;
+  PetscReal *PCv10 = fmm->PCv10, *PCv11 = fmm->PCv11, *PCv12 = fmm->PCv12;
+  PetscReal *PCv20 = fmm->PCv20, *PCv21 = fmm->PCv21, *PCv22 = fmm->PCv22;
+
+  dl = 1.0/(PetscReal)(2*N0);
+  ie = -2;
+
+  for(j=start[1]; j<start[1]+width[1]; j++){
+   if (start[1]+j*width[1] == j*sizes[1]) {
+    for (i=start[0]; i<start[0]+width[0]; i++) {
+      ie += 2;
+      s1  = (PetscReal)(2*j+1)*dl;
+      s2  = (PetscReal)(2*i+1)*dl;
+      xb0[ie] = 0.0; 
+        yb0[ie] = s1-dl;
+        zb0[ie] = s2-dl;
+      xb1[ie] = 0.0;
+        yb1[ie] = s1+dl;
+        zb1[ie] = s2+dl;
+      xb2[ie] = 0.0;
+        yb2[ie] = s1+dl;
+        zb2[ie] = s2-dl;
+      xb0[ie+1] = 0.0;
+        yb0[ie+1] = s1-dl;
+        zb0[ie+1] = s2-dl;
+      xb1[ie+1] = 0.0;
+        yb1[ie+1] = s1-dl;
+        zb1[ie+1] = s2+dl;
+      xb2[ie+1] = 0.0;
+        yb2[ie+1] = s1+dl;
+        zb2[ie+1] = s2+dl;
+      xb0[2*N0*N0+ie] = 1.0;
+        yb0[2*N0*N0+ie] = s1-dl;
+        zb0[2*N0*N0+ie] = s2-dl;
+      xb1[2*N0*N0+ie] = 1.0;
+        yb1[2*N0*N0+ie] = s1+dl;
+        zb1[2*N0*N0+ie] = s2-dl;
+      xb2[2*N0*N0+ie] = 1.0;
+        yb2[2*N0*N0+ie] = s1+dl;
+        zb2[2*N0*N0+ie] = s2+dl;
+      xb0[2*N0*N0+ie+1] = 1.0;
+        yb0[2*N0*N0+ie+1] = s1-dl;
+        zb0[2*N0*N0+ie+1] = s2-dl;
+      xb1[2*N0*N0+ie+1] = 1.0;
+        yb1[2*N0*N0+ie+1] = s1+dl;
+        zb1[2*N0*N0+ie+1] = s2+dl;
+      xb2[2*N0*N0+ie+1] = 1.0;
+        yb2[2*N0*N0+ie+1] = s1-dl;
+        zb2[2*N0*N0+ie+1] = s2+dl;
+      xb0[4*N0*N0+ie] = s1+dl;
+        yb0[4*N0*N0+ie] = 0.0;
+	zb0[4*N0*N0+ie] = s2-dl;
+	xb1[4*N0*N0+ie] = s1-dl;
+        yb1[4*N0*N0+ie] = 0.0;
+        zb1[4*N0*N0+ie] = s2+dl;
+      xb2[4*N0*N0+ie] = s1-dl;
+        yb2[4*N0*N0+ie] = 0.0;
+        zb2[4*N0*N0+ie] = s2-dl;
+      xb0[4*N0*N0+ie+1] = s1+dl;
+        yb0[4*N0*N0+ie+1] = 0.0;
+        zb0[4*N0*N0+ie+1] = s2-dl;
+      xb1[4*N0*N0+ie+1] = s1+dl;
+        yb1[4*N0*N0+ie+1] = 0.0;
+        zb1[4*N0*N0+ie+1] = s2+dl;
+      xb2[4*N0*N0+ie+1] = s1-dl;
+        yb2[4*N0*N0+ie+1] = 0.0;
+        zb2[4*N0*N0+ie+1] = s2+dl;
+      xb0[6*N0*N0+ie] = s1-dl;
+        yb0[6*N0*N0+ie] = 1.0;
+        zb0[6*N0*N0+ie] = s2-dl;
+      xb1[6*N0*N0+ie] = s1-dl;
+        yb1[6*N0*N0+ie] = 1.0;
+        zb1[6*N0*N0+ie] = s2+dl;
+      xb2[6*N0*N0+ie] = s1+dl;
+        yb2[6*N0*N0+ie] = 1.0;
+        zb2[6*N0*N0+ie] = s2-dl;
+      xb0[6*N0*N0+ie+1] = s1-dl;
+        yb0[6*N0*N0+ie+1] = 1.0;
+        zb0[6*N0*N0+ie+1] = s2+dl;
+      xb1[6*N0*N0+ie+1] = s1+dl;
+        yb1[6*N0*N0+ie+1] = 1.0;
+        zb1[6*N0*N0+ie+1] = s2+dl;
+      xb2[6*N0*N0+ie+1] = s1+dl;
+        yb2[6*N0*N0+ie+1] = 1.0;
+        zb2[6*N0*N0+ie+1] = s2-dl;
+      xb0[8*N0*N0+ie] = s1-dl;
+        yb0[8*N0*N0+ie] = s2-dl;
+        zb0[8*N0*N0+ie] = 0.0;
+      xb1[8*N0*N0+ie] = s1+dl;
+        yb1[8*N0*N0+ie] = s2+dl;
+        zb1[8*N0*N0+ie] = 0.0;
+      xb2[8*N0*N0+ie] = s1+dl;
+        yb2[8*N0*N0+ie] = s2-dl;
+        zb2[8*N0*N0+ie] = 0.0 ;
+      xb0[8*N0*N0+ie+1] = s1-dl;
+        yb0[8*N0*N0+ie+1] = s2-dl;
+        zb0[8*N0*N0+ie+1] = 0.0;
+      xb1[8*N0*N0+ie+1] = s1-dl;
+        yb1[8*N0*N0+ie+1] = s2+dl;
+        zb1[8*N0*N0+ie+1] = 0.0;
+      xb2[8*N0*N0+ie+1] = s1+dl;
+        yb2[8*N0*N0+ie+1] = s2+dl;
+        zb2[8*N0*N0+ie+1] = 0.0;
+      xb0[10*N0*N0+ie] = s1-dl;
+        yb0[10*N0*N0+ie] = s2-dl;
+        zb0[10*N0*N0+ie] = 1.0;
+      xb1[10*N0*N0+ie] = s1+dl;
+        yb1[10*N0*N0+ie] = s2-dl;
+        zb1[10*N0*N0+ie] = 1.0;
+      xb2[10*N0*N0+ie] = s1+dl;
+        yb2[10*N0*N0+ie] = s2+dl;
+        zb2[10*N0*N0+ie] = 1.0;
+      xb0[10*N0*N0+ie+1] = s1-dl;
+        yb0[10*N0*N0+ie+1] = s2-dl;
+        zb0[10*N0*N0+ie+1] = 1.0;
+      xb1[10*N0*N0+ie+1] = s1+dl;
+        yb1[10*N0*N0+ie+1] = s2+dl;
+        zb1[10*N0*N0+ie+1] = 1.0;
+      xb2[10*N0*N0+ie+1] = s1-dl;
+        yb2[10*N0*N0+ie+1] = s2+dl;
+        zb2[10*N0*N0+ie+1] = 1.0;
+    }
+   }
+  }
+
+  tg[0]  = 0.25*(1.0+s3);
+  vg[0]  = tg[0];
+  tg[1]  = tg[0];
+  vg[1]  = 0.25*(1.0-s3);
+  tg[2]  = vg[1];
+  vg[2]  = tg[0];
+  tg[3]  = vg[1];
+  vg[3]  = vg[1];
+  tg[4]  = 0.25*(3.0+s3);
+  vg[4]  = tg[0];
+  tg[5]  = tg[4];
+  vg[5]  = vg[1];
+  tg[6]  = 0.25*(3.0-s3);
+  vg[6]  = tg[0];
+  tg[7]  = tg[6];
+  vg[7]  = vg[1];
+  tg[8]  = tg[4];
+  vg[8]  = tg[4];
+  tg[9]  = tg[4];
+  vg[9]  = tg[6];
+  tg[10] = tg[6];
+  vg[10] = tg[4];
+  tg[11] = tg[6];
+  vg[11] = tg[6];
+  tg[12] = tg[0];
+  vg[12] = tg[4];
+  tg[13] = tg[0];
+  vg[13] = tg[6];
+  tg[14] = vg[1];
+  vg[14] = tg[4];
+  tg[15] = vg[1];
+  vg[15] = tg[6];
+
+  for(i = 0; i < 12*width[0]*width[1]; i++){
+    x21 = xb1[i]-xb0[i];
+    x31 = xb2[i]-xb0[i];
+    y21 = yb1[i]-yb0[i];
+    y31 = yb2[i]-yb0[i];
+    z21 = zb1[i]-zb0[i];
+    z31 = zb2[i]-zb0[i];
+    dxe[i] = y21*z31-z21*y31;
+    dye[i] = z21*x31-x21*z31;
+    dze[i] = x21*y31-y21*x31;
+  
+    dd = sqrt(dxe[i]*dxe[i]+dye[i]*dye[i]+dze[i]*dze[i]);
+    dxe[i] = dxe[i]/dd;
+    dye[i] = dye[i]/dd;
+    dze[i] = dze[i]/dd;
+
+    aa = sqrt(x21*x21+y21*y21+z21*z21);
+    bb = sqrt(((xb1[i]-xb2[i])*(xb1[i]-xb2[i]))+((yb1[i]-yb2[i])*(yb1[i]-yb2[i]))+((zb1[i]-zb2[i])*(zb1[i]-zb2[i])));
+    cc = sqrt(x31*x31+y31*y31+z31*z31);
+    dd = 0.5*(aa+bb+cc);
+    jac[i] = 2.0*sqrt(dd*(dd-aa)*(dd-bb)*(dd-cc));
+
+    if (fabs(dze[i]) >= s3){
+      PCv00[i] = x21;
+      PCv01[i] = x31;
+      PCv02[i] = xb0[i];
+      PCv10[i] = y21;
+      PCv11[i] = y31;
+      PCv12[i] = yb0[i];
+      PCv20[i] = -(dxe[i]*PCv00[i]+dye[i]*PCv10[i])/dze[i];
+      PCv21[i] = -(dxe[i]*PCv01[i]+dye[i]*PCv11[i])/dze[i];
+      PCv22[i] = zb0[i];
+    }
+    else if (fabs(dye[i]) >= s3){
+      PCv00[i] = x21;
+      PCv01[i] = x31;
+      PCv02[i] = xb0[i];
+      PCv20[i] = z21;
+      PCv21[i] = z31;
+      PCv22[i] = zb0[i];
+      PCv10[i] =- (dxe[i]*PCv00[i]+dze[i]*PCv20[i])/dye[i];
+      PCv11[i] =- (dxe[i]*PCv01[i]+dze[i]*PCv21[i])/dye[i];
+      PCv12[i] = yb0[i];
+     }
+    else{
+      PCv10[i] = y21;
+      PCv11[i] = y31;
+      PCv12[i] = yb0[i];
+      PCv20[i] = z21;
+      PCv21[i] = z31;
+      PCv22[i] = zb0[i];
+      PCv00[i] =- (dze[i]*PCv20[i]+dye[i]*PCv10[i])/dxe[i];
+      PCv01[i] =- (dze[i]*PCv21[i]+dye[i]*PCv11[i])/dxe[i];
+      PCv02[i] = xb0[i];
+    }
+    xm[i] = PCv00[i]*0.25+PCv01[i]*0.5+PCv02[i];
+    ym[i] = PCv10[i]*0.25+PCv11[i]*0.5+PCv12[i];
+    zm[i] = PCv20[i]*0.25+PCv21[i]*0.5+PCv22[i];
+    ic += 1;
+  }
+
+  for(i = 0; i < nb; i++){
+    if (i < (2*N0*N0)){
+      BCT[i] = 0;
+      BCV[i] = nbodyG(xm[i],ym[i],zm[i],0.0,0.0,1.5,eps);
+    }
+    else{
+      BCT[i] = 1;
+      BCV[i] = nbodyGn(xm[i],ym[i],zm[i],0.0,0.0,1.5,dxe[i],dye[i],dze[i],eps);
+    }
+  }
+
+  nb = ic; // this is ok, but we should improve memory allocation
+  assert( nb == ic );
+
+  for (m = 0; m < nb; m++){
+    rhs[m] = 0.0;
+    for (k = 0; k < nb; k++){
+      D1 = 0.0;
+      D2 = 0.0;
+      for (i = 0; i < 16; i++){
+        ut  = tg[i]*(1.0-vg[i]);
+        xk  = PCv00[k]*ut+PCv01[k]*vg[i]+PCv02[k];
+        yk  = PCv10[k]*ut+PCv11[k]*vg[i]+PCv12[k];
+        zk  = PCv20[k]*ut+PCv21[k]*vg[i]+PCv22[k];
+        D1 += nbodyG(xk,yk,zk,xm[m],ym[m],zm[m],eps)*(1.0-vg[i]);
+        D2 += nbodyGn(xk,yk,zk,xm[m],ym[m],zm[m],dxe[k],dye[k],dze[k],eps)*(1.0-vg[i]);
+      }
+      D1 = jac[k]*D1/16.0;
+      D2 = jac[k]*D2/16.0;
+
+      if (k == m)
+        del = 1.0;
+      else
+        del = 0.0;
+
+      if (BCT[k] == 0){
+        rhs[m] += BCV[k]*(-D2+0.5*del);
+        #if DENSE
+          PetscReal value = -D1;
+          MatSetValue(G,m,k,value,INSERT_VALUES);
+        #endif
+      }
+      else{
+        rhs[m] += BCV[k]*D1;
+        #if DENSE
+          PetscReal value = D2-0.5*del;
+          MatSetValue(G,m,k,value,INSERT_VALUES);
+        #endif
+      }
+    }
+    VecSetValue(f,m,rhs[m],INSERT_VALUES);
+  }
+
+  ierr = VecAssemblyBegin(f);CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(f);CHKERRQ(ierr);
+  #if DENSE
+    ierr = MatAssemblyBegin(G,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+    ierr = MatAssemblyEnd(G,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  #endif
+   
+  KSP ksp;
+  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
+  ierr = KSPSetOperators(ksp,G,G,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
+  ierr = KSPSetOptionsPrefix(ksp,"pc_fmm_");CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+
+  ierr = KSPSolve(ksp,f,u);CHKERRQ(ierr);
+  ierr = VecGetArray(u,&un);CHKERRQ(ierr);
+
+  for (m = 0; m < nb; m++){
+    if (BCT[m]==0){
+      phi[m]  = BCV[m];
+      dphi[m] = un[m];
+    }
+    else{
+      phi[m]  = un[m];
+      dphi[m] = BCV[m];
+    }
+  }
+
+  for (j = 0; j < ni; j++){
+    for (k = 0; k < nb; k++){
+      D1 = 0.0;
+      D2 = 0.0;
+      for (i = 0; i < 16; i++){
+        ut  = tg[i]*(1.0-vg[i]);
+        xk  = PCv00[k]*ut+PCv01[k]*vg[i]+PCv02[k];
+        yk  = PCv10[k]*ut+PCv11[k]*vg[i]+PCv12[k];
+        zk  = PCv20[k]*ut+PCv21[k]*vg[i]+PCv22[k];
+        D1 += nbodyG(xk,yk,zk,xi[j],yi[j],zi[j],eps)*(1.0-vg[i]);
+        D2 += nbodyGn(xk,yk,zk,xi[j],yi[j],zi[j],dxe[k],dye[k],dze[k],eps)*(1.0-vg[i]);
+      }
+      D1 = jac[k]*D1/16.0;
+      D2 = jac[k]*D2/16.0;
+      ui[j] += phi[k]*D2-dphi[k]*D1;
+    }
+  }
+
+  ic = 0;
+  for (i=info.xs; i<info.xs+info.xm; i++) { // Loop over interior points
+    for (j=info.ys; j<info.ys+info.ym; j++) {
+      for (k=info.zs; k<info.zs+info.zm; k++) {
+        if (i != 0 && i != info.mx-1 && j != 0 && j != info.my-1 && k != 0 && k != info.mz-1) {
+          Y[k][j][i] = ui[ic];
+          ic++;
+        } else {
+          Y[k][j][i] = 0;
+        }
+      }
+    }
+  }
+  
+  ierr = DMDAVecRestoreArray(da,x,&X);CHKERRQ(ierr);
+  ierr = DMDAVecRestoreArray(da,y,&Y);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCView_FMM"
+static PetscErrorCode PCView_FMM(PC pc, PetscViewer viewer)
+{
+  PC_FMM *fmm = (PC_FMM*)pc->data;
+  PetscBool      isascii;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
+  if (!isascii) PetscFunctionReturn(0);
+  if (!fmm->da) PetscFunctionReturn(0);
+  ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
+  ierr = PetscViewerASCIIPrintf(viewer,"DAView for dof grid:\n");CHKERRQ(ierr);
+  ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
+  ierr = DMView(fmm->da,viewer);CHKERRQ(ierr);
+  ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
+  ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCReset_FMM"
+static PetscErrorCode PCReset_FMM(PC pc)
+{
+  PC_FMM *fmm = (PC_FMM*)pc->data;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = DMDestroy(&(fmm->da));CHKERRQ(ierr);
+  ierr = VecDestroy(&(fmm->u));CHKERRQ(ierr);
+  ierr = VecDestroy(&(fmm->f));CHKERRQ(ierr);
+  ierr = MatDestroy(&(fmm->G));CHKERRQ(ierr);
+  ierr = PetscFree(fmm->tg);CHKERRQ(ierr);
+  ierr = PetscFree(fmm->vg);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->xi,fmm->yi,fmm->zi);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->xb0,fmm->yb0,fmm->zb0);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->xb1,fmm->yb1,fmm->zb1);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->xb2,fmm->yb2,fmm->zb2);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->xm,fmm->ym,fmm->zm);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->dxe,fmm->dye,fmm->dze);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->rhs,fmm->un,fmm->ui);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->PCv00,fmm->PCv01,fmm->PCv02);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->PCv10,fmm->PCv11,fmm->PCv12);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->PCv20,fmm->PCv21,fmm->PCv22);CHKERRQ(ierr);
+  ierr = PetscFree3(fmm->ri,fmm->jac,fmm->BCV);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "PCDestroy_FMM"
+static PetscErrorCode PCDestroy_FMM(PC pc)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = PCReset_FMM(pc);CHKERRQ(ierr);
+  ierr = PetscFree(pc->data);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+EXTERN_C_BEGIN
+#undef  __FUNCT__
+#define __FUNCT__ "PCCreate_FMM"
+PetscErrorCode PCCreate_FMM(PC pc)
+{
+  PC_FMM *fmm = 0;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = PetscNewLog(pc,PC_FMM,&fmm);CHKERRQ(ierr);
+  pc->data = (void*)fmm;
+  pc->ops->setup               = PCSetUp_FMM;
+  pc->ops->reset               = PCReset_FMM;
+  pc->ops->destroy             = PCDestroy_FMM;
+  pc->ops->setfromoptions      = 0;// PCSetFromOptions_FMM;
+  pc->ops->view                = PCView_FMM;
+  pc->ops->apply               = PCApply_FMM;
+  pc->ops->applytranspose      = 0;// PCApplyTranspose_FMM;
+
+  PetscFunctionReturn(0);
+}
+EXTERN_C_END
+#ifndef __pcfmm3d_h
+#define __pcfmm3d_h
+
+typedef struct {
+  DM  da;
+  IGA iga;
+  Vec u,f;
+  Mat G;
+  PetscInt ni,nb,*BCT;
+  PetscReal *xi,*yi,*zi,*ri,*ub,*xm,*ym,*zm,*um,*dxe,*dye,*dze,*re,*rhs,*un,*ui,*jac,eps;
+  PetscReal *BCV, *xb0,*yb0,*zb0,*xb1,*yb1,*zb1,*xb2,*yb2,*zb2, *tg, *vg, *phi, *dphi;
+  PetscReal *PCv00, *PCv01, *PCv02, *PCv10, *PCv11, *PCv12, *PCv20, *PCv21, *PCv22;
+} PC_FMM;
+
+#define PCFMM "fmm"
+
+PETSC_EXTERN PetscErrorCode PCFMMSetIGA(PC pc,IGA iga);
+
+#ifdef __cplusplus
+extern "C" void FMM_Init(double eps2, int verbose);
+extern "C" void FMM(int ni, double * xi, double * yi, double * zi, double * vi,
+  	            int nj, double * xj, double * yj, double * zj, double * vj);
+#endif
+#endif
-TARGETS = Poisson2D
+2d: Poisson2D.o fmm2D.o
+	${CLINKER} $^ ${PETIGA_LIB} -L../../exafmm2d/wrappers -lfmm -ltbb -lm
+	mpirun -np 1 ./a.out -ksp_monitor -ksp_norm_type PRECONDITIONED -ksp_max_it 100 -pc_type hypre 
 
-all: ${TARGETS}
+3d: Poisson3D.o fmm3D.o
+	${CLINKER} $^ ${PETIGA_LIB} -L../../exafmm2d/wrappers -lfmm -ltbb -lm
+	mpirun -np 1 ./a.out -ksp_monitor -ksp_norm_type PRECONDITIONED -ksp_max_it 15 -pc_type hypre
 
-Poisson2D: Poisson2D.o fmm.o
-	${CLINKER} $^ -o $@ ${PETIGA_LIB} -L../../exafmm2d/wrappers -lfmm -ltbb
-#	${RM} -f $<
-#	mpirun -np 1 ./Poisson2D -ksp_monitor -ksp_norm_type UNPRECONDITIONED -ksp_max_it 100 -pc_type fmm -log_summary
-	mpirun -np 1 ./Poisson2D -ksp_monitor -ksp_norm_type UNPRECONDITIONED -draw_pause -1 -iga_elements 64
 clean::
-	-@${RM} ${TARGETS}
+	-@${RM} *.o *.out
 
 include ${PETIGA_DIR}/conf/petigavariables
 include ${PETIGA_DIR}/conf/petigarules
-

fmm/nbody.h

-void nbodyG(int ni, PetscReal *xi, PetscReal *yi, PetscReal *vi,
-            int nj, PetscReal *xj, PetscReal *yj, PetscReal *vj, PetscReal eps) {
-  PetscInt i, j;
-  PetscReal eps2 = eps * eps;
-  for (i=0; i<ni; i++) {
-    PetscReal v = 0;
-    for (j=0; j<nj; j++) {
-      PetscReal dx = xi[i] - xj[j];
-      PetscReal dy = yi[i] - yj[j];
-      PetscReal r = sqrt(dx * dx + dy * dy + eps2);
-      v += vj[j] * log(r);
-    }
-    vi[i] -= v;
-  }
-}
-
-void nbodyGn(int ni, PetscReal *xi, PetscReal *yi, PetscReal *vi,
-             int nj, PetscReal *xj, PetscReal *yj, PetscReal *vj,
-	     PetscReal *dxe, PetscReal *dye, PetscReal *re, PetscReal eps) {
-  PetscInt i, j;
-  PetscReal eps2 = eps * eps;
-  for (i=0; i<ni; i++) {
-    PetscReal v = 0;
-    for (j=0; j<nj; j++) {
-      PetscReal dx = xi[i] - xj[j];
-      PetscReal dy = yi[i] - yj[j];
-      PetscReal rdn = dx * dye[j] - dy * dxe[j];
-      PetscReal r2 = dx * dx + dy * dy + eps2;
-      v += vj[j] * rdn / re[j] / r2;
-    }
-    vi[i] -= v;
-  }
-}
+void nbodyG(int ni, PetscReal *xi, PetscReal *yi, PetscReal *vi,
+            int nj, PetscReal *xj, PetscReal *yj, PetscReal *vj, PetscReal eps) {
+  PetscInt i, j;
+  PetscReal eps2 = eps * eps;
+  for (i=0; i<ni; i++) {
+    PetscReal v = 0;
+    for (j=0; j<nj; j++) {
+      PetscReal dx = xi[i] - xj[j];
+      PetscReal dy = yi[i] - yj[j];
+      PetscReal r = sqrt(dx * dx + dy * dy + eps2);
+      v += vj[j] * log(r);
+    }
+    vi[i] -= v;
+  }
+}
+
+void nbodyGn(int ni, PetscReal *xi, PetscReal *yi, PetscReal *vi,
+             int nj, PetscReal *xj, PetscReal *yj, PetscReal *vj,
+	     PetscReal *dxe, PetscReal *dye, PetscReal *re, PetscReal eps) {
+  PetscInt i, j;
+  PetscReal eps2 = eps * eps;
+  for (i=0; i<ni; i++) {
+    PetscReal v = 0;
+    for (j=0; j<nj; j++) {
+      PetscReal dx = xi[i] - xj[j];
+      PetscReal dy = yi[i] - yj[j];
+      PetscReal rdn = dx * dye[j] - dy * dxe[j];
+      PetscReal r2 = dx * dx + dy * dy + eps2;
+      v += vj[j] * rdn / re[j] / r2;
+    }
+    vi[i] -= v;
+  }
+}
+double nbodyG(PetscReal xi, PetscReal yi, PetscReal zi,
+            PetscReal xj, PetscReal yj, PetscReal zj, PetscReal eps) {
+  PetscReal eps2 = eps * eps;
+  PetscReal dx = xi - xj;
+  PetscReal dy = yi - yj;
+  PetscReal dz = zi - zj;
+  PetscReal r  = sqrt(dx * dx + dy * dy + dz * dz + eps2);
+  return (-0.25/r/M_PI);
+}
+
+double nbodyGn(PetscReal xi, PetscReal yi, PetscReal zi,
+             PetscReal xj, PetscReal yj, PetscReal zj,
+             PetscReal dxe, PetscReal dye, PetscReal dze, PetscReal eps) {
+  PetscReal eps2 = eps * eps;
+  PetscReal dx = xi - xj;
+  PetscReal dy = yi - yj;
+  PetscReal dz = zi - zj;
+  PetscReal rdn = dx * dxe + dy * dye + dz * dze;
+  PetscReal r2 = dx * dx + dy * dy + dz * dz + eps2;
+   
+  return 0.25 * rdn / pow(r2,1.5) / M_PI;
+}
+/*void nbodyG(int ni, PetscReal *xi, PetscReal *yi, PetscReal *zi, PetscReal *vi,
+            int nj, PetscReal *xj, PetscReal *yj, PetscReal *zj, PetscReal *vj, PetscReal eps) {
+  PetscInt i, j;
+  PetscReal eps2 = eps * eps;
+  for (i=0; i<ni; i++) {
+    PetscReal v = 0;
+    for (j=0; j<nj; j++) {
+      PetscReal dx = xi[i] - xj[j];
+      PetscReal dy = yi[i] - yj[j];
+      PetscReal dz = zi[i] - zj[j];
+      PetscReal r = sqrt(dx * dx + dy * dy + dz * dz + eps2);
+      v -= vj[j] * (0.25/r/M_PI);
+    }
+    vi[i] += v;
+  }
+}
+
+void nbodyGn(int ni, PetscReal *xi, PetscReal *yi, PetscReal *zi, PetscReal *vi,
+             int nj, PetscReal *xj, PetscReal *yj, PetscReal *zj, PetscReal *vj,
+	     PetscReal *dxe, PetscReal *dye, PetscReal *dze, PetscReal *re, PetscReal eps) {
+  PetscInt i, j;
+  PetscReal eps2 = eps * eps;
+  for (i=0; i<ni; i++) {
+    PetscReal v = 0;
+    for (j=0; j<nj; j++) {
+      PetscReal dx = xi[i] - xj[j];
+      PetscReal dy = yi[i] - yj[j];
+      PetscReal dz = zi[i] - zj[j];
+      PetscReal rdn = dx * dxe[j] + dy * dye[j] + dz * dze[j];
+      PetscReal r2 = dx * dx + dy * dy + dz * dz + eps2;
+      v += vj[j] * 0.25 * rdn / pow(re[j],1.5) / M_PI;
+    }
+    vi[i] += v;
+  }
+}*/
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.