LoAdSG
MatrixVectorInhomogen.h
1//
2// Created by to35jepo on 12/7/22.
3//
4
5#ifndef SGRUN_MATRIXVECTORINHOMOGEN_H_H
6#define SGRUN_MATRIXVECTORINHOMOGEN_H_H
7#include "../extemp/vector.h"
8#include "../tests/old_versions/MatrixVectorMultiplicationPrewavelets.h"
9#include "../sgrid/multiDepthHashGrid.h"
10#include "../BasisTransformations/BasisTransformations.h"
11#include "../stencils/PoissonStencil.h"
12#include "../iterator/depthIterator.h"
13
14
15class MatrixVectorInhomogen { ;
16public:
17 MatrixVectorInhomogen(AdaptiveSparseGrid& grid, MultiLevelAdaptiveSparseGrid& mgrid) : gM(mgrid), g(grid), z(grid), u(grid),
18 nodal(mgrid),
19 u_new(grid),
20 u_old(grid), Ax_neu(grid), list(grid), P(mgrid){
21
22
23
24
25 };
26
27
28 template<class Problem>
29 void multiplication(VectorSparseG &prew, VectorSparseG &Ax, Problem &matrixProblem);
30
31
41 template<class Problem>
42 void multiplicationRHSHomogen(VectorSparseG &prew, VectorSparseG &Ax, Problem &matrixProblem);
43
44
45
46
47 // hier wird if(type==irgendwas)
48 //void multiplication(VectorSparseG &prew, VectorSparseG &Ax, StencilType type);
49
50
51
52
53private:
54
55 MultiLevelVector gM;
56 MultiLevelVector nodal;
57 MultiLevelVector P;
58 VectorSparseG Ax_neu;
59 VectorSparseG g;
60 VectorSparseG z;
61 VectorSparseG u;
62 VectorSparseG u_new;
63 VectorSparseG u_old;
64
65 ListOfDepthOrderedSubgrids list;
66
67 template<class Problem>
68 void CaseFunction(VectorSparseG &prew, bool *restrictions, Problem &matrixProblem);
69
70 template<class Problem>
71 void CaseFunctionTransformation(VectorSparseG &prew, bool *restrictions, Problem &matrixProblem);
72
73 inline void calcU(VectorSparseG &prew, Depth &Tiefe, bool *restrictions);
74
75 inline void calcNodal(VectorSparseG &prew, Depth &Tiefe);
76
77 template<class Problem>
78 void applyProblemGlobal(VectorSparseG &input, VectorSparseG &output, Problem &matrixProblem, Depth &T) {
79 AdaptiveSparseGrid_Base *grid = input.getSparseGrid();
80 int maxocc = grid->getMaximalOccupiedSecondTable();
81 for (unsigned long i = 0; i < maxocc; i++) {
82 IndexDimension Index = grid->getIndexOfTable(i);
83 Depth Tlocal(Index);
84 if (Tlocal <= T) {
85
86 double val = applyProblemLocal(Index, input, matrixProblem, T);
87 //output = val | Index;
88 output.setValue(i, val);
89 }
90 }
91 };
92
93 template<class Problem>
94 double applyProblemLocal(IndexDimension &Index, VectorSparseG &input, Problem &matrixProblem, Depth &T) {
95
96 double val = 0.0;
97 IndexDimension NextTest;
98 bool exist_index;
99 for (MultiDimCompass mc; mc.goon(); ++mc) {
100 NextTest = Index.nextThree_boundary(&mc, T.returnTiefen());
101 //exist_index = getNextIndex(Index, mc, T, NextTest);
102 if(mc.goon()) {
103 double value = 0.0;
104 value = matrixProblem.returnValue(Index, mc);
105 val = val + value * input.getValue(NextTest);
106 }
107 }
108 return val;
109 };
110
111 inline void prolongation1D(VectorSparseG &coarse, VectorSparseG &fine, int &t, int &d) {
112 AdaptiveSparseGrid_Base *grid = fine.getSparseGrid();
113 unsigned long maxocc = grid->getMaximalOccupiedSecondTable();
114
115
116 for (unsigned long i = 0; i < maxocc; i++) {
117 IndexDimension Index = grid->getIndexOfTable(i);
118
119 if (Index.getDepth(d) == t) {
120 double value = fine.getValue(i);
121/*
122 if (Index.isAtLeftBoundary(d))
123 value = fine.getValue(Index);
124
125 if (Index.isAtRightBoundary(d))
126 value = fine.getValue(Index);
127*/
128
129 if ((!Index.isAtRightBoundary(d)) && (!Index.isAtLeftBoundary(d))) {
130 value = value + 0.5 * (fine.getValue(Index.nextRight(d)) + fine.getValue(Index.nextLeft(d)));
131
132 }
133 //coarse.setValue(Index, value);
134 coarse = value | Index;
135
136
137 }
138
139
140 }
141
142
143 }
144
145 void prolongation1D_inplace(VectorSparseG &vec, int &t, int &d) {
146 AdaptiveSparseGrid_Base *grid = vec.getSparseGrid();
147
148
149 ListOfDepthOrderedSubgrids::iterator iter(list);
150 iter.gotoBegin();
151 do {
152 Depth Tlocal = iter.getDepth();
153 if (Tlocal.at(d) == t) {
154 SubgridFixedDepth::iterator inneriter(*iter.getSubgrid());
155 inneriter.gotobegin();
156 do {
157 IndexDimension Index = inneriter.getPoint();
158 unsigned long i = inneriter.geti();
159 double value = vec.getValue(i);
160 if ((!Index.isAtRightBoundary(d)) && (!Index.isAtLeftBoundary(d)))
161 value = value + 0.5 * (vec.getValue(Index.nextRight(d)) + vec.getValue(Index.nextLeft(d)));
162 //coarse = value | Index;
163 vec.setValue(i, value);
164
165 } while (inneriter.next());
166 }
167
168 } while (iter.next());
169 }
170
171
172
173 inline void prolongation(VectorSparseG &coarse, VectorSparseG &fine, Depth &Tcoarse, Depth &Tfine) {
174 fine = coarse;
175 int tcoarse;
176 int tfine;
177
178 VectorSparseG testvec(*coarse.getSparseGrid());
179 testvec = coarse;
180 VectorSparseG testvec2(*coarse.getSparseGrid());
181
182 for (int d = 0; d < DimensionSparseGrid; d++) {
183
184 tcoarse = Tcoarse.at(d) + 1;
185 tfine = Tfine.at(d);
186 for (int t = tcoarse; t <= tfine; ++t) {
187 //prolongation1D(fine, fine, t, d);
188 //prolongation1D_inplace(fine,t,d);
189 prolongation1D_inplace_singleHash(list,fine,t,d);
190 /* testvec2 = testvec-fine;
191 if(L_infty(testvec2)>1e-10) {
192 for(unsigned long k=0; k < testvec.getSparseGrid()->getMaximalOccupiedSecondTable(); k++){
193 if(abs(testvec2.getValue(k))>1e-10) {
194 testvec2.getSparseGrid()->getIndexOfTable(k).Print();
195 cout << testvec2.getValue(k) << endl;
196 cout << fine.getValue(k) << endl;
197 cout << testvec.getValue(k) << endl;
198 exit(EXIT_FAILURE);
199 }
200 }
201
202 }*/
203 }
204 }
205 };
206
207 static void prolongation1D_inplace_singleHash(ListOfDepthOrderedSubgrids& list, VectorSparseG &vec,const int &t,const int &d) {
208 AdaptiveSparseGrid_Base *grid = vec.getSparseGrid();
209
210
211 ListOfDepthOrderedSubgrids::iterator iter(list);
212 iter.gotoBegin();
213 do {
214
215 Depth Tlocal = iter.getDepth();
216
217
218 if (Tlocal.at(d) == t) {
219
220
221 vector<SingleDepthHashGrid*> fineDepthGrids = grid->getMultiDepthHashGrid()->getGridsForDepthInDirection(Tlocal,d);
222 SubgridFixedDepth::iterator inneriter(*iter.getSubgrid());
223 inneriter.gotobegin();
224 do {
225 IndexDimension Index = inneriter.getPoint();
226 unsigned long i = inneriter.geti();
227 double value = vec.getValue(i);
228 if ((!Index.isAtRightBoundary(d)) && (!Index.isAtLeftBoundary(d))){
229 auto right = Index.nextRight(d);
230 auto left = Index.nextLeft(d);
231 // auto* gridRight = fineDepthGrids[right.getDepth(d)];
232 // auto* gridLeft = fineDepthGrids[left.getDepth(d)];
233 // auto valRight = vec.getValue(right,gridRight);
234 // auto valLeft = vec.getValue(left,gridLeft);
235 value = value + 0.5 * (vec.getValue(right,fineDepthGrids[right.getDepth(d)]) + vec.getValue(left,fineDepthGrids[left.getDepth(d)]));
236
237 // value = value + 0.5 * (vec.getValue(Index.nextRight(d)) + vec.getValue(Index.nextLeft(d)));
238 // auto p= vec.getValue(right,fineDepthGrids[right.getDepth(d)]);
239 }
240 //coarse = value | Index;
241 vec.setValue(i, value);
242 // cout << Depth(Index.nextRight(d)).at(d)<<", \t" <<Depth(Index.nextLeft(d)).at(d) << endl;
243
244 // cout << "Left: " ;
245 // Depth(Index.nextLeft(d)).Print();
246 // cout << "Right: " ;
247 // Depth(Index.nextRight(d)).Print();
248 } while (inneriter.next());
249
250
251 }
252
253 // cout << "iter" << endl;
254 } while (iter.next());
255 // cout << "end" << endl;
256
257
258 }
259
260 inline void restriction1D(VectorSparseG &fine, VectorSparseG &coarse, int t, int d) {
261
262 AdaptiveSparseGrid_Base *grid = fine.getSparseGrid();
263 unsigned long maxocc = grid->getMaximalOccupiedSecondTable();
264 double value;
265 double valueL;
266 double valueR;
267 for (unsigned long i = 0; i < maxocc; i++) {
268 IndexDimension Index = grid->getIndexOfTable(i);
269
270 if (Index.getDepth(d) <= t) {
271
272 value = fine.getValue(i);
273 valueL = 0.0;
274 valueR = 0.0;
275 if (t >= 0) {
276 if (!(Index.isAtRightBoundary(d))) valueL = 0.5 * fine.getValue(Index.nextRight(d, t + 1));
277 if (!(Index.isAtLeftBoundary(d))) valueR = 0.5 * fine.getValue(Index.nextLeft(d, t + 1));
278 }
279
280
281 value = value + valueL + valueR;
282 //coarse.setValue(Index, value);
283 coarse.setValue(i, value);
284 //coarse = value | Index;
285 }
286
287
288 }
289
290 };
291
292
293 inline void restriction1D_inverted_inplace( const VectorSparseG &fine,const VectorSparseG &coarse, int t, int d) {
294
295 double value;
296
297 ListOfDepthOrderedSubgrids::iterator iter(list);
298 iter.gotoBegin();
299 do {
300 Depth Tlocal = iter.getDepth();
301 // SingleDepthHashGrid& coarseDepthGrid = grid->getMultiDepthHashGrid()->getGridForDepth(Tlocal);
302 if(Tlocal.at(d)<=t){
303
304 SubgridFixedDepth::iterator inneriter(*iter.getSubgrid());
305 inneriter.gotobegin();
306 // Depth fineDepth = Depth(inneriter.getPoint());
307 // SingleDepthHashGrid& fineDepthGrid = grid->getMultiDepthHashGrid()->getGridForDepth(Depth(inneriter.getPoint().nextRight(d,t+1)));
308 do {
309
310 unsigned long i = inneriter.geti();
311 double value = fine.getValue(i);
312 coarse.setValue(i,value);
313
314 } while (inneriter.next());
315
316 }
317
318 if (Tlocal.at(d) == t + 1) {
319
320 SubgridFixedDepth::iterator inneriter(*iter.getSubgrid());
321 inneriter.gotobegin();
322 // Depth fineDepth = Depth(inneriter.getPoint());
323 // SingleDepthHashGrid& fineDepthGrid = grid->getMultiDepthHashGrid()->getGridForDepth(Depth(inneriter.getPoint().nextRight(d,t+1)));
324 do {
325
326 IndexDimension Index = inneriter.getPoint();
327 unsigned long i = inneriter.geti();
328 value = 0.5 * fine.getValue(i);
329
330 if (t >= 0) {
331
332 IndexDimension rightIndex = Index.nextRight(d);
333 if (!Index.isAtRightBoundary(d)){
334 // vec.getValue(rightIndex,value,singleDepthHashGrid)
335 coarse.addToValue(rightIndex,value);
336 // value += vec.getValue(rightIndex);
337 }
338
339
340 IndexDimension leftIndex = Index.nextLeft(d);
341 if (!Index.isAtLeftBoundary(d)){
342 // vec.getValue(rightIndex,value,singleDepthHashGrid)
343 coarse.addToValue(leftIndex,value);
344 // value += vec.getValue(rightIndex);
345 }
346 }
347 } while (inneriter.next());
348 }
349 } while (iter.next());
350
351 };
352
353
354 inline void restriction1D_inverted_inplace( const VectorSparseG &vec, int t, int d) {
355
356 double value;
357
358 ListOfDepthOrderedSubgrids::iterator iter(list);
359 iter.gotoBegin();
360 do {
361 Depth Tlocal = iter.getDepth();
362 if (Tlocal.at(d) == t + 1) {
363
364 SubgridFixedDepth::iterator inneriter(*iter.getSubgrid());
365 inneriter.gotobegin();
366 // Depth fineDepth = Depth(inneriter.getPoint());
367 // SingleDepthHashGrid& fineDepthGrid = grid->getMultiDepthHashGrid()->getGridForDepth(Depth(inneriter.getPoint().nextRight(d,t+1)));
368 do {
369
370 IndexDimension Index = inneriter.getPoint();
371 unsigned long i = inneriter.geti();
372 value = 0.5 * vec.getValue(i);
373
374 if (t >= 0) {
375
376 IndexDimension rightIndex = Index.nextRight(d);
377 if (!Index.isAtRightBoundary(d)){
378 // vec.getValue(rightIndex,value,singleDepthHashGrid)
379 vec.addToValue(rightIndex,value);
380 // value += vec.getValue(rightIndex);
381 }
382
383
384 IndexDimension leftIndex = Index.nextLeft(d);
385 if (!Index.isAtLeftBoundary(d)){
386 // vec.getValue(rightIndex,value,singleDepthHashGrid)
387 vec.addToValue(leftIndex,value);
388 // value += vec.getValue(rightIndex);
389 }
390 }
391 } while (inneriter.next());
392 }
393 } while (iter.next());
394
395 };
396
397
398 inline void DualSpaceTransformationNeumann(VectorSparseG &Ax_hier, VectorSparseG &Ax, Depth &T) {
399 SingleDepthHashGrid& depthGrid = Ax_hier.getSparseGrid()->getMultiDepthHashGrid()->getGridForDepth(T);
400 const auto& mapping = depthGrid._mapPosToGridPos;
401
402 for (size_t i = 0; i < mapping.size(); i++)
403 {
404
405 IndexDimension I = depthGrid._map.getIndexOfTable(i);
406 double coeff = 0.0;
407 IndexDimension J;
408 if(Ax_hier.getSparseGrid()->getActiveTable()[mapping[i]]) {
409 for (MultiDimFiveCompass mc; mc.goon(); ++mc) {
410 double basis_coeff = 1.0;
411 J = I.nextFive_Neumann(&mc, T, &basis_coeff);
412 if (mc.goon())
413 coeff = coeff + Ax_hier.getValue(J) * basis_coeff;
414 }
415
416 Ax.setValue(mapping[i], coeff);
417 }
418 }
419
420 }
421
422 inline void DualSpaceTransformationDirichlet(VectorSparseG &g, VectorSparseG &Ax, Depth &T) {
423 SingleDepthHashGrid& depthGrid = g.getSparseGrid()->getMultiDepthHashGrid()->getGridForDepth(T);
424 const auto& mapping = depthGrid._mapPosToGridPos;
425
426 for (size_t i = 0; i < mapping.size(); i++)
427 {
428
429 IndexDimension I = depthGrid._map.getIndexOfTable(i);
430 Richtung richtung[DimensionSparseGrid];
431 if(I.isRandnah(richtung)){
432
433 }
434 double coeff = 0.0;
435 IndexDimension J;
436 if(g.getSparseGrid()->getActiveTable()[mapping[i]]) {
437 for (MultiDimCompass mc; mc.goon(); ++mc) {
438 double basis_coeff = 1.0;
439 J = I.nextThree_boundary(&mc, T);
440 if (mc.goon())
441 coeff = coeff + g.getValue(J) * basis_coeff;
442 }
443
444 Ax.setValue(mapping[i], coeff);
445 }
446 }
447
448 }
449
450 class iterator{
451 public:
452 iterator(Richtung* richtung_){
453 maxShift = 1;
454 for(int d=0; d<DimensionSparseGrid; d++){
455 richtung[d]=richtung_[d];
456 if(richtung[d]!=Mitte)
457 maxShift *=2;
458
459 }
460 }
461 bool goon(){
462 return false;
463 }
464 private:
465 Richtung richtung[DimensionSparseGrid];
466 int maxShift;
467 int shiftnumber;
468 };
469
470};
471
473
474template<class Problem>
475void MatrixVectorInhomogen::multiplication(VectorSparseG &prew, VectorSparseG &Ax, Problem &matrixProblem) {
476
477 Ax = 0.0;
478 nodal = 0.0;
479
480 AdaptiveSparseGrid_Base *grid = prew.getSparseGrid();
481
482 int k = 1;
483 grid->WorkOnHangingNodes = true;
484
485 ListOfDepthOrderedSubgrids::iterator iter(list);
486 iter.gotoBegin();
487 do {
488 Depth T = iter.getDepth();
489 u = 0.0;
490 calcNodal(prew, T);
491 nodal.setMultiLevelValues(u, T);
492 } while (iter.next());
493
494 grid->WorkOnHangingNodes = false;
495
496
497 for (CasesIterator iter; iter.goon(); ++iter) {
498 Ax_neu = 0.0;
499 bool *restrictions = iter.getcase();
500
501 grid->WorkOnHangingNodes = true;
502
503 Ax_neu = 0.0;
504 gM = 0.0;
505
506
507 CaseFunction(prew, restrictions, matrixProblem);
508
509
510
511 grid->WorkOnHangingNodes = false;
512 Ax = Ax + Ax_neu;
513
514 k++;
515 }
516
517}
518
519
520template<class Problem>
521void MatrixVectorInhomogen::multiplicationRHSHomogen(VectorSparseG &prew, VectorSparseG &Ax, Problem &matrixProblem) {
522
523 Ax = 0.0;
524 nodal = 0.0;
525 AdaptiveSparseGrid *grid = prew.getSparseGrid();
526
527 int k = 1;
528 grid->WorkOnHangingNodes = true;
529
530
531 ListOfDepthOrderedSubgrids::iterator iter(list);
532 iter.gotoBegin();
533 do {
534 Depth T = iter.getDepth();
535 u = 0.0;
536 calcNodal(prew, T);
537 nodal.setMultiLevelValues2(u, T);
538 } while (iter.next());
539
540 grid->WorkOnHangingNodes = false;
541
542
543 for (CasesIterator iter; iter.goon(); ++iter) {
544 Ax_neu = 0.0;
545 bool *restrictions = iter.getcase();
546
547 grid->WorkOnHangingNodes = true;
548
549 u = 0.0;
550 g = 0.0;
551 z = 0.0;
552 Ax_neu = 0.0;
553 gM = 0.0;
554
555
556 CaseFunctionTransformation(prew, restrictions, matrixProblem);
557
558 grid->WorkOnHangingNodes = false;
559 Ax = Ax + Ax_neu;
560
561 k++;
562 }
563
564}
565
566template<class Problem>
567void MatrixVectorInhomogen::CaseFunction(VectorSparseG &prew, bool *restrictions, Problem &matrixProblem) {
568
569 list.SortiereTiefenBoundary(restrictions);
570 ListOfDepthOrderedSubgrids::iterator iter(list);
571
572 std::list<Depth>* sortedDepths = list.getSortierteTiefen(); // TODO maybe not deference to prevent useless copy of all elements
573
574 if (sortedDepths->size() == 0) return;
575
576 P+=nodal;
577
578 AdaptiveSparseGrid_Base *grid = prew.getSparseGrid();
579
580 for(int d=0; d < DimensionSparseGrid; d++) {
581 if(restrictions[d]==0)
582 for (Depth T: *sortedDepths) {
583 int t = T.at(d)+1;
584 Depth T_fine = T;
585 T_fine.set(t,d);
586
587 if(list.isIncluded(T_fine)) {
588 u_new = 0.0;
589 u_new.setMultiLevelValues2(P, T);
590
591 prolongation1D_inplace(u_new, t, d);
592 P.addMultiLevelValues(u_new, T_fine);
593 }
594
595 }
596 }
597
598
599 for (Depth T: *sortedDepths) {
600 grid->WorkOnHangingNodes = true;
601 u = 0.0;
602
603
604 u.setMultiLevelValues2(P,T);
605
606
607
608 z = 0.0;
609
610 matrixProblem.initialize(T);
611 applyProblemGlobal(u, z, matrixProblem, T);
612
613 g = 0.0;
614 g.setMultiLevelValues(gM, T);
615
616 z = z + g;
617 g = z;
618
619
620
621
622
623
624 Depth Tcoarse(T);
625 for (int d = 0; d < DimensionSparseGrid; d++) {
626 if (restrictions[d]) {
627
628
629 Tcoarse.set(T.at(d) - 1, d);
630 if (list.isIncluded(Tcoarse)) {
631 g = 0.0;
632
633
634 restriction1D_inverted_inplace(z, g, Tcoarse.at(d), d);
635
636
637 z = 0.0;
638 z.setMultiLevelValues(gM, Tcoarse);
639 z = z + g;
640 gM.setMultiLevelValues(g, Tcoarse);
641 }
642
643
644 }
645
646
647 }
648
649
650 if (list.isIncluded(Tcoarse)) {
651 grid->WorkOnHangingNodes = false;
652
653 DualSpaceTransformationNeumann(g, Ax_neu, Tcoarse);
654 }
655
656 }
657
658};
659
660
661template<class Problem>
662void MatrixVectorInhomogen::CaseFunctionTransformation(VectorSparseG &prew, bool *restrictions, Problem &matrixProblem) {
663
664 list.SortiereTiefenBoundary(restrictions);
665 ListOfDepthOrderedSubgrids::iterator iter(list);
666
667 std::list<Depth>* sortedDepths = list.getSortierteTiefen(); // TODO maybe not deference to prevent useless copy of all elements
668
669 if (sortedDepths->size() == 0) return;
670
671 P+=nodal;
672
673 AdaptiveSparseGrid *grid = prew.getSparseGrid();
674 for (Depth T: *sortedDepths) {
675
676 for (int d = 0; d < DimensionSparseGrid; d++) {
677 if (restrictions[d] == 0) {
678 int t = T.at(d)+1;
679 Depth T_fine = T;
680 T_fine.set(t,d);
681 }
682 }
683 }
684
685 for(int d=0; d < DimensionSparseGrid; d++) {
686 if(restrictions[d]==0)
687 for (Depth T: *sortedDepths) {
688 int t = T.at(d)+1;
689 Depth T_fine = T;
690 T_fine.set(t,d);
691
692 if(list.isIncluded(T_fine)) {
693 u_new = 0.0;
694 u_new.setMultiLevelValues2(P, T);
695
696 prolongation1D_inplace(u_new, t, d);
697 P.addMultiLevelValues(u_new, T_fine);
698 }
699
700 }
701 }
702
703
704 for (Depth T: *sortedDepths) {
705 grid->WorkOnHangingNodes = true;
706 u = 0.0;
707
708
709 u.setMultiLevelValues2(P,T);
710
711
712
713 z = 0.0;
714
715 matrixProblem.initialize(T);
716 applyProblemGlobal(u, z, matrixProblem, T);
717
718 z.addMultiLevelValues2(gM,T);
719
720
721
722
723
724
725 Depth Tcoarse(T);
726 for (int d = 0; d < DimensionSparseGrid; d++) {
727 if (restrictions[d]) {
728
729 Tcoarse.set(T.at(d) - 1, d);
730 if (list.isIncluded(Tcoarse)) {
731
732 restriction1D_inverted_inplace(z,Tcoarse.at(d), d);
733 g=z;
734 z.addMultiLevelValues2(gM,Tcoarse);
735 gM.setMultiLevelValues2(g, Tcoarse);
736
737
738 }
739
740 }
741 }
742
743
744 if (list.isIncluded(Tcoarse)) {
745 grid->WorkOnHangingNodes = false;
746 g = 0.0;
747 DualSpaceTransformationNeumann(z, g, Tcoarse);
748 DualSpaceTransformationDirichlet(g,Ax_neu,Tcoarse);
749 }
750
751 }
752
753};
754
755
756void MatrixVectorInhomogen::calcNodal(VectorSparseG &prew, Depth &Tiefe) {
757 SubgridFixedDepth::iterator iter(*list.getSubgrid(Tiefe));
758 iter.gotobegin();
759 do {
760
761 IndexDimension Index = iter.getPoint();
762 double coeff = prew.getValue(Index);
763 // Depth Tneu(Index);
764 double basis_coeff = 1.0;
765 for (MultiDimFiveCompass mc; mc.hasNext(); ++mc) {
766
767 basis_coeff = 1.0;
768 IndexDimension J = Index.nextFive_Neumann(&mc, Tiefe.returnTiefen(), &basis_coeff);
769
770 double val = u.getValue(J) + coeff * basis_coeff;
771
772 //hier.setValue(J, val);
773 //u = val | J;
774 u.setValue(J, val);
775
776
777 }
778
779
780 } while (iter.next());
781}
782
783#endif //SGRUN_MATRIXVECTORINHOMOGEN_H_H
AdaptiveSparseGrid_Base
Definition sparseGrid.h:86
AdaptiveSparseGrid_Base::getIndexOfTable
IndexDimension getIndexOfTable(unsigned long i)
Definition sparseGrid.h:566
AdaptiveSparseGrid
Definition sparseGrid.h:277
ListOfDepthOrderedSubgrids
Definition ListOfDepthOrderedGrids.h:115
MultiDimCompass
Definition index.h:356
MultiDimFiveCompass
Definition index.h:436
SubgridFixedDepth::iterator
Definition ListOfDepthOrderedGrids.h:82