Commits

Nathan Collier committed 4679e64

added interface for boundary integrals

Comments (0)

Files changed (3)

   IGAUserSystem     System;
   void              *SysCtx;  
   /**/
+  IGAUserSystem     BndSystem[3][2]; 
+  void              *BndSysCtx[3][2];  
+  /**/
   IGAUserFunction   Function;
   void              *FunCtx;
   IGAUserJacobian   Jacobian;
 PETSC_EXTERN PetscErrorCode IGAElementAssembleVec(IGAElement element,const PetscScalar F[],Vec vec);
 PETSC_EXTERN PetscErrorCode IGAElementAssembleMat(IGAElement element,const PetscScalar K[],Mat mat);
 
+PETSC_EXTERN PetscErrorCode IGABeginBoundaryElement(IGA iga,PetscInt dir, PetscInt side,IGAElement *element);
+PETSC_EXTERN PetscBool      IGANextBoundaryElement(IGA iga,PetscInt dir, PetscInt side,IGAElement element);
+PETSC_EXTERN PetscErrorCode IGABoundaryElementBeginPoint(IGAElement element,PetscInt dir, PetscInt side,IGAPoint *point);
+PETSC_EXTERN PetscBool IGABoundaryElementNextPoint(IGAElement element,PetscInt dir, PetscInt side,IGAPoint point);
+
+
 /* ---------------------------------------------------------------- */
 
 struct _n_IGAPoint {
                                               Mat J,
                                               IGAUserIEJacobian,void *);
 
+PETSC_EXTERN PetscErrorCode IGABoundaryFormSystem(IGA iga,Mat matA,Vec vecB,
+						  IGAUserSystem System,
+						  PetscInt dir, PetscInt side, void *ctx);
+
 /* ---------------------------------------------------------------- */
 
 #ifndef PetscMalloc1
   }
   PetscFunctionReturn(0);
 }
+
+
+
+
+
+
+
+
+
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABeginBoundaryElement"
+PetscErrorCode IGABeginBoundaryElement(IGA iga,PetscInt dir, PetscInt side,IGAElement *element)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidPointer(element,2);
+  IGACheckSetUp(iga,1);
+  ierr = IGAGetElement(iga,element);CHKERRQ(ierr);
+  (*element)->index = -1;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGANextBoundaryElement"
+PetscBool IGANextBoundaryElement(IGA iga,PetscInt dir, PetscInt side,IGAElement element)
+{
+  PetscInt i,dim  = element->dim;
+  PetscInt *start = element->start;
+  PetscInt *width = element->width;
+  PetscInt *ID    = element->ID;
+  PetscInt index,coord;
+  PetscErrorCode ierr;
+  /* */
+  element->nval = 0;
+  element->nvec = 0;
+  element->nmat = 0;
+  /* */
+  index = ++element->index;
+  if (PetscUnlikely(index >= element->count)) {
+    element->index = -1;
+    return PETSC_FALSE;
+  }
+  /* Loops through boundary elements on my partition */
+  for (i=0; i<dim; i++) {
+    coord = index % width[i];
+    index = (index - coord) / width[i];
+    ID[i] = coord + start[i];
+  }
+#undef  CHKERRRETURN
+#define CHKERRRETURN(n,r) do{if(PetscUnlikely(n)){CHKERRCONTINUE(n);return(r);}}while(0)
+  ierr = IGAElementBuildMapping(element);  CHKERRRETURN(ierr,PETSC_FALSE);
+  ierr = IGAElementBuildGeometry(element); CHKERRRETURN(ierr,PETSC_FALSE);
+  ierr = IGAElementBuildFix(element);      CHKERRRETURN(ierr,PETSC_FALSE);
+#undef  CHKERRRETURN
+  return PETSC_TRUE;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundaryElementBeginPoint"
+PetscErrorCode IGABoundaryElementBeginPoint(IGAElement element,PetscInt dir, PetscInt side,IGAPoint *point)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(element,1);
+  PetscValidPointer(point,2);
+  if (PetscUnlikely(element->index < 0))
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
+  ierr = IGAElementBuildQuadrature(element);CHKERRQ(ierr);
+  ierr = IGAElementBuildShapeFuns(element);CHKERRQ(ierr);
+  *point = element->iterator;
+
+  (*point)->count = element->nqp;
+  (*point)->index = -1;
+
+  (*point)->neq = element->neq;
+  (*point)->nen = element->nen;
+  (*point)->dof = element->dof;
+  (*point)->dim = element->dim;
+  (*point)->nsd = element->nsd;
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundaryElementNextPoint"
+PetscBool IGABoundaryElementNextPoint(IGAElement element,PetscInt dir, PetscInt side,IGAPoint point)
+{
+  PetscInt nen = point->nen;
+  PetscInt dim = point->dim;
+  PetscInt nsd = point->nsd;
+  PetscInt index;
+  /*PetscValidPointer(element,1);*/
+  /*PetscValidPointer(point,2);*/
+  /* */
+  point->nvec = 0;
+  point->nmat = 0;
+  /* */
+  index = ++point->index;
+  if (PetscUnlikely(index == 0))            goto start;
+  if (PetscUnlikely(index >= point->count)) goto stop;
+
+  point->weight   += 1;
+  point->detJac   += 1;
+
+  point->point    += dim;
+  point->scale    += dim;
+  point->basis[0] += nen;
+  point->basis[1] += nen*dim;
+  point->basis[2] += nen*dim*dim;
+  point->basis[3] += nen*dim*dim*dim;
+
+  point->detX     += 1;
+  point->gradX[0] += dim*dim;
+  point->gradX[1] += dim*dim;
+  point->shape[0] += nen;
+  point->shape[1] += nen*dim;
+  point->shape[2] += nen*dim*dim;
+  point->shape[3] += nen*dim*dim*dim;
+
+  return PETSC_TRUE;
+
+ start:
+
+  point->weight   = element->weight;
+  point->detJac   = element->detJac;
+
+  point->point    = element->point;
+  point->scale    = element->scale;
+  point->basis[0] = element->basis[0];
+  point->basis[1] = element->basis[1];
+  point->basis[2] = element->basis[2];
+  point->basis[3] = element->basis[3];
+
+  if (element->geometry)
+    point->geometry = element->geometryX;
+  else
+    point->geometry = PETSC_NULL;
+  if (element->geometry && dim == nsd) { /* XXX */
+    point->detX     = element->detX;
+    point->gradX[0] = element->gradX[0];
+    point->gradX[1] = element->gradX[1];
+    point->shape[0] = element->shape[0];
+    point->shape[1] = element->shape[1];
+    point->shape[2] = element->shape[2];
+    point->shape[3] = element->shape[3];
+  } else {
+    point->shape[0] = element->basis[0];
+    point->shape[1] = element->basis[1];
+    point->shape[2] = element->basis[2];
+    point->shape[3] = element->basis[3];
+  }
+
+  return PETSC_TRUE;
+
+ stop:
+
+  point->index = -1;
+  return PETSC_FALSE;
+
+}
 PetscErrorCode IGAComputeSystem(IGA iga,Mat matA,Vec vecB)
 {
   IGAUserSystem  System;
-  void           *SysCtx;
+  IGAUserSystem  BndSystem;
+  void           *SysCtx,*BndSysCtx;
   PetscErrorCode ierr;
+  PetscInt       dir,side;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidHeaderSpecific(matA,MAT_CLASSID,2);
   System = iga->userops->System;
   SysCtx = iga->userops->SysCtx;
   ierr = IGAFormSystem(iga,matA,vecB,System,SysCtx);CHKERRQ(ierr);
+  for(dir=0;dir<iga->dim;dir++){
+    for(side=0;side<2;side++){
+      BndSystem = iga->userops->BndSystem[dir][side];
+      BndSysCtx = iga->userops->BndSysCtx[dir][side];
+      if(BndSystem){
+	ierr = IGABoundaryFormSystem(iga,matA,vecB,BndSystem,dir,side,BndSysCtx);CHKERRQ(ierr);
+      }
+    }
+  }      
   PetscFunctionReturn(0);
 }
 
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAFormSystem"
+PetscErrorCode IGABoundaryFormSystem(IGA iga,Mat matA,Vec vecB,
+				     IGAUserSystem System,
+				     PetscInt dir, PetscInt side, void *ctx)
+{
+  IGAElement     element;
+  IGAPoint       point;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(matA,MAT_CLASSID,2);
+  PetscValidHeaderSpecific(vecB,VEC_CLASSID,3);
+  IGACheckSetUp(iga,1);
+
+  /* Element loop */
+  ierr = PetscLogEventBegin(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
+  ierr = IGABeginBoundaryElement(iga,dir,side,&element);CHKERRQ(ierr);
+  while (IGANextBoundaryElement(iga,dir,side,element)) {
+    PetscScalar *A, *B;
+    ierr = IGAElementGetWorkMat(element,&A);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&B);CHKERRQ(ierr);
+    /* Quadrature loop */
+    ierr = IGABoundaryElementBeginPoint(element,dir,side,&point);CHKERRQ(ierr);
+    while (IGABoundaryElementNextPoint(element,dir,side,point)) {
+      PetscScalar *K, *F;
+      ierr = IGAPointGetWorkMat(point,&K);CHKERRQ(ierr);
+      ierr = IGAPointGetWorkVec(point,&F);CHKERRQ(ierr);
+      ierr = System(point,K,F,ctx);CHKERRQ(ierr);
+      ierr = IGAPointAddMat(point,K,A);CHKERRQ(ierr);
+      ierr = IGAPointAddVec(point,F,B);CHKERRQ(ierr);
+    }
+    ierr = IGAElementEndPoint(element,&point);CHKERRQ(ierr);
+    /* */
+    ierr = IGAElementFixSystem(element,A,B);CHKERRQ(ierr);
+    ierr = IGAElementAssembleMat(element,A,matA);CHKERRQ(ierr);
+    ierr = IGAElementAssembleVec(element,B,vecB);CHKERRQ(ierr);
+  }
+  ierr = IGAEndElement(iga,&element);CHKERRQ(ierr);
+  ierr = PetscLogEventEnd(IGA_FormSystem,iga,matA,vecB,0);CHKERRQ(ierr);
+
+  ierr = MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd  (matA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = VecAssemblyBegin(vecB);CHKERRQ(ierr);
+  ierr = VecAssemblyEnd  (vecB);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGACreateKSP"
 /*@
    IGACreateKSP - Creates a KSP (linear solver) which uses the same