Newer
Older
Tomáš Oberhuber
committed
/***************************************************************************
tnlGrid.h - description
-------------------
begin : Dec 12, 2010
copyright : (C) 2010 by Tomas Oberhuber
email : tomas.oberhuber@fjfi.cvut.cz
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#ifndef TNLGRID_H_
#define TNLGRID_H_
#include <iomanip>
#include <fstream>
#include <core/tnlAssert.h>
#include <core/tnlArray.h>
#include <core/tnlVector.h>
Tomáš Oberhuber
committed
using namespace std;
template< int Dimensions, typename Real = double, typename Device = tnlHost, typename Index = int >
Tomáš Oberhuber
committed
class tnlGrid : public tnlArray< Dimensions, Real, Device, Index >
{
//! We do not allow constructor without parameters.
tnlGrid();
//! We do not allow copy constructor without object name.
tnlGrid( const tnlGrid< Dimensions, Real, Device, Index >& a );
public:
tnlGrid( const tnlString& name );
tnlGrid( const tnlString& name,
const tnlGrid< Dimensions, Real, tnlHost, Index >& grid );
tnlGrid( const tnlString& name,
const tnlGrid< Dimensions, Real, tnlCuda, Index >& grid );
const tnlTuple< Dimensions, Index >& getDimensions() const;
Tomáš Oberhuber
committed
//! Sets the dimensions
/***
* This method also must recompute space steps. It is save to call setDimensions and
* setDomain in any order. Both recompute the space steps.
*/
bool setDimensions( const tnlTuple< Dimensions, Index >& dimensions );
Tomáš Oberhuber
committed
//! Sets the computation domain in form of "rectangle".
/***
* This method also must recompute space steps. It is save to call setDimensions and
* setDomain in any order. Both recompute the space steps.
*/
bool setDomain( const tnlTuple< Dimensions, Real >& lowerCorner,
const tnlTuple< Dimensions, Real >& upperCorner );
Tomáš Oberhuber
committed
//! Set dimensions and domain of the grid using another grid as a template
bool setLike( const tnlGrid< Dimensions, Real, tnlHost, Index >& v );
//! Set dimensions and domain of the grid using another grid as a template
bool setLike( const tnlGrid< Dimensions, Real, tnlCuda, Index >& v );
const tnlTuple< Dimensions, Real >& getDomainLowerCorner() const;
Tomáš Oberhuber
committed
const tnlTuple< Dimensions, Real >& getDomainUpperCorner() const;
Tomáš Oberhuber
committed
const tnlTuple< Dimensions, Real >& getSpaceSteps() const;
Tomáš Oberhuber
committed
tnlString getType() const;
bool operator == ( const tnlGrid< Dimensions, Real, Device, Index >& array ) const;
bool operator != ( const tnlGrid< Dimensions, Real, Device, Index >& array ) const;
template< typename Real2, typename Device2, typename Index2 >
Tomáš Oberhuber
committed
tnlGrid< Dimensions, Real, Device, Index >& operator = ( const tnlGrid< Dimensions, Real2, Device2, Index2 >& array );
//! This method interpolates value at given point.
Real getValue( const tnlTuple< Dimensions, Real >& point ) const;
Tomáš Oberhuber
committed
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
//! Interpolation for 1D grid.
Real getValue( const Real& x ) const;
//! Interpolation for 2D grid.
Real getValue( const Real& x,
const Real& y ) const;
//! Interpolation for 3D grid.
Real getValue( const Real& x,
const Real& y,
const Real& z ) const;
//! Forward difference w.r.t x
Real Partial_x_f( const Index i1 ) const;
//! Forward difference w.r.t x in two dimensions
Real Partial_x_f( const Index i1,
const Index i2 ) const;
//! Forward difference w.r.t x in three dimensions
Real Partial_x_f( const Index i1,
const Index i2,
const Index i3 ) const;
//! Backward difference w.r.t x
Real Partial_x_b( const Index i1 ) const;
//! Backward difference w.r.t x in two dimensions
Real Partial_x_b( const Index i1,
const Index i2 ) const;
//! Backward difference w.r.t x in three dimensions
Real Partial_x_b( const Index i1,
const Index i2,
const Index i3 ) const;
//! Central difference w.r.t. x
Real Partial_x( const Index i1 ) const;
//! Central difference w.r.t. x in two dimensions
Real Partial_x( const Index i1,
const Index i2 ) const;
//! Central difference w.r.t. x
Real Partial_x( const Index i1,
const Index i2,
const Index i3 ) const;
//! Second order difference w.r.t. x
Real Partial_xx( const Index i1 ) const;
//! Second order difference w.r.t. x in two dimensions
Real Partial_xx( const Index i1,
const Index i2 ) const;
//! Second order difference w.r.t. x
Real Partial_xx( const Index i1,
const Index i2,
const Index i3 ) const;
//! Forward difference w.r.t y
Real Partial_y_f( const Index i1,
const Index i2 ) const;
//! Forward difference w.r.t y in three dimensions
Real Partial_y_f( const Index i1,
const Index i2,
const Index i3 ) const;
//! Backward difference w.r.t y
Real Partial_y_b( const Index i1,
const Index i2 ) const;
//! Backward difference w.r.t y in three dimensions
Real Partial_y_b( const Index i1,
const Index i2,
const Index i3 ) const;
//! Central difference w.r.t y
Real Partial_y( const Index i1,
const Index i2 ) const;
//! Central difference w.r.t y
Real Partial_y( const Index i1,
const Index i2,
const Index i3 ) const;
//! Second order difference w.r.t. y
Real Partial_yy( const Index i1,
const Index i2 ) const;
//! Second order difference w.r.t. y in three dimensions
Real Partial_yy( const Index i1,
const Index i2,
const Index i3 ) const;
//! Forward difference w.r.t z
Real Partial_z_f( const Index i1,
const Index i2,
const Index i3 ) const;
//! Backward difference w.r.t z
Real Partial_z_b( const Index i1,
const Index i2,
const Index i3 ) const;
//! Central difference w.r.t z
Real Partial_z( const Index i1,
const Index i2,
const Index i3 ) const;
//! Second order difference w.r.t. z
Real Partial_zz( const Index i1,
const Index i2,
const Index i3 ) const;
//! Set space dependent Dirichlet boundary conditions
void setDirichletBC( const tnlGrid< Dimensions, Real, Device, Index >&bc,
const tnlTuple< Dimensions, bool >& lowerBC,
const tnlTuple< Dimensions, bool >& upperBC );
Tomáš Oberhuber
committed
//! Set constant Dirichlet boundary conditions
void setDirichletBC( const Real& bc,
const tnlTuple< Dimensions, bool >& lowerBC,
const tnlTuple< Dimensions, bool >& upperBC );
Tomáš Oberhuber
committed
//! Set space dependent Neumann boundary conditions
void setNeumannBC( const tnlGrid< Dimensions, Real, Device, Index >&bc,
const tnlTuple< Dimensions, bool >& lowerBC,
const tnlTuple< Dimensions, bool >& upperBC );
Tomáš Oberhuber
committed
//! Set constant Neumann boundary conditions
void setNeumannBC( const Real& bc,
const tnlTuple< Dimensions, bool >& lowerBC,
const tnlTuple< Dimensions, bool >& upperBC );
Tomáš Oberhuber
committed
Real getMax() const;
Real getMin() const;
Real getAbsMax() const;
Real getAbsMin() const;
Real getLpNorm() const;
Real getSum() const;
Real getDifferenceMax( const tnlVector< Real, tnlHost, Index >& v ) const;
Tomáš Oberhuber
committed
Real getDifferenceMin( const tnlVector< Real, tnlHost, Index >& v ) const;
Tomáš Oberhuber
committed
Real getDifferenceAbsMax( const tnlVector< Real, tnlHost, Index >& v ) const;
Tomáš Oberhuber
committed
Real getDifferenceAbsMin( const tnlVector< Real, tnlHost, Index >& v ) const;
Tomáš Oberhuber
committed
Real getDifferenceLpNorm( const tnlVector< Real, tnlHost, Index >& v, const Real& p ) const;
Tomáš Oberhuber
committed
Real getDifferenceSum( const tnlVector< Real, tnlHost, Index >& v ) const;
Tomáš Oberhuber
committed
//! Method for saving the object to a file as a binary data
bool save( tnlFile& file ) const;
//! Method for restoring the object from a file
bool load( tnlFile& file );
bool save( const tnlString& fileName ) const;
bool load( const tnlString& fileName );
//! This method writes the grid in some format suitable for some other preprocessing.
/*! Possible formats are:
* 1) Gnuplot format (gnuplot)
* 2) VTI format (vti)
* 3) Povray format (povray) - only for 3D.
*/
bool draw( const tnlString& fileName,
const tnlString& format,
const tnlTuple< Dimensions, Index > steps = ( tnlTuple< Dimensions, Index > ) 1 ) const;
Tomáš Oberhuber
committed
protected:
tnlTuple< Dimensions, Real > domainLowerCorner, domainUpperCorner, spaceSteps;
Tomáš Oberhuber
committed
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
};
#ifdef HAVE_CUDA
template< int Dimensions, typename Real, typename Index >
__global__ void setConstantDirichletBC( const Index xSize,
const Index ySize,
const Real bc,
Real* u );
template< int Dimensions, typename Real, typename Index >
__global__ void setDirichletBC( const Index xSize,
const Index ySize,
const Real* bc,
Real* u );
template< int Dimensions, typename Real, typename Index >
__global__ void setConstantNeumannBC( const Index xSize,
const Index ySize,
const Real hx,
const Real hy,
const Real bc,
Real* u );
template< int Dimensions, typename Real, typename Index >
__global__ void setNeumannBC( const Index xSize,
const Index ySize,
const Real hx,
const Real hy,
const Real* bc,
Real* u );
#endif
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
tnlGrid< Dimensions, Real, Device, Index > :: tnlGrid( const tnlString& name )
: tnlArray< Dimensions, Real, Device, Index >( name )
{
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
tnlGrid< Dimensions, Real, Device, Index > :: tnlGrid( const tnlString& name,
const tnlGrid< Dimensions, Real, tnlHost, Index >& grid )
: tnlArray< Dimensions, Real, Device, Index >( name, grid )
{
this -> setDomain( grid. getDomainLowerCorner(),
grid. getDomainUpperCorner() );
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
tnlGrid< Dimensions, Real, Device, Index > :: tnlGrid( const tnlString& name,
const tnlGrid< Dimensions, Real, tnlCuda, Index >& grid )
: tnlArray< Dimensions, Real, Device, Index >( name, grid )
{
this -> setDomain( grid. getDomainLowerCorner(),
grid. getDomainUpperCorner() );
}
template< int Dimensions, typename Real, typename Device, typename Index >
const tnlTuple< Dimensions, Index >& tnlGrid< Dimensions, Real, Device, Index > :: getDimensions() const
Tomáš Oberhuber
committed
{
return tnlArray< Dimensions, Real, Device, Index > :: getDimensions();
}
template< int Dimensions, typename Real, typename Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: setDimensions( const tnlTuple< Dimensions, Index >& dimensions )
Tomáš Oberhuber
committed
{
if( ! tnlArray< Dimensions, Real, Device, Index > :: setDimensions( dimensions ) )
return false;
for( int i = 0; i < Dimensions; i ++ )
spaceSteps[ i ] = ( domainUpperCorner[ i ] - domainLowerCorner[ i ] ) / ( Real ) ( this -> getDimensions()[ i ] - 1 );
return true;
}
template< int Dimensions, typename Real, typename Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: setDomain( const tnlTuple< Dimensions, Real >& lowerCorner,
const tnlTuple< Dimensions, Real >& upperCorner )
Tomáš Oberhuber
committed
{
if( lowerCorner >= upperCorner )
{
cerr << "Wrong parameters for the grid domain of " << this -> getName() << ". The low corner must by smaller than the high corner." << endl
<< "lowerCorner = " << lowerCorner << endl << "upperCorner = " << upperCorner << endl;
return false;
}
domainLowerCorner = lowerCorner;
domainUpperCorner = upperCorner;
for( int i = 0; i < Dimensions; i ++ )
spaceSteps[ i ] = ( domainUpperCorner[ i ] - domainLowerCorner[ i ] ) / ( Real ) ( this -> getDimensions()[ i ] - 1 );
return true;
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
bool tnlGrid< Dimensions, Real, Device, Index > :: setLike( const tnlGrid< Dimensions, Real, tnlHost, Index >& v )
{
return tnlArray< Dimensions, Real, Device, Index > :: setLike( v ) &&
this -> setDomain( v. getDomainLowerCorner(), v. getDomainUpperCorner() );
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
bool tnlGrid< Dimensions, Real, Device, Index > :: setLike( const tnlGrid< Dimensions, Real, tnlCuda, Index >& v )
{
return tnlArray< Dimensions, Real, Device, Index > :: setLike( v ) &&
this -> setDomain( v. getDomainLowerCorner(), v. getDomainUpperCorner() );
}
template< int Dimensions, typename Real, typename Device, typename Index >
const tnlTuple< Dimensions, Real >& tnlGrid< Dimensions, Real, Device, Index > :: getDomainLowerCorner() const
Tomáš Oberhuber
committed
{
return this -> domainLowerCorner;
}
template< int Dimensions, typename Real, typename Device, typename Index >
const tnlTuple< Dimensions, Real >& tnlGrid< Dimensions, Real, Device, Index > :: getDomainUpperCorner() const
Tomáš Oberhuber
committed
{
return this -> domainUpperCorner;
}
template< int Dimensions, typename Real, typename Device, typename Index >
const tnlTuple< Dimensions, Real >& tnlGrid< Dimensions, Real, Device, Index > :: getSpaceSteps() const
Tomáš Oberhuber
committed
{
return spaceSteps;
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
tnlString tnlGrid< Dimensions, Real, Device, Index > :: getType() const
{
return tnlString( "tnlGrid< ") +
tnlString( Dimensions ) +
tnlString( ", " ) +
tnlString( getParameterType< Real >() ) +
Tomáš Oberhuber
committed
tnlString( ", " ) +
Tomáš Oberhuber
committed
tnlString( ", " ) +
tnlString( getParameterType< Index >() ) +
Tomáš Oberhuber
committed
tnlString( " >" );
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
bool tnlGrid< Dimensions, Real, Device, Index > :: operator == ( const tnlGrid< Dimensions, Real, Device, Index >& grid ) const
{
tnlAssert( this -> getDomainLowerCorner() == grid. getDomainLowerCorner() &&
this -> getDomainUpperCorner() == grid. getDomainUpperCorner(),
cerr << "You are attempting to compare two grids with different domains." << endl
<< "First grid name is " << this -> getName()
<< " domain is ( " << this -> getDomainLowerCorner() << " )- ("
<< this -> getDomainUpperCorner() << ")" << endl
<< "Second grid is " << grid. getName()
<< " domain is ( " << grid. getDomainLowerCorner() << " )- ("
<< grid. getDomainUpperCorner() << ")" << endl; );
return tnlArray< Dimensions, Real, Device, Index > :: operator == ( grid );
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
bool tnlGrid< Dimensions, Real, Device, Index > :: operator != ( const tnlGrid< Dimensions, Real, Device, Index >& grid ) const
{
return ! ( (* this ) == grid );
}
template< int Dimensions, typename Real, typename Device, typename Index >
template< typename Real2, typename Device2, typename Index2 >
Tomáš Oberhuber
committed
tnlGrid< Dimensions, Real, Device, Index >&
tnlGrid< Dimensions, Real, Device, Index >
:: operator = ( const tnlGrid< Dimensions, Real2, Device2, Index2 >& grid )
{
tnlAssert( this -> getDimensions() == grid. getDimensions(),
cerr << "You are attempting to assign two arrays with different dimensions." << endl
<< "First array name is " << this -> getName()
<< " dimensions are ( " << this -> getDimensions() << " )" << endl
<< "Second array is " << grid. getName()
<< " dimensions are ( " << grid. getDimensions() << " )" << endl; );
tnlVector< Real, Device, Index > :: operator = ( grid );
Tomáš Oberhuber
committed
return ( *this );
}
template< int Dimensions, typename Real, typename Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getValue( const tnlTuple< Dimensions, Real >& point ) const
Tomáš Oberhuber
committed
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
{
if( Dimensions == 1)
{
tnlAssert( 0, cerr << "Interpolation is not implemented yet.");
}
if( Dimensions == 2 )
{
Real x = ( point[ 0 ] - domainLowerCorner[ 0 ] ) / spaceSteps[ 0 ];
Real y = ( point[ 1 ] - domainLowerCorner[ 1 ] ) / spaceSteps[ 1 ];
Index ix = ( Index ) ( x );
Index iy = ( Index ) ( y );
Real dx = x - ( Real ) ix;
Real dy = y - ( Real ) iy;
if( iy >= this -> getDimensions()[ tnlY ] - 1 )
{
if( ix >= this -> getDimensions()[ tnlX ] - 1 )
return this -> getElement( this -> getDimensions()[ tnlX ] - 1,
this -> getDimensions()[ tnlY ] - 1 );
return ( Real( 1.0 ) - dx ) * this -> getElement( ix,
this -> getDimensions()[ tnlY ] - 1 ) +
dx * this -> getElement( ix + 1,
this -> getDimensions()[ tnlY ] - 1 );
}
if( ix >= this -> getDimensions()[ tnlX ] - 1 )
return ( Real( 1.0 ) - dy ) * this -> getElement( this -> getDimensions()[ tnlX ] - 1,
iy ) +
dy * this -> getElement( this -> getDimensions()[ tnlX ] - 1,
iy + 1 );
Real a1, a2;
a1 = ( Real( 1.0 ) - dx ) * this -> getElement( ix, iy ) +
dx * this -> getElement( ix + 1, iy );
a2 = ( Real( 1.0 ) - dx ) * this -> getElement( ix, iy + 1 ) +
dx * this -> getElement( ix + 1, iy + 1 );
return ( Real( 1.0 ) - dy ) * a1 + dy * a2;
}
if( Dimensions == 3 )
{
tnlAssert( 0, cerr << "Interpolation is not implemented yet.");
}
if( Dimensions > 3 )
{
cerr << "Interpolation for 4D grids is not implemented yet." << endl;
abort();
}
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: getValue( const Real& x ) const
{
return this -> getValue( tnlTuple< 1, Real >( x ) );
Tomáš Oberhuber
committed
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: getValue( const Real& x,
const Real& y ) const
{
return this -> getValue( tnlTuple< 2, Real >( x, y ) );
Tomáš Oberhuber
committed
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: getValue( const Real& x,
const Real& y,
const Real& z ) const
{
return this -> getValue( tnlTuple< 3, Real >( x, y, z ) );
Tomáš Oberhuber
committed
}
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_f( const Index i1 ) const
{
tnlAssert( Dimensions == 1,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 1 is expected." << endl; );
tnlAssert( i1 >= 0 &&
i1 < ( this -> getDimensions(). x() - 1 ),
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions(). y() - 1 ) << " )" << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 + 1 ) -
tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 ) ) / Hx;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_f( const Index i1,
const Index i2 ) const
{
tnlAssert( Dimensions == 2,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 2 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 >= 0 &&
i1 < this -> getDimensions(). x() - 1 && i2 <= this -> getDimensions(). y() - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions() .x() - 1 ) << " ) " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions(). y() - 1 ) << " ] " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( this -> getElement( i1 + 1, i2 ) -
this -> getElement( i1, i2 ) ) / Hx;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_f( const Index i1,
const Index i2,
const Index i3 ) const
{
tnlAssert( Dimensions == 3,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 3 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 >= 0 && i3 >= 0 &&
i1 < this -> getDimensions()[ tnlX ] - 1 &&
i2 <= this -> getDimensions()[ tnlY ] - 1 &&
i3 <= this -> getDimensions()[ tnlZ ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ) " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl
<< " i3 = " << i3 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlZ ] - 1 ) << " ] " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( this -> getElement( i1 + 1, i2, i3 ) -
this -> getElement( i1, i2, i3 ) ) / Hx;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_b( const Index i1 ) const
{
tnlAssert( Dimensions == 1,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 1 is expected." << endl; );
tnlAssert( i1 > 0 &&
i1 <= ( tnlArray< 1, Real, tnlHost, Index > :: getDimensions()[ tnlX ] - 1 ),
cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
( tnlArray< 1, Real, tnlHost, Index > :: getDimensions()[ tnlX ] - 1 ) << " ] " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 ) -
tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 - 1 ) ) / Hx;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_b( const Index i1,
const Index i2 ) const
{
tnlAssert( Dimensions == 2,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 2 are expected." << endl; );
tnlAssert( i1 > 0 && i2 >= 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 && i2 <= this -> getDimensions()[ tnlY ] - 1,
cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( this -> getElement( i1, i2 ) -
this -> getElement( i1 - 1, i2 ) ) / Hx;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_b( const Index i1,
const Index i2,
const Index i3 ) const
{
tnlAssert( Dimensions == 3,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 3 are expected." << endl; );
tnlAssert( i1 > 0 && i2 >= 0 && i3 >= 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 &&
i2 <= this -> getDimensions()[ tnlY ] - 1 &&
i3 <= this -> getDimensions()[ tnlZ ] - 1,
cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl
<< " i3 = " << i3 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlZ ] - 1 ) << " ] " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( this -> getElement( i1, i2, i3 ) -
this -> getElement( i1 - 1, i2, i3 ) ) / Hx;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x( const Index i1 ) const
{
tnlAssert( Dimensions == 1,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 1 is expected." << endl; );
tnlAssert( i1 > 0 &&
i1 < ( this -> getDimensions()[ tnlX ] - 1 ),
cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ) " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 + 1 ) -
tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 - 1 ) ) / ( 2.0 * Hx );
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x( const Index i1,
const Index i2 ) const
{
tnlAssert( Dimensions == 2,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 2 are expected." << endl; );
tnlAssert( i1 > 0 && i2 >= 0 &&
i1 < this -> getDimensions()[ tnlX ] - 1 && i2 <= this -> getDimensions()[ tnlY ] - 1,
cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ) " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( this -> getElement( i1 + 1, i2 ) -
this -> getElement( i1 - 1, i2 ) ) / ( 2.0 * Hx );
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x( const Index i1,
const Index i2,
const Index i3 ) const
{
tnlAssert( Dimensions == 3,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 3 are expected." << endl; );
tnlAssert( i1 > 0 && i2 >= 0 && i3 >= 0 &&
i1 < this -> getDimensions()[ tnlX ] - 1 &&
i2 <= this -> getDimensions()[ tnlY ] - 1 &&
i3 <= this -> getDimensions()[ tnlZ ] - 1,
cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ) " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl
<< " i3 = " << i3 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlZ ] - 1 ) << " ] " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( this -> getElement( i1 + 1, i2, i3 ) -
this -> getElement( i1 - 1, i2, i3 ) ) / ( 2.0 * Hx );
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_xx( const Index i1 ) const
{
tnlAssert( Dimensions == 1,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 1 is expected." << endl; );
tnlAssert( i1 > 0 &&
i1 < ( tnlArray< 1, Real, tnlHost, Index > :: getDimensions()[ tnlX ] - 1 ),
cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
( tnlArray< 1, Real, tnlHost, Index > :: getDimensions()[ tnlX ] - 1 ) << " ) " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 + 1 ) -
2.0 * tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 ) +
tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 - 1 ) ) / ( Hx * Hx );
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_xx( const Index i1,
const Index i2 ) const
{
tnlAssert( Dimensions == 2,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 2 are expected." << endl; );
tnlAssert( i1 > 0 && i2 >= 0 &&
i1 < this -> getDimensions()[ tnlX ] - 1 && i2 <= this -> getDimensions()[ tnlY ] - 1,
cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ) " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( this -> getElement( i1 + 1, i2 ) -
2.0 * this -> getElement( i1, i2 ) +
this -> getElement( i1 - 1, i2 ) ) / ( Hx * Hx );
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_xx( const Index i1,
const Index i2,
const Index i3 ) const
{
tnlAssert( Dimensions == 3,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 3 are expected." << endl; );
tnlAssert( i1 > 0 && i2 >= 0 && i3 >= 0 &&
i1 < this -> getDimensions()[ tnlX ] - 1 &&
i2 <= this -> getDimensions()[ tnlY ] - 1 &&
i3 <= this -> getDimensions()[ tnlZ ] - 1,
cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ) " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl
<< " i3 = " << i3 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlZ ] - 1 ) << " ] " << endl; );
const Real& Hx = spaceSteps[ 0 ];
tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
return ( this -> getElement( i1 + 1, i2, i3 ) -
2.0 * this -> getElement( i1, i2, i3 ) +
this -> getElement( i1 - 1, i2, i3 ) ) / ( Hx * Hx );
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y_f( const Index i1,
const Index i2 ) const
{
tnlAssert( Dimensions == 2,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 2 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 >= 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 && i2 < this -> getDimensions()[ tnlY ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ) " << endl; );
const Real& Hy = spaceSteps[ 1 ];
tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
return ( this -> getElement( i1, i2 + 1 ) -
this -> getElement( i1, i2 ) ) / Hy;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y_f( const Index i1,
const Index i2,
const Index i3 ) const
{
tnlAssert( Dimensions == 3,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 3 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 >= 0 && i3 >= 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 &&
i2 < this -> getDimensions()[ tnlY ] - 1 &&
i3 <= this -> getDimensions()[ tnlZ ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ) " << endl
<< " i3 = " << i3 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlZ ] - 1 ) << " ] " << endl; );
const Real& Hy = spaceSteps[ 1 ];
tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
return ( this -> getElement( i1, i2 + 1, i3 ) -
this -> getElement( i1, i2, i3 ) ) / Hy;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y_b( const Index i1,
const Index i2 ) const
{
tnlAssert( Dimensions == 2,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 2 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 > 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 && i2 <= this -> getDimensions()[ tnlY ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl; );
const Real& Hy = spaceSteps[ 1 ];
tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
return ( this -> getElement( i1, i2 ) -
this -> getElement( i1, i2 - 1 ) ) / Hy;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y_b( const Index i1,
const Index i2,
const Index i3 ) const
{
tnlAssert( Dimensions == 3,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 3 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 > 0 && i3 >= 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 &&
i2 <= this -> getDimensions()[ tnlY ] - 1 &&
i3 <= this -> getDimensions()[ tnlZ ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl
<< " i3 = " << i3 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlZ ] - 1 ) << " ] " << endl; );
const Real& Hy = spaceSteps[ 1 ];
tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
return ( this -> getElement( i1, i2, i3 ) -
this -> getElement( i1, i2 - 1, i3 ) ) / Hy;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y( const Index i1,
const Index i2 ) const
{
tnlAssert( Dimensions == 2,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 2 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 > 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 && i2 < this -> getDimensions()[ tnlY ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ) " << endl; );
const Real& Hy = spaceSteps[ 1 ];
tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
return ( this -> getElement( i1, i2 + 1 ) -
this -> getElement( i1, i2 - 1 ) ) / ( 2.0 * Hy );
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y( const Index i1,
const Index i2,
const Index i3 ) const
{
tnlAssert( Dimensions == 3,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 3 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 > 0 && i3 >= 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 &&
i2 < this -> getDimensions()[ tnlY ] - 1 &&
i3 <= this -> getDimensions()[ tnlZ ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ) " << endl
<< " i3 = " << i3 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlZ ] - 1 ) << " ] " << endl; );
const Real& Hy = spaceSteps[ 1 ];
tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
return ( this -> getElement( i1, i2 + 1, i3 ) -
this -> getElement( i1, i2 - 1, i3 ) ) / ( 2.0 * Hy );
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_yy( const Index i1,
const Index i2 ) const
{
tnlAssert( Dimensions == 2,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 2 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 > 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 && i2 < this -> getDimensions()[ tnlY ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ) " << endl; );
const Real& Hy = spaceSteps[ 1 ];
tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
return ( this -> getElement( i1, i2 + 1 ) -
2.0 * this -> getElement( i1, i2 ) +
this -> getElement( i1, i2 - 1 ) ) / ( Hy * Hy );
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_yy( const Index i1,
const Index i2,
const Index i3 ) const
{
tnlAssert( Dimensions == 3,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 3 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 > 0 && i3 >= 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 &&
i2 < this -> getDimensions()[ tnlY ] - 1 &&
i3 <= this -> getDimensions()[ tnlZ ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ) " << endl
<< " i3 = " << i3 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlZ ] - 1 ) << " ] " << endl; );
const Real& Hy = spaceSteps[ 1 ];
tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
return ( this -> getElement( i1, i2 + 1, i3 ) -
2.0 * this -> getElement( i1, i2, i3 ) +
this -> getElement( i1, i2 - 1, i3 ) ) / ( Hy * Hy );
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_z_f( const Index i1,
const Index i2,
const Index i3 ) const
{
tnlAssert( Dimensions == 3,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 3 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 >= 0 && i3 >= 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 &&
i2 <= this -> getDimensions()[ tnlY ] - 1 &&
i3 < this -> getDimensions()[ tnlZ ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl
<< " i3 = " << i3 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlZ ] - 1 ) << " ) " << endl; );
const Real& Hz = spaceSteps[ 2 ];
tnlAssert( Hz > 0, cerr << "Hz = " << Hz << endl; );
return ( this -> getElement( i1, i2, i3 + 1 ) -
this -> getElement( i1, i2, i3 ) ) / Hz;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_z_b( const Index i1,
const Index i2,
const Index i3 ) const
{
tnlAssert( Dimensions == 3,
cerr << "The array " << this -> getName()
<< " has " << Dimensions << " but 3 are expected." << endl; );
tnlAssert( i1 >= 0 && i2 >= 0 && i3 > 0 &&
i1 <= this -> getDimensions()[ tnlX ] - 1 &&
i2 <= this -> getDimensions()[ tnlY ] - 1 &&
i3 <= this -> getDimensions()[ tnlZ ] - 1,
cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlX ] - 1 ) << " ] " << endl
<< " i2 = " << i2 << " and it should be in [ 0, " <<
( this -> getDimensions()[ tnlY ] - 1 ) << " ] " << endl
<< " i3 = " << i3 << " and it should be in ( 0, " <<
( this -> getDimensions()[ tnlZ ] - 1 ) << " ] " << endl; );
const Real& Hz = spaceSteps[ 2 ];
tnlAssert( Hz > 0, cerr << "Hz = " << Hz << endl; );
return ( this -> getElement( i1, i2, i3 ) -
this -> getElement( i1, i2, i3 - 1 ) ) / Hz;
};
template< int Dimensions, typename Real, typename Device, typename Index >
Tomáš Oberhuber
committed
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_z( const Index i1,
const Index i2,
const Index i3 ) const