meshconv v3.11

src/backgrid.cc

@@ -26,9 +26,9 @@ double BackGrid::search_distance = -1.0;
 void BackGrid::debug_output_dimensions(){
   int x, y, z, sum=0, cnt=0;
 
-  for(x=0; x<segments[0]; ++x){
-    for(y=0; y<segments[1]; ++y){
-      for(z=0; z<segments[2]; ++z){
+  for(x = 0; x < segments[0]; ++x){
+    for(y = 0; y < segments[1]; ++y){
+      for(z = 0; z < segments[2]; ++z){
         sum += grid[x][y][z].size();
         ++cnt;
       }
@@ -44,7 +44,6 @@ void BackGrid::debug_output_dimensions(){
 void BackGrid::init(Mesh *m, Mesh *n, uchar operation_index_mode, bool no_scaling){
   int x, y;
   Vertex dist;
-  double distpower = 1.0;
   double max_spread = 0.0;
   int seg_deg;
 
@@ -76,18 +75,18 @@ void BackGrid::init(Mesh *m, Mesh *n, uchar operation_index_mode, bool no_scalin
 
   //calculate minimal and maximal coordinates of backgrid
   determine_coordinate_extrema(*m, minc, maxc);
-  for(int p=0; p<m->dim_of_world; ++p){
+  for(int p = 0; p < m->dim_of_world; ++p){
     if(maxc[p] - minc[p] > max_spread) max_spread = maxc[p] - minc[p];
   }
   if(n != NULL){
-    for(x=0; x<n->numvertices; ++x){
-      for(y=0; y<m->dim_of_world; ++y){
+    for(x = 0; x < n->numvertices; ++x){
+      for(y = 0; y < m->dim_of_world; ++y){
         if(n->vertices[x][y] < minc[y]) minc[y] = n->vertices[x][y];
         if(n->vertices[x][y] > maxc[y]) maxc[y] = n->vertices[x][y];
       }
     }
   }
-  for(y=0; y<m->dim_of_world; ++y){
+  for(y = 0; y < m->dim_of_world; ++y){
     minc[y] = minc[y] - coarse_epsilon;
     maxc[y] = maxc[y] + coarse_epsilon;
   }
@@ -95,30 +94,31 @@ void BackGrid::init(Mesh *m, Mesh *n, uchar operation_index_mode, bool no_scalin
   //calculate sizes of segments
   if(no_scaling){
     //boxes must be cubes (all edges of same length)
-    distpower = 0.0;
+    double maxdist = 0.0;
     y = 0;
-    for(x=0; x<m->dim_of_world; ++x){
+    for(x = 0; x < m->dim_of_world; ++x){
       dist[x] = maxc[x] - minc[x];
-      if(dist[x] > distpower){
-        distpower = dist[x];
+      if(dist[x] > maxdist){
+        maxdist = dist[x];
         y = x;
       }
     }
-    for(x=0; x<m->dim_of_world; ++x){
-      seg_size[x] = distpower / (double)(seg_deg);
-      segments[x] = int((dist[x] / distpower) * seg_deg);
+    for(x = 0; x < m->dim_of_world; ++x){
+      seg_size[x] = maxdist / (double)(seg_deg);
+      segments[x] = int((dist[x] / maxdist) * seg_deg);
       if(dist[x] < dist[y]) ++segments[x];
     }
     segments[y] = seg_deg;
     min_seg_size = seg_size[0];
   }else{
     //boxes need not be cubes
+    double distpower = 1.0;
     min_seg_size = 9e100;
-    for(x=0; x<m->dim_of_world; ++x){
+    for(x = 0; x < m->dim_of_world; ++x){
       dist[x] = maxc[x] - minc[x];
       if(dist[x] > epsilon && dist[x] < min_seg_size) min_seg_size = dist[x];
     }
-    for(x=0; x<m->dim_of_world; ++x){
+    for(x = 0; x < m->dim_of_world; ++x){
       if(dist[x] <= epsilon) dist[x] = min_seg_size;
       else dist[x] = 1.0 / dist[x];
       distpower *= dist[x];
@@ -127,7 +127,7 @@ void BackGrid::init(Mesh *m, Mesh *n, uchar operation_index_mode, bool no_scalin
 
     double tmp_segments[3];
     int latest_zero_dimension = 0, latest_non_zero_dimension = 0, num_zero_dimensions = 0;
-    for(x=0; x<m->dim_of_world; ++x){
+    for(x = 0; x < m->dim_of_world; ++x){
       tmp_segments[x] = int(distpower / dist[x]);
       if(tmp_segments[x] < 1.0){
         latest_zero_dimension = x;
@@ -166,13 +166,13 @@ void BackGrid::init(Mesh *m, Mesh *n, uchar operation_index_mode, bool no_scalin
                              * tmp_segments[latest_non_zero_dimension];
       }
     }
-    for(x=0; x<m->dim_of_world; ++x){
+    for(x = 0; x < m->dim_of_world; ++x){
       segments[x] = tmp_segments[x];
       if(segments[x] <= 0) segments[x] = 1;
     }
 
     min_seg_size = 9e100;
-    for(x=0; x<m->dim_of_world; ++x){
+    for(x = 0; x < m->dim_of_world; ++x){
       seg_size[x] = (maxc[x] - minc[x]) / (double)segments[x];
       if(min_seg_size > seg_size[x]) min_seg_size = seg_size[x];
     }
@@ -195,9 +195,9 @@ void BackGrid::init(Mesh *m, Mesh *n, uchar operation_index_mode, bool no_scalin
   //create data structure for backgrid
   if(m->dim_of_world == 2) segments[2] = 1;
   grid.resize(segments[0]);
-  for(x=0; x<segments[0]; ++x){
+  for(x = 0; x < segments[0]; ++x){
     grid[x].resize(segments[1]);
-    for(y=0; y<segments[1]; ++y){
+    for(y = 0; y < segments[1]; ++y){
       grid[x][y].resize(segments[2]);
     }
   }
@@ -870,20 +870,20 @@ void BackGrid::attach_mesh_3d(){
   #define DEBUG_REFINE_BACKGRID_3D 0
 
   #if DEBUG_REFINE_BACKGRID_3D != 0
-  rebuild_element_indices(*n);
-  int debug_element = 11;
-  vector<vector<vector<short int> > > z_bt;
-  z_bt.resize(segments[0]);
-  for(x=0; x<segments[0]; ++x){
-    z_bt[x].resize(segments[1]);
-    for(y=0; y<segments[1]; ++y){
-      z_bt[x][y].resize(segments[2]);
-      for(z=0; z<segments[2]; ++z) z_bt[x][y][z] = 0;
+    rebuild_element_indices(*n);
+    int debug_element = 11;
+    vector<vector<vector<short int> > > z_bt;
+    z_bt.resize(segments[0]);
+    for(x=0; x<segments[0]; ++x){
+      z_bt[x].resize(segments[1]);
+      for(y=0; y<segments[1]; ++y){
+        z_bt[x][y].resize(segments[2]);
+        for(z=0; z<segments[2]; ++z) z_bt[x][y][z] = 0;
+      }
     }
-  }
   #endif
 
-  for(it=n->elements.begin(); it!=n->elements.end(); ++it){
+  for(it = n->elements.begin(); it != n->elements.end(); ++it){
     //create operation_index
     if(operation_index_mode == oim_create){
       if(it->operation_index != NULL){
@@ -913,9 +913,9 @@ void BackGrid::attach_mesh_3d(){
     if(cz2 >= segments[2]) cz2 = segments[2] - 1;
 
     //test for all these segments if an intersection with 'it' occurs
-    for(z=cz1; z<=cz2; ++z){
-      for(y=cy1; y<=cy2; ++y){
-        for(x=cx1; x<=cx2; ++x){
+    for(z = cz1; z <= cz2; ++z){
+      for(y = cy1; y <= cy2; ++y){
+        for(x = cx1; x <= cx2; ++x){
           //get coordinates of the current segment
           sx1 = x * seg_size[0] + minc[0];
           sx2 = sx1 + seg_size[0];

src/boundary_operations.cc

@@ -240,7 +240,6 @@ vector<vector<vector<int> > >* box_boundary_periodic(Mesh &m, int *dimensions){
         bfaces[e].insert(new_face, new_value);
       }else{
         vector<int> &other_value = find_result[0]->second;
-        //! correct ordering of local indices
         for(int i = 1; i <= m.dim_of_elements; ++i){
           for(int j = 1; j <= m.dim_of_elements; ++j){
             int w = 0;

src/commands.cc

@@ -101,6 +101,28 @@ void call_assign(vector<string> &args, bool test){
   }
 }
 
+//----------------------------------------< call_arh_compression >
+void call_arh_compression(vector<string> &args, bool test){
+  if(!test){
+    string task_description = "arh_compression";
+    int prompt_length = task_description.size() + 5;
+    if(undefined_variable_warning != ""){
+      cout << task_description << fill_dots(prompt_length, task_description_length)
+           << undefined_variable_warning << endl;
+    }
+
+    check_argument_count(args, 1);
+    use_compression = ival(args[0]);
+    #ifdef NO_COMPRESSION
+      if(use_compression != 0){
+        error("meshconv was not compiled with compression-support!"
+                       " Try compilation without -DNO_COMPRESSION.");
+      }
+    #endif
+    if(use_compression > 2) error("Unknown compression-type.");
+  }
+}
+
 //----------------------------------------< call_arh_info >
 void call_arh_info(vector<string> &args, bool test){
   string task_description = (args.empty())? "" : "arh info for \"" + file_name(args[0]) + "\"";
@@ -1171,20 +1193,20 @@ void call_detect_holes(vector<string> &args, bool test){
       cout << task_description + fill_dots(prompt_length, task_description_length) << flush;
     }
     if(undefined_variable_warning != "") cout << undefined_variable_warning << flush;
-    check_argument_count(args, 2);
-    Mesh *m = lookup_mesh(args[1]);
+    check_argument_count(args, 1, 2);
+    Mesh *m = lookup_mesh(args[args.size() - 1]);
 
     int cnt = 0;
     if(m->dim_of_world == 3 && m->dim_of_elements == 2){
       map<int, vector<pair<int, int> > > holes;
       cnt = detect_holes(*m, holes);
-      output_holes(args[0], *m, holes);
+      if(args.size() == 2) output_holes(args[0], *m, holes);
     }else{
       error("Hole-detection is implemented for 3D-surfaces only.");
     }
 
     if(cnt == 0 || verbosity == 1) cout << "ok" << endl;
-    else if(verbosity >= 2) cout << cnt << " holes" << endl;
+    else if(verbosity >= 2 || args.size() == 1) cout << cnt << " holes" << endl;
   }else{
     if(task_description_length < prompt_length) task_description_length = prompt_length;
   }
@@ -3265,28 +3287,6 @@ void call_set_max_threads(vector<string> &args, bool test){
   }
 }
 
-//----------------------------------------< call_set_output_compression >
-void call_set_output_compression(vector<string> &args, bool test){
-  if(!test){
-    string task_description = "set_output_compression";
-    int prompt_length = task_description.size() + 5;
-    if(undefined_variable_warning != ""){
-      cout << task_description << fill_dots(prompt_length, task_description_length)
-           << undefined_variable_warning << endl;
-    }
-
-    check_argument_count(args, 1);
-    use_compression = ival(args[0]);
-    #ifdef NO_COMPRESSION
-      if(use_compression != 0){
-        error("meshconv was not compiled with compression-support!"
-                       " Try compilation without -DNO_COMPRESSION.");
-      }
-    #endif
-    if(use_compression > 2) error("Unknown compression-type.");
-  }
-}
-
 //----------------------------------------< call_set_precision >
 void call_set_precision(vector<string> &args, bool test){
   if(!test){
@@ -3563,12 +3563,18 @@ void call_split_edges_element_count(vector<string> &args, bool test){
       cout << task_description + fill_dots(prompt_length, task_description_length) << flush;
     }
     if(undefined_variable_warning != "") cout << undefined_variable_warning << flush;
-    check_argument_count(args, 2);
+    check_argument_count_at_least(args, 2);
     Mesh *m = lookup_mesh(args[0]);
     int min_element_count = ival(args[1]);
+    vector<DataVector*> vecs, el_vecs;
+    for(int i = 2, i_end = args.size(); i < i_end; ++i){
+      if(args[i][0] == '^') el_vecs.push_back(lookup_vector(args[i]));
+      else vecs.push_back(lookup_vector(args[i]));
+    }
+
 
     if((m->dim_of_world == 3 || m->dim_of_world == 2) && m->dim_of_elements == 2){
-      split_edges_3d(*m, -1.0, min_element_count);
+      split_edges_3d(*m, vecs, el_vecs, -1.0, min_element_count);
     }else{
       error("Invalid dimensionality of input mesh.");
     }
@@ -3589,12 +3595,17 @@ void call_split_edges_max_length(vector<string> &args, bool test){
       cout << task_description + fill_dots(prompt_length, task_description_length) << flush;
     }
     if(undefined_variable_warning != "") cout << undefined_variable_warning << flush;
-    check_argument_count(args, 2);
+    check_argument_count_at_least(args, 2);
     Mesh *m = lookup_mesh(args[0]);
     double max_edge_length = dval(args[1]);
+    vector<DataVector*> vecs, el_vecs;
+    for(int i = 2, i_end = args.size(); i < i_end; ++i){
+      if(args[i][0] == '^') el_vecs.push_back(lookup_vector(args[i]));
+      else vecs.push_back(lookup_vector(args[i]));
+    }
 
     if((m->dim_of_world == 3 || m->dim_of_world == 2) && m->dim_of_elements == 2){
-      split_edges_3d(*m, max_edge_length, INT_MAX);
+      split_edges_3d(*m, vecs, el_vecs, max_edge_length, INT_MAX);
     }else{
       error("Invalid dimensionality of input mesh.");
     }
@@ -4430,12 +4441,14 @@ void call_voxels_to_mesh(vector<string> &args, bool test){
       cout << task_description + fill_dots(prompt_length, task_description_length) << flush;
     }
     if(undefined_variable_warning != "") cout << undefined_variable_warning << flush;
-    check_argument_count(args, 1);
+    check_argument_count(args, 1, 2);
     if(!voxel_manager.has_data()) error("No voxel-data loaded!");
     Mesh *m = lookup_mesh(args[0], true);
     delete_refinement_manager(args[0]);
+    double voxel_size = -1.0;
+    if(args.size() == 2) voxel_size = dval(args[1]);
 
-    voxel_manager.voxels_to_mesh(*m);
+    voxel_manager.voxels_to_mesh(*m, voxel_size);
 
     if(verbosity >= 1) cout << "ok" << endl;
   }else{
@@ -4955,6 +4968,7 @@ void execute_next_command(string &cmdline, bool test){
   if(test) ++command_count_test; else ++command_count_exe;
   if(skip_mode == skip_off){
     if(srcmd == "about") call_about(test);
+    else if(srcmd == "arh_compression") call_arh_compression(cmdargs, test);
     else if(srcmd == "arh_info") call_arh_info(cmdargs, test);
     else if(srcmd == "arh_type") call_arh_type(cmdargs, test);
     else if(srcmd == "assign") call_assign(cmdargs, test);
@@ -5067,7 +5081,6 @@ void execute_next_command(string &cmdline, bool test){
     else if(srcmd == "set_line_triangle_intersection")
                                        call_set_line_triangle_intersection(cmdargs, test);
     else if(srcmd == "set_max_threads") call_set_max_threads(cmdargs, test);
-    else if(srcmd == "set_output_compression") call_set_output_compression(cmdargs, test);
     else if(srcmd == "set_precision") call_set_precision(cmdargs, test);
     else if(srcmd == "set_refine_mode") call_set_refine_mode(cmdargs, test);
     else if(srcmd == "set_temporary_file_folder") call_set_temporary_file_folder(cmdargs, test);

src/filereaders.cc

@@ -2051,7 +2051,7 @@ void recurse_refine(Mesh &m, list<Element>::iterator &e_it, vector<uint64_t> &ve
     m.hierarchy.insert(make_pair(fi, entry_value));
 
     //setup right child element
-    type = level % 3 + RefinementManager::base_level_type;
+    type = (level + RefinementManager::base_level_type) % 3;
     if(m.dim_of_world == 2){
       v[0] = e_it->v[1];
       v[1] = e_it->v[2];
@@ -4125,11 +4125,12 @@ void create_bt_mesh(Mesh &m, int dim, int size_x, int size_y, int size_z,
         for(zb = 0; zb <= 1; ++zb){
           for(yb = 0; yb <= 1; ++yb){
             for(xb = 0; xb <= 1; ++xb){
-              if(x-xb>=0 && x-xb<size_x && y-yb>=0 && y-yb<size_y && z-zb>=0 && z-zb<size_z){
-                if(boxes[x-xb][y-yb][z-zb] != '0'){
-                  current_box_type = box_types.find(boxes[x-xb][y-yb][z-zb]);
+              if(x - xb >= 0 && x - xb < size_x && y - yb >= 0 && y - yb < size_y
+                                                    && z - zb >= 0 && z-zb<size_z){
+                if(boxes[x - xb][y - yb][z - zb] != '0'){
+                  current_box_type = box_types.find(boxes[x - xb][y - yb][z - zb]);
                   if(current_box_type == box_types.end()){
-                    string errtmp(1, boxes[x-xb][y-yb][z-zb]);
+                    string errtmp(1, boxes[x - xb][y - yb][z - zb]);
                     error("Undefined box-type '" + errtmp +
                           "'. (Check for misspelled keywords in header!)");
                   }

src/interface.cc

@@ -738,7 +738,7 @@ void replace_constants(string &str){
   constants.insert(make_pair("NON_HIERARCHICAL", "0")); //make_periodic
   constants.insert(make_pair("HIERARCHICAL", "1"));
 
-  for(map<string, string>::iterator it=constants.begin(); it!=constants.end(); ++it){
+  for(map<string, string>::iterator it = constants.begin(); it != constants.end(); ++it){
     pos = 0;
     do{
       pos = str.find(it->first, pos);

src/mesh_maintenance.cc

@@ -241,7 +241,7 @@ int64_t forest_index(Mesh &m, int64_t parent_index){//$ACL_exclude
     tmp >>= 1;
     ++level;
   }while(tmp > 0);
-  next_index = int64_t(m.numbaseelements) * ((1 << (level+1)) - 1);
+  next_index = int64_t(m.numbaseelements) * ((1 << (level + 1)) - 1);
   diff = parent_index - m.numbaseelements * ((1 << level) - 1);
   return next_index + 2 * diff;
 }
@@ -253,7 +253,7 @@ int forest_index(Mesh &m, int parent_index){//$ACL_exclude
     tmp >>= 1;
     ++level;
   }while(tmp > 0);
-  next_index = m.numbaseelements * ((1 << (level+1)) - 1);
+  next_index = m.numbaseelements * ((1 << (level + 1)) - 1);
   diff = parent_index - m.numbaseelements * ((1 << level) - 1);
   return next_index + 2 * diff;
 }

src/mesh_operations.cc

@@ -1474,7 +1474,8 @@ struct SplitEdge{
 };
 
 //----------------------------------------< split_edges_3d >
-void split_edges_3d(Mesh &m, double max_edge_length, int min_elements){
+void split_edges_3d(Mesh &m, vector<DataVector*> &vecs, vector<DataVector*> &el_vecs,
+                                             double max_edge_length, int min_elements){
   //prepare progress-output
   double progress_step = 0.0;
   int progress_index = 1;
@@ -1483,6 +1484,8 @@ void split_edges_3d(Mesh &m, double max_edge_length, int min_elements){
   bool do_progress = (progress_units > 0 && min_elements < INT_MAX && verbosity >= 2);
   if(do_progress) cout << "analysing: " << flush;
 
+  if(!el_vecs.empty()) rebuild_element_indices(m);
+
   //collect edges
   unsigned j;
   SplitEdge c_edge;
@@ -1526,33 +1529,48 @@ void split_edges_3d(Mesh &m, double max_edge_length, int min_elements){
     current_progress_units = (progress_units == 1)? 10 : progress_units;
     progress_step = double(min_elements - original_elements) / current_progress_units;
     if(progress_step < 1.0) progress_step = 1.0;
-    cout << "done. collapsing: " << flush;
+    cout << "done. splitting: " << flush;
+  }
+
+  //determine component-counts of multi-component-vectors
+  vector<int> vertex_vector_component_sizes(vecs.size());
+  for(size_t vi = 0; vi < vecs.size(); ++vi){
+    if(vecs[vi]->size() == 0 || vecs[vi]->size() % m.numvertices != 0){
+      error("Vertex-vector " + ft(uint64_t(vi + 1)) + " does not fit the mesh ("
+                   + ft(vecs[vi]->size()) + " / " + ft(m.numvertices) + ").");
+    }
+    vertex_vector_component_sizes[vi] = vecs[vi]->size() / m.numvertices;
+  }
+  vector<int> element_vector_component_sizes(el_vecs.size());
+  for(size_t vi = 0; vi < el_vecs.size(); ++vi){
+    if(el_vecs[vi]->size() == 0 || el_vecs[vi]->size() % m.numelements != 0){
+      error("Element-vector " + ft(uint64_t(vi + 1)) + " does not fit the mesh ("
+                   + ft(el_vecs[vi]->size()) + " / " + ft(m.numelements) + ").");
+    }
+    element_vector_component_sizes[vi] = el_vecs[vi]->size() / m.numelements;
   }
 
   //split longest edge until all edges are no longer than 'max_edge_length'
+  map<int, list<Element>::iterator> deleted_elements;
   Vertex nv;
   int sv[3];
   SplitEdge child_edge1, child_edge2;
+  int current_num_elements = m.numelements;
   while(true){
-    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    //! debug (output distance-map)
-    /*!zzz("\ndistance-map:");
-    for(multimap<double, SplitEdge*>::iterator it = distance_map.begin(),
-                                                 end = distance_map.end(); it != end; ++it){
-      zzz(ft(it->first, 7, 4) + "(" + ft(it->second->v[0]) + ", " + ft(it->second->v[1]) + ")");
-    }
-    zzz();*/
-    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
     multimap<double, SplitEdge*>::iterator s_it = distance_map.end();
     --s_it;
-    if(s_it->first <= max_edge_length || m.numelements >= min_elements) break;
+    if(s_it->first <= max_edge_length || current_num_elements >= min_elements) break;
     SplitEdge &edge = *s_it->second;
-    //!zzz("split-edge: (" + ft(edge.v[0]) + ", " + ft(edge.v[1]) + ")"); //!!!!!!!!
 
     //add new vertex
     for(int i = 0; i < 3; ++i) nv[i] = (m.vertices[edge.v[0]][i] + m.vertices[edge.v[1]][i]) / 2.0;
     add_vertex(m, nv);
+    for(size_t vi = 0; vi < vecs.size(); ++vi){
+      for(size_t ci = 0, c_end = vertex_vector_component_sizes[vi]; ci < c_end; ++ci){
+        vecs[vi]->push_back(((*vecs[vi])[edge.v[0] * c_end + ci]
+                           + (*vecs[vi])[edge.v[1] * c_end + ci]) / 2.0);
+      }
+    }
 
     //run through elements at split-edge
     for(int i = 0; i < edge.el_cnt; ++i){
@@ -1567,6 +1585,12 @@ void split_edges_3d(Mesh &m, double max_edge_length, int min_elements){
       }
       sv[2] = edge.v[0];
       add_element(m, sv);
+      ++current_num_elements;
+      for(size_t vi = 0; vi < el_vecs.size(); ++vi){
+        for(size_t ci = 0, c_end = element_vector_component_sizes[vi]; ci < c_end; ++ci){
+          el_vecs[vi]->push_back((*el_vecs[vi])[edge.elements[i]->index]);
+        }
+      }
       //update neighbour-edge (sv[1], sv[2])
       if(sv[1] < sv[2]){
         c_edge.v[0] = sv[1];
@@ -1575,13 +1599,11 @@ void split_edges_3d(Mesh &m, double max_edge_length, int min_elements){
         c_edge.v[0] = sv[2];
         c_edge.v[1] = sv[1];
       }
-      //!zzz("  find (" + ft(c_edge.v[0]) + ", " + ft(c_edge.v[1]) + ")"); //!!!!!!
       set<SplitEdge>::iterator n_it = edges.find(c_edge);
       for(int e = 0; e < n_it->el_cnt; ++e){
         if(n_it->elements[e] == edge.elements[i]){
           list<Element>::iterator t_it = m.elements.end();
           --t_it;
-          //!zzz("  set (" + ft(t_it->v[0]) + ", " + ft(t_it->v[1]) + ", " + ft(t_it->v[2]) + ")"); //!
           *const_cast<list<Element>::iterator*>(&n_it->elements[e]) = t_it;
           break;
         }
@@ -1602,6 +1624,12 @@ void split_edges_3d(Mesh &m, double max_edge_length, int min_elements){
       //add second child-element
       sv[2] = edge.v[1];
       add_element(m, sv);
+      ++current_num_elements;
+      for(size_t vi = 0; vi < el_vecs.size(); ++vi){
+        for(size_t ci = 0, c_end = element_vector_component_sizes[vi]; ci < c_end; ++ci){
+          el_vecs[vi]->push_back((*el_vecs[vi])[edge.elements[i]->index]);
+        }
+      }
       //update neighbour-edge (sv[1], sv[2])
       if(sv[1] < sv[2]){
         c_edge.v[0] = sv[1];
@@ -1610,13 +1638,11 @@ void split_edges_3d(Mesh &m, double max_edge_length, int min_elements){
         c_edge.v[0] = sv[2];
         c_edge.v[1] = sv[1];
       }
-      //!zzz("  find (" + ft(c_edge.v[0]) + ", " + ft(c_edge.v[1]) + ")"); //!!!!!!
       n_it = edges.find(c_edge);
       for(int e = 0; e < n_it->el_cnt; ++e){
         if(n_it->elements[e] == edge.elements[i]){
           list<Element>::iterator t_it = m.elements.end();
           --t_it;
-          //!zzz("  set (" + ft(t_it->v[0]) + ", " + ft(t_it->v[1]) + ", " + ft(t_it->v[2]) + ")"); //!
           *const_cast<list<Element>::iterator*>(&n_it->elements[e]) = t_it;
           break;
         }
@@ -1634,10 +1660,8 @@ void split_edges_3d(Mesh &m, double max_edge_length, int min_elements){
         --child_edge2.elements[1];
       }
 
-      //!zzz("delete element(" + ft(edge.elements[i]->v[0]) + ", " + ft(edge.elements[i]->v[1])
-      //!                                          + ", " + ft(edge.elements[i]->v[2]) + ")"); //!!!
-      m.elements.erase(edge.elements[i]);
-      --m.numelements;
+      deleted_elements[edge.elements[i]->index] = edge.elements[i];
+      --current_num_elements;
 
       //add new edge bisecting deleted element
       c_edge.v[0] = sv[1];
@@ -1646,18 +1670,15 @@ void split_edges_3d(Mesh &m, double max_edge_length, int min_elements){
       c_edge.elements[1] = m.elements.end();
       c_edge.elements[0] = --c_edge.elements[1];
       --c_edge.elements[0];
-      //!zzz("  new inter-edge: (" + ft(c_edge.v[0]) + ", " + ft(c_edge.v[1])); //!!!
       distance_map.insert(pair<double, SplitEdge*>(
                           edge_length_3d(m.vertices[sv[1]], m.vertices[sv[0]]),
                           const_cast<SplitEdge*>(&(*edges.insert(c_edge).first))));
     }
 
     //add child-edges of split-edge
-    //!zzz("  new split-edge1: (" + ft(child_edge1.v[0]) + ", " + ft(child_edge1.v[1])); //!!!
     distance_map.insert(pair<double, SplitEdge*>(
              edge_length_3d(m.vertices[child_edge1.v[0]], m.vertices[child_edge1.v[1]]),
                              const_cast<SplitEdge*>(&(*edges.insert(child_edge1).first))));
-    //!zzz("  new split-edge2: (" + ft(child_edge2.v[0]) + ", " + ft(child_edge2.v[1])); //!!!
     distance_map.insert(pair<double, SplitEdge*>(
              edge_length_3d(m.vertices[child_edge2.v[0]], m.vertices[child_edge2.v[1]]),
                              const_cast<SplitEdge*>(&(*edges.insert(child_edge2).first))));
@@ -1673,12 +1694,16 @@ void split_edges_3d(Mesh &m, double max_edge_length, int min_elements){
     }
   }
 
+  //update element-vectors
+  update_element_vectors_after_deletions(m, el_vecs, deleted_elements);
+
   m.has_valid_element_indices = false;
   m.boundary_checked = false;
 }
 
 //----------------------------------------< eliminate_flat_triangles >
-//Eliminates triangles which have at least one inner angle below the specified 'angle_threshold' by //swapping edges. If element-vectors are specified they are updated according to the mesh-changes
+//Eliminates triangles which have at least one inner angle below the specified 'angle_threshold' by
+//swapping edges. If element-vectors are specified they are updated according to the mesh-changes
 //(i.e. the element-values of the triangles to be eliminated are dropped in favor of the element-
 //values of the neighbour which is involved in the edge-swap).
 void eliminate_flat_triangles(Mesh &m, vector<DataVector*> &vecs, double angle_threshold){
@@ -2452,7 +2477,6 @@ void debug_output_collapsibility(Mesh &m, multimap<double, DecimateEdge*> &dista
 //'max_feature_angle' or lower are preserved if at least 'min_feature_net_size' many such edges are
 //connected to each other. Vertex- and element-vectors contained in 'vecs' and 'el_vecs', respec-
 //tively are updated according to the mesh-changes.
-#define ZIB_COND if(collapse_counter == 33261) //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! $zib
 void decimate_3d(Mesh &m, double min_edge_length, int max_elements, double max_inner_angle,
                  double max_feature_angle, int min_feature_net_size, vector<DataVector*> &vecs,
                                                                    vector<DataVector*> &el_vecs){
@@ -2791,50 +2815,6 @@ void decimate_3d(Mesh &m, double min_edge_length, int max_elements, double max_i
     int v0 = del_edge.v[0];
     int v1 = del_edge.v[1];
 
-    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    //! debug (disallow junctions)
-    //determine vertices of elements adajacent to 'del_edge' that are not part of 'del_edge'
-    /*int del_element_vertices[2] = {-1, -1};
-    int vi;
-    int16_t idj = 0;
-    for(; idj < del_edge.el_cnt; ++idj){
-      for(int i = 0; i < 3; ++i){
-        vi = del_edge.elements[idj]->v[i];