#include #include "scalarfield2d.hpp" #include "ibsimu.hpp" ScalarFieldOnPlane::ScalarFieldOnPlane() : _F(NULL) { } ScalarFieldOnPlane::ScalarFieldOnPlane( const ScalarFieldOnPlane &p ) : _size(p._size), _origo(p._origo), _max(p._max), _h(p._h), _div_h(p._div_h) { _vb[0] = p._vb[0]; _vb[1] = p._vb[1]; _vb[2] = p._vb[2]; _psize[0] = p._psize[0]; _psize[1] = p._psize[1]; _porigo[0] = p._porigo[0]; _porigo[1] = p._porigo[1]; _F = new double[_psize[0]*_psize[1]]; memcpy( _F, p._F, _psize[0]*_psize[1]*sizeof(double) ); } ScalarFieldOnPlane &ScalarFieldOnPlane::operator=( const ScalarFieldOnPlane &p ) { _size = p._size; _origo = p._origo; _max = p._max; _h = p._h; _div_h = p._div_h; _vb[0] = p._vb[0]; _vb[1] = p._vb[1]; _vb[2] = p._vb[2]; _psize[0] = p._psize[0]; _psize[1] = p._psize[1]; _porigo[0] = p._porigo[0]; _porigo[1] = p._porigo[1]; if( _F ) delete [] _F; _F = new double[_psize[0]*_psize[1]]; memcpy( _F, p._F, _psize[0]*_psize[1]*sizeof(double) ); return( *this ); } // Construct as a copy of a partial plane from MeshScalarField defined // by vb, level and index range, where range is (imin, jmin, imax, jmax) void ScalarFieldOnPlane::copy_partial_plane( const MeshScalarField &f, uint32_t vb[3], uint32_t level, uint32_t range[4] ) { ibsimu.message(1) << "copy_partial_plane():\n"; ibsimu.message(1) << "vb = {" << vb[0] << ", " << vb[1] << ", " << vb[2] << "}\n"; ibsimu.message(1) << "level = " << level << "\n"; ibsimu.message(1) << "range = {" << range[0] << ", " << range[1] << ", " << range[2] << ", " << range[3] << "}\n"; // Check validity of vb if( vb[0] == vb[1] || vb[0] == vb[2] || vb[1] == vb[2] ) throw( Error( ERROR_LOCATION, "illegal index base vectors, dual use of index" ) ); for( int a = 0; a < 3; a++ ) { if( vb[a] > 2 ) throw( Error( ERROR_LOCATION, "illegal index base vectors, larger than 2" ) ); } // Check level if( level >= f.size(vb[2]) ) throw( Error( ERROR_LOCATION, "level request non-existent in field" ) ); // Check range for( int a = 0; a < 2; a++ ) { if( range[2*a+0] >= range[2*a+1] ) throw( Error( ERROR_LOCATION, "illegal range" ) ); if( range[2*a+1] >= f.size(vb[a]) ) throw( Error( ERROR_LOCATION, "range request non-existent in field" ) ); } // Delete old field if( _F ) delete [] _F; // Set parameters _h = f.h(); _div_h = 1.0/_h; for( int a = 0; a < 3; a++ ) _vb[a] = vb[a]; for( int a = 0; a < 2; a++ ) { _psize[a] = _size(_vb[a]) = range[2*a+1]-range[2*a+0]+1; _porigo[a] = _origo(_vb[a]) = f.origo(_vb[a])+range[2*a+0]*f.h(); } _size(_vb[2]) = 1; _origo(_vb[2]) = f.origo(_vb[2])+level*f.h(); _max = Vec3D( _origo(0)+_h*(_size[0]-1), _origo(1)+_h*(_size[1]-1), _origo(2)+_h*(_size[2]-1) ); ibsimu.message(1) << "origo = {" << _origo[0] << ", " << _origo[1] << ", " << _origo[2] << "}\n"; ibsimu.message(1) << "max = {" << _max[0] << ", " << _max[1] << ", " << _max[2] << "}\n"; // Allocate field _F = new double[_psize[0]*_psize[1]]; // Copy the field uint32_t input[3]; input[_vb[2]] = level; for( uint32_t j = 0; j < _psize[1]; j++ ) { input[_vb[1]] = j + range[2]; for( uint32_t i = 0; i < _psize[0]; i++ ) { input[_vb[0]] = i + range[0]; _F[i+j*_psize[0]] = f( input[0], input[1], input[2] ); } } } ScalarFieldOnPlane::ScalarFieldOnPlane( const MeshScalarField &f, uint32_t vb[3], uint32_t level, uint32_t range[4] ) : _F(NULL) { copy_partial_plane( f, vb, level, range ); } ScalarFieldOnPlane::~ScalarFieldOnPlane() { if( _F ) delete [] _F; } // Interpolating evaluator. Returns field value on a point on the // plane, where x projects to along a plane normal double ScalarFieldOnPlane::operator()( const Vec3D &x ) const { double y[2] = { x[_vb[0]], x[_vb[1]] }; int32_t i = (int32_t)floor( (y[0]-_porigo[0])*_div_h ); int32_t j = (int32_t)floor( (y[1]-_porigo[1])*_div_h ); if( i < 0 ) i = 0; else if( i >= (int32_t)_psize[0]-1 ) i = _psize[0]-2; if( j < 0 ) j = 0; else if( j >= (int32_t)_psize[1]-1 ) j = _psize[1]-2; double t = _div_h*( y[0]-(i*_h+_porigo[0]) ); double u = _div_h*( y[1]-(j*_h+_porigo[1]) ); int32_t ptr = _psize[0]*j + i; return( (1.0-t)*(1.0-u)*_F[ptr] + ( t)*(1.0-u)*_F[ptr+1] + (1.0-t)*( u)*_F[ptr+_psize[0]] + ( t)*( u)*_F[ptr+1+_psize[0]] ); } void ScalarFieldOnBox::copy_plane( uint32_t index, uint32_t vb0, uint32_t vb1, uint32_t vb2, uint32_t levelindex, const MeshScalarField &f ) { uint32_t vb[3] = { vb0, vb1, vb2 }; uint32_t level = _range[2*vb2+levelindex]; uint32_t planerange[4] = { _range[2*vb0+0], _range[2*vb0+1], _range[2*vb1+0], _range[2*vb1+1] }; _plane[index].copy_partial_plane( f, vb, level, planerange ); } // Box range (xmin, xmax, ymin, ymax, zmin, zmax) rounded to closest grid level ScalarFieldOnBox::ScalarFieldOnBox( const MeshScalarField &f, double range[6] ) { // Calculate range indices for( int a = 0; a < 3; a++ ) { _range[2*a+0] = (uint32_t)floor( (range[2*a+0]-f.origo(a))*f.div_h() + 0.5 ); _range[2*a+1] = (uint32_t)floor( (range[2*a+1]-f.origo(a))*f.div_h() + 0.5 ); } ibsimu.message(1) << "Copying field on a boundary of box from (" << _range[0] << ", " << _range[2] << ", " << _range[4] << ") to (" << _range[1] << ", " << _range[3] << ", " << _range[5] << ")\n"; // Check ranges for( int a = 0; a < 3; a++ ) { if( _range[2*a+0] >= _range[2*a+1] ) throw( Error( ERROR_LOCATION, "illegal range" ) ); if( _range[2*a+0] >= f.size(a) ) throw( Error( ERROR_LOCATION, "range request non-existent in field" ) ); if( _range[2*a+1] >= f.size(a) ) throw( Error( ERROR_LOCATION, "range request non-existent in field" ) ); } // xmin, xmax copy_plane( 0, 1, 2, 0, 0, f ); copy_plane( 1, 1, 2, 0, 1, f ); // ymin, ymax copy_plane( 2, 0, 2, 1, 0, f ); copy_plane( 3, 0, 2, 1, 1, f ); // zmin, zmax copy_plane( 4, 0, 1, 2, 0, f ); copy_plane( 5, 0, 1, 2, 1, f ); } ScalarFieldOnBox::~ScalarFieldOnBox() { }