diff --git a/Adaptive_benchmark/ArctanSphericalWaveFront.xinp b/Adaptive_benchmark/ArctanSphericalWaveFront.xinp
new file mode 100644
index 0000000..0371b00
--- /dev/null
+++ b/Adaptive_benchmark/ArctanSphericalWaveFront.xinp
@@ -0,0 +1,81 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<!--*************************************************************
+
+               Arctan Spherical Wave Front
+
+
+(http://math.nist.gov/cgi-bin/amr-display-problem.cgi)
+
+-nabla^2(u) = f,   in Omega
+         u  = g,   on Gamma
+
+Omega : (0,1)^3
+
+Exact solution: u = atan(alpha*(r-r0))
+                r = sqrt((x-x0)^2 + (y-y0)^2 + (z-z0)^2)
+Parameters:     (x0,y0,z0) center of the spherical wave front
+                r0         radius of the cirular wave front
+                alpha      controls the steepness of the wave front
+
+Parameter Reccomendations:
+   - Mild wave front
+        (x0,y0,z0) = (-.05 ,-.05, -.05)
+               r0  = 0.7
+            alpha  = 20
+
+   - Steep wave front
+        (x0,y0,z0) = (-.05 ,-.05, -.05)
+               r0  = 0.7
+            alpha  = 1000
+
+**********************************************************************-->
+
+
+<simulation>
+
+  <!-- General - geometry definitions !-->
+  <geometry sets="true">
+    <raiseorder patch="1" u="1" v="1" w="1"/>
+    <refine     patch="1" u="7" v="7" w="7"/>
+  </geometry>
+
+  <!-- General - boundary conditions !-->
+  <boundaryconditions>
+    <dirichlet set="Boundary" comp="1" type="anasol"/>
+  </boundaryconditions>
+
+  <postprocessing>
+    <vtfformat nu="2" nv="2" nw="2">BINARY</vtfformat>
+  </postprocessing>
+
+  <!-- Problem-specific block !-->
+  <poisson>
+    <source type="expression"> x0=-.05; y0=-.05; z0=-.05;  r=sqrt(pow(x-x0,2)+pow(y-y0,2)+pow(z-z0,2)); alpha=20; r0=0.7;
+
+(alpha*(2*x - 2*x0)^2)/(4*(alpha^2*(r0 - ((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2))^2 + 1)*((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(3/2)) - (3*alpha)/((alpha^2*(r0 - ((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2))^2 + 1)*((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2)) + (alpha*(2*y - 2*y0)^2)/(4*(alpha^2*(r0 - ((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2))^2 + 1)*((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(3/2)) + (alpha*(2*z - 2*z0)^2)/(4*(alpha^2*(r0 - ((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2))^2 + 1)*((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(3/2)) - (alpha^3*(2*z - 2*z0)^2*(r0 - ((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2)))/(2*(alpha^2*(r0 - ((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2))^2 + 1)^2*((x - x0)^2 + (y - y0)^2 + (z - z0)^2)) - (alpha^3*(2*x - 2*x0)^2*(r0 - ((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2)))/(2*(alpha^2*(r0 - ((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2))^2 + 1)^2*((x - x0)^2 + (y - y0)^2 + (z - z0)^2)) - (alpha^3*(2*y - 2*y0)^2*(r0 - ((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2)))/(2*(alpha^2*(r0 - ((x - x0)^2 + (y - y0)^2 + (z - z0)^2)^(1/2))^2 + 1)^2*((x - x0)^2 + (y - y0)^2 + (z - z0)^2))
+
+    </source>
+    <anasol type="expression">
+      <variables> x0=-.05; y0=-.05; z0=-.05; r=sqrt(pow(x-x0,2)+pow(y-y0,2)+pow(z-z0,2)); alpha=20; r0=0.7 </variables>
+
+      <primary>   atan(alpha*(r-r0)) </primary>
+      <secondary> -(alpha*(2*x - 2*x0))*pow(2*(alpha*alpha*pow(r - r0,2)+ 1)*r*r, -0.5) |
+                  -(alpha*(2*y - 2*y0))*pow(2*(alpha*alpha*pow(r - r0,2)+ 1)*r*r, -0.5) |
+                  -(alpha*(2*z - 2*z0))*pow(2*(alpha*alpha*pow(r - r0,2)+ 1)*r*r, -0.5)
+      </secondary>
+    </anasol>
+  </poisson>
+
+
+  <!--General - adaptive refinement parameters -->
+  <adaptive>
+    <maxstep>3</maxstep>
+    <beta type="AVERAGE">100</beta>
+    <maxdof>10000</maxdof>
+    <errtol>0.000001</errtol>
+    <knot_mult>1</knot_mult>
+    <scheme>isotropic_function</scheme>
+  </adaptive>
+
+</simulation>
diff --git a/Adaptive_benchmark/Fischera.g2 b/Adaptive_benchmark/Fischera.g2
new file mode 100644
index 0000000..b6d6013
--- /dev/null
+++ b/Adaptive_benchmark/Fischera.g2
@@ -0,0 +1,112 @@
+700 1 0 0
+3 0
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+-1.000000 -1.000000 -1.000000
+0.000000 -1.000000 -1.000000
+-1.000000 0.000000 -1.000000
+0.000000 0.000000 -1.000000
+-1.000000 -1.000000 0.000000
+0.000000 -1.000000 0.000000
+-1.000000 0.000000 0.000000
+0.000000 0.000000 0.000000
+700 1 0 0
+3 0
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+-1.000000 -1.000000 0.000000
+0.000000 -1.000000 0.000000
+-1.000000 0.000000 0.000000
+0.000000 0.000000 0.000000
+-1.000000 -1.000000 1.000000
+0.000000 -1.000000 1.000000
+-1.000000 0.000000 1.000000
+0.000000 0.000000 1.000000
+700 1 0 0
+3 0
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+-1.000000 0.000000 -1.000000
+0.000000 0.000000 -1.000000
+-1.000000 1.000000 -1.000000
+0.000000 1.000000 -1.000000
+-1.000000 0.000000 0.000000
+0.000000 0.000000 0.000000
+-1.000000 1.000000 0.000000
+0.000000 1.000000 0.000000
+700 1 0 0
+3 0
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+-1.000000 0.000000 0.000000
+0.000000 0.000000 0.000000
+-1.000000 1.000000 0.000000
+0.000000 1.000000 0.000000
+-1.000000 0.000000 1.000000
+0.000000 0.000000 1.000000
+-1.000000 1.000000 1.000000
+0.000000 1.000000 1.000000
+700 1 0 0
+3 0
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+0.000000 -1.000000 -1.000000
+1.000000 -1.000000 -1.000000
+0.000000 0.000000 -1.000000
+1.000000 0.000000 -1.000000
+0.000000 -1.000000 0.000000
+1.000000 -1.000000 0.000000
+0.000000 0.000000 0.000000
+1.000000 0.000000 0.000000
+700 1 0 0
+3 0
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+0.000000 -1.000000 0.000000
+1.000000 -1.000000 0.000000
+0.000000 0.000000 0.000000
+1.000000 0.000000 0.000000
+0.000000 -1.000000 1.000000
+1.000000 -1.000000 1.000000
+0.000000 0.000000 1.000000
+1.000000 0.000000 1.000000
+700 1 0 0
+3 0
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+2 2
+0.000000 0.000000 1.000000 1.000000
+0.000000 0.000000 -1.000000
+1.000000 0.000000 -1.000000
+0.000000 1.000000 -1.000000
+1.000000 1.000000 -1.000000
+0.000000 0.000000 0.000000
+1.000000 0.000000 0.000000
+0.000000 1.000000 0.000000
+1.000000 1.000000 0.000000
diff --git a/Adaptive_benchmark/FischeraCorner.xinp b/Adaptive_benchmark/FischeraCorner.xinp
new file mode 100644
index 0000000..5e7f90d
--- /dev/null
+++ b/Adaptive_benchmark/FischeraCorner.xinp
@@ -0,0 +1,108 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<!--*************************************************************
+
+               Fischera Corner with Vertex Singularity
+
+
+(http://math.nist.gov/cgi-bin/amr-display-problem.cgi)
+
+-nabla^2(u) = f,   in Omega
+         u  = g,   on Gamma
+
+Omega : (-1,1)^3 \ [0,1)^3
+
+Exact solution: u = r^q
+Parameters:     q > -0.5 determines the strength of the singularity
+
+**********************************************************************-->
+
+<simulation>
+
+  <!-- General - geometry definitions !-->
+  <geometry>
+    <patchfile>Fischera.g2</patchfile>
+<!--
+-->
+    <raiseorder lowerpatch="1" upperpatch="7" u="1" v="1" w="1"/>
+    <refine     lowerpatch="1" upperpatch="7" u="3" v="3" w="3"/>
+
+    <topology>
+      <connection master="5" mface="6" slave="6" sface="5"/>
+      <connection master="5" mface="4" slave="7" sface="3"/>
+      <connection master="3" mface="6" slave="4" sface="5"/>
+      <connection master="3" mface="2" slave="7" sface="1"/>
+      <connection master="2" mface="4" slave="4" sface="3"/>
+      <connection master="2" mface="2" slave="6" sface="1"/>
+      <connection master="1" mface="6" slave="2" sface="5"/>
+      <connection master="1" mface="4" slave="3" sface="3"/>
+      <connection master="1" mface="2" slave="5" sface="1"/>
+    </topology>
+
+    <topologysets>
+      <set name="Boundary" type="face">
+        <item patch="7">6</item>
+        <item patch="7">5</item>
+        <item patch="7">4</item>
+        <item patch="7">2</item>
+        <item patch="6">6</item>
+        <item patch="6">4</item>
+        <item patch="6">3</item>
+        <item patch="6">2</item>
+        <item patch="5">5</item>
+        <item patch="5">4</item>
+        <item patch="5">3</item>
+        <item patch="5">2</item>
+        <item patch="4">6</item>
+        <item patch="4">4</item>
+        <item patch="4">2</item>
+        <item patch="4">1</item>
+        <item patch="3">5</item>
+        <item patch="3">4</item>
+        <item patch="3">1</item>
+        <item patch="2">6</item>
+        <item patch="2">3</item>
+        <item patch="2">1</item>
+        <item patch="1">5</item>
+        <item patch="1">3</item>
+        <item patch="1">2</item>
+        <item patch="1">1</item>
+      </set>
+    </topologysets>
+
+  </geometry>
+
+  <!-- General - boundary conditions !-->
+  <boundaryconditions>
+    <dirichlet set="Boundary" comp="1" type="anasol"/>
+  </boundaryconditions>
+
+  <postprocessing>
+    <vtfformat nu="2" nv="2" nw="2">BINARY</vtfformat>
+  </postprocessing>
+
+  <!-- Problem-specific block !-->
+  <poisson>
+    <source type="expression">r2=(x*x+y*y+z*z); q=0.5; -q*(q+1)*pow(r2,q/2-1)</source>
+    <anasol type="expression">
+      <primary>  r =sqrt(x*x + y*y + z*z); q=0.5;  pow(r,q)  </primary>
+      <secondary>r2=    (x*x + y*y + z*z); q=0.5; -q*x*pow(r2,q/2 - 1) |
+                 r2=    (x*x + y*y + z*z); q=0.5; -q*y*pow(r2,q/2 - 1) |
+                 r2=    (x*x + y*y + z*z); q=0.5; -q*z*pow(r2,q/2 - 1)
+      </secondary>
+    </anasol>
+
+  </poisson>
+
+
+  <!--General - adaptive refinement parameters -->
+  <adaptive>
+    <maxstep>4</maxstep>
+    <beta type="AVERAGE">100</beta>
+    <maxdof>100000</maxdof>
+    <errtol>0.000001</errtol>
+    <knot_mult>1</knot_mult>
+    <scheme>isotropic_function</scheme>
+  </adaptive>
+
+</simulation>
diff --git a/Adaptive_benchmark/FischeraLine.xinp b/Adaptive_benchmark/FischeraLine.xinp
new file mode 100644
index 0000000..7426244
--- /dev/null
+++ b/Adaptive_benchmark/FischeraLine.xinp
@@ -0,0 +1,105 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<!--*************************************************************
+
+          Fichera Corner with Vertex and Edge Singularities
+
+
+(http://math.nist.gov/cgi-bin/amr-display-problem.cgi)
+
+-nabla^2(u) = 1/r,   in Omega
+         u  = 0,     on Gamma
+
+Omega : (-1,1)^3 \ [0,1)^3
+
+Exact solution: (refence solution available online)
+
+**********************************************************************-->
+
+<simulation>
+
+  <!-- General - geometry definitions !-->
+  <geometry>
+    <patchfile>Fischera.g2</patchfile>
+<!--
+    <raiseorder patch="1" u="1" v="1" w="1"/>
+-->
+
+    <topology>
+      <connection master="5" mface="6" slave="6" sface="5"/>
+      <connection master="5" mface="4" slave="7" sface="3"/>
+      <connection master="3" mface="6" slave="4" sface="5"/>
+      <connection master="3" mface="2" slave="7" sface="1"/>
+      <connection master="2" mface="4" slave="4" sface="3"/>
+      <connection master="2" mface="2" slave="6" sface="1"/>
+      <connection master="1" mface="6" slave="2" sface="5"/>
+      <connection master="1" mface="4" slave="3" sface="3"/>
+      <connection master="1" mface="2" slave="5" sface="1"/>
+    </topology>
+
+    <topologysets>
+      <set name="Boundary" type="face">
+        <item patch="7">6</item>
+        <item patch="7">5</item>
+        <item patch="7">4</item>
+        <item patch="7">2</item>
+        <item patch="6">6</item>
+        <item patch="6">4</item>
+        <item patch="6">3</item>
+        <item patch="6">2</item>
+        <item patch="5">5</item>
+        <item patch="5">4</item>
+        <item patch="5">3</item>
+        <item patch="5">2</item>
+        <item patch="4">6</item>
+        <item patch="4">4</item>
+        <item patch="4">2</item>
+        <item patch="4">1</item>
+        <item patch="3">5</item>
+        <item patch="3">4</item>
+        <item patch="3">1</item>
+        <item patch="2">6</item>
+        <item patch="2">3</item>
+        <item patch="2">1</item>
+        <item patch="1">5</item>
+        <item patch="1">3</item>
+        <item patch="1">2</item>
+        <item patch="1">1</item>
+      </set>
+    </topologysets>
+
+  </geometry>
+
+  <!-- General - boundary conditions !-->
+  <boundaryconditions>
+    <dirichlet set="Boundary" comp="1"/>
+  </boundaryconditions>
+
+  <postprocessing>
+    <vtfformat nu="2" nv="2" nw="2">BINARY</vtfformat>
+    <projection type="cgl2"/>
+  </postprocessing>
+
+  <!-- Problem-specific block !-->
+  <poisson>
+    <source type="expression">r=sqrt(x*x+y*y+z*z); 1/r</source>
+  </poisson>
+
+
+  <!--General - adaptive refinement parameters -->
+  <adaptive>
+    <maxstep>4</maxstep>
+    <beta type="AVERAGE">100</beta>
+    <maxdof>10000</maxdof>
+    <errtol>0.000001</errtol>
+    <knot_mult>1</knot_mult>
+    <scheme>isotropic_function</scheme>
+    <!--
+      <use_norm>1</use_norm>
+      <use_sub_norm>1</use_sub_norm>
+      <store_eps_mesh/>
+      <test_linear_independence/>
+    -->
+  </adaptive>
+
+</simulation>
diff --git a/Adaptive_benchmark/GenFischera.py b/Adaptive_benchmark/GenFischera.py
new file mode 100644
index 0000000..f61d951
--- /dev/null
+++ b/Adaptive_benchmark/GenFischera.py
@@ -0,0 +1,16 @@
+from splipy import *
+from splipy.IO import *
+import splipy.curve_factory as cf
+import splipy.surface_factory as sf
+import splipy.volume_factory as vf
+import numpy as np
+
+v = Volume() - [1,1,1];
+
+with G2("Fischera.g2") as f:
+    for i in range(2):
+        for j in range(2):
+            for k in range(2):
+                if i==j==k==1: # skip first octant
+                    continue
+                f.write(v + [i,j,k])
diff --git a/Adaptive_benchmark/Peak3D.xinp b/Adaptive_benchmark/Peak3D.xinp
new file mode 100644
index 0000000..5a67649
--- /dev/null
+++ b/Adaptive_benchmark/Peak3D.xinp
@@ -0,0 +1,76 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<!--*************************************************************
+
+                           Peak 3D
+
+
+(http://math.nist.gov/cgi-bin/amr-display-problem.cgi)
+
+-nabla^2(u) = f,   in Omega
+         u  = g,   on Gamma
+
+Omega : (0,1)^3
+
+Exact solution: u = exp(-alpha*r^2)
+                r = sqrt((x-x0)^2 + (y-y0)^2 + (z-z0)^2)
+Parameters:     (x0,y0,z0) center of the peak
+                alpha      the strength of the peak
+
+Parameter Reccomendations:
+   - Mild peak
+        (x0,y0,z0) = (0.5 ,0.5, 0.5)
+            alpha  = 1e3
+
+   - sharp peak
+        (x0,y0,z0) = (0.5 ,0.5, 0.5)
+            alpha  = 1e5
+
+**********************************************************************-->
+
+
+<simulation>
+
+  <!-- General - geometry definitions !-->
+  <geometry sets="true">
+    <raiseorder patch="1" u="1" v="1" w="1"/>
+    <refine     patch="1" u="7" v="7" w="7"/>
+  </geometry>
+
+  <!-- General - boundary conditions !-->
+  <boundaryconditions>
+    <dirichlet set="Boundary" comp="1" type="anasol"/>
+  </boundaryconditions>
+
+  <postprocessing>
+    <vtfformat nu="2" nv="2" nw="2">BINARY</vtfformat>
+  </postprocessing>
+
+  <!-- Problem-specific block !-->
+  <poisson>
+    <source type="expression">  x0=0.5; y0=0.5; z0=0.5; r2=pow(x-x0,2)+pow(y-y0,2)+pow(z-z0,2); alpha=1e3;
+6*alpha*exp(-alpha*r2) - alpha^2*exp(-alpha*r2)*(2*x - 2*x0)^2 - alpha^2*exp(-alpha*r2)*(2*y - 2*y0)^2 - alpha^2*exp(-alpha*r2)*(2*z - 2*z0)^2
+    </source>
+    <anasol type="expression">
+      <variables> x0=0.5; y0=0.5; z0=0.5; r2=pow(x-x0,2)+pow(y-y0,2)+pow(z-z0,2); alpha=1e3; </variables>
+
+      <primary>    exp(-alpha*r2) </primary>
+      <secondary>  alpha*exp(-alpha*r2)*(2*x - 2*x0) |
+                   alpha*exp(-alpha*r2)*(2*y - 2*y0) |
+                   alpha*exp(-alpha*r2)*(2*z - 2*z0)
+   </secondary>
+    </anasol>
+  </poisson>
+
+
+  <!--General - adaptive refinement parameters -->
+  <adaptive>
+    <maxstep>3</maxstep>
+    <beta type="AVERAGE">100</beta>
+    <maxdof>10000</maxdof>
+    <errtol>0.000001</errtol>
+    <knot_mult>1</knot_mult>
+    <scheme>isotropic_function</scheme>
+  </adaptive>
+
+</simulation>
diff --git a/SIMPoisson.h b/SIMPoisson.h
index 37c554c..7b25a09 100644
--- a/SIMPoisson.h
+++ b/SIMPoisson.h
@@ -130,6 +130,8 @@ template<class Dim> class SIMPoisson : public SIMMultiPatchModelGen<Dim>
       // Project the FE stresses onto the splines basis
       this->setMode(SIM::RECOVERY);
       size_t i = 0;
+      if (projections == &myProj)
+        myProj.resize(this->opt.project.size());
       for (auto pit  = this->opt.project.begin();
                 pit != this->opt.project.end(); ++i, ++pit) {
         Matrix ssol;