diff --git a/src/geometry/block.cc b/src/geometry/block.cc
index 76795bd..0f6f72c 100644
--- a/src/geometry/block.cc
+++ b/src/geometry/block.cc
@@ -1,365 +1,367 @@
 #include "block.hh"
 #include "block_surface_manipulations.hh"
 #include <limits>
 #include <cmath>
 #include <vector>
 #include <deque>
 #include <utility>
 #include <tuple>
 #include <algorithm>
 
 
 #include <cadd_mesh.hh>
 #include <exception>
 #include <sstream>
 
 
 class BlockException: public std::exception {
 public:
   BlockException(std::stringstream & _what):what_stream(_what.str()) {};
   BlockException(std::string &_what):what_stream(_what) {};
   ~BlockException() throw() {};
   virtual const char* what() const throw() {
     return this->what_stream.c_str();
   }
 private:
   std::string what_stream;
 };
 
 template<Uint DIM>
 Block<DIM>::Block(const std::vector<Real> & _min_max, std::string _name)
   :Geometry<DIM>(_name) {
   if (_min_max.size() != 2*DIM) {
     std::stringstream error;
     error << "The length of 'min_max' should be = " << 2*DIM;
 
     throw BlockException(error);
   }
   for (Uint i = 0 ;  i < 2*DIM ; ++i) {
     this->min_max[i] = _min_max[i];
   }
 }
 
 template<Uint DIM>
 Block<DIM>::Block(const Real * min_max, std::string _name)
   :Geometry<DIM>(_name) {
   for (Uint i = 0; i < 2 * DIM; ++i) {
     this->min_max[i] = min_max[i];
   }
 }
 
 template<Uint DIM>
 Block<DIM>::Block(const Block & other)
   :Geometry<DIM>(other.name) {
   for (Uint i = 0 ;  i < 2 * DIM ; ++i) {
     this->min_max[i] = other.min_max[i];
   }
 }
 
 
 template <Uint DIM>
 void Block<DIM>::shift(const PointRef<DIM> & offset) {
   for (Uint i = 0 ;  i < DIM ; ++i) {
     this->min_max[2 * i + 0] += offset.getComponent(i);
     this->min_max[2 * i + 1] += offset.getComponent(i);
   }
 }
 
 template<Uint DIM>
 InsideObject Block<DIM>::from_border(const Real * point, Real tol) const {
   Real lower_violation = -std::numeric_limits<Real>::max();
   Real upper_violation = -std::numeric_limits<Real>::max();
   for (Uint i = 0; i < DIM; ++i) {
     keep_max(lower_violation, this->min_max[2*i] - point[i]);
     keep_max(upper_violation, point[i] - this->min_max[2*i+1]);
   }
 
   return {max(lower_violation, upper_violation),
       max(lower_violation-tol, upper_violation+tol) <= 0};
 }
 
 template<Uint DIM, bool eval_min, bool eval_max>
 inline void in_direction_helper(Real & min_val, Real & max_val,
                                 const Real * direction, Real * point) {
   Real dot_prod = 0;
   Real norm = 0;
   for (Uint i = 0; i < DIM; ++i) {
     dot_prod += point[i] * direction[i];
     norm += direction[i] * direction[i];
   }
   dot_prod /= sqrt(norm);
 
   if (eval_min) {
     keep_min(min_val, dot_prod);
   }
   if (eval_max) {
     keep_max(max_val, dot_prod);
   }
 }
 
 template<Uint DIM>
 Real Block<DIM>::min_in_direction(const Real * dir, const Real & unit) const {
   Real point[DIM];
   Real min_val = std::numeric_limits<Real>::max();
   for (Uint i = 0; i < 2; ++i) {
     point[0] =  min_max[i];
     for (Uint j = 0; j < 2; ++j) {
       point[1] = min_max[2 + j];
       if (DIM == 3) {
         for (Uint k = 0; k < 2; ++k) {
           point[2] = min_max[4 + k];
           in_direction_helper<DIM, true, false>(min_val, min_val, dir, point);
         }
       } else if (DIM == 2) {
         in_direction_helper<DIM, true, false>(min_val, min_val, dir, point);
       }
     }
   }
   return min_val/unit;
 }
 
 
 template<Uint DIM>
 Real Block<DIM>::max_in_direction(const Real * dir, const Real & unit) const {
   Real point[DIM];
   Real max_val = -std::numeric_limits<Real>::max();
   for (Uint i = 0; i < 2; ++i) {
     point[0] =  min_max[i];
     for (Uint j = 0; j < 2; ++j) {
       point[1] = min_max[2 + j];
       if (DIM == 3) {
         for (Uint k = 0; k < 2; ++k) {
           point[2] = min_max[4 + k];
           in_direction_helper<DIM, false, true>(max_val, max_val, dir, point);
         }
       } else if (DIM == 2) {
         in_direction_helper<DIM, false, true>(max_val, max_val, dir, point);
       }
     }
   }
   return max_val/unit;
 }
 
 template<Uint DIM>
 Real Block<DIM>::size_in_direction(const Real * dir, const Real & unit) const {
   Real point[DIM];
   Real min_val = std::numeric_limits<Real>::max();
   Real max_val = -std::numeric_limits<Real>::max();
   for (Uint i = 0; i < 2; ++i) {
     point[0] =  min_max[i];
     for (Uint j = 0; j < 2; ++j) {
       point[1] = min_max[2 + j];
       if (DIM == 3) {
         for (Uint k = 0; k < 2; ++k) {
           point[2] = min_max[4 + k];
           in_direction_helper<DIM, true, true>(min_val, max_val, dir, point);
         }
       } else if (DIM == 2) {
         in_direction_helper<DIM, true, true>(min_val, max_val, dir, point);
       }
     }
   }
   return (max_val - min_val)/unit;
 }
 
 template<Uint DIM>
 void Block<DIM>::printself(std::ostream & stream, int indent) const {
   indent_space(stream, indent);
   stream << "Block " << this->get_name() << ": " << std::endl;
   indent_space(stream, indent);
   stream << "x_min " << this->min_max[0] << ", x_max = " << this->min_max[1] << std::endl;
   indent_space(stream, indent);
   stream << "y_min " << this->min_max[2] << ", y_max = " << this->min_max[3] << std::endl;
   if (DIM == 3) {
     indent_space(stream, indent);
     stream << "z_min " << this->min_max[4] << ", z_max = " << this->min_max[5] << std::endl;
   }
 }
 
 template<Uint DIM>
 inline void enforce_slave(const PointRef<DIM> & point,
                           const Mesh<DIM> & auxiliary,
                           const std::vector<Real> & refinement,
                           const Real & step_size,
                           PointContainer<DIM> & points) {
   Real ref_sq = step_size * auxiliary.interpolate(point, refinement);
   ref_sq *= ref_sq;
   std::deque<Uint> delete_points;
   for (Uint point_id = 0 ;  point_id < points.getSize() ; ++point_id) {
     if (ref_sq > distance2(point, points.getPoint(point_id))) {
       delete_points.push_front(point_id);
     }
   }
   auto point_it = delete_points.begin();
   const auto point_end = delete_points.end();
   for(; point_it != point_end; ++point_it) {
     points.deletePoint(*point_it);
   }
   points.addPoint(point);
 }
 
 template <Uint DIM>
 inline bool masterdom(const PointRef<DIM> & point, const Real * min_max,
                       const bool periodicity[DIM], bool master[DIM], bool & slave) {
   bool is_master = false;
   bool is_not_master = false;
   slave = false;
   for (Uint i = 0 ;  i < DIM ; ++i) {
     master[i] = (point[i] == min_max[2*i]) && periodicity[i];
     // not a master if i'm at the max of any periodic direction
     is_not_master |= (point[i] == min_max[2*i + 1]) && periodicity[i];
     is_master |= master[i];
     slave |= (point[i] == min_max[2*i + 1]) && periodicity[i];
   }
   return is_master && !is_not_master;
 }
 
 template <Uint DIM>
 inline void drive_slaves(const PointRef<DIM> & point, const Real * min_max,
                          const bool periodicity[DIM],
                          const Mesh<DIM> & auxiliary,
                          const std::vector<Real> & refinement,
                          const Real & step_size,
                          PointContainer<DIM> & points) {
   bool master[DIM];
   bool slave;
   if (masterdom(point, min_max, periodicity, master, slave)) {
     PointContainer<DIM> slaves;
     slaves.addPoint(point);
     PointRef<DIM> tmp;
     for (Uint i = 0 ;  i < DIM ; ++i) {
       if ( master[i]) {
         Uint current_nb_slaves = slaves.getSize();
         for (Uint slave_num = 0 ;  slave_num < current_nb_slaves ; ++slave_num) {
           tmp = slaves.getConstPoint(slave_num);
           tmp[i] += min_max[2*i+1]-min_max[2*i];
           slaves.addPoint(tmp);
         }
       }
     }
     for (Uint slave_num = 0 ;  slave_num < slaves.getSize() ; ++slave_num) {
       enforce_slave(slaves.getConstPoint(slave_num),
                     auxiliary,
                     refinement,
                     step_size,
                     points);
     }
   } else if (!slave) {
     points.addPoint(point);
   }
 }
 template <Uint DIM>
 inline void micro_add_surface_point(PointRef<DIM> & point_hi,
                                     PointRef<DIM> & point_lo,
                                     const Uint xi,
                                     const Real * min_max,
                                     const Mesh<DIM> & auxiliary,
                                     const std::vector<Real> & refinement,
                                     const bool periodicity[DIM],
                                     const Real & repr_const,
                                     PointContainer<DIM> & points) {
   point_hi[xi] = point_lo[xi] = min_max[2*xi];
   if (check_if_ok(point_lo, auxiliary, refinement, repr_const, points)) {
     drive_slaves(point_lo, min_max, periodicity, auxiliary, refinement,
                  repr_const,  points);
   }
   if (check_if_ok(point_hi, auxiliary, refinement, repr_const, points)) {
     drive_slaves(point_hi, min_max, periodicity, auxiliary, refinement,
                  repr_const, points);
   }
   point_hi[xi] = point_lo[xi] = min_max[2*xi+1];
   if (!periodicity[xi]) {
     add_if_ok(point_hi, auxiliary, refinement, repr_const, points);
     add_if_ok(point_lo, auxiliary, refinement, repr_const, points);
   }
 }
 
 template <Uint DIM>
 void recurse_through_dimensions(PointRef<DIM> & point,
                                 const Uint xi,
                                 const Uint my_offset_dim,
                                 const Real * min_max,
                                 const Mesh<DIM> & auxiliary,
                                 const std::vector<Real> & refinement,
                                 const bool periodicity[DIM],
                                 const Real step_size[DIM],
                                 const Real & repr_const,
                                 PointContainer<DIM> & points,
                                 Uint * starts) {
   Uint my_dim = (xi+my_offset_dim) % DIM;
   Real length = min_max[2*my_dim + 1] - min_max[2*my_dim];
   Uint nb_steps = length/step_size[my_dim];
   PointRef<DIM> point_lo = point;
   PointRef<DIM> point_hi = point;
 
   // We mesh from the outside inwards, to make sure that corners/edges get
   // points
   for (Uint step = starts[my_offset_dim-1]; step < nb_steps/2+1; ++step) {
     point_lo[my_dim] = min_max[2*my_dim]     + step * step_size[my_dim];
     point_hi[my_dim] = min_max[2*my_dim + 1] - step * step_size[my_dim];
     if (DIM == my_offset_dim + 1) {
       micro_add_surface_point(point_hi, point_lo, xi, min_max, auxiliary,
                               refinement, periodicity, repr_const, points);
     } else {
       recurse_through_dimensions(point_lo, xi, my_offset_dim + 1, min_max,
                                  auxiliary, refinement, periodicity, step_size,
                                  repr_const, points, starts);
       recurse_through_dimensions(point_hi, xi, my_offset_dim + 1, min_max,
                                  auxiliary, refinement, periodicity, step_size,
                                  repr_const, points, starts);
     }
   }
 }
 
 
 template <Uint DIM>
 void Block<DIM>::generate_surface_points(const Mesh<DIM> & auxiliary,
                                          const std::vector<Real> & refinement,
                                          const Real step_size[DIM],
                                          const Real & repr_const,
                                          const bool periodicity[DIM],
                                          PointContainer<DIM> & points) const {
   PointRef<DIM> point;
   Uint starts[DIM] = {0};
 
   for (Uint xi = 0 ;  xi < DIM ; ++xi) {
     starts[DIM-1-xi] = 1;
     recurse_through_dimensions(point, xi, 1, this->min_max, auxiliary,
                                refinement, periodicity, step_size,
                                repr_const, points, starts);
   }
   points.cleanDuplicates(1.);
 }
 
 
 
 template<Uint DIM>
 void Block<DIM>::complement_periodically(const Mesh<DIM> & auxiliary,
                                          const std::vector<Real> & refinement,
                                          const PointContainer<DIM> & points,
                                          PointContainer<DIM> & complement,
                                          int direction,
                                          Uint dim) const {
   // relevant coordinate of the good and the bad side
   Real ok_coord, bad_coord;
   if (direction > 0) {
     ok_coord  = this->min_max[2*dim + 0];
     bad_coord = this->min_max[2*dim + 1];
   } else {
     ok_coord  = this->min_max[2*dim + 1];
     bad_coord = this->min_max[2*dim + 0];
   }
   Real tol = std::max(fabs(ok_coord), 1.);
   tol = 2000*tol*std::numeric_limits<Real>::epsilon();
   PointRef<DIM> point;
-  for (Uint i = 0; i < points.getSize() ; ++i) {
-    point = points.getConstPoint(i);
+
+  PointContainer<DIM> container = points + complement;
+  for (Uint i = 0; i < container.getSize() ; ++i) {
+    point = container.getConstPoint(i);
     if (fabs(point.getComponent(dim) - ok_coord) < tol) {
       point.getArray()[dim] = bad_coord;
       add_periodic_complement(auxiliary, refinement, point, points, complement);
     }
   }
 }
 
 template class Block<twoD>;
 template class Block<threeD>;
diff --git a/src/geometry/geometry.hh b/src/geometry/geometry.hh
index eb30156..f16f85f 100644
--- a/src/geometry/geometry.hh
+++ b/src/geometry/geometry.hh
@@ -1,212 +1,212 @@
 /**
  * @file   geometry.hh
  * @author Till Junge <junge@lsmspc42.epfl.ch>
  * @date   Wed Apr  4 15:15:10 2012
  *
  * @brief  Geometry motherclass (interface) for CADD mesher
  *
  * @section LICENSE
  *
  * <insert lisence here>
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #ifndef __GEOMETRY_HH__
 #define __GEOMETRY_HH__
 
 #include "common.hh"
 #include "point_container.hh"
 #include <sstream>
 
 //#include "cadd_mesh.hh"
 
 template <Uint DIM>
 class Mesh;
 
 struct InsideObject {
   Real distance;
   bool is_inside;
 };
 
 /*! \brief Geometry interface for CADD mesher
  * You will never use this class.*/
 template <Uint DIM>
 class Geometry {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   /*! Constructor
    *\param _name a string identifier for the geometry */
   Geometry(std::string _name = ""):name(_name){};
   virtual ~Geometry(){};
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
 
   /// function to print the content of the class
   virtual void printself(std::ostream & stream, int indent = 0) const
   {
     stream << "This is just a interface class and this function should never "
            << "ever be called" <<std::endl;
   }
 
   virtual void shift(const PointRef<DIM> & offset) = 0;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 
   /// returns the string identifier
   const std::string & get_name() const
   {
     return this->name;
   }
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 public:
   /*! creates a copy of the object (calling new) and returns a Geometry
    *  pointer*/
   virtual Geometry * resolveType() const  = 0;
   /*! test whether a point is inside the domain
    * \param point coordinates of the point to test*/
   virtual bool is_inside(const Real * point, Real tol) const = 0;
   bool py_is_inside(const PointRef<DIM> & point, Real tol) const {
     return this->is_inside(point, tol);
   }
   virtual inline bool is_inside(const PointRef<DIM> & point, Real tol) const {
     return this->is_inside(&(point.getComponent(0)), tol);}
   /*! Returns the minimum value of all constraint violations (negative return
    * value means point is inside and return value from the nearest boundary)
    * \param point coordintates of the point to test*/
   virtual InsideObject from_border(const Real * point, Real tol) const = 0;
   inline InsideObject from_border(const PointRef<DIM> & point, Real tol) const {
     return this->from_border(&(point.getComponent(0)), tol);}
 
   /*! projects the geometry on a line of direction \a dir and returns the
    * minimum of the line segment spanned
    * \param dir direction vector (will be scaled to unity by the method
    * \param unit unit in which the distance should be expressed, by default unity*/
   virtual Real  min_in_direction(const Real * dir, const Real & unit = 1) const = 0;
 
   /*! projects the geometry on a line of direction \a dir and returns the
    * maximum of the line segment spanned
    * \param dir direction vector (will be scaled to unity by the method
    * \param unit unit in which the distance should be expressed, by default unity*/
   virtual Real  max_in_direction(const Real * dir, const Real & unit = 1) const = 0;
   /*! projects the geometry on a line of direction \a dir and returns the length
    * of the projection
    * \param dir direction vector (will be scaled to unity by the method
    * \param unit unit in which the distance should be expressed, by default unity*/
   virtual Real size_in_direction(const Real * dir, const Real & unit = 1) const ;
 
   /*! put nodes onto the surface of the geometry
    *  \param auxiliary mesh for refinement interpolation
    *  \param refinement vector of refinements used with the auxiliary mesh
    *  \param step_size array of periodic cell sizes as returned by
    *         lattice::getSpatialConstants()
    *  \param repr_const representative constant of the periodic cell
    *  \param periodicity self-explanatory
    *  \param points container of nodal points for meshing
    */
   virtual void generate_surface_points(const Mesh<DIM> & auxiliary,
                                        const std::vector<Real> & refinement,
                                        const Real step_size[DIM],
                                        const Real & repr_const,
                                        const bool periodicity[DIM],
                                        PointContainer<DIM> & points) const;
 
   /*! When the dimensions of the overall geometry are not a power of two times
    *  the lattice size, the domain can be not filled correctly. This method 
    *  corrects the fucked up side by periodically repeating the ok side.
    *  \param direction positive values mean that the min side is used to augment
    *                   the max size, negative values the inverse.
    *  \param dim dimension to fix: 0:x, 1:y, 2:z
    */
   virtual void complement_periodically(const Mesh<DIM> & auxiliary,
                                        const std::vector<Real> & refinement,
                                        const PointContainer<DIM> & points,
                                        PointContainer<DIM> & complement,
                                        int direction,
                                        Uint dim) const;
 
 protected:
   std::string name;
 };
 
 /// standard output stream operator
 template<Uint DIM>
 inline std::ostream & operator <<(std::ostream & stream, const Geometry<DIM> & _this)
 {
   _this.printself(stream);
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 template <Uint DIM, typename subGeom>
 Geometry<DIM> * newGeom(const subGeom & geometry) {
   return static_cast<Geometry<DIM>*>(new subGeom(geometry));
 }
 
 
 /* -------------------------------------------------------------------------- */
 template<Uint DIM>
 inline bool check_if_ok(const PointRef<DIM> & point,
                         const Mesh<DIM> & auxiliary,
                         const std::vector<Real> & refinement,
                         const Real & step_size,
                         PointContainer<DIM> & points) {
   Real ref_sq = square(step_size * auxiliary.interpolate(point,
                                                          refinement));
   Real step_sq = step_size*step_size;
   // if refinement is smaller than 1 => full resolution
   if (ref_sq < step_sq) {
     return false;
   }
   // check whether the new point in too close to an existing point
   for (Uint point_id = 0; point_id < points.getSize() ; ++point_id) {
     if (ref_sq > distance2<DIM>(point, points.getPoint(point_id))) {
       return false;
     }
   }
   return true;
 }
 
 /* -------------------------------------------------------------------------- */
 template<Uint DIM>
 inline bool add_if_ok(const PointRef<DIM> & point,
                       const Mesh<DIM> & auxiliary,
                       const std::vector<Real> & refinement,
                       const Real & step_size,
                       PointContainer<DIM> & points) {
   if (check_if_ok(point, auxiliary, refinement, step_size, points)) {
     points.addPoint(point);
     return true;
   } else {
     return false;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template<Uint DIM>
 bool add_periodic_complement(const Mesh<DIM> & auxiliary,
                              const std::vector<Real> & refinement,
-                             const PointRef<DIM> point,
+                             const PointRef<DIM> & point,
                              const PointContainer<DIM> & points,
                              PointContainer<DIM> & complement) {
   Real local_refinement = auxiliary.interpolate(point, refinement);
   if (local_refinement <= 1.) {
     complement.addPoint(point);
     return true;
   } else {
     return false;
   }
 }
 
 
 #endif /* __GEOMETRY_HH__ */