Allow reuse of iterator between iterations

jedbrown · Mar 13, 2011 · 0ba4ee9 · 0ba4ee9
1 parent e3f409f
commit 0ba4ee9
Show file tree

Hide file tree

Showing 2 changed files with 42 additions and 13 deletions.
diff --git a/src/fs/interface/fsrulesetit.c b/src/fs/interface/fsrulesetit.c
@@ -23,6 +23,7 @@ static dErr ValueCachePhysicalReset(struct ValueCachePhysical *phys);
 static dErr ValueCachePhysicalDestroy(struct ValueCachePhysical *phys);
 
 struct ExplicitCache {
+  dInt maxnpatches,maxQ,maxP;
   dInt  P;                     /**< number of finite element modes on each patch */
   dInt  *bidx;                 /**< For each patch, indices of modes with support on patch. Indices are relative to expanded space. */
   dReal *interp;               /**< Interpolation: u_quad[patch*Q+i] = sum_j interp[patch*Q*P+i*P+j] * u_patch[eidx[patch*P+j]] */
@@ -761,10 +762,14 @@ static dErr ValueCacheSetUp(struct ValueCache *vc,dInt bs,dInt nmax)
   dErr err;
 
   dFunctionBegin;
-  err = dMemzero(vc,sizeof(*vc));dCHK(err);
-  vc->bs = bs;
-  vc->n_alloc = nmax;
-  err = dMallocA4(bs*nmax,&vc->u_alloc,bs*nmax*3,&vc->du_alloc,bs*nmax,&vc->v_alloc,bs*nmax*3,&vc->dv_alloc);dCHK(err);
+  if (bs != vc->bs || nmax > vc->n_alloc) {
+    err = dFree4(vc->u_alloc,vc->du_alloc,vc->v_alloc,vc->dv_alloc);dCHK(err);
+    err = dMemzero(vc,sizeof(*vc));dCHK(err);
+    vc->bs = bs;
+    vc->n_alloc = nmax;
+    err = dMallocA4(bs*nmax,&vc->u_alloc,bs*nmax*3,&vc->du_alloc,bs*nmax,&vc->v_alloc,bs*nmax*3,&vc->dv_alloc);dCHK(err);
+  }
+  err = ValueCacheReset(vc,dFALSE);dCHK(err);
   dFunctionReturn(0);
 }
 
@@ -922,8 +927,14 @@ static dErr ExplicitCacheSetUp(struct ExplicitCache *explicit,dInt maxnpatches,d
   dErr err;
 
   dFunctionBegin;
-  explicit->P = -1;
-  err = dMallocA3(maxnpatches*maxP,&explicit->bidx,maxnpatches*maxQ*maxP,&explicit->interp,maxnpatches*maxQ*maxP*3,&explicit->deriv);dCHK(err);
+  if (!(maxnpatches <= explicit->maxnpatches && maxQ <= explicit->maxQ && maxP <= explicit->maxP)) {
+    err = dFree3(explicit->bidx,explicit->interp,explicit->deriv);dCHK(err);
+    explicit->maxnpatches = maxnpatches;
+    explicit->maxQ        = maxQ;
+    explicit->maxP        = maxP;
+    explicit->P = -1;
+    err = dMallocA3(maxnpatches*maxP,&explicit->bidx,maxnpatches*maxQ*maxP,&explicit->interp,maxnpatches*maxQ*maxP*3,&explicit->deriv);dCHK(err);
+  }
   err = dMemzero(explicit->bidx,maxnpatches*maxP*sizeof(explicit->bidx[0]));dCHK(err);
   err = dMemzero(explicit->interp,maxnpatches*maxQ*maxP*sizeof(explicit->interp[0]));dCHK(err);
   err = dMemzero(explicit->deriv,maxnpatches*maxQ*maxP*3*sizeof(explicit->deriv[0]));dCHK(err);

diff --git a/src/fs/tests/ex1.c b/src/fs/tests/ex1.c
@@ -254,6 +254,7 @@ struct ProjContext {
   dFS fs;
   Vec x,y;
   dRuleset ruleset[dQUADRATURE_METHOD_INVALID];
+  dRulesetIterator regioniter[dQUADRATURE_METHOD_INVALID];
 };
 
 static dErr ProjGetRuleset(struct ProjContext *proj,dQuadratureMethod qmethod,dRuleset *ruleset)
@@ -270,6 +271,25 @@ static dErr ProjGetRuleset(struct ProjContext *proj,dQuadratureMethod qmethod,dR
   dFunctionReturn(0);
 }
 
+static dErr ProjGetRegionIterator(struct ProjContext *proj,dQuadratureMethod qmethod,dRulesetIterator *riter)
+{
+  dErr err;
+
+  dFunctionBegin;
+  if (!proj->regioniter[qmethod]) {
+    dRulesetIterator iter;
+    dRuleset ruleset;
+    dFS cfs;
+    err = ProjGetRuleset(proj,qmethod,&ruleset);dCHK(err);
+    err = dFSGetCoordinateFS(proj->fs,&cfs);dCHK(err);
+    err = dRulesetCreateIterator(ruleset,cfs,&iter);dCHK(err);
+    err = dRulesetIteratorAddFS(iter,proj->fs);dCHK(err);
+    proj->regioniter[qmethod] = iter;
+  }
+  *riter = proj->regioniter[qmethod];
+  dFunctionReturn(0);
+}
+
 static dErr ProjResidual(dUNUSED SNES snes,Vec gx,Vec gy,void *ctx)
 {
   dErr err;
@@ -385,18 +405,14 @@ static dErr ProjResidual3(dUNUSED SNES snes,Vec gx,Vec gy,void *ctx)
 {
   dErr err;
   struct ProjContext *proj = ctx;
-  dFS cfs,fs = proj->fs;
-  dRuleset ruleset;
+  dFS fs = proj->fs;
   Vec Coords;
   dRulesetIterator iter;
 
   dFunctionBegin;
-  err = ProjGetRuleset(proj,gopt.proj_qmethod,&ruleset);dCHK(err);
   err = VecZeroEntries(gy);dCHK(err);
-  err = dFSGetCoordinateFS(fs,&cfs);dCHK(err);
   err = dFSGetGeometryVectorExpanded(fs,&Coords);dCHK(err);
-  err = dRulesetCreateIterator(ruleset,cfs,&iter);dCHK(err);
-  err = dRulesetIteratorAddFS(iter,fs);dCHK(err);
+  err = ProjGetRegionIterator(proj,gopt.proj_qmethod,&iter);dCHK(err);
   err = dRulesetIteratorStart(iter,Coords,NULL,gx,gy);dCHK(err);
 
   while (dRulesetIteratorHasPatch(iter)) {
@@ -417,7 +433,6 @@ static dErr ProjResidual3(dUNUSED SNES snes,Vec gx,Vec gy,void *ctx)
   }
 
   err = dRulesetIteratorFinish(iter);dCHK(err);
-  err = dRulesetIteratorDestroy(iter);dCHK(err);
   dFunctionReturn(0);
 }
 
@@ -660,6 +675,9 @@ static dErr doProjection(dFS fs)
   for (dInt i=0; i<ALEN(proj.ruleset); i++) {
     if (proj.ruleset[i]) {err = dRulesetDestroy(proj.ruleset[i]);dCHK(err);}
   }
+  for (dInt i=0; i<ALEN(proj.regioniter); i++) {
+    if (proj.regioniter[i]) {err = dRulesetIteratorDestroy(proj.regioniter[i]);dCHK(err);}
+  }
   dFunctionReturn(0);
 }