LoAdSG
Stencil.h
1//
2// Created by to35jepo on 3/7/23.
3//
4
5#ifndef SGRUN_STENCIL_H
6#define SGRUN_STENCIL_H
7
8#include "../indices/index.h"
9
10
11
12
13#include "InterfaceIntegration/interfaceMatrices.h"
14#include "InterfaceIntegration/constantIntegrators.h"
15#include "InterfaceIntegration/interatorBasisFunction.h"
16#include "InterfaceIntegration/HelmholtzIntegrator/helmIntegrator.h"
17
18#include "../cells/celldimension.h"
19#include "../cells/iterator_neu.h"
20#include "../applications/norms.h"
21#include "../iterator/RectangularIterator.h"
22#include "PoissonStencil.h"
23
24
25#include <omp.h>
26
27#include <sstream>
28
29
30
31
32
33
34
35
36template <class Stencilclass>
37double getLocalStencilBoundary(IndexDimension Index, Depth &T, VectorSparseG &u, Stencilclass stencilclass) {
38
39
40 stencilclass.initialize(T);
41
42 double val = 0.0;
43 IndexDimension NextTest;
44
45 for (MultiDimCompass mc; mc.goon(); ++mc) {
46 NextTest = Index.nextThree_boundary(&mc, T.returnTiefen());
47
48 if(mc.goon()) {
49 double value = 0.0;
50 value = stencilclass.returnValue(Index, mc);
51 val = val + value * u.getValue(NextTest);
52 }
53 }
54 return val;
55
56}
57
58
59class CellIterator{
60public:
61 CellIterator(MultiDimCompass mc, Depth T, IndexDimension center){
62 int val=1;
63
64 // Berechne Anzahl der Zellen, die MC berühren
65 for(int d=0; d<DimensionSparseGrid; d++){
66 if(mc.getRichtung(d)==Mitte&&center.isNotAtBoundary(d))val=2*val;
67 if(mc.getRichtung(d)==Links&&center.isAtLeftBoundary(d)) val=0.0;
68 if(mc.getRichtung(d)==Rechts&&center.isAtRightBoundary(d)) val=0.0;
69 }
70 maxShift = val;
71
72 Indices = new IndexDimension[maxShift];
73 recursive(center,0,mc,T);
74 }
75 ~CellIterator(){
76 delete[] Indices;
77 }
78
79 void PrintIndices(){
80 for(int j=0; j<maxShift;j++){
81 Indices[j].PrintCoord();
82 }
83 }
84
85 int maxShift;
86
87 bool goon(){
88 return (shiftNumber<maxShift);
89 }
90
91 void operator++(){
92 shiftNumber++;
93 }
94
95 IndexDimension getCell(){
96 return Indices[shiftNumber];
97 }
98
99 void start(){
100 shiftNumber=0;
101 }
102
103private:
104
105 IndexDimension* Indices;
106 int k = 0;
107 int shiftNumber=0;
108
109 void recursive(IndexDimension I, int d, MultiDimCompass mc, Depth T){
110 if(d == DimensionSparseGrid) {
111 Indices[k]=I;
112 k++;
113 }else {
114 int D = d+1;
115 if (mc.getRichtung(d) == Mitte) {
116 // rechts
117 recursive(I, D, mc,T);
118 //links
119 if(!I.isAtLeftBoundary(d)){
120 I = I.nextLeft(d, T.at(d));
121 recursive(I, D, mc, T);
122 }
123 }
124 if (mc.getRichtung(d) == Links && !I.isAtLeftBoundary(d)) {
125 I = I.nextLeft(d, T.at(d));
126 recursive(I, D, mc,T);
127 }
128 if (mc.getRichtung(d) == Rechts &&!I.isAtRightBoundary(d)) {
129 recursive(I, D, mc,T);
130 }
131 }
132 }
133};
134
135
136
137
138class StencilTemplate{
139
140 template<class StencilA, class StencilB>
141 friend class ADD;
142
143
144 public:
145 // ghost functions and members
146 virtual void resetActiveWorkers(){};
147 virtual void receiveApplySendOnActiveWorkers(){};
148 virtual int getNumberProcesses(){return 0;};
149 virtual int size(){return 0;};
150 std::vector<int> active_worker;
151 virtual void applyLocalStiffnessMatricesDepth(VectorSparseG& u,VectorSparseG& z,Depth& depth){};
152 virtual double applyLocalStiffnessMatricesFixedDepthIndex_onNode(VectorSparseG& input, VectorSparseG& output, Depth& D, IndexDimension& Index){return 0.0;};
153 virtual int getNodeForDepth(Depth& T){return 0;};
154 void applyLocalStiffnessMatricesOnIndices_onNode(VectorSparseG& input, VectorSparseG& output, Depth& D, std::vector<IndexDimension>& Indices){};
155
156
157 StencilTemplate(AdaptiveSparseGrid& grid_): grid(grid_){};
158
159 virtual inline void initialize(Depth &T_) {
160 T = T_;
161
162 for (int d = 0; d < DimensionSparseGrid; d++) {
163 int exp = T_.at(d);
164 if (exp == 0) meshwidth[d] = 1.0;
165 else {
166
167 int value2 = 1 << exp;
168 exponent[d] = double(value2);
169 meshwidth[d] = 1.0 / double(value2);
170 }
171 }
172
173 };
174
175 inline void copyInitialization(StencilTemplate& stencil){
176 T = stencil.T;
177
178 for (int d = 0; d < DimensionSparseGrid; d++) {
179 meshwidth[d]=stencil.meshwidth[d];
180 exponent[d]=stencil.exponent[d];
181
182 }
183 }
184
185
186
187 virtual inline void applyStencilOnCell(CellDimension& cell,VectorSparseG& input, VectorSparseG& output){
188
189 for(CellIndexIterator outerIter(&cell); outerIter.goon(); ++outerIter) {
190
191
192 IndexDimension p = outerIter.getIndex();
193
194 unsigned long kp;
195 bool occP = output.getSparseGrid()->occupied(kp, p);
196 CellIndexDirection dirP = outerIter.getCellIndexDirection();
197
198
199 if (occP) {
200
201 for (CellIndexIterator innerIter(outerIter); innerIter.goon(); ++innerIter) {
202
203
204 IndexDimension q = innerIter.getIndex();
205 if (q.isNotAtBoundary()) {
206
207 CellIndexDirection dirQ = innerIter.getCellIndexDirection();
208
209 unsigned long kq;
210
211 bool occQ = output.getSparseGrid()->occupied(kq, q);
212 if (occQ) {
213
214 double val = integration(cell, dirP, dirQ);
215 if (p == q) {
216 val = val * input.getValue(kp);
217
218
219 output.addToValue(kq, val);
220 } else {
221 double val_p = val * input.getValue(kp);
222
223 output.addToValue(kq, val_p);
224
225 double val_q = val * input.getValue(kq);
226
227 output.addToValue(kp, val_q);
228 }
229 }
230
231 }
232
233 }
234
235 }
236 }
237
238
239
240
241
242
243
244
245
246 }
247
248 virtual inline void applyStencilOnCell_BinaryIterator(CellDimension& cell,VectorSparseG& input, VectorSparseG& output){
249
250
251
252 int iterend = PowerOfTwo<DimensionSparseGrid>::value;
253
254 CellIndexIterator outerIter(&cell);
255 for(int iteri=0; iteri<iterend ; iteri++) {
256 IndexDimension p = outerIter.getIndexByBinary(iteri);
257 unsigned long kp;
258 bool occP = output.getSparseGrid()->occupied(kp, p);
259 CellIndexDirection dirP = outerIter.getCellIndexDirection();
260 if (occP) {
261 for (int iterj = iteri; iterj < iterend; iterj++) {
262 IndexDimension q = outerIter.getIndexByBinary(iterj);
263 if (q.isNotAtBoundary()) {
264 CellIndexDirection dirQ = outerIter.getCellIndexDirection();
265 unsigned long kq;
266
267 bool occQ = output.getSparseGrid()->occupied(kq, q);
268 if (occQ) {
269
270 double val = integration(cell, dirP, dirQ);
271 if (p == q) {
272 val = val * input.getValue(kp);
273
274 //#pragma omp critical
275 output.addToValue(kq, val);
276 } else {
277 double val_p = val * input.getValue(kp);
278 //#pragma omp critical
279 output.addToValue(kq, val_p);
280
281 double val_q = val * input.getValue(kq);
282 //#pragma omp critical
283 output.addToValue(kp, val_q);
284 }
285 }
286
287 }
288 }
289 }
290 }
291
292 }
293
294 virtual inline void applyStencilOnCell_nosymmetry(CellDimension& cell,VectorSparseG& input, VectorSparseG& output) {
295
296
297 for (CellIndexIterator outerIter(&cell); outerIter.goon(); ++outerIter) {
298
299
300 IndexDimension p = outerIter.getIndex();
301
302 unsigned long kp;
303 bool occP = output.getSparseGrid()->occupied(kp, p);
304 CellIndexDirection dirP = outerIter.getCellIndexDirection();
305
306
307 if (occP) {
308
309
310
311 for (CellIndexIterator innerIter(&cell); innerIter.goon(); ++innerIter) {
312
313 IndexDimension q = innerIter.getIndex();
314 if(q.isNotAtBoundary()){
315 CellIndexDirection dirQ = innerIter.getCellIndexDirection();
316
317 unsigned long kq;
318
319 bool occQ= output.getSparseGrid()->occupied(kq,q);
320 if(occQ){
321
322
323 double val = integration(cell, dirP, dirQ);
324 val =val*input.getValue(kp);
325 output.addToValue(kq, val);
326
327
328 }
329 }
330 } }
331 }
332 }
333
334 virtual inline void applyStencilOnCellMPI(CellDimension& cell,VectorSparseG& input, VectorSparseG& output) {
335
336
337 int iterend = POW2(DimensionSparseGrid);
338
339
340 {
341 CellIndexIterator outerIter(&cell);
342
343 for (int iteri = 0; iteri < iterend; iteri++) {
344
345
346 IndexDimension p = outerIter.getIndex();
347
348 unsigned long kp;
349 bool occP = output.getSparseGrid()->occupied(kp, p);
350 CellIndexDirection dirP = outerIter.getCellIndexDirection();
351
352
353 if (occP) {
354 for (CellIndexIterator innerIter(outerIter); innerIter.goon(); ++innerIter) {
355 IndexDimension q = innerIter.getIndex();
356 if (q.isNotAtBoundary()) {
357 CellIndexDirection dirQ = innerIter.getCellIndexDirection();
358
359 unsigned long kq;
360
361 bool occQ = output.getSparseGrid()->occupied(kq, q);
362 if(occQ){
363 double val = integration(cell, dirP, dirQ);
364 if(p==q){
365 val =val*input.getValue(kp);
366
367 output.addToValue(kq, val);
368 }else{
369 double val_p =val*input.getValue(kp);
370
371 output.addToValue(kq, val_p);
372
373 double val_q =val*input.getValue(kq);
374
375 output.addToValue(kp, val_q);
376 }
377 }
378
379 }
380
381 }
382
383
384 }
385 ++outerIter;
386 }
387 }
388 }
389
390
391 virtual inline void applyStencilOnCell_MPI_OMP(CellDimension& cell,VectorSparseG& input, VectorSparseG& output) {
392
393
394 int iterend = POW2(DimensionSparseGrid);
395 CellIndexIterator outerIter(&cell);
396
397 int entries = iterend*(iterend+1)/2;
398 CellIndexDirection pCellIndexDirection[2][entries];
399 bool calcIntgral[entries];
400 bool diagonalEntry[entries];
401 unsigned long occupied[2][entries];
402 double values[2][entries];
403 int j=0;
404
405
406 for (CellIndexIterator outerIter2(&cell); outerIter2.goon(); ++outerIter2){
407
408 IndexDimension p = outerIter2.getIndex();
409 unsigned long kp;
410 bool occP = output.getSparseGrid()->occupied(kp, p);
411 CellIndexDirection dirP = outerIter2.getCellIndexDirection();
412
413
414 for (CellIndexIterator innerIter(outerIter2); innerIter.goon(); ++innerIter) {
415 IndexDimension q = innerIter.getIndex();
416
417 if (q.isNotAtBoundary()) {
418
419 CellIndexDirection dirQ = innerIter.getCellIndexDirection();
420
421 unsigned long kq;
422
423 bool occQ = output.getSparseGrid()->occupied(kq, q);
424 if (occQ && occP) {
425 calcIntgral[j] = true;
426 pCellIndexDirection[0][j] = dirP;
427 pCellIndexDirection[1][j] = dirQ;
428 occupied[0][j]=kp;
429 occupied[1][j]=kq;
430 if (p == q) {
431 diagonalEntry[j] = true;
432
433
434 } else {
435 diagonalEntry[j] = false;
436 }
437 } else {
438 calcIntgral[j] = false;
439 }
440
441 } else {
442 calcIntgral[j] = false;
443 }
444 j++;
445 }
446 }
447
448
449#pragma omp parallel for schedule(dynamic)
450 for (int i = 0; i < entries; i++) {
451 if (calcIntgral[i]) {
452 double val = integration(cell, pCellIndexDirection[0][i], pCellIndexDirection[1][i]);
453 unsigned long kp = occupied[0][i];
454 unsigned long kq = occupied[1][i];
455 if (diagonalEntry[i]) {
456 val = val * input.getValue(kp);
457 values[0][i] = val;
458
459 } else {
460 double val_p = val * input.getValue(kp);
461 values[1][i] = val_p;
462
463
464 double val_q = val * input.getValue(kq);
465
466 values[0][i] = val_q;
467
468 }
469
470 }
471 }
472
473
474
475
476
477 for (int i = 0; i < entries; i++) {
478 if (calcIntgral[i]) {
479 unsigned long kp = occupied[0][i];
480 unsigned long kq = occupied[1][i];
481 if (diagonalEntry[i]) {
482
483 double val =values[0][i];
484
485 output.addToValue(kq, val);
486 } else {
487 //double val_p = val * input.getValue(kp);
488 double val_p =values[1][i];
489
490
491 output.addToValue(kq, val_p);
492
493 double val_q= values[0][i];
494 output.addToValue(kp, val_q);
495 }
496
497 }
498 }
499
500 }
501
502
503 virtual inline void applyStencilOnCell_MPI_OMP_nosymmetry(CellDimension& cell,VectorSparseG& input, VectorSparseG& output) {
504
505
506
507 int iterend = POW2(DimensionSparseGrid);
508
509 CellIndexIterator outerIter(&cell);
510
511 int entries = iterend*iterend;
512 CellIndexDirection pCellIndexDirection[2][entries];
513 bool calcIntgral[entries];
514
515 unsigned long occupied[2][entries];
516 double values[entries];
517 int j=0;
518
519
520 for (CellIndexIterator outerIter2(&cell); outerIter2.goon(); ++outerIter2){
521
522 IndexDimension p = outerIter2.getIndex();
523 unsigned long kp;
524 bool occP = output.getSparseGrid()->occupied(kp, p);
525 CellIndexDirection dirP = outerIter2.getCellIndexDirection();
526
527
528 /* for (CellIndexIterator innerIter(&cell); innerIter.goon(); ++innerIter) {
529 IndexDimension q = innerIter.getIndex();
530 if (q.isNotAtBoundary()) {
531 CellIndexDirection dirQ = innerIter.getCellIndexDirection();
532
533 unsigned long kq;
534
535 bool occQ = output.getSparseGrid()->occupied(kq, q);
536 if (occQ) {
537
538
539 double val = integration(cell, dirP, dirQ);
540
541 val = val * input.getValue(kp);
542 output.addToValue(kp, val);
543
544
545 }
546
547 }
548 }*/
549 for (CellIndexIterator innerIter(&cell); innerIter.goon(); ++innerIter) {
550 IndexDimension q = innerIter.getIndex();
551
552 if (q.isNotAtBoundary() && p.isNotAtBoundary()){
553 CellIndexDirection dirQ = innerIter.getCellIndexDirection();
554
555 unsigned long kq;
556
557 bool occQ = output.getSparseGrid()->occupied(kq, q);
558 if (occQ && occP) {
559 calcIntgral[j] = true;
560 pCellIndexDirection[0][j] = dirP;
561 pCellIndexDirection[1][j] = dirQ;
562 occupied[0][j]=kp;
563 occupied[1][j]=kq;
564 } else {
565 calcIntgral[j] = false;
566 }
567
568 } else {
569 calcIntgral[j] = false;
570 }
571 j++;
572 }
573
574 }
575
576
577#pragma omp parallel for schedule(runtime)
578 for (int i = 0; i < entries; i++) {
579 if (calcIntgral[i]) {
580 double val = integration(cell, pCellIndexDirection[0][i], pCellIndexDirection[1][i]);
581 unsigned long kp = occupied[0][i];
582 //unsigned long kq = occupied[1][i];
583
584 double val_p = val * input.getValue(kp);
585
586 values[i] = val_p;
587 }
588 }
589
590
591 for (int i = 0; i < entries; i++) {
592 if (calcIntgral[i]) {
593 //unsigned long kp = occupied[0][i];
594 unsigned long kq = occupied[1][i];
595 double val_p =values[i];
596 output.addToValue(kq, val_p);
597 }
598 }
599
600/*
601 //#pragma omp parallel for schedule(runtime)
602 for (int iteri=0; iteri < iterend; iteri++) {
603 int thread_num = omp_get_thread_num();
604
605
606 IndexDimension p = outerIter.getIndex(iteri);
607
608
609
610 unsigned long kp;
611 bool occP = output.getSparseGrid()->occupied(kp, p);
612 CellIndexDirection dirP = outerIter.getCellIndexDirection(iteri);
613
614
615
616 if (occP) {
617 for (CellIndexIterator innerIter(outerIter,iteri); innerIter.goon(); ++innerIter) {
618 IndexDimension q = innerIter.getIndex();
619
620
621 if (q.isNotAtBoundary()) {
622
623 CellIndexDirection dirQ = innerIter.getCellIndexDirection();
624
625 unsigned long kq;
626
627 bool occQ = output.getSparseGrid()->occupied(kq, q);
628 if(occQ){
629 double val = integration(cell, dirP, dirQ);
630 if(p==q){
631 val =val*input.getValue(kp);
632
633
634 output.addToValue(kq, val);
635 }else{
636 double val_p =val*input.getValue(kp);
637
638
639 output.addToValue(kq, val_p);
640
641 double val_q =val*input.getValue(kq);
642
643
644 output.addToValue(kp, val_q);
645 }
646 }
647
648 }
649
650 }
651
652
653 }
654
655 ++outerIter;
656
657 }
658 output_test=output_test-output;
659 if(L_infty(output_test)>1e-20){
660
661 exit(0);}*/
662
663
664 }
665
666
667 virtual inline double returnValue(const IndexDimension &Index, const MultiDimCompass &mc){
668
669
670 double val = 0.0;
671
672
673
674 CellIterator cells(mc,T,Index);
675
676
677 for(cells.start(); cells.goon();++cells){
678 IndexDimension cell = cells.getCell();
679 if(cellInGrid(cell)) {
680
681
682 val += localValue(Index, cell, mc);
683
684
685
686
687 }
688
689 }
690
691
692
693
694
695 return val;
696
697
698
699 };
700
701
702
703 bool cellInGrid(IndexDimension cell) const{
704 IndexDimension Imin = cell;
705 IndexDimension Imax = cell;
706
707 for(int d=0; d<DimensionSparseGrid; d++){
708 Imax = Imax.nextRight(d,T.at(d));
709 }
710
711 RectangularIterator iterator(Imin,Imax,T);
712 for(;iterator.goon();++iterator){
713 IndexDimension I = iterator.getIndex();
714 unsigned long k;
715 if(I.isNotAtBoundary() && !(grid.occupied(k,I))) return false;
716 }
717
718 return true;
719 }
735 virtual inline double localValue(const IndexDimension &center ,const IndexDimension& cell, const MultiDimCompass &mc)const{
736 //cout << "value " << endl;
737 return 0.0;
738 };
739
740 virtual inline double integration(CellDimension& cell, CellIndexDirection& dirP, CellIndexDirection& dirQ){
741
742 cout << "not implemented yet! " << endl;
743 exit(1);
744 return 0.0;
745 }
746 Depth T;
747
748 virtual TypeMatrixVectorMultiplication getTypeMatrixVectorMultiplication(){return typeMatrixVectorMultiplication;};
749
750
751protected:
752
753
754
755
756 // integrate (-x+p)*(x-(p-h))
757 inline double integral_left(const IndexDimension& center,const IndexDimension& cell, int d) const{
758
759 double p=center.coordinate(d);
760 double h=meshwidth[d];
761 double x = cell.coordinate(d);
762 double val1=-(1.0/3.0)*x*x*x+0.5*(p-h)*x*x+0.5*p*x*x-p*(p-h)*x;
763
764 x = x+h;
765 double val2=-(1.0/3.0)*x*x*x+0.5*(p-h)*x*x+0.5*p*x*x-p*(p-h)*x;
766 return val2-val1;
767 }
768
769 // integrate (-x+(p+h))*(x-p)
770 inline double integral_right(const IndexDimension& center,const IndexDimension& cell, int d) const{
771 double p=center.coordinate(d);
772 double h=meshwidth[d];
773 double x = cell.coordinate(d);
774
775 double val1=-(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
776
777 x = x+h;
778 double val2=-(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
779
780 return val2-val1;
781 //return -(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
782 }
783
784
785 inline double integral_CenterR(const IndexDimension& center,const IndexDimension& cell, int d) const{
786
787
788 double p =center.coordinate(d);
789 double h =meshwidth[d];
790 double x = cell.coordinate(d);
791
792
793
794 double val1= (1.0/3.0)*x*x*x-(p+h)*x*x+(p+h)*(p+h)*x;
795
796
797 x = x+h;
798 double val2= (1.0/3.0)*x*x*x-(p+h)*x*x+(p+h)*(p+h)*x;
799
800 return val2-val1;
801
802 }
803
804
805
806 // integrate (x-(p-h))²
807 inline double integral_CenterL(const IndexDimension& center,const IndexDimension& cell, int d) const{
808 double p = center.coordinate(d);
809 double h=meshwidth[d];
810 double x = cell.coordinate(d);
811 double val1=(1.0/3.0)*x*x*x-(p-h)*x*x+(p-h)*(p-h)*x;
812
813
814 x = x+h;
815 double val2=(1.0/3.0)*x*x*x-(p-h)*x*x+(p-h)*(p-h)*x;
816
817
818 return val2-val1;
819
820 }
821
822
823 virtual inline double integral_RR(CellDimension& cell,int d){
824 double b = cell.getRightBoundary(d);
825 double a = cell.getLeftBoundary(d);
826 double h= b-a;
827
828 //integrate (x-b)² from a to b
829 double val1=(1.0/3.0)*h*h*h;
830
831 return val1;
832 }
833
834
835 virtual inline double integral_LL(CellDimension& cell,int d){
836 double b = cell.getRightBoundary(d);
837 double a = cell.getLeftBoundary(d);
838 double h= b-a;
839
840
841 //integrate (a-x)² from a to b
842 double val1=(1.0/3.0)*h*h*h;
843
844 return val1;
845 }
846
847 virtual inline double integral_LR(CellDimension& cell,int d){
848 double b = cell.getRightBoundary(d);
849 double a = cell.getLeftBoundary(d);
850 double h= b-a;
851
852 double val1=(1.0/6.0)*h*h*h;
853 return val1;
854 }
855
856
857
858
859 double meshwidth[DimensionSparseGrid];
860 double exponent[DimensionSparseGrid];
861
862
863
864 AdaptiveSparseGrid& grid;
865
866 TypeMatrixVectorMultiplication typeMatrixVectorMultiplication = StencilOnTheFly;
867
868
869};
870
871
872
873class HelmHoltz: public StencilTemplate{
874public:
875 HelmHoltz(AdaptiveSparseGrid& grid_): StencilTemplate(grid_){};
876
877
878 inline double localValue(const IndexDimension &center ,const IndexDimension& cell, const MultiDimCompass &mc)const{
879
880
881 double value = 1.0;
882 for(int d=0; d<DimensionSparseGrid; d++) {
883 if (mc.getRichtung(d) == Mitte && cell.getIndex(d) == center.getIndex(d))
884 value *= integral_CenterR(center,cell, d);
885 if (mc.getRichtung(d) == Mitte && !(cell.getIndex(d) == center.getIndex(d)))
886 value *= integral_CenterL(center,cell, d);
887 if (mc.getRichtung(d) == Links && !(cell.getIndex(d) == center.getIndex(d)))
888 value *= integral_left(center,cell, d);
889 if (mc.getRichtung(d) == Rechts && cell.getIndex(d) == center.getIndex(d))
890 value *= integral_right(center,cell, d);
891
892 value = exponent[d]*exponent[d]*value;
893 }
894
895
896
897
898 return value;}
899
900 inline double integration(CellDimension& cell, CellIndexDirection& dirP, CellIndexDirection& dirQ){
901 double value = 1.0;
902
903 for(int d=0; d<DimensionSparseGrid; d++) {
904
905
906 if (dirP.getDirection(d)==dirQ.getDirection(d)){
907 value *= integral_RR(cell, d);
908 value = exponent[d] * exponent[d] * value;
909
910 }
911 else {
912 value *= integral_LR(cell,d);
913 value = exponent[d] * exponent[d] * value;
914 }
915
916 }
917 return value;
918
919 }
920
921protected:
922
923 // integrate (-x+p)*(x-(p-h))
924 inline double integral_left(const IndexDimension& center,const IndexDimension& cell, int d) const{
925
926 double p =center.coordinate(d);
927 double h=meshwidth[d];
928 double x = cell.coordinate(d);
929 double val1=-(1.0/3.0)*x*x*x+0.5*(p-h)*x*x+0.5*p*x*x-p*(p-h)*x;
930
931 x = x+h;
932 double val2=-(1.0/3.0)*x*x*x+0.5*(p-h)*x*x+0.5*p*x*x-p*(p-h)*x;
933 return val2-val1;
934 }
935
936 // integrate (-x+(p+h))*(x-p)
937 inline double integral_right(const IndexDimension& center,const IndexDimension& cell, int d) const{
938 double p =center.coordinate(d);
939 double h=meshwidth[d];
940 double x = cell.coordinate(d);
941
942 double val1=-(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
943
944 x = x+h;
945 double val2=-(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
946
947 return val2-val1;
948 //return -(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
949 }
950
951 inline double integral_CenterR(const IndexDimension& center,const IndexDimension& cell, int d) const{
952
953 double p =center.coordinate(d);
954 double h=meshwidth[d];
955 double x = cell.coordinate(d);
956
957
958
959 double val1= (1.0/3.0)*x*x*x-(p+h)*x*x+(p+h)*(p+h)*x;
960
961
962 x = x+h;
963 double val2= (1.0/3.0)*x*x*x-(p+h)*x*x+(p+h)*(p+h)*x;
964
965 return val2-val1;
966
967 }
968
969 // integrate (x-(p-h))²
970 inline double integral_CenterL(const IndexDimension& center,const IndexDimension& cell, int d) const{
971 double p = center.coordinate(d);
972 double h=meshwidth[d];
973 double x = cell.coordinate(d);
974 double val1=(1.0/3.0)*x*x*x-(p-h)*x*x+(p-h)*(p-h)*x;
975
976
977 x = x+h;
978 double val2=(1.0/3.0)*x*x*x-(p-h)*x*x+(p-h)*(p-h)*x;
979
980
981 return val2-val1;
982
983 }
984
985
986
987
988};
989
990template <typename F>
991class StencilInterface: public StencilTemplate {
992public:
993
994 StencilInterface(AdaptiveSparseGrid& grid_, const F &varCoeff) : StencilTemplate(grid_), localStiffnessHelm(
995 varCoeff, 1e-8, 0, 8, 1.0e8) {
996 };
997
998 inline double
999 localValue(const IndexDimension &center, const IndexDimension &cell, const MultiDimCompass &mc) const {
1000
1001 BasisFunctionType u[DimensionSparseGrid];
1002 BasisFunctionType v[DimensionSparseGrid];
1003 double p_left[DimensionSparseGrid];
1004 double p_right[DimensionSparseGrid];
1005
1006 //IntegratorHelm<double (*) (const std::array<double, DimensionSparseGrid>&), DimensionSparseGrid> localStiffnessHelm(&var_coeff, 1e-10, 0, 10, 1.0e9 );
1007
1008 double value = 0.0;
1009 bool integrate = true;
1010 for (int d = 0; d < DimensionSparseGrid; d++) {
1011 p_left[d] = cell.coordinate(d);
1012 p_right[d] = cell.coordinate(d) + meshwidth[d];
1013
1014 if (mc.getRichtung(d) == Mitte && cell.getIndex(d) == center.getIndex(d)) {
1015 u[d] = leftBasis;
1016 v[d] = leftBasis;
1017 } else if (mc.getRichtung(d) == Mitte && cell.getIndex(d) != center.getIndex(d)) {
1018
1019 u[d] = rightBasis;
1020 v[d] = rightBasis;
1021 } else if (mc.getRichtung(d) == Links && cell.getIndex(d) != center.getIndex(d)) {
1022 u[d] = leftBasis;
1023 v[d] = rightBasis;
1024 } else if (mc.getRichtung(d) == Rechts && cell.getIndex(d) == center.getIndex(d)) {
1025 v[d] = leftBasis;
1026 u[d] = rightBasis;
1027 } else {
1028 integrate = false;
1029 break;
1030 }
1031
1032
1033 }
1034
1035 if (integrate) {
1036
1037 value = localStiffnessHelm.stencil_integration(p_left, p_right, u, v);
1038
1039
1040 }
1041
1042 stencil_count++;
1043 return value;
1044 };
1045
1046 static int stencil_count;
1047
1048 inline double integration(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
1049
1050 BasisFunctionType u[DimensionSparseGrid];
1051 BasisFunctionType v[DimensionSparseGrid];
1052 double p_left[DimensionSparseGrid];
1053 double p_right[DimensionSparseGrid];
1054
1055 double value = 0.0;
1056
1057 for (int d = 0; d < DimensionSparseGrid; d++) {
1058 p_left[d] = cell.getLeftBoundary(d);
1059 p_right[d] = cell.getRightBoundary(d);
1060
1061
1062 if (dirP.getDirection(d) == Left) {
1063 u[d] = leftBasis;
1064 } else {
1065 u[d] = rightBasis;
1066 }
1067
1068 if (dirQ.getDirection(d) == Left) {
1069 v[d] = leftBasis;
1070 } else {
1071 v[d] = rightBasis;
1072 }
1073
1074
1075 }
1076 value = localStiffnessHelm.stencil_integration(p_left, p_right, u, v);
1077
1078 return value;
1079 }
1080protected:
1081
1082 IntegratorHelm<const F, DimensionSparseGrid> localStiffnessHelm;
1083
1084
1085/* inline double integration(CellDimension &cell, CellIndexDirection &dirP, CellIndexDirection &dirQ) {
1086
1087 BasisFunctionType u[DimensionSparseGrid];
1088 BasisFunctionType v[DimensionSparseGrid];
1089 double p_left[DimensionSparseGrid];
1090 double p_right[DimensionSparseGrid];
1091
1092 double value = 0.0;
1093
1094 for (int d = 0; d < DimensionSparseGrid; d++) {
1095 p_left[d] = cell.getLeftBoundary(d);
1096 p_right[d] = cell.getRightBoundary(d);
1097
1098
1099 if (dirP.getDirection(d) == Left) {
1100 u[d] = leftBasis;
1101 } else {
1102 u[d] = rightBasis;
1103 }
1104
1105 if (dirQ.getDirection(d) == Left) {
1106 v[d] = leftBasis;
1107 } else {
1108 v[d] = rightBasis;
1109 }
1110
1111
1112 }
1113 value = localStiffnessHelm.stencil_integration(p_left, p_right, u, v);
1114 stencil_count++;
1115 return value;
1116 }*/
1117};
1118
1119template <class F>
1120int StencilInterface<F>::stencil_count=0;
1121
1122
1123
1127class VariableCoefficientMass: public StencilTemplate{
1128
1129public:
1130 VariableCoefficientMass(AdaptiveSparseGrid& grid_): StencilTemplate(grid_){};
1131
1132
1133
1134
1135 inline double localValue(const IndexDimension &center ,const IndexDimension& cell, const MultiDimCompass &mc)const{
1136
1137
1138 double value = 1.0;
1139 for(int d=0; d<DimensionSparseGrid; d++) {
1140 if (mc.getRichtung(d) == Mitte && cell.getIndex(d) == center.getIndex(d))
1141 value *= integral_CenterR(center,cell, d);
1142 if (mc.getRichtung(d) == Mitte && !(cell.getIndex(d) == center.getIndex(d)))
1143 value *= integral_CenterL(center,cell, d);
1144 if (mc.getRichtung(d) == Links && !(cell.getIndex(d) == center.getIndex(d)))
1145 value *= integral_left(center,cell, d);
1146 if (mc.getRichtung(d) == Rechts && cell.getIndex(d) == center.getIndex(d))
1147 value *= integral_right(center,cell, d);
1148
1149 value = exponent[d]*exponent[d]*value;
1150 }
1151
1152
1153 return value;};
1154
1155 inline double integration(CellDimension& cell, CellIndexDirection& dirP, CellIndexDirection& dirQ){
1156 double value = 1.0;
1157
1158 for(int d=0; d<DimensionSparseGrid; d++) {
1159
1160
1161 if (dirP.getDirection(d)==dirQ.getDirection(d)){
1162 if(dirP.getDirection(d)==Right) {
1163 value *= integral_LL(cell, d);
1164
1165
1166 value = exponent[d] * exponent[d] * value;
1167 }
1168 if(dirP.getDirection(d)==Left) {
1169 value *= integral_RR(cell, d);
1170
1171 value = exponent[d] * exponent[d] * value;
1172 }
1173 }
1174 else {
1175 value *= integral_LR(cell,d);
1176
1177 value = exponent[d] * exponent[d] * value;
1178 }
1179
1180
1181
1182 }
1183 return value;
1184
1185 }
1186
1187 protected:
1188
1189
1190 // integrate (-x+p)*(x-(p-h))*(1-x²)
1191 inline double integral_left(const IndexDimension& center,const IndexDimension& cell, int d) const{
1192
1193 double p =center.coordinate(d);
1194 double h=meshwidth[d];
1195
1196 double x = cell.coordinate(d);
1197 double val1=-(1.0/3.0)*x*x*x+0.5*(p-h)*x*x+0.5*p*x*x-p*(p-h)*x;
1198 double val12=-(1.0/5.0)*x*x*x*x*x+(1.0/4.0)*p*x*x*x*x+(1.0/4.0)*(p-h)*x*x*x*x-(1.0/3.0)*p*(p-h)*x*x*x;
1199
1200 val1 = val1 - val12;
1201 x = x+h;
1202 double val2=-(1.0/3.0)*x*x*x+0.5*(p-h)*x*x+0.5*p*x*x-p*(p-h)*x;
1203 double val22=-(1.0/5.0)*x*x*x*x*x+(1.0/4.0)*p*x*x*x*x+(1.0/4.0)*(p-h)*x*x*x*x-(1.0/3.0)*p*(p-h)*x*x*x;
1204
1205 val2= val2 - val22;
1206
1207 return val2-val1;
1208 }
1209
1210 // integrate (x-(p-h))²
1211 inline double integral_CenterL(const IndexDimension& center,const IndexDimension& cell, int d) const{
1212 double p =center.coordinate(d);
1213 double h=meshwidth[d];
1214 double x = cell.coordinate(d);
1215
1216 double val1=(1.0/3.0)*x*x*x-(p-h)*x*x+(p-h)*(p-h)*x;
1217 double val12=(1.0/5.0)*x*x*x*x*x-(2.0/4.0)*(p-h)*x*x*x*x+(1.0/3.0)*(p-h)*(p-h)*x*x*x;
1218 val1 = val1 - val12;
1219
1220
1221 x = x+h;
1222 double val2=(1.0/3.0)*x*x*x-(p-h)*x*x+(p-h)*(p-h)*x;
1223 double val22=(1.0/5.0)*x*x*x*x*x-(2.0/4.0)*(p-h)*x*x*x*x+(1.0/3.0)*(p-h)*(p-h)*x*x*x;
1224 val2= val2 - val22;
1225
1226 return val2-val1;
1227
1228 }
1229
1230 // integrate (-x+(p+h))²*(1-x²)
1231 inline double integral_CenterR(const IndexDimension& center,const IndexDimension& cell, int d) const{
1232
1233 double p =center.coordinate(d);
1234 double h=meshwidth[d];
1235 double x = cell.coordinate(d);
1236
1237
1238 double val1= (1.0/3.0)*x*x*x-(p+h)*x*x+(p+h)*(p+h)*x;
1239 double val12=(1.0/5.0)*x*x*x*x*x-(2.0/4.0)*(p+h)*x*x*x*x+(1.0/3.0)*(p+h)*(p+h)*x*x*x;
1240 val1 = val1 - val12;
1241
1242 x = x+h;
1243 double val2= (1.0/3.0)*x*x*x-(p+h)*x*x+(p+h)*(p+h)*x;
1244 double val22=(1.0/5.0)*x*x*x*x*x-(2.0/4.0)*(p+h)*x*x*x*x+(1.0/3.0)*(p+h)*(p+h)*x*x*x;
1245 val2= val2 - val22;
1246 return val2-val1;
1247
1248 }
1249
1250
1255 inline double integral_right(const IndexDimension& center,const IndexDimension& cell, int d) const{
1256 double p =center.coordinate(d);
1257 double h=meshwidth[d];
1258 double x = cell.coordinate(d);
1259
1260 double val1=-(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
1261 double val12=(1.0/5.0)*x*x*x*x*x-(1.0/4.0)*p*x*x*x*x-(1.0/4.0)*(p+h)*x*x*x*x+(1.0/3.0)*p*(p+h)*x*x*x;
1262 val1 = val1 + val12;
1263
1264 x = x+h;
1265 double val2=-(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
1266 double val22=(1.0/5.0)*x*x*x*x*x-(1.0/4.0)*p*x*x*x*x-(1.0/4.0)*(p+h)*x*x*x*x+(1.0/3.0)*p*(p+h)*x*x*x;
1267 val2= val2 + val22;
1268 return val2-val1;
1269
1270 }
1271
1272
1273 inline double integral_RR(CellDimension& cell,int d){
1274 double b = cell.getRightBoundary(d);
1275 double a = cell.getLeftBoundary(d);
1276 double h= b-a;
1277
1278 //integrate (x-b)²(1-x²) from a to b
1279 double val1=(-1.0*h*h*h*(-10.0 + 6.0*a*a + 3.0*a*b + b*b))/30.0;
1280 return val1;
1281
1282 }
1283
1284 inline double integral_LL(CellDimension& cell,int d){
1285 double b = cell.getRightBoundary(d);
1286 double a = cell.getLeftBoundary(d);
1287 double h= b-a;
1288
1289 //integrate (a-x)²(1-x²) from a to b
1290 double val1=(-1.0*h*h*h*(-10.0 + 6.0*b*b + 3.0*a*b + a*a))/30.0;
1291 return val1;
1292
1293 }
1294 inline double integral_LR(CellDimension& cell,int d){
1295 double b = cell.getRightBoundary(d);
1296 double a = cell.getLeftBoundary(d);
1297 double h= b-a;
1298
1299 //integrate (x-b)(a-x)(1-x²) from a to b
1300 double val1=(-1.0*h*h*h*(-10.0 + 3.0*a*a + 4.0*a*b + 3.0*b*b))/60.0;
1301
1302 return val1;
1303
1304 }
1305
1306
1307};
1308
1309
1310
1311
1315class ConstantCoefficientMass: public StencilTemplate{
1316public:
1317 ConstantCoefficientMass(AdaptiveSparseGrid& grid_): StencilTemplate(grid_){};
1318
1319 double c=1.0;
1320
1321 inline double localValue(const IndexDimension &center ,const IndexDimension& cell, const MultiDimCompass &mc)const{
1322
1323
1324 double value = 1.0;
1325 for(int d=0; d<DimensionSparseGrid; d++) {
1326 if (mc.getRichtung(d) == Mitte && cell.getIndex(d) == center.getIndex(d))
1327 value *= integral_CenterR(center,cell, d);
1328 if (mc.getRichtung(d) == Mitte && !(cell.getIndex(d) == center.getIndex(d)))
1329 value *= integral_CenterL(center,cell, d);
1330 if (mc.getRichtung(d) == Links && !(cell.getIndex(d) == center.getIndex(d)))
1331 value *= integral_left(center,cell, d);
1332 if (mc.getRichtung(d) == Rechts && cell.getIndex(d) == center.getIndex(d))
1333 value *= integral_right(center,cell, d);
1334
1335 value = exponent[d]*exponent[d]*value;
1336 }
1337
1338
1339
1340
1341 return value;}
1342
1343 inline double integration(CellDimension& cell, CellIndexDirection& dirP, CellIndexDirection& dirQ){
1344 double value = 1.0;
1345
1346 for(int d=0; d<DimensionSparseGrid; d++) {
1347
1348
1349 if (dirP.getDirection(d)==dirQ.getDirection(d)){
1350 value *= integral_RR(cell, d);
1351 value = exponent[d] * exponent[d] * value;
1352
1353 }
1354 else {
1355 value *= integral_LR(cell,d);
1356 value = exponent[d] * exponent[d] * value;
1357 }
1358
1359 }
1360 return c*value;
1361
1362 }
1363
1364protected:
1365
1366 // integrate (-x+p)*(x-(p-h))
1367 inline double integral_left(const IndexDimension& center,const IndexDimension& cell, int d) const{
1368
1369 double p =center.coordinate(d);
1370 double h=meshwidth[d];
1371 double x = cell.coordinate(d);
1372 double val1=-(1.0/3.0)*x*x*x+0.5*(p-h)*x*x+0.5*p*x*x-p*(p-h)*x;
1373
1374 x = x+h;
1375 double val2=-(1.0/3.0)*x*x*x+0.5*(p-h)*x*x+0.5*p*x*x-p*(p-h)*x;
1376 return val2-val1;
1377 }
1378
1379 // integrate (-x+(p+h))*(x-p)
1380 inline double integral_right(const IndexDimension& center,const IndexDimension& cell, int d) const{
1381 double p =center.coordinate(d);
1382 double h=meshwidth[d];
1383 double x = cell.coordinate(d);
1384
1385 double val1=-(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
1386
1387 x = x+h;
1388 double val2=-(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
1389
1390 return val2-val1;
1391 //return -(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
1392 }
1393
1394 inline double integral_CenterR(const IndexDimension& center,const IndexDimension& cell, int d) const{
1395
1396 double p =center.coordinate(d);
1397 double h=meshwidth[d];
1398 double x = cell.coordinate(d);
1399
1400
1401
1402 double val1= (1.0/3.0)*x*x*x-(p+h)*x*x+(p+h)*(p+h)*x;
1403
1404
1405 x = x+h;
1406 double val2= (1.0/3.0)*x*x*x-(p+h)*x*x+(p+h)*(p+h)*x;
1407
1408 return val2-val1;
1409
1410 }
1411
1412 // integrate (x-(p-h))²
1413 inline double integral_CenterL(const IndexDimension& center,const IndexDimension& cell, int d) const{
1414 double p = center.coordinate(d);
1415 double h=meshwidth[d];
1416 double x = cell.coordinate(d);
1417 double val1=(1.0/3.0)*x*x*x-(p-h)*x*x+(p-h)*(p-h)*x;
1418
1419
1420 x = x+h;
1421 double val2=(1.0/3.0)*x*x*x-(p-h)*x*x+(p-h)*(p-h)*x;
1422
1423
1424 return val2-val1;
1425
1426 }
1427
1428
1429
1430
1431};
1432
1433class StencilVariableCoeff{
1434public:
1435 inline void initialize(Depth &T_) {
1436
1437 T = T_;
1438 poisson.initialize(T);
1439
1440 for (int d = 0; d < DimensionSparseGrid; d++) {
1441 int exp = T_.at(d);
1442 if (exp == 0) meshwidth[d] = 1.0;
1443 else {
1444
1445 int value2 = 1 << exp;
1446 exponent[d] = value2;
1447 meshwidth[d] = 1 / double(value2);
1448
1449 }
1450 }
1451
1452 };
1453
1454
1455
1456 inline double returnValue(const IndexDimension &Index, const MultiDimCompass &mc){
1457 double val = 0.0;
1458
1459
1460 for(int d=0; d<DimensionSparseGrid; d++)
1461 coord[d]=Index.coordinate(d);
1462
1463 CellIterator cells(mc,T,Index);
1464
1465
1466 for(cells.start(); cells.goon();++cells){
1467 IndexDimension cell = cells.getCell();
1468 val+= localValue(Index,cell,mc);
1469
1470 }
1471
1472
1473
1474
1475 val += poisson.returnValue(Index,mc);
1476 return val;
1477
1478
1479
1480 };
1481 // integrate (-x+p)*(x-(p-h))*(1-x²)
1482 double integral_left(IndexDimension& cell, int d) const{
1483
1484 double p = coord[d];
1485 double h=meshwidth[d];
1486
1487 double x = cell.coordinate(d);
1488 double val1=-(1.0/3.0)*x*x*x+0.5*(p-h)*x*x+0.5*p*x*x-p*(p-h)*x;
1489 double val12=-(1.0/5.0)*x*x*x*x*x+(1.0/4.0)*p*x*x*x*x+(1.0/4.0)*(p-h)*x*x*x*x-(1.0/3.0)*p*(p-h)*x*x*x;
1490
1491 val1 = val1 - val12;
1492 x = x+h;
1493 double val2=-(1.0/3.0)*x*x*x+0.5*(p-h)*x*x+0.5*p*x*x-p*(p-h)*x;
1494 double val22=-(1.0/5.0)*x*x*x*x*x+(1.0/4.0)*p*x*x*x*x+(1.0/4.0)*(p-h)*x*x*x*x-(1.0/3.0)*p*(p-h)*x*x*x;
1495
1496 val2= val2 - val22;
1497
1498 return val2-val1;
1499 }
1500
1501 // integrate (x-(p-h))²
1502 double integral_CenterL(IndexDimension& cell, int d) const{
1503 double p = coord[d];
1504 double h=meshwidth[d];
1505 double x = cell.coordinate(d);
1506
1507 double val1=(1.0/3.0)*x*x*x-(p-h)*x*x+(p-h)*(p-h)*x;
1508 double val12=(1.0/5.0)*x*x*x*x*x-(2.0/4.0)*(p-h)*x*x*x*x+(1.0/3.0)*(p-h)*(p-h)*x*x*x;
1509 val1 = val1 - val12;
1510
1511
1512 x = x+h;
1513 double val2=(1.0/3.0)*x*x*x-(p-h)*x*x+(p-h)*(p-h)*x;
1514 double val22=(1.0/5.0)*x*x*x*x*x-(2.0/4.0)*(p-h)*x*x*x*x+(1.0/3.0)*(p-h)*(p-h)*x*x*x;
1515 val2= val2 - val22;
1516
1517 return val2-val1;
1518
1519 }
1520
1521 // integrate (-x+(p+h))²*(1-x²)
1522 double integral_CenterR(IndexDimension& cell, int d) const{
1523
1524 double p = coord[d];
1525 double h=meshwidth[d];
1526 double x = cell.coordinate(d);
1527
1528
1529 double val1= (1.0/3.0)*x*x*x-(p+h)*x*x+(p+h)*(p+h)*x;
1530 double val12=(1.0/5.0)*x*x*x*x*x-(2.0/4.0)*(p+h)*x*x*x*x+(1.0/3.0)*(p+h)*(p+h)*x*x*x;
1531 val1 = val1 - val12;
1532
1533 x = x+h;
1534 double val2= (1.0/3.0)*x*x*x-(p+h)*x*x+(p+h)*(p+h)*x;
1535 double val22=(1.0/5.0)*x*x*x*x*x-(2.0/4.0)*(p+h)*x*x*x*x+(1.0/3.0)*(p+h)*(p+h)*x*x*x;
1536 val2= val2 - val22;
1537 return val2-val1;
1538
1539 }
1540
1541
1542 // integrate (-x+(p+h))*(x-p)*(1-x²)
1543 double integral_right(IndexDimension& cell, int d) const{
1544 double p = coord[d];
1545 double h=meshwidth[d];
1546 double x = cell.coordinate(d);
1547
1548 double val1=-(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
1549 double val12=(1.0/5.0)*x*x*x*x*x-(1.0/4.0)*p*x*x*x*x-(1.0/4.0)*(p+h)*x*x*x*x+(1.0/3.0)*p*(p+h)*x*x*x;
1550 val1 = val1 + val12;
1551
1552 x = x+h;
1553 double val2=-(1.0/3.0)*x*x*x+(1.0/2.0)*p*x*x+(1.0/2.0)*(p+h)*x*x-p*(p+h)*x;
1554 double val22=(1.0/5.0)*x*x*x*x*x-(1.0/4.0)*p*x*x*x*x-(1.0/4.0)*(p+h)*x*x*x*x+(1.0/3.0)*p*(p+h)*x*x*x;
1555 val2= val2 + val22;
1556 return val2-val1;
1557
1558 }
1559
1560
1561private:
1562 double localValue(const IndexDimension &center ,IndexDimension& cell,const MultiDimCompass& mc)const{
1563
1564
1565 double value = 1.0;
1566 for(int d=0; d<DimensionSparseGrid; d++) {
1567 if (mc.getRichtung(d) == Mitte && cell.getIndex(d) == center.getIndex(d))
1568 value *= integral_CenterR(cell, d);
1569 if (mc.getRichtung(d) == Mitte && !(cell.getIndex(d) == center.getIndex(d)))
1570 value *= integral_CenterL(cell, d);
1571 if (mc.getRichtung(d) == Links && !(cell.getIndex(d) == center.getIndex(d)))
1572 value *= integral_left(cell, d);
1573 if (mc.getRichtung(d) == Rechts && cell.getIndex(d) == center.getIndex(d))
1574 value *= integral_right(cell, d);
1575
1576 value = exponent[d]*exponent[d]*value;
1577 }
1578
1579
1580
1581
1582 return value;};
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594 Depth T;
1595 double meshwidth[DimensionSparseGrid];
1596 double exponent[DimensionSparseGrid];
1597 double coord[DimensionSparseGrid];
1598
1599 PoissonStencil poisson;
1600
1601};
1602
1603
1604
1605class Poisson : public StencilTemplate{
1606public:
1607 Poisson(AdaptiveSparseGrid& grid_): StencilTemplate(grid_){};
1608
1609 TypeMatrixVectorMultiplication getTypeMatrixVectorMultiplication(){
1610 return StencilOnTheFly;
1611 }
1612
1613
1614 inline double localValue(const IndexDimension &center ,const IndexDimension& cell, const MultiDimCompass &mc)const{
1615
1616
1617 double value =0.0;
1618 double ddx;
1619 for(int dx=0; dx < DimensionSparseGrid; dx++) {
1620 ddx=1.0;
1621 if (mc.getRichtung(dx) == Mitte && cell.getIndex(dx) == center.getIndex(dx)){
1622
1623 ddx *= exponent[dx];
1624 ddx *= Integral_Mass(center,cell,mc,dx);
1625
1626 value += ddx;
1627
1628 }
1629
1630 if (mc.getRichtung(dx) == Mitte && !(cell.getIndex(dx) == center.getIndex(dx))){
1631
1632 ddx *= exponent[dx];
1633 ddx *= Integral_Mass(center,cell,mc,dx);
1634
1635 value += ddx;
1636
1637 }
1638
1639 if (mc.getRichtung(dx) == Links && !(cell.getIndex(dx) == center.getIndex(dx))){
1640
1641 ddx *= -1.0*exponent[dx];
1642 ddx *= Integral_Mass(center,cell,mc,dx);
1643
1644 value += ddx;
1645
1646 }
1647
1648 if (mc.getRichtung(dx) == Rechts && cell.getIndex(dx) == center.getIndex(dx)){
1649
1650 ddx *= -1.0*exponent[dx];
1651 ddx *= Integral_Mass(center,cell,mc,dx);
1652
1653 value += ddx;
1654
1655 }
1656
1657 }
1658
1659
1660 return value;}
1661
1662 inline double Integral_Mass(const IndexDimension &center ,const IndexDimension& cell, const MultiDimCompass &mc, int D)const{
1663
1664
1665 double value = 1.0;
1666 for(int d=0; d<DimensionSparseGrid; d++) {
1667
1668 if(D!=d) {
1669 if (mc.getRichtung(d) == Mitte && cell.getIndex(d) == center.getIndex(d)){
1670
1671
1672 value *= integral_CenterR(center,cell, d);
1673
1674 // links links
1675 value = exponent[d] * exponent[d] * value;
1676
1677 }
1678 else if (mc.getRichtung(d) == Mitte && (cell.getIndex(d) != center.getIndex(d))){
1679 value *= integral_CenterL(center,cell, d);
1680
1681 // rechts rechts
1682 value = exponent[d] * exponent[d] * value;
1683
1684 }
1685
1686 else if (mc.getRichtung(d) == Links && (cell.getIndex(d) != center.getIndex(d))){
1687 value *= integral_left(center,cell, d);
1688
1689 //
1690 value = exponent[d] * exponent[d] * value;
1691
1692 }
1693
1694 else if (mc.getRichtung(d) == Rechts && cell.getIndex(d) == center.getIndex(d)){
1695 value *= integral_right(center,cell, d);
1696 value = exponent[d] * exponent[d] * value;
1697
1698 }
1699
1700
1701
1702
1703 }
1704 }
1705
1706
1707
1708
1709 return value;
1710 }
1711
1712
1713 inline double integration(CellDimension& cell, CellIndexDirection& dirP, CellIndexDirection& dirQ){
1714
1715 double h[DimensionSparseGrid];
1716
1717 for(int i=0;i<DimensionSparseGrid;++i) h[i] =cell.getRightBoundary(i) - cell.getLeftBoundary(i);
1718
1719/* double sum=0.0;
1720 for(int i=0;i<DimensionSparseGrid;++i) { // int du/dxi * dv/dxi d(x1...xd)
1721 double prod=1.0;
1722 for(int j=0;j<DimensionSparseGrid;++j) { // integral factor in direction j
1723 if(i==j) { // factor: int du/dxi * dv/dxi dxi
1724 if(dirP.getDirection(j)==dirQ.getDirection(j)) prod = prod / h[j];
1725 else prod = - prod / h[j];
1726 }
1727 else { // factor: int du/dxi * dv/dxi dxj == int u*v dxj
1728 if(dirP.getDirection(j)==dirQ.getDirection(j)) prod = prod * (1.0/3.0) * h[j];
1729 else prod = prod * (1.0/6.0) * h[j];
1730 }
1731 }
1732 sum = sum + prod;
1733 }
1734 return sum;*/
1735 double value =0.0;
1736 double ddx;
1737 for(int dx=0; dx < DimensionSparseGrid; dx++) {
1738 ddx=1.0;
1739 if (dirP.getDirection(dx)==dirQ.getDirection(dx)){
1740
1741
1742 ddx *= exponent[dx];
1743 ddx *= Integral_Mass(cell,dirP,dirQ,dx);
1744
1745
1746 value += ddx;
1747
1748
1749 } else {
1750
1751
1752 ddx *= -1.0*exponent[dx];
1753 ddx *= Integral_Mass(cell,dirP,dirQ,dx);
1754
1755 value += ddx;
1756
1757
1758 }
1759
1760 }
1761
1762
1763 return value ; // sum = sum_i int du/dxi * dv/dxi d(x1...xd)
1764 }
1765
1766
1767 inline double Integral_Mass(CellDimension& cell, CellIndexDirection& dirP, CellIndexDirection& dirQ, int D) {
1768
1769
1770 double value = 1.0;
1771 for(int d=0; d<DimensionSparseGrid; d++) {
1772
1773 if(D!=d) {
1774 if (dirP.getDirection(d)==dirQ.getDirection(d)){
1775 if(dirP.getDirection(d)==Right) value *= integral_RR(cell, d);
1776 if(dirP.getDirection(d)==Left) value *= integral_LL(cell,d);
1777
1778
1779 value = exponent[d] * exponent[d] * value;
1780
1781 }
1782 else {
1783 value *= integral_LR(cell,d);
1784 value = exponent[d] * exponent[d] * value;
1785 }
1786 }
1787 }
1788 return value;
1789 }
1790
1791
1792
1793};
1794
1795
1796
1797
1798template<class StencilA,class StencilB>
1799class ADD:public StencilTemplate{
1800public:
1801 ADD(StencilA stencilA, StencilB stencilB,AdaptiveSparseGrid& grid_):StencilTemplate(grid_),stencilP(stencilA), stencilM(stencilB) {};
1802 virtual inline double returnValue(const IndexDimension &Index, const MultiDimCompass &mc){
1803
1804 double val = 0.0;
1805
1806 CellIterator cells(mc,stencilM.T,Index);
1807
1808
1809 for(cells.start(); cells.goon();++cells){
1810 IndexDimension cell = cells.getCell();
1811 if(stencilM.cellInGrid(cell)) {
1812
1813
1814 val += stencilM.localValue(Index, cell, mc);
1815 val += stencilP.localValue(Index, cell, mc);
1816
1817
1818
1819 }
1820
1821 }
1822
1823
1824
1825
1826
1827 return val;
1828
1829
1830
1831 };
1832
1833 inline void initialize(Depth &T_) {
1834
1835 stencilP.initialize(T_);
1836 stencilM.copyInitialization(stencilP);
1837
1838
1839
1840 };
1841
1842/*
1843 inline void applyStencilOnCell(CellDimension& cell,VectorSparseG& input, VectorSparseG& output){
1844 stencilM.applyStencilOnCell(cell,input,output);
1845 stencilP.applyStencilOnCell(cell,input,output);
1846
1847
1848 }
1849*/
1850 inline double integration(CellDimension& cell, CellIndexDirection& dirP, CellIndexDirection& dirQ){
1851
1852 return stencilP.integration(cell,dirP,dirQ)+stencilM.integration(cell,dirP,dirQ);
1853 }
1854/* inline void applyStencilOnCell(CellDimension& cell,VectorSparseG& input, VectorSparseG& output){
1855
1856 for(CellIndexIterator outerIter(&cell); outerIter.goon(); ++outerIter) {
1857
1858
1859 IndexDimension p = outerIter.getIndex();
1860
1861 unsigned long kp;
1862 bool occP = output.getSparseGrid()->occupied(kp, p);
1863 CellIndexDirection dirP = outerIter.getCellIndexDirection();
1864
1865
1866 if (occP) {
1867
1868 for (CellIndexIterator innerIter(outerIter); innerIter.goon(); ++innerIter) {
1869
1870
1871 IndexDimension q = innerIter.getIndex();
1872 if (q.isNotAtBoundary()) {
1873
1874 CellIndexDirection dirQ = innerIter.getCellIndexDirection();
1875
1876 unsigned long kq;
1877
1878 bool occQ = output.getSparseGrid()->occupied(kq, q);
1879 if (occQ) {
1880
1881 double val = stencilM.integration(cell, dirP, dirQ)+stencilP.integration(cell, dirP, dirQ);
1882 if (p == q) {
1883 val = val * input.getValue(kp);
1884
1885 //#pragma omp critical
1886 output.addToValue(kq, val);
1887 } else {
1888 double val_p = val * input.getValue(kp);
1889 //#pragma omp critical
1890 output.addToValue(kq, val_p);
1891
1892 double val_q = val * input.getValue(kq);
1893 //#pragma omp critical
1894 output.addToValue(kp, val_q);
1895 }
1896 }
1897
1898 }
1899
1900 }
1901
1902 }
1903 }
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913 }*/
1914/* inline void applyStencilOnCell_BinaryIterator(CellDimension& cell,VectorSparseG& input, VectorSparseG& output){
1915 stencilM.applyStencilOnCell_BinaryIterator(cell,input,output);
1916 stencilP.applyStencilOnCell_BinaryIterator(cell,input,output);
1917 }*/
1918/* inline void applyStencilOnCell_MPI_OMP(CellDimension& cell,VectorSparseG& input, VectorSparseG& output){
1919 stencilM.applyStencilOnCell_MPI_OMP(cell,input,output);
1920 stencilP.applyStencilOnCell_MPI_OMP(cell,input,output);
1921
1922
1923 }*/
1924/* inline void applyStencilOnCell_MPI_OMP_nosymmetry(CellDimension& cell,VectorSparseG& input, VectorSparseG& output){
1925 stencilM.applyStencilOnCell_MPI_OMP_nosymmetry(cell,input,output);
1926 stencilP.applyStencilOnCell_MPI_OMP_nosymmetry(cell,input,output);
1927
1928
1929 }*/
1930private:
1931
1932 StencilA stencilP;
1933 StencilB stencilM;
1934
1935};
1936
1937
1938
1939#endif //SGRUN_STENCIL_H
AdaptiveSparseGrid
Definition sparseGrid.h:277
ConstantCoefficientMass
Definition Stencil.h:1315
IntegratorHelm
Integrator for the Helmholtz operator with variable coefficients.
Definition helmIntegrator.h:27
MultiDimCompass
Definition index.h:356
RectangularIterator
Definition RectangularIterator.h:13
VariableCoefficientMass
Definition Stencil.h:1127
VariableCoefficientMass::integral_right
double integral_right(const IndexDimension &center, const IndexDimension &cell, int d) const
Definition Stencil.h:1255