Source

PetIGA / src / petigapc.c

#include "petiga.h"
/* ---------------------------------------------------------------- */

PETSC_STATIC_INLINE
PetscInt Index(const PetscInt shape[],
               PetscInt i,PetscInt j,PetscInt k)
{
  return i + j * shape[0] + k * shape[0] * shape[1];
}

PETSC_STATIC_INLINE
PetscInt Color(const PetscInt shape[3],
               const PetscInt start[3],
               const PetscInt width[3],
               PetscInt i,PetscInt j,PetscInt k)
{
  PetscInt r=0,g=0,b=0;
  PetscInt L[3],R[3],C[2];

  if (i<0||i>=shape[0]) return PETSC_MIN_INT;
  if (j<0||j>=shape[1]) return PETSC_MIN_INT;
  if (k<0||k>=shape[2]) return PETSC_MIN_INT;

  L[0] = start[0]; R[0] = start[0] + width[0] - 1;
  L[1] = start[1]; R[1] = start[1] + width[1] - 1;
  L[2] = start[2]; R[2] = start[2] + width[2] - 1;
  if (i<L[0]) r = i - L[0]; if (i>R[0]) r = i - R[0];
  if (j<L[1]) g = j - L[1]; if (j>R[1]) g = j - R[1];
  if (k<L[2]) b = k - L[2]; if (k>R[2]) b = k - R[2];

  C[0] = shape[0] - width[0] + 1;
  C[1] = shape[1] - width[1] + 1;
  return Index(C,r,g,b);
}

static const
PetscInt STENCIL[7][3] = {{ 0,  0, -1},
                          { 0, -1,  0},
                          {-1,  0,  0},
                          { 0,  0,  0},
                          { 0,  0, +1},
                          { 0, +1,  0},
                          {+1,  0,  0}};

PETSC_STATIC_INLINE
#undef  __FUNCT__
#define __FUNCT__ "ComputeBDDCGraph"
PetscErrorCode ComputeBDDCGraph(PetscInt bs,const PetscInt shape[3],
                                const PetscInt start[3],const PetscInt width[3],
                                PetscInt *_nvtx,PetscInt *_xadj[],PetscInt *_adjy[])
{
  PetscInt       c,i,j,k,s,v,pos;
  PetscInt       nvtx=0,*xadj;
  PetscInt       nadj=0,*adjy;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidIntPointer(_nvtx,5);
  PetscValidPointer(_xadj,6);
  PetscValidPointer(_adjy,7);

  /* Compute the number of vertices and adjacencies */
  nvtx = shape[0]*shape[1]*shape[2];
  for (k=0; k<shape[2]; k++)
    for (j=0; j<shape[1]; j++)
      for (i=0; i<shape[0]; i++) {
        PetscInt color = Color(shape,start,width,i,j,k);
        for (s=0; s<7; s++) {
          PetscInt ii = i + STENCIL[s][0];
          PetscInt jj = j + STENCIL[s][1];
          PetscInt kk = k + STENCIL[s][2];
          PetscInt cc = Color(shape,start,width,ii,jj,kk);
          if (cc == color) nadj++;
        }
      }

  /* Allocate arrays to store the adjacency graph */
  nvtx *= bs; nadj *= bs; /* adjust for block size */
  ierr = PetscMalloc((nvtx+1)*sizeof(PetscInt),&xadj);CHKERRQ(ierr);
  ierr = PetscMalloc(nadj*sizeof(PetscInt),&adjy);CHKERRQ(ierr);

  /* Fill the adjacency graph */
  pos = 0; xadj[pos++] = 0;
  for (k=0; k<shape[2]; k++)
    for (j=0; j<shape[1]; j++)
      for (i=0; i<shape[0]; i++) {
        /* Compute the color of this vertex */
        PetscInt color = Color(shape,start,width,i,j,k);
        /* Compute the list of neighbor vertices
           having the same color */
        PetscInt nv = 0, vertices[7];
        for (s=0; s<7; s++) { /* loop over neighbors */
          PetscInt ii = i + STENCIL[s][0];
          PetscInt jj = j + STENCIL[s][1];
          PetscInt kk = k + STENCIL[s][2];
          PetscInt cc = Color(shape,start,width,ii,jj,kk);
          if (cc == color)
            vertices[nv++] = Index(shape,ii,jj,kk);
        }
        for (c=0; c<bs; c++) {
          xadj[pos] = xadj[pos-1];
          for (v=0; v<nv; v++)
            adjy[xadj[pos]++] = c + bs*vertices[v];
          pos++;
        }
      }

  *_nvtx = nvtx;
  *_xadj = xadj;
  *_adjy = adjy;

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "IGAPreparePCBDDC"
PetscErrorCode IGAPreparePCBDDC(IGA iga,PC pc)
{
  Mat            mat;
  void           (*f)(void);
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
  PetscValidHeaderSpecific(pc,PC_CLASSID,2);
  IGACheckSetUp(iga,1);

  ierr = PCGetOperators(pc,PETSC_NULL,&mat,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscObjectQueryFunction((PetscObject)mat,"MatISGetLocalMat_C",&f);CHKERRQ(ierr);
  if (!f) PetscFunctionReturn(0);
  ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr);
  if (!f) PetscFunctionReturn(0);

#if defined(PETSC_HAVE_PCBDDC)  /* XXX */
  {
    PetscInt i,dim,dof;
    PetscInt shape[3] = {1,1,1};
    PetscInt start[3] = {0,0,0};
    PetscInt width[3] = {1,1,1};
    PetscInt nvtx=0,*xadj=0,*adjy=0;
    ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
    ierr = IGAGetDof(iga,&dof);CHKERRQ(ierr);
    for (i=0; i<dim; i++) {
      shape[i] = iga->node_gwidth[i];
      start[i] = iga->node_lstart[i] - iga->node_gstart[i];
      width[i] = iga->node_lwidth[i];
    }
    ierr = ComputeBDDCGraph(dof,shape,start,width,&nvtx,&xadj,&adjy);CHKERRQ(ierr);
    ierr = PCBDDCSetLocalAdjacencyGraph(pc,nvtx,xadj,adjy,PETSC_OWN_POINTER);CHKERRQ(ierr);
  }
#endif

  PetscFunctionReturn(0);
}

/* ---------------------------------------------------------------- */

#undef  __FUNCT__
#define __FUNCT__ "IGA_OptionsHandler_PC"
static PetscErrorCode IGA_OptionsHandler_PC(PetscObject obj,void *ctx)
{
  PC             pc = (PC)obj;
  Mat            mat;
  IGA            iga;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc,PC_CLASSID,1);
  if (PetscOptionsPublishCount != 1) PetscFunctionReturn(0);
  ierr = PCGetOperators(pc,PETSC_NULL,&mat,PETSC_NULL);CHKERRQ(ierr);
  if (!mat) PetscFunctionReturn(0);
  PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
  ierr = PetscObjectQuery((PetscObject)mat,"IGA",(PetscObject*)&iga);CHKERRQ(ierr);
  if (!iga) PetscFunctionReturn(0);
  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
  /* */
  ierr = IGAPreparePCBDDC(iga,pc);CHKERRQ(ierr);
  /* */
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "IGASetOptionsHandlerPC"
PetscErrorCode IGASetOptionsHandlerPC(PC pc)
{
  PetscErrorCode ierr;
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc,PC_CLASSID,1);
  ierr = PetscObjectAddOptionsHandler((PetscObject)pc,IGA_OptionsHandler_PC,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/* ---------------------------------------------------------------- */