Discontinuous Galerkin Library
#include "dg/algorithm.h"
interpolation.h
Go to the documentation of this file.
1#pragma once
2//#include <iomanip>
3
4#include <cusp/coo_matrix.h>
5#include <cusp/csr_matrix.h>
6#include "dg/backend/typedefs.h"
7#include "grid.h"
8#include "evaluation.h"
9#include "functions.h"
10#include "operator_tensor.h"
11#include "xspacelib.h"
12
18namespace dg{
19
20namespace create{
22namespace detail{
23
32template<class real_type>
33std::vector<real_type> coefficients( real_type xn, unsigned n)
34{
35 assert( xn <= 1. && xn >= -1.);
36 std::vector<real_type> px(n);
37 if( xn == -1)
38 {
39 for( unsigned u=0; u<n; u++)
40 px[u] = (u%2 == 0) ? +1. : -1.;
41 }
42 else if( xn == 1)
43 {
44 for( unsigned i=0; i<n; i++)
45 px[i] = 1.;
46 }
47 else
48 {
49 px[0] = 1.;
50 if( n > 1)
51 {
52 px[1] = xn;
53 for( unsigned i=1; i<n-1; i++)
54 px[i+1] = ((real_type)(2*i+1)*xn*px[i]-(real_type)i*px[i-1])/(real_type)(i+1);
55 }
56 }
57 return px;
58}
59
60// evaluate n base polynomials for n given abscissas
61template<class real_type>
62std::vector<real_type> lagrange( real_type x, const std::vector<real_type>& xi)
63{
64 unsigned n = xi.size();
65 std::vector<real_type> l( n , 1.);
66 for( unsigned i=0; i<n; i++)
67 for( unsigned k=0; k<n; k++)
68 {
69 if ( k != i)
70 l[i] *= (x-xi[k])/(xi[i]-xi[k]);
71 }
72 return l;
73}
74
75//THERE IS A BUG FOR PERIODIC BC !!
76template<class real_type>
77std::vector<real_type> choose_1d_abscissas( real_type X,
78 unsigned points_per_line, const RealGrid1d<real_type>& g,
79 const thrust::host_vector<real_type>& abs,
80 thrust::host_vector<unsigned>& cols)
81{
82 assert( abs.size() >= points_per_line && "There must be more points to interpolate\n");
83 dg::bc bcx = g.bcx();
84 //determine which cell (X) lies in
85 real_type xnn = (X-g.x0())/g.h();
86 unsigned n = (unsigned)floor(xnn);
87 //intervall correction
88 if (n==g.N() && bcx != dg::PER) {
89 n-=1;
90 }
91 // look for closest abscissa
92 std::vector<real_type> xs( points_per_line, 0);
93 // X <= *it
94 auto it = std::lower_bound( abs.begin()+n*g.n(), abs.begin() + (n+1)*g.n(),
95 X);
96 cols.resize( points_per_line, 0);
97 switch( points_per_line)
98 {
99 case 1: xs[0] = 1.;
100 if( it == abs.begin())
101 cols[0] = 0;
102 else if( it == abs.end())
103 cols[0] = it - abs.begin() - 1;
104 else
105 {
106 if ( fabs(X - *it) < fabs( X - *(it-1)))
107 cols[0] = it - abs.begin();
108 else
109 cols[0] = it - abs.begin() -1;
110 }
111 break;
112 case 2: if( it == abs.begin())
113 {
114 if( bcx == dg::PER)
115 {
116 xs[0] = *it;
117 xs[1] = *(abs.end() -1)-g.lx();;
118 cols[0] = 0, cols[1] = abs.end()-abs.begin()-1;
119 }
120 else
121 {
122 //xs[0] = *it;
123 //xs[1] = *(it+1);
124 //cols[0] = 0, cols[1] = 1;
125 // This makes it consistent with fem_t
126 xs.resize(1);
127 xs[0] = *it;
128 cols[0] = 0;
129 }
130 }
131 else if( it == abs.end())
132 {
133 if( bcx == dg::PER)
134 {
135 xs[0] = *(abs.begin())+g.lx();
136 xs[1] = *(it-1);
137 cols[0] = 0, cols[1] = it-abs.begin()-1;
138 }
139 else
140 {
141 //xs[0] = *(it-2);
142 //xs[0] = *(it-1);
143 //cols[0] = it - abs.begin() - 2;
144 //cols[1] = it - abs.begin() - 1;
145 // This makes it consistent with fem_t
146 xs.resize(1);
147 xs[0] = *(it-1);
148 cols[0] = it-abs.begin()-1;
149 }
150 }
151 else
152 {
153 xs[0] = *(it-1);
154 xs[1] = *it;
155 cols[0] = it - abs.begin() - 1;
156 cols[1] = cols[0]+1;
157 }
158 break;
159 case 4: if( it <= abs.begin() +1)
160 {
161 if( bcx == dg::PER)
162 {
163 xs[0] = *abs.begin(), cols[0] = 0;
164 xs[1] = *(abs.begin()+1), cols[1] = 1;
165 xs[2] = it == abs.begin() ? *(abs.end() -2) : *(abs.begin()+2);
166 cols[2] = it == abs.begin() ? abs.end()-abs.begin() -2 : 2;
167 xs[3] = *(abs.end() -1);
168 cols[3] = abs.end()-abs.begin() -1;
169 }
170 else
171 {
172 it = abs.begin();
173 xs[0] = *it, xs[1] = *(it+1);
174 xs[2] = *(it+2), xs[3] = *(it+3);
175 cols[0] = 0, cols[1] = 1;
176 cols[2] = 2, cols[3] = 3;
177 }
178 }
179 else if( it >= abs.end() -1)
180 {
181 if( bcx == dg::PER)
182 {
183 xs[0] = *abs.begin(), cols[0] = 0;
184 xs[1] = it == abs.end() ? *(abs.begin()+1) : *(abs.end() -3) ;
185 cols[1] = it == abs.end() ? 1 : abs.end()-abs.begin()-3 ;
186 xs[2] = *(abs.end() - 2), cols[2] = abs.end()-abs.begin()-2;
187 xs[3] = *(abs.end() - 1), cols[3] = abs.end()-abs.begin()-1;
188 }
189 else
190 {
191 it = abs.end();
192 xs[0] = *(it-4), xs[1] = *(it-3);
193 xs[2] = *(it-2), xs[3] = *(it-1);
194 cols[0] = it - abs.begin() - 4;
195 cols[1] = cols[0]+1;
196 cols[2] = cols[1]+1;
197 cols[3] = cols[2]+1;
198 }
199 }
200 else
201 {
202 xs[0] = *(it-2), xs[1] = *(it-1);
203 xs[2] = *(it ), xs[3] = *(it+1);
204 cols[0] = it - abs.begin() - 2;
205 cols[1] = cols[0]+1;
206 cols[2] = cols[1]+1;
207 cols[3] = cols[2]+1;
208 }
209 break;
210 }
211 return xs;
212}
213
214}//namespace detail
218
253template<class real_type>
254cusp::coo_matrix<int, real_type, cusp::host_memory> interpolation(
255 const thrust::host_vector<real_type>& x,
256 const RealGrid1d<real_type>& g,
257 dg::bc bcx = dg::NEU,
258 std::string method = "dg")
259{
260 cusp::array1d<real_type, cusp::host_memory> values;
261 cusp::array1d<int, cusp::host_memory> row_indices;
262 cusp::array1d<int, cusp::host_memory> column_indices;
263 if( method == "dg")
264 {
265 dg::Operator<real_type> forward( g.dlt().forward());
266 for( unsigned i=0; i<x.size(); i++)
267 {
268 real_type X = x[i];
269 bool negative = false;
270 g.shift( negative, X, bcx);
271
272 //determine which cell (x) lies in
273 real_type xnn = (X-g.x0())/g.h();
274 unsigned n = (unsigned)floor(xnn);
275 //determine normalized coordinates
276 real_type xn = 2.*xnn - (real_type)(2*n+1);
277 //intervall correction
278 if (n==g.N()) {
279 n-=1;
280 xn = 1.;
281 }
282 //evaluate 2d Legendre polynomials at (xn, yn)...
283 std::vector<real_type> px = detail::coefficients( xn, g.n());
284 //...these are the matrix coefficients with which to multiply
285 std::vector<real_type> pxF(px.size(),0);
286 for( unsigned l=0; l<g.n(); l++)
287 for( unsigned k=0; k<g.n(); k++)
288 pxF[l]+= px[k]*forward(k,l);
289 unsigned cols = n*g.n();
290 for ( unsigned l=0; l<g.n(); l++)
291 {
292 row_indices.push_back(i);
293 column_indices.push_back( cols + l);
294 values.push_back(negative ? -pxF[l] : pxF[l]);
295 }
296 }
297 }
298 else
299 {
300 unsigned points_per_line = 1;
301 if( method == "nearest")
302 points_per_line = 1;
303 else if( method == "linear")
304 points_per_line = 2;
305 else if( method == "cubic")
306 points_per_line = 4;
307 else
308 throw std::runtime_error( "Interpolation method "+method+" not recognized!\n");
309 thrust::host_vector<real_type> abs = dg::create::abscissas( g);
310 dg::RealGrid1d<real_type> gx( g.x0(), g.x1(), g.n(), g.N(), bcx);
311 for( unsigned i=0; i<x.size(); i++)
312 {
313 real_type X = x[i];
314 bool negative = false;
315 g.shift( negative, X, bcx);
316
317 thrust::host_vector<unsigned> cols;
318 std::vector<real_type> xs = detail::choose_1d_abscissas( X,
319 points_per_line, gx, abs, cols);
320
321 std::vector<real_type> px = detail::lagrange( X, xs);
322 // px may have size != points_per_line (at boundary)
323 for ( unsigned l=0; l<px.size(); l++)
324 {
325 row_indices.push_back(i);
326 column_indices.push_back( cols[l]);
327 values.push_back(negative ? -px[l] : px[l]);
328 }
329 }
330 }
331 cusp::coo_matrix<int, real_type, cusp::host_memory> A(
332 x.size(), g.size(), values.size());
333 A.row_indices = row_indices;
334 A.column_indices = column_indices;
335 A.values = values;
336 return A;
337}
338
360template<class real_type>
361cusp::coo_matrix<int, real_type, cusp::host_memory> interpolation(
362 const thrust::host_vector<real_type>& x,
363 const thrust::host_vector<real_type>& y,
364 const aRealTopology2d<real_type>& g,
365 dg::bc bcx = dg::NEU, dg::bc bcy = dg::NEU,
366 std::string method = "dg")
367{
368 assert( x.size() == y.size());
369 cusp::array1d<real_type, cusp::host_memory> values;
370 cusp::array1d<int, cusp::host_memory> row_indices;
371 cusp::array1d<int, cusp::host_memory> column_indices;
372 if( method == "dg")
373 {
374 std::vector<real_type> gauss_nodesx = g.dltx().abscissas();
375 std::vector<real_type> gauss_nodesy = g.dlty().abscissas();
376 dg::Operator<real_type> forwardx( g.dltx().forward());
377 dg::Operator<real_type> forwardy( g.dlty().forward());
378
379
380 for( int i=0; i<(int)x.size(); i++)
381 {
382 real_type X = x[i], Y = y[i];
383 bool negative=false;
384 g.shift( negative,X,Y, bcx, bcy);
385
386 //determine which cell (x,y) lies in
387 real_type xnn = (X-g.x0())/g.hx();
388 real_type ynn = (Y-g.y0())/g.hy();
389 unsigned nn = (unsigned)floor(xnn);
390 unsigned mm = (unsigned)floor(ynn);
391 //determine normalized coordinates
392 real_type xn = 2.*xnn - (real_type)(2*nn+1);
393 real_type yn = 2.*ynn - (real_type)(2*mm+1);
394 //interval correction
395 if (nn==g.Nx()) {
396 nn-=1;
397 xn =1.;
398 }
399 if (mm==g.Ny()) {
400 mm-=1;
401 yn =1.;
402 }
403 //Test if the point is a Gauss point since then no interpolation is needed
404 int idxX =-1, idxY = -1;
405 for( unsigned k=0; k<g.nx(); k++)
406 {
407 if( fabs( xn - gauss_nodesx[k]) < 1e-14)
408 idxX = nn*g.nx() + k; //determine which grid column it is
409 }
410 for( unsigned k=0; k<g.ny(); k++)
411 {
412 if( fabs( yn - gauss_nodesy[k]) < 1e-14)
413 idxY = mm*g.ny() + k; //determine grid line
414 }
415 if( idxX < 0 && idxY < 0 ) //there is no corresponding point
416 {
417 //evaluate 2d Legendre polynomials at (xn, yn)...
418 std::vector<real_type> px = detail::coefficients( xn, g.nx()),
419 py = detail::coefficients( yn, g.ny());
420 std::vector<real_type> pxF(g.nx(),0), pyF(g.ny(), 0);
421 for( unsigned l=0; l<g.nx(); l++)
422 for( unsigned k=0; k<g.nx(); k++)
423 pxF[l]+= px[k]*forwardx(k,l);
424 for( unsigned l=0; l<g.ny(); l++)
425 for( unsigned k=0; k<g.ny(); k++)
426 pyF[l]+= py[k]*forwardy(k,l);
427 //these are the matrix coefficients with which to multiply
428 for(unsigned k=0; k<g.ny(); k++)
429 for( unsigned l=0; l<g.nx(); l++)
430 {
431 row_indices.push_back( i);
432 column_indices.push_back( ((mm*g.ny()+k)*g.Nx()+nn)*g.nx() + l);
433 real_type pxy = pyF[k]*pxF[l];
434 if( !negative)
435 values.push_back( pxy);
436 else
437 values.push_back( -pxy);
438 }
439 }
440 else if ( idxX < 0 && idxY >=0) //there is a corresponding line
441 {
442 std::vector<real_type> px = detail::coefficients( xn, g.nx());
443 std::vector<real_type> pxF(g.nx(),0);
444 for( unsigned l=0; l<g.nx(); l++)
445 for( unsigned k=0; k<g.nx(); k++)
446 pxF[l]+= px[k]*forwardx(k,l);
447 for( unsigned l=0; l<g.nx(); l++)
448 {
449 row_indices.push_back( i);
450 column_indices.push_back( ((idxY)*g.Nx() + nn)*g.nx() + l);
451 if( !negative)
452 values.push_back( pxF[l]);
453 else
454 values.push_back(-pxF[l]);
455
456 }
457 }
458 else if ( idxX >= 0 && idxY < 0) //there is a corresponding column
459 {
460 std::vector<real_type> py = detail::coefficients( yn, g.ny());
461 std::vector<real_type> pyF(g.ny(),0);
462 for( unsigned l=0; l<g.ny(); l++)
463 for( unsigned k=0; k<g.ny(); k++)
464 pyF[l]+= py[k]*forwardy(k,l);
465 for( unsigned k=0; k<g.ny(); k++)
466 {
467 row_indices.push_back(i);
468 column_indices.push_back((mm*g.ny()+k)*g.Nx()*g.nx() + idxX);
469 if( !negative)
470 values.push_back( pyF[k]);
471 else
472 values.push_back(-pyF[k]);
473
474 }
475 }
476 else //the point already exists
477 {
478 row_indices.push_back(i);
479 column_indices.push_back(idxY*g.Nx()*g.nx() + idxX);
480 if( !negative)
481 values.push_back( 1.);
482 else
483 values.push_back(-1.);
484 }
485
486 }
487 }
488 else
489 {
490 unsigned points_per_line = 1;
491 if( method == "nearest")
492 points_per_line = 1;
493 else if( method == "linear")
494 points_per_line = 2;
495 else if( method == "cubic")
496 points_per_line = 4;
497 else
498 throw std::runtime_error( "Interpolation method "+method+" not recognized!\n");
499 RealGrid1d<real_type> gx(g.x0(), g.x1(), g.nx(), g.Nx(), bcx);
500 RealGrid1d<real_type> gy(g.y0(), g.y1(), g.ny(), g.Ny(), bcy);
501 thrust::host_vector<real_type> absX = dg::create::abscissas( gx);
502 thrust::host_vector<real_type> absY = dg::create::abscissas( gy);
503
504 for( unsigned i=0; i<x.size(); i++)
505 {
506 real_type X = x[i], Y = y[i];
507 bool negative = false;
508 g.shift( negative, X, Y, bcx, bcy);
509
510 thrust::host_vector<unsigned> colsX, colsY;
511 std::vector<real_type> xs = detail::choose_1d_abscissas( X,
512 points_per_line, gx, absX, colsX);
513 std::vector<real_type> ys = detail::choose_1d_abscissas( Y,
514 points_per_line, gy, absY, colsY);
515
516 //evaluate 2d Legendre polynomials at (xn, yn)...
517 std::vector<real_type> pxy( points_per_line*points_per_line);
518 std::vector<real_type> px = detail::lagrange( X, xs),
519 py = detail::lagrange( Y, ys);
520 // note: px , py may have size != points_per_line at boundary
521 for(unsigned k=0; k<py.size(); k++)
522 for( unsigned l=0; l<px.size(); l++)
523 pxy[k*px.size()+l]= py[k]*px[l];
524 for( unsigned k=0; k<py.size(); k++)
525 for( unsigned l=0; l<px.size(); l++)
526 {
527 if( fabs(pxy[k*px.size() +l]) > 1e-14)
528 {
529 row_indices.push_back( i);
530 column_indices.push_back( (colsY[k])*g.nx()*g.Nx() +
531 colsX[l]);
532 values.push_back( negative ? - pxy[k*px.size()+l]
533 : pxy[k*px.size()+l]);
534 }
535 }
536 }
537 }
538 cusp::coo_matrix<int, real_type, cusp::host_memory> A( x.size(),
539 g.size(), values.size());
540 A.row_indices = row_indices;
541 A.column_indices = column_indices;
542 A.values = values;
543
544 return A;
545}
546
547
548
572template<class real_type>
573cusp::coo_matrix<int, real_type, cusp::host_memory> interpolation(
574 const thrust::host_vector<real_type>& x,
575 const thrust::host_vector<real_type>& y,
576 const thrust::host_vector<real_type>& z,
577 const aRealTopology3d<real_type>& g,
578 dg::bc bcx = dg::NEU, dg::bc bcy = dg::NEU, dg::bc bcz = dg::PER,
579 std::string method = "dg")
580{
581 assert( x.size() == y.size());
582 assert( y.size() == z.size());
583 cusp::array1d<real_type, cusp::host_memory> values;
584 cusp::array1d<int, cusp::host_memory> row_indices;
585 cusp::array1d<int, cusp::host_memory> column_indices;
586
587 if( method == "dg")
588 {
589 std::vector<real_type> gauss_nodesx = g.dltx().abscissas();
590 std::vector<real_type> gauss_nodesy = g.dlty().abscissas();
591 std::vector<real_type> gauss_nodesz = g.dltz().abscissas();
592 dg::Operator<real_type> forwardx( g.dltx().forward());
593 dg::Operator<real_type> forwardy( g.dlty().forward());
594 dg::Operator<real_type> forwardz( g.dltz().forward());
595 for( int i=0; i<(int)x.size(); i++)
596 {
597 real_type X = x[i], Y = y[i], Z = z[i];
598 bool negative = false;
599 g.shift( negative,X,Y,Z, bcx, bcy, bcz);
600
601 //determine which cell (x,y,z) lies in
602 real_type xnn = (X-g.x0())/g.hx();
603 real_type ynn = (Y-g.y0())/g.hy();
604 real_type znn = (Z-g.z0())/g.hz();
605 unsigned nn = (unsigned)floor(xnn);
606 unsigned mm = (unsigned)floor(ynn);
607 unsigned ll = (unsigned)floor(znn);
608 //determine normalized coordinates
609 real_type xn = 2.*xnn - (real_type)(2*nn+1);
610 real_type yn = 2.*ynn - (real_type)(2*mm+1);
611 real_type zn = 2.*znn - (real_type)(2*ll+1);
612 //interval correction
613 if (nn==g.Nx()) {
614 nn-=1;
615 xn =1.;
616 }
617 if (mm==g.Ny()) {
618 mm-=1;
619 yn =1.;
620 }
621 if (ll==g.Nz()) {
622 ll-=1;
623 zn =1.;
624 }
625 //Test if the point is a Gauss point since then no interpolation is needed
626 int idxX =-1, idxY = -1, idxZ = -1;
627 for( unsigned k=0; k<g.nx(); k++)
628 {
629 if( fabs( xn - gauss_nodesx[k]) < 1e-14)
630 idxX = nn*g.nx() + k; //determine which grid column it is
631 }
632 for( unsigned k=0; k<g.ny(); k++)
633 {
634 if( fabs( yn - gauss_nodesy[k]) < 1e-14)
635 idxY = mm*g.ny() + k; //determine grid line
636 }
637 for( unsigned k=0; k<g.nz(); k++)
638 {
639 if( fabs( zn - gauss_nodesz[k]) < 1e-14)
640 idxZ = ll*g.nz() + k; //determine grid line
641 }
642 if( idxX >= 0 && idxY >= 0 && idxZ >= 0) //the point already exists
643 {
644 row_indices.push_back(i);
645 column_indices.push_back((idxZ*g.Ny()*g.ny()+idxY)*g.Nx()*g.nx() + idxX);
646 if( !negative)
647 values.push_back( 1.);
648 else
649 values.push_back(-1.);
650 }
651 else if ( idxX < 0 && idxY >=0 && idxZ >= 0)
652 {
653 std::vector<real_type> px = detail::coefficients( xn, g.nx());
654 std::vector<real_type> pxF(g.nx(),0);
655 for( unsigned l=0; l<g.nx(); l++)
656 for( unsigned k=0; k<g.nx(); k++)
657 pxF[l]+= px[k]*forwardx(k,l);
658 for( unsigned l=0; l<g.nx(); l++)
659 {
660 row_indices.push_back( i);
661 column_indices.push_back( (idxZ*g.Ny()*g.ny() +
662 idxY)*g.Nx()*g.nx() + nn*g.nx() + l);
663 if( !negative)
664 values.push_back( pxF[l]);
665 else
666 values.push_back(-pxF[l]);
667 }
668 }
669 else if ( idxX >= 0 && idxY < 0 && idxZ >= 0)
670 {
671 std::vector<real_type> py = detail::coefficients( yn, g.ny());
672 std::vector<real_type> pyF(g.ny(),0);
673 for( unsigned l=0; l<g.ny(); l++)
674 for( unsigned k=0; k<g.ny(); k++)
675 pyF[l]+= py[k]*forwardy(k,l);
676 for( unsigned k=0; k<g.ny(); k++)
677 {
678 row_indices.push_back(i);
679 column_indices.push_back(((idxZ*g.Ny()+mm)*g.ny()+k)*g.Nx()*g.nx() + idxX);
680 if(!negative)
681 values.push_back( pyF[k]);
682 else
683 values.push_back(-pyF[k]);
684 }
685 }
686 else
687 {
688 //evaluate 3d Legendre polynomials at (xn, yn, zn)...
689 std::vector<real_type> px = detail::coefficients( xn, g.nx()),
690 py = detail::coefficients( yn, g.ny()),
691 pz = detail::coefficients( zn, g.nz());
692 std::vector<real_type> pxF(g.nx(),0), pyF(g.ny(), 0), pzF( g.nz(), 0);
693 for( unsigned l=0; l<g.nx(); l++)
694 for( unsigned k=0; k<g.nx(); k++)
695 pxF[l]+= px[k]*forwardx(k,l);
696 for( unsigned l=0; l<g.ny(); l++)
697 for( unsigned k=0; k<g.ny(); k++)
698 pyF[l]+= py[k]*forwardy(k,l);
699 for( unsigned l=0; l<g.nz(); l++)
700 for( unsigned k=0; k<g.nz(); k++)
701 pzF[l]+= pz[k]*forwardz(k,l);
702 //these are the matrix coefficients with which to multiply
703 for( unsigned s=0; s<g.nz(); s++)
704 for( unsigned k=0; k<g.ny(); k++)
705 for( unsigned l=0; l<g.nx(); l++)
706 {
707 row_indices.push_back( i);
708 column_indices.push_back(
709 ((((ll*g.nz()+s)*g.Ny()+mm)*g.ny()+k)*g.Nx()+nn)*g.nx()+l);
710 real_type pxyz = pzF[s]*pyF[k]*pxF[l];
711 if( !negative)
712 values.push_back( pxyz);
713 else
714 values.push_back(-pxyz);
715 }
716 }
717 }
718 }
719 else
720 {
721 unsigned points_per_line = 1;
722 if( method == "nearest")
723 points_per_line = 1;
724 else if( method == "linear")
725 points_per_line = 2;
726 else if( method == "cubic")
727 points_per_line = 4;
728 else
729 throw std::runtime_error( "Interpolation method "+method+" not recognized!\n");
730 RealGrid1d<real_type> gx(g.x0(), g.x1(), g.nx(), g.Nx(), bcx);
731 RealGrid1d<real_type> gy(g.y0(), g.y1(), g.ny(), g.Ny(), bcy);
732 RealGrid1d<real_type> gz(g.z0(), g.z1(), g.nz(), g.Nz(), bcz);
733 thrust::host_vector<real_type> absX = dg::create::abscissas( gx);
734 thrust::host_vector<real_type> absY = dg::create::abscissas( gy);
735 thrust::host_vector<real_type> absZ = dg::create::abscissas( gz);
736 for( unsigned i=0; i<x.size(); i++)
737 {
738 real_type X = x[i], Y = y[i], Z = z[i];
739 bool negative = false;
740 g.shift( negative, X, Y, Z, bcx, bcy, bcz);
741
742 thrust::host_vector<unsigned> colsX, colsY, colsZ;
743 std::vector<real_type> xs = detail::choose_1d_abscissas( X,
744 points_per_line, gx, absX, colsX);
745 std::vector<real_type> ys = detail::choose_1d_abscissas( Y,
746 points_per_line, gy, absY, colsY);
747 std::vector<real_type> zs = detail::choose_1d_abscissas( Z,
748 points_per_line, gz, absZ, colsZ);
749
750 //evaluate 3d Legendre polynomials at (xn, yn, zn)...
751 std::vector<real_type> pxyz( points_per_line*points_per_line
752 *points_per_line);
753 std::vector<real_type> px = detail::lagrange( X, xs),
754 py = detail::lagrange( Y, ys),
755 pz = detail::lagrange( Z, zs);
756 // note: px, py, pz may have size != points_per_line at boundary
757 for( unsigned m=0; m<pz.size(); m++)
758 for( unsigned k=0; k<py.size(); k++)
759 for( unsigned l=0; l<px.size(); l++)
760 pxyz[(m*py.size()+k)*px.size()+l]= pz[m]*py[k]*px[l];
761 for( unsigned m=0; m<pz.size(); m++)
762 for( unsigned k=0; k<py.size(); k++)
763 for( unsigned l=0; l<px.size(); l++)
764 {
765 if( fabs(pxyz[(m*py.size()+k)*px.size() +l]) > 1e-14)
766 {
767 row_indices.push_back( i);
768 column_indices.push_back( ((colsZ[m])*g.ny()*g.Ny() +
769 colsY[k])*g.nx()*g.Nx() + colsX[l]);
770 values.push_back( negative ?
771 -pxyz[(m*py.size()+k)*px.size()+l]
772 : pxyz[(m*py.size()+k)*px.size()+l] );
773 }
774 }
775 }
776 }
777 cusp::coo_matrix<int, real_type, cusp::host_memory> A( x.size(), g.size(),
778 values.size());
779 A.row_indices = row_indices;
780 A.column_indices = column_indices;
781 A.values = values;
782
783 return A;
784}
803template<class real_type>
804cusp::coo_matrix<int, real_type, cusp::host_memory> interpolation( const RealGrid1d<real_type>& g_new, const RealGrid1d<real_type>& g_old, std::string method = "dg")
805{
806 //assert both grids are on the same box
807 assert( g_new.x0() >= g_old.x0());
808 assert( g_new.x1() <= g_old.x1());
809 thrust::host_vector<real_type> pointsX = dg::evaluate( dg::cooX1d, g_new);
810 return interpolation( pointsX, g_old, g_old.bcx(), method);
811
812}
814template<class real_type>
815cusp::coo_matrix<int, real_type, cusp::host_memory> interpolation( const aRealTopology2d<real_type>& g_new, const aRealTopology2d<real_type>& g_old, std::string method = "dg")
816{
817 //assert both grids are on the same box
818 assert( g_new.x0() >= g_old.x0());
819 assert( g_new.x1() <= g_old.x1());
820 assert( g_new.y0() >= g_old.y0());
821 assert( g_new.y1() <= g_old.y1());
822 thrust::host_vector<real_type> pointsX = dg::evaluate( dg::cooX2d, g_new);
823
824 thrust::host_vector<real_type> pointsY = dg::evaluate( dg::cooY2d, g_new);
825 return interpolation( pointsX, pointsY, g_old, g_old.bcx(), g_old.bcy(), method);
826
827}
828
830template<class real_type>
831cusp::coo_matrix<int, real_type, cusp::host_memory> interpolation( const aRealTopology3d<real_type>& g_new, const aRealTopology3d<real_type>& g_old, std::string method = "dg")
832{
833 //assert both grids are on the same box
834 assert( g_new.x0() >= g_old.x0());
835 assert( g_new.x1() <= g_old.x1());
836 assert( g_new.y0() >= g_old.y0());
837 assert( g_new.y1() <= g_old.y1());
838 assert( g_new.z0() >= g_old.z0());
839 assert( g_new.z1() <= g_old.z1());
840 thrust::host_vector<real_type> pointsX = dg::evaluate( dg::cooX3d, g_new);
841 thrust::host_vector<real_type> pointsY = dg::evaluate( dg::cooY3d, g_new);
842 thrust::host_vector<real_type> pointsZ = dg::evaluate( dg::cooZ3d, g_new);
843 return interpolation( pointsX, pointsY, pointsZ, g_old, g_old.bcx(), g_old.bcy(), g_old.bcz(), method);
844
845}
847template<class real_type>
848cusp::coo_matrix<int, real_type, cusp::host_memory> interpolation( const aRealTopology3d<real_type>& g_new, const aRealTopology2d<real_type>& g_old, std::string method = "dg")
849{
850 //assert both grids are on the same box
851 assert( g_new.x0() >= g_old.x0());
852 assert( g_new.x1() <= g_old.x1());
853 assert( g_new.y0() >= g_old.y0());
854 assert( g_new.y1() <= g_old.y1());
855 thrust::host_vector<real_type> pointsX = dg::evaluate( dg::cooX3d, g_new);
856 thrust::host_vector<real_type> pointsY = dg::evaluate( dg::cooY3d, g_new);
857 return interpolation( pointsX, pointsY, g_old, g_old.bcx(), g_old.bcy(), method);
858
859}
861
862
863}//namespace create
864
865
883template<class real_type>
884real_type interpolate(
885 dg::space sp,
886 const thrust::host_vector<real_type>& v,
887 real_type x,
888 const RealGrid1d<real_type>& g,
889 dg::bc bcx = dg::NEU)
890{
891 assert( v.size() == g.size());
892 bool negative = false;
893 g.shift( negative, x, bcx);
894
895 //determine which cell (x) lies in
896
897 real_type xnn = (x-g.x0())/g.h();
898 unsigned n = (unsigned)floor(xnn);
899 //determine normalized coordinates
900
901 real_type xn = 2.*xnn - (real_type)(2*n+1);
902 //interval correction
903 if (n==g.N()) {
904 n-=1;
905 xn = 1.;
906 }
907 //evaluate 1d Legendre polynomials at (xn)...
908 std::vector<real_type> px = create::detail::coefficients( xn, g.n());
909 if( sp == dg::xspace)
910 {
911 dg::Operator<real_type> forward( g.dlt().forward());
912 std::vector<real_type> pxF(g.n(),0);
913 for( unsigned l=0; l<g.n(); l++)
914 for( unsigned k=0; k<g.n(); k++)
915 pxF[l]+= px[k]*forward(k,l);
916 for( unsigned k=0; k<g.n(); k++)
917 px[k] = pxF[k];
918 }
919 //these are the matrix coefficients with which to multiply
920 unsigned cols = (n)*g.n();
921 //multiply x
922 real_type value = 0;
923 for( unsigned j=0; j<g.n(); j++)
924 {
925 if(negative)
926 value -= v[cols + j]*px[j];
927 else
928 value += v[cols + j]*px[j];
929 }
930 return value;
931}
932
952template<class real_type>
953real_type interpolate(
954 dg::space sp,
955 const thrust::host_vector<real_type>& v,
956 real_type x, real_type y,
957 const aRealTopology2d<real_type>& g,
958 dg::bc bcx = dg::NEU, dg::bc bcy = dg::NEU )
959{
960 assert( v.size() == g.size());
961 bool negative = false;
962 g.shift( negative, x,y, bcx, bcy);
963
964 //determine which cell (x,y) lies in
965
966 real_type xnn = (x-g.x0())/g.hx();
967 real_type ynn = (y-g.y0())/g.hy();
968 unsigned n = (unsigned)floor(xnn);
969 unsigned m = (unsigned)floor(ynn);
970 //determine normalized coordinates
971
972 real_type xn = 2.*xnn - (real_type)(2*n+1);
973 real_type yn = 2.*ynn - (real_type)(2*m+1);
974 //interval correction
975 if (n==g.Nx()) {
976 n-=1;
977 xn = 1.;
978 }
979 if (m==g.Ny()) {
980 m-=1;
981 yn =1.;
982 }
983 //evaluate 2d Legendre polynomials at (xn, yn)...
984 std::vector<real_type> px = create::detail::coefficients( xn, g.nx()),
985 py = create::detail::coefficients( yn, g.ny());
986 if( sp == dg::xspace)
987 {
988 dg::Operator<real_type> forwardx( g.dltx().forward());
989 dg::Operator<real_type> forwardy( g.dlty().forward());
990 std::vector<real_type> pxF(g.nx(),0), pyF(g.ny(), 0);
991 for( unsigned l=0; l<g.nx(); l++)
992 for( unsigned k=0; k<g.nx(); k++)
993 pxF[l]+= px[k]*forwardx(k,l);
994 for( unsigned l=0; l<g.ny(); l++)
995 for( unsigned k=0; k<g.ny(); k++)
996 pyF[l]+= py[k]*forwardy(k,l);
997 px = pxF, py = pyF;
998 }
999 //these are the matrix coefficients with which to multiply
1000 unsigned cols = (m)*g.Nx()*g.ny()*g.nx() + (n)*g.nx();
1001 //multiply x
1002 real_type value = 0;
1003 for( unsigned i=0; i<g.ny(); i++)
1004 for( unsigned j=0; j<g.nx(); j++)
1005 {
1006 if(negative)
1007 value -= v[cols + i*g.Nx()*g.nx() + j]*px[j]*py[i];
1008 else
1009 value += v[cols + i*g.Nx()*g.nx() + j]*px[j]*py[i];
1010 }
1011 return value;
1012}
1013
1014} //namespace dg
evaluation.h
Function discretization routines.
functions.h
Some utility functions for the dg::evaluate routines.
grid.h
base topology classes
dg::cooX1d
static DG_DEVICE double cooX1d(double x)
Definition: functions.h:38
dg::cooY2d
static DG_DEVICE double cooY2d(double x, double y)
Definition: functions.h:45
dg::cooZ3d
static DG_DEVICE double cooZ3d(double x, double y, double z)
Definition: functions.h:49
dg::cooY3d
static DG_DEVICE double cooY3d(double x, double y, double z)
Definition: functions.h:47
dg::cooX2d
static DG_DEVICE double cooX2d(double x, double y)
Definition: functions.h:40
dg::cooX3d
static DG_DEVICE double cooX3d(double x, double y, double z)
Definition: functions.h:42
dg::bc
bc
Switch between boundary conditions.
Definition: enums.h:15
dg::space
space
Space of DG coefficients.
Definition: enums.h:164
dg::coo3d::z
@ z
z direction
dg::PER
@ PER
periodic boundaries
Definition: enums.h:16
dg::NEU
@ NEU
Neumann on both boundaries.
Definition: enums.h:20
dg::xspace
@ xspace
Configuration space "nodal values".
Definition: enums.h:166
dg::forward
@ forward
forward derivative (cell to the right and current cell)
Definition: enums.h:98
dg::coo2d::y
@ y
y direction
dg::coo2d::x
@ x
x direction
dg::evaluate
thrust::host_vector< real_type > evaluate(UnaryOp f, const RealGrid1d< real_type > &g)
Evaluate a 1d function on grid coordinates.
Definition: evaluation.h:67
dg::create::interpolation
cusp::coo_matrix< int, real_type, cusp::host_memory > interpolation(const thrust::host_vector< real_type > &x, const RealGrid1d< real_type > &g, dg::bc bcx=dg::NEU, std::string method="dg")
Create interpolation matrix.
Definition: interpolation.h:254
dg::interpolate
real_type interpolate(dg::space sp, const thrust::host_vector< real_type > &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.
Definition: interpolation.h:884
dg
This is the namespace for all functions and classes defined and used by the discontinuous Galerkin li...
operator_tensor.h
dg::RealGrid1d
1D grid
Definition: grid.h:80
dg::RealGrid1d::h
real_type h() const
cell size
Definition: grid.h:129
dg::RealGrid1d::bcx
bc bcx() const
boundary conditions
Definition: grid.h:147
dg::RealGrid1d::n
unsigned n() const
number of polynomial coefficients
Definition: grid.h:141
dg::RealGrid1d::dlt
const DLT< real_type > & dlt() const
the discrete legendre transformation
Definition: grid.h:197
dg::RealGrid1d::shift
void shift(bool &negative, real_type &x) const
Shift any point coordinate to a corresponding grid coordinate according to the boundary condition.
Definition: grid.h:214
dg::RealGrid1d::size
unsigned size() const
n()*N() (Total number of grid points)
Definition: grid.h:191
dg::RealGrid1d::x1
real_type x1() const
right boundary
Definition: grid.h:117
dg::RealGrid1d::N
unsigned N() const
number of cells
Definition: grid.h:135
dg::RealGrid1d::x0
real_type x0() const
left boundary
Definition: grid.h:111
dg::aRealTopology2d
An abstract base class for two-dimensional grids.
Definition: grid.h:277
dg::aRealTopology2d::dltx
const DLT< real_type > & dltx() const
discrete legendre trafo
Definition: grid.h:372
dg::aRealTopology2d::shift
void shift(bool &negative, real_type &x, real_type &y) const
Shift any point coordinate to a corresponding grid coordinate according to the boundary condition.
Definition: grid.h:454
dg::aRealTopology2d::x0
real_type x0() const
Left boundary in x.
Definition: grid.h:288
dg::aRealTopology2d::dlty
const DLT< real_type > & dlty() const
discrete legendre transformation in y
Definition: grid.h:374
dg::aRealTopology2d::y1
real_type y1() const
Right boundary in y.
Definition: grid.h:306
dg::aRealTopology2d::hy
real_type hy() const
cell size in y
Definition: grid.h:330
dg::aRealTopology2d::ny
unsigned ny() const
number of polynomial coefficients in y
Definition: grid.h:340
dg::aRealTopology2d::bcx
bc bcx() const
boundary conditions in x
Definition: grid.h:358
dg::aRealTopology2d::size
unsigned size() const
The total number of points.
Definition: grid.h:424
dg::aRealTopology2d::hx
real_type hx() const
cell size in x
Definition: grid.h:324
dg::aRealTopology2d::y0
real_type y0() const
left boundary in y
Definition: grid.h:300
dg::aRealTopology2d::Nx
unsigned Nx() const
number of cells in x
Definition: grid.h:346
dg::aRealTopology2d::bcy
bc bcy() const
boundary conditions in y
Definition: grid.h:364
dg::aRealTopology2d::x1
real_type x1() const
Right boundary in x.
Definition: grid.h:294
dg::aRealTopology2d::nx
unsigned nx() const
number of polynomial coefficients in x
Definition: grid.h:338
dg::aRealTopology2d::Ny
unsigned Ny() const
number of cells in y
Definition: grid.h:352
dg::aRealTopology3d
An abstract base class for three-dimensional grids.
Definition: grid.h:523
dg::aRealTopology3d::y0
real_type y0() const
left boundary in y
Definition: grid.h:547
dg::aRealTopology3d::size
unsigned size() const
The total number of points.
Definition: grid.h:670
dg::aRealTopology3d::dltz
const DLT< real_type > & dltz() const
discrete legendre transformation in z
Definition: grid.h:658
dg::aRealTopology3d::dltx
const DLT< real_type > & dltx() const
discrete legendre transformation in x
Definition: grid.h:654
dg::aRealTopology3d::nz
unsigned nz() const
number of polynomial coefficients in z
Definition: grid.h:616
dg::aRealTopology3d::bcz
bc bcz() const
boundary conditions in z
Definition: grid.h:652
dg::aRealTopology3d::Nx
unsigned Nx() const
number of points in x
Definition: grid.h:622
dg::aRealTopology3d::ny
unsigned ny() const
number of polynomial coefficients in y
Definition: grid.h:614
dg::aRealTopology3d::shift
void shift(bool &negative, real_type &x, real_type &y, real_type &z) const
Shift any point coordinate to a corresponding grid coordinate according to the boundary condition.
Definition: grid.h:710
dg::aRealTopology3d::dlty
const DLT< real_type > & dlty() const
discrete legendre transformation in y
Definition: grid.h:656
dg::aRealTopology3d::x0
real_type x0() const
left boundary in x
Definition: grid.h:534
dg::aRealTopology3d::hy
real_type hy() const
cell size in y
Definition: grid.h:598
dg::aRealTopology3d::y1
real_type y1() const
right boundary in y
Definition: grid.h:553
dg::aRealTopology3d::hx
real_type hx() const
cell size in x
Definition: grid.h:592
dg::aRealTopology3d::z0
real_type z0() const
left boundary in z
Definition: grid.h:560
dg::aRealTopology3d::x1
real_type x1() const
right boundary in x
Definition: grid.h:540
dg::aRealTopology3d::Ny
unsigned Ny() const
number of points in y
Definition: grid.h:628
dg::aRealTopology3d::bcy
bc bcy() const
boundary conditions in y
Definition: grid.h:646
dg::aRealTopology3d::bcx
bc bcx() const
boundary conditions in x
Definition: grid.h:640
dg::aRealTopology3d::hz
real_type hz() const
cell size in z
Definition: grid.h:604
dg::aRealTopology3d::z1
real_type z1() const
right boundary in z
Definition: grid.h:566
dg::aRealTopology3d::Nz
unsigned Nz() const
number of points in z
Definition: grid.h:634
dg::aRealTopology3d::nx
unsigned nx() const
number of polynomial coefficients in x
Definition: grid.h:612
typedefs.h
Useful typedefs of commonly used types.
xspacelib.h
utility functions