diff --git a/src/domain/domain.cc b/src/domain/domain.cc
index 1ce43d0..f69bd46 100644
--- a/src/domain/domain.cc
+++ b/src/domain/domain.cc
@@ -1,276 +1,276 @@
 
 #include "domain.hh"
 #include "tree.hh"
 #include <limits>
 #include <cmath>
 #include <iterator>
 #include <stdexcept>
 
 
 template <Uint DIM>
 Domain<DIM>::Domain(const Geometry<DIM> & overall_,
                     const Geometry<DIM> & refined_,
                     const Geometry<DIM> & atomistic_,
                     const Geometry<DIM> & atomistic_skinned_,
                     const Lattice<DIM> & lattice_,
                     bool boundaries_special_treatment_,
                     const PointRef<DIM> & centre,
                     Real tolerance,
                     const std::string & qh_flags_)
   :auxiliary_mesh(NULL), tree(NULL), main_mesh(NULL),
   filled(false), main_meshed(false),
   boundaries_special_treatment(boundaries_special_treatment_),
   mesh_quality(-std::numeric_limits<Real>::max()),
   tol(tolerance),
   qh_flags(qh_flags_)
 {
   this->overall = overall_.resolveType();
   this->refined = refined_.resolveType();
   this->atomistic = atomistic_.resolveType();
   this->atomistic_skinned = atomistic_skinned_.resolveType();
   this->lattice = lattice_.resolveType();
 
   Real max_size = -std::numeric_limits<Real>::max();
   for (Uint dir = 0; dir < DIM ; ++dir) {
     Real size = this->overall->size_in_direction
       (this->lattice->getLatticeVector(dir),
        this->lattice->getConstant(dir));
     keep_max(max_size, size);
   }
   // compute the next larger power of 2
   int exponent = log(max_size)/log(2) +1;
   Uint pow2size = pow(2, exponent);
   if (&centre == NULL) {
     this->tree = new Tree<DIM>(NULL, 2*pow2size, this->lattice, this);
   } else {
     this->tree = new Tree<DIM>(centre.getArray(), 2*pow2size, this->lattice, this);
   }
 
 }
 
 
 template <Uint DIM>
 Domain<DIM>::~Domain() {
   delete this->auxiliary_mesh;
   delete this->main_mesh;
   delete this->tree;
 }
 
 
 
 template <Uint DIM>
 void Domain<DIM>::setPointRefinement(const Real * coords, Real & refinement ) {
   this->refinement_points.addPoint(coords);
   this->refinement.push_back(refinement);
 }
 
 template <Uint DIM>
 void Domain<DIM>::fill_points(bool periodicity[DIM]) throw(std::string){
   this->initialiseAuxiliaryMesh();
-  this->auxiliary_mesh->dump_paraview("test");
+  //this->auxiliary_mesh->dump_paraview("test");
   if (this->boundaries_special_treatment) {
     this->overall->generate_surface_points(*this->auxiliary_mesh,
                                            this->refinement,
                                            this->lattice->getSpatialConstants(),
                                            this->lattice->getRepresentativeConstant(),
                                            periodicity,
                                            this->points);
   }
   this->tree->fill(this->points);
   this->filled = true;
 
   this->points.reportDuplicates(1);
 }
 
 template<Uint DIM>
 void Domain<DIM>::complement_periodically(int direction, Uint dim) throw(std::string){
   if (!this->filled) {
     throw("Mesh has got to be filled before it can be complemented");
   }
   this->overall->complement_periodically(*this->auxiliary_mesh,
                                          this->refinement,
                                          this->points,
                                          direction,
                                          dim);
 }
 
 template <Uint DIM>
 void Domain<DIM>::printself(std::ostream & stream, int indent) const{
   indent_space(stream, indent);
   stream << "PRINTSELF IS NOT YET IMOPLEMENTED FOR CLASS DOMAIN";
 }
 
 template <Uint DIM>
 void Domain<DIM>::initialiseAuxiliaryMesh() {
   if (this->auxiliary_mesh != NULL) {
     delete this->auxiliary_mesh;
   }
   this->auxiliary_mesh = new Mesh<DIM>(this->refinement_points.getVector(), 0);
   this->auxiliary_mesh->delaunay(std::numeric_limits<Real>::max(), this->qh_flags);
   this->auxiliary_mesh->initInterpolation();
 }
 
 template <Uint DIM>
 Mesh<DIM> & Domain<DIM>::getMesh(Real quality) throw(std::string) {
   if ((!this->main_meshed) || (quality!=this->mesh_quality)) {
     this->main_mesh = new Mesh<DIM>(this->points.getVector(),
                                     this->tol,
                                     this->atomistic_skinned);
     try {
       this->main_mesh->delaunay(quality, this->qh_flags);
     } catch (std::string & error) {
       this->tree->dump();
       this->points.dump_paraview();
       throw("affen_fehler " + error);
     }
     this->main_mesh->cleanup(*this->atomistic);
     this->main_meshed = true;
     this->fill_atoms();
     this->fill_coupling_atoms();
     this->mesh_quality = quality;
   }
 
   return *this->main_mesh;
 }
 
 template <Uint DIM>
 void Domain<DIM>::splitDegenerates(const Real & radius_threshold) {
   this->main_mesh->splitDegenerates(radius_threshold, this->points);
   delete this->main_mesh;
   this->main_meshed = false;
 }
 
 template <Uint DIM>
 Mesh<DIM> & Domain<DIM>::getAuxiliaryMesh() {
   if (this->auxiliary_mesh == NULL) {
     this->initialiseAuxiliaryMesh();
   }
   return *this->auxiliary_mesh;
 }
 
 template <Uint DIM>
 Tree<DIM> & Domain<DIM>::getTree() {
   return *this->tree;
 }
 
 template <Uint DIM>
 bool Domain<DIM>::inFemDomain(const PointRef<DIM> & point) {
   return this->main_mesh->containsExactNode(point);
 }
 
 template <Uint DIM>
 void Domain<DIM>::fill_atoms() {
   this->atoms.clear();
   for (Uint point_id = 0; point_id < this->points.getSize() ; ++point_id) {
     PointRef<DIM> point = this->points.getPoint(point_id);
     if (this->refined->is_inside(point, this->tol)) {
       this->atoms.addPoint(point);
     }
   }
 }
 template <Uint DIM>
 void Domain<DIM>::fill_coupling_atoms() {
   this->coupling_atom_ids.clear();
   for (Uint atom_id = 0 ;  atom_id < this->atoms.getSize() ; ++atom_id) {
     const PointRef<DIM> & point = this->atoms.getConstPoint(atom_id);
     if (!this->atomistic_skinned->is_inside(point, this->tol) &&
       this->inFemDomain(point)) {
       this->coupling_atom_ids.push_front(atom_id);
     }
   }
 }
 
 template <Uint DIM>
 void Domain<DIM>::dump_atoms_lammps(std::string filename) {
   if (this->atoms.getSize() == 0) {
     throw(std::runtime_error("there are no atoms"));
   }
   this->atoms.dump_lammps(filename);
 }
 
 template <Uint DIM>
 void Domain<DIM>::dump_atoms_paraview(std::string filename) {
   if (this->atoms.getSize() == 0) {
     throw(std::runtime_error("there are no atoms"));
   }
   this->atoms.dump_paraview(filename);
 }
 
 template <Uint DIM>
 const PointContainer<DIM>  & Domain<DIM>::getAtoms() {
   if (this->atoms.getSize() == 0) {
     throw(std::runtime_error("there are no atoms. Did you forget to call Domain<DIM>::fill_atoms() on this domain?"));
   }
   return this->atoms;
 }
 
 template <Uint DIM>
 bool Domain<DIM>::inDomain(PointRef<DIM> & point) {
   InsideObject eval = this->overall->from_border(point, this->tol);
   if (!eval.is_inside) {
     // we're out of the domain
     return false;
   } else {
     if (this->boundaries_special_treatment) {
       int int_refinement = this->interpolate(point);
       Real refinement = int_refinement
         * this->lattice->getRepresentativeConstant();
       return max(-eval.distance, 0.) >= refinement*.75; // heuristic value
     } else {
       return true;
     }
   }
 }
 
 template<Uint DIM>
 void Domain<DIM>::compute_coupling_atoms(PointContainer<DIM> & interface,
                                          PointContainer<DIM> & pad,
                                          const Geometry<DIM> & atomistic) {
   if (this->atoms.getSize() == 0) {
     throw("there are no atoms");
   }
   if (!this->main_meshed) {
     throw std::runtime_error("this domain has not been meshed. Did you call getMesh(quality)?");
   }
   if (&atomistic == NULL) {
     this->main_mesh->compute_interface_atoms(interface, *this->atomistic,
                                             this->tol);
   } else {
     this->main_mesh->compute_interface_atoms(interface, atomistic,
                                             this->tol);
   }
   // now, clean the pad atoms to get rid of the interface atoms
   for (Uint interface_index = 0 ;  interface_index < interface.getSize() ;
        ++interface_index) {
     PointRef<DIM> interface_point = interface.getPoint(interface_index);
     auto pad_index = this->coupling_atom_ids.begin();
     auto pad_before_index = this->coupling_atom_ids.before_begin();
     auto pad_end = this->coupling_atom_ids.end();
     bool searching = true;
     for (; pad_index != pad_end && searching; ++pad_index, ++pad_before_index) {
       PointRef<DIM> pad_point = this->atoms.getPoint(*pad_index);
       if (pad_point == interface_point) {
         this->coupling_atom_ids.erase_after(pad_before_index);
         searching=false;
       }
     }
   }
   auto pad_index = this->coupling_atom_ids.begin();
   auto pad_end = this->coupling_atom_ids.end();
   for(;pad_index != pad_end; ++pad_index) {
     pad.addPoint(this->atoms.getPoint(*pad_index));
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template<Uint DIM>
 void Domain<DIM>::insertCustomPoint(const PointRef<DIM> & point) {
   this->points.addPoint(point);
   if (this->main_meshed) {
     throw std::runtime_error("This domain was aleady meshed. The new point was inserted, but may just have fucked up the domain"
                              );
   }
 }
 
 template class Domain<twoD>;
 template class Domain<threeD>;