Remove the block size attribute from dFS, use the one that was promot…

…ed to DM

include/dohpfsimpl.h

@@ -101,7 +101,6 @@ struct _p_dFS {
   dBool        spacebuilt;
   dBool        assemblefull;    /**< Use full order constraints for assembly */
   dBool        assemblereduced; /**< Assemble only diagonal part of blocks, only matters for bs>1 and MATAIJ */
-  dInt         bs;              /**< Block size (number of dofs per node) */
   dInt         nelem;
   dInt        *off;             /**< Offset of element dofs in expanded vector */
   struct _dFSIntegrationLink *integration;

src/fs/impls/cont/cont.c

@@ -388,7 +388,7 @@ static dErr dFSGetSubElementMeshSize_Cont(dFS fs,dInt *nelem,dInt *nvert,dInt *n
 {
   /* dFS_Cont *cont = fs->data; */
   dErr    err;
-  dInt    nents,nsub,n;
+  dInt    nents,nsub,n,bs;
   dMeshEH *ents;
   dPolynomialOrder *degree;
   Vec     X;
@@ -407,8 +407,9 @@ static dErr dFSGetSubElementMeshSize_Cont(dFS fs,dInt *nelem,dInt *nvert,dInt *n
   err = VecGetLocalSize(X,&n);dCHK(err);
   err = VecDohpRestoreClosure(fs->gvec,&X);dCHK(err);
   err = dFree2(ents,degree);dCHK(err);
+  err = dFSGetBlockSize(fs,&bs);dCHK(err);
   *nelem = nsub;
-  *nvert = n/fs->bs;
+  *nvert = n/bs;
   *nconn = nsub*8;              /* all hex */
   dFunctionReturn(0);
 }

src/fs/impls/cont/contviewdhm.c

@@ -110,7 +110,7 @@ dErr dFSView_Cont_DHM(dFS fs,dViewer viewer)
     dht_FS     fs5;
     dht_Field *field5;
     hid_t      h5t_fs,fsdset,fsspace;
-    dInt       i;
+    dInt       i,bs;
     PetscMPIInt size;
 
     err = dViewerDHMGetFSType(viewer,&h5t_fs);dCHK(err);
@@ -121,12 +121,13 @@ dErr dFSView_Cont_DHM(dFS fs,dViewer viewer)
     err = dMeshGetTagName(fs->mesh,fs->tag.bstatus,&fs5.bstatus);dCHK(err);
     herr = H5Rcreate(&fs5.mesh,meshgrp,mstatestr,H5R_OBJECT,-1);dH5CHK(herr,H5Rcreate);
     fs5.time = dhm->time;
+    err = dFSGetBlockSize(fs,&bs);dCHK(err);
     err = PetscObjectStateQuery((PetscObject)fs,&fs5.internal_state);dCHK(err);
     err = MPI_Comm_size(((dObject)fs)->comm,&size);dCHK(err);
     fs5.number_of_subdomains = size;
-    fs5.fields.len = fs->bs;
+    fs5.fields.len = bs;
     err = dMallocA(fs5.fields.len,&field5);dCHK(err);
-    for (i=0; i<fs->bs; i++) {
+    for (i=0; i<bs; i++) {
       field5[i].name = fs->fieldname[i];
       field5[i].units.dimensions = (char*)"m s-1"; /* we only use it as \c const */
       field5[i].units.scale = exp(1.0);
@@ -323,7 +324,7 @@ dErr dFSLoadIntoFS_Cont_DHM(PetscViewer viewer,const char name[],dFS fs)
         dMeshTag tag;
         dJacobi jac;
         dMeshESH set,*sets;
-        dInt nsets;
+        dInt nsets,bs;
         dIInt readrank = 0;         /* Hard-code the rank for now */
         err = dMeshGetTag(mesh,fs5.partition,&tag);dCHK(err);
         err = dMeshGetTaggedSet(mesh,tag,&readrank,&set);dCHK(err);
@@ -345,7 +346,8 @@ dErr dFSLoadIntoFS_Cont_DHM(PetscViewer viewer,const char name[],dFS fs)
         err = dFree(sets);dCHK(err);
 
         err = dFSSetBlockSize(fs,(dInt)fs5.fields.len);dCHK(err);
-        for (dInt i=0; i<fs->bs; i++) {
+        err = dFSGetBlockSize(fs,&bs);dCHK(err);
+        for (dInt i=0; i<bs; i++) {
           const dht_Field *field5 = fs5.fields.p;
           err = dFSSetFieldName(fs,i,field5[i].name);dCHK(err);
         }

src/fs/interface/fs.c

@@ -95,22 +95,25 @@ dErr dFSSetDegree(dFS fs,dJacobi jac,dMeshTag deg)
 dErr dFSSetBlockSize(dFS fs,dInt bs)
 {
   dErr err;
+  dInt obs;
 
   dFunctionBegin;
   dValidHeader(fs,DM_CLASSID,1);
-  for (dInt i=0; i<fs->bs; i++) {err = dFree(fs->fieldname[i]);dCHK(err);}
+  err = dFSGetBlockSize(fs,&obs);dCHK(err);
+  for (dInt i=0; i<obs; i++) {err = dFree(fs->fieldname[i]);dCHK(err);}
   err = dFree(fs->fieldname);dCHK(err);
   err = dCallocA(bs,&fs->fieldname);dCHK(err);
-  fs->bs = bs;
+  fs->dm.bs = bs;
   dFunctionReturn(0);
 }
 
 dErr dFSGetBlockSize(dFS fs,dInt *bs)
 {
+  dErr err;
   dFunctionBegin;
   dValidHeader(fs,DM_CLASSID,1);
   dValidIntPointer(bs,2);
-  *bs = fs->bs;
+  err = DMGetBlockSize((DM)fs,bs);dCHK(err);
   dFunctionReturn(0);
 }
 
@@ -125,10 +128,12 @@ dErr dFSGetBlockSize(dFS fs,dInt *bs)
 dErr dFSSetFieldName(dFS fs,dInt fn,const char *fname)
 {
   dErr err;
+  dInt bs;
 
   dFunctionBegin;
   dValidHeader(fs,DM_CLASSID,1);
-  if (fn < 0 || fs->bs <= fn) dERROR(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field number %d out of range",fn);
+  err = dFSGetBlockSize(fs,&bs);dCHK(err);
+  if (fn < 0 || bs <= fn) dERROR(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field number %d out of range",fn);
   if (fs->fieldname[fn]) {err = dFree(fs->fieldname[fn]);dCHK(err);}
   err = PetscStrallocpy(fname,&fs->fieldname[fn]);dCHK(err);
   dFunctionReturn(0);
@@ -169,11 +174,13 @@ dErr dFSRegisterBoundarySet(dFS fs,dMeshESH bset,dFSBStatus bstat,dFSConstraintF
   iMesh_Instance   mi;
   dErr             err;
   dIInt            ierr;
+  dInt             bs;
 
   dFunctionBegin;
   dValidHeader(fs,DM_CLASSID,1);
+  err = dFSGetBlockSize(fs,&bs);dCHK(err);
   if (!dFSBStatusValid(bstat)) dERROR(PETSC_COMM_SELF,1,"Boundary status %x invalid",bstat);
-  if (dFSBStatusStrongCount(bstat) > fs->bs) dERROR(PETSC_COMM_SELF,1,"Cannot impose strong conditions on more dofs than the block size");
+  if (dFSBStatusStrongCount(bstat) > bs) dERROR(PETSC_COMM_SELF,1,"Cannot impose strong conditions on more dofs than the block size");
   err = dMeshTagSSetData(fs->mesh,fs->tag.bstatus,&bset,1,&bstat,sizeof(bstat),dDATA_BYTE);dCHK(err);
   if (cfunc) {
     struct dFSConstraintCtx ctx;
@@ -225,19 +232,20 @@ dErr dFSView(dFS fs,dViewer viewer)
       err = PetscViewerASCIIPrintf(viewer,"number of vertices=%d edges=%d faces=%d regions=%d\n",nents[0],nents[1],nents[2],nents[3]);dCHK(err);
     }
     {                           /* print aggregate sizes */
-      dInt lm[4],gm[4];
+      dInt lm[4],gm[4],bs;
       err = MatGetSize(fs->E,&lm[0],&lm[1]);dCHK(err);
-      if (lm[0]%fs->bs || lm[1]%fs->bs) dERROR(PETSC_COMM_SELF,1,"Constraint matrix not a multiple of block size, should not happen");
-      lm[0] /= fs->bs;
-      lm[1] /= fs->bs;
+      err = dFSGetBlockSize(fs,&bs);dCHK(err);
+      if (lm[0]%bs || lm[1]%bs) dERROR(PETSC_COMM_SELF,1,"Constraint matrix not a multiple of block size, should not happen");
+      lm[0] /= bs;
+      lm[1] /= bs;
       if (lm[1] != fs->nc) dERROR(PETSC_COMM_SELF,1,"Inconsistent number of closure nodes");
       lm[2] = fs->n;
       lm[3] = fs->ngh;
       err = MPI_Reduce(lm,gm,4,MPIU_INT,MPI_SUM,0,((dObject)fs)->comm);dCHK(err);
       err = PetscViewerASCIIPrintf(viewer,"On rank 0: %d/%d expanded nodes constrained against %d+%d / %d+%d real nodes, %d / %d closure\n",
                                    lm[0],gm[0], lm[2],lm[3], gm[2],gm[3], lm[1],gm[1]);dCHK(err);
       err = PetscViewerASCIIPrintf(viewer,"Block size %d: global dofs %d, ghost dofs %d, closure dofs %d\n",
-                                   fs->bs,fs->bs*gm[2],fs->bs*gm[3],fs->bs*gm[1]);dCHK(err);
+                                   bs,bs*gm[2],bs*gm[3],bs*gm[1]);dCHK(err);
     }
     if (fs->ops->view) {
       err = (*fs->ops->view)(fs,viewer);dCHK(err);
@@ -282,13 +290,15 @@ dErr DMDestroy_dFS(DM dm)
 {
   dFS fs = (dFS)dm;
   dErr err;
+  dInt bs;
 
   dFunctionBegin;
   dValidHeader(fs,DM_CLASSID,1);
   if (fs->ops->impldestroy) {
     err = (*fs->ops->impldestroy)(fs);dCHK(err);
   }
-  for (dInt i=0; i<fs->bs; i++) {err = dFree(fs->fieldname[i]);dCHK(err);}
+  err = dFSGetBlockSize(fs,&bs);dCHK(err);
+  for (dInt i=0; i<bs; i++) {err = dFree(fs->fieldname[i]);dCHK(err);}
   err = dFree(fs->fieldname);dCHK(err);
   err = VecDestroy(&fs->gvec);dCHK(err);
   err = VecDestroy(&fs->dcache);dCHK(err);
@@ -655,7 +665,8 @@ dErr dFSGetMatrix(dFS fs,const MatType mtype,Mat *inJ)
   dValidCharPointer(mtype,2);
   dValidPointer(inJ,3);
   *inJ = 0;
-  bs = fs->bs; n = fs->n;
+  n = fs->n;
+  err = dFSGetBlockSize(fs,&bs);dCHK(err);
   err = MatCreate(((dObject)fs)->comm,&J);dCHK(err);
   err = MatSetSizes(J,bs*n,bs*n,PETSC_DETERMINE,PETSC_DETERMINE);dCHK(err);
   err = MatSetType(J,mtype);dCHK(err);
@@ -706,7 +717,7 @@ dErr dFSMatSetValuesBlockedExpanded(dFS fs,Mat A,dInt m,const dInt idxm[],dInt n
   dValidPointer(idxm,4);
   dValidPointer(idxn,6);
   dValidPointer(v,7);
-  bs = fs->bs;
+  err = dFSGetBlockSize(fs,&bs);dCHK(err);
   err = PetscLogEventBegin(dLOG_FSMatSetValuesExpanded,fs,A,0,0);dCHK(err);
   err = MatMAIJGetAIJ(fs->assemblefull?fs->E:fs->Ep,&E);dCHK(err); /* Does not reference so do not destroy or return E */
   err = MatGetRowIJ(E,0,dFALSE,dFALSE,&cn,&ci,&cj,&done);dCHK(err);

src/fs/interface/fsbuild.c

@@ -164,12 +164,13 @@ dErr dFSBuildSpaceVectors_Private(dFS fs,dMeshTag indexTag,const dInt inodes[],d
 
   /* Set ghost indices of every node using \a ghidx, create global vector. */
   {
-    dInt gh=0;
+    dInt gh=0,bs;
     for (dInt i=0; i<ghents_s; i++) {
       for (dInt j=0; j<inodes[idx[i]]; j++) ghidx[gh++] = gcoffset[i] + j;
     }
     if (gh != fs->ngh) dERROR(PETSC_COMM_SELF,1,"Ghost count inconsistent");
-    err = VecCreateDohp(((dObject)fs)->comm,fs->bs,fs->n,fs->nc,fs->ngh,ghidx,&fs->gvec);dCHK(err);
+    err = dFSGetBlockSize(fs,&bs);dCHK(err);
+    err = VecCreateDohp(((dObject)fs)->comm,bs,fs->n,fs->nc,fs->ngh,ghidx,&fs->gvec);dCHK(err);
   }
   err = dFree3(gcoffset,ghidx,idx);dCHK(err);
 

src/fs/interface/fsreg.c

@@ -58,7 +58,7 @@ dErr dFSCreate(MPI_Comm comm,dFS *infs)
   err = PetscNewLog(fs,struct _dFSOps,&fs->ops);dCHK(err);
 
 /* Defaults */
-  fs->bs = 1;
+  fs->dm.bs = 1;
   err = dCallocA(1,&fs->fieldname);dCHK(err);
 
   *infs = fs;

src/fs/interface/fsrot.c

@@ -51,7 +51,8 @@ dErr dFSRotationCreate(dFS fs,IS is,dReal rmat[],dInt ns[],Vec v,dFSRotation *in
   *inrot = 0;
   err = PetscHeaderCreate(rot,_p_dFSRotation,struct _dFSRotationOps,dFSROT_CLASSID,0,"dFSRotation",PETSC_COMM_SELF,dFSRotationDestroy,dFSRotationView);dCHK(err);
 
-  bs = rot->bs = fs->bs;
+  err = dFSGetBlockSize(fs,&bs);dCHK(err);
+  rot->bs = bs;
   err = ISGetSize(is,&n);dCHK(err);
   rot->n = n;
   err = PetscObjectReference((PetscObject)is);dCHK(err);