# PetIGA

committed ee7de47

Code for Anders Logg's challenge in Google+

# demo/LoggChallenge.c

`+#include "petiga.h"`
`+`
`+/*`
`+`
`+Solve the partial diferential equation`
`+`
`+   -Laplacian(u) = f`
`+`
`+with homogeneous Dirichlet boundary conditions on the unit square for`
`+`
`+   f(x,y) = 2 pi^2 sin(pi*x) * sin(pi*y).`
`+`
`++ Who can obtain the smallest error?`
`++ Who can compute a solution with an error smaller than 10^-6?`
`+`
`+*/`
`+`
`+PetscReal Forcing(PetscReal x, PetscReal y)`
`+{`
`+  PetscReal pi = M_PI;`
`+  return 2*pi*pi * sin(pi*x) * sin(pi*y);`
`+}`
`+`
`+#undef  __FUNCT__`
`+#define __FUNCT__ "System"`
`+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)`
`+{`
`+  PetscReal x = p->point[0];`
`+  PetscReal y = p->point[1];`
`+  PetscReal *N0 = p->shape[0];`
`+  PetscReal (*N1)[2] = (PetscReal (*)[2]) p->shape[1];`
`+  PetscInt a,b,nen=p->nen;`
`+  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;`
`+    }`
`+    F[a] = Na * Forcing(x,y);`
`+  }`
`+  return 0;`
`+}`
`+`
`+PetscReal Exact(PetscReal x, PetscReal y)`
`+{`
`+  PetscReal pi = M_PI;`
`+  return sin(pi*x) * sin(pi*y);`
`+}`
`+`
`+#undef  __FUNCT__`
`+#define __FUNCT__ "Error"`
`+PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)`
`+{`
`+  PetscReal x = p->point[0];`
`+  PetscReal y = p->point[1];`
`+  PetscReal u_exact = Exact(x,y);`
`+  PetscScalar u;`
`+  IGAPointFormValue(p,U,&u);`
`+  PetscReal e = PetscAbsScalar(u - u_exact);`
`+  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);`
`+`
`+  PetscLogDouble tic,toc;`
`+  ierr = PetscGetTime(&tic);`
`+`
`+  IGA iga;`
`+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);`
`+  ierr = IGASetDim(iga,2);CHKERRQ(ierr);`
`+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);`
`+`
`+  IGABoundary bnd;`
`+  PetscInt dir,side;`
`+  for (dir=0; dir<2; dir++) {`
`+    for (side=0; side<2; side++) {`
`+      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);`
`+      ierr = IGABoundarySetValue(bnd,0,0.0);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);`
`+`
`+  PetscLogDouble ta1,ta2,ta;`
`+  ierr = PetscGetTime(&ta1);`
`+  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);`
`+  ierr = PetscGetTime(&ta2);`
`+  ta = ta2-ta1;`
`+`
`+  KSP ksp;`
`+  ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);`
`+  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);`
`+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);`
`+`
`+  PetscLogDouble ts1,ts2,ts;`
`+  ierr = PetscGetTime(&ts1);`
`+  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);`
`+  ierr = PetscGetTime(&ts2);`
`+  ts = ts2-ts1;`
`+`
`+  ierr = PetscGetTime(&toc);`
`+  PetscLogDouble tt = toc-tic;`
`+`
`+  Vec r;`
`+  PetscReal rnorm;`
`+  ierr = VecDuplicate(b,&r);CHKERRQ(ierr);`
`+  ierr = MatMult(A,x,r);CHKERRQ(ierr);`
`+  ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);`
`+  ierr = VecNorm(r,NORM_2,&rnorm);CHKERRQ(ierr);`
`+  ierr = VecDestroy(&r);CHKERRQ(ierr);`
`+`
`+  PetscInt its;`
`+  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);`
`+`
`+  ierr = IGAReset(iga);CHKERRQ(ierr);`
`+  for (dir=0; dir<2; dir++) {`
`+    IGARule rule;`
`+    ierr = IGAGetRule(iga,dir,&rule);CHKERRQ(ierr);`
`+    ierr = IGARuleInit(rule,9);CHKERRQ(ierr);`
`+  }`
`+  ierr = IGASetUp(iga);CHKERRQ(ierr);`
`+`
`+  PetscScalar error = 0;`
`+  ierr = IGAFormScalar(iga,x,1,&error,Error,PETSC_NULL);CHKERRQ(ierr);`
`+  error = PetscSqrtReal(PetscRealPart(error));`
`+`
`+  ierr = PetscPrintf(PETSC_COMM_WORLD,`
`+                     "Error=%E , ||b-Ax||=%E in %D its , "`
`+                     "Time: %f [assembly: %f (%.0f\%), solve: %f (%.0f\%)]\n",`
`+                     error,rnorm,its,tt,ta,ta/tt*100,ts,ts/tt*100);CHKERRQ(ierr);`
`+`
`+  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;`
`+}`

# demo/makefile

` Poisson: Poisson.o chkopts`
` 	\${CLINKER} -o \$@ \$< \${PETIGA_LIB}`
` 	\${RM} -f \$<`
`+LoggChallenge: LoggChallenge.o chkopts`
`+	\${CLINKER} -o \$@ \$< \${PETIGA_LIB}`
`+	\${RM} -f \$<`
` Bratu: Bratu.o BratuFJ.o chkopts`
` 	\${CLINKER} -o \$@ \$< BratuFJ.o \${PETIGA_LIB}`
` 	\${RM} -f \$< BratuFJ.o bratufj.mod`
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.