Commits

Lisandro Dalcin  committed cc2d452

Implement constant Neumann boundary conditions

  • Participants
  • Parent commits f9566d6

Comments (0)

Files changed (8)

File demo/Neumann.c

+#include "petiga.h"
+
+#define pi M_PI
+
+
+PetscReal Exact(PetscReal x[3])
+{
+  return sin(pi*x[0]) + sin(pi*x[1]) + sin(pi*x[2]);
+}
+
+PetscReal Forcing(PetscReal x[3])
+{
+  return pi*pi * Exact(x);
+}
+
+#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);
+
+  PetscReal x[3] = {0,0,0};
+  IGAPointFormPoint(p,x);
+  PetscReal f = Forcing(x);
+
+  const PetscReal *N0,*N1;
+  IGAPointGetShapeFuns(p,0,&N0);
+  IGAPointGetShapeFuns(p,1,&N1);
+
+  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] = N0[a]*f;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "Error"
+PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
+{
+  PetscReal x[3] = {0,0,0};
+  IGAPointFormPoint(p,x);
+  PetscScalar u;
+  IGAPointFormValue(p,U,&u);
+  PetscReal e = PetscAbsScalar(u) -  Exact(x);
+  S[0] = e*e;
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  PetscInt dim = 3;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Neumann Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+
+  IGABoundary bnd;
+  PetscInt dir,side;
+  for (dir=0; dir<dim; dir++) {
+    for (side=0; side<2; side++) {
+      PetscScalar load = -pi;
+      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+      ierr = IGABoundarySetLoad(bnd,0,load);CHKERRQ(ierr);
+    }
+  }
+
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  Mat A;
+  Vec x,b;
+  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);
+  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+
+  KSP ksp;
+  ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
+  ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
+  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+  {
+    MPI_Comm comm;
+    MatNullSpace nsp;
+    ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
+    ierr = MatNullSpaceCreate(comm,PETSC_TRUE,0,PETSC_NULL,&nsp);CHKERRQ(ierr);
+    ierr = KSPSetNullSpace(ksp,nsp);CHKERRQ(ierr);
+    ierr = MatNullSpaceDestroy(&nsp);CHKERRQ(ierr);
+  }
+  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+  {
+    PetscReal vmin;
+    ierr = VecMin(x,0,&vmin);CHKERRQ(ierr);
+    ierr = VecShift(x,-vmin);CHKERRQ(ierr);
+  }
+
+  PetscScalar se = 0;
+  ierr = IGAFormScalar(iga,x,1,&se,Error,PETSC_NULL);CHKERRQ(ierr);
+  PetscReal e = PetscSqrtReal(PetscRealPart(se));
+
+  PetscBool error = PETSC_FALSE;
+  ierr = PetscOptionsGetBool(0,"-error",&error,0);CHKERRQ(ierr);
+  if (error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"Error=%G\n",e);CHKERRQ(ierr);}
+  else if (e>1e-4) SETERRQ1(PETSC_COMM_WORLD,1,"Error=%G\n",e);
+
+  if (dim < 3) {
+    PetscBool draw = PETSC_FALSE;
+    ierr = PetscOptionsGetBool(0,"-draw",&draw,0);CHKERRQ(ierr);
+    if (draw) {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);
+
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}

File demo/makefile

 	  AdvectionDiffusion \
 	  Laplace \
 	  Poisson \
+	  Neumann \
 	  Bratu \
 	  PatternFormation \
 	  CahnHilliard2D \
 Poisson: Poisson.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
 	${RM} -f $<
+Neumann: Neumann.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
 LoggChallenge: LoggChallenge.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
 	${RM} -f $<
 	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 2 -steady false -ts_max_steps 2
 runex6d_1:
 	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 2 -steady -iga_collocation
+runex7a_1:
+	-@${MPIEXEC} -n 1 ./Neumann ${OPTS} -dim 1
+runex7a_4:
+	-@${MPIEXEC} -n 4 ./Neumann ${OPTS} -dim 1
+runex7b_1:
+	-@${MPIEXEC} -n 1 ./Neumann ${OPTS} -dim 2
+runex7b_4:
+	-@${MPIEXEC} -n 4 ./Neumann ${OPTS} -dim 2
+runex7c_4:
+	-@${MPIEXEC} -n 4 ./Neumann ${OPTS} -dim 3
+
 
 InputOutput := InputOutput.PETSc runex0a_1 runex0a_2 runex0a_3 runex0a_4 runex0a_8 runex0a.rm InputOutput.rm
 L2Proj1D := L2Proj1D.PETSc runex1a_1 runex1a_4 L2Proj1D.rm
 Poisson1D := Poisson1D.PETSc runex2a_1 runex2a_4 Poisson1D.rm
 Poisson2D := Poisson2D.PETSc runex2b_1 runex2b_4 Poisson2D.rm
 Poisson3D := Poisson3D.PETSc runex2c_1 runex2c_4 Poisson3D.rm
+Neumann := Neumann.PETSc runex7a_1 runex7a_4 runex7b_1 runex7b_4 runex7c_4 Neumann.rm
 Bratu := Bratu.PETSc runex6a_1 runex6a_2 runex6b_1 runex6b_4 runex6c_1 runex6c_4 runex6d_1 Bratu.rm
 CahnHilliard2D := CahnHilliard2D.PETSc runex4_1 runex4_4 CahnHilliard2D.rm
 PatternFormation := PatternFormation.PETSc runex5a_1 runex5a_4 runex5b_1 runex5b_4 PatternFormation.rm
 Poisson      := $(Poisson1D) $(Poisson2D) $(Poisson3D)
 CahnHilliard := $(CahnHilliard2D) $(CahnHilliard3D)
 
-TESTEXAMPLES_C := $(InputOutput) $(L2Proj) $(Poisson) $(Bratu) $(CahnHilliard) $(PatternFormation)
+TESTEXAMPLES_C := $(InputOutput) $(L2Proj) $(Poisson) $(Neumann) $(Bratu) $(CahnHilliard) $(PatternFormation)
 TESTEXAMPLES_F :=
 TESTEXAMPLES_FORTRAN:=$(TESTEXAMPLES_F)
 testexamples:

File include/petiga.h

   PetscInt refct;
   /**/
   PetscInt    dof;
-  PetscInt    nbc;
+  /**/
+  PetscInt    count;
   PetscInt    *field;
   PetscScalar *value;
+  /**/
+  PetscInt    nload;
+  PetscInt    *iload;
+  PetscScalar *vload;
 };
 PETSC_EXTERN PetscErrorCode IGABoundaryCreate(IGABoundary *boundary);
 PETSC_EXTERN PetscErrorCode IGABoundaryDestroy(IGABoundary *boundary);
 PETSC_EXTERN PetscErrorCode IGABoundaryReset(IGABoundary boundary);
 PETSC_EXTERN PetscErrorCode IGABoundaryReference(IGABoundary boundary);
 PETSC_EXTERN PetscErrorCode IGABoundaryInit(IGABoundary boundary,PetscInt dof);
+PETSC_EXTERN PetscErrorCode IGABoundaryClear(IGABoundary boundary);
 PETSC_EXTERN PetscErrorCode IGABoundarySetValue(IGABoundary boundary,PetscInt field,PetscScalar value);
+PETSC_EXTERN PetscErrorCode IGABoundarySetLoad (IGABoundary boundary,PetscInt field,PetscScalar value);
 
 /* ---------------------------------------------------------------- */
 
   PetscInt     nfix;
   PetscInt    *ifix;
   PetscScalar *vfix;
-  PetscScalar *xfix;
+  PetscScalar *ufix;
+
+  PetscInt     nflux;
+  PetscInt    *iflux;
+  PetscScalar *vflux;
 
   PetscInt    nval;
   PetscScalar *wval[8];

File src/petiga1d.f90

 pure subroutine IGA_Quadrature_1D(&
-     inq,iX,iW,iL,           &
-     W,J,X,L)                &
+     inq,iX,iW,iL,                &
+     W,J,X,L)                     &
   bind(C, name="IGA_Quadrature_1D")
   use PetIGA
   implicit none
 
 
 pure subroutine IGA_BasisFuns_1D(&
-     order,                 &
-     rational,W,            &
-     inq,ina,ind,iN,        &
-     N0,N1,N2,N3)           &
+     order,                      &
+     rational,W,                 &
+     inq,ina,ind,iN,             &
+     N0,N1,N2,N3)                &
   bind(C, name="IGA_BasisFuns_1D")
   use PetIGA
   implicit none
 
 
 pure subroutine IGA_ShapeFuns_1D(&
-     order,                 &
-     nqp,nen,X,             &
-     M0,M1,M2,M3,           &
-     N0,N1,N2,N3,           &
-     DetF,F,G)              &
+     order,                      &
+     nqp,nen,X,                  &
+     M0,M1,M2,M3,                &
+     N0,N1,N2,N3,                &
+     DetF,F,G)                   &
   bind(C, name="IGA_ShapeFuns_1D")
   use PetIGA
   implicit none

File src/petiga2d.f90

 pure subroutine IGA_Quadrature_2D(&
-     inq,iX,iW,iL,           &
-     jnq,jX,jW,jL,           &
-     W,J,X,L)                &
+     inq,iX,iW,iL,                &
+     jnq,jX,jW,jL,                &
+     W,J,X,L)                     &
   bind(C, name="IGA_Quadrature_2D")
   use PetIGA
   implicit none
 
 
 pure subroutine IGA_BasisFuns_2D(&
-     order,                 &
-     rational,W,            &
-     inq,ina,ind,iN,        &
-     jnq,jna,jnd,jN,        &
-     N0,N1,N2,N3)           &
+     order,                      &
+     rational,W,                 &
+     inq,ina,ind,iN,             &
+     jnq,jna,jnd,jN,             &
+     N0,N1,N2,N3)                &
   bind(C, name="IGA_BasisFuns_2D")
   use PetIGA
   implicit none
 
 
 pure subroutine IGA_ShapeFuns_2D(&
-     order,                 &
-     nqp,nen,X,             &
-     M0,M1,M2,M3,           &
-     N0,N1,N2,N3,           &
-     DetF,F,G)              &
+     order,                      &
+     nqp,nen,X,                  &
+     M0,M1,M2,M3,                &
+     N0,N1,N2,N3,                &
+     DetF,F,G)                   &
   bind(C, name="IGA_ShapeFuns_2D")
   use PetIGA
   implicit none
 contains
 include 'petigageo.f90.in'
 end subroutine IGA_ShapeFuns_2D
+
+
+subroutine IGA_BoundaryArea_2D(&
+     m,axis,side,              &
+     geometry,Cx,              &
+     rational,Cw,              &
+     nqp,W,nen,ndr,N,          &
+     dS)                       &
+  bind(C, name="IGA_BoundaryArea_2D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter         :: nsd = 2
+  integer(kind=IGA_INTEGER_KIND), parameter         :: dim = 1
+  integer(kind=IGA_INTEGER_KIND), intent(in)        :: m(2)
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: axis, side
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: geometry, rational
+  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cx(nsd,m(1),m(2))
+  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cw(    m(1),m(2))
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: nqp, nen, ndr
+  real   (kind=IGA_REAL_KIND   ), intent(in)        :: W(nqp), N(0:ndr,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)       :: dS
+  integer(kind=IGA_INTEGER_KIND)  :: k, q
+  real   (kind=IGA_REAL_KIND   )  :: N0(nen), N1(dim,nen), detJ
+  real   (kind=IGA_REAL_KIND   ), pointer :: Xx(:,:), Xw(:)
+  select case (axis)
+  case (0)
+     if (side==0) k=1
+     if (side==1) k=m(1)
+     Xx => Cx(:,k,:); Xw => Cw(k,:)
+  case (1)
+     if (side==0) k=1
+     if (side==1) k=m(2)
+     Xx => Cx(:,:,k); Xw => Cw(:,k)
+  end select
+  detJ = 1.0
+  dS = 0.0
+  do q=1,nqp
+     N0(  :) = N(0,:,q)
+     N1(1,:) = N(1,:,q)
+     if (rational /= 0) then
+        call Rationalize(nen,Xw,N0,N1)
+     end if
+     if (geometry /= 0) then
+        call Jacobian(nen,N1,Xx,detJ)
+     end if
+     dS = dS + detJ * W(q)
+  end do
+contains
+pure subroutine Rationalize(nen,W,R0,R1)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter     :: dim = 1
+  integer(kind=IGA_INTEGER_KIND), intent(in)    :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: W(nen)
+  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R0(    nen)
+  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R1(dim,nen)
+  integer(kind=IGA_INTEGER_KIND)  :: i
+  real   (kind=IGA_REAL_KIND   )  :: W0
+  R0 = W * R0
+  W0 = sum(R0)
+  R0 = R0 / W0
+  forall(i=1:dim) &
+  R1(i,:) = W*R1(i,:) - R0 * sum(W*R1(i,:))
+  R1 = R1 / W0
+end subroutine Rationalize
+pure subroutine Jacobian(nen,N,X,J)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: nsd = 2
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: X(nsd,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: J
+  real   (kind=IGA_REAL_KIND   )  :: F(dim,nsd), M(dim,dim)
+  F = matmul(N,transpose(X))
+  M = matmul(F,transpose(F))
+  J = sqrt(abs(Determinant(dim,M)))
+end subroutine Jacobian
+include 'petigainv.f90.in'
+end subroutine IGA_BoundaryArea_2D

File src/petiga3d.f90

 pure subroutine IGA_Quadrature_3D(&
-     inq,iX,iW,iL,           &
-     jnq,jX,jW,jL,           &
-     knq,kX,kW,kL,           &
-     W,J,X,L)                &
+     inq,iX,iW,iL,                &
+     jnq,jX,jW,jL,                &
+     knq,kX,kW,kL,                &
+     W,J,X,L)                     &
   bind(C, name="IGA_Quadrature_3D")
   use PetIGA
   implicit none
 
 
 pure subroutine IGA_BasisFuns_3D(&
-     order,                 &
-     rational,W,            &
-     inq,ina,ind,iN,        &
-     jnq,jna,jnd,jN,        &
-     knq,kna,knd,kN,        &
-     N0,N1,N2,N3)           &
+     order,                      &
+     rational,W,                 &
+     inq,ina,ind,iN,             &
+     jnq,jna,jnd,jN,             &
+     knq,kna,knd,kN,             &
+     N0,N1,N2,N3)                &
   bind(C, name="IGA_BasisFuns_3D")
   use PetIGA
   implicit none
 
 
 pure subroutine IGA_ShapeFuns_3D(&
-     order,                 &
-     nqp,nen,X,             &
-     M0,M1,M2,M3,           &
-     N0,N1,N2,N3,           &
-     DetF,F,G)              &
+     order,                      &
+     nqp,nen,X,                  &
+     M0,M1,M2,M3,                &
+     N0,N1,N2,N3,                &
+     DetF,F,G)                   &
   bind(C, name="IGA_ShapeFuns_3D")
   use PetIGA
   implicit none
 contains
 include 'petigageo.f90.in'
 end subroutine IGA_ShapeFuns_3D
+
+
+subroutine IGA_BoundaryArea_3D(&
+     m,axis,side,              &
+     geometry,Cx,              &
+     rational,Cw,              &
+     inqp,iW,inen,indr,iN,     &
+     jnqp,jW,jnen,jndr,jN,     &
+     dS)                       &
+  bind(C, name="IGA_BoundaryArea_3D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter         :: nsd = 3
+  integer(kind=IGA_INTEGER_KIND), parameter         :: dim = 2
+  integer(kind=IGA_INTEGER_KIND), intent(in)        :: m(3)
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: axis, side
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: geometry, rational
+  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cx(nsd,m(1),m(2),m(3))
+  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cw(    m(1),m(2),m(3))
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: inqp, inen, indr
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: jnqp, jnen, jndr
+  real   (kind=IGA_REAL_KIND   ), intent(in)        :: iW(inqp), iN(0:indr,inen,inqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)        :: jW(jnqp), jN(0:jndr,jnen,jnqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)       :: dS
+  integer(kind=IGA_INTEGER_KIND)  :: k, nen, iq, jq, ia, ja
+  real   (kind=IGA_REAL_KIND   )  :: N0(inen,jnen), N1(dim,inen,jnen), detJ
+  real   (kind=IGA_REAL_KIND   ), pointer :: Xx(:,:,:), Xw(:,:)
+  nen = inen*jnen
+  select case (axis)
+  case (0)
+     if (side==0) k=1
+     if (side==1) k=m(1)
+     Xx => Cx(:,k,:,:); Xw => Cw(k,:,:)
+  case (1)
+     if (side==0) k=1
+     if (side==1) k=m(2)
+     Xx => Cx(:,:,k,:); Xw => Cw(:,k,:)
+  case (2)
+     if (side==0) k=1
+     if (side==1) k=m(3)
+     Xx => Cx(:,:,:,k); Xw => Cw(:,:,k)
+  end select
+  detJ = 1.0
+  dS = 0.0
+  do jq=1,jnqp
+     do iq=1,inqp
+        forall (ia=1:inen, ja=1:jnen)
+           N0(ia,ja) = iN(0,ia,iq) * jN(0,ja,jq)
+        end forall
+        forall (ia=1:inen, ja=1:jnen)
+           N1(1,ia,ja) = iN(1,ia,iq) * jN(0,ja,jq)
+           N1(2,ia,ja) = iN(0,ia,iq) * jN(1,ja,jq)
+        end forall
+        if (rational /= 0) then
+           call Rationalize(nen,Xw,N0,N1)
+        end if
+        if (geometry /= 0) then
+           call Jacobian(nen,N1,Xx,detJ)
+        end if
+        dS = dS + detJ * (iW(iq)*jW(jq))
+     end do
+  end do
+contains
+pure subroutine Rationalize(nen,W,R0,R1)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter     :: dim = 2
+  integer(kind=IGA_INTEGER_KIND), intent(in)    :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: W(nen)
+  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R0(    nen)
+  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R1(dim,nen)
+  integer(kind=IGA_INTEGER_KIND)  :: i
+  real   (kind=IGA_REAL_KIND   )  :: W0
+  R0 = W * R0
+  W0 = sum(R0)
+  R0 = R0 / W0
+  forall(i=1:dim) &
+  R1(i,:) = W*R1(i,:) - R0 * sum(W*R1(i,:))
+  R1 = R1 / W0
+end subroutine Rationalize
+pure subroutine Jacobian(nen,N,X,J)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: nsd = 3
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: X(nsd,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: J
+  real   (kind=IGA_REAL_KIND   )  :: F(dim,nsd), M(dim,dim)
+  F = matmul(N,transpose(X))
+  M = matmul(F,transpose(F))
+  J = sqrt(abs(Determinant(dim,M)))
+end subroutine Jacobian
+include 'petigainv.f90.in'
+end subroutine IGA_BoundaryArea_3D

File src/petigabound.c

   if (!boundary) PetscFunctionReturn(0);
   PetscValidPointer(boundary,1);
   boundary->dof = 0;
-  boundary->nbc = 0;
+  boundary->count = 0;
   ierr = PetscFree(boundary->field);CHKERRQ(ierr);
   ierr = PetscFree(boundary->value);CHKERRQ(ierr);
+  boundary->nload = 0;
+  ierr = PetscFree(boundary->iload);CHKERRQ(ierr);
+  ierr = PetscFree(boundary->vload);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
   PetscValidPointer(boundary,1);
   ierr = IGABoundaryReset(boundary);CHKERRQ(ierr);
   boundary->dof = dof;
-  boundary->nbc = 0;
+  boundary->count = 0;
   ierr = PetscMalloc1(dof,PetscInt,   &boundary->field);CHKERRQ(ierr);
   ierr = PetscMalloc1(dof,PetscScalar,&boundary->value);CHKERRQ(ierr);
+  boundary->nload = 0;
+  ierr = PetscMalloc1(dof,PetscInt,   &boundary->iload);CHKERRQ(ierr);
+  ierr = PetscMalloc1(dof,PetscScalar,&boundary->vload);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundaryClear"
+PetscErrorCode IGABoundaryClear(IGABoundary boundary)
+{
+  PetscFunctionBegin;
+  PetscValidPointer(boundary,1);
+  boundary->count   = 0;
+  boundary->nload = 0;
   PetscFunctionReturn(0);
 }
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGABoundarySetValue"
 /*@
-   IGABoundarySetValue - Used to set a constant Dirichlet condition on the given boundary.
+   IGABoundarySetValue - Used to set a constant Dirichlet boundary
+   condition on the given boundary.
    
    Logically Collective on IGABoundary
 
   if (field >= dof) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D, but dof %D",field,dof);
   { /**/
     PetscInt pos;
-    for (pos=0; pos<boundary->nbc; pos++)
+    for (pos=0; pos<boundary->count; pos++)
       if (boundary->field[pos] == field) break;
-    if (pos==boundary->nbc) boundary->nbc++;
+    if (pos==boundary->count) boundary->count++;
     boundary->field[pos] = field;
     boundary->value[pos] = value;
   }
   PetscFunctionReturn(0);
 }
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGABoundarySetLoad"
+/*@
+   IGABoundarySetLoad - Used to set a constant Neumann boundary
+   condition on the given boundary.
+   
+   Logically Collective on IGABoundary
+
+   Input Parameters:
++  boundary - the IGAAxis context
+.  field - the index of the field on which to enforce the condition
+-  value - the value to set
+
+   Level: normal
+
+.keywords: IGA, boundary, Neumann
+@*/
+PetscErrorCode IGABoundarySetLoad(IGABoundary boundary,PetscInt field,PetscScalar value)
+{
+  PetscInt dof;
+  PetscFunctionBegin;
+  PetscValidPointer(boundary,1);
+  dof = boundary->dof;
+  if (field <  0)   SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D must be nonnegative",field);
+  if (field >= dof) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D, but dof %D",field,dof);
+  { /**/
+    PetscInt pos;
+    for (pos=0; pos<boundary->nload; pos++)
+      if (boundary->iload[pos] == field) break;
+    if (pos==boundary->nload) boundary->nload++;
+    boundary->iload[pos] = field;
+    boundary->vload[pos] = value;
+  }
+  PetscFunctionReturn(0);
+}

File src/petigaelem.c

   element->nfix = 0;
   ierr = PetscFree(element->ifix);CHKERRQ(ierr);
   ierr = PetscFree(element->vfix);CHKERRQ(ierr);
-  ierr = PetscFree(element->xfix);CHKERRQ(ierr);
+  ierr = PetscFree(element->ufix);CHKERRQ(ierr);
+  element->nflux = 0;
+  ierr = PetscFree(element->iflux);CHKERRQ(ierr);
+  ierr = PetscFree(element->vflux);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
     element->nfix = 0;
     ierr = PetscMalloc1(nen*dof,PetscInt,   &element->ifix);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*dof,PetscScalar,&element->vfix);CHKERRQ(ierr);
-    ierr = PetscMalloc1(nen*dof,PetscScalar,&element->xfix);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nen*dof,PetscScalar,&element->ufix);CHKERRQ(ierr);
+    element->nflux = 0;
+    ierr = PetscMalloc1(nen*dof,PetscInt,   &element->iflux);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nen*dof,PetscScalar,&element->vflux);CHKERRQ(ierr);
   }
   ierr = IGAPointInit(element->iterator,element);CHKERRQ(ierr);
   PetscFunctionReturn(0);
   PetscFunctionReturn(0);
 }
 
-static void AddFix(IGAElement element,PetscInt dir,PetscInt side,PetscInt a)
+/* -------------------------------------------------------------------------- */
+
+static void AddFixa(IGAElement element,IGABoundary b,PetscInt a)
+{
+  if (b->count) {
+    PetscInt dof = element->dof;
+    PetscInt count = element->nfix;
+    PetscInt *index = element->ifix;
+    PetscScalar *value = element->vfix;
+    PetscInt j,k,n = b->count;
+    for (k=0; k<n; k++) {
+      PetscInt idx = a*dof + b->field[k];
+      PetscScalar val = b->value[k];
+      for (j=0; j<count; j++)
+        if (index[j] == idx) break;
+      if (j == count) count++;
+      index[j] = idx;
+      value[j] = val;
+    }
+    element->nfix = count;
+  }
+}
+
+EXTERN_C_BEGIN
+extern void IGA_BoundaryArea_2D(const PetscInt[],PetscInt,PetscInt,
+                                PetscInt,const PetscReal[],
+                                PetscInt,const PetscReal[],
+                                PetscInt,const PetscReal[],PetscInt,PetscInt,const PetscReal[],
+                                PetscReal*);
+extern void IGA_BoundaryArea_3D(const PetscInt[],PetscInt,PetscInt,
+                                PetscInt,const PetscReal[],
+                                PetscInt,const PetscReal[],
+                                PetscInt,const PetscReal[],PetscInt,PetscInt,const PetscReal[],
+                                PetscInt,const PetscReal[],PetscInt,PetscInt,const PetscReal[],
+                                PetscReal*);
+EXTERN_C_END
+
+static PetscReal BoundaryArea(IGAElement element,PetscInt dir,PetscInt side)
+{
+  PetscReal A = 1.0;
+  PetscInt *ID = element->ID;
+  IGABasis *BD = element->BD;
+  PetscInt i,dim = element->dim;
+  if (dim == 1) return A;
+  for (i=0; i<dim; i++)
+    if (i != dir) {
+      PetscReal L = BD[i]->detJ[ID[i]];
+      PetscInt  n = BD[i]->nen;
+      A *= L/n;
+    }
+  if (!element->geometry) {
+    A *= (dim==2) ? 2 : 4; /* sum(W) = 2 */
+  } else {
+    PetscInt shape[3] = {1,1,1};
+    PetscInt k,nqp[3],nen[3],ndr[3];
+    PetscReal *W[3],*N[3],dS = 1.0;
+    for (i=0; i<dim; i++)
+      shape[i] = BD[i]->nen;
+    for (k=0,i=0; i<dim; i++) {
+      if (i == dir) continue;
+      nqp[k] = BD[i]->nqp;
+      nen[k] = BD[i]->nen;
+      ndr[k] = BD[i]->d;
+      W[k]   = BD[i]->weight;
+      N[k]   = BD[i]->value+ID[i]*nqp[k]*nen[k]*(ndr[k]+1);
+      k++;
+    }
+    switch (dim) {
+    case 2: IGA_BoundaryArea_2D(shape,dir,side,
+                                element->geometry,element->geometryX,
+                                element->rational,element->geometryW,
+                                nqp[0],W[0],nen[0],ndr[0],N[0],
+                                &dS); break;
+    case 3: IGA_BoundaryArea_3D(shape,dir,side,
+                                element->geometry,element->geometryX,
+                                element->rational,element->geometryW,
+                                nqp[0],W[0],nen[0],ndr[0],N[0],
+                                nqp[1],W[1],nen[1],ndr[1],N[1],
+                                &dS);break;
+    }
+    A *= dS;
+  }
+  return A;
+}
+
+static void AddFlux(IGAElement element,IGABoundary b,PetscInt a,PetscReal A)
+{
+  if (b->nload) {
+    PetscInt dof = element->dof;
+    PetscInt count = element->nflux;
+    PetscInt *index = element->iflux;
+    PetscScalar *value = element->vflux;
+    PetscInt j,k,n = b->nload;
+    for (k=0; k<n; k++) {
+      PetscInt idx = a*dof + b->iload[k];
+      PetscScalar val = b->vload[k];
+      for (j=0; j<count; j++)
+        if (index[j] == idx) break;
+      if (j == count) value[count++] = 0.0;
+      index[j] = idx;
+      value[j]+= val*A;
+    }
+    element->nflux = count;
+  }
+}
+
+static void BuildFix(IGAElement element,PetscInt dir,PetscInt side)
 {
   IGABoundary b = element->parent->boundary[dir][side];
-  PetscInt dof = element->dof;
-  PetscInt nfix = element->nfix;
-  PetscInt *ifix = element->ifix;
-  PetscScalar *vfix = element->vfix;
-  PetscInt j,k,n = b->nbc;
-  for (k=0; k<n; k++) {
-    PetscInt idx = a*dof + b->field[k];
-    PetscScalar val = b->value[k];
-    for (j=0; j<nfix; j++)
-      if (ifix[j] == idx) break;
-    if (j==nfix) nfix++;
-    ifix[j] = idx;
-    vfix[j] = val;
+  if (b->count || b->nload) {
+    PetscReal Area = b->nload ? BoundaryArea(element,dir,side) : 1.0;
+    IGABasis *BD = element->BD;
+    PetscInt S[3]={0,0,0},E[3]={1,1,1};
+    PetscInt ia,ja,ka,jstride,kstride,a;
+    PetscInt i,dim = element->dim;
+    for (i=0; i<dim; i++) E[i] = BD[i]->nen;
+    jstride = E[0]; kstride = E[0]*E[1];
+    if (side) S[dir] = E[dir]-1;
+    else      E[dir] = S[dir]+1;
+    for (ka=S[2]; ka<E[2]; ka++)
+      for (ja=S[1]; ja<E[1]; ja++)
+        for (ia=S[0]; ia<E[0]; ia++)
+          {
+            a = ia + ja*jstride + ka*kstride;
+            AddFixa(element,b,a);
+            AddFlux(element,b,a,Area);
+          }
   }
-  element->nfix = nfix;
 }
 
 #undef  __FUNCT__
   PetscValidPointer(element,1);
   if (PetscUnlikely(element->index < 0))
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call during element loop");
-  {
+  element->nfix  = 0;
+  element->nflux = 0;
+  if (!element->collocation) {
+    IGAAxis  *AX = element->parent->axis;
     IGABasis *BD = element->BD;
     PetscInt *ID = element->ID;
+    PetscInt i,dim = element->dim;
+    for (i=0; i<dim; i++) {
+      PetscBool w = AX[i]->periodic;
+      PetscInt  e = BD[i]->nel-1; /* last element */
+      if (ID[i] == 0 && !w) BuildFix(element,i,0);
+      if (ID[i] == e && !w) BuildFix(element,i,1);
+    }
+  } else {
     PetscInt A0[3] = {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT};
     PetscInt A1[3] = {PETSC_MAX_INT,PETSC_MAX_INT,PETSC_MAX_INT};
-    PetscBool onboundary = PETSC_FALSE;
-    PetscInt i,dim = element->dim;
-    for (i=0; i<dim; i++) {
-      PetscInt  e = BD[i]->nel-1; /* last element */
-      PetscInt  n = BD[i]->nnp-1; /* last node    */
-      PetscBool b = element->collocation;
-      if (element->parent->axis[i]->periodic) continue; /* XXX */
-      if (ID[i] == 0 || b) { A0[i] = 0; onboundary = PETSC_TRUE; }
-      if (ID[i] == e || b) { A1[i] = n; onboundary = PETSC_TRUE; }
+    {
+      IGAAxis  *AX = element->parent->axis;
+      PetscInt i,dim = element->dim;
+      for (i=0; i<dim; i++) {
+        PetscBool w = AX[i]->periodic;
+        PetscInt  n = AX[i]->nnp-1; /* last node */
+        A0[i] = 0; if (!w) A1[i] = n;
+      }
     }
-    element->nfix = 0;
-    if (onboundary) {
+    {
+      IGABoundary (*b)[2] = element->parent->boundary;
+      IGABasis *BD = element->BD;
+      PetscInt *ID = element->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 a = 0;
-      for (ka=0; ka<knen; ka++) {
-        for (ja=0; ja<jnen; ja++) {
-          for (ia=0; ia<inen; ia++) {
-            PetscInt iA = ioffset + ia;
-            PetscInt jA = joffset + ja;
-            PetscInt kA = koffset + ka;
-            /**/ if (iA == A0[0]) AddFix(element,0,0,a);
-            else if (iA == A1[0]) AddFix(element,0,1,a);
-            /**/ if (jA == A0[1]) AddFix(element,1,0,a);
-            else if (jA == A1[1]) AddFix(element,1,1,a);
-            /**/ if (kA == A0[2]) AddFix(element,2,0,a);
-            else if (kA == A1[2]) AddFix(element,2,1,a);
-            a++;
-          }
-        }
-      }
+      for (ka=0; ka<knen; ka++)
+        for (ja=0; ja<jnen; ja++)
+          for (ia=0; ia<inen; ia++)
+            {
+              PetscInt iA = ioffset + ia;
+              PetscInt jA = joffset + ja;
+              PetscInt kA = koffset + ka;
+              /**/ if (iA == A0[0]) AddFixa(element,b[0][0],a);
+              else if (iA == A1[0]) AddFixa(element,b[0][1],a);
+              /**/ if (jA == A0[1]) AddFixa(element,b[1][0],a);
+              else if (jA == A1[1]) AddFixa(element,b[1][1],a);
+              /**/ if (kA == A0[2]) AddFixa(element,b[2][0],a);
+              else if (kA == A1[2]) AddFixa(element,b[2][1],a);
+              a++;
+            }
     }
   }
   PetscFunctionReturn(0);
   PetscValidPointer(element,1);
   PetscValidScalarPointer(U,2);
   {
-    PetscInt f,nfix=element->nfix;
-    for (f=0; f<nfix; f++) {
+    PetscInt f,n;
+    n = element->nfix;
+    for (f=0; f<n; f++) {
       PetscInt k = element->ifix[f];
-      element->xfix[f] = U[k];
+      element->ufix[f] = U[k];
       U[k] = element->vfix[f];
     }
   }
   if (!element->collocation) {
     PetscInt M = element->neq * element->dof;
     PetscInt N = element->nen * element->dof;
-    PetscInt f,nfix=element->nfix;
-    for (f=0; f<nfix; f++) {
+    PetscInt f,n;
+    n = element->nflux;
+    for (f=0; f<n; f++) {
+      PetscInt    k = element->iflux[f];
+      PetscScalar v = element->vflux[f];
+      F[k] += v;
+    }
+    n = element->nfix;
+    for (f=0; f<n; f++) {
       PetscInt    k = element->ifix[f];
       PetscScalar v = element->vfix[f];
       PetscInt i,j;
-      for (i=0; i<M; i++)
-        F[i] -= K[i*N+k] * v;
+      for (i=0; i<M; i++) F[i] -= K[i*N+k] * v;
       for (i=0; i<M; i++) K[i*N+k] = 0.0;
       for (j=0; j<N; j++) K[k*N+j] = 0.0;
       K[k*N+k] = 1.0;
   } else {
     PetscInt M = element->neq * element->dof;
     PetscInt N = element->nen * element->dof;
-    PetscInt f,nfix=element->nfix;
-    for (f=0; f<nfix; f++) {
+    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]) {
   PetscValidPointer(element,1);
   PetscValidScalarPointer(F,2);
   if (!element->collocation) {
-    PetscInt f,nfix=element->nfix;
-    for (f=0; f<nfix; f++) {
+    PetscInt f,n;
+    n = element->nflux;
+    for (f=0; f<n; f++) {
+      PetscInt    k = element->iflux[f];
+      PetscScalar v = element->vflux[f];
+      F[k] -= v;
+    }
+    n = element->nfix;
+    for (f=0; f<n; f++) {
       PetscInt    k = element->ifix[f];
-      PetscScalar u = element->xfix[f];
       PetscScalar v = element->vfix[f];
+      PetscScalar u = element->ufix[f];
       F[k] = u - v;
     }
   } else {
-    PetscInt f,nfix=element->nfix;
-    for (f=0; f<nfix; f++) {
+    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 u = element->xfix[f];
         PetscScalar v = element->vfix[f];
+        PetscScalar u = element->ufix[f];
         F[c] = u - v;
       }
     }
   if (!element->collocation) {
     PetscInt M = element->neq * element->dof;
     PetscInt N = element->nen * element->dof;
-    PetscInt f,nfix=element->nfix;
-    for (f=0; f<nfix; f++) {
+    PetscInt f,n;
+    n = element->nfix;
+    for (f=0; f<n; f++) {
       PetscInt k = element->ifix[f];
       PetscInt i,j;
       for (i=0; i<M; i++) J[i*N+k] = 0.0;
   } else {
     PetscInt M = element->neq * element->dof;
     PetscInt N = element->nen * element->dof;
-    PetscInt f,nfix=element->nfix;
-    for (f=0; f<nfix; f++) {
+    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]) {