Source

PetIGA / src / petigabasis.c

Full commit
Lisandro Dalcin f3d7eb8 

Lisandro Dalcin ccbe5f2 






Lisandro Dalcin f3d7eb8 






















Lisandro Dalcin 67ed8c1 
Lisandro Dalcin f3d7eb8 




Lisandro Dalcin 67ed8c1 











Lisandro Dalcin 7dc1b6c 
Lisandro Dalcin 67ed8c1 



Lisandro Dalcin 82de52a 

Lisandro Dalcin 67ed8c1 



Lisandro Dalcin f3d7eb8 








Lisandro Dalcin ccbe5f2 











Lisandro Dalcin f3d7eb8 
Lisandro Dalcin ccbe5f2 

Lisandro Dalcin f3d7eb8 


Lisandro Dalcin f0d255a 

Lisandro Dalcin f3d7eb8 
Lisandro Dalcin ccbe5f2 
Lisandro Dalcin af24fc2 

Lisandro Dalcin f3d7eb8 


Lisandro Dalcin 7dc1b6c 
Lisandro Dalcin f3d7eb8 



Lisandro Dalcin ccbe5f2 
Lisandro Dalcin f3d7eb8 




Lisandro Dalcin 790d381 

Lisandro Dalcin f3d7eb8 
Lisandro Dalcin ccbe5f2 
Lisandro Dalcin f3d7eb8 


Lisandro Dalcin ccbe5f2 












Lisandro Dalcin f3d7eb8 



Lisandro Dalcin af24fc2 



Lisandro Dalcin dd7d956 
Lisandro Dalcin ccbe5f2 







Lisandro Dalcin 7dc1b6c 
Lisandro Dalcin f3d7eb8 



Lisandro Dalcin ccbe5f2 
Lisandro Dalcin 7dc1b6c 
Lisandro Dalcin f3d7eb8 











Lisandro Dalcin ccbe5f2 
Lisandro Dalcin f3d7eb8 
Lisandro Dalcin af24fc2 
Lisandro Dalcin f3d7eb8 

Lisandro Dalcin 67ed8c1 
Lisandro Dalcin f3d7eb8 





Lisandro Dalcin 7dc1b6c 

Lisandro Dalcin f3d7eb8 



Lisandro Dalcin 67ed8c1 
Lisandro Dalcin 82de52a 
Lisandro Dalcin 8f8bb27 


Lisandro Dalcin 82de52a 

Lisandro Dalcin 8f8bb27 
Lisandro Dalcin 82de52a 

Lisandro Dalcin 8f8bb27 
Lisandro Dalcin ccbe5f2 

Lisandro Dalcin 82de52a 

Lisandro Dalcin f3d7eb8 

Lisandro Dalcin f0d255a 









Lisandro Dalcin 05b4b56 
Lisandro Dalcin f0d255a 








































Lisandro Dalcin 87105e5 
Lisandro Dalcin f0d255a 











































#include "petiga.h"

const char *const IGABasisTypes[] = {
  "BSPLINE",
  "BERNSTEIN",
  "LAGRANGE",
  /* */
  "IGABasisType","IGA_BASIS_",0};

#undef  __FUNCT__
#define __FUNCT__ "IGABasisCreate"
PetscErrorCode IGABasisCreate(IGABasis *basis)
{
  PetscErrorCode ierr;
  PetscFunctionBegin;
  PetscValidPointer(basis,3);
  ierr = PetscNew(struct _n_IGABasis,basis);CHKERRQ(ierr);
  (*basis)->refct = 1;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "IGABasisDestroy"
PetscErrorCode IGABasisDestroy(IGABasis *_basis)
{
  IGABasis       basis;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  PetscValidPointer(_basis,1);
  basis = *_basis; *_basis = 0;
  if (!basis) PetscFunctionReturn(0);
  if (--basis->refct > 0) PetscFunctionReturn(0);
  ierr = IGABasisReset(basis);CHKERRQ(ierr);
  ierr = PetscFree(basis);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "IGABasisReset"
PetscErrorCode IGABasisReset(IGABasis basis)
{
  PetscErrorCode ierr;
  PetscFunctionBegin;
  if (!basis) PetscFunctionReturn(0);
  PetscValidPointer(basis,1);
  basis->nel = 0;
  basis->nqp = 0;
  basis->nen = 0;
  basis->p   = 0;
  basis->d   = 0;
  ierr = PetscFree(basis->offset);CHKERRQ(ierr);
  ierr = PetscFree(basis->detJ);CHKERRQ(ierr);
  ierr = PetscFree(basis->weight);CHKERRQ(ierr);
  ierr = PetscFree(basis->point);CHKERRQ(ierr);
  ierr = PetscFree(basis->value);CHKERRQ(ierr);
  ierr = PetscFree(basis->bnd_value[0]);CHKERRQ(ierr);
  ierr = PetscFree(basis->bnd_value[1]);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "IGABasisReference"
PetscErrorCode IGABasisReference(IGABasis basis)
{
  PetscFunctionBegin;
  PetscValidPointer(basis,1);
  basis->refct++;
  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 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       m,p;
  const PetscInt *span;
  const PetscReal*U,*X,*W;
  PetscInt       iel,nel;
  PetscInt       iqp,nqp;
  PetscInt       nen,ndr;
  PetscInt       *offset;
  PetscReal      *detJ;
  PetscReal      *weight;
  PetscReal      *point;
  PetscReal      *value;
  void          (*ComputeBasis)(PetscInt,PetscReal,PetscInt,PetscInt,const PetscReal[],PetscReal[]);
  PetscErrorCode ierr;
  PetscFunctionBegin;
  PetscValidPointer(basis,1);
  PetscValidPointer(axis,2);
  PetscValidPointer(rule,3);
  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;

  nel  = axis->nel;
  span = axis->span;
  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];
  }
  for (iel=0; iel<nel; iel++) {
    PetscInt  k = span[iel];
    PetscReal u0 = U[k], u1 = U[k+1];
    PetscReal J = (u1-u0)/2;
    PetscReal *u = &point[iel*nqp];
    PetscReal *N = &value[iel*nqp*nen*ndr];
    detJ[iel] = J;
    for (iqp=0; iqp<nqp; iqp++) {
      u[iqp] = (X[iqp] + 1) * J + u0;
      ComputeBasis(k,u[iqp],p,d,U,&N[iqp*nen*ndr]);
    }
    offset[iel] = k-p;
  }

  ierr = IGABasisReset(basis);CHKERRQ(ierr);

  basis->nel    = nel;
  basis->nqp    = nqp;
  basis->nen    = nen;
  basis->p      = p;
  basis->d      = d;
  basis->offset = offset;

  basis->detJ   = detJ;
  basis->weight = weight;
  basis->point  = point;
  basis->value  = value;

  {
    PetscInt  o0 = offset[0], o1 = offset[nel-1];
    PetscInt  k0 = span[0],   k1 = span[nel-1];
    PetscReal u0 = U[k0],     u1 = U[k1+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] =  o0; basis->bnd_offset[1] =  o1;
    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;
    ComputeBasis(k0,u0,p,d,U,basis->bnd_value[0]);
    ComputeBasis(k1,u1,p,d,U,basis->bnd_value[1]);
  }

  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
extern PetscInt  IGA_FindSpan(PetscInt n,PetscInt p,PetscReal u, const PetscReal *U);
extern PetscReal IGA_Greville(PetscInt i,PetscInt p,const PetscReal U[]);
EXTERN_C_END

#undef  __FUNCT__
#define __FUNCT__ "IGABasisInitCollocation"
PetscErrorCode IGABasisInitCollocation(IGABasis basis,IGAAxis axis,PetscInt d)
{
  PetscInt       p,n;
  const PetscReal*U;
  PetscInt       iel,nel;
  PetscInt       iqp,nqp;
  PetscInt       nen,ndr;
  PetscInt       *offset;
  PetscReal      *detJ;
  PetscReal      *weight;
  PetscReal      *point;
  PetscReal      *value;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  PetscValidPointer(basis,1);
  PetscValidPointer(axis,2);
  if (d < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
                      "Derivative order must be grather than zero, got %D",d);

  p = axis->p;
  n = axis->m - p -1;
  U = axis->U;

  nel  = axis->nnp;
  nqp  = 1;
  nen  = p+1;
  ndr  = d+1;

  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);

  for (iqp=0; iqp<nqp; iqp++) {
    weight[iqp] = 1.0;
  }
  for (iel=0; iel<nel; iel++) {
    PetscReal u = IGA_Greville(iel,p,U);
    PetscInt  k = IGA_FindSpan(n,p,u,U);
    PetscReal *N = &value[iel*nen*ndr];
    offset[iel] = k-p;
    point[iel]  = u;
    detJ[iel]   = U[k+1]-U[k];
    IGA_Basis_BSpline(k,u,p,d,U,N);
  }

  ierr = IGABasisReset(basis);CHKERRQ(ierr);

  basis->nel    = nel;
  basis->nqp    = nqp;
  basis->nen    = nen;
  basis->p      = p;
  basis->d      = d;
  basis->offset = offset;

  basis->detJ   = detJ;
  basis->weight = weight;
  basis->point  = point;
  basis->value  = value;

  PetscFunctionReturn(0);
}

PetscInt IGA_FindSpan(PetscInt n,PetscInt p,PetscReal u, const PetscReal U[])
{
  PetscInt low,high,span;
  if(u >= U[n+1]) return n;
  if(u <= U[p])   return p;
  low  = p;
  high = n+1;
  span = (high+low)/2;
  while(u < U[span] || u >= U[span+1]){
    if(u < U[span])
      high = span;
    else
      low = span;
    span = (high+low)/2;
  }
  return span;
}

PetscReal IGA_Greville(PetscInt i,PetscInt p,const PetscReal U[])
{
  PetscInt j;
  PetscReal u = 0.0;
  for(j=0;j<p;j++) u += U[i+j+1];
  u /= p;
  return u;
}