Lisandro Dalcin avatar Lisandro Dalcin committed ccbe5f2

Add support for different basis types (with Gabriel Espinosa)

Comments (0)

Files changed (4)

 PETSC_EXTERN PetscErrorCode IGARuleSetRule(IGARule rule,PetscInt q,const PetscReal x[],const PetscReal w[]);
 PETSC_EXTERN PetscErrorCode IGARuleGetRule(IGARule rule,PetscInt *q,PetscReal *x[],PetscReal *w[]);
 
+typedef enum {
+  IGA_BASIS_BSPLINE=0,
+  IGA_BASIS_BERNSTEIN,
+  IGA_BASIS_LAGRANGE
+} IGABasisType;
+
+PETSC_EXTERN const char *const IGABasisTypes[];
+
 struct _n_IGABasis {
   PetscInt refct;
   /**/
+  IGABasisType type;  /* basis type */
+  /**/
   PetscInt  nel;      /* number of elements */
   PetscInt  nqp;      /* number of quadrature points */
   PetscInt  nen;      /* number of local basis functions */
 PETSC_EXTERN PetscErrorCode IGABasisDestroy(IGABasis *basis);
 PETSC_EXTERN PetscErrorCode IGABasisReset(IGABasis basis);
 PETSC_EXTERN PetscErrorCode IGABasisReference(IGABasis basis);
+PETSC_EXTERN PetscErrorCode IGABasisSetType(IGABasis basis,IGABasisType type);
 PETSC_EXTERN PetscErrorCode IGABasisInitQuadrature (IGABasis basis,IGAAxis axis,IGARule rule,PetscInt order);
 PETSC_EXTERN PetscErrorCode IGABasisInitCollocation(IGABasis basis,IGAAxis axis,PetscInt order);
 
   {
     PetscBool flg;
     PetscInt  i,nw,nb;
+    IGABasisType btype[3] = {IGA_BASIS_BSPLINE,IGA_BASIS_BSPLINE,IGA_BASIS_BSPLINE};
     PetscBool wraps[3]    = {PETSC_FALSE,PETSC_FALSE,PETSC_FALSE};
     PetscInt  np,procs[3] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
     PetscInt  nd,degrs[3] = {2,2,2};
     for (i=0; i<dim; i++) wraps[i] = iga->axis[i]->periodic;
     for (i=0; i<dim; i++) if (iga->axis[i]->p   > 0) degrs[i] = iga->axis[i]->p;
     for (i=0; i<dim; i++) if (iga->rule[i]->nqp > 0) quadr[i] = iga->rule[i]->nqp;
+    for (i=0; i<dim; i++) btype[i] = iga->basis[i]->type;
 
     ierr = PetscObjectOptionsBegin((PetscObject)iga);CHKERRQ(ierr);
 
 	ierr = IGAAxisSetPeriodic(iga->axis[i],w);CHKERRQ(ierr);
       }
 
+    /* Basis */
+    ierr = PetscOptionsEnum("-iga_basis_type","Basis type","IGABasisSetType",IGABasisTypes,(PetscEnum)btype[0],(PetscEnum*)&btype[0],&flg);CHKERRQ(ierr);
+    for (i=0; i<dim; i++) btype[i] = btype[0]; /* XXX */
+    if (flg) for (i=0; i<dim; i++) {
+        ierr = IGABasisSetType(iga->basis[i],btype[i]);CHKERRQ(ierr);
+      }
+    for (i=0; i<dim; i++) if (btype[i] != IGA_BASIS_BSPLINE) conts[i] = 0;
+
     /* Geometry */
     ierr = PetscOptionsString("-iga_geometry","Specify IGA geometry file","IGARead",filename,filename,sizeof(filename),&flg);CHKERRQ(ierr);
     if (flg) { /* load from file */
     if (flg) { ierr = IGASetOrder(iga,order);CHKERRQ(ierr);}
 
     /* Quadrature */
-    ierr = PetscOptionsIntArray ("-iga_quadrature","Quadrature points","IGARuleInit",quadr,(nq=dim,&nq),&flg);CHKERRQ(ierr);
+    ierr = PetscOptionsIntArray("-iga_quadrature","Quadrature points","IGARuleInit",quadr,(nq=dim,&nq),&flg);CHKERRQ(ierr);
     if (flg) for (i=0; i<dim; i++) {
 	PetscInt q = (i<nq) ? quadr[i] : quadr[0];
 	if (q > 0) {ierr = IGARuleInit(iga->rule[i],q);CHKERRQ(ierr);}

src/petigabasis.c

 #include "petiga.h"
 
+const char *const IGABasisTypes[] = {
+  "BSPLINE",
+  "BERNSTEIN",
+  "LAGRANGE",
+  /* */
+  "IGABasisType","IGA_BASIS_",0};
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGABasisCreate"
 PetscErrorCode IGABasisCreate(IGABasis *basis)
   PetscFunctionReturn(0);
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "IGABasisSetType"
+PetscErrorCode IGABasisSetType(IGABasis basis,IGABasisType type)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(basis,1);
+  if (basis->type != type) {ierr = IGABasisReset(basis);CHKERRQ(ierr);}
+  basis->type = type;
+  PetscFunctionReturn(0);
+}
+
 EXTERN_C_BEGIN
-extern void IGA_Basis_BSpline(PetscInt i,PetscReal u,PetscInt p,PetscInt d,const PetscReal U[],PetscReal N[]);
+extern void IGA_Basis_BSpline (PetscInt i,PetscReal u,PetscInt p,PetscInt d,const PetscReal U[],PetscReal B[]);
+extern void IGA_Basis_Lagrange(PetscInt i,PetscReal u,PetscInt p,PetscInt d,const PetscReal U[],PetscReal L[]);
 EXTERN_C_END
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGABasisInitQuadrature"
 PetscErrorCode IGABasisInitQuadrature(IGABasis basis,IGAAxis axis,IGARule rule,PetscInt d)
 {
-  PetscInt       p;
+  PetscInt       m,p;
   const PetscInt *span;
   const PetscReal*U,*X,*W;
   PetscInt       iel,nel;
   PetscReal      *weight;
   PetscReal      *point;
   PetscReal      *value;
+  void          (*ComputeBasis)(PetscInt,PetscReal,PetscInt,PetscInt,const PetscReal[],PetscReal[]);
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(basis,1);
   if (d < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
                       "Derivative order must be grather than zero, got %D",d);
 
+  m = axis->m;
   p = axis->p;
   U = axis->U;
 
+  if (basis->type != IGA_BASIS_BSPLINE) {
+    PetscInt s,j,k=1;
+    while (k < m) {
+      j = k; s = 1; while (++k < m && U[j] == U[k]) s++;
+      if (s < p) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
+                          "Basis type %s requires C^0 continuity, "
+                          "Knot U[%D]=%G has multiplicity %D "
+                          "less than polynomial degree %D",
+                          IGABasisTypes[basis->type],j,U[j],s,p);
+    }
+    /* XXX */printf("using %s basis\n",IGABasisTypes[basis->type]);
+  }
+
   nqp = rule->nqp;
   X   = rule->point;
   W   = rule->weight;
   nen  = p+1;
   ndr  = d+1;
 
+  switch (basis->type) {
+  case IGA_BASIS_BSPLINE:
+  case IGA_BASIS_BERNSTEIN:
+    ComputeBasis = IGA_Basis_BSpline; break;
+  case IGA_BASIS_LAGRANGE:
+    ComputeBasis = IGA_Basis_Lagrange; break;
+  }
+
   ierr = PetscMalloc1(nel,PetscInt,&offset);CHKERRQ(ierr);
   ierr = PetscMalloc1(nel,PetscReal,&detJ);CHKERRQ(ierr);
   ierr = PetscMalloc1(nqp,PetscReal,&weight);CHKERRQ(ierr);
   ierr = PetscMalloc1(nel*nqp,PetscReal,&point);CHKERRQ(ierr);
   ierr = PetscMalloc1(nel*nqp*nen*ndr,PetscReal,&value);CHKERRQ(ierr);
+  ierr = PetscMemzero(value,nel*nqp*nen*ndr*sizeof(PetscReal));CHKERRQ(ierr);
 
   for (iqp=0; iqp<nqp; iqp++) {
     weight[iqp] = W[iqp];
     detJ[iel] = J;
     for (iqp=0; iqp<nqp; iqp++) {
       u[iqp] = (X[iqp] + 1) * J + u0;
-      IGA_Basis_BSpline(k,u[iqp],p,d,U,&N[iqp*nen*ndr]);
+      ComputeBasis(k,u[iqp],p,d,U,&N[iqp*nen*ndr]);
     }
     offset[iel] = k-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]);
+    ComputeBasis(k0,u0,p,d,U,basis->bnd_value[0]);
+    ComputeBasis(k1,u1,p,d,U,basis->bnd_value[1]);
   }
 
   PetscFunctionReturn(0);

src/petigabsp.f90

 ! -*- f90 -*-
 
-subroutine IGA_Basis_BSpline(i,uu,p,d,U,N) &
+subroutine IGA_Basis_BSpline(k,uu,p,d,U,B) &
   bind(C, name="IGA_Basis_BSpline")
   use PetIGA
   implicit none
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: i, p, d
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: k, p, d
   real   (kind=IGA_REAL_KIND   ), intent(in),value :: uu
-  real   (kind=IGA_REAL_KIND   ), intent(in)       :: U(0:i+p)
-  real   (kind=IGA_REAL_KIND   ), intent(out)      :: N(0:d,0:p)
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: U(0:k+p)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: B(0:d,0:p)
   real   (kind=IGA_REAL_KIND   )  :: ders(0:p,0:d)
-  call BasisFunsDers(i,uu,p,d,U,ders)
-  N = transpose(ders)
+  call BasisFunsDers(k,uu,p,d,U,ders)
+  B = transpose(ders)
 contains
 pure subroutine BasisFunsDers(i,uu,p,n,U,ders)
   use PetIGA
   end do
 end subroutine BasisFunsDers
 end subroutine IGA_Basis_BSpline
+
+
+subroutine IGA_Basis_Lagrange(kk,uu,p,d,U,B) &
+  bind(C, name="IGA_Basis_Lagrange")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: kk, p, d
+  real   (kind=IGA_REAL_KIND   ), intent(in),value :: uu
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: U(0:kk+p)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: B(0:d,0:p)
+  integer(kind=IGA_INTEGER_KIND)  :: m, i, j, k, l
+  real   (kind=IGA_REAL_KIND   )  :: Lp, Ls1, Ls2, Ls3
+  real   (kind=IGA_REAL_KIND   )  :: X(0:p)
+
+  forall (m=0:p) X(m) = U(kk) + m * (U(kk+1) - U(kk)) / p
+
+  do m = 0, p
+     Lp = 1
+     do i = 0, p
+        if (i == m) cycle
+        Lp = Lp * (uu-X(i))/(X(m)-X(i))
+     end do
+     B(0,m) = Lp
+  end do
+
+  if (d < 1) return
+  do m = 0, p
+     Ls1 = 0
+     do j = 0, p
+        if (j == m) cycle
+        Lp = 1
+        do i = 0, p
+           if (i == m .or. i == j) cycle
+           Lp = Lp * (uu-X(i))/(X(m)-X(i))
+        end do
+        Ls1 = Ls1 + Lp/(X(m)-X(j))
+     end do
+     B(1,m) = Ls1
+  end do
+
+  if (d < 2) return
+  do m = 0, p
+     Ls2 = 0
+     do k = 0, p
+        if (k == m) cycle
+        Ls1 = 0
+        do j = 0, p
+           if (j == m .or. j == k) cycle
+           Lp = 1
+           do i = 0, p
+              if (i == m .or. i == k .or. i == j) cycle
+              Lp = Lp * (uu-X(i))/(X(m)-X(i))
+           end do
+           Ls1 = Ls1 + Lp/(X(m)-X(j))
+        end do
+        Ls2 = Ls2 + Ls1/(X(m)-X(k))
+     end do
+     B(2,m) = Ls2
+  end do
+
+
+  if (d < 3) return
+  do m = 0, p
+     Ls3 = 0
+     do l = 0, p
+        if (l == m) cycle
+        Ls2 = 0
+        do k = 0, p
+           if (k == m .or. k == l) cycle
+           Ls1 = 0
+           do j = 0, p
+              if (j == m .or. j == l .or. j == k) cycle
+              Lp = 1
+              do i = 0, p
+                 if (i == m .or. i == l .or. i == k .or. i == j) cycle
+                 Lp = Lp * (uu-X(i))/(X(m)-X(i))
+              end do
+              Ls1 = Ls1 + Lp/(X(m)-X(j))
+           end do
+           Ls2 = Ls2 + Ls1/(X(m)-X(k))
+        end do
+        Ls3 = Ls3 + Ls2/(X(m)-X(l))
+     end do
+     B(3,m) = Ls3
+  end do
+
+end subroutine IGA_Basis_Lagrange
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.