diff --git a/src/findpts_imp.h b/src/findpts_imp.h
index b54acd9..6677c99 100644
--- a/src/findpts_imp.h
+++ b/src/findpts_imp.h
@@ -1,4 +1,3 @@
-
 #define obbox           TOKEN_PASTE(obbox_,D)
 #define local_hash_data TOKEN_PASTE(findpts_local_hash_data_,D)
 #define hash_data       TOKEN_PASTE(findpts_hash_data_,D)
@@ -23,6 +22,7 @@
 #define eval_src_pt         TOKEN_PASTE(eval_src_pt_ ,D)
 #define eval_out_pt         TOKEN_PASTE(eval_out_pt_ ,D)
 #define setup_aux           TOKEN_PASTE(setup_aux_,D)
+#define setup_fev_aux       TOKEN_PASTE(setup_fev_aux_,D)
 #define findpts_setup       TOKEN_PASTE(PREFIXED_NAME(findpts_setup_),D)
 #define findpts_free        TOKEN_PASTE(PREFIXED_NAME(findpts_free_ ),D)
 #define findpts             TOKEN_PASTE(PREFIXED_NAME(findpts_      ),D)
@@ -200,6 +200,8 @@ struct findpts_data {
   struct crystal cr;
   struct findpts_local_data local;
   struct hash_data hash;
+  struct array savpt;
+  uint         fevsetup;
 };
 
 static void setup_aux(
@@ -214,6 +216,7 @@ static void setup_aux(
                       npt_max, newt_tol);
   hash_build(&fd->hash,&fd->local.hd,fd->local.obb,nel,
              global_hash_size,&fd->cr);
+  fd->fevsetup = 0;
 }
 
 struct findpts_data *findpts_setup(
@@ -236,6 +239,7 @@ void findpts_free(struct findpts_data *fd)
   hash_free(&fd->hash);
   findpts_local_free(&fd->local);
   crystal_free(&fd->cr);
+  if (fd->fevsetup==1) array_free(&fd->savpt);
   free(fd);
 }
 
@@ -250,6 +254,7 @@ void findpts(      uint   *const  code_base   , const unsigned  code_stride   ,
              const double *const     x_base[D], const unsigned     x_stride[D],
              const uint npt, struct findpts_data *const fd)
 {
+  if (fd->fevsetup==1) array_free(&fd->savpt); fd->fevsetup=0;
   const uint np = fd->cr.comm.np, id=fd->cr.comm.id;
   struct array hash_pt, src_pt, out_pt;
   /* look locally first */
@@ -377,16 +382,14 @@ void findpts(      uint   *const  code_base   , const unsigned  code_stride   ,
 struct eval_src_pt { double r[D]; uint index, proc, el; };
 struct eval_out_pt { double out; uint index, proc; };
 
-void findpts_eval(
-        double *const  out_base, const unsigned  out_stride,
+static void setup_fev_aux(
   const uint   *const code_base, const unsigned code_stride,
   const uint   *const proc_base, const unsigned proc_stride,
   const uint   *const   el_base, const unsigned   el_stride,
   const double *const    r_base, const unsigned    r_stride,
-  const uint npt,
-  const double *const in, struct findpts_data *const fd)
+  const uint npt, struct findpts_data *const fd)
 {
-  struct array src, outpt;
+  struct array src;
   /* copy user data, weed out unfound points, send out */
   {
     uint index;
@@ -411,22 +414,59 @@ void findpts_eval(
     src.n = pt - (struct eval_src_pt*)src.ptr;
     sarray_transfer(struct eval_src_pt,&src,proc,1,&fd->cr);
   }
-  /* evaluate points, send back */
+  /* setup space for source points*/
   {
     uint n=src.n;
+    uint d;
     const struct eval_src_pt *spt;
-    struct eval_out_pt *opt;
+    struct eval_src_pt *opt;
     /* group points by element */
     sarray_sort(struct eval_src_pt,src.ptr,n, el,0, &fd->cr.data);
-    array_init(struct eval_out_pt,&outpt,n), outpt.n=n;
-    spt=src.ptr, opt=outpt.ptr;
-    for(;n;--n,++spt,++opt) opt->index=spt->index,opt->proc=spt->proc;
-    spt=src.ptr, opt=outpt.ptr;
-    findpts_local_eval(&opt->out ,sizeof(struct eval_out_pt),
-                       &spt->el  ,sizeof(struct eval_src_pt),
-                        spt->r   ,sizeof(struct eval_src_pt),
-                       src.n, in,&fd->local);
+    array_init(struct eval_src_pt,&fd->savpt,n), fd->savpt.n=n;
+    spt=src.ptr, opt=fd->savpt.ptr;
+    for(;n;--n,++spt,++opt) {
+        opt->index= spt->index;
+        opt->proc = spt->proc;
+        opt->el   = spt->el  ;
+        for(d=0;d<D;++d) opt->r[d] = spt->r[d];
+    }
     array_free(&src);
+  }
+}
+
+void findpts_eval(
+        double *const  out_base, const unsigned  out_stride,
+  const uint   *const code_base, const unsigned code_stride,
+  const uint   *const proc_base, const unsigned proc_stride,
+  const uint   *const   el_base, const unsigned   el_stride,
+  const double *const    r_base, const unsigned    r_stride,
+  const uint npt,
+  const double *const in, struct findpts_data *const fd)
+{ 
+  if (fd->fevsetup==0) {
+    setup_fev_aux(code_base,code_stride,
+                  proc_base,proc_stride,
+                    el_base,  el_stride,
+                     r_base,   r_stride,
+                        npt,         fd);
+    fd->fevsetup=1;
+  }
+  struct array outpt;
+  /* evaluate points, send back */
+  { 
+    uint n = fd->savpt.n;
+    struct eval_src_pt *opt;
+    struct eval_out_pt *opto;
+    array_init(struct eval_out_pt,&outpt,n), outpt.n=n;
+    opto=outpt.ptr;
+    opt=fd->savpt.ptr;
+    for(;n;--n,++opto,++opt) opto->index=opt->index,opto->proc=opt->proc;
+    opto=outpt.ptr;
+    opt=fd->savpt.ptr;
+    findpts_local_eval(&opto->out ,sizeof(struct eval_out_pt),
+                       &opt->el  ,sizeof(struct eval_src_pt),
+                       opt->r   ,sizeof(struct eval_src_pt),
+                       fd->savpt.n, in,&fd->local);
     sarray_transfer(struct eval_out_pt,&outpt,proc,1,&fd->cr);
   }
   /* copy results to user data */
@@ -445,6 +485,7 @@ void findpts_eval(
 #undef findpts_free
 #undef findpts_setup
 #undef setup_aux
+#undef setup_fev_aux
 #undef eval_out_pt
 #undef eval_src_pt
 #undef out_pt