LoAdSG
StencilSGpp.h
1//
2// Created by to35jepo on 5/3/24.
3//
4
5#ifndef RUN_STENCILSGPP_H
6#define RUN_STENCILSGPP_H
7#include "Stencil.h"
8
9#include "../../../SGpp/base/src/sgpp_base.hpp"
10
11
12enum QuadratureType {
13 SGadaptive,
14 SGregular,
15 MonteCarlo,
16 MonteCarlo_Hierarchichal,
17 MonteCarloFixed
18};
19
20
21template <typename F>
22class StencilSGpp: public StencilTemplate {
23public:
24
25
26 int level;
27 bool adaptive;
28 int numberofgridpoints=0;
29 CellDimension cellDimension;
30 QuadratureType quadrature_type;
31 int gridpoints=2147483647;
32 std::unique_ptr<sgpp::base::Grid> grid_MC;
33 int maxlevel=0;
34 int getMinimumGridPointsOfACell(){
35 int localvalue = gridpoints;
36 int global_min;
37 MPI_Reduce(&localvalue, &global_min, 1, MPI_INT, MPI_MIN, 0, MPI_COMM_WORLD);
38 return global_min;
39 };
40
41 int mcPoints=1;
42
43 StencilSGpp(AdaptiveSparseGrid& grid_, const F &varCoeff_, int level_=1,QuadratureType quadrature_type_=SGadaptive) : StencilTemplate(grid_), varCoeff(varCoeff_), level(level_), quadrature_type(quadrature_type_){
44
45 if(quadrature_type==MonteCarloFixed) {
46 grid_MC = std::unique_ptr<sgpp::base::Grid>(sgpp::base::Grid::createLinearGrid(DimensionSparseGrid));
47 grid_MC->getGenerator().regular(level);
48
49
50 }
51 //gridpoints = grid_MC->getSize();
52 mcPoints=grid.getDOFS();
53
54 };
55
56
57
58
59 inline double returnValue(const IndexDimension &Index, const MultiDimCompass &mc){cout << "StencilSGPP returnValue not implmented" << endl; exit(1);}
60
61 inline double integration(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
62 double value = 0.0;
63
64
65 if(quadrature_type == SGadaptive) {
66 value = integration_adaptive(cell,dirP,dirQ);
67 } else if(quadrature_type == SGregular){
68 value = integration_regular(cell,dirP,dirQ);
69 } else if(quadrature_type == MonteCarlo){
70 value = integration_MC(cell,dirP,dirQ);
71 } else if(quadrature_type == MonteCarlo_Hierarchichal){
72 value = integration_MC_AD(cell,dirP,dirQ);
73 }else if(quadrature_type == MonteCarloFixed) {
74 value = integration_MC_fixedCell(cell, dirP, dirQ);
75 }
76
77
78 return value;
79 }
80 inline double integration_regular(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ);
81
82 inline double integrationLocallyAdaptive(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ);
83
84
85 inline double integration_adaptive(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ);
86
87 inline double integration_gausslegendre(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ);
88
89
90
91 inline double integration_MC(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ);
92
93
94 inline double integration_MC_fixedCell(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ);
95
96 inline double integration_MC_fixedCell2(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ);
97
98
99
100 inline double integration_MC_fixed(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ);
101
102
103
104 inline double integration_MC_AD(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ);
105protected:
106
107 F varCoeff;
108
109
110
111 //in dieser Klasse soll die Funktion erstellt werden, die dann über das Einheitsquadrat integriert werden soll. Stichwort: Substitution
112 class FunctionClass{
113
114
115 public:
116 FunctionClass (CellDimension& cell_, CellIndexDirection& dirP_, CellIndexDirection &dirQ_,F varCoeff_ ): cell(cell_),dirP(dirP_),dirQ(dirQ_),varCoeff(varCoeff_){
117 // transformation of variable xi gives
118 // J : xi -> (bi-ai)xi+ai
119 // det(J')=(b1-a1)*...*(bd-ad)
120 det_J_prime=1.0;
121 for(int i=0; i<DimensionSparseGrid; i++){
122 double h = cell.getRightBoundary(i)-cell.getLeftBoundary(i);
123 det_J_prime*=h;
124 meshwidth[i]=h;
125
126 }
127
128 }
129 double evaluation(int dim, double* x, void* clientdata){
130 // transform variable
131 // J : xi -> (bi-ai)xi+ai
132 // det(J')=(b1-a1)*...*(bd-ad)
133
134 for(int i=0; i<DimensionSparseGrid; i++){
135 x[i]=meshwidth[i]*x[i]+cell.getLeftBoundary(i);
136 }
137
138
139
140
141
142
143 double value = 1.0;
144
145 for(int i = 0; i < DimensionSparseGrid; ++i )
146 {
147
148 // === Muliply by the first base function ===
149 if(dirP.getDirection(i)==Left){
150 // p(x) = (x-b)/(a-b)=(b-x)/h
151 value *= (cell.getRightBoundary(i)-x[i]) / (meshwidth[i]);
152
153 }else{
154 // p(x) = (x-a)/(b-a)
155 value *= (x[i] - cell.getLeftBoundary(i)) / (meshwidth[i]);
156 }
157
158 // === Muliply by the second base function ===
159 if(dirQ.getDirection(i)==Left){
160 // p(x) = (x-b)/(a-b)=(b-x)/h
161 value *= (cell.getRightBoundary(i)-x[i]) / (meshwidth[i]);
162
163
164 }else{
165 // p(x) = (x-a)/(b-a)
166 value *= (x[i] - cell.getLeftBoundary(i)) / (meshwidth[i]);
167 }
168
169
170 }
171
172
173 value *= varCoeff(x);
174
175
176 return value;
177 };
178
179
180
181
182 double detJPrime(){return det_J_prime;};
183
184 static double evaluationWrapper(int dim, double* x, void* clientdata) {
185 return static_cast<FunctionClass*>(clientdata)->evaluation(dim, x, nullptr);
186 }
187 private:
188 double det_J_prime = 1.0;
189 double meshwidth[DimensionSparseGrid];
190 CellDimension cell;
191 CellIndexDirection dirP;
192 CellIndexDirection dirQ;
193 F varCoeff;
194 };
195
196
197
198};
199
200template<typename F>
201double StencilSGpp<F>::integration_MC_AD(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
202
203
204 FunctionClass function(cell,dirP,dirQ,varCoeff);
205
206 int levels[DimensionSparseGrid];
207 int maxlevel=0;
208 for(int d=0; d<DimensionSparseGrid;d++){
209 levels[d]=grid.getMaxDepth(d)-(cell.getDepth(d));
210 maxlevel<levels[d]?maxlevel=levels[d]:maxlevel;
211 }
212
213
214
215
216
217 int dim = DimensionSparseGrid;
218 std::unique_ptr<sgpp::base::Grid> grid_sgpp(sgpp::base::Grid::createLinearGrid(dim));
219 sgpp::base::GridStorage& gridStorage = grid_sgpp->getStorage();
220
221 // create regular grid, level
222 grid_sgpp->getGenerator().regular(1);
223
224
225 sgpp::base::DataVector alpha(gridStorage.getSize());
226 sgpp::base::DataVector refined(gridStorage.getSize());
227
228
229 std::vector<size_t> addedPoints;
230
231
232 for (int step = 0; step < maxlevel; step++) {
233
234 int alpha_c=0;
235 for (size_t i = 0; i < gridStorage.getSize(); i++) {
236 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
237 alpha[i]=1.0;
238
239
240 // do not allow refinement if max level in one direction is reached
241 for(size_t k=0; k< DimensionSparseGrid; k++){
242 if(gp.getLevel(k)> levels[k]){
243 alpha[i]=0.0;
244 break;
245 }
246
247 }
248
249
250
251 // ensure that grid point can only be refined once
252 if(refined[i]>1.0) alpha[i]=0.0;
253
254
255 if(alpha[i]>0.0) {
256 alpha_c++;
257
258 // this means that this grid point can't be refined in the next refinement step
259 refined[i] = 2.0;
260 }
261
262
263 }
264
265 if (alpha_c<1)break;
266
267 sgpp::base::SurplusRefinementFunctor functor(alpha,alpha_c, 0.5);
268 grid_sgpp->getGenerator().refine(functor, &addedPoints);
269
270
271 alpha.resize(gridStorage.getSize());
272 refined.resize(gridStorage.getSize());
273
274 if(addedPoints.size()<1)break;
275 addedPoints.clear();
276 }
277
278 double p[DimensionSparseGrid];
279 for (size_t i = 0; i < gridStorage.getSize(); i++) {
280 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
281 for(int j=0; j<DimensionSparseGrid;j++){
282 p[j] = gp.getStandardCoordinate(j);
283 }
284 alpha[i] = function.evaluation(DimensionSparseGrid, p, NULL);
285 }
286
287 if(numberofgridpoints<gridStorage.getSize()){
288 numberofgridpoints=gridStorage.getSize();
289 cellDimension = cell;
290 }
291
292
293
294 std::unique_ptr<sgpp::base::OperationHierarchisation>(sgpp::op_factory::createOperationHierarchisation(*grid_sgpp))
295 ->doHierarchisation(alpha);
296
297 sgpp::base::OperationQuadratureMC opMC(*grid_sgpp,gridStorage.getSize());
298
299
300
301
302 //cout << " do mc" << endl;
303 double res = opMC.doQuadrature(alpha);
304 //double res = opMC.doQuadratureFunc(FunctionClass::evaluationWrapper, &function);
305
306
307
308 res *= function.detJPrime();
309
310 return res;
311}
312
313template<typename F>
314double StencilSGpp<F>::integration_MC_fixed(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
315
316 FunctionClass function(cell,dirP,dirQ,varCoeff);
317
318 int dim = DimensionSparseGrid;
319
320
321
322 std::unique_ptr<sgpp::base::Grid> grid2(sgpp::base::Grid::createLinearBoundaryGrid(DimensionSparseGrid));
323
324
325
326
327 sgpp::base::OperationQuadratureMC opMC(*grid2, grid.getDOFS());
328
329 double res = opMC.doQuadratureFunc(FunctionClass::evaluationWrapper, &function);
330 res *= function.detJPrime();
331
332 return res;
333}
334
335template<typename F>
336double StencilSGpp<F>::integration_MC_fixedCell(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
337
338
339 FunctionClass function(cell,dirP,dirQ,varCoeff);
340
341 int l=0;
342 for(int d=0; d<DimensionSparseGrid;d++){
343 l+=cell.getDepth(d)-1;
344 }
345
346 int vol_cell= POW2(unsigned (l));
347
348 if(vol_cell>grid_MC->getSize()) mcPoints=1;
349 else mcPoints = grid_MC->getSize()/vol_cell;
350
351
352
353 sgpp::base::OperationQuadratureMC opMC(*grid_MC,mcPoints);
354
355 double res = opMC.doQuadratureFunc(FunctionClass::evaluationWrapper, &function);
356
357
358
359 res *= function.detJPrime();
360
361 return res;
362}
363
364template<typename F>
365double StencilSGpp<F>::integration_MC_fixedCell2(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
366
367
368 FunctionClass function(cell,dirP,dirQ,varCoeff);
369
370
371
372 mcPoints = grid_MC->getSize();
373 sgpp::base::OperationQuadratureMC opMC(*grid_MC,mcPoints);
374
375 double res = opMC.doQuadratureFunc(FunctionClass::evaluationWrapper, &function);
376
377
378
379 res *= function.detJPrime();
380
381 return res;
382}
383
384
385
386template<typename F>
387double StencilSGpp<F>::integration_MC(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
388
389
390 FunctionClass function(cell,dirP,dirQ,varCoeff);
391
392 int levels[DimensionSparseGrid];
393 int maxlevel=0;
394 for(int d=0; d<DimensionSparseGrid;d++){
395 levels[d]=grid.getMaxDepth(d)-(cell.getDepth(d));
396 maxlevel<levels[d]?maxlevel=levels[d]:maxlevel;
397 }
398
399
400
401
402
403 int dim = DimensionSparseGrid;
404 std::unique_ptr<sgpp::base::Grid> grid_sgpp(sgpp::base::Grid::createLinearGrid(dim));
405 sgpp::base::GridStorage& gridStorage = grid_sgpp->getStorage();
406
407 // create regular grid, level
408 grid_sgpp->getGenerator().regular(1);
409
410
411 sgpp::base::DataVector alpha(gridStorage.getSize());
412 sgpp::base::DataVector refined(gridStorage.getSize());
413
414
415 std::vector<size_t> addedPoints;
416
417
418 for (int step = 0; step < maxlevel; step++) {
419
420 int alpha_c=0;
421 for (size_t i = 0; i < gridStorage.getSize(); i++) {
422 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
423 alpha[i]=1.0;
424
425
426 // do not allow refinement if max level in one direction is reached
427 for(size_t k=0; k< DimensionSparseGrid; k++){
428 if(gp.getLevel(k)> levels[k]){
429 alpha[i]=0.0;
430 break;
431 }
432
433 }
434
435
436
437 // ensure that grid point can only be refined once
438 if(refined[i]>1.0) alpha[i]=0.0;
439
440
441 if(alpha[i]>0.0) {
442 alpha_c++;
443
444 // this means that this grid point can't be refined in the next refinement step
445 refined[i] = 2.0;
446 }
447
448
449 }
450
451 if (alpha_c<1)break;
452
453 sgpp::base::SurplusRefinementFunctor functor(alpha,alpha_c, 0.5);
454 grid_sgpp->getGenerator().refine(functor, &addedPoints);
455
456
457 alpha.resize(gridStorage.getSize());
458 refined.resize(gridStorage.getSize());
459
460 if(addedPoints.size()<1)break;
461 addedPoints.clear();
462 }
463
464/* double p[DimensionSparseGrid];
465 for (size_t i = 0; i < gridStorage.getSize(); i++) {
466 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
467 for(int j=0; j<DimensionSparseGrid;j++){
468 p[j] = gp.getStandardCoordinate(j);
469 }
470 alpha[i] = function.evaluation(DimensionSparseGrid, p, NULL);
471 }
472
473 if(numberofgridpoints<gridStorage.getSize()){
474 numberofgridpoints=gridStorage.getSize();
475 cellDimension = cell;
476 }
477
478
479
480 std::unique_ptr<sgpp::base::OperationHierarchisation>(sgpp::op_factory::createOperationHierarchisation(*grid_sgpp))
481 ->doHierarchisation(alpha);*/
482
483 sgpp::base::OperationQuadratureMC opMC(*grid_sgpp,gridStorage.getSize());
484/* std::unique_ptr<sgpp::base::OperationHierarchisation>(sgpp::op_factory::createOperationHierarchisation(*grid_sgpp))
485 ->doHierarchisation(alpha);*/
486
487
488
489 //cout << " do mc" << endl;
490 //double res = opMC.doQuadrature(alpha);
491 double res = opMC.doQuadratureFunc(FunctionClass::evaluationWrapper, &function);
492
493
494
495 res *= function.detJPrime();
496
497 return res;
498}
499
500template<typename F>
501double
502StencilSGpp<F>::integration_gausslegendre(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
503
504
505 FunctionClass function(cell,dirP,dirQ,varCoeff);
506
507 int levels[DimensionSparseGrid];
508 int maxlevel=0;
509 for(int d=0; d<DimensionSparseGrid;d++){
510 levels[d]=grid.getMaxDepth(d)-(cell.getDepth(d));
511 maxlevel<levels[d]?maxlevel=levels[d]:maxlevel;
512 }
513
514
515
516
517
518 int dim = DimensionSparseGrid;
519 std::unique_ptr<sgpp::base::Grid> grid_sgpp(sgpp::base::Grid::createLinearGrid(dim));
520 sgpp::base::GridStorage& gridStorage = grid_sgpp->getStorage();
521
522 // create regular grid, level
523 grid_sgpp->getGenerator().regular(1);
524
525
526 sgpp::base::DataVector alpha(gridStorage.getSize());
527 sgpp::base::DataVector refined(gridStorage.getSize());
528
529
530 std::vector<size_t> addedPoints;
531
532
533 for (int step = 0; step < maxlevel; step++) {
534
535 int alpha_c=0;
536 for (size_t i = 0; i < gridStorage.getSize(); i++) {
537 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
538 alpha[i]=1.0;
539
540
541 // do not allow refinement if max level in one direction is reached
542 for(size_t k=0; k< DimensionSparseGrid; k++){
543 if(gp.getLevel(k)> levels[k]){
544 alpha[i]=0.0;
545 break;
546 }
547
548 }
549
550
551
552 // ensure that grid point can only be refined once
553 if(refined[i]>1.0) alpha[i]=0.0;
554
555
556 if(alpha[i]>0.0) {
557 alpha_c++;
558
559 // this means that this grid point can't be refined in the next refinement step
560 refined[i] = 2.0;
561 }
562
563
564 }
565
566 if (alpha_c<1)break;
567
568 sgpp::base::SurplusRefinementFunctor functor(alpha,alpha_c, 0.5);
569 grid_sgpp->getGenerator().refine(functor, &addedPoints);
570
571
572 alpha.resize(gridStorage.getSize());
573 refined.resize(gridStorage.getSize());
574
575 if(addedPoints.size()<1)break;
576 addedPoints.clear();
577 }
578
579
580 double p[DimensionSparseGrid];
581 for (size_t i = 0; i < gridStorage.getSize(); i++) {
582 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
583 for(int j=0; j<DimensionSparseGrid;j++){
584 p[j] = gp.getStandardCoordinate(j);
585 }
586 alpha[i] = function.evaluation(DimensionSparseGrid, p, NULL);
587 }
588
589
590 double res;
591
592
593
594 res*= function.detJPrime();
595
596 return res;
597}
598
599template<typename F>
600double StencilSGpp<F>::integration_adaptive(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
601
602
603 FunctionClass function(cell,dirP,dirQ,varCoeff);
604
605 int levels[DimensionSparseGrid];
606 int maxlevel=0;
607 for(int d=0; d<DimensionSparseGrid;d++){
608 levels[d]=grid.getMaxDepth(d)-(cell.getDepth(d));
609 maxlevel<levels[d]?maxlevel=levels[d]:maxlevel;
610 }
611
612
613
614
615
616 int dim = DimensionSparseGrid;
617 std::unique_ptr<sgpp::base::Grid> grid_sgpp(sgpp::base::Grid::createLinearGrid(dim));
618 sgpp::base::GridStorage& gridStorage = grid_sgpp->getStorage();
619
620 // create regular grid, level
621 grid_sgpp->getGenerator().regular(1);
622
623
624 sgpp::base::DataVector alpha(gridStorage.getSize());
625 sgpp::base::DataVector refined(gridStorage.getSize());
626
627
628 std::vector<size_t> addedPoints;
629
630
631 for (int step = 0; step < maxlevel; step++) {
632
633 int alpha_c=0;
634 for (size_t i = 0; i < gridStorage.getSize(); i++) {
635 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
636 alpha[i]=1.0;
637
638
639 // do not allow refinement if max level in one direction is reached
640 for(size_t k=0; k< DimensionSparseGrid; k++){
641 if(gp.getLevel(k)> levels[k]){
642 alpha[i]=0.0;
643 break;
644 }
645
646 }
647
648
649
650 // ensure that grid point can only be refined once
651 if(refined[i]>1.0) alpha[i]=0.0;
652
653
654 if(alpha[i]>0.0) {
655 alpha_c++;
656
657 // this means that this grid point can't be refined in the next refinement step
658 refined[i] = 2.0;
659 }
660
661
662 }
663
664 if (alpha_c<1)break;
665
666 sgpp::base::SurplusRefinementFunctor functor(alpha,alpha_c, 0.5);
667 grid_sgpp->getGenerator().refine(functor, &addedPoints);
668
669
670 alpha.resize(gridStorage.getSize());
671 refined.resize(gridStorage.getSize());
672
673 if(addedPoints.size()<1)break;
674 addedPoints.clear();
675 }
676
677
678 double p[DimensionSparseGrid];
679 for (size_t i = 0; i < gridStorage.getSize(); i++) {
680 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
681 for(int j=0; j<DimensionSparseGrid;j++){
682 p[j] = gp.getStandardCoordinate(j);
683 }
684 alpha[i] = function.evaluation(DimensionSparseGrid, p, NULL);
685 }
686
687 if(numberofgridpoints<gridStorage.getSize()){
688 numberofgridpoints=gridStorage.getSize();
689 cellDimension = cell;
690 }
691
692
693 // Prints Grid
694/*
695 sgpp::base::GridPrinter gridPrinter(*grid_sgpp);
696 string filename = "sgpp_grid.gnu";
697 gridPrinter.printSparseGrid(alpha,filename,false);
698*/
699
700 std::unique_ptr<sgpp::base::OperationHierarchisation>(sgpp::op_factory::createOperationHierarchisation(*grid_sgpp))
701 ->doHierarchisation(alpha);
702
703 // direct quadrature
704 std::unique_ptr<sgpp::base::OperationQuadrature> opQ(
705 sgpp::op_factory::createOperationQuadrature(*grid_sgpp));
706 double res = opQ->doQuadrature(alpha);
707
708 res *= function.detJPrime();
709
710 return res;
711}
712
713template<typename F>
714double
715StencilSGpp<F>::integrationLocallyAdaptive(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
716
717
718 FunctionClass function(cell,dirP,dirQ,varCoeff);
719
720
721 int dim = DimensionSparseGrid;
722 std::unique_ptr<sgpp::base::Grid> grid_sgpp(sgpp::base::Grid::createLinearBoundaryGrid(dim));
723 sgpp::base::GridStorage& gridStorage = grid_sgpp->getStorage();
724
725 // create regular grid, level
726 grid_sgpp->getGenerator().regular(level);
727
728
729 sgpp::base::DataVector alpha(gridStorage.getSize());
730 sgpp::base::DataVector refined(gridStorage.getSize());
731
732 double p[DimensionSparseGrid];
733 sgpp::base::DataVector funEvals(gridStorage.getSize());
734 for (size_t i = 0; i < gridStorage.getSize(); i++) {
735 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
736 for(size_t k=0; k<DimensionSparseGrid; k++){
737 p[k]=gp.getStandardCoordinate(k);
738 }
739 funEvals[i] = function.evaluation(DimensionSparseGrid,p,NULL);
740 refined[i]=0.0;
741
742 }
743
744 std::vector<size_t> addedPoints;
745
746 alpha.copyFrom(funEvals);
747 sgpp::op_factory::createOperationHierarchisation(*grid_sgpp)->doHierarchisation(alpha);
748
749
750
751
752
753 int number_refinements=10;
754 for (int step = 0; step < number_refinements-1; step++) {
755
756
757 int alpha_c=0;
758 for (size_t i = 0; i < gridStorage.getSize(); i++) {
759 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
760
761 for(size_t k=0; k< DimensionSparseGrid; k++){
762 if(gp.getLevel(k)>= grid.getMaxDepth(int(k))-cell.getDepth(int(k))){
763
764 alpha[i]=0.0;
765 break;
766 }
767
768 }
769
770 if(refined[i]>1.0) alpha[i]=0.0;
771
772
773 if(alpha[i]>0.0) {
774 alpha_c++;
775 refined[i] = 2.0;
776 }
777
778
779
780
781
782 }
783
784 cout << " alpha > 0 " << alpha_c << endl;
785
786
787
788 sgpp::base::SurplusRefinementFunctor functor(alpha,alpha_c, 0.0);
789 grid_sgpp->getGenerator().refine(functor, &addedPoints);
790
791 cout << " added points " << addedPoints.size() << endl;
792
793
794
795
796
797
798
799 alpha.resize(gridStorage.getSize());
800 funEvals.resize(gridStorage.getSize());
801 refined.resize(gridStorage.getSize());
802
803
804 for (size_t i = 0; i < addedPoints.size(); i++) {
805 size_t seq = addedPoints[i];
806 sgpp::base::GridPoint& gp = gridStorage.getPoint(seq);
807 for(size_t k=0; k<DimensionSparseGrid; k++){
808 p[k]=gp.getStandardCoordinate(k);
809 }
810 funEvals[seq] = function.evaluation(DimensionSparseGrid,p,NULL);
811 }
812
813
814
815 alpha.copyFrom(funEvals);
816 sgpp::op_factory::createOperationHierarchisation(*grid_sgpp)->doHierarchisation(alpha);
817
818 int max=0;
819 int sum=0;
820 for(size_t i=0; i<gridStorage.getSize(); i++){
821 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
822 for(int d=0; d<DimensionSparseGrid;d++){
823 sum+=gp.getLevel(d);
824 }
825 if(max<sum)max=sum;
826 sum=0;
827 }
828
829 cout << "max " << max << endl;
830
831
832
833
834
835 if(addedPoints.size()<1)break;
836 addedPoints.clear();
837 }
838
839
840 cout << " final points " << gridStorage.getSize() << endl;
841
842
843
844
845 // direct quadrature
846 std::unique_ptr<sgpp::base::OperationQuadrature> opQ(
847 sgpp::op_factory::createOperationQuadrature(*grid_sgpp));
848 double res = opQ->doQuadrature(alpha);
849
850 res *= function.detJPrime();
851
852 return res;
853}
854
855template<typename F>
856double StencilSGpp<F>::integration_regular(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
857
858
859 FunctionClass function(cell,dirP,dirQ,varCoeff);
860
861 int dim = DimensionSparseGrid;
862 std::unique_ptr<sgpp::base::Grid> grid_sgpp(sgpp::base::Grid::createLinearBoundaryGrid(dim));
863 sgpp::base::GridStorage& gridStorage = grid_sgpp->getStorage();
864
865 // create regular grid, level 3
866
867 grid_sgpp->getGenerator().regular(level);
868
869
870 //ab hier soll SGpp anfangen, function.evaluation() ist dann die Funktion die von SGpp integriert werden soll
871
872 sgpp::base::DataVector alpha(gridStorage.getSize());
873 double p[DimensionSparseGrid];
874 for (size_t i = 0; i < gridStorage.getSize(); i++) {
875 sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
876 for(int j=0; j<DimensionSparseGrid;j++){
877 p[j] = gp.getStandardCoordinate(j);
878 }
879 alpha[i] = function.evaluation(DimensionSparseGrid, p, NULL);
880 }
881
882 std::unique_ptr<sgpp::base::OperationHierarchisation>(sgpp::op_factory::createOperationHierarchisation(*grid_sgpp))
883 ->doHierarchisation(alpha);
884
885 // direct quadrature
886 std::unique_ptr<sgpp::base::OperationQuadrature> opQ(
887 sgpp::op_factory::createOperationQuadrature(*grid_sgpp));
888 double res = opQ->doQuadrature(alpha);
889
890 res *= function.detJPrime();
891
892
893
894
895 return res;
896}
897
898
899#endif //RUN_STENCILSGPP_H
AdaptiveSparseGrid
Definition sparseGrid.h:277
MultiDimCompass
Definition index.h:356