1: #define PETSCDM_DLL
  2: 
 3:  #include petscda.h
 4:  #include petscmat.h


  7: typedef struct _SlicedOps *SlicedOps;
  8: struct _SlicedOps {
  9:   PetscErrorCode (*view)(Sliced,PetscViewer);
 10:   PetscErrorCode (*createglobalvector)(Sliced,Vec*);
 11:   PetscErrorCode (*getcoloring)(Sliced,ISColoringType,ISColoring*);
 12:   PetscErrorCode (*getmatrix)(Sliced,MatType,Mat*);
 13:   PetscErrorCode (*getinterpolation)(Sliced,Sliced,Mat*,Vec*);
 14:   PetscErrorCode (*refine)(Sliced,MPI_Comm,Sliced*);
 15: };

 17: struct _p_Sliced {
 18:   PETSCHEADER(struct _SlicedOps);
 19:   Vec      globalvector;
 20:   PetscInt bs,n,N,Nghosts,*ghosts;
 21:   PetscInt d_nz,o_nz,*d_nnz,*o_nnz;
 22: };

 26: /*@C
 27:     SlicedGetMatrix - Creates a matrix with the correct parallel layout required for 
 28:       computing the Jacobian on a function defined using the informatin in Sliced.

 30:     Collective on Sliced

 32:     Input Parameter:
 33: +   slice - the slice object
 34: -   mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ,
 35:             or any type which inherits from one of these (such as MATAIJ, MATLUSOL, etc.).

 37:     Output Parameters:
 38: .   J  - matrix with the correct nonzero preallocation
 39:         (obviously without the correct Jacobian values)

 41:     Level: advanced

 43:     Notes: This properly preallocates the number of nonzeros in the sparse matrix so you 
 44:        do not need to do it yourself.

 46: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), DASetBlockFills()

 48: @*/
 49: PetscErrorCode  SlicedGetMatrix(Sliced slice, MatType mtype,Mat *J)
 50: {
 51:   PetscErrorCode         ierr,*globals,rstart,i;
 52:   ISLocalToGlobalMapping lmap;

 55:   MatCreate(slice->comm,J);
 56:   MatSetSizes(*J,slice->n,slice->n,PETSC_DETERMINE,PETSC_DETERMINE);
 57:   MatSetType(*J,mtype);
 58:   MatSetBlockSize(*J,slice->bs);
 59:   MatSeqAIJSetPreallocation(*J,slice->d_nz,slice->d_nnz);
 60:   MatMPIAIJSetPreallocation(*J,slice->d_nz,slice->d_nnz,slice->o_nz,slice->o_nnz);
 61:   MatSeqBAIJSetPreallocation(*J,slice->bs,slice->d_nz,slice->d_nnz);
 62:   MatMPIBAIJSetPreallocation(*J,slice->bs,slice->d_nz,slice->d_nnz,slice->o_nz,slice->o_nnz);

 64:   PetscMalloc((slice->n+slice->Nghosts+1)*sizeof(PetscInt),&globals);
 65:   MatGetOwnershipRange(*J,&rstart,PETSC_NULL);
 66:   for (i=0; i<slice->n; i++) {
 67:     globals[i] = rstart + i;
 68:   }
 69:   PetscMemcpy(globals+slice->n,slice->ghosts,slice->Nghosts*sizeof(PetscInt));
 70:   ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,slice->n+slice->Nghosts,globals,&lmap);
 71:   PetscFree(globals);
 72:   MatSetLocalToGlobalMapping(*J,lmap);
 73:   ISLocalToGlobalMappingDestroy(lmap);
 74:   return(0);
 75: }

 79: /*@C
 80:     SlicedSetGhosts - Sets the global indices of other processes elements that will
 81:       be ghosts on this process

 83:     Not Collective

 85:     Input Parameters:
 86: +    slice - the Sliced object
 87: .    bs - block size
 88: .    nlocal - number of local (non-ghost) entries
 89: .    Nghosts - number of ghosts on this process
 90: -    ghosts - indices of all the ghost points

 92:     Level: advanced

 94: .seealso SlicedDestroy(), SlicedCreateGlobalVector(), SlicedGetGlobalIndices()

 96: @*/
 97: PetscErrorCode  SlicedSetGhosts(Sliced slice,PetscInt bs,PetscInt nlocal,PetscInt Nghosts,const PetscInt ghosts[])
 98: {

103:   PetscFree(slice->ghosts);
104:   PetscMalloc((1+Nghosts)*sizeof(PetscInt),&slice->ghosts);
105:   PetscMemcpy(slice->ghosts,ghosts,Nghosts*sizeof(PetscInt));
106:   slice->bs      = bs;
107:   slice->n       = nlocal;
108:   slice->Nghosts = Nghosts;
109:   return(0);
110: }

114: /*@C
115:     SlicedSetPreallocation - sets the matrix memory preallocation for matrices computed by Sliced

117:     Not Collective

119:     Input Parameters:
120: +    slice - the Sliced object
121: .    d_nz - maximum number of nonzeros in any row of diagonal block
122: .    d_nnz - number of nonzeros in each row of diagonal block
123: .    o_nz - maximum number of nonzeros in any row of off-diagonal block
124: .    o_nnz - number of nonzeros in each row of off-diagonal block


127:     Level: advanced

129: .seealso SlicedDestroy(), SlicedCreateGlobalVector(), SlicedGetGlobalIndices(), MatMPIAIJSetPreallocation(),
130:          MatMPIBAIJSetPreallocation()

132: @*/
133: PetscErrorCode  SlicedSetPreallocation(Sliced slice,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
134: {
137:   slice->d_nz  = d_nz;
138:   slice->d_nnz = (PetscInt*)d_nnz;
139:   slice->o_nz  = o_nz;
140:   slice->o_nnz = (PetscInt*)o_nnz;
141:   return(0);
142: }

146: /*@C
147:     SlicedCreate - Creates a DM object, used to manage data for a unstructured problem

149:     Collective on MPI_Comm

151:     Input Parameter:
152: .   comm - the processors that will share the global vector

154:     Output Parameters:
155: .   slice - the slice object

157:     Level: advanced

159: .seealso SlicedDestroy(), SlicedCreateGlobalVector(), SlicedGetGlobalIndices()

161: @*/
162: PetscErrorCode  SlicedCreate(MPI_Comm comm,Sliced *slice)
163: {
165:   Sliced         p;

169:   *slice = PETSC_NULL;
170: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
171:   DMInitializePackage(PETSC_NULL);
172: #endif

174:   PetscHeaderCreate(p,_p_Sliced,struct _SlicedOps,DA_COOKIE,0,"Sliced",comm,SlicedDestroy,0);
175:   p->ops->createglobalvector = SlicedCreateGlobalVector;
176:   p->ops->getmatrix          = SlicedGetMatrix;
177:   *slice = p;
178:   return(0);
179: }

183: /*@C
184:     SlicedDestroy - Destroys a vector slice.

186:     Collective on Sliced

188:     Input Parameter:
189: .   slice - the slice object

191:     Level: advanced

193: .seealso SlicedCreate(), SlicedCreateGlobalVector(), SlicedGetGlobalIndices()

195: @*/
196: PetscErrorCode  SlicedDestroy(Sliced slice)
197: {
198:   PetscErrorCode     ierr;

201:   if (--slice->refct > 0) return(0);
202:   if (slice->globalvector) {VecDestroy(slice->globalvector);}
203:   PetscHeaderDestroy(slice);
204:   return(0);
205: }

209: /*@C
210:     SlicedCreateGlobalVector - Creates a vector of the correct size to be gathered into 
211:         by the slice.

213:     Collective on Sliced

215:     Input Parameter:
216: .    slice - the slice object

218:     Output Parameters:
219: .   gvec - the global vector

221:     Level: advanced

223:     Notes: Once this has been created you cannot add additional arrays or vectors to be packed.

225: .seealso SlicedDestroy(), SlicedCreate(), SlicedGetGlobalIndices()

227: @*/
228: PetscErrorCode  SlicedCreateGlobalVector(Sliced slice,Vec *gvec)
229: {
230:   PetscErrorCode     ierr;


234:   if (slice->globalvector) {
235:     VecDuplicate(slice->globalvector,gvec);
236:   } else {
237:     VecCreateGhostBlock(slice->comm,slice->bs,slice->n,PETSC_DETERMINE,slice->Nghosts,slice->ghosts,&slice->globalvector);
238:     *gvec = slice->globalvector;
239:     PetscObjectReference((PetscObject)*gvec);
240:   }
241:   return(0);
242: }

246: /*@C
247:     SlicedGetGlobalIndices - Gets the global indices for all the local entries

249:     Collective on Sliced

251:     Input Parameter:
252: .    slice - the slice object

254:     Output Parameters:
255: .    idx - the individual indices for each packed vector/array
256:  
257:     Level: advanced

259:     Notes:
260:        The idx parameters should be freed by the calling routine with PetscFree()

262: .seealso SlicedDestroy(), SlicedCreateGlobalVector(), SlicedCreate()

264: @*/
265: PetscErrorCode  SlicedGetGlobalIndices(Sliced slice,PetscInt *idx[])
266: {
267:   return(0);
268: }