Commits

Lisandro Dalcin committed f9e155e

Improve handling of BCs for collocation

Comments (0)

Files changed (4)

+/*
+  This code solves the Laplace problem in one of the following ways:
+
+  1) On the parametric unit domain [0,1]^dim (default)
+
+    To solve on the parametric domain, do not specify a geometry
+    file. You may change the discretization by altering the dimension
+    of the space (-iga_dim), the number of uniform elements in each
+    direction (-iga_elements), the polynomial order (-iga_degree),
+    and the continuity (-iga_continuity).
+
+  2) On a geometry
+
+    If a geometry file is specified (-iga_geometry), the
+    discretization will be what is read in from the geometry and is
+    not editable from the command line.
+
+  Note that the boundary conditions for this problem are such that
+  the solution is always u(x)=1 (unit Dirichlet on the left side and
+  free Neumann on the right). The error in the solution may be
+  computed by using the -print_error command.
+
+*/
+
+
 #include "petiga.h"
 
 #undef  __FUNCT__
-#define __FUNCT__ "SystemLaplace"
-PetscErrorCode SystemLaplace(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+#define __FUNCT__ "SystemGalerkin"
+PetscErrorCode SystemGalerkin(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
 
   const PetscReal *N1;
   IGAPointGetShapeFuns(p,1,&N1);
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "SystemPoisson"
-PetscErrorCode SystemPoisson(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+#define __FUNCT__ "SystemCollocation"
+PetscErrorCode SystemCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
 
-  const PetscReal *N;
-  IGAPointGetShapeFuns(p,0,&N);
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,1,&N1);
+  const PetscReal (*N2)[dim][dim];
+  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
 
-  PetscInt a,b,i;
-  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;
-    }
-    F[a] = 1.0*N[a];
-  }
+  PetscInt a,i;
+  for (a=0; a<nen; a++)
+    for (i=0; i<dim; i++)
+      K[a] += -N2[a][i][i];
+
+  F[0] = 0.0;
+
   return 0;
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "SystemCollocation"
-PetscErrorCode SystemCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
-{
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
-
-  PetscInt Nb[3] = {0,0,0};
-  Nb[0] = p->parent->parent->axis[0]->nnp;
-  Nb[1] = p->parent->parent->axis[1]->nnp;
-  Nb[2] = p->parent->parent->axis[2]->nnp;
-
-  const PetscReal *N0,(*N1)[dim],(*N2)[dim][dim];
-  IGAPointGetBasisFuns(p,0,(const PetscReal**)&N0);
-  IGAPointGetBasisFuns(p,1,(const PetscReal**)&N1);
-  IGAPointGetBasisFuns(p,2,(const PetscReal**)&N2);
-  
-  PetscInt a,i;
-  PetscReal n[3] = {0,0,0};
-  PetscBool Dirichlet=PETSC_FALSE,Neumann=PETSC_FALSE;
-  for (i=0; i<dim; i++) if (p->parent->ID[i] == 0) Dirichlet=PETSC_TRUE;
-  if (!Dirichlet) { 
-    for (i=0; i<dim; i++) {
-      if (p->parent->ID[i] == Nb[i]-1) {
-	Neumann=PETSC_TRUE;
-	n[i] = 1;
-	i=dim;
-      }
-    }
-  }
-
-  if(Dirichlet){
-    for (a=0; a<nen; a++) K[a] += N0[a];
-    F[0] = 1.0;
-  }else if(Neumann){
-    for (a=0; a<nen; a++) 
-      for (i=0; i<dim; i++) 
-	K[a] += N1[a][i]*n[i];    
-    F[0] = 0.0;
-  }else{
-    for (a=0; a<nen; a++) 
-      for (i=0; i<dim; i++) 
-	K[a] += -N2[a][i][i];
-  }
-  return 0;
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "ErrorLaplace"
-PetscErrorCode ErrorLaplace(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
+#define __FUNCT__ "Error"
+PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
 {
   PetscScalar u;
   IGAPointFormValue(p,U,&u);
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
 
-  /*
-    This code solves the Laplace problem in one of the following ways:
-    
-    1) On the parametric unit domain [0,1]^dim (default)
-    
-      To solve on the parametric domain, do not specify a geometry
-      file. You may change the discretization by altering the
-      dimension of the space (-dim), the number of uniform elements in
-      each direction (-iga_elements), the polynomial order
-      (-iga_degree), and the continuity (-iga_continuity). Note that
-      the boundary conditions for this problem are such that the
-      solution is always u(x)=1 (unit Dirichlet on the left side and
-      free Neumann on the right). The error in the solution may be
-      computed by using the -print_error command.
-
-    2) On a geometry
-
-      If a geometry file is specified (-iga_geometry), then the code
-      will solve the Poisson problem on this geometry. The forcing is
-      set to 1 and we have 0 Dirichlet conditions everywhere. We use
-      this mode to test geometries as the solution should display the
-      same symmetry as the geometry. The discretization will be what
-      is read in from the geometry and is not editable from the
-      commandline.
-
-   */
-
   PetscErrorCode ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
   // Setup options
 
-  PetscInt  i; 
-  PetscInt  dim = 3; 
-  PetscBool print_error = PETSC_FALSE; 
-  PetscBool draw = PETSC_FALSE; 
-  PetscBool Collocation = PETSC_FALSE;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Laplace Options","IGA");CHKERRQ(ierr); 
-  ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr); 
-  ierr = PetscOptionsBool("-print_error","Prints the L2 error of the solution",__FILE__,print_error,&print_error,PETSC_NULL);CHKERRQ(ierr); 
-  ierr = PetscOptionsBool("-draw","If dim <= 2, then draw the solution to the screen",__FILE__,draw,&draw,PETSC_NULL);CHKERRQ(ierr); 
-  ierr = PetscOptionsBool("-collocation","Enable to use collocation",__FILE__,Collocation,&Collocation,PETSC_NULL);CHKERRQ(ierr); 
-  ierr = PetscOptionsEnd();CHKERRQ(ierr); 
+  PetscInt  i;
+  PetscBool collocation = PETSC_FALSE;
+  PetscBool print_error = PETSC_FALSE;
+  PetscBool save = PETSC_FALSE;
+  PetscBool draw = PETSC_FALSE;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Laplace Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-collocation","Enable to use collocation",__FILE__,collocation,&collocation,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-print_error","Prints the L2 error of the solution",__FILE__,print_error,&print_error,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-save","Save the solution to file",__FILE__,save,&save,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-draw","If dim <= 2, then draw the solution to the screen",__FILE__,draw,&draw,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
   // Initialize the discretization
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
-  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  if (iga->dim < 1) {ierr = IGASetDim(iga,3);CHKERRQ(ierr);}
   ierr = IGASetUp(iga);CHKERRQ(ierr);
+  if (collocation) {ierr = IGASetUseCollocation(iga,PETSC_TRUE);CHKERRQ(ierr);}
 
   // Set boundary conditions
 
-  if (iga->geometry) {
-    for(i=0; i<dim; i++) {
-      IGABoundary bnd;
-      ierr = IGAGetBoundary(iga,i,0,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-      ierr = IGAGetBoundary(iga,i,1,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-    }
-  }else{
-    for (i=0; i<dim; i++) {
-      IGABoundary bnd;
-      ierr = IGAGetBoundary(iga,i,0,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,1.0);CHKERRQ(ierr);
-    }
+  PetscInt dim;
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
+  for (i=0; i<dim; i++) {
+    IGABoundary bnd;
+    ierr = IGAGetBoundary(iga,i,0,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetValue(bnd,0,1.0);CHKERRQ(ierr);
+    ierr = IGAGetBoundary(iga,i,1,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetLoad(bnd,0,0.0);CHKERRQ(ierr);
   }
-  
+
   // Assemble
 
   Mat A;
-  Vec x,b;
+  Vec b,x;
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
-  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  if (Collocation){
-    ierr = IGASetUseCollocation(iga,PETSC_TRUE);CHKERRQ(ierr); 
-    ierr = IGASetUserSystem(iga,SystemCollocation,PETSC_NULL);CHKERRQ(ierr); 
-    ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
-  }else{
-    if (iga->geometry){ 
-      ierr = IGASetUserSystem(iga,SystemPoisson,PETSC_NULL);CHKERRQ(ierr); 
-    }else{
-      ierr = IGASetUserSystem(iga,SystemLaplace,PETSC_NULL);CHKERRQ(ierr); 
-    }
-    ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+  if (!iga->collocation) {
+    ierr = IGASetUserSystem(iga,SystemGalerkin,PETSC_NULL);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
+  } else {
+    ierr = IGASetUserSystem(iga,SystemCollocation,PETSC_NULL);CHKERRQ(ierr);
+    ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_FALSE);CHKERRQ(ierr);
   }
+  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   // Solve
 
 
   // Various post-processing options
 
-  if (iga->geometry) {
-    MPI_Comm        comm;
-    PetscViewer     viewer;
-    ierr = PetscObjectGetComm((PetscObject)x,&comm);CHKERRQ(ierr);
-    ierr = PetscViewerBinaryOpen(comm,"solution.dat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
-    ierr = VecView(x,viewer);CHKERRQ(ierr);
-    ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
-    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
-  }
-
-  if (print_error && !iga->geometry) {
+  if (print_error) {
     PetscScalar error = 0;
-    ierr = IGAFormScalar(iga,x,1,&error,ErrorLaplace,PETSC_NULL);CHKERRQ(ierr);
+    ierr = IGAFormScalar(iga,x,1,&error,Error,PETSC_NULL);CHKERRQ(ierr);
     error = PetscSqrtReal(PetscRealPart(error));
     ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);
   }
-
-  if (draw && dim <= 2) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
+  if (save) {
+    ierr = IGAWrite(iga,"LaplaceGeometry.dat");CHKERRQ(ierr);
+    ierr = IGAWriteVec(iga,x,"LaplaceSolution.dat");CHKERRQ(ierr);
+  }
+  if (draw && dim <= 2) {
+    ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
+  }
 
   // Cleanup
 
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidLogicalCollectiveBool(iga,collocation,2);
-  if (collocation && !iga->collocation) {
+  if (collocation) {
     PetscMPIInt size = 1;
-    PetscInt i, dim = (iga->dim > 0) ? iga->dim : 3;
-    PetscBool periodic = PETSC_FALSE;
     ierr = MPI_Comm_size(((PetscObject)iga)->comm,&size);CHKERRQ(ierr);
-    for (i=0; i<dim; i++) if(iga->axis[i]->periodic) periodic = PETSC_TRUE;
     if (size > 1) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_SUP,
                           "Collocation not supported in parallel");
+  }
+  if (collocation && iga->setup) {
+    PetscInt i, dim = iga->dim;
+    PetscBool periodic = PETSC_FALSE;
+    for (i=0; i<dim; i++) if(iga->axis[i]->periodic) periodic = PETSC_TRUE;
     if (periodic) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_SUP,
                           "Collocation not supported with periodicity");
   }
-  if (collocation && iga->setup) { /* collocation */
+  if (collocation && iga->setup) {
     PetscInt i;
     for (i=0; i<3; i++) {
       ierr = IGABasisDestroy(&iga->node_basis[i]);CHKERRQ(ierr);

src/petigabasis.c

                       "Derivative order must be grather than zero, got %D",d);
 
   p = axis->p;
-  n = axis->m - p -1;
+  n = axis->m - p - 1;
   U = axis->U;
 
   nel  = axis->nnp;
     PetscInt  k = IGA_FindSpan(n,p,u,U);
     PetscReal *N = &value[iel*nen*ndr];
     offset[iel] = k-p;
+    detJ[iel]   = U[k+1]-U[k];
     point[iel]  = u;
-    detJ[iel]   = U[k+1]-U[k];
     IGA_Basis_BSpline(k,u,p,d,U,N);
   }
 
   basis->point  = point;
   basis->value  = value;
 
+  {
+    PetscInt  k0 = p,    k1 = n;
+    PetscReal u0 = U[p], u1 = U[n+1];
+    ierr = PetscMalloc1(nen*ndr,PetscReal,&basis->bnd_value[0]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nen*ndr,PetscReal,&basis->bnd_value[1]);CHKERRQ(ierr);
+    basis->bnd_offset[0] = k0-p; basis->bnd_offset[1] =  k1-p;
+    basis->bnd_detJ  [0] =  1.0; basis->bnd_detJ  [1] =   1.0;
+    basis->bnd_weight[0] =  1.0; basis->bnd_weight[1] =   1.0;
+    basis->bnd_point [0] =   u0; basis->bnd_point [1] =    u1;
+    IGA_Basis_BSpline(k0,u0,p,d,U,basis->bnd_value[0]);
+    IGA_Basis_BSpline(k1,u1,p,d,U,basis->bnd_value[1]);
+  }
+
   PetscFunctionReturn(0);
 }
 
   }
 }
 
+PETSC_STATIC_INLINE
+IGABoundary AtBoundary(IGAElement element,PetscInt dir,PetscInt side)
+{
+  IGABoundary b = element->parent->boundary[dir][side];
+  PetscInt e = side ? element->sizes[dir]-1 : 0;
+  return (element->ID[dir] == e) ? b : PETSC_NULL;
+}
+
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim, PetscReal a[], PetscReal b[])
+{
+  PetscInt i; PetscReal s;
+  for (s=0.0, i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementBuildFix"
 PetscErrorCode IGAElementBuildFix(IGAElement element)
   PetscFunctionReturn(0);
  collocation:
   {
-    PetscInt M = element->neq * element->dof;
-    PetscInt N = element->nen * element->dof;
-    PetscInt f,n;
-    n = element->nfix;
-    for (f=0; f<n; f++) {
-      PetscInt k = element->ifix[f];
-      PetscInt a = k / element->dof;
-      if (element->index == element->colmap[a]) {
-        PetscInt    c = k % element->dof;
-        PetscScalar v = element->vfix[f];
-        PetscInt i,j;
-        for (i=0; i<M; i++) K[i*N+k] = 0.0;
-        for (j=0; j<N; j++) K[c*N+j] = 0.0;
-        K[c*N+k] = 1.0;
-        F[c] = v;
+    PetscInt dim = element->dim;
+    PetscInt dof = element->dof;
+    PetscInt nen = element->nen;
+    PetscInt N = nen * dof;
+    PetscInt dir,side;
+    for (dir=0; dir<dim; dir++) {
+      for (side=0; side<2; side++) {
+        IGABoundary b = AtBoundary(element,dir,side);
+        if(b && b->nload) {
+          PetscInt  f,n = b->nload;
+          PetscReal normal[3] = {0.0,0.0,0.0};
+          PetscReal *dshape;
+          if (!element->geometry) {
+            normal[dir] = side ? 1.0 : -1.0;
+            dshape = element->basis[1];
+          } else {
+            PetscReal dS, *F = element->gradX[0];
+            IGA_GetNormal(dim,dir,side,F,&dS,normal);
+            dshape = element->shape[1];
+          }
+          for (f=0; f<n; f++) {
+            PetscInt  c = b->iload[f];
+            PetscReal v = b->vload[f];
+            PetscInt  a,j;
+            for (j=0; j<N; j++) K[c*N+j] = 0.0;
+            for (a=0; a<nen; a++)
+              K[c*N+a*dof+c] = DOT(dim,&dshape[a*dim],normal);
+            F[c] = v;
+          }
+        }
+        if (b && b->count) {
+          PetscInt  f,n = b->count;
+          PetscReal *shape;
+          if (!element->geometry) {
+            shape = element->basis[0];
+          } else {
+            shape = element->shape[0];
+          }
+          for (f=0; f<n; f++) {
+            PetscInt  c = b->field[f];
+            PetscReal v = b->value[f];
+            PetscInt  a,j;
+            for (j=0; j<N; j++) K[c*N+j] = 0.0;
+            for (a=0; a<nen; a++)
+              K[c*N+a*dof+c] = shape[a];
+            F[c] = v;
+          }
+        }
       }
     }
   }
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.