Commits

Nathan Collier  committed c996fbd

added a crude and simple adaptive L2 projection

  • Participants
  • Parent commits a0e316e

Comments (0)

Files changed (2)

File demo/AdaptiveL2Projection.c

+#include "petiga.h"
+
+#define SQ(A) (A)*(A)
+
+PetscScalar Function(PetscReal x, PetscReal y, PetscReal z)
+{
+  PetscScalar x0=0.85,y0=0.45;
+  PetscScalar sigx=0.125,sigy=0.2;
+  return exp( -( SQ(x-x0) / (2.*SQ(sigx)) + SQ(y-y0) / (2.*SQ(sigy) ) ) );
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt nen;
+  IGAPointGetSizes(p,0,&nen,0);
+  
+  PetscReal x[3] = {0,0,0};
+  IGAPointFormPoint(p,x);
+  PetscScalar f = Function(x[0],x[1],x[2]);
+  
+  const PetscReal *N;
+  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N);
+  
+  PetscInt a,b;
+  for (a=0; a<nen; a++) {
+    for (b=0; b<nen; b++) {
+      K[a*nen+b] = N[a] * N[b];
+    }
+    F[a] = N[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 f = Function(x[0],x[1],x[2]);
+  
+  PetscScalar u;
+  IGAPointFormValue(p,U,&u);
+  
+  PetscReal e = PetscAbsScalar(u - f);
+  S[0] = e*e;
+  S[1] = f;
+
+  return 0;
+}
+
+#define MAX_BREAKS 1000
+
+#undef  __FUNCT__
+#define __FUNCT__ "ComputeError"
+PetscErrorCode ComputeError(PetscInt dim,PetscInt p,PetscInt C,PetscReal (*U)[MAX_BREAKS],PetscInt *N,PetscReal *error)
+{
+  
+  PetscErrorCode  ierr;
+  PetscInt i;
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  for (i=0; i<dim; i++) {
+    IGAAxis axis;
+    ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
+    ierr = IGAAxisSetDegree(axis,p);CHKERRQ(ierr);
+    ierr = IGAAxisInitBreaks(axis,N[i],U[i],C);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 = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+
+  PetscReal errors[2];
+  ierr  = IGAFormScalar(iga,x,2,&errors[0],Error,PETSC_NULL);CHKERRQ(ierr);
+  errors[0] = PetscSqrtReal(PetscRealPart(errors[0]));
+  errors[1] = PetscSqrtReal(PetscRealPart(errors[1]));
+  *error = errors[0]/errors[1];
+
+  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);
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "FindBestKnotToInsert"
+PetscErrorCode FindBestKnotToInsert(PetscInt dim,PetscInt idim,PetscInt p,PetscInt C,PetscReal (*U)[MAX_BREAKS],
+				    PetscInt *N,PetscInt *insert,PetscReal *min_error) 
+{
+  PetscErrorCode  ierr;
+  PetscInt i,j;
+
+  // Initialize other dimensions to values in U
+  PetscReal Ur[3][MAX_BREAKS];
+  PetscInt  Nr[3];
+  for(i=0;i<dim;i++) {
+    if(i != idim) for(j=0;j<N[i];j++) Ur[i][j] = U[i][j];
+    Nr[i] = N[i];
+  }
+  Nr[idim] += 1;
+
+  *insert    = -1;
+  *min_error = 1e20;
+
+  // Try refining each nonzero span
+  for(i=0;i<N[idim]-1;i++){  
+    for(j=0;j<N[idim]+1;j++){
+      if(j<=i)   Ur[idim][j] = U[idim][j]; 
+      if(j==i+1) Ur[idim][j] = 0.5*(U[idim][j-1]+U[idim][j]);
+      if(j>i+1)  Ur[idim][j] = U[idim][j-1]; 
+    }
+    PetscReal error;
+    ierr = ComputeError(dim,p,C,&Ur[0],&Nr[0],&error);CHKERRQ(ierr);
+    if(error < *min_error){
+      *min_error = error;
+      *insert    = i;
+    }
+  }
+
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  PetscReal U[3][MAX_BREAKS];
+  PetscInt  N[3];
+
+  PetscInt i,j;
+  PetscInt dim = 2, p = 2, C = 1;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-p","polynomial order",       __FILE__,p,&p,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-C","global continuity order",__FILE__,C,&C,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+
+  // initial breaks (2x2 mesh)
+  for (i=0; i<dim; i++) { U[i][0] = 0.; U[i][1] = .5; U[i][2] = 1.; N[i] = 3; }
+
+  PetscInt step = 0;
+  PetscReal min_error = 1.e20;
+  while(step<50 && min_error > 1e-9){
+
+    PetscInt insert[3],imin=-1;
+    PetscReal error[3];
+    min_error = 1.e20;
+
+    for(i=0;i<dim;i++){
+      ierr = FindBestKnotToInsert(dim,i,p,C,&U[0],&N[0],&insert[i],&error[i]);CHKERRQ(ierr);
+      if(error[i] < min_error){
+	min_error = error[i];
+	imin = i;
+      }
+    }
+
+    if(imin < 0) {
+      SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"No minimum found.");
+    }
+
+    PetscScalar in = 0.5*(U[imin][insert[imin]]+U[imin][insert[imin]+1]);
+    for(j=N[imin];j>insert[imin];j--) U[imin][j] = U[imin][j-1];
+    U[imin][insert[imin]+1] = in;
+    N[imin] += 1;
+
+    printf("Step %d: Error %.16e\n",step,min_error);
+
+    step += 1;
+  }
+  
+  
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  for (i=0; i<dim; i++) {
+    IGAAxis axis;
+    ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
+    ierr = IGAAxisSetDegree(axis,p);CHKERRQ(ierr);
+    ierr = IGAAxisInitBreaks(axis,N[i],U[i],C);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 = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+  
+  ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
+
+  ierr = IGAWrite(iga,"iga.dat");CHKERRQ(ierr);
+  ierr = IGAWriteVec(iga,x,"sol.dat");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

 	  L2Proj1D L2Proj2D \
 	  Poisson1D Poisson2D Poisson3D \
 	  L2Projection \
+	  AdaptiveL2Projection \
 	  AdvectionDiffusion \
 	  Laplace \
 	  Poisson \
 L2Projection: L2Projection.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
 	${RM} -f $<
+AdaptiveL2Projection: AdaptiveL2Projection.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
 Laplace: Laplace.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
 	${RM} -f $<