Viewing fully implemented for dRule.

Signed-off-by: Jed Brown <jed@59A2.org>

include/dohpjacobi.h

@@ -77,6 +77,11 @@ EXTERN dErr dJacobiSetDegrees(dJacobi,dInt,dInt);
 EXTERN dErr dJacobiGetRule(dJacobi jac,dTopology top,const dInt rsize[],dRule *rule,void **base,dInt *index);
 EXTERN dErr dJacobiGetEFS(dJacobi jac,dTopology top,const dInt bsize[],dRule *rule,dEFS *efs,void **base,dInt *index);
 
+EXTERN dErr dRuleView(dRule *rule,PetscViewer);
+EXTERN dErr dRuleGetSize(dRule *rule,dInt *dim,dInt *nnodes);
+EXTERN dErr dRuleGetNodeWeight(dRule *rule,dReal *coord,dReal *weight);
+EXTERN dErr dRuleGetTensorNodeWeight(dRule *rule,dInt *dim,dInt *nnodes,const dReal **coord,const dReal **weight);
+
 PETSC_EXTERN_CXX_END
 
 #endif /* _DOHPJACOBI_H */

include/private/fsimpl.h

@@ -22,10 +22,6 @@ struct p_dRule {
   struct v_dRuleOps *ops;
   void              *data;
 };
-#define dRuleView(rule,view)                                    (*rule->ops->view)(rule,view)
-#define dRuleGetSize(rule,dim,nnodes)                           (*rule->ops->getSize)(rule,dim,nnodes)
-#define dRuleGetNodeWeight(rule,coord,weight)                   (*rule->ops->getNodeWeight)(rule,coord,weight)
-#define dRuleGetTensorNodeWeight(rule,dim,nnodes,coord,weights) (*rule->ops->getTensorNodeWeight)(rule,dim,nnodes,coord,weights)
 
 /**
 * Operations required for an EFS.  Defined here so that these function calls can be inlined.

src/jacobi/impls/tensor/efstopo.c

@@ -1,5 +1,6 @@
 #include "tensor.h"
 
+static dErr dEFSView_Tensor_Private(const char *,dRule*,dInt,TensorBasis*,PetscViewer);
 #ifdef _F
 # undef _F
 #endif
@@ -14,23 +15,25 @@ _F(dEFSGetSizes_Tensor_Quad);
 _F(dEFSGetSizes_Tensor_Hex);
 #undef _F
 #define _F(f) static dErr f(dEFS*,dInt,dInt*,dScalar**restrict,const dScalar[],dScalar[],dApplyMode,InsertMode) /* dEFSApply */
-_F(dEFSApply_Tensor_Line);
-_F(dEFSApply_Tensor_Quad);
+//_F(dEFSApply_Tensor_Line);
+//_F(dEFSApply_Tensor_Quad);
 _F(dEFSApply_Tensor_Hex);
 #undef _F
-#define _F(f) static dErr f(dEFS*,dInt,dInt,const dScalar[],dScalar[],InsertMode,ScatterMode) /* dEFSScatterInt */
+#if 0
+# define _F(f) static dErr f(dEFS*,dInt,dInt,const dScalar[],dScalar[],InsertMode,ScatterMode) /* dEFSScatterInt */
 _F(dEFSScatterInt_Tensor_Line);
 _F(dEFSScatterInt_Tensor_Quad);
 _F(dEFSScatterInt_Tensor_Hex);
-#undef _F
-#define _F(f) static dErr f(dEFS,dEFS,dInt*,dScalar**restrict,const dScalar[],dScalar[],InsertMode,ScatterMode) /* dEFSScatterFacet */
+# undef _F
+# define _F(f) static dErr f(dEFS,dEFS,dInt*,dScalar**restrict,const dScalar[],dScalar[],InsertMode,ScatterMode) /* dEFSScatterFacet */
 _F(dEFSScatterFacet_Tensor_Line);
 _F(dEFSScatterFacet_Tensor_Quad);
 _F(dEFSScatterFacet_Tensor_Hex);
-#undef _F
-#define _F(f) static dErr f(dInt D,const dInt P[],const dInt Q[],dInt *wlen,dScalar *restrict* work,dReal *A[],dTransposeMode tpose[],const dScalar f[],dScalar g[restrict],InsertMode imode)
-_F(TensorMult_Line);
-_F(TensorMult_Quad);
+# undef _F
+#endif
+#define _F(f) static dErr f(dInt D,const dInt P[],const dInt Q[],dInt *wlen,dScalar **restrict work,dReal *A[],dTransposeMode tpose[],const dScalar f[],dScalar g[restrict],InsertMode imode)
+//_F(TensorMult_Line);
+//_F(TensorMult_Quad);
 _F(TensorMult_Hex);
 #undef _F
 
@@ -46,12 +49,12 @@ dErr dJacobiEFSOpsSetUp_Tensor(dJacobi jac)
 {
   static const struct v_dEFSOps efsOpsLine = { .view = dEFSView_Tensor_Line,
                                                .getSizes = dEFSGetSizes_Tensor_Line,
-                                               .apply = dEFSApply_Tensor_Line,
+                                               .apply = 0, /* dEFSApply_Tensor_Line, */
                                                .scatterInt = 0,
                                                .scatterFacet = 0 };
   static const struct v_dEFSOps efsOpsQuad = { .view = dEFSView_Tensor_Quad,
                                                .getSizes = dEFSGetSizes_Tensor_Quad,
-                                               .apply = dEFSApply_Tensor_Quad,
+                                               .apply = 0, /* dEFSApply_Tensor_Quad, */
                                                .scatterInt = 0,
                                                .scatterFacet = 0 };
   static const struct v_dEFSOps efsOpsHex  = { .view = dEFSView_Tensor_Hex,
@@ -78,6 +81,130 @@ dErr dJacobiEFSOpsSetUp_Tensor(dJacobi jac)
   dFunctionReturn(0);
 }
 
+static dErr dEFSView_Tensor_Private(const char *name,dRule *rule,dInt n,TensorBasis *b,PetscViewer viewer)
+{
+  dBool ascii;
+  dErr err;
+
+  dFunctionBegin;
+  err = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&ascii);dCHK(err);
+  if (ascii) {
+    err = PetscViewerASCIIPrintf(viewer,"dEFS type %s\n",name);dCHK(err);
+    err = PetscViewerASCIIPrintf(viewer,"based on dRule:\n");dCHK(err);
+    err = PetscViewerASCIIPushTab(viewer);dCHK(err);
+    err = dRuleView(rule,viewer);dCHK(err);
+    err = PetscViewerASCIIPopTab(viewer);
+    for (dInt i=0; i<n; i++) {
+      err = TensorBasisView(b[i],viewer);
+    }
+  }
+  dFunctionReturn(0);
+}
+
+
+static dErr dEFSView_Tensor_Line(dEFS *efs,PetscViewer viewer)
+{
+  dErr err;
+
+  dFunctionBegin;
+  err = dEFSView_Tensor_Private("Tensor_Line",efs->rule,1,((struct s_dEFS_Tensor_Line*)efs->data)->basis,viewer);dCHK(err);
+  dFunctionReturn(0);
+}
+static dErr dEFSView_Tensor_Quad(dEFS *efs,PetscViewer viewer)
+{
+  dErr err;
+
+  dFunctionBegin;
+  err = dEFSView_Tensor_Private("Tensor_Quad",efs->rule,2,((struct s_dEFS_Tensor_Quad*)efs->data)->basis,viewer);dCHK(err);
+  dFunctionReturn(0);
+}
+static dErr dEFSView_Tensor_Hex(dEFS *efs,PetscViewer viewer)
+{
+  dErr err;
+
+  dFunctionBegin;
+  err = dEFSView_Tensor_Private("Tensor_Hex",efs->rule,3,((struct s_dEFS_Tensor_Hex*)efs->data)->basis,viewer);dCHK(err);
+  dFunctionReturn(0);
+}
+
+static dErr dEFSGetSizes_Tensor_Line(dEFS *efs,dInt *dim,dInt *inodes,dInt *total)
+{
+  TensorBasis *b = ((struct s_dEFS_Tensor_Line*)efs->data)->basis;
+
+  dFunctionBegin;
+  if (dim) *dim = 1;
+  if (inodes) *inodes = b[0]->P - 2;
+  if (total) *total = b[0]->P;
+  dFunctionReturn(0);
+}
+static dErr dEFSGetSizes_Tensor_Quad(dEFS *efs,dInt *dim,dInt *inodes,dInt *total)
+{
+  TensorBasis *b = ((struct s_dEFS_Tensor_Quad*)efs->data)->basis;
+
+  dFunctionBegin;
+  if (dim) *dim = 2;
+  if (inodes) *inodes = (b[0]->P - 2) * (b[1]->P - 2);
+  if (total) *total = b[0]->P * b[1]->P;
+  dFunctionReturn(0);
+}
+static dErr dEFSGetSizes_Tensor_Hex(dEFS *efs,dInt *dim,dInt *inodes,dInt *total)
+{
+  TensorBasis *b = ((struct s_dEFS_Tensor_Hex*)efs->data)->basis;
+
+  dFunctionBegin;
+  if (dim) *dim = 3;
+  if (inodes) *inodes = (b[0]->P - 2) * (b[1]->P - 2) * (b[2]->P - 2);
+  if (total) *total = b[0]->P * b[1]->P * b[2]->P;
+  dFunctionReturn(0);
+}
+
+static dErr dEFSApply_Tensor_Hex(dEFS *efs,dInt D,dInt *wlen,dScalar**restrict work,const dScalar in[],dScalar out[],dApplyMode amode,InsertMode imode)
+{
+  TensorBasis *b = ((struct s_dEFS_Tensor_Hex*)efs->data)->basis;
+  /* const dInt *P=&b->P,*Q=&b->Q; */
+  dScalar *A[3];
+  dInt P[3],Q[3],chunk;
+  dTransposeMode tpose[3];
+  dErr err;
+
+  dFunctionBegin;
+  for (dInt i=0; i<3; i++) {
+    P[i] = b[i]->P;
+    Q[i] = b[i]->Q;
+  }
+  switch (amode) {
+    case dAPPLY_INTERP:
+      for (dInt i=0; i<3; i++) {
+        A[i] = b[i]->interp;
+        tpose[i] = dTRANSPOSE_NO;
+      }
+      err = TensorMult_Hex(D,P,Q,wlen,work,A,tpose,in,out,imode);dCHK(err);
+      break;
+    case dAPPLY_INTERP_TRANSPOSE:
+      for (dInt i=0; i<3; i++) {
+        A[i] = b[i]->interp;
+        tpose[i] = dTRANSPOSE_YES;
+      }
+      err = TensorMult_Hex(D,Q,P,wlen,work,A,tpose,in,out,imode);dCHK(err);
+      break;
+    case dAPPLY_GRAD:
+      chunk = D*Q[0]*Q[1]*Q[2]; /* number of dofs for each component of the gradient */
+      for (dInt i=0; i<3; i++) {
+        for (dInt j=0; j<3; j++) {
+          if (i == j) A[j] = b[j]->deriv;
+          else A[j] = b[j]->interp;
+          tpose[j] = dTRANSPOSE_NO;
+        }
+        err = TensorMult_Hex(D,P,Q,wlen,work,A,tpose,in,&out[i*chunk],imode);dCHK(err);
+      }
+      break;
+    case dAPPLY_GRAD_TRANSPOSE:
+    default:
+      dERROR(1,"invalid dApplyMode %d specified",amode);
+  }
+  dFunctionReturn(0);
+}
+
 
 /** 
 * The core computational kernel.  Performs a tensor product operation with the matrices A,B,C.
@@ -95,7 +222,7 @@ dErr dJacobiEFSOpsSetUp_Tensor(dJacobi jac)
 * 
 * @return err
 */
-static dErr TensorMult_Hex(dInt D,const dInt P[3],const dInt Q[3],dInt *wlen,dScalar *restrict* work,
+static dErr TensorMult_Hex(dInt D,const dInt P[3],const dInt Q[3],dInt *wlen,dScalar **restrict work,
                            dReal *A[3],dTransposeMode tpose[3],const dScalar f[],dScalar g[restrict],InsertMode imode)
 {
   dInt i,j,k,l,d,idx;

src/jacobi/impls/tensor/tensor.c

@@ -189,7 +189,7 @@ static dErr dRealTableView(dInt m,dInt n,const dReal mat[],const char *name,dVie
   if (!ascii) dFunctionReturn(0);
   for (dInt i=0; i<m; i++) {
     if (name) {
-      err = PetscViewerASCIIPrintf(viewer,"%10s[%2d][%2d:%2d] ",name,i,0,n-1);dCHK(err);
+      err = PetscViewerASCIIPrintf(viewer,"%10s[%2d][%2d:%2d] ",name,i,0,n);dCHK(err);
     }
     err = PetscViewerASCIIUseTabs(viewer,PETSC_NO);dCHK(err);
     for (dInt j=0; j<n; j++) {
@@ -246,7 +246,7 @@ static dErr dJacobiGetRule_Tensor(dJacobi jac,dTopology top,const dInt rsize[],d
       if (base) {
         rule->ops = this->ruleOpsQuad;
         rule->data = start;
-        quad = (struct s_dRule_Tensor_Quad*)&start;
+        quad = (struct s_dRule_Tensor_Quad*)start;
         err = TensorGetRule(this,rsize[0],(TensorRule*)&quad->rule[0]);dCHK(err);
         err = TensorGetRule(this,rsize[1],(TensorRule*)&quad->rule[1]);dCHK(err);
       }

src/jacobi/interface/jacobi.c

@@ -322,3 +322,44 @@ dErr dJacobiGetEFS(dJacobi jac,dTopology top,const dInt bsize[],dRule *rule,dEFS
   err = jac->ops->getefs(jac,top,bsize,rule,efs,base,index);dCHK(err);
   dFunctionReturn(0);
 }
+
+dErr dRuleView(dRule *rule,PetscViewer viewer)
+{
+  dErr err;
+
+  dFunctionBegin;
+  dValidPointer(rule,1);
+  dValidHeader(viewer,PETSC_VIEWER_COOKIE,2);
+  err = (*rule->ops->view)(rule,viewer);dCHK(err);
+  dFunctionReturn(0);
+}
+
+dErr dRuleGetSize(dRule *rule,dInt *dim,dInt *nnodes)
+{
+  dErr err;
+
+  dFunctionBegin;
+  dValidPointer(rule,1);
+  err = (*rule->ops->getSize)(rule,dim,nnodes);dCHK(err);
+  dFunctionReturn(0);
+}
+
+dErr dRuleGetNodeWeight(dRule *rule,dReal *coord,dReal *weight)
+{
+  dErr err;
+
+  dFunctionBegin;
+  dValidPointer(rule,1);
+  err = (*rule->ops->getNodeWeight)(rule,coord,weight);dCHK(err);
+  dFunctionReturn(0);
+}
+
+dErr dRuleGetTensorNodeWeight(dRule *rule,dInt *dim,dInt *nnodes,const dReal **coord,const dReal **weight)
+{
+  dErr err;
+
+  dFunctionBegin;
+  dValidPointer(rule,1);
+  err = (*rule->ops->getTensorNodeWeight)(rule,dim,nnodes,coord,weight);dCHK(err);
+  dFunctionReturn(0);
+}

src/jacobi/tests/ex1.c

@@ -1,6 +1,7 @@
 static const char help[] = "Tests the dJacobi object.";
 
 #include "dohpjacobi.h"
+#include "private/fsimpl.h"
 
 dErr checkRulesAndEFS(dJacobi);
 
@@ -21,9 +22,9 @@ int main(int argc,char *argv[])
   err = dJacobiSetDegrees(jac,8,4);dCHK(err);
   err = dJacobiSetFromOptions(jac);dCHK(err);
   err = dJacobiSetUp(jac);dCHK(err);
-  err = dJacobiView(jac,viewer);dCHK(err);
+  /* err = dJacobiView(jac,viewer);dCHK(err); */
   err = dJacobiSetUp(jac);dCHK(err);
-  err = getRules(jac);dCHK(err);  
+  err = checkRulesAndEFS(jac);dCHK(err);  
   err = dJacobiDestroy(jac);dCHK(err);
   err = PetscFinalize();dCHK(err);
   dFunctionReturn(0);
@@ -39,6 +40,7 @@ dErr checkRulesAndEFS(dJacobi jac)
   MPI_Comm comm = ((PetscObject)jac)->comm;
   dRule *rule;
   dEFS *efs;
+  dInt index;
   void **rbase,**ebase;
   dErr err;
 
@@ -58,7 +60,7 @@ dErr checkRulesAndEFS(dJacobi jac)
 
   for (dInt i=0; i<N; i++) {
     err = dPrintf(comm,"Rule for element %d\n",i);dCHK(err);
-    err = dRuleView(&rule,viewer);dCHK(err);
+    err = dRuleView(&rule[i],viewer);dCHK(err);
   }
   dFunctionReturn(0);
 }