Discontinuous Galerkin Library
#include "dg/algorithm.h"
gridX.h
Go to the documentation of this file.
1#pragma once
2#include <cassert>
3#include <cmath> //floor function
4#include "dlt.h"
5#include "grid.h"
6#include "../enums.h"
7
45namespace dg{
47template<class real_type>
48struct RealGridX2d;
49template<class real_type>
50struct RealGridX3d;
52
66template<class real_type>
68{
69 using value_type = real_type;
71 using host_vector = thrust::host_vector<real_type>;
83 RealGridX1d( real_type x0, real_type x1, real_type f, unsigned n, unsigned N, bc bcx = NEU):
84 x0_(x0), x1_(x1), f_(f),
85 n_(n), Nx_(N), bcx_(bcx), dlt_(n)
86 {
87 assert( (f >= 0) && (f < 0.5) );
88 assert( fabs(outer_N() - f*(real_type)N) < 1e-14);
89 assert( x1 > x0 );
90 assert( N > 0 );
91 assert( n != 0 );
92 assert( bcx != PER);
93 }
99 real_type x0() const {return x0_;}
105 real_type x1() const {return x1_;}
111 real_type f() const {return f_;}
117 real_type lx() const {return x1_-x0_;}
123 real_type h() const {return lx()/(real_type)Nx_;}
129 unsigned N() const {return Nx_;}
135 unsigned outer_N() const {return (unsigned)(round(f_*(real_type)Nx_));}
141 unsigned inner_N() const {return N()-2*outer_N();}
147 unsigned n() const {return n_;}
153 bc bcx() const {return bcx_;}
159 unsigned size() const { return n_*Nx_;}
165 void display( std::ostream& os = std::cout) const
166 {
167 os << "RealGrid parameters are: \n"
168 <<" n = "<<n_<<"\n"
169 <<" N = "<<Nx_<<"\n"
170 <<" inner N = "<<inner_N()<<"\n"
171 <<" outer N = "<<outer_N()<<"\n"
172 <<" h = "<<h()<<"\n"
173 <<" x0 = "<<x0_<<"\n"
174 <<" x1 = "<<x1_<<"\n"
175 <<" lx = "<<lx()<<"\n"
176 <<"Boundary conditions in x are: \n"
177 <<" "<<bc2str(bcx_)<<"\n";
178 }
184 const DLT<real_type>& dlt() const {return dlt_;}
185 RealGrid1d<real_type> grid() const{return RealGrid1d<real_type>( x0_, x1_, n_, Nx_, bcx_);}
186
195 void shift_topologic( real_type x0, real_type& x1) const
196 {
197 assert( contains(x0));
198 real_type deltaX;
199 real_type xleft = x0_ + f_*lx();
200 real_type xright = x1_ - f_*lx();
201 if( x0 >= xleft && x0<xright)
202 {
203 if( x1 > xleft) deltaX = (x1 -xleft);
204 else deltaX = xright - x1;
205 unsigned N = floor(deltaX/(xright-xleft));
206 if( x1 > xright ) x1 -= N*lx();
207 if( x1 < xleft ) x1 += N*lx();
208 }
209 else if( x0 < xleft && x1 >=xleft)
210 x1 += (xright-xleft);
211 else if( x0 >= xright && x1 < xright)
212 x1 -= (xright-xleft);
213
214 }
215
224 bool contains( real_type x) const
225 {
226 if( (x>=x0_ && x <= x1_)) return true;
227 return false;
228 }
229 private:
230 real_type x0_, x1_, f_;
231 unsigned n_, Nx_;
232 bc bcx_;
233 DLT<real_type> dlt_;
234};
235
236//template<class real_type>
237//struct aRealTopologyX3d; //forward declare 3d version
238
255template<class real_type>
257{
258 using value_type = real_type;
260 using host_vector = thrust::host_vector<real_type>;
262
268 real_type x0() const {return x0_;}
274 real_type x1() const {return x1_;}
280 real_type y0() const {return y0_;}
286 real_type y1() const {return y1_;}
292 real_type lx() const {return x1_-x0_;}
298 real_type ly() const {return y1_-y0_;}
304 real_type hx() const {return lx()/(real_type)Nx_;}
310 real_type hy() const {return ly()/(real_type)Ny_;}
316 real_type fx() const {return fx_;}
322 real_type fy() const {return fy_;}
328 unsigned n() const {return n_;}
334 unsigned Nx() const {return Nx_;}
340 unsigned inner_Nx() const {return Nx_ - outer_Nx();}
346 unsigned outer_Nx() const {return (unsigned)round(fx_*(real_type)Nx_);}
352 unsigned Ny() const {return Ny_;}
358 unsigned inner_Ny() const {return Ny_-2*outer_Ny();}
364 unsigned outer_Ny() const {return (unsigned)round(fy_*(real_type)Ny_);}
370 bc bcx() const {return bcx_;}
376 bc bcy() const {return bcy_;}
382 RealGrid2d<real_type> grid() const {return RealGrid2d<real_type>( x0_,x1_,y0_,y1_,n_,Nx_,Ny_,bcx_,bcy_);}
388 const DLT<real_type>& dlt() const{return dlt_;}
394 unsigned size() const { return n_*n_*Nx_*Ny_;}
400 void display( std::ostream& os = std::cout) const
401 {
402 os << "Grid parameters are: \n"
403 <<" n = "<<n_<<"\n"
404 <<" Nx = "<<Nx_<<"\n"
405 <<" inner Nx = "<<inner_Nx()<<"\n"
406 <<" outer Nx = "<<outer_Nx()<<"\n"
407 <<" Ny = "<<Ny_<<"\n"
408 <<" inner Ny = "<<inner_Ny()<<"\n"
409 <<" outer Ny = "<<outer_Ny()<<"\n"
410 <<" hx = "<<hx()<<"\n"
411 <<" hy = "<<hy()<<"\n"
412 <<" x0 = "<<x0_<<"\n"
413 <<" x1 = "<<x1_<<"\n"
414 <<" y0 = "<<y0_<<"\n"
415 <<" y1 = "<<y1_<<"\n"
416 <<" lx = "<<lx()<<"\n"
417 <<" ly = "<<ly()<<"\n"
418 <<"Boundary conditions in x are: \n"
419 <<" "<<bc2str(bcx_)<<"\n"
420 <<"Boundary conditions in y are: \n"
421 <<" "<<bc2str(bcy_)<<"\n";
422 }
423
434 void shift_topologic( real_type x0, real_type y0, real_type& x1, real_type& y1) const
435 {
436 assert( contains(x0, y0));
437 real_type deltaX;
438 if( x1 > x0_) deltaX = (x1 -x0_);
439 else deltaX = x1_ - x1;
440 unsigned N = floor(deltaX/lx());
441 if( x1 > x1_ && bcx_ == dg::PER) x1 -= N*lx();
442 if( x1 < x0_ && bcx_ == dg::PER) x1 += N*lx();
443
444 if( x0 < x1_ - fx_*(x1_-x0_) ) //if x0 is one of the inner points
445 {
446 real_type deltaY;
447 real_type yleft = y0_ + fy_*ly();
448 real_type yright = y1_ - fy_*ly();
449 if( y0 >= yleft && y0<yright)
450 {
451 if( y1 > yleft) deltaY = (y1 -yleft);
452 else deltaY = yright - y1;
453 unsigned N = floor(deltaY/(yright-yleft));
454 if( y1 > yright ) y1 -= N*ly();
455 if( y1 < yleft ) y1 += N*ly();
456 }
457 else if( y0 < yleft && y1 >=yleft)
458 y1 += (yright-yleft);
459 else if( y0 >= yright && y1 < yright)
460 y1 -= (yright-yleft);
461 }
462
463 }
464
474 bool contains( real_type x, real_type y)const
475 {
476 if( (x>=x0_ && x <= x1_) && (y>=y0_ && y <= y1_)) return true;
477 return false;
478 }
479 protected:
481 ~aRealTopologyX2d() = default;
484 aRealTopologyX2d( real_type x0, real_type x1, real_type y0, real_type y1, real_type fx, real_type fy, unsigned n, unsigned Nx, unsigned Ny, bc bcx, bc bcy):
485 x0_(x0), x1_(x1), y0_(y0), y1_(y1), fx_(fx), fy_(fy),
486 n_(n), Nx_(Nx), Ny_(Ny), bcx_(bcx), bcy_( bcy), dlt_(n)
487 {
488 assert( (fy_ >= 0.) && (fy_ < 0.5) );
489 assert( (fx_ >= 0.) && (fx_ < 1.) );
490 assert( fabs(outer_Nx() - fx_*(real_type)Nx) < 1e-14);
491 assert( fabs(outer_Ny() - fy_*(real_type)Ny) < 1e-14);
492 assert( n != 0);
493 assert( x1 > x0 && y1 > y0);
494 assert( Nx_ > 0 && Ny > 0 );
495 assert( bcy != PER);
496 }
498 aRealTopologyX2d(const aRealTopologyX2d& src) = default;
501 private:
502 real_type x0_, x1_, y0_, y1_;
503 real_type fx_, fy_;
504 unsigned n_, Nx_, Ny_;
505 bc bcx_, bcy_;
506 DLT<real_type> dlt_;
507};
512template<class real_type>
513struct RealGridX2d : public aRealTopologyX2d<real_type>
514{
517 RealGridX2d( real_type x0, real_type x1, real_type y0, real_type y1, real_type fx, real_type fy, unsigned n, unsigned Nx, unsigned Ny, bc bcx=PER, bc bcy=NEU):
518 aRealTopologyX2d<real_type>(x0,x1,y0,y1,fx,fy,n,Nx,Ny,bcx,bcy) { }
520 explicit RealGridX2d( const aRealTopologyX2d<real_type>& src): aRealTopologyX2d<real_type>(src){}
521};
522
539template<class real_type>
541{
542 using value_type = real_type;
544 using host_vector = thrust::host_vector<real_type>;
551 real_type x0() const {return x0_;}
557 real_type x1() const {return x1_;}
558
564 real_type y0() const {return y0_;}
570 real_type y1() const {return y1_;}
571
577 real_type z0() const {return z0_;}
583 real_type z1() const {return z1_;}
584
590 real_type lx() const {return x1_-x0_;}
596 real_type ly() const {return y1_-y0_;}
602 real_type lz() const {return z1_-z0_;}
603
609 real_type hx() const {return lx()/(real_type)Nx_;}
615 real_type hy() const {return ly()/(real_type)Ny_;}
621 real_type hz() const {return lz()/(real_type)Nz_;}
627 real_type fx() const {return fx_;}
633 real_type fy() const {return fy_;}
639 unsigned n() const {return n_;}
645 unsigned Nx() const {return Nx_;}
651 unsigned inner_Nx() const {return Nx_ - outer_Nx();}
657 unsigned outer_Nx() const {return (unsigned)round(fx_*(real_type)Nx_);}
663 unsigned Ny() const {return Ny_;}
669 unsigned inner_Ny() const {return Ny_-2*outer_Ny();}
675 unsigned outer_Ny() const {return (unsigned)round(fy_*(real_type)Ny_);}
681 unsigned Nz() const {return Nz_;}
687 bc bcx() const {return bcx_;}
693 bc bcy() const {return bcy_;}
699 bc bcz() const {return bcz_;}
706 return RealGrid3d<real_type>( x0_,x1_,y0_,y1_,z0_,z1_,n_,Nx_,Ny_,Nz_,bcx_,bcy_,bcz_);
707 }
713 const DLT<real_type>& dlt() const{return dlt_;}
719 unsigned size() const { return n_*n_*Nx_*Ny_*Nz_;}
725 void display( std::ostream& os = std::cout) const
726 {
727 os << "Grid parameters are: \n"
728 <<" n = "<<n_<<"\n"
729 <<" Nx = "<<Nx_<<"\n"
730 <<" inner Nx = "<<inner_Nx()<<"\n"
731 <<" outer Nx = "<<outer_Nx()<<"\n"
732 <<" Ny = "<<Ny_<<"\n"
733 <<" inner Ny = "<<inner_Ny()<<"\n"
734 <<" outer Ny = "<<outer_Ny()<<"\n"
735 <<" Nz = "<<Nz_<<"\n"
736 <<" hx = "<<hx()<<"\n"
737 <<" hy = "<<hy()<<"\n"
738 <<" hz = "<<hz()<<"\n"
739 <<" x0 = "<<x0_<<"\n"
740 <<" x1 = "<<x1_<<"\n"
741 <<" y0 = "<<y0_<<"\n"
742 <<" y1 = "<<y1_<<"\n"
743 <<" z0 = "<<z0_<<"\n"
744 <<" z1 = "<<z1_<<"\n"
745 <<" lx = "<<lx()<<"\n"
746 <<" ly = "<<ly()<<"\n"
747 <<" lz = "<<lz()<<"\n"
748 <<"Boundary conditions in x are: \n"
749 <<" "<<bc2str(bcx_)<<"\n"
750 <<"Boundary conditions in y are: \n"
751 <<" "<<bc2str(bcy_)<<"\n"
752 <<"Boundary conditions in z are: \n"
753 <<" "<<bc2str(bcz_)<<"\n";
754 }
765 bool contains( real_type x, real_type y, real_type z)const
766 {
767 if( (x>=x0_ && x <= x1_) && (y>=y0_ && y <= y1_) && (z>=z0_ && z<=z1_))
768 return true;
769 return false;
770 }
771 protected:
773 ~aRealTopologyX3d() = default;
776 aRealTopologyX3d( real_type x0, real_type x1, real_type y0, real_type y1, real_type z0, real_type z1, real_type fx, real_type fy, unsigned n, unsigned Nx, unsigned Ny, unsigned Nz, bc bcx, bc bcy, bc bcz):
777 x0_(x0), x1_(x1), y0_(y0), y1_(y1), z0_(z0), z1_(z1), fx_(fx), fy_(fy),
778 n_(n), Nx_(Nx), Ny_(Ny), Nz_(Nz), bcx_(bcx), bcy_( bcy), bcz_( bcz), dlt_(n)
779 {
780 assert( (fy_ >= 0.) && (fy_ < 0.5) );
781 assert( (fx_ >= 0.) && (fx_ < 1.) );
782 assert( fabs(outer_Nx() - fx_*(real_type)Nx) < 1e-14);
783 assert( fabs(outer_Ny() - fy_*(real_type)Ny) < 1e-14);
784 assert( n != 0);
785 assert( x1 > x0 && y1 > y0 ); assert( z1 > z0 );
786 assert( Nx_ > 0 && Ny > 0); assert( Nz > 0);
787 }
789 aRealTopologyX3d(const aRealTopologyX3d& src) = default;
792 private:
793 real_type x0_, x1_, y0_, y1_, z0_, z1_;
794 real_type fx_,fy_;
795 unsigned n_, Nx_, Ny_, Nz_;
796 bc bcx_, bcy_, bcz_;
797 DLT<real_type> dlt_;
798};
799
804template<class real_type>
805struct RealGridX3d : public aRealTopologyX3d<real_type>
806{
809 RealGridX3d( real_type x0, real_type x1, real_type y0, real_type y1, real_type z0, real_type z1, real_type fx, real_type fy, unsigned n, unsigned Nx, unsigned Ny, unsigned Nz, bc bcx=PER, bc bcy=NEU, bc bcz=PER):
810 aRealTopologyX3d<real_type>(x0,x1,y0,y1,z0,z1,fx,fy,n,Nx,Ny,Nz,bcx,bcy,bcz) { }
812 explicit RealGridX3d( const aRealTopologyX3d<real_type>& src): aRealTopologyX3d<real_type>(src){}
813};
814
823
824}// namespace dg
Struct holding coefficients for Discrete Legendre Transformation (DLT) related operations.
Definition: dlt.h:23
The discrete legendre trafo class.
base topology classes
bc
Switch between boundary conditions.
Definition: enums.h:15
static std::string bc2str(bc bcx)
write a string describing boundary condition to an output stream
Definition: enums.h:37
@ z
z direction
@ PER
periodic boundaries
Definition: enums.h:16
@ NEU
Neumann on both boundaries.
Definition: enums.h:20
@ y
y direction
@ x
x direction
This is the namespace for all functions and classes defined and used by the discontinuous Galerkin li...
1D grid
Definition: grid.h:80
The simplest implementation of aRealTopology2d.
Definition: grid.h:818
The simplest implementation of aRealTopology3d.
Definition: grid.h:844
1D grid for X-point topology
Definition: gridX.h:68
RealGridX1d(real_type x0, real_type x1, real_type f, unsigned n, unsigned N, bc bcx=NEU)
1D X-point grid
Definition: gridX.h:83
unsigned inner_N() const
number of cells in the inner region
Definition: gridX.h:141
const DLT< real_type > & dlt() const
the discrete legendre transformation
Definition: gridX.h:184
bool contains(real_type x) const
Check if the grid contains a point.
Definition: gridX.h:224
bc bcx() const
boundary conditions
Definition: gridX.h:153
real_type value_type
Definition: gridX.h:69
unsigned size() const
the total number of points
Definition: gridX.h:159
unsigned outer_N() const
number of cells in one of the outer regions
Definition: gridX.h:135
unsigned N() const
number of cells
Definition: gridX.h:129
void shift_topologic(real_type x0, real_type &x1) const
Shifts a point coordinate due to topology.
Definition: gridX.h:195
void display(std::ostream &os=std::cout) const
Display.
Definition: gridX.h:165
real_type x0() const
left boundary
Definition: gridX.h:99
unsigned n() const
number of polynomial coefficients
Definition: gridX.h:147
real_type x1() const
right boundary
Definition: gridX.h:105
real_type f() const
Factor.
Definition: gridX.h:111
thrust::host_vector< real_type > host_vector
The host vector type used by host functions like evaluate.
Definition: gridX.h:71
real_type lx() const
total length of interval
Definition: gridX.h:117
real_type h() const
cell size
Definition: gridX.h:123
RealGrid1d< real_type > grid() const
Definition: gridX.h:185
The simplest implementation of aRealTopologyX2d.
Definition: gridX.h:514
RealGridX2d(const aRealTopologyX2d< real_type > &src)
allow explicit type conversion from any other topology
Definition: gridX.h:520
RealGridX2d(real_type x0, real_type x1, real_type y0, real_type y1, real_type fx, real_type fy, unsigned n, unsigned Nx, unsigned Ny, bc bcx=PER, bc bcy=NEU)
Construct a 2D X-point grid.
Definition: gridX.h:517
The simplest implementation of aRealTopologyX3d.
Definition: gridX.h:806
RealGridX3d(real_type x0, real_type x1, real_type y0, real_type y1, real_type z0, real_type z1, real_type fx, real_type fy, unsigned n, unsigned Nx, unsigned Ny, unsigned Nz, bc bcx=PER, bc bcy=NEU, bc bcz=PER)
Construct a 3D X-point grid.
Definition: gridX.h:809
RealGridX3d(const aRealTopologyX3d< real_type > &src)
allow explicit type conversion from any other topology
Definition: gridX.h:812
A 2D grid class with X-point topology.
Definition: gridX.h:257
thrust::host_vector< real_type > host_vector
The host vector type used by host functions like evaluate.
Definition: gridX.h:260
unsigned size() const
real_typehe total number of points
Definition: gridX.h:394
bool contains(real_type x, real_type y) const
Check if the grid contains a point.
Definition: gridX.h:474
real_type hx() const
cell size in x
Definition: gridX.h:304
real_type y0() const
left boundary in y
Definition: gridX.h:280
bc bcy() const
boundary conditions in y
Definition: gridX.h:376
real_type lx() const
length of x
Definition: gridX.h:292
RealGrid2d< real_type > grid() const
Return a copy without topology.
Definition: gridX.h:382
real_type y1() const
Right boundary in y.
Definition: gridX.h:286
aRealTopologyX2d & operator=(const aRealTopologyX2d &src)=default
unsigned Nx() const
number of cells in x
Definition: gridX.h:334
real_type x1() const
Right boundary in x.
Definition: gridX.h:274
real_type x0() const
Left boundary in x.
Definition: gridX.h:268
void display(std::ostream &os=std::cout) const
Display.
Definition: gridX.h:400
real_type fx() const
partition factor in x
Definition: gridX.h:316
aRealTopologyX2d(real_type x0, real_type x1, real_type y0, real_type y1, real_type fx, real_type fy, unsigned n, unsigned Nx, unsigned Ny, bc bcx, bc bcy)
Construct a 2D X-point grid.
Definition: gridX.h:484
bc bcx() const
boundary conditions in x
Definition: gridX.h:370
unsigned inner_Nx() const
number of topological cells in x
Definition: gridX.h:340
unsigned Ny() const
number of cells in y
Definition: gridX.h:352
real_type ly() const
length of y
Definition: gridX.h:298
real_type fy() const
partition factor in y
Definition: gridX.h:322
unsigned outer_Ny() const
number of cells in one of the outer regions of y
Definition: gridX.h:364
real_type hy() const
cell size in y
Definition: gridX.h:310
unsigned n() const
number of polynomial coefficients in x and y
Definition: gridX.h:328
unsigned outer_Nx() const
number of smooth rows in x
Definition: gridX.h:346
void shift_topologic(real_type x0, real_type y0, real_type &x1, real_type &y1) const
Shifts a point coordinate due to topology.
Definition: gridX.h:434
unsigned inner_Ny() const
number of cells in the inner region of y
Definition: gridX.h:358
aRealTopologyX2d(const aRealTopologyX2d &src)=default
~aRealTopologyX2d()=default
disallow destruction through base class pointer
real_type value_type
Definition: gridX.h:258
const DLT< real_type > & dlt() const
discrete legendre trafo
Definition: gridX.h:388
A 3D grid class with X-point topology.
Definition: gridX.h:541
bc bcx() const
boundary conditions in x
Definition: gridX.h:687
RealGrid3d< real_type > grid() const
Return a copy without topology.
Definition: gridX.h:705
aRealTopologyX3d & operator=(const aRealTopologyX3d &src)=default
unsigned n() const
number of polynomial coefficients in x and y
Definition: gridX.h:639
real_type y1() const
right boundary in y
Definition: gridX.h:570
real_type z1() const
right boundary in z
Definition: gridX.h:583
real_type lz() const
length in z
Definition: gridX.h:602
real_type hx() const
cell size in x
Definition: gridX.h:609
aRealTopologyX3d(real_type x0, real_type x1, real_type y0, real_type y1, real_type z0, real_type z1, real_type fx, real_type fy, unsigned n, unsigned Nx, unsigned Ny, unsigned Nz, bc bcx, bc bcy, bc bcz)
Construct a 3D X-point grid.
Definition: gridX.h:776
void display(std::ostream &os=std::cout) const
Display.
Definition: gridX.h:725
const DLT< real_type > & dlt() const
discrete legendre transformation
Definition: gridX.h:713
real_type fx() const
partition factor in x
Definition: gridX.h:627
unsigned Nz() const
number of points in z
Definition: gridX.h:681
unsigned inner_Ny() const
number of cells in the inner region of y
Definition: gridX.h:669
real_type hy() const
cell size in y
Definition: gridX.h:615
unsigned Ny() const
number of cells in y
Definition: gridX.h:663
real_type hz() const
cell size in z
Definition: gridX.h:621
unsigned Nx() const
number of points in x
Definition: gridX.h:645
bc bcy() const
boundary conditions in y
Definition: gridX.h:693
unsigned size() const
real_typehe total number of points
Definition: gridX.h:719
~aRealTopologyX3d()=default
disallow destruction through base class pointer
real_type z0() const
left boundary in z
Definition: gridX.h:577
real_type y0() const
left boundary in y
Definition: gridX.h:564
thrust::host_vector< real_type > host_vector
The host vector type used by host functions like evaluate.
Definition: gridX.h:544
unsigned inner_Nx() const
number of topological cells in x
Definition: gridX.h:651
real_type value_type
Definition: gridX.h:542
real_type x0() const
left boundary in x
Definition: gridX.h:551
unsigned outer_Ny() const
number of cells in one of the outer regions of y
Definition: gridX.h:675
real_type x1() const
right boundary in x
Definition: gridX.h:557
real_type ly() const
length in y
Definition: gridX.h:596
real_type fy() const
partition factor in y
Definition: gridX.h:633
bc bcz() const
boundary conditions in z
Definition: gridX.h:699
bool contains(real_type x, real_type y, real_type z) const
Check if the grid contains a point.
Definition: gridX.h:765
aRealTopologyX3d(const aRealTopologyX3d &src)=default
real_type lx() const
length in x
Definition: gridX.h:590
unsigned outer_Nx() const
number of smooth rows in x
Definition: gridX.h:657