Commits

Lisandro Dalcin committed a686184

A bunch of fixes here and there to fix boundary integrals

Comments (0)

Files changed (5)

demo/BoundaryIntegral.c

 #define __FUNCT__ "System"
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscReal *N0 = p->shape[0];
-  PetscReal (*N1)[2] = (PetscReal (*)[2]) p->shape[1];
-  PetscInt a,b,nen=p->nen;
+  PetscInt nen,dim;
+  IGAPointGetSizes(p,0,&nen,0);
+  IGAPointGetDims(p,&dim,0);
+
+  const PetscReal *N1;
+  IGAPointGetShapeFuns(p,1,&N1);
+
+  PetscInt a,b,i;
   for (a=0; a<nen; a++) {
-    PetscReal Na   = N0[a];
-    PetscReal Na_x = N1[a][0];
-    PetscReal Na_y = N1[a][1];
     for (b=0; b<nen; b++) {
-      PetscReal Nb_x = N1[b][0];
-      PetscReal Nb_y = N1[b][1];
-      K[a*nen+b] = Na_x*Nb_x + Na_y*Nb_y;
+      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] = Na * 0.0;
   }
   return 0;
 }
 
+
+
 #undef  __FUNCT__
 #define __FUNCT__ "Neumann"
 PetscErrorCode Neumann(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 
 typedef struct { 
   PetscInt dir;
+  PetscInt side;
 } AppCtx;
 
 #undef  __FUNCT__
   AppCtx *user = (AppCtx *)ctx;
   PetscScalar u;
   IGAPointFormValue(p,U,&u);
-  PetscReal e = u - p->point[user->dir];
+  PetscReal x;
+  if (user->side == 0)
+    x = 1 - p->point[user->dir];
+  else
+    x = p->point[user->dir];
+  PetscReal e = u - x;
   S[0] = e*e;
   return 0;
 }
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
   AppCtx user;
-  user.dir = 0;
+  user.dir  = 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("-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 = IGASetDof(iga,1);CHKERRQ(ierr);
 
   IGABoundary bnd;
-  ierr = IGAGetBoundary(iga,user.dir,0,&bnd);CHKERRQ(ierr);
+  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,1,&bnd);CHKERRQ(ierr); 
+  ierr = IGAGetBoundary(iga,user.dir,n,&bnd);CHKERRQ(ierr); 
   ierr = IGABoundarySetUserSystem(bnd,Neumann,PETSC_NULL);CHKERRQ(ierr);
 
   ierr = IGASetFromOptions(iga);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));
                        /*1: [nqp][nen][nsd]           */
                        /*2: [nqp][nen][nsd][nsd]      */
                        /*3: [nqp][nen][nsd][nsd][nsd] */
+  
+  PetscReal *normal;   /*   [nqp][dim]                */
+  
 
   IGA      parent;
   IGAPoint iterator;
                        /*2: [nen][nsd][nsd] */
                        /*3: [nen][nsd][nsd][nsd] */
 
+  PetscReal *normal;   /*   [dim] */
+
   IGAElement parent;
 
   PetscInt    nvec;
   ierr = PetscFree(element->shape[2]);CHKERRQ(ierr);
   ierr = PetscFree(element->shape[3]);CHKERRQ(ierr);
 
+  ierr = PetscFree(element->normal);CHKERRQ(ierr);
+
   ierr = IGAElementFreeFix(element);CHKERRQ(ierr);
   ierr = IGAElementFreeWork(element);CHKERRQ(ierr);
   ierr = IGAPointReset(element->iterator);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*nen*dim*dim,PetscReal,&element->shape[2]);CHKERRQ(ierr);
     ierr = PetscMalloc1(nqp*nen*dim*dim*dim,PetscReal,&element->shape[3]);CHKERRQ(ierr);
 
+
     ierr = PetscMemzero(element->weight,  sizeof(PetscReal)*nqp);CHKERRQ(ierr);
     ierr = PetscMemzero(element->detJac,  sizeof(PetscReal)*nqp);CHKERRQ(ierr);
 
     ierr = PetscMemzero(element->shape[1],sizeof(PetscReal)*nqp*nen*dim);CHKERRQ(ierr);
     ierr = PetscMemzero(element->shape[2],sizeof(PetscReal)*nqp*nen*dim*dim);CHKERRQ(ierr);
     ierr = PetscMemzero(element->shape[3],sizeof(PetscReal)*nqp*nen*dim*dim*dim);CHKERRQ(ierr);
+
+    ierr = PetscMalloc1(nqp*dim,PetscReal,&element->normal);CHKERRQ(ierr);
+    ierr = PetscMemzero(element->normal,sizeof(PetscReal)*nqp*dim);CHKERRQ(ierr);
   }
   if (!element->collocation) {
     element->neq    = element->nen;
   point->shape[2] += nen*dim*dim;
   point->shape[3] += nen*dim*dim*dim;
 
+  point->normal   += dim;
+
   return PETSC_TRUE;
 
  start:
     point->shape[3] = element->basis[3];
   }
 
+  point->normal = element->normal;
+
   return PETSC_TRUE;
 
  stop:
     IGABasis *BD = element->BD;
     PetscInt *ID = element->ID;
 
-    PetscInt ione=1;
     PetscReal U[2],one=1;
     IGAAxisGetLimits(element->parent->axis[dir],&U[0],&U[1]);
 
     switch (element->dim) {
     case 3:
-      if(dir==0){
-	IGA_Quadrature_3D(ione,&U[side],&one,&one,
-			  IGA_Quadrature_ARGS(ID,BD,1),
-			  IGA_Quadrature_ARGS(ID,BD,2),
-			  element->weight,element->detJac,
-			  element->point,element->scale);
-      }else if(dir==1){
-	IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
-			  ione,&U[side],&one,&one,
-			  IGA_Quadrature_ARGS(ID,BD,2),
-			  element->weight,element->detJac,
-			  element->point,element->scale);
-      }else{
-	IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
-			  IGA_Quadrature_ARGS(ID,BD,1),
-			  ione,&U[side],&one,&one,
-			  element->weight,element->detJac,
-			  element->point,element->scale);
+      switch (dir) {
+      case 0: IGA_Quadrature_3D(1,&U[side],&one,&one,
+                                IGA_Quadrature_ARGS(ID,BD,1),
+                                IGA_Quadrature_ARGS(ID,BD,2),
+                                element->weight,element->detJac,
+                                element->point,element->scale); break;
+      case 1: IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
+                                1,&U[side],&one,&one,
+                                IGA_Quadrature_ARGS(ID,BD,2),
+                                element->weight,element->detJac,
+                                element->point,element->scale); break;
+      case 2: IGA_Quadrature_3D(IGA_Quadrature_ARGS(ID,BD,0),
+                                IGA_Quadrature_ARGS(ID,BD,1),
+                                1,&U[side],&one,&one,
+                                element->weight,element->detJac,
+                                element->point,element->scale); break;
       }
       break;
     case 2:
-      if(dir==0){
-	IGA_Quadrature_2D(ione,&U[side],&one,&one,
-			  IGA_Quadrature_ARGS(ID,BD,1),
-			  element->weight,element->detJac,
-			  element->point,element->scale);
-      }else{
-	IGA_Quadrature_2D(IGA_Quadrature_ARGS(ID,BD,0),
-			  ione,&U[side],&one,&one,
-			  element->weight,element->detJac,
-			  element->point,element->scale);
+      switch (dir) {
+      case 0: IGA_Quadrature_2D(1,&U[side],&one,&one,
+                                IGA_Quadrature_ARGS(ID,BD,1),
+                                element->weight,element->detJac,
+                                element->point,element->scale); break;
+      case 1: IGA_Quadrature_2D(IGA_Quadrature_ARGS(ID,BD,0),
+                                1,&U[side],&one,&one,
+                                element->weight,element->detJac,
+                                element->point,element->scale); break;
       }
       break;
-    case 1: IGA_Quadrature_1D(IGA_Quadrature_ARGS(ID,BD,0),
-			      element->weight,element->detJac,
-			      element->point,element->scale); break;
+    case 1:
+      switch (dir) {
+      case 0: IGA_Quadrature_1D(1,&U[side],&one,&one,
+                                element->weight,element->detJac,
+                                element->point,element->scale); break;
+      }
+      break;
     }
   }
   PetscFunctionReturn(0);
   PetscFunctionReturn(0);
 }
 
+EXTERN_C_BEGIN
 extern void IGA_DersBasisFuns(PetscInt i,PetscReal u,PetscInt p,PetscInt d,const PetscReal U[],PetscReal N[]);
+EXTERN_C_END
+
+EXTERN_C_BEGIN
+extern void IGA_GetNormal(PetscInt dim,PetscInt dir,PetscInt side,const PetscReal F[],PetscReal *dS,PetscReal N[]);
+EXTERN_C_END
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementBuildBoundaryShapeFuns"
     PetscReal **N = element->basis;
 
     PetscInt k = element->parent->axis[dir]->span[element->ID[dir]];
-    PetscReal u = element->point[dir];
     PetscInt p = element->parent->axis[dir]->p;
     PetscInt d = element->parent->order;
-    PetscInt ione = 1;
     PetscReal *U = element->parent->axis[dir]->U;
     PetscInt nen = BD[dir]->nen;
     PetscReal value[nen*(d+1)];
-    IGA_DersBasisFuns(k,u,p,d,U,&value[0]);
+
+    PetscReal u[2];
+    IGAAxisGetLimits(element->parent->axis[dir],&u[0],&u[1]);
+    IGA_DersBasisFuns(k,u[side],p,d,U,value);
 
     switch (element->dim) {
     case 3:
-      if(dir == 0){
-	IGA_BasisFuns_3D(order,element->rational,
-			 element->rationalW,
-			 ione,BD[0]->nen,BD[0]->d,&value[0],
-			 IGA_BasisFuns_ARGS(ID,BD,1),
-			 IGA_BasisFuns_ARGS(ID,BD,2),
-			 N[0],N[1],N[2],N[3]);
-      }else if(dir == 1){
-	IGA_BasisFuns_3D(order,element->rational,
-			 element->rationalW,
-			 IGA_BasisFuns_ARGS(ID,BD,0),
-			 ione,BD[1]->nen,BD[1]->d,&value[0],
-			 IGA_BasisFuns_ARGS(ID,BD,2),
-			 N[0],N[1],N[2],N[3]);
-      }else{
-	IGA_BasisFuns_3D(order,element->rational,
-			 element->rationalW,
-			 IGA_BasisFuns_ARGS(ID,BD,0),
-			 IGA_BasisFuns_ARGS(ID,BD,1),
-			 ione,BD[2]->nen,BD[2]->d,&value[0],
-			 N[0],N[1],N[2],N[3]);
+      switch (dir) {
+      case 0: IGA_BasisFuns_3D(order,
+                               element->rational, element->rationalW,
+                               1,BD[0]->nen,BD[0]->d,value,
+                               IGA_BasisFuns_ARGS(ID,BD,1),
+                               IGA_BasisFuns_ARGS(ID,BD,2),
+                               N[0],N[1],N[2],N[3]); break;
+      case 1: IGA_BasisFuns_3D(order,
+                               element->rational,element->rationalW,
+                               IGA_BasisFuns_ARGS(ID,BD,0),
+                               1,BD[1]->nen,BD[1]->d,value,
+                               IGA_BasisFuns_ARGS(ID,BD,2),
+                               N[0],N[1],N[2],N[3]); break;
+      case 2: IGA_BasisFuns_3D(order,
+                               element->rational,element->rationalW,
+                               IGA_BasisFuns_ARGS(ID,BD,0),
+                               IGA_BasisFuns_ARGS(ID,BD,1),
+                               1,BD[2]->nen,BD[2]->d,value,
+                               N[0],N[1],N[2],N[3]); break;
       }
       break;
     case 2:
-      if(dir == 0){
-	IGA_BasisFuns_2D(order,element->rational,
-			 element->rationalW,
-			 ione,BD[0]->nen,BD[0]->d,&value[0],
-			 IGA_BasisFuns_ARGS(ID,BD,1),
-			 N[0],N[1],N[2],N[3]);
-      }else{
-	IGA_BasisFuns_2D(order,element->rational,
-			 element->rationalW,
-			 IGA_BasisFuns_ARGS(ID,BD,0),
-			 ione,BD[1]->nen,BD[1]->d,&value[0],
-			 N[0],N[1],N[2],N[3]);
+      switch (dir) {
+      case 0: IGA_BasisFuns_2D(order,
+                               element->rational,element->rationalW,
+                               1,BD[0]->nen,BD[0]->d,value,
+                               IGA_BasisFuns_ARGS(ID,BD,1),
+                               N[0],N[1],N[2],N[3]); break;
+      case 1: IGA_BasisFuns_2D(order,
+                               element->rational,element->rationalW,
+                               IGA_BasisFuns_ARGS(ID,BD,0),
+                               1,BD[1]->nen,BD[1]->d,value,
+                               N[0],N[1],N[2],N[3]); break;
       }
       break;
-    case 1: IGA_BasisFuns_1D(order,element->rational,
-			     element->rationalW,
-			     IGA_BasisFuns_ARGS(ID,BD,0),
-			     N[0],N[1],N[2],N[3]); break;
+    case 1: 
+      switch (dir) {
+      case 0: IGA_BasisFuns_1D(order,
+                               element->rational,element->rationalW,
+                               1,BD[0]->nen,BD[0]->d,value,
+                               N[0],N[1],N[2],N[3]); break;
+      }
+      break;
     }
   }
   if (element->dim == element->nsd) /* XXX */
 			     N[0],N[1],N[2],N[3],
 			     dX,gX[0],gX[1]); break;
     }
+    {
+      PetscInt  dim = element->dim;
+      PetscReal *F  = element->gradX[0];
+      PetscReal dS, *n = element->normal;
+      for (q=0; q<nqp; q++) {
+        IGA_GetNormal(dim,dir,side,F,&dS,&n[q*dim]);
+        element->detJac[q] *= dS;
+      }
+    }
+  }
+  if (!element->geometry) {
+    PetscInt q, nqp = element->nqp;
+    PetscInt i, dim = element->dim;
+    PetscReal *n  = element->normal;
     for (q=0; q<nqp; q++)
-      element->detJac[q] *= dX[q];
+      for (i=0; i<dim; i++)
+        n[q*dim+i] = (dir==i) ? ((side==0)?-1.0:1.0) : 0.0;
   }
   PetscFunctionReturn(0);
 }

src/petigaftn.F90

      type(C_PTR) :: detX
      type(C_PTR) :: gradX(0:1)
      type(C_PTR) :: shape(0:3)
+     type(C_PTR) :: normal
   end type IGAPoint
 
 
       call IGAPoint_InvGradGeomMap(p,G)
     end function IGA_InvGradGeomMap
 
+    function IGA_Normal(p) result (N)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      real(kind=IGA_REAL_KIND), pointer :: N(:)
+      call c2f(p%normal,N,(/p%dim/))
+    end function IGA_Normal
+
     function IGA_Shape0(p) result(N)
       use ISO_C_BINDING, only: c2f => C_F_POINTER
       implicit none

src/petigaval.F90

 include 'petigainv.f90.in'
 end subroutine IGA_GetInvGradGeomMap
 
+subroutine IGA_GetNormal(dim,axis,side,F,dS,N) &
+  bind(C, name="IGA_GetNormal")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: dim
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: axis,side
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: F(dim,dim)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: dS
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: N(dim)
+  select case (dim)
+  case (3)
+     select case (axis)
+     case (0)
+        N = normal3(F(2,:),F(3,:))
+     case (1)
+        N = normal3(F(3,:),F(1,:))
+     case (2)
+        N = normal3(F(1,:),F(2,:))
+     end select
+  case (2)
+     select case (axis)
+     case (0); N = normal2(F(2,:)) 
+     case (1); N = normal2(F(1,:)) 
+     end select
+  case (1)
+     N = (/ 1.0 /)
+  end select
+  dS = sqrt(sum(N*N))
+  if (side == 0) then
+     N = -N/dS
+  else
+     N = +N/dS
+  endif
+contains
+function normal2(t) result(n)
+  implicit none
+  real(kind=IGA_REAL_KIND)             :: n(2)
+  real(kind=IGA_REAL_KIND), intent(in) :: t(2)
+  ! n_i = eps_ij n_j
+  n(1) = + t(2)
+  n(2) = - t(1)
+end function normal2
+function normal3(s,t) result(n)
+  implicit none
+  real(kind=IGA_REAL_KIND)             :: n(3)
+  real(kind=IGA_REAL_KIND), intent(in) :: s(3)
+  real(kind=IGA_REAL_KIND), intent(in) :: t(3)
+  ! c_i = eps_ijk s_j tj
+  n(1) = s(2) * t(3) - s(3) * t(2)
+  n(2) = s(3) * t(1) - s(1) * t(3)
+  n(3) = s(1) * t(2) - s(2) * t(1)
+end function normal3
+end subroutine IGA_GetNormal
+
 
 subroutine IGA_GetValue(nen,dof,N,U,V) &
   bind(C, name="IGA_GetValue")