Permalink
Browse files

Added polynomial interpolation.

  • Loading branch information...
Jed Brown
Jed Brown committed Mar 28, 2008
1 parent f2e7037 commit 06cb1d46cea45bc6e493ef75e3399a5feb6b6a8e
Showing with 59 additions and 0 deletions.
  1. +2 −0 .gitignore
  2. +3 −0 Makefile
  3. +54 −0 util.C
View
@@ -8,3 +8,5 @@ poisson
elliptic
istest
core
+stokes
+TAGS
View
@@ -30,5 +30,8 @@ stokes : stokes.o chebyshev.o chkopts
istest : istest.o chkopts
${CLINKER} -o istest istest.o ${PETSC_LIB}
+util : util.o chkopts
+ ${CLINKER} -o util util.o ${PETSC_LIB}
+
etags : chebyshev.c stokes.C util.C
${ETAGS} $^
View
54 util.C
@@ -1,3 +1,7 @@
+#ifdef MAIN
+#include <stdio.h>
+#include <petscvec.h>
+#endif
class BlockIt {
public:
@@ -66,3 +70,53 @@ int productInt(int d, int dim[]) {
void zeroInt(int d, int v[]) {
for (int i=0; i < d; i++) v[i] = 0;
}
+
+#undef __FUNCT__
+#define __FUNCT__ "polyInterp"
+PetscErrorCode polyInterp(PetscInt n, PetscReal *x, PetscScalar *w, PetscReal x0, PetscReal x1, PetscScalar *f0, PetscScalar *f1)
+{
+ PetscScalar *tmp;
+ PetscInt o, e;
+ PetscReal y;
+
+ PetscFunctionBegin;
+ for (int di=1; di < n; di++) { // offset (column of table)
+ o = di % 2; e = (o+1) % 2;
+ for (int i=0; i < n-di; i++) { // nodes
+ w[i*4+2*o] = ((x0-x[i+di])*w[i*4+2*e ] + (x[i]-x0)*w[(i+1)*4+2*e ]) / (x[i] - x[i+di]);
+ w[i*4+2*o+1] = ((x1-x[i+di])*w[i*4+2*e+1] + (x[i]-x1)*w[(i+1)*4+2*e+1]) / (x[i] - x[i+di]);
+ }
+ }
+ *f0 = w[2*o];
+ *f1 = w[2*o+1];
+ PetscFunctionReturn(0);
+}
+
+#ifdef MAIN
+#define N 20
+PetscScalar func(PetscReal x) {
+ //return pow(x,6) + 3.1*pow(x,4) + 2.7*pow(x,3);
+ return cos(x);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[])
+{
+ PetscReal x[N], w[N*4], x0, x1, f0, f1;
+ PetscErrorCode ierr;
+
+ for (int order = 2; order < N; order++) {
+ for (int i=0; i < order; i++) {
+ x[i] = 1.0 + 1.0 * i;
+ w[i*4] = func(x[i]);
+ w[i*4+1] = w[i*4];
+ }
+ x0 = 1.43; x1 = 3.1;
+ ierr = polyInterp(order, x, w, x0, x1, &f0, &f1);CHKERRQ(ierr);
+ printf("[%d], %f ~ %f %f ~ %f\n", order, f0, func(x0), f1, func(x1));
+ }
+ return 0;
+}
+
+#endif

0 comments on commit 06cb1d4

Please sign in to comment.