From 15788747b65b12f26d838545dd4ff0497623a39f Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Mon, 4 Nov 2013 19:02:48 +0100 Subject: [PATCH] regions are now properly nested, i.e. ellipsoids are nested (not becoming cuboids at coarser levels), convex hulls are linearly scaled versions of one another... --- mesh.hh | 78 ++++++++++++++-------- plugins/convex_hull.hh | 8 ++- plugins/region_convex_hull.cc | 14 +++- plugins/region_ellipsoid.cc | 118 ++++++++++++++++++++++++++-------- region_generator.cc | 6 +- region_generator.hh | 6 +- 6 files changed, 169 insertions(+), 61 deletions(-) diff --git a/mesh.hh b/mesh.hh index 19eede6..72635ba 100644 --- a/mesh.hh +++ b/mesh.hh @@ -36,10 +36,10 @@ public: : nx_( 0 ), ny_ ( 0 ), nz_( 0 ) { } - refinement_mask( size_t nx, size_t ny, size_t nz ) + refinement_mask( size_t nx, size_t ny, size_t nz, short value = 0. ) : nx_( nx ), ny_( ny ), nz_( nz ) { - mask_.assign( nx_*ny_*nz_, 0 ); + mask_.assign( nx_*ny_*nz_, value ); } refinement_mask( const refinement_mask& r ) @@ -60,12 +60,12 @@ public: return *this; } - void init( size_t nx, size_t ny, size_t nz ) + void init( size_t nx, size_t ny, size_t nz, short value = 0. ) { nx_ = nx; ny_ = ny; nz_ = nz; - mask_.assign( nx_*ny_*nz_, 0 ); + mask_.assign( nx_*ny_*nz_, value ); } const short& operator()( size_t i, size_t j, size_t k ) const @@ -537,7 +537,7 @@ public: m_yoffabs, //!< vector of x-offsets of a level mesh relative to the coarser level m_zoffabs; //!< vector of x-offsets of a level mesh relative to the coarser level - refinement_mask ref_mask; + std::vector< refinement_mask* > m_ref_masks; bool bhave_refmask; protected: @@ -614,15 +614,20 @@ public: m_yoffabs = gh.m_yoffabs; m_zoffabs = gh.m_zoffabs; - ref_mask = gh.ref_mask; + //ref_mask = gh.ref_mask; bhave_refmask = gh.bhave_refmask; + + if( bhave_refmask ) + { + for( size_t i=0; ideallocate(); - } //! free all memory occupied by the grid hierarchy @@ -633,11 +638,14 @@ public: m_pgrids.clear(); std::vector< MeshvarBnd* >().swap( m_pgrids ); - m_xoffabs.clear(); m_yoffabs.clear(); m_zoffabs.clear(); m_levelmin = 0; + + for( size_t i=0; i= levelmin(); --ilevel ) { - xq[0] = (offset_abs(ilevel,0) + i)*dx + 0.5*dx + shift[0]; - for( size_t j=0; jinit( size(ilevel,0), size(ilevel,1), size(ilevel,2), 0.0 ); + + for( size_t i=0; iquery_point( xq ) ) - ref_mask(i,j,k) = true; + xq[1] = (offset_abs(ilevel,1) + j)*dx + 0.5*dx + shift[1]; + for( size_t k=0; kquery_point( xq, ilevel ) ) + (*m_ref_masks[ilevel])(i,j,k) = 1.0;//true; + else + (*m_ref_masks[ilevel])(i,j,k) = 0.0; + } } } } bhave_refmask = true; + } } @@ -762,16 +775,19 @@ public: */ bool is_refined( unsigned ilevel, int i, int j, int k ) const { - if( ilevel == levelmax() ){ + if( bhave_refmask ) + return !(*m_ref_masks[ilevel-1])(offset(ilevel,0)+i/2,offset(ilevel,1)+j/2,offset(ilevel,2)+k/2); + + /*if( ilevel == levelmax() ){ if( !bhave_refmask ) return false; else if( ref_mask(offset(levelmax(),0)+i/2,offset(levelmax(),1)+j/2,offset(levelmax(),2)+k/2) ) return false; else return true; - } + }*/ if( ilevel == levelmax()-1 && bhave_refmask ) - return ref_mask(i,j,k); + return (*m_ref_masks[ilevel])(i,j,k); if( i < offset(ilevel+1,0) || i >= offset(ilevel+1, 0)+(int)size(ilevel+1,0)/2 || j < offset(ilevel+1,1) || j >= offset(ilevel+1, 1)+(int)size(ilevel+1,1)/2 || @@ -849,6 +865,9 @@ public: } m_levelmin = lmax; + + for( unsigned i=0; i<= lmax; ++i ) + m_ref_masks.push_back( new refinement_mask(size(i,0),size(i,1),size(i,2),(short)(i!=lmax)) ); } //! multiply entire grid hierarchy by a constant @@ -937,7 +956,12 @@ public: GridHierarchy& operator=( const GridHierarchy& gh ) { bhave_refmask = gh.bhave_refmask; - ref_mask = gh.ref_mask; + + if( bhave_refmask ) + { + for( unsigned i=0; i<=gh.levelmax(); ++i ) + m_ref_masks.push_back( new refinement_mask( *(gh.m_ref_masks[i]) ) ); + } if( !is_consistent(gh) ) { @@ -1002,6 +1026,8 @@ public: m_xoffabs.push_back( 2*(m_xoffabs.back() + xoff) ); m_yoffabs.push_back( 2*(m_yoffabs.back() + yoff) ); m_zoffabs.push_back( 2*(m_zoffabs.back() + zoff) ); + + m_ref_masks.push_back( new refinement_mask(nx,ny,nz,0.0) ); } /*! cut a refinement patch to the specified size diff --git a/plugins/convex_hull.hh b/plugins/convex_hull.hh index 7785860..30a5227 100644 --- a/plugins/convex_hull.hh +++ b/plugins/convex_hull.hh @@ -216,18 +216,20 @@ struct convex_hull{ } template< typename T > - bool check_point( const T* x ) const + bool check_point( const T* x, double dist = 0.0 ) const { + dist *= -1.0; + for( size_t i=0; i level_dist_; + void apply_shift( size_t Np, double *p, int *shift, int levelmin ) { double dx = 1.0/(1< 0; --ilevel ) + { + dx = 1.0/(1ul<<(ilevel)); + level_dist_[ilevel-1] = level_dist_[ilevel] + padding_ * dx; + } } ~region_convex_hull_plugin() @@ -114,8 +124,8 @@ public: // it might have enlarged it, but who cares... } - bool query_point( double *x ) - { return phull_->check_point( x ); } + bool query_point( double *x, int ilevel ) + { return phull_->check_point( x, level_dist_[ilevel] ); } bool is_grid_dim_forced( size_t* ndims ) { return false; } diff --git a/plugins/region_ellipsoid.cc b/plugins/region_ellipsoid.cc index ad4a974..33c4ff5 100644 --- a/plugins/region_ellipsoid.cc +++ b/plugins/region_ellipsoid.cc @@ -164,11 +164,13 @@ protected: float *Q; float *u; + float detA, detA13; + float v1[3],v2[3],v3[3],r1,r2,r3; float V[9], mu[3]; - bool axes_computed; + bool axes_computed, hold_point_data; void compute_axes( void ) { @@ -286,7 +288,7 @@ protected: public: min_ellipsoid( size_t N_, double* P ) - : N( N_ ), axes_computed( false ) + : N( N_ ), axes_computed( false ), hold_point_data( true ) { // --- initialize --- LOGINFO("computing minimum bounding ellipsoid from %lld points",N); @@ -339,10 +341,13 @@ public: Inverse_3x3( Ainv, A ); for( size_t i=0; i<9; ++i ){ A[i] /= 3.0; Ainv[i] *= 3.0; } + + detA = Determinant_3x3( A ); + detA13 = pow( detA, 1.0/3.0 ); } min_ellipsoid( const double* A_, const double *c_ ) - : N( 0 ), axes_computed( false ) + : N( 0 ), axes_computed( false ), hold_point_data( false ) { for( int i=0; i<9; ++i ) { A[i] = A_[i]; Ainv[i] = 0.0; } @@ -350,15 +355,48 @@ public: c[i] = c_[i]; } + min_ellipsoid( const min_ellipsoid& e ) + : N( 0 ), hold_point_data( false ) + { + for( int i=0; i<16; ++i ) + X[i] = e.X[i]; + + for( int i=0; i<3; ++i ) + { + c[i] = e.c[i]; + v1[i] = e.v1[i]; + v2[i] = e.v2[i]; + v3[i] = e.v3[i]; + mu[i] = e.mu[i]; + } + + for( int i=0; i<9; ++i ) + { + A[i] = e.A[i]; + Ainv[i] = e.Ainv[i]; + V[i] = e.V[i]; + } + + N = e.N; + detA = e.detA; + detA13 = e.detA13; + axes_computed = e.axes_computed; + } + ~min_ellipsoid() { - delete[] u; - delete[] Q; + if( hold_point_data ) + { + delete[] u; + delete[] Q; + } } template - bool check_point( const T *x ) + bool check_point( const T *x, double dist = 0.0 ) { + dist = (dist + 1.0) * detA13; + float q[3] = {x[0]-c[0],x[1]-c[1],x[2]-c[2]}; double r = 0.0; @@ -416,17 +454,12 @@ public: void expand_ellipsoid( float dr ) { - //print(); - - - - if( !axes_computed ) { - std::cerr << "computing axes.....\n"; + LOGUSER("computing ellipsoid axes....."); compute_axes(); } - + float muold[3] = {mu[0],mu[1],mu[2]}; float munew[3]; for( int i=0; i<3; ++i ) munew[i] = sgn(mu[i])/sqr(1.0/sqrt(fabs(mu[i]))+dr); @@ -446,7 +479,11 @@ public: Inverse_3x3( A, Ainv ); - //print(); + LOGUSER("computing ellipsoid axes....."); + compute_axes(); + + LOGINFO("mu = %f %f %f -> %f %f %f", muold[0], muold[1], muold[2], mu[0], mu[1], mu[2]); + //print(); } }; @@ -458,12 +495,10 @@ public: class region_ellipsoid_plugin : public region_generator_plugin{ private: - min_ellipsoid *pellip_; + std::vector< min_ellipsoid * > pellip_; int shift[3], shift_level, padding_; double vfac_; bool do_extra_padding_; - - void apply_shift( size_t Np, double *p, int *shift, int levelmin ) { @@ -482,6 +517,10 @@ public: { std::vector pp; + for( int i=0; i<=levelmax_; ++i ) + pellip_.push_back( NULL ); + + // sanity check if( !cf.containsKey("setup", "region_point_file") && !( cf.containsKey("setup","region_ellipsoid_matrix[0]") && @@ -533,7 +572,9 @@ public: shift_level = point_levelmin; } - pellip_ = new min_ellipsoid( pp.size()/3, &pp[0] ); + + + pellip_[levelmax_-1] = new min_ellipsoid( pp.size()/3, &pp[0] ); } else { @@ -550,7 +591,7 @@ public: strtmp = cf.getValue("setup","region_ellipsoid_center"); sscanf( strtmp.c_str(), "%lf,%lf,%lf", &c[0],&c[1],&c[2] ); - pellip_ = new min_ellipsoid( A, c ); + pellip_[levelmax_-1] = new min_ellipsoid( A, c ); } @@ -586,8 +627,8 @@ public: // output the center float c[3], A[9]; - pellip_->get_center( c ); - pellip_->get_matrix( A ); + pellip_[levelmax_-1]->get_center( c ); + pellip_[levelmax_-1]->get_matrix( A ); LOGINFO("Region center for ellipsoid determined at\n\t xc = ( %f %f %f )",c[0],c[1],c[2]); LOGINFO("Ellipsoid matrix determined as\n\t ( %f %f %f )\n\t A = ( %f %f %f )\n\t ( %f %f %f )", @@ -598,8 +639,18 @@ public: //expand the ellipsoid by one grid cell unsigned levelmax = cf.getValue("setup","levelmax"); + unsigned npad = cf.getValue("setup","padding"); double dx = 1.0/(1ul<expand_ellipsoid( dx ); + pellip_[levelmax_-1]->expand_ellipsoid( dx ); + + + // generate the higher level ellipsoids + for( int ilevel = levelmax_-1; ilevel > 0; --ilevel ) + { + dx = 1.0/(1ul<<(ilevel)); + pellip_[ ilevel-1 ] = new min_ellipsoid( *pellip_[ ilevel ] ); + pellip_[ ilevel-1 ]->expand_ellipsoid( npad*dx ); + } @@ -617,12 +668,21 @@ public: ~region_ellipsoid_plugin() { - delete pellip_; + for( int i=0; i<=levelmax_; ++i ) + if( pellip_[i] != NULL ) + delete pellip_[i]; } void get_AABB( double *left, double *right, unsigned level ) { - pellip_->get_AABB( left, right ); + if( level <= levelmin_ ) + { + left[0] = left[1] = left[2] = 0.0; + right[0] = right[1] = right[2] = 1.0; + return; + } + + pellip_[level-1]->get_AABB( left, right ); double dx = 1.0/(1ul<check_point( x ); } + bool query_point( double *x, int level ) + { + return pellip_[level]->check_point( x ); + } bool is_grid_dim_forced( size_t* ndims ) { return false; } @@ -654,7 +716,7 @@ public: void get_center( double *xc ) { float c[3]; - pellip_->get_center( c ); + pellip_[levelmax_]->get_center( c ); xc[0] = c[0]; xc[1] = c[1]; @@ -666,7 +728,7 @@ public: double dx = 1.0/(1<get_center( c ); + pellip_[levelmax_]->get_center( c ); xc[0] = c[0]+shift[0]*dx; xc[1] = c[1]+shift[1]*dx; diff --git a/region_generator.cc b/region_generator.cc index 54e7d4d..381e9f7 100644 --- a/region_generator.cc +++ b/region_generator.cc @@ -173,8 +173,10 @@ public: //fprintf(stderr,"left = %f,%f,%f - right = %f,%f,%f\n",left[0],left[1],left[2],right[0],right[1],right[2]); } - bool query_point( double *x ) + bool query_point( double *x, int ilevel ) { + return true; +/* bool check = true; double dx; for( int i=0; i<3; ++i ) @@ -186,6 +188,8 @@ public: check &= ((dx >= padding_fine_) & (dx <= lxref_[i]-padding_fine_)); } return check; + */ + } bool is_grid_dim_forced( size_t* ndims ) diff --git a/region_generator.hh b/region_generator.hh index 8173adb..ebdc6e0 100644 --- a/region_generator.hh +++ b/region_generator.hh @@ -12,10 +12,14 @@ class region_generator_plugin{ public: config_file *pcf_; + int levelmin_, levelmax_; + public: region_generator_plugin( config_file& cf ) : pcf_( &cf ) { + levelmin_ = cf.getValue("setup","levelmin"); + levelmax_ = cf.getValue("setup","levelmax"); } //! destructor @@ -25,7 +29,7 @@ public: virtual void get_AABB( double *left, double *right, unsigned level) = 0; //! query whether a point intersects the region - virtual bool query_point( double *x ) = 0; + virtual bool query_point( double *x, int level ) = 0; //! query whether the region generator explicitly forces the grid dimensions virtual bool is_grid_dim_forced( size_t *ndims ) = 0;