Add dRuleGetPatches to enable integration over separate patches withi…

…n an element

Includes simple test program

include/dohpjacimpl.h

@@ -20,6 +20,7 @@ struct _dRuleOps {
   dErr (*view)(dRule,dViewer);
   dErr (*getSize)(dRule,dInt*,dInt*); /**< topological dimension of the space, total number of nodes */
   dErr (*getNodeWeight)(dRule,dReal[],dReal[]); /**< nodes and weights in interlaced ordering, arrays must be large enough */
+  dErr (*getPatches)(dRule,dInt*,dInt*,const dInt **,const dReal**);
   dErr (*getTensorNodeWeight)(dRule,dInt*,dInt[],const dReal**,const dReal**); /**< topological dimension, number of
                                                                                * nodes in each direction, weights in
                                                                                * each direction.  Does not copy, may not

include/dohpjacobi.h

@@ -170,6 +170,7 @@ extern dErr dJacobiAddConstraints(dJacobi,dInt,const dInt[],const dInt[],const d
 
 extern dErr dRuleView(dRule rule,dViewer);
 extern dErr dRuleGetSize(dRule rule,dInt *dim,dInt *nnodes);
+extern dErr dRuleGetPatches(dRule rule,dInt *npatches,dInt *off,const dInt **ind,const dReal **weight);
 extern dErr dRuleGetNodeWeight(dRule rule,dReal *coord,dReal *weight);
 extern dErr dRuleGetTensorNodeWeight(dRule rule,dInt *dim,dInt *nnodes,const dReal *coord[],const dReal *weight[]);
 extern dErr dRuleComputeGeometry(dRule rule,const dReal vtx[restrict][3],dReal[restrict][3],dReal jinv[restrict][3][3],dReal jdet[restrict]);

src/jacobi/interface/jacobi.c

@@ -290,6 +290,8 @@ dErr dRuleView(dRule rule,PetscViewer viewer)
   dFunctionReturn(0);
 }
 
+/** Get topological dimension of rule and total number of nodes
+ */
 dErr dRuleGetSize(dRule rule,dInt *dim,dInt *nnodes)
 {
   dErr err;
@@ -300,6 +302,28 @@ dErr dRuleGetSize(dRule rule,dInt *dim,dInt *nnodes)
   dFunctionReturn(0);
 }
 
+/** Decompose a composite rule into patches that preserve additional sparsity.
+ *
+ * @note Rationale: There is a one-many relationship between elements and patches.  High-order quadratures cannot be
+ * decomposed and are used during traversals that require it.  When a high-order rule has embedded low-order rules,
+ * additional sparsity can be preserved by using the low-order rules instead of the single high-order rule.  However,
+ * such low-order quadrature involves smaller inner loops and may not be the fastest way to evaluate residuals (where
+ * tensor products can be used and sparsity does not appear explicitly).
+ */
+dErr dRuleGetPatches(dRule rule,dInt *npatches,dInt *patchsize,const dInt **ind,const dReal **weights)
+{
+  dErr err;
+
+  dFunctionBegin;
+  dValidPointer(rule,1);
+  if (rule->ops.getPatches) {
+    err = (*rule->ops.getPatches)(rule,npatches,patchsize,ind,weights);dCHK(err);
+  } else dERROR(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Rule does not have patches.  Maybe trivial patches should be implemented?");
+  dFunctionReturn(0);
+}
+
+/** Get location of nodes and associated weights on reference element
+ */
 dErr dRuleGetNodeWeight(dRule rule,dReal *coord,dReal *weight)
 {
   dErr err;
@@ -337,7 +361,9 @@ dErr dRuleGetNodeWeight(dRule rule,dReal *coord,dReal *weight)
   dFunctionReturn(0);
 }
 
-dErr dRuleGetTensorNodeWeight(dRule rule,dInt *dim,dInt *nnodes,const dReal *coord[],const dReal *weight[])
+/** Get tensor-product representation of coordinate nodes and weights on reference element
+ */
+dErr dRuleGetTensorNodeWeight(dRule rule,dInt *dim,dInt nnodes[3],const dReal *coord[3],const dReal *weight[3])
 {
   dErr err;
 

src/jacobi/quadrature/impls/tensor/ruletopo.c

@@ -2,6 +2,8 @@
 #include "tensorquad.h"
 
 static dErr dRuleView_Tensor_Private(const char*,dInt,TensorRule*,PetscViewer);
+static dErr dRuleGetPatches_Tensor_All(dRule,dInt*,dInt*,const dInt **,const dReal**);
+
 #define _F(f) static dErr f(dRule,PetscViewer)
 _F(dRuleView_Tensor_Line);
 _F(dRuleView_Tensor_Quad);
@@ -231,27 +233,41 @@ static dErr dRuleComputeGeometry_Tensor_Hex(dRule rule,const dReal x[restrict][3
   dFunctionReturn(0);
 }
 
+static dErr dRuleGetPatches_Tensor_All(dRule grule,dInt *npatches,dInt *patchsize,const dInt **ind,const dReal **weight)
+{
+  dRule_Tensor *rule = (dRule_Tensor*)grule;
+  dFunctionBegin;
+  *npatches  = rule->npatches;
+  *patchsize = rule->patchsize;
+  *ind       = rule->patchind;
+  *weight    = rule->patchweight;
+  dFunctionReturn(0);
+}
+
 dErr dQuadratureRuleOpsSetUp_Tensor(dQuadrature quad)
 {
   static const struct _dRuleOps ruleOpsLine = {
     .view                = dRuleView_Tensor_Line,
     .getSize             = dRuleGetSize_Tensor_Line,
     .getNodeWeight       = NULL, /* dRuleGetNodeWeight_Tensor_Line, */
     .getTensorNodeWeight = dRuleGetTensorNodeWeight_Tensor_Line,
+    .getPatches          = dRuleGetPatches_Tensor_All,
     .computeGeometry     = NULL, /* Not implemented */
   };
   static const struct _dRuleOps ruleOpsQuad = {
     .view                = dRuleView_Tensor_Quad,
     .getSize             = dRuleGetSize_Tensor_Quad,
     .getNodeWeight       = NULL, /* dRuleGetNodeWeight_Tensor_Quad, */
     .getTensorNodeWeight = dRuleGetTensorNodeWeight_Tensor_Quad,
+    .getPatches          = dRuleGetPatches_Tensor_All,
     .computeGeometry     = NULL, /* Not implemented */
   };
   static const struct _dRuleOps ruleOpsHex  = {
     .view                = dRuleView_Tensor_Hex,
     .getSize             = dRuleGetSize_Tensor_Hex,
     .getNodeWeight       = NULL, /* dRuleGetNodeWeight_Tensor_Hex, */
     .getTensorNodeWeight = dRuleGetTensorNodeWeight_Tensor_Hex,
+    .getPatches          = dRuleGetPatches_Tensor_All,
     .computeGeometry     = dRuleComputeGeometry_Tensor_Hex,
   };
   dQuadrature_Tensor *tnsr = quad->data;

src/jacobi/quadrature/impls/tensor/tensorquad.c

@@ -68,9 +68,9 @@ static dErr TensorGetRule(dQuadrature_Tensor *tnsr,dQuadratureMethod method,dInt
       case dQUADRATURE_METHOD_SPARSE:
         /* The meaning of this option is somewhat unfortunate, but I don't see a clearly better API at the moment.  Here
          * we just assume a Legendre-Gauss-Lobatto tensor product basis, and construct a quadrature (either Gauss or
-         * Lobatto) on the subelements.  A better system would somehow have the subelement details available explicitly.
+         * Lobatto) on the patches.  A better system would somehow have the subelement details available explicitly.
          */
-        size = order;
+        size = 2*order;
         if (tnsr->alpha != 0 || tnsr->beta != 0) dERROR(PETSC_COMM_SELF,PETSC_ERR_SUP,"only alpha=0, beta=0 (Legendre)");
         switch (tnsr->family) {
           case dGAUSS_GAUSS: nodes_and_weights = two_point_gauss; break;
@@ -115,6 +115,57 @@ static dErr TensorGetRule(dQuadrature_Tensor *tnsr,dQuadratureMethod method,dInt
   dFunctionReturn(0);
 }
 
+static dErr dRulePatchSetup_Tensor(dRule_Tensor *rule,dQuadratureMethod method)
+{
+  dErr err;
+  dInt i,j,k,ii,jj,kk,dim,N[3],Q[3],P[3],NN,PP;
+  const dReal *W[3];
+
+  dFunctionBegin;
+  err = dRuleGetTensorNodeWeight((dRule)rule,&dim,Q,NULL,W);dCHK(err);
+  switch (method) {
+  case dQUADRATURE_METHOD_SPARSE: /* Tile the element with 2*2*2 quadrature point patches */
+    for (i=0; i<3; i++) {
+      P[i] = Q[i] > 1 ? 2 : 1;
+      N[i] = Q[i] / P[i];
+      if (N[i]*P[i] != Q[i]) dERROR(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Trying to use sparse quadrature, but the number of points %D is not even",Q[i]);
+    }
+    break;
+  default:                      /* Use one big patch for the whole element */
+    N[0] = N[1] = N[2] = 1;
+    err = dMemcpy(P,Q,sizeof(Q));dCHK(err);
+  }
+  rule->npatches = NN = N[0]*N[1]*N[2];          /* Number of patches in this element */
+  rule->patchsize = PP = P[0]*P[1]*P[2];          /* Number of nodes per patch */
+  err = dMallocA2(NN*PP,&rule->patchind,NN*PP,&rule->patchweight);dCHK(err);
+  /* 3 nested loops over the patches */
+  for (i=0; i<N[0]; i++) {
+    for (j=0; j<N[1]; j++) {
+      for (k=0; k<N[2]; k++) {
+        /* 3 nested loops over quadrature points on each patch */
+        for (ii=0; ii<P[0]; ii++) {
+          for (jj=0; jj<P[1]; jj++) {
+            for (kk=0; kk<P[2]; kk++) {
+              const dInt
+                patch    = (i*N[1]+j)*N[2]+k,    /* Patch number within element */
+                patchidx = (ii*P[1]+jj)*P[2]+kk, /* Index of node within current patch */
+                ielem    = i*P[0]+ii,
+                jelem    = j*P[1]+jj,
+                kelem    = k*P[2]+kk,
+                elemidx  = (ielem*Q[1] + jelem)*Q[2] + kelem; /* Index of node within current element */
+              rule->patchind[patch*PP+patchidx] = elemidx;
+              rule->patchweight[patch*PP+patchidx] = (rule->trule[0]->weight[ielem]
+                                                     * (dim < 2 ? 1 : rule->trule[1]->weight[jelem])
+                                                     * (dim < 3 ? 1 : rule->trule[2]->weight[jelem]));
+            }
+          }
+        }
+      }
+    }
+  }
+  dFunctionReturn(0);
+}
+
 static dErr dQuadratureGetRules_Tensor_Private(dQuadrature quad,dInt n,const dEntTopology topo[],const dPolynomialOrder order[],dRule rules[],dQuadratureMethod method)
 {
   dQuadrature_Tensor *tnsr = quad->data;
@@ -147,6 +198,7 @@ static dErr dQuadratureGetRules_Tensor_Private(dQuadrature quad,dInt n,const dEn
           break;
         default: dERROR(PETSC_COMM_SELF,1,"no rule available for given topology %s",dMeshEntTopologyName(topo[i]));
       }
+      err = dRulePatchSetup_Tensor(newrule,method);dCHK(err);
       kh_val(tnsr->rules,kiter) = newrule;
     }
     rules[i] = (dRule)kh_val(tnsr->rules,kiter);
@@ -323,8 +375,11 @@ static dErr dQuadratureDestroy_Tensor(dQuadrature quad)
   kh_destroy_facetrule(tnsr->facetrules);
   /* Destroy all the cached (n-dimensional) rules */
   for (k=kh_begin(tnsr->rules); k!=kh_end(tnsr->rules); k++) {
+    dRule_Tensor *rule;
     if (!kh_exist(tnsr->rules,k)) continue;
-    err = dFree(kh_val(tnsr->rules,k));dCHK(err);
+    rule = kh_val(tnsr->rules,k);
+    err = dFree2(rule->patchind,rule->patchweight);dCHK(err);
+    err = dFree(rule);dCHK(err);
   }
   kh_destroy_rule(tnsr->rules);
   /* Destroy all cached 1D tensor rules */

src/jacobi/quadrature/impls/tensor/tensorquad.h

@@ -17,6 +17,9 @@ struct s_TensorRule {
 typedef struct {
   dRuleHEADER;
   dEntTopology topo;
+  dInt npatches,patchsize;
+  dInt *patchind;
+  dReal *patchweight;
   TensorRule trule[3];
 } dRule_Tensor;
 

src/jacobi/tests/CMakeLists.txt

@@ -1,4 +1,4 @@
-foreach (i RANGE 1 3)
+foreach (i RANGE 1 4)
   dohp_link_executable("jacobi-ex${i}" "ex${i}.c")
 endforeach ()
 
@@ -11,3 +11,5 @@ dohp_add_test (jacobi-3-modal-gauss 1 jacobi-ex3 -djac_type modal -rule_degree 6
 dohp_add_test (jacobi-3-tensor-gauss 1 jacobi-ex3 -djac_type tensor -rule_degree 8 -dquad_tensor_gauss_family gauss)
 dohp_add_test (jacobi-3-tensor-lobatto 1 jacobi-ex3 -djac_type tensor -rule_degree 6 -dquad_tensor_gauss_family lobatto)
 dohp_add_test (jacobi-3-tensor-sparse-gauss 1 jacobi-ex3 -djac_type tensor -rule_degree 4 -dquad_tensor_gauss_family gauss -qmethod sparse -just_view)
+dohp_add_test (jacobi-4-tensor-sparse 1 jacobi-ex4 -quadrature_method sparse)
+dohp_add_test (jacobi-4-tensor-fast 1 jacobi-ex4 -quadrature_method fast)

src/jacobi/tests/ex4.c

@@ -0,0 +1,66 @@
+static const char help[] = "Tests modal bases for dJacobi.";
+
+#include <dohpjacobi.h>
+#include <dohpsys.h>
+#include <dohp.h>
+
+static dErr TestPatches(dJacobi jac,PetscViewer viewer)
+{
+  dEntTopology topo[] = {dTOPO_HEX,dTOPO_HEX,dTOPO_HEX,dTOPO_HEX};
+  dPolynomialOrder rdegree[4],bdegree[4];
+  dRule *rules;
+  dEFS *efs;
+  dQuadrature quad;
+  dQuadratureMethod method = dQUADRATURE_METHOD_FAST;
+  dInt bp = 5;
+  dErr err;
+
+  dFunctionBegin;
+  err = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TestPatches Options",NULL);dCHK(err);
+  err = PetscOptionsInt("-basis_degree","Degree of tensor-product basis polynomial",NULL,bp,&bp,NULL);dCHK(err);
+  err = PetscOptionsEnum("-quadrature_method","Method to use for sample quadrature",NULL,dQuadratureMethods,(PetscEnum)method,(PetscEnum*)&method,NULL);dCHK(err);
+  err = PetscOptionsEnd();dCHK(err);
+  for (dInt i=0; i<4; i++) {
+    rdegree[i] = dPolynomialOrderCreate(0,i+1,i+1,i+1);
+    bdegree[i] = dPolynomialOrderCreate(0,bp,bp,bp);
+  }
+  err = dJacobiGetQuadrature(jac,method,&quad);dCHK(err);
+  err = dQuadratureGetRules(quad,4,topo,rdegree,&rules);dCHK(err);
+  err = dJacobiGetEFS(jac,4,topo,bdegree,rules,&efs);dCHK(err);
+
+  for (dInt i=0; i<4; i++) {
+    dInt npatches,patchsize,dim,tsize[3];
+    const dInt *ind;
+    const dReal *weight,*tcoord[3],*tweight[3];
+    dRule rule;
+    err = dEFSGetRule(efs[i],&rule);dCHK(err);
+    err = dRuleGetTensorNodeWeight(rule,&dim,tsize,tcoord,tweight);dCHK(err);
+    err = dRuleGetPatches(rule,&npatches,&patchsize,&ind,&weight);dCHK(err);
+    if (npatches*patchsize != tsize[0]*tsize[1]*tsize[2]) dERROR(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent patch sizes");
+    err = PetscViewerASCIIPrintf(viewer,"Element %D: npatches %D  patchsize %D\n",i,npatches,patchsize);dCHK(err);
+    err = dIntTableView(npatches,patchsize,ind,"indices",viewer);dCHK(err);
+    err = dRealTableView(npatches,patchsize,weight,"weights",viewer);dCHK(err);
+  }
+  err = dFree(rules);dCHK(err);
+  err = dFree(efs);dCHK(err);
+  dFunctionReturn(0);
+}
+
+int main(int argc,char *argv[])
+{
+  dJacobi jac;
+  MPI_Comm comm;
+  PetscViewer viewer;
+  dErr err;
+
+  err = dInitialize(&argc,&argv,0,help);dCHK(err);
+  comm = PETSC_COMM_WORLD;
+  viewer = PETSC_VIEWER_STDOUT_WORLD;
+
+  err = dJacobiCreate(comm,&jac);dCHK(err);
+  err = dJacobiSetFromOptions(jac);dCHK(err);
+  err = TestPatches(jac,viewer);dCHK(err);
+  err = dJacobiDestroy(jac);dCHK(err);
+  err = dFinalize();dCHK(err);
+  return 0;
+}

src/jacobi/tests/refout/jacobi-4-tensor-fast.refout

@@ -0,0 +1,12 @@
+Element 0: npatches 1  patchsize 1
+   indices[ 0][ 0: 1]    0
+   weights[ 0][ 0: 1]    8.00000
+Element 1: npatches 1  patchsize 8
+   indices[ 0][ 0: 8]    0   1   2   3   4   5   6   7
+   weights[ 0][ 0: 8]    1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000
+Element 2: npatches 1  patchsize 8
+   indices[ 0][ 0: 8]    0   1   2   3   4   5   6   7
+   weights[ 0][ 0: 8]    1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000
+Element 3: npatches 1  patchsize 27
+   indices[ 0][ 0:27]    0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26
+   weights[ 0][ 0:27]    0.17147   0.17147   0.17147   0.43896   0.43896   0.43896   0.17147   0.17147   0.17147   0.27435   0.27435   0.27435   0.70233   0.70233   0.70233   0.27435   0.27435   0.27435   0.17147   0.17147   0.17147   0.43896   0.43896   0.43896   0.17147   0.17147   0.17147