Peano
Loading...
Searching...
No Matches
StencilFactory.cpp
Go to the documentation of this file.
1#include "StencilFactory.h"
2#include "ElementMatrix.h"
3
4#include "peano4/utils/Loop.h"
5#include "tarch/Assertions.h"
6
7
8void toolbox::finiteelements::preprocessBoundaryStencil( Stencil& stencil, const std::bitset<Dimensions*2>& boundaryFaceNormals ) {
9 assertion2( boundaryFaceNormals!=0, stencil, boundaryFaceNormals );
10 for (int d=0; d<Dimensions; d++) {
11 if (boundaryFaceNormals[d] || boundaryFaceNormals[d+Dimensions]) {
12 dfore(k,3,d,1) {
13 int entry = peano4::utils::dLinearised(k,3);
14 stencil(entry) *= 2;
15 }
16 }
17 }
18}
19
20
22toolbox::finiteelements::exchangeCoordinates( const Stencil& stencil, int coord0, int coord1 ) {
23 Stencil result;
24
25#ifdef Dim3
26 d3for3(destination)
28 source = destination;
29 source(coord0) = destination(coord1);
30 source(coord1) = destination(coord0);
31 result( peano::utils::dLinearised(destination,3) ) = stencil( peano::utils::dLinearised(source,3) );
33 #else
34 assertionMsg( false, "operation only supported for d=3" );
35#endif
36
37 return result;
38}
39
40
43 result(0) = 1.0 / 3.0;
44 result(1) = 2.0 / 3.0;
45 result(2) = 3.0 / 3.0;
46 result(3) = 2.0 / 3.0;
47 result(4) = 1.0 / 3.0;
48 return result;
49}
50
51
54 result(0) = 1.0 / 6.0;
55 result(1) = 4.0 / 6.0;
56 result(2) = 1.0 / 6.0;
57 return result;
58}
59
60
63 result(0) = 0.0;
64 result(1) = 1.0;
65 result(2) = 0.0;
66 return result;
67}
68
69
72 result(0) = 1.0;
73 result(1) = 0.0;
74 result(2) = 0.0;
75 return result;
76}
77
78
81 result(0) = 0.0;
82 result(1) = 0.0;
83 result(2) = 1.0;
84 return result;
85}
86
87
90 result(0) = -1.0;
91 result(1) = 1.0;
92 result(2) = 0.0;
93 return result;
94}
95
96
99 result(0) = 0.0;
100 result(1) = -1.0;
101 result(2) = 1.0;
102 return result;
103}
104
105
108 result(0) = -1.0;
109 result(1) = 0.0;
110 result(2) = 1.0;
111 return result;
112}
113
114
117 result(0) = 1.0 / 2.0;
118 result(1) = 0.0;
119 result(2) = 1.0 / 2.0;
120 return result;
121}
122
123
126 result(0) = -1.0;
127 result(1) = 2.0;
128 result(2) = -1.0;
129 return result;
130}
131
132
136 ) {
138
139 for (int i=0; i<3; i++) {
140 for (int j=0; j<3; j++) {
141 result(i+j*3) = a(i) * b(j);
142 }
143 }
144
145 return result;
146}
147
148
152 ) {
154
155 for (int i=0; i<5; i++) {
156 for (int j=0; j<5; j++) {
157 result(i+j*5) = a(i) * b(j);
158 }
159 }
160
161 return result;
162}
163
164
169 ) {
171
172 for (int i=0; i<3; i++) {
173 for (int j=0; j<3; j++) {
174 for (int k=0; k<3; k++) {
175 result(i+j*3+k*3*3) = a(i) * b(j) * c(k);
176 }
177 }
178 }
179
180 return result;
181}
182
183
188 ) {
190
191 for (int i=0; i<5; i++) {
192 for (int j=0; j<5; j++) {
193 for (int k=0; k<5; k++) {
194 result(i+j*5+k*5*5) = a(i) * b(j) * c(k);
195 }
196 }
197 }
198
199 return result;
200}
201
202
208) {
210
211 for (int ai=0; ai<3; ai++) {
212 for (int bi=0; bi<3; bi++) {
213 for (int ci=0; ci<3; ci++) {
214 for (int di=0; di<3; di++) {
215 result(ai+bi*3+ci*3*3+di*3*3*3) = a(ai) * b(bi) * c(ci) * d(di);
216 }
217 }
218 }
219 }
220
221 return result;
222}
223
224
230) {
232
233 for (int ai=0; ai<5; ai++) {
234 for (int bi=0; bi<5; bi++) {
235 for (int ci=0; ci<5; ci++) {
236 for (int di=0; di<5; di++) {
237 result(ai+bi*5+ci*5*5+di*5*5*5) = a(ai) * b(bi) * c(ci) * d(di);
238 }
239 }
240 }
241 }
242
243 return result;
244}
245
246
253) {
255
256 for (int ai=0; ai<3; ai++) {
257 for (int bi=0; bi<3; bi++) {
258 for (int ci=0; ci<3; ci++) {
259 for (int di=0; di<3; di++) {
260 for (int ei=0; ei<3; ei++) {
261 result(ai+bi*3+ci*3*3+di*3*3*3+ei*3*3*3*3) = a(ai) * b(bi) * c(ci) * d(di) * e(ei);
262 }
263 }
264 }
265 }
266 }
267
268 return result;
269}
270
271
278) {
280
281 for (int ai=0; ai<5; ai++) {
282 for (int bi=0; bi<5; bi++) {
283 for (int ci=0; ci<5; ci++) {
284 for (int di=0; di<5; di++) {
285 for (int ei=0; ei<5; ei++) {
286 result(ai+bi*5+ci*5*5+di*5*5*5+ei*5*5*5*5) = a(ai) * b(bi) * c(ci) * d(di) * e(ei);
287 }
288 }
289 }
290 }
291 }
292
293 return result;
294}
295
296
339
340
342 const tarch::la::Vector<Dimensions,std::complex<double> >& h,
343 const tarch::la::Vector<Dimensions,std::complex<double> >& scaling
344) {
346
347 #if Dimensions==2
348 result = scaling(0) * h(1)/h(0) *
349 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
352 ))
353 + scaling(1) * h(0)/h(1) *
354 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
357 ));
358 #elif Dimensions==3
359 result = scaling(0) * h(1)*h(2)/h(0) *
360 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
364 ))
365 + scaling(1) * h(0)*h(2)/h(1) *
366 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
370 ))
371 + scaling(2) * h(0)*h(1)/h(2) *
372 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
376 ));
377 #elif defined(Dim4)
378 result = scaling(0) * h(1)*h(2)*h(3)/h(0) *
379 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
384 ))
385 + scaling(1) * h(0)*h(2)*h(3)/h(1) *
386 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
391 ))
392 + scaling(2) * h(0)*h(1)*h(3)/h(2) *
393 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
398 ))
399 + scaling(2) * h(0)*h(1)*h(2)/h(3) *
400 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
405 ));
406 #elif defined(Dim5)
407 result = scaling(0) * h(1)*h(2)*h(3)*h(4)/h(0) *
408 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
414 ))
415 + scaling(1) * h(0)*h(2)*h(3)*h(4)/h(1) *
416 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
422 ))
423 + scaling(2) * h(0)*h(1)*h(3)*h(4)/h(2) *
424 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
430 ))
431 + scaling(2) * h(0)*h(1)*h(2)*h(4)/h(3) *
432 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
438 ))
439 + scaling(2) * h(0)*h(1)*h(2)*h(3)/h(4) *
440 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
446 ));
447 #else
448 assertionMsg( false, "dimension not supported" );
449 #endif
450
451 return result;
452}
453
456 const tarch::la::Vector<Dimensions,std::complex<double> >& scaling
457) {
459
460 #if Dimensions==2
461 result = scaling(0) * h(1)/h(0) *
462 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
465 ))
466 + scaling(1) * h(0)/h(1) *
467 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
470 ));
471 #elif Dimensions==3
472 result = scaling(0) * h(1)*h(2)/h(0) *
473 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
477 ))
478 + scaling(1) * h(0)*h(2)/h(1) *
479 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
483 ))
484 + scaling(2) * h(0)*h(1)/h(2) *
485 tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
489 ));
490 #else
491 assertionMsg( false, "dimension not supported" );
492 #endif
493
494 return result;
495}
496
497
501) {
503
504 #if Dimensions==2
505 result(0) = 0.0;
506 result(1) = - convCoeff(1) - fabs(convCoeff(1));
507 result(2) = 0.0;
508 result(3) = - convCoeff(0) - fabs(convCoeff(0));
509 result(4) = 2.0*(fabs(convCoeff(0)) + fabs(convCoeff(1)));
510 result(5) = convCoeff(0) - fabs(convCoeff(0));
511 result(6) = 0.0;
512 result(7) = convCoeff(1) - fabs(convCoeff(1));
513 result(8) = 0.0;
514
515 assertionNumericalEquals(h(0), h(1));
516 result *= h(0)*0.5;
517 #elif Dimensions==3
518 result(0) = 0.0;
519 result(1) = 0.0;
520 result(2) = 0.0;
521 result(3) = 0.0;
522 result(4) = -convCoeff(2) - fabs(convCoeff(2));
523 result(5) = 0.0;
524 result(6) = 0.0;
525 result(7) = 0.0;
526 result(8) = 0.0;
527
528 result(9) = 0.0;
529 result(10) = - convCoeff(1) - fabs(convCoeff(1));
530 result(11) = 0.0;
531 result(12) = - convCoeff(0) - fabs(convCoeff(0));
532 result(13) = 2.0*(fabs(convCoeff(0)) + fabs(convCoeff(1)) + fabs(convCoeff(2)));
533 result(14) = convCoeff(0) - fabs(convCoeff(0));
534 result(15) = 0.0;
535 result(16) = convCoeff(1) - fabs(convCoeff(1));
536 result(17) = 0.0;
537
538 result(18) = 0.0;
539 result(19) = 0.0;
540 result(20) = 0.0;
541 result(21) = 0.0;
542 result(22) = convCoeff(2) - fabs(convCoeff(2));
543 result(23) = 0.0;
544 result(24) = 0.0;
545 result(25) = 0.0;
546 result(26) = 0.0;
547
548 assertionNumericalEquals(h(0), h(1));
549 result *= h(0)*h(1)*0.5;
550 #else
551 assertionMsg( false, "dimension not supported" );
552 #endif
553
554 return result;
555}
556
557
560 double& rhs,
561 const std::bitset<Dimensions*2>& boundaryFaceNormals,
563 double derivative
564) {
565
566 #if Dimensions==2
567
568 if ( boundaryFaceNormals[0] ) { // left boundary
569 stencil(4) += stencil(0);
570 stencil(0) = 0.0;
571 stencil(4) += stencil(3);
572 stencil(3) = 0.0;
573 stencil(4) += stencil(6);
574 stencil(6) = 0.0;
575 }
576
577 if ( boundaryFaceNormals[1] ) { // bottom boundary
578 stencil(4) += stencil(0);
579 stencil(0) = 0.0;
580 stencil(4) += stencil(1);
581 stencil(1) = 0.0;
582 stencil(4) += stencil(2);
583 stencil(2) = 0.0;
584 }
585
586 if ( boundaryFaceNormals[2] ) { // right boundary
587 stencil(4) += stencil(2);
588 stencil(2) = 0.0;
589 stencil(4) += stencil(5);
590 stencil(5) = 0.0;
591 stencil(4) += stencil(8);
592 stencil(8) = 0.0;
593
594/*
595 stencil(0) += stencil(2);
596 stencil(2) = 0.0;
597 stencil(3) += stencil(5);
598 stencil(5) = 0.0;
599 stencil(6) += stencil(8);
600 stencil(8) = 0.0;
601 rhs += h(0)*derivative*(stencil(2)+stencil(5)+stencil(8) );
602*/
603 }
604
605 if ( boundaryFaceNormals[3] ) { // top boundary
606 stencil(4) += stencil(6);
607 stencil(6) = 0.0;
608 stencil(4) += stencil(7);
609 stencil(7) = 0.0;
610 stencil(4) += stencil(8);
611 stencil(8) = 0.0;
612 }
613
614 #endif
615}
616
617
638
639
642
643 #if Dimensions==2
644 result = tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
647 ))
648 * h(0) * h(1);
649 #elif Dimensions==3
650 result = tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
654 )) * h(0) * h(1) * h(2);
655 #else
656 assertionMsg( false, "dimension not supported" );
657 #endif
658
659 return result;
660}
661
662
665 const std::complex<double>& phi ) {
667
668 #if Dimensions==2
669 result = tarch::la::convertScalar< std::complex<double> > ( h(0) * h(1)* toolbox::finiteelements::stencilProduct(
672 ));
673 #elif Dimensions==3
674 result = tarch::la::convertScalar< std::complex<double> >( h(0) * h(1) * h(2) * toolbox::finiteelements::stencilProduct(
678 ));
679 #else
680 assertionMsg( false, "dimension not supported" );
681 #endif
682
683 return phi * result;
684}
685
687 const tarch::la::Vector<Dimensions,std::complex<double> >& h,
688 const std::complex<double>& phi
689) {
691
692 #if Dimensions==2
693 result = h(0) * h(1)* ( tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
696 ));
697 #elif Dimensions==3
698 result = h(0) * h(1) * h(2) * ( tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
702 ));
703 #elif defined(Dim4)
704 result = h(0) * h(1) * h(2) * h(3) * ( tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
709 ));
710 #elif defined(Dim5)
711 result = h(0) * h(1) * h(2) * h(3) * h(4) * ( tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
717 ));
718 #else
719 assertionMsg( false, "dimension not supported" );
720 #endif
721
722 return phi * result;
723}
724
725
728
729 #if Dimensions==2
733 )
734 * h(0) * h(1);
735 #elif Dimensions==3
740 )
741 * h(0) * h(1) * h(2);
742 #else
743 assertionMsg( false, "dimension not supported" );
744 #endif
745
746 return result;
747}
748
749
752
753 #if Dimensions==2
754 result = tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
757 ))
758 * h(0) * h(1);
759 #elif Dimensions==3
760 result = tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
764 ))
765 * h(0) * h(1) * h(2);
766 #elif defined(Dim4)
767 result = tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
772 ))
773 * h(0) * h(1) * h(2) * h(3);
774 #elif defined(Dim4)
775 result = tarch::la::convertScalar< std::complex<double> >( toolbox::finiteelements::stencilProduct(
781 ))
782 * h(0) * h(1) * h(2) * h(3) * h(4);
783 #else
784 assertionMsg( false, "dimension not supported" );
785 #endif
786
787 return result;
788}
789
790
829
830
833
834 dfor2(i)
836 double commonFacesPowerTwo = 1.0;
837 for (int d=0; d<Dimensions; d++) {
838 stencilEntry(d) = i(d)+adjacentCell(d);
839 if (stencilEntry(d)==1) commonFacesPowerTwo *= 2.0;
840 }
841 result(peano4::utils::dLinearised(stencilEntry,3)) = stencil(peano4::utils::dLinearised(stencilEntry,3)) / commonFacesPowerTwo;
843
844 return result;
845}
846
847
849 const tarch::la::Vector<Dimensions,double>& cellCentre,
851 int integrationPointsPerAxis,
852 std::function<double(const tarch::la::Vector<Dimensions,double>&)> epsilon
853) {
854 static toolbox::finiteelements::ElementWiseAssemblyMatrix referenceStiffnessMatrix =
857 );
858
860
861 const double termScaling = 1.0 / std::pow(integrationPointsPerAxis,Dimensions);
862
863 dfor( subVolume, integrationPointsPerAxis ) {
865 cellCentre - h/2.0 +
866 tarch::la::multiplyComponents( tarch::la::convertScalar<double>( subVolume )+0.5, h/static_cast<double>(integrationPointsPerAxis));
867
868 double epsilonScaling = epsilon(x) * tarch::la::volume(h) / h(0) / h(0);
869
870 dfor2( targetVertex )
871 double cellScaling = epsilonScaling;
872 for (int d=0; d<Dimensions; d++) {
873 int mirroredIndex = targetVertex(d)==0 ? subVolume(d) : integrationPointsPerAxis-subVolume(d)-1;
874 cellScaling *= termScaling * ( std::pow(mirroredIndex+1,Dimensions) - std::pow(mirroredIndex,Dimensions) );
875 }
876 dfor2( testVertex )
877 result(targetVertexScalar,testVertexScalar) += cellScaling * referenceStiffnessMatrix(targetVertexScalar,testVertexScalar);
880 }
881
882 return result;
883}
#define assertion2(expr, param0, param1)
#define assertionMsg(expr, message)
#define assertionNumericalEquals(lhs, rhs)
#define d3for3(counter)
If Dimensions is not set to three, we might nevertheless need three-dimensional loops.
Definition Loop.h:838
#define dfor2(counter)
Shortcut For dfor(counter,2)
Definition Loop.h:911
#define dfore(counter, max, dim, value)
This is an exclusive d-dimensional for loop.
Definition Loop.h:942
#define dfor(counter, max)
d-dimensional Loop
Definition Loop.h:313
#define enddforx
I prefer to use this macro for dforx instead of a closing bracket as many syntax parser fail otherwis...
Definition Loop.h:933
Static (i.e.
Definition Matrix.h:31
const float const float const float struct part *restrict struct part *restrict const float a
Definition hydro_iact.h:55
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
Definition Loop.cpp:106
Scalar volume(const Vector< Size, Scalar > &vector)
Computes the volume of the tetrahedron spanned by the Cartesian unit vectors scaled by the correspond...
Matrix< Rows, Cols, Scalar > multiplyComponents(const Matrix< Rows, X, Scalar > &lhs, const Matrix< X, Cols, Scalar > &rhs)
ElementWiseAssemblyMatrix getPoissonMatrixWithJumpingCoefficient(const tarch::la::Vector< Dimensions, double > &cellCentre, const tarch::la::Vector< Dimensions, double > &h, int integrationPointsPerAxis, std::function< double(const tarch::la::Vector< Dimensions, double > &)>)
Integrate over element with a subvoxel scheme (cmp.
tarch::la::Vector< 3, double > get1DIdentity()
Stencil is not scaled at all with any mesh width.
void applyNeumannBC(tarch::la::Vector< ThreePowerD, double > &stencil, double &rhs, const std::bitset< Dimensions *2 > &boundaryFaceNormals, const tarch::la::Vector< Dimensions, double > &h, double derivative)
tarch::la::Vector< FivePowerD, double > getDLinearInterpolation()
Computes the stencil for a d-linear interpolation.
tarch::la::Vector< 5, double > get1DLinearInterpolationStencil()
Stencil is not scaled at all with any mesh width.
tarch::la::Vector< ThreePowerD, double > getUpwindDiscretisedConvection(const tarch::la::Vector< Dimensions, double > &h, const tarch::la::Vector< Dimensions, double > &convCoeff)
Computes a convection operator, defined by the convection coefficients, discretised with the first or...
Stencil exchangeCoordinates(const Stencil &stencil, int coord0, int coord1)
Exchanges the coordinates of a stencil.
tarch::la::Vector< ThreePowerD, double > getMassMatrix(const tarch::la::Vector< Dimensions, double > &h)
Computes the mass matrix.
tarch::la::Vector< ThreePowerD, std::complex< double > > getHelmholtzShiftMatrix(const tarch::la::Vector< Dimensions, double > &h, const std::complex< double > &phi)
Computes the Helmholtz shift term matrix.
tarch::la::Vector< 3, double > get1DCentralDifferencesStencil()
tarch::la::Vector< ThreePowerD, double > getLaplacian(const tarch::la::Vector< Dimensions, double > &h, const tarch::la::Vector< Dimensions, double > &scaling=1.0)
Computes the Laplacian.
tarch::la::Vector< 3, double > get1DRightIdentity()
void preprocessBoundaryStencil(Stencil &stencil, const std::bitset< Dimensions *2 > &boundaryFaceNormals)
Preprocess boundary stencils held by boundary vertices.
tarch::la::Vector< 3, double > get1DMeanValueStencil()
Stencil is not scaled at all with any mesh width.
tarch::la::Vector< StencilSize *StencilSize, double > stencilProduct(const tarch::la::Vector< StencilSize, double > &a, const tarch::la::Vector< StencilSize, double > &b)
tarch::la::Vector< 3, double > get1DDownwindStencil()
Stencil is not scaled at all with any mesh width.
tarch::la::Vector< 3, double > get1DLaplaceStencil()
Stencil is not scaled at all with any mesh width.
tarch::la::Vector< 3, double > get1DUpwindStencil()
tarch::la::Vector< 3, double > get1DLeftIdentity()
tarch::la::Vector< 3, double > get1DMassStencil()
ElementWiseAssemblyMatrix getElementWiseAssemblyMatrix(const Stencil &stencil)
Derive element-wise matrix from stencil.
tarch::la::Vector< ThreePowerD, double > getIdentity(const tarch::la::Vector< Dimensions, double > &h)
Computes the mass matrix.
Stencil extractElementStencil(const Stencil &stencil, const tarch::la::Vector< Dimensions, int > &adjacentCell)
The standard stencil spans 2^d cells adjacent to a vertex.
Simple vector class.
Definition Vector.h:134