Commits

Lisandro Dalcin  committed 1f98deb

Improve handling of BCs for collocation

  • Participants
  • Parent commits f9e155e

Comments (0)

Files changed (4)

File demo/BoundaryIntegral.c

 #include "petiga.h"
 
+typedef struct {
+  PetscInt dir;
+  PetscInt side;
+} AppCtx;
+
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
+
+PETSC_STATIC_INLINE
+PetscReal DEL2(PetscInt dim,const PetscReal a[dim][dim])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i][i];
+  return s;
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "System"
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
-
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,1,&N1);
-
-  PetscInt a,b,i;
+  PetscInt dim = p->dim;
+  PetscInt nen = p->nen;
+  const PetscReal (*N1)[dim]; //*((PetscReal**)&N1) = p->shape[1];
+  IGAPointGetShapeFuns(p,1,(const PetscReal **)&N1);
+  PetscInt a,b;
   for (a=0; a<nen; a++) {
-    for (b=0; b<nen; b++) {
-      PetscScalar Kab = 0.0;
-      for (i=0; i<dim; i++)
-        Kab += N1[a*dim+i]*N1[b*dim+i];
-      K[a*nen+b] = Kab;
-    }
+    for (b=0; b<nen; b++)
+      K[a*nen+b] = DOT(dim,N1[a],N1[b]);
+    F[a] = 0.0;
   }
   return 0;
 }
 
-
-
 #undef  __FUNCT__
 #define __FUNCT__ "Neumann"
 PetscErrorCode Neumann(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscReal *N0 = p->shape[0];
-  PetscInt a,nen=p->nen;
-  for (a=0; a<nen; a++) {
-    PetscReal Na   = N0[a];
-    F[a] = Na * 1.0;
-  }
+  PetscInt nen = p->nen;
+  const PetscReal *N0 = p->shape[0];
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    F[a] = N0[a] * 1.0;
   return 0;
 }
 
-typedef struct { 
-  PetscInt dir;
-  PetscInt side;
-} AppCtx;
+#undef  __FUNCT__
+#define __FUNCT__ "SystemCollocation"
+PetscErrorCode SystemCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+
+  const PetscReal (*N2)[dim][dim]; //*((PetscReal**)&N2) = p->shape[2];
+  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
+
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] = -DEL2(dim,N2[a]);
+  F[0] = 0.0;
+
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "DirichletCollocation"
+PetscErrorCode DirichletCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt nen = p->nen;
+  const PetscReal *N0 = p->shape[0];
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] = N0[a];
+  F[0] = 1.0;
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "NeumannCollocation"
+PetscErrorCode NeumannCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+  const PetscReal (*N1)[dim];
+  const PetscReal *normal = p->normal;
+  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] = DOT(dim,N1[a],normal);
+  F[0] = 1.0;
+  return 0;
+}
+PetscErrorCode Neumann0Collocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{(void)NeumannCollocation(p,K,F,ctx); F[0] = 0.0; return 0;}
 
 #undef  __FUNCT__
 #define __FUNCT__ "Error"
   IGAPointFormValue(p,U,&u);
   PetscReal x;
   if (user->side == 0)
-    x = 1 - p->point[user->dir];
+    x = 1 - p->point[user->dir] + 1;
   else
-    x = p->point[user->dir];
+    x = p->point[user->dir] + 1;
   PetscReal e = u - x;
   S[0] = e*e;
   return 0;
   user.side = 1;
 
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-dir", "direction",__FILE__,user.dir, &user.dir, PETSC_NULL);CHKERRQ(ierr); 
-  ierr = PetscOptionsInt("-side","side",     __FILE__,user.side,&user.side,PETSC_NULL);CHKERRQ(ierr); 
+  ierr = PetscOptionsInt("-dir", "direction",__FILE__,user.dir, &user.dir, PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-side","side",     __FILE__,user.side,&user.side,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  if (iga->dim < 1) {ierr = IGASetDim(iga,3);CHKERRQ(ierr);}
+
+  PetscInt dim;
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
 
   IGABoundary bnd;
-  PetscInt d = !user.side; 
+  PetscInt d =  !user.side;
   PetscInt n = !!user.side;
-  ierr = IGAGetBoundary(iga,user.dir,d,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-  ierr = IGAGetBoundary(iga,user.dir,n,&bnd);CHKERRQ(ierr); 
-  ierr = IGABoundarySetUserSystem(bnd,Neumann,PETSC_NULL);CHKERRQ(ierr);
-
-  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  if (!iga->collocation) {
+    ierr = IGAGetBoundary(iga,user.dir,d,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetValue(bnd,0,1.0);CHKERRQ(ierr);
+    ierr = IGAGetBoundary(iga,user.dir,n,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetUserSystem(bnd,Neumann,PETSC_NULL);CHKERRQ(ierr);
+  } else {
+    PetscInt dir,side;
+    for (dir=0; dir<dim; dir++) {
+      for (side=0; side<2; side++) {
+        ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+        ierr = IGABoundarySetUserSystem(bnd,Neumann0Collocation,PETSC_NULL);CHKERRQ(ierr);
+      }
+    }
+    ierr = IGAGetBoundary(iga,user.dir,d,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetUserSystem(bnd,DirichletCollocation,PETSC_NULL);CHKERRQ(ierr);
+    ierr = IGAGetBoundary(iga,user.dir,n,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetUserSystem(bnd,NeumannCollocation,PETSC_NULL);CHKERRQ(ierr);
+  }
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   Mat A;
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,PETSC_NULL);CHKERRQ(ierr);
+  if (!iga->collocation) {
+    ierr = IGASetUserSystem(iga,System,PETSC_NULL);CHKERRQ(ierr);
+  } else {
+    ierr = IGASetUserSystem(iga,SystemCollocation,PETSC_NULL);CHKERRQ(ierr);
+  }
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
-  
+
   KSP ksp;
   ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
   ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
-  PetscInt dim;
-  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
-  if (dim <= 2) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
-
   PetscScalar error = 0;
   ierr = IGAFormScalar(iga,x,1,&error,Error,&user);CHKERRQ(ierr);
   error = PetscSqrtReal(PetscRealPart(error));
   ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);
-  
+
+  PetscBool draw = PETSC_TRUE;
+  if (draw && dim <= 2) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
+
   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&x);CHKERRQ(ierr);
   ierr = VecDestroy(&b);CHKERRQ(ierr);
   ierr = IGADestroy(&iga);CHKERRQ(ierr);
 
-  PetscBool flag = PETSC_FALSE;
-  PetscReal secs = -1;
-  ierr = PetscOptionsHasName(0,"-sleep",&flag);CHKERRQ(ierr);
-  ierr = PetscOptionsGetReal(0,"-sleep",&secs,0);CHKERRQ(ierr);
-  if (flag) {ierr = PetscSleep(secs);CHKERRQ(ierr);}
-
   ierr = PetscFinalize();CHKERRQ(ierr);
   return 0;
 }

File demo/Laplace.c

 
 */
 
+#include "petiga.h"
 
-#include "petiga.h"
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
+
+PETSC_STATIC_INLINE
+PetscReal DEL2(PetscInt dim,const PetscReal a[dim][dim])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i][i];
+  return s;
+}
 
 #undef  __FUNCT__
 #define __FUNCT__ "SystemGalerkin"
   PetscInt nen = p->nen;
   PetscInt dim = p->dim;
 
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,1,&N1);
+  const PetscReal (*N1)[dim];
+  IGAPointGetShapeFuns(p,1,(const PetscReal **)&N1);
 
-  PetscInt a,b,i;
+  PetscInt a,b;
   for (a=0; a<nen; a++) {
-    for (b=0; b<nen; b++) {
-      PetscScalar Kab = 0.0;
-      for (i=0; i<dim; i++)
-        Kab += N1[a*dim+i]*N1[b*dim+i];
-      K[a*nen+b] = Kab;
-    }
+    for (b=0; b<nen; b++)
+      K[a*nen+b] = DOT(dim,N1[a],N1[b]);
+    F[a] = 0.0;
   }
   return 0;
 }
   const PetscReal (*N2)[dim][dim];
   IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
 
-  PetscInt a,i;
+  PetscInt a;
   for (a=0; a<nen; a++)
-    for (i=0; i<dim; i++)
-      K[a] += -N2[a][i][i];
-
+    K[a] += -DEL2(dim,N2[a]);
   F[0] = 0.0;
 
   return 0;

File include/petiga.h

   PetscInt  sizes[3];
   PetscInt  ID[3];
   PetscBool atboundary;
-  PetscInt  atboundary_id;
+  PetscInt  boundary_id;
   /**/
   PetscInt count;
   PetscInt index;

File src/petigaelem.c

   element->nval = 0;
   element->nvec = 0;
   element->nmat = 0;
+  element->boundary_id = -1;
   /* */
   index = ++element->index;
   if (PetscUnlikely(index >= element->count)) {
 #define __FUNCT__ "IGAElementNextUserOps"
 PetscBool IGAElementNextUserOps(IGAElement element,IGAUserOps *userops)
 {
-  static PetscInt counter = 0; /* XXX */
-  IGA iga = element->parent;
+  IGA      iga = element->parent;
   PetscInt dim = element->dim;
-  for (; PetscUnlikely(counter < 2*dim); counter++) {
-    PetscInt dir  = counter / 2;
-    PetscInt side = counter % 2;
-    PetscInt iel  = side ? element->sizes[dir]-1 : 0;
-    if (element->ID[dir] != iel) continue;
+  while (++element->boundary_id < 2*dim) {
+    PetscInt i = element->boundary_id / 2;
+    PetscInt s = element->boundary_id % 2;
+    PetscInt e = s ? element->sizes[i]-1 : 0;
+    if (element->ID[i] != e) continue;
+    if (!iga->boundary[i][s]->userops) continue;
+    *userops = iga->boundary[i][s]->userops;
     element->atboundary = PETSC_TRUE;
-    element->atboundary_id = counter;
-    *userops = iga->boundary[dir][side]->userops;
-    if (!(*userops)) continue;
-    counter++;
     return PETSC_TRUE;
   }
-  if (PetscLikely(counter == 2*dim)) {
+  if (element->boundary_id++ == 2*dim) {
+    *userops = iga->userops;
     element->atboundary = PETSC_FALSE;
-    element->atboundary_id = 0;
-    *userops = iga->userops;
-    counter++;
     return PETSC_TRUE;
   }
   *userops = 0;
-  counter = 0;
+  element->atboundary  = PETSC_FALSE;
+  element->boundary_id = -1;
   return PETSC_FALSE;
 }
 
   (*point)->nsd = element->nsd;
   (*point)->npd = element->npd;
 
-  if (PetscUnlikely(!element->atboundary)) {
+  if (PetscLikely(!element->atboundary)) {
     ierr = IGAElementBuildQuadrature(element);CHKERRQ(ierr);
     ierr = IGAElementBuildShapeFuns(element);CHKERRQ(ierr);
   } else {
-    PetscInt i = element->atboundary_id / 2;
-    PetscInt s = element->atboundary_id % 2;
+    PetscInt i = element->boundary_id / 2;
+    PetscInt s = element->boundary_id % 2;
     (*point)->count = element->nqp / element->BD[i]->nqp;
     ierr = IGAElementBuildQuadratureAtBoundary(element,i,s);CHKERRQ(ierr);
     ierr = IGAElementBuildShapeFunsAtBoundary (element,i,s);CHKERRQ(ierr);
-    PetscFunctionReturn(0);
   }
   PetscFunctionReturn(0);
 }
     PetscFunctionReturn(PETSC_ERR_PLIB);
   }
   *point = PETSC_NULL;
-  PetscFunctionReturn(0);
+  /* XXX */
+  if (PetscLikely(!element->collocation)) PetscFunctionReturn(0);
+  if (PetscLikely(!element->atboundary))  PetscFunctionReturn(0);
+  element->atboundary  = PETSC_FALSE;
+  element->boundary_id = 2*element->dim;
+  /* XXX */
+ PetscFunctionReturn(0);
 }
 
 #undef  __FUNCT__
 }
 
 PETSC_STATIC_INLINE
-PetscReal DOT(PetscInt dim, PetscReal a[], PetscReal b[])
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
 {
-  PetscInt i; PetscReal s;
-  for (s=0.0, i=0; i<dim; i++) s += a[i]*b[i];
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
   return s;
 }