Interpolation, Projection, Transformation

Collaboration diagram for Interpolation, Projection, Transformation:

Topics
	Fast interpolation

Functions
template<class Topology >
Topology::host_vector	dg::forward_transform (const typename Topology::host_vector &in, const Topology &g)
	Transform a vector from dg::xspace (nodal values) to dg::lspace (modal values)
template<class RecursiveHostVector , class real_type , size_t Nd>
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const RecursiveHostVector &x, const aRealTopology< real_type, Nd > &g, std::array< dg::bc, Nd > bcx, std::string method="dg")
	Create interpolation matrix of a list of points in given grid.
template<class host_vector , class real_type , typename = std::enable_if_t<dg::is_vector_v<host_vector>>>
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const host_vector &x, const RealGrid1d< real_type > &g, dg::bc bcx=dg::NEU, std::string method="dg")
	Create interpolation matrix.
template<class host_vector , class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const host_vector &x, const host_vector &y, const aRealTopology2d< real_type > &g, dg::bc bcx=dg::NEU, dg::bc bcy=dg::NEU, std::string method="dg")
	Create interpolation matrix.
template<class host_vector , class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const host_vector &x, const host_vector &y, const host_vector &z, const aRealTopology3d< real_type > &g, dg::bc bcx=dg::NEU, dg::bc bcy=dg::NEU, dg::bc bcz=dg::PER, std::string method="dg")
	Create interpolation matrix.
template<class real_type , size_t Nd>
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const aRealTopology< real_type, Nd > &g_new, const aRealTopology< real_type, Nd > &g_old, std::string method="dg")
	Create interpolation between two grids.
template<class host_vector , class real_type >
real_type	dg::interpolate (dg::space sp, const host_vector &v, real_type x, const RealGrid1d< real_type > &g, dg::bc bcx=dg::NEU)
	Interpolate a vector on a single point on a 1d Grid.
template<class host_vector , class real_type >
real_type	dg::interpolate (dg::space sp, const host_vector &v, real_type x, real_type y, const aRealTopology< real_type, 2 > &g, dg::bc bcx=dg::NEU, dg::bc bcy=dg::NEU)
	Interpolate a vector on a single point on a 2d Grid.
template<class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const thrust::host_vector< real_type > &x, const RealGridX1d< real_type > &g)
	Create interpolation matrix.
template<class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const thrust::host_vector< real_type > &x, const thrust::host_vector< real_type > &y, const aRealTopologyX2d< real_type > &g, dg::bc globalbcz=dg::NEU)
	Create interpolation matrix.
template<class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const thrust::host_vector< real_type > &x, const thrust::host_vector< real_type > &y, const thrust::host_vector< real_type > &z, const aRealTopologyX3d< real_type > &g, dg::bc globalbcz=dg::NEU)
	Create interpolation matrix.
template<class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const RealGridX1d< real_type > &g_new, const RealGridX1d< real_type > &g_old)
	Create interpolation between two grids.
template<class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const aRealTopologyX2d< real_type > &g_new, const aRealTopologyX2d< real_type > &g_old)
	Create interpolation between two grids.
template<class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::interpolation (const aRealTopologyX3d< real_type > &g_new, const aRealTopologyX3d< real_type > &g_old)
	Create interpolation between two grids.
template<class real_type >
real_type	dg::interpolate (dg::space sp, const thrust::host_vector< real_type > &v, real_type x, real_type y, const aRealTopologyX2d< real_type > &g, dg::bc bcx=dg::NEU, dg::bc bcy=dg::NEU)
	Interpolate a vector on a single point on a 2d Grid.
template<class MPITopology , typename = std::enable_if_t<dg::is_vector_v< typename MPITopology::host_vector, MPIVectorTag>>>
dg::MIHMatrix_t< typename MPITopology::value_type >	dg::create::interpolation (const MPITopology &g_new, const MPITopology &g_old, std::string method="dg")
	Create interpolation between two grids.
template<class MPITopology , typename = std::enable_if_t<dg::is_vector_v< typename MPITopology::host_vector, MPIVectorTag>>>
dg::MIHMatrix_t< typename MPITopology::value_type >	dg::create::projection (const MPITopology &g_new, const MPITopology &g_old, std::string method="dg")
	Create a projection between two grids.
template<class RecursiveHostVector , class real_type , size_t Nd>
dg::MIHMatrix_t< real_type >	dg::create::interpolation (const RecursiveHostVector &x, const aRealMPITopology< real_type, Nd > &g, std::array< dg::bc, Nd > bcx, std::string method="dg")
	Create interpolation matrix of a list of points in given grid.
template<class host_vector , class real_type >
dg::MIHMatrix_t< real_type >	dg::create::interpolation (const host_vector &x, const RealMPIGrid1d< real_type > &g, dg::bc bcx=dg::NEU, std::string method="dg")
	Create an MPI row distributed interpolation 1d matrix.
template<class host_vector , class real_type >
dg::MIHMatrix_t< real_type >	dg::create::interpolation (const host_vector &x, const host_vector &y, const aRealMPITopology2d< real_type > &g, dg::bc bcx=dg::NEU, dg::bc bcy=dg::NEU, std::string method="dg")
	Create an MPI row distributed interpolation 2d matrix.
template<class host_vector , class real_type >
dg::MIHMatrix_t< real_type >	dg::create::interpolation (const host_vector &x, const host_vector &y, const host_vector &z, const aRealMPITopology3d< real_type > &g, dg::bc bcx=dg::NEU, dg::bc bcy=dg::NEU, dg::bc bcz=dg::PER, std::string method="dg")
	Create an MPI row distributed interpolation 3d matrix.
template<class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::diagonal (const thrust::host_vector< real_type > &diagonal)
	Create a diagonal matrix.
template<class real_type , size_t Nd>
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::projection (const aRealTopology< real_type, Nd > &g_new, const aRealTopology< real_type, Nd > &g_old, std::string method="dg")
	Create a projection between two grids.
template<class real_type , size_t Nd>
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::transformation (const aRealTopology< real_type, Nd > &g_new, const aRealTopology< real_type, Nd > &g_old)
	Create a transformation matrix between two grids.
template<class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::projection (const RealGridX1d< real_type > &g_new, const RealGridX1d< real_type > &g_old, std::string method="dg")
template<class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::projection (const aRealTopologyX2d< real_type > &g_new, const aRealTopologyX2d< real_type > &g_old, std::string method="dg")
template<class real_type >
dg::SparseMatrix< int, real_type, thrust::host_vector >	dg::create::projection (const aRealTopologyX3d< real_type > &g_new, const aRealTopologyX3d< real_type > &g_old, std::string method="dg")

Detailed Description

Function Documentation

diagonal()

template<class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::diagonal

(

const thrust::host_vector< real_type > &

diagonal

)

Create a diagonal matrix.

This matrix is given by \( D_{ij} = d_i \delta_{ij}\)

Parameters

diagonal

The diagonal elements d_i

Returns

diagonal matrix

forward_transform()

template<class Topology >

Topology::host_vector dg::forward_transform	(	const typename Topology::host_vector &	in,
		const Topology &	g )

Transform a vector from dg::xspace (nodal values) to dg::lspace (modal values)

This can speedup the dg::interpolate function

Parameters

in	input
g	A grid

Returns

the vector in dg::lspace

See also

dg::interpolate

interpolate() [1/3]

template<class host_vector , class real_type >

real_type dg::interpolate	(	dg::space	sp,
		const host_vector &	v,
		real_type	x,
		const RealGrid1d< real_type > &	g,
		dg::bc	bcx = dg::NEU )

Interpolate a vector on a single point on a 1d Grid.

Parameters

sp	Indicate whether the elements of the vector v are in xspace (nodal values) or lspace (modal values) (choose dg::xspace if you don't know what is going on here, It is faster to interpolate in dg::lspace so consider transforming v using dg::forward_transform( ) if you do it very many times)
v	The vector to interpolate
x	X-coordinate of interpolation point
g	The Grid on which to operate

Parameters

bcx

determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.

Returns

interpolated point

interpolate() [2/3]

template<class host_vector , class real_type >

real_type dg::interpolate	(	dg::space	sp,
		const host_vector &	v,
		real_type	x,
		real_type	y,
		const aRealTopology< real_type, 2 > &	g,
		dg::bc	bcx = dg::NEU,
		dg::bc	bcy = dg::NEU )

Interpolate a vector on a single point on a 2d Grid.

Parameters

sp	Indicate whether the elements of the vector v are in xspace (nodal values) or lspace (modal values) (choose dg::xspace if you don't know what is going on here, It is faster to interpolate in dg::lspace so consider transforming v using dg::forward_transform( ) if you do it very many times)
v	The vector to interpolate in dg::xspace, or dg::lspace s.a. dg::forward_transform( )
x	X-coordinate of interpolation point
y	Y-coordinate of interpolation point
g	The Grid on which to operate

Parameters

bcx

determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.

bcy

analogous to bcx, applies to y direction

Returns

interpolated point

interpolate() [3/3]

template<class real_type >

real_type dg::interpolate	(	dg::space	sp,
		const thrust::host_vector< real_type > &	v,
		real_type	x,
		real_type	y,
		const aRealTopologyX2d< real_type > &	g,
		dg::bc	bcx = dg::NEU,
		dg::bc	bcy = dg::NEU )

Interpolate a vector on a single point on a 2d Grid.

Parameters

sp	Indicate whether the elements of the vector v are in xspace or lspace (choose dg::xspace if you don't know what is going on here, It is faster to interpolate in dg::lspace so consider transforming v using dg::forward_transform( ) if you do it very many times)
v	The vector to interpolate in dg::xspace, or dg::lspace s.a. dg::forward_transform( )
x	X-coordinate of interpolation point
y	Y-coordinate of interpolation point
g	The Grid on which to operate

Parameters

bcx

determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.

bcy

analogous to bcx, applies to y direction

Returns

interpolated point

Note

g.contains(x,y) must return true

interpolation() [1/16]

template<class real_type , size_t Nd>

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const aRealTopology< real_type, Nd > &	g_new,
		const aRealTopology< real_type, Nd > &	g_old,
		std::string	method = "dg" )

Create interpolation between two grids.

This matrix interpolates vectors on the old grid g_old to the Gaussian nodes of the new grid g_new. The interpolation is of the order g_old.n()

See also

Introduction to dg methods

for integer multiples between old and new grid you may want to consider the dg::create::fast_interpolation functions

Parameters

g_new	The new grid
g_old	The old grid. The boundaries of the new grid must lie within the boundaries of the old grid.

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

Returns

Interpolation matrix with g_old.size() columns and g_new.size() rows

Attention

Explicit zeros in the returned matrix are removed

interpolation() [2/16]

template<class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const aRealTopologyX2d< real_type > &	g_new,
		const aRealTopologyX2d< real_type > &	g_old )

Create interpolation between two grids.

This matrix can be applied to vectors defined on the old grid to obtain its values on the new grid.

Parameters

g_new	The new points
g_old	The old grid

Returns

Interpolation matrix

Note

The boundaries of the old grid must lie within the boundaries of the new grid

interpolation() [3/16]

template<class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const aRealTopologyX3d< real_type > &	g_new,
		const aRealTopologyX3d< real_type > &	g_old )

Create interpolation between two grids.

This matrix can be applied to vectors defined on the old grid to obtain its values on the new grid.

Parameters

g_new	The new points
g_old	The old grid

Returns

Interpolation matrix

Note

The boundaries of the old grid must lie within the boundaries of the new grid

interpolation() [4/16]

template<class host_vector , class real_type >

dg::MIHMatrix_t< real_type > dg::create::interpolation	(	const host_vector &	x,
		const host_vector &	y,
		const aRealMPITopology2d< real_type > &	g,
		dg::bc	bcx = dg::NEU,
		dg::bc	bcy = dg::NEU,
		std::string	method = "dg" )

Create an MPI row distributed interpolation 2d matrix.

Note

In the MPI version each process creates a local interpolation matrix with local row and global column indices using the given points and

auto mat = dg::create::interpolation ( x,y, g.global(), bcx, bcy, method);
return dg::make_mpi_matrix( mat, g);

The created matrix has g.size() columns and x.size() rows. Per default it uses polynomial interpolation given by the dG polynomials, i.e. the interpolation has order g.n() . When applied to a vector the result contains the interpolated values at the given interpolation points. The given boundary conditions determine how interpolation points outside the grid domain are treated.

        //create equidistant values
        thrust::host_vector<double> x( g.size()), y(x);
        for( unsigned i=0; i<g.Ny()*g.ny(); i++)
            for( unsigned j=0; j<g.Nx()*g.nx(); j++)
            {
                //intentionally set values outside the grid domain
                x[i*g.Nx()*g.nx() + j] =
                        g.x0() + g.lx() + (j+0.5)*g.hx()/(double)(g.nx());
                y[i*g.Nx()*g.ny() + j] =
                        g.y0() + 2*g.ly() + (i+0.5)*g.hy()/(double)(g.ny());
            }
        //use DIR because the coo.2d is zero on the right boundary
        Matrix B = dg::create::interpolation( x, y, g, dg::DIR, dg::DIR, "dg");
        //values outside the grid are mirrored back in
 
        const thrust::host_vector<double> vec = dg::evaluate( function, g);
        thrust::host_vector<double> inter(vec);
        dg::blas2::symv( B, vec, inter);

See also

Introduction to dg methods

Parameters

x	X-coordinates of interpolation points
y	Y-coordinates of interpolation points (y.size() must equal x.size())
g	The Grid on which to operate

Parameters

bcx

determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.

bcy

analogous to bcx, applies to y direction

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

Returns

interpolation matrix

Attention

removes explicit zeros in the interpolation matrix

interpolation() [5/16]

template<class host_vector , class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const host_vector &	x,
		const host_vector &	y,
		const aRealTopology2d< real_type > &	g,
		dg::bc	bcx = dg::NEU,
		dg::bc	bcy = dg::NEU,
		std::string	method = "dg" )

Create interpolation matrix.

The created matrix has g.size() columns and x.size() rows. Per default it uses polynomial interpolation given by the dG polynomials, i.e. the interpolation has order g.n() . When applied to a vector the result contains the interpolated values at the given interpolation points. The given boundary conditions determine how interpolation points outside the grid domain are treated.

        //create equidistant values
        thrust::host_vector<double> x( g.size()), y(x);
        for( unsigned i=0; i<g.Ny()*g.ny(); i++)
            for( unsigned j=0; j<g.Nx()*g.nx(); j++)
            {
                //intentionally set values outside the grid domain
                x[i*g.Nx()*g.nx() + j] =
                        g.x0() + g.lx() + (j+0.5)*g.hx()/(double)(g.nx());
                y[i*g.Nx()*g.ny() + j] =
                        g.y0() + 2*g.ly() + (i+0.5)*g.hy()/(double)(g.ny());
            }
        //use DIR because the coo.2d is zero on the right boundary
        Matrix B = dg::create::interpolation( x, y, g, dg::DIR, dg::DIR, "dg");
        //values outside the grid are mirrored back in
 
        const thrust::host_vector<double> vec = dg::evaluate( function, g);
        thrust::host_vector<double> inter(vec);
        dg::blas2::symv( B, vec, inter);

See also

Introduction to dg methods

Parameters

x	X-coordinates of interpolation points
y	Y-coordinates of interpolation points (y.size() must equal x.size())
g	The Grid on which to operate

Parameters

bcx

determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.

bcy

analogous to bcx, applies to y direction

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

Returns

interpolation matrix

Attention

removes explicit zeros in the interpolation matrix

interpolation() [6/16]

template<class host_vector , class real_type >

dg::MIHMatrix_t< real_type > dg::create::interpolation	(	const host_vector &	x,
		const host_vector &	y,
		const host_vector &	z,
		const aRealMPITopology3d< real_type > &	g,
		dg::bc	bcx = dg::NEU,
		dg::bc	bcy = dg::NEU,
		dg::bc	bcz = dg::PER,
		std::string	method = "dg" )

Create an MPI row distributed interpolation 3d matrix.

Note

In the MPI version each process creates a local interpolation matrix with local row and global column indices using the given points

auto mat = dg::create::interpolation ( x,y,z, g.global(), bcx, bcy, bcz, method);
return dg::make_mpi_matrix( mat, g);

The created matrix has g.size() columns and x.size() rows. Per default it uses polynomial interpolation given by the dG polynomials, i.e. the interpolation has order g.n() . When applied to a vector the result contains the interpolated values at the given interpolation points.

        //create equidistant values
        thrust::host_vector<double> x( g.size()), y(x), z(x);
        for( unsigned k=0; k<g.nz()*g.Nz(); k++)
            for( unsigned i=0; i<g.Ny()*g.ny(); i++)
                for( unsigned j=0; j<g.Nx()*g.nx(); j++)
                {
                    x[(k*g.Ny()*g.ny() + i)*g.Nx()*g.nx() + j] =
                            g.x0() + (j+0.5)*g.hx()/(double)(g.nx());
                    y[(k*g.Ny()*g.ny() + i)*g.Nx()*g.nx() + j] =
                            g.y0() + (i+0.5)*g.hy()/(double)(g.ny());
                    z[(k*g.Ny()*g.ny() + i)*g.Nx()*g.nx() + j] =
                            g.z0() + (k+0.5)*g.hz()/(double)(g.nz());
                }
        auto method = GENERATE( as<std::string>{}, "nearest", "linear", "dg", "cubic");
        Matrix B = dg::create::interpolation( x, y, z, g, dg::DIR, dg::DIR, dg::PER,
                method);
        const thrust::host_vector<double> vec = dg::evaluate( function, g);
        thrust::host_vector<double> inter(vec);
        dg::blas2::symv( B, vec, inter);

See also

Introduction to dg methods

Parameters

x	X-coordinates of interpolation points
y	Y-coordinates of interpolation points (y.size() must equal x.size())
z	Z-coordinates of interpolation points (z.size() must equal x.size())
g	The Grid on which to operate

Parameters

bcx	determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.
bcy	analogous to bcx, applies to y direction
bcz	analogous to bcx, applies to z direction

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

Returns

interpolation matrix

Attention

removes explicit zeros from the interpolation matrix

all points (x, y, z) must lie within or on the boundaries of g

interpolation() [7/16]

template<class host_vector , class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const host_vector &	x,
		const host_vector &	y,
		const host_vector &	z,
		const aRealTopology3d< real_type > &	g,
		dg::bc	bcx = dg::NEU,
		dg::bc	bcy = dg::NEU,
		dg::bc	bcz = dg::PER,
		std::string	method = "dg" )

Create interpolation matrix.

The created matrix has g.size() columns and x.size() rows. Per default it uses polynomial interpolation given by the dG polynomials, i.e. the interpolation has order g.n() . When applied to a vector the result contains the interpolated values at the given interpolation points.

        //create equidistant values
        thrust::host_vector<double> x( g.size()), y(x), z(x);
        for( unsigned k=0; k<g.nz()*g.Nz(); k++)
            for( unsigned i=0; i<g.Ny()*g.ny(); i++)
                for( unsigned j=0; j<g.Nx()*g.nx(); j++)
                {
                    x[(k*g.Ny()*g.ny() + i)*g.Nx()*g.nx() + j] =
                            g.x0() + (j+0.5)*g.hx()/(double)(g.nx());
                    y[(k*g.Ny()*g.ny() + i)*g.Nx()*g.nx() + j] =
                            g.y0() + (i+0.5)*g.hy()/(double)(g.ny());
                    z[(k*g.Ny()*g.ny() + i)*g.Nx()*g.nx() + j] =
                            g.z0() + (k+0.5)*g.hz()/(double)(g.nz());
                }
        auto method = GENERATE( as<std::string>{}, "nearest", "linear", "dg", "cubic");
        Matrix B = dg::create::interpolation( x, y, z, g, dg::DIR, dg::DIR, dg::PER,
                method);
        const thrust::host_vector<double> vec = dg::evaluate( function, g);
        thrust::host_vector<double> inter(vec);
        dg::blas2::symv( B, vec, inter);

See also

Introduction to dg methods

Parameters

x	X-coordinates of interpolation points
y	Y-coordinates of interpolation points (y.size() must equal x.size())
z	Z-coordinates of interpolation points (z.size() must equal x.size())
g	The Grid on which to operate

Parameters

bcx	determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.
bcy	analogous to bcx, applies to y direction
bcz	analogous to bcx, applies to z direction

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

Returns

interpolation matrix

Attention

removes explicit zeros from the interpolation matrix

all points (x, y, z) must lie within or on the boundaries of g

interpolation() [8/16]

template<class host_vector , class real_type , typename = std::enable_if_t<dg::is_vector_v<host_vector>>>

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const host_vector &	x,
		const RealGrid1d< real_type > &	g,
		dg::bc	bcx = dg::NEU,
		std::string	method = "dg" )

Create interpolation matrix.

The created matrix has g.size() columns and x.size() rows. Per default it uses polynomial interpolation given by the dG polynomials, i.e. the interpolation has order g.n() . When applied to a vector the result contains the interpolated values at the given interpolation points. The given boundary conditions determine how interpolation points outside the grid domain are treated.

See also

Introduction to dg methods

Parameters

x	X-coordinates of interpolation points
g	The Grid on which to operate

Parameters

bcx

determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

Returns

interpolation matrix

Attention

removes explicit zeros in the interpolation matrix

interpolation() [9/16]

template<class host_vector , class real_type >

dg::MIHMatrix_t< real_type > dg::create::interpolation	(	const host_vector &	x,
		const RealMPIGrid1d< real_type > &	g,
		dg::bc	bcx = dg::NEU,
		std::string	method = "dg" )

Create an MPI row distributed interpolation 1d matrix.

Note

In the MPI version each process creates a local interpolation matrix with local row and global column indices using the given points and

auto mat = dg::create::interpolation ( x, g.global(), bcx, method);
return dg::make_mpi_matrix( mat, g);

The created matrix has g.size() columns and x.size() rows. Per default it uses polynomial interpolation given by the dG polynomials, i.e. the interpolation has order g.n() . When applied to a vector the result contains the interpolated values at the given interpolation points. The given boundary conditions determine how interpolation points outside the grid domain are treated.

See also

Introduction to dg methods

Parameters

x	X-coordinates of interpolation points
g	The Grid on which to operate

Parameters

bcx

determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

Returns

interpolation matrix

Attention

removes explicit zeros in the interpolation matrix

interpolation() [10/16]

template<class MPITopology , typename = std::enable_if_t<dg::is_vector_v< typename MPITopology::host_vector, MPIVectorTag>>>

dg::MIHMatrix_t< typename MPITopology::value_type > dg::create::interpolation	(	const MPITopology &	g_new,
		const MPITopology &	g_old,
		std::string	method = "dg" )

Create interpolation between two grids.

This matrix interpolates vectors on the old grid g_old to the Gaussian nodes of the new grid g_new. The interpolation is of the order g_old.n()

See also

Introduction to dg methods

for integer multiples between old and new grid you may want to consider the dg::create::fast_interpolation functions

Parameters

g_new	The new grid
g_old	The old grid. The boundaries of the new grid must lie within the boundaries of the old grid.

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

Returns

Interpolation matrix with g_old.size() columns and g_new.size() rows

Attention

Explicit zeros in the returned matrix are removed

interpolation() [11/16]

template<class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const RealGridX1d< real_type > &	g_new,
		const RealGridX1d< real_type > &	g_old )

Create interpolation between two grids.

This matrix can be applied to vectors defined on the old grid to obtain its values on the new grid.

Parameters

g_new	The new points
g_old	The old grid

Returns

Interpolation matrix

Note

The boundaries of the old grid must lie within the boundaries of the new grid

interpolation() [12/16]

template<class RecursiveHostVector , class real_type , size_t Nd>

dg::MIHMatrix_t< real_type > dg::create::interpolation	(	const RecursiveHostVector &	x,
		const aRealMPITopology< real_type, Nd > &	g,
		std::array< dg::bc, Nd >	bcx,
		std::string	method = "dg" )

Create interpolation matrix of a list of points in given grid.

The created matrix has g.size() columns and x[0].size() rows. Per default it uses polynomial interpolation given by the dG polynomials, i.e. the interpolation has order g.n(u) in each direction. When applied to a vector the result contains the interpolated values at the given interpolation points. The given boundary conditions determine how interpolation points outside the grid domain are treated.

See also

Introduction to dg methods

Template Parameters

Nd	Number of dimensions

Parameters

x	Must be of size Nd coordinates of interpolation points (x[0] is the list of x-coordinates, x[1] is the list of y-coordinates, etc.
g	The Grid on which to operate

Parameters

bcx

determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

interpolation() [13/16]

template<class RecursiveHostVector , class real_type , size_t Nd>

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const RecursiveHostVector &	x,
		const aRealTopology< real_type, Nd > &	g,
		std::array< dg::bc, Nd >	bcx,
		std::string	method = "dg" )

Create interpolation matrix of a list of points in given grid.

The created matrix has g.size() columns and x[0].size() rows. Per default it uses polynomial interpolation given by the dG polynomials, i.e. the interpolation has order g.n(u) in each direction. When applied to a vector the result contains the interpolated values at the given interpolation points. The given boundary conditions determine how interpolation points outside the grid domain are treated.

See also

Introduction to dg methods

Template Parameters

Nd	Number of dimensions

Parameters

x	Must be of size Nd coordinates of interpolation points (x[0] is the list of x-coordinates, x[1] is the list of y-coordinates, etc.
g	The Grid on which to operate

Parameters

bcx

determines what to do when a point lies outside the boundary in x. If dg::PER, the point will be shifted topologically back onto the domain. Else the point will be mirrored at the boundary: dg::NEU will then simply interpolate at the resulting point, dg::DIR will take the negative of the interpolation. (dg::DIR_NEU and dg::NEU_DIR apply dg::NEU / dg::DIR to the respective left or right boundary ) This means the result of the interpolation is as if the interpolated function were Fourier transformed with the correct boundary condition and thus extended beyond the grid boundaries. Note that if a point lies directly on the boundary between two grid cells, the value of the polynomial to the right is taken.

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

interpolation() [14/16]

template<class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const thrust::host_vector< real_type > &	x,
		const RealGridX1d< real_type > &	g )

Create interpolation matrix.

The matrix, when applied to a vector, interpolates its values to the given coordinates

Parameters

x	X-coordinates of interpolation points
g	The Grid on which to operate

Returns

interpolation matrix

interpolation() [15/16]

template<class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const thrust::host_vector< real_type > &	x,
		const thrust::host_vector< real_type > &	y,
		const aRealTopologyX2d< real_type > &	g,
		dg::bc	globalbcz = dg::NEU )

Create interpolation matrix.

The matrix, when applied to a vector, interpolates its values to the given coordinates

Parameters

x	X-coordinates of interpolation points
y	Y-coordinates of interpolation points
g	The Grid on which to operate
globalbcz	NEU for common interpolation. DIR for zeros at Box

Returns

interpolation matrix

interpolation() [16/16]

template<class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::interpolation	(	const thrust::host_vector< real_type > &	x,
		const thrust::host_vector< real_type > &	y,
		const thrust::host_vector< real_type > &	z,
		const aRealTopologyX3d< real_type > &	g,
		dg::bc	globalbcz = dg::NEU )

Create interpolation matrix.

The matrix, when applied to a vector, interpolates its values to the given coordinates. In z-direction only a nearest neighbor interpolation is used

Parameters

x	X-coordinates of interpolation points
y	Y-coordinates of interpolation points
z	Z-coordinates of interpolation points
g	The Grid on which to operate
globalbcz	determines what to do if values lie exactly on the boundary

Returns

interpolation matrix

Note

The values of x, y and z must lie within the boundaries of g

projection() [1/5]

template<class real_type , size_t Nd>

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::projection	(	const aRealTopology< real_type, Nd > &	g_new,
		const aRealTopology< real_type, Nd > &	g_old,
		std::string	method = "dg" )

Create a projection between two grids.

This matrix can be applied to vectors defined on the old (fine) grid to obtain its values projected on the new (coarse) grid. (Projection means that the projection integrals over the base polynomials are computed). If the fine grid is a multiple of the coarse grid, the integral value of the projected vector will be conserved and the difference in the L2 norm between old and new vector small. The projection matrix is the adjoint of the interpolation matrix

See also

Introduction to dg methods

for integer multiples between old and new Topological grids and operations you may want to consider the dg::create::fast_projection A large collection

Parameters

g_new	The new (coarse) grid
g_old	The old (fine) grid

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

Returns

Projection matrix

Note

The boundaries of the old grid must lie within the boundaries of the new grid

also check dg::create::transformation, which is the more general solution

Attention

Projection only works if the number of cells in the fine grid is a multiple of the number of cells in the coarse grid and if the number of polynomial coefficients is lower or the same in the new grid

projection() [2/5]

template<class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::projection	(	const aRealTopologyX2d< real_type > &	g_new,
		const aRealTopologyX2d< real_type > &	g_old,
		std::string	method = "dg" )

projection() [3/5]

template<class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::projection	(	const aRealTopologyX3d< real_type > &	g_new,
		const aRealTopologyX3d< real_type > &	g_old,
		std::string	method = "dg" )

projection() [4/5]

template<class MPITopology , typename = std::enable_if_t<dg::is_vector_v< typename MPITopology::host_vector, MPIVectorTag>>>

dg::MIHMatrix_t< typename MPITopology::value_type > dg::create::projection	(	const MPITopology &	g_new,
		const MPITopology &	g_old,
		std::string	method = "dg" )

Create a projection between two grids.

This matrix can be applied to vectors defined on the old (fine) grid to obtain its values projected on the new (coarse) grid. (Projection means that the projection integrals over the base polynomials are computed). If the fine grid is a multiple of the coarse grid, the integral value of the projected vector will be conserved and the difference in the L2 norm between old and new vector small. The projection matrix is the adjoint of the interpolation matrix

See also

Introduction to dg methods

for integer multiples between old and new Topological grids and operations you may want to consider the dg::create::fast_projection A large collection

Parameters

g_new	The new (coarse) grid
g_old	The old (fine) grid

Parameters

method

Several interpolation methods are available: dg uses the native dG interpolation scheme given by the grid, nearest searches for the nearest point and copies its value, linear searches for the two (in 2d four, etc.) closest points and linearly interpolates their values, cubic searches for the four (in 2d 16, etc) closest points and interpolates a cubic polynomial. Pay attention that linear and cubic entail nearest neighbor communication in mpi.

Returns

Projection matrix

Note

The boundaries of the old grid must lie within the boundaries of the new grid

also check dg::create::transformation, which is the more general solution

Attention

Projection only works if the number of cells in the fine grid is a multiple of the number of cells in the coarse grid and if the number of polynomial coefficients is lower or the same in the new grid

projection() [5/5]

template<class real_type >

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::projection	(	const RealGridX1d< real_type > &	g_new,
		const RealGridX1d< real_type > &	g_old,
		std::string	method = "dg" )

transformation()

template<class real_type , size_t Nd>

dg::SparseMatrix< int, real_type, thrust::host_vector > dg::create::transformation	(	const aRealTopology< real_type, Nd > &	g_new,
		const aRealTopology< real_type, Nd > &	g_old )

Create a transformation matrix between two grids.

The transformation matrix is probably the most correct way of transforming dG vectors between any two grids of different resolution. It first finds the least common multiple grid (lcm) of the old and the new grid. Then it interpolates the values to the lcm grid and finally projects them back to the new grid. In total we have

\[ \mathcal T = P Q \]

where \( Q\) is the interpolation matrix and \( P \) the projection. If either new or old grid is already the lcm grid this function reduces to the interpolation/projection function.

See also

Introduction to dg methods

Parameters

g_new	The new grid
g_old	The old grid

Returns

transformation matrix

Note

The boundaries of the old grid must lie within the boundaries of the new grid

If the grid are very incompatible the matrix-matrix multiplication can take a while