src/tetmesh/looseoctree.h

Go to the documentation of this file.
00001 /***************************************************************************
00002                           looseoctree.h  -  description
00003                              -------------------
00004     begin                : Wed Feb 25 2004
00005     copyright            : (C) 2004 by Roman Geus
00006     email                : roman.geus@psi.ch
00007 ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 
00018 #ifndef LOOSEOCTREE_H
00019 #define LOOSEOCTREE_H
00020 
00021 #include <vector>
00022 #include "octant.h"
00023 
00024 using namespace std;
00025 
00026 namespace mesh {
00027 
00037 template <typename NodeType>
00038 class LooseOctree {
00039 public: // typedefs
00042     typedef bool (NodeType::*filter)(const Vector3&) const;
00043 public:
00047     LooseOctree(const Vector3& center, double length, int max_depth);
00050     ~LooseOctree();
00055     Octant<NodeType>* add_node(NodeType* node, Octant<NodeType>* start=0, int depth=0);
00064     void find_nodes_by_point(const Vector3& coord,
00065                              std::vector<NodeType *>& node_list,
00066                              int& nof_visited,
00067                              LooseOctree::filter filter,
00068                              Octant<NodeType>* octant=0);
00077     NodeType* find_one_node_by_point(const Vector3& coord,
00078                                        int& nof_visited,
00079                                        LooseOctree::filter filter,
00080                                        Octant<NodeType>* octant=0);
00081      
00085     void get_statistics(int& max_depth,
00086                         double& avg_node_depth,
00087                         int& nof_octants,
00088                         int& nof_nodes);
00089 
00095     void export_vtk();
00098     void export_vtk_level(int level, std::string filename);
00099 
00100 protected:
00101     void get_statistics_recursive(Octant<NodeType>* octant,
00102                                   int depth,
00103                                   int& max_depth,
00104                                   double& sum_depth,
00105                                   int& nof_octants,
00106                                   int& nof_nodes);
00107 
00108     void export_vtk_recursive(Octant<NodeType>* octant,int depth,
00109                               int x,
00110                               int y,
00111                               int z,
00112                               int target_depth,
00113                               std::vector<int>& cells);
00114 
00115 private:
00117     Octant<NodeType>* _root;
00119     int _max_depth;
00120 };
00121   
00122 
00123 template <typename NodeType>
00124 LooseOctree<NodeType>::LooseOctree(const Vector3& center, double length, int max_depth)
00125     : _max_depth(max_depth)
00126 {
00127     _root = new Octant<NodeType>(center, length);
00128 }
00129 
00130 template <typename NodeType>
00131 LooseOctree<NodeType>::~LooseOctree() {
00132     delete _root;
00133     _root = 0;
00134 }
00135 
00136 template <typename NodeType>
00137 Octant<NodeType>* LooseOctree<NodeType>::add_node(NodeType* node, Octant<NodeType>* octant, int depth) {
00138     // root is the default octant
00139     if (octant == 0)
00140         octant = _root;
00141       
00142     Vector3 node_extent, node_center;
00143     node->get_bounding_box(node_center, node_extent);
00144 
00145     // if the octant is twice as big as the node,
00146     // we will add it to a child.
00147     if ((depth < _max_depth ) &&
00148         octant->_length > 2*node_extent.x &&
00149         octant->_length > 2*node_extent.y &&
00150         octant->_length > 2*node_extent.z
00151         ) {
00152         // find child octant containing the mid-point of box
00153         int child_idx = octant->get_child_index_for_point(node_center);
00154 
00155         // If Octant does not yet exist: create new one
00156         if (octant->_children[child_idx] == 0)
00157             octant->_children[child_idx] = new Octant<NodeType>(octant, child_idx);
00158 
00159         // call myself recursively
00160         return add_node(node, octant->_children[child_idx], depth + 1);
00161 
00162     } else {
00163         // found destination octant: add node
00164         octant->add_node(node);
00165         return octant;
00166     }
00167 }
00168 
00169 template <typename NodeType>
00170 void LooseOctree<NodeType>::find_nodes_by_point(const Vector3& coord,
00171                                       std::vector<NodeType *>& nodes,
00172                                       int& nof_visited,
00173                                       LooseOctree<NodeType>::filter filter,
00174                                       Octant<NodeType>* octant) {
00175     // root is the default octant
00176     if (octant == 0)
00177         octant = _root;
00178     // return if coord is outside octant (including overlapping region)
00179     Vector3 dist = coord - octant->_center;
00180     double ext_length = 2.0 * octant->_length;
00181     if (fabs(dist.x) > ext_length ||
00182         fabs(dist.y) > ext_length ||
00183         fabs(dist.z) > ext_length)
00184         return;
00185     // check nodes stored in current Octant
00186     typename std::vector<NodeType *>::iterator it = octant->_nodes.begin();
00187     while (it != octant->_nodes.end()) {
00188         NodeType* node = *it;
00189         if ((node->*filter)(coord))
00190             nodes.push_back(node);
00191         nof_visited ++;
00192         it ++;
00193     }
00194     // descend to children
00195     for (int i = 0; i < 8; i ++)
00196         if (octant->_children[i])
00197             find_nodes_by_point(coord, nodes, nof_visited, filter, octant->_children[i]);
00198 }
00199 
00200 template <typename NodeType>
00201 NodeType* LooseOctree<NodeType>::find_one_node_by_point(const Vector3& coord,
00202                                                 int& nof_visited,
00203                                                 LooseOctree<NodeType>::filter filter,
00204                                                 Octant<NodeType>* octant) {
00205     NodeType* node;
00206     
00207     // root is the default octant
00208     if (octant == 0)
00209         octant = _root;
00210     // return if coord is outside octant (including overlapping region)
00211     Vector3 dist = coord - octant->_center;
00212     double ext_length = 2.0 * octant->_length;
00213     if (fabs(dist.x) > ext_length ||
00214         fabs(dist.y) > ext_length ||
00215         fabs(dist.z) > ext_length)
00216         return 0;
00217     // check nodes stored in current Octant
00218     typename std::vector<NodeType *>::iterator it = octant->_nodes.begin();
00219     while (it != octant->_nodes.end()) {
00220         node = *it;
00221         nof_visited ++;
00222         it ++;
00223         if ((node->*filter)(coord))
00224             return node;
00225     }
00226     // descend to children
00227     for (int i = 0; i < 8; i ++) {
00228         if (octant->_children[i]) {
00229             node = find_one_node_by_point(coord, nof_visited, filter, octant->_children[i]);
00230             // if a node was found abort descend to other children
00231             if (node)
00232                 return node;
00233         }
00234     }
00235     // Nothing found, return 0
00236     return 0;
00237 }
00238 
00239 template <typename NodeType>
00240 void LooseOctree<NodeType>::get_statistics(int& max_depth,
00241                                  double& avg_node_depth,
00242                                  int& nof_octants,
00243                                  int& nof_nodes) {
00244     Octant<NodeType>* octant = _root;
00245     int depth = 0;
00246     double sum_depth = 0.0;
00247 
00248     max_depth = 0;
00249     nof_octants = 0;
00250     nof_nodes = 0;
00251     get_statistics_recursive(octant, depth, max_depth, sum_depth,
00252                              nof_octants, nof_nodes);
00253     avg_node_depth = sum_depth / nof_nodes;
00254 }
00255   
00256 template <typename NodeType>
00257 void LooseOctree<NodeType>::get_statistics_recursive(Octant<NodeType>* octant,
00258                                            int depth,
00259                                            int& max_depth,
00260                                            double& sum_depth,
00261                                            int& nof_octants,
00262                                            int& nof_nodes) {
00263     // update statistics
00264     if (depth > max_depth)
00265         max_depth = depth;
00266     sum_depth += octant->_nodes.size() * depth;
00267     nof_octants ++;
00268     nof_nodes += octant->_nodes.size();
00269 
00270     // descend to children
00271     for (int i = 0; i < 8; i ++)
00272         if (octant->_children[i])
00273             get_statistics_recursive(octant->_children[i], depth+1, max_depth, sum_depth,
00274                                      nof_octants, nof_nodes);
00275 }
00276 
00277 template <typename NodeType>
00278 void LooseOctree<NodeType>::export_vtk() {
00279     // get max depth
00280     int max_depth, nof_octants, nof_nodes;
00281     double avg_node_depth;
00282     get_statistics(max_depth, avg_node_depth, nof_octants, nof_nodes);
00283 
00284     // export all levels to separate VTK files
00285     for (int level = 0; level <= max_depth; ++ level) {
00286         std::ostringstream buf;
00287         buf << "octree_l" << level << ".vtk";
00288         export_vtk_level(level, buf.str());
00289     }
00290 }
00291 
00292 template <typename NodeType>
00293 void LooseOctree<NodeType>::export_vtk_level(int level, std::string filename) {
00294     std::vector<int> cells;
00295     // recursive descend
00296     export_vtk_recursive(_root, 0, 0, 0, 0, level, cells);
00297     int nof_octants = cells.size() / 9;
00298 
00299     std::ofstream of;
00300     of.open(filename.c_str());
00301     rAssert(of.is_open());
00302 
00303     // VTK header
00304     of.precision(6);
00305     of << scientific;
00306     of << "# vtk DataFile Version 2.0" << endl;
00307     of << "generated using VtkExport::export_eigenfields" << endl;
00308     of << "ASCII" << endl << endl;
00309     of << "DATASET UNSTRUCTURED_GRID" << endl;
00310 
00311     // export all points in level, (2**depth + 1)**3 points
00312     int n = (1 << level) + 1;
00313     Vector3 origin = _root->_center;
00314     origin -= _root->_length;
00315     double side = 2.0*_root->_length / (1 << level);
00316     of << "POINTS " << n*n*n << " float" << endl;
00317 
00318     for (int z = 0; z < n; ++ z)
00319         for (int y = 0; y < n; ++ y)
00320             for (int x = 0; x < n; ++ x)
00321                 of << origin.x + x*side << " "
00322                    << origin.y + y*side << " "
00323                    << origin.z + z*side << endl;
00324     of << endl;
00325 
00326     // cells
00327     of << "CELLS " << nof_octants << " " << cells.size() << endl;
00328     std::vector<int>::iterator it = cells.begin();
00329     for (int i = 0; it != cells.end(); ++ it, ++ i) {
00330         of << *it;
00331         if (i % 9 == 8)
00332             of << endl;
00333         else
00334             of << " ";
00335     }
00336     of << endl;
00337 
00338     // cell types
00339     of << "CELL_TYPES " << nof_octants << endl;
00340     for (int i = 0; i < nof_octants; ++ i)
00341         of << "11" << endl;
00342     of << endl;
00343 
00344     // cell data
00345     of << "CELL_DATA " << nof_octants << endl;
00346     of << "SCALARS cell_scalars int 1" << endl;
00347     of << "LOOKUP_TABLE default" << endl;
00348     for (int i = 0; i < nof_octants; ++ i)
00349         of << "0" << endl;
00350     of << endl;
00351 }
00352 
00353 template <typename NodeType>
00354 void LooseOctree<NodeType>::export_vtk_recursive(Octant<NodeType>* octant,
00355                                        int depth,
00356                                        int x,
00357                                        int y,
00358                                        int z,
00359                                        int target_depth,
00360                                        std::vector<int>& cells) {
00361     if (depth == target_depth) {
00362         // export octant
00363         int n = (1 << depth) + 1;
00364         cells.push_back(8);
00365         cells.push_back((x+0) + (y+0)*n + (z+0)*n*n);
00366         cells.push_back((x+1) + (y+0)*n + (z+0)*n*n);
00367         cells.push_back((x+0) + (y+1)*n + (z+0)*n*n);
00368         cells.push_back((x+1) + (y+1)*n + (z+0)*n*n);
00369         cells.push_back((x+0) + (y+0)*n + (z+1)*n*n);
00370         cells.push_back((x+1) + (y+0)*n + (z+1)*n*n);
00371         cells.push_back((x+0) + (y+1)*n + (z+1)*n*n);
00372         cells.push_back((x+1) + (y+1)*n + (z+1)*n*n);
00373     } else {
00374         // descend to children
00375         for (int i = 0; i < 8; i ++) {
00376             if (octant->_children[i])
00377                 export_vtk_recursive(octant->_children[i],
00378                                      depth+1,
00379                                      x + ((i & 0x1) == 0x1) * (1 << (target_depth - depth - 1)),
00380                                      y + ((i & 0x2) == 0x2) * (1 << (target_depth - depth - 1)),
00381                                      z + ((i & 0x4) == 0x4) * (1 << (target_depth - depth - 1)),
00382                                      target_depth,
00383                                      cells);
00384         }
00385     }
00386 }
00387 
00388 } // namespace mesh
00389 
00390 #endif

Generated on Fri Oct 26 13:35:12 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7