Commits

Nathan Collier committed f45d44b

added support for precomputation of 1d basis functions for collocation

Comments (0)

Files changed (3)

 
 typedef struct _n_IGAElement  *IGAElement;
 typedef struct _n_IGAPoint    *IGAPoint;
+
 typedef struct _n_IGAColPoint *IGAColPoint;
+typedef struct _n_IGAColBasis *IGAColBasis;
 
 /* ---------------------------------------------------------------- */
 
   IGABasis basis[3];
   IGABoundary boundary[3][2];
   IGAElement  iterator;
-  IGAColPoint point_iterator;
+
+  IGAColPoint point_iterator; /* stuff added for collocation */
+  IGAColBasis colbasis[3];
 
   PetscInt  proc_sizes[3];
   PetscInt  proc_ranks[3];
 PETSC_EXTERN PetscErrorCode IGAColFormSystem(IGA iga,Mat matA,Vec vecB,IGAColUserSystem System,void *ctx);
 PETSC_EXTERN PetscErrorCode IGAColComputeSystem(IGA iga,Mat A,Vec B);
 
+struct _n_IGAColBasis {
+  PetscInt refct;
+  /**/
+  PetscInt  ncp;      /* number of collocation points */
+  PetscInt  nen;      /* number of local basis functions */
+  PetscInt  p,d;      /* polynomial order, last derivative index */
+
+  PetscInt  *offset;  /* [ncp] basis offset   */
+  PetscReal *detJ;    /* [ncp]                */
+  PetscReal *point;   /* [ncp]                */
+  PetscReal *value;   /* [ncp][nen][d+1]      */
+};
+PETSC_EXTERN PetscErrorCode IGAColBasisCreate(IGAColBasis *basis);
+PETSC_EXTERN PetscErrorCode IGAColBasisDestroy(IGAColBasis *basis);
+PETSC_EXTERN PetscErrorCode IGAColBasisReset(IGAColBasis basis);
+PETSC_EXTERN PetscErrorCode IGAColBasisReference(IGAColBasis basis);
+PETSC_EXTERN PetscErrorCode IGAColBasisInit(IGAColBasis basis,IGAAxis axis,PetscInt d);
+
 /* ---------------------------------------------------------------- */
 
 #ifndef PetscMalloc1
     ierr = IGAAxisCreate(&iga->axis[i]);CHKERRQ(ierr);
     ierr = IGARuleCreate(&iga->rule[i]);CHKERRQ(ierr);
     ierr = IGABasisCreate(&iga->basis[i]);CHKERRQ(ierr);
+    ierr = IGAColBasisCreate(&iga->colbasis[i]);CHKERRQ(ierr); /* collocation */
     ierr = IGABoundaryCreate(&iga->boundary[i][0]);CHKERRQ(ierr);
     ierr = IGABoundaryCreate(&iga->boundary[i][1]);CHKERRQ(ierr);
     iga->proc_sizes[i] = -1;
     ierr = IGAAxisDestroy(&iga->axis[i]);CHKERRQ(ierr);
     ierr = IGARuleDestroy(&iga->rule[i]);CHKERRQ(ierr);
     ierr = IGABasisDestroy(&iga->basis[i]);CHKERRQ(ierr);
+    ierr = IGAColBasisDestroy(&iga->colbasis[i]);CHKERRQ(ierr); /* collocation */
     ierr = IGABoundaryDestroy(&iga->boundary[i][0]);CHKERRQ(ierr);
     ierr = IGABoundaryDestroy(&iga->boundary[i][1]);CHKERRQ(ierr);
   }
       ierr = IGARuleInit(iga->rule[i],q);CHKERRQ(ierr);
     }
     ierr = IGABasisInit(iga->basis[i],iga->axis[i],iga->rule[i],iga->order);CHKERRQ(ierr);
+    ierr = IGAColBasisInit(iga->colbasis[i],iga->axis[i],iga->order);CHKERRQ(ierr);
   }
   ierr = IGAElementInit(iga->iterator,iga);CHKERRQ(ierr);
   ierr = IGAColPointInit(iga->point_iterator,iga);CHKERRQ(ierr);
 
 #undef  __FUNCT__
 #define __FUNCT__ "_FindSpan"
-PetscInt _FindSpan(PetscInt n,PetscInt p,PetscReal u, PetscReal *U)
+PetscInt _FindSpan(PetscInt n,PetscInt p,PetscReal u, const PetscReal *U)
 {
   /* n is the index of the last basis */
   PetscInt low,high,span;
 
 #undef  __FUNCT__
 #define __FUNCT__ "_Greville"
-PetscReal _Greville(PetscInt i,PetscInt p,PetscInt m,PetscReal *U)
+PetscReal _Greville(PetscInt i,PetscInt p,PetscInt m,const PetscReal *U)
 {
   PetscInt j;
   PetscReal X = 0.0;
 #define __FUNCT__ "IGAColPointNext"
 PetscBool IGAColPointNext(IGAColPoint point)
 {
-  IGA      iga = point->parent;
   PetscInt i,dim  = point->dim;
   PetscInt *start = point->start;
   PetscInt *width = point->width;
   PetscInt *ID    = point->ID;
-  PetscInt *span  = point->span;
-  PetscReal *pnt  = point->point;
   PetscInt index,coord;
- 
+  IGAColBasis *BD = point->parent->colbasis;
+
   point->nvec = 0;
   point->nmat = 0;
 
     coord   = index % width[i];
     index   = (index - coord) / width[i];
     ID[i]   = coord + start[i]; 
-    pnt[i]  = _Greville(ID[i],iga->axis[i]->p,iga->axis[i]->m,iga->axis[i]->U);
-    span[i] = _FindSpan(iga->axis[i]->m-iga->axis[i]->p-1,iga->axis[i]->p,pnt[i],iga->axis[i]->U);
+    point->point[i] = BD[i]->point[ID[i]];
   }
-  for (i=dim; i<3; i++) {
-    span[i] = 0;
-  }
-  //printf("ID: {%d,%d,%d}  pnt: {%.2f,%.2f,%.2f}  span: {%d,%d,%d}\n",ID[0],ID[1],ID[2],pnt[0],pnt[1],pnt[2],span[0],span[1],span[2]);
   IGAColPointBuildMapping(point);
   IGAColPointBuildGeometry(point);
   IGAColPointBuildShapeFuns(point);
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during collocation point loop");
   { /* */
     IGA      iga = point->parent;
-    IGABasis *BD = iga->basis;
-    PetscInt *span = point->span;
-    PetscInt ia, inen = BD[0]->nen, ioffset = span[0]-iga->axis[0]->p;
-    PetscInt ja, jnen = BD[1]->nen, joffset = span[1]-iga->axis[1]->p;
-    PetscInt ka, knen = BD[2]->nen, koffset = span[2]-iga->axis[2]->p;
+    IGAColBasis *BD = iga->colbasis;
+    PetscInt *ID = point->ID;
+    
+    PetscInt ia, inen = BD[0]->nen, ioffset = BD[0]->offset[ID[0]];
+    PetscInt ja, jnen = BD[1]->nen, joffset = BD[1]->offset[ID[1]];
+    PetscInt ka, knen = BD[2]->nen, koffset = BD[2]->offset[ID[2]];
     PetscInt *start = iga->node_gstart, *width = iga->node_gwidth;
     PetscInt istart = start[0]/*istride = 1*/;
     PetscInt jstart = start[1], jstride = width[0];
 extern void IGA_DersBasisFuns(PetscInt i,PetscReal u,PetscInt p,PetscInt d,const PetscReal U[],PetscReal N[]);
 EXTERN_C_END
 
-#define IGA_BasisFuns_ARGS(BD,i) 1,BD[i]->nen,BD[i]->d,point->basis1d[i]
+#define IGA_BasisFuns_ARGS(ID,BD,i) 1,BD[i]->nen,BD[i]->d,BD[i]->value+ID[i]*BD[i]->nen*(BD[i]->d+1)
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGAColPointBuildShapeFuns"
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during collocation point loop");
   order = point->parent->order;
   {
-    IGABasis *BD  = point->parent->basis;
+    IGAColBasis *BD  = point->parent->colbasis;
+    PetscInt *ID  = point->ID;
     PetscReal **N = point->basis;
     switch (point->dim) {
-    case 3: 
-      IGA_DersBasisFuns(point->span[0],point->point[0],point->parent->axis[0]->p,point->parent->order,point->parent->axis[0]->U,point->basis1d[0]);
-      IGA_DersBasisFuns(point->span[1],point->point[1],point->parent->axis[1]->p,point->parent->order,point->parent->axis[1]->U,point->basis1d[1]);
-      IGA_DersBasisFuns(point->span[2],point->point[2],point->parent->axis[2]->p,point->parent->order,point->parent->axis[2]->U,point->basis1d[2]);
-      IGA_BasisFuns_3D(order,point->rational,
-		       point->geometryW,
-		       IGA_BasisFuns_ARGS(BD,0),
-		       IGA_BasisFuns_ARGS(BD,1),
-		       IGA_BasisFuns_ARGS(BD,2),
-		       N[0],N[1],N[2],N[3]); break;
-    case 2: 
-      IGA_DersBasisFuns(point->span[0],point->point[0],point->parent->axis[0]->p,point->parent->order,point->parent->axis[0]->U,point->basis1d[0]);
-      IGA_DersBasisFuns(point->span[1],point->point[1],point->parent->axis[1]->p,point->parent->order,point->parent->axis[1]->U,point->basis1d[1]);
-      IGA_BasisFuns_2D(order,point->rational,
-		       point->geometryW,
-		       IGA_BasisFuns_ARGS(BD,0),
-		       IGA_BasisFuns_ARGS(BD,1),
-		       N[0],N[1],N[2],N[3]); break;
-    case 1: 
-      IGA_DersBasisFuns(point->span[0],point->point[0],point->parent->axis[0]->p,point->parent->order,point->parent->axis[0]->U,point->basis1d[0]);
-      IGA_BasisFuns_1D(order,point->rational,
-		       point->geometryW,
-		       IGA_BasisFuns_ARGS(BD,0),
-		       N[0],N[1],N[2],N[3]); break;
+    case 3: IGA_BasisFuns_3D(order,point->rational,
+                             point->geometryW,
+                             IGA_BasisFuns_ARGS(ID,BD,0),
+                             IGA_BasisFuns_ARGS(ID,BD,1),
+                             IGA_BasisFuns_ARGS(ID,BD,2),
+                             N[0],N[1],N[2],N[3]); break;
+    case 2: IGA_BasisFuns_2D(order,point->rational,
+                             point->geometryW,
+                             IGA_BasisFuns_ARGS(ID,BD,0),
+                             IGA_BasisFuns_ARGS(ID,BD,1),
+                             N[0],N[1],N[2],N[3]); break;
+    case 1: IGA_BasisFuns_1D(order,point->rational,
+                             point->geometryW,
+                             IGA_BasisFuns_ARGS(ID,BD,0),
+                             N[0],N[1],N[2],N[3]); break;
     }
   }
   if (point->dim == point->nsd) /* XXX */
   PetscFunctionReturn(0);
 }
 
+/* The following routines are a slight alteration of the petigaksp.c routines */
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAColComputeSystem"
 PetscErrorCode IGAColComputeSystem(IGA iga,Mat matA,Vec vecB)
   if (SysCtx) iga->userops->ColSysCtx = SysCtx;
   PetscFunctionReturn(0);
 }
+
+/* The following routines are a slight alteration of the petigabasis.c routines */
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColBasisCreate"
+PetscErrorCode IGAColBasisCreate(IGAColBasis *basis)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(basis,3);
+  ierr = PetscNew(struct _n_IGAColBasis,basis);CHKERRQ(ierr);
+  (*basis)->refct = 1;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColBasisDestroy"
+PetscErrorCode IGAColBasisDestroy(IGAColBasis *_basis)
+{
+  IGAColBasis       basis;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(_basis,1);
+  basis = *_basis; *_basis = 0;
+  if (!basis) PetscFunctionReturn(0);
+  if (--basis->refct > 0) PetscFunctionReturn(0);
+  ierr = IGAColBasisReset(basis);CHKERRQ(ierr);
+  ierr = PetscFree(basis);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColBasisReset"
+PetscErrorCode IGAColBasisReset(IGAColBasis basis)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  if (!basis) PetscFunctionReturn(0);
+  PetscValidPointer(basis,1);
+  basis->ncp = 0;
+  basis->nen = 0;
+  basis->p   = 0;
+  basis->d   = 0;
+  ierr = PetscFree(basis->offset);CHKERRQ(ierr);
+  ierr = PetscFree(basis->detJ);CHKERRQ(ierr);
+  ierr = PetscFree(basis->point);CHKERRQ(ierr);
+  ierr = PetscFree(basis->value);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColBasisReference"
+PetscErrorCode IGAColBasisReference(IGAColBasis basis)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(basis,1);
+  basis->refct++;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAColBasisInit"
+PetscErrorCode IGAColBasisInit(IGAColBasis basis,IGAAxis axis,PetscInt d)
+{
+  PetscInt       m,p;
+  const PetscReal*U;
+  PetscInt       icp,ncp;
+  PetscInt       nen,ndr;
+  PetscInt       *offset;
+  PetscReal      *detJ;
+  PetscReal      *point;
+  PetscReal      *value;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(basis,1);
+  PetscValidPointer(axis,2);
+  if (d < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
+                      "Derivative order must be grather than zero, got %D",d);
+
+  m = axis->m;
+  p = axis->p;
+  U = axis->U;
+
+  ncp  = axis->nnp;
+  nen  = p+1;
+  ndr  = d+1;
+
+  ierr = PetscMalloc1(ncp,PetscInt,&offset);CHKERRQ(ierr);
+  ierr = PetscMalloc1(ncp,PetscReal,&detJ);CHKERRQ(ierr);
+  ierr = PetscMalloc1(ncp,PetscReal,&point);CHKERRQ(ierr);
+  ierr = PetscMalloc1(ncp*nen*ndr,PetscReal,&value);CHKERRQ(ierr);
+  
+  for (icp=0; icp<ncp; icp++) {
+    point[icp]   = _Greville(icp,p,m,U); /* XXX Later fix to be a choice */
+    PetscInt k   = _FindSpan(ncp-1,p,point[icp],U);
+    IGA_DersBasisFuns(k,point[icp],p,d,U,&value[icp*nen*ndr]);
+    offset[icp] = k-p;
+    detJ[icp]   = 0.5*(U[k+1]-U[k]);
+  }
+
+  ierr = IGAColBasisReset(basis);CHKERRQ(ierr);
+
+  basis->ncp    = ncp;
+  basis->nen    = nen;
+  basis->p      = p;
+  basis->d      = d;
+  basis->offset = offset;
+
+  basis->detJ   = detJ;
+  basis->point  = point;
+  basis->value  = value;
+
+  PetscFunctionReturn(0);
+}
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.