9 assertion2( boundaryFaceNormals!=0, stencil, boundaryFaceNormals );
10 for (
int d=0; d<Dimensions; d++) {
11 if (boundaryFaceNormals[d] || boundaryFaceNormals[d+Dimensions]) {
29 source(coord0) = destination(coord1);
30 source(coord1) = destination(coord0);
31 result( peano::utils::dLinearised(destination,3) ) = stencil( peano::utils::dLinearised(source,3) );
34 assertionMsg(
false,
"operation only supported for d=3" );
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;
54 result(0) = 1.0 / 6.0;
55 result(1) = 4.0 / 6.0;
56 result(2) = 1.0 / 6.0;
117 result(0) = 1.0 / 2.0;
119 result(2) = 1.0 / 2.0;
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);
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);
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);
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);
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);
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);
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);
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);
304 result = scaling(0) * h(1)/h(0) *
309 + scaling(1) * h(0)/h(1) *
315 result = scaling(0) * h(1)*h(2)/h(0) *
321 + scaling(1) * h(0)*h(2)/h(1) *
327 + scaling(2) * h(0)*h(1)/h(2) *
348 result = scaling(0) * h(1)/h(0) *
353 + scaling(1) * h(0)/h(1) *
359 result = scaling(0) * h(1)*h(2)/h(0) *
365 + scaling(1) * h(0)*h(2)/h(1) *
371 + scaling(2) * h(0)*h(1)/h(2) *
378 result = scaling(0) * h(1)*h(2)*h(3)/h(0) *
385 + scaling(1) * h(0)*h(2)*h(3)/h(1) *
392 + scaling(2) * h(0)*h(1)*h(3)/h(2) *
399 + scaling(2) * h(0)*h(1)*h(2)/h(3) *
407 result = scaling(0) * h(1)*h(2)*h(3)*h(4)/h(0) *
415 + scaling(1) * h(0)*h(2)*h(3)*h(4)/h(1) *
423 + scaling(2) * h(0)*h(1)*h(3)*h(4)/h(2) *
431 + scaling(2) * h(0)*h(1)*h(2)*h(4)/h(3) *
439 + scaling(2) * h(0)*h(1)*h(2)*h(3)/h(4) *
461 result = scaling(0) * h(1)/h(0) *
466 + scaling(1) * h(0)/h(1) *
472 result = scaling(0) * h(1)*h(2)/h(0) *
478 + scaling(1) * h(0)*h(2)/h(1) *
484 + scaling(2) * h(0)*h(1)/h(2) *
506 result(1) = - convCoeff(1) - fabs(convCoeff(1));
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));
512 result(7) = convCoeff(1) - fabs(convCoeff(1));
522 result(4) = -convCoeff(2) - fabs(convCoeff(2));
529 result(10) = - convCoeff(1) - fabs(convCoeff(1));
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));
535 result(16) = convCoeff(1) - fabs(convCoeff(1));
542 result(22) = convCoeff(2) - fabs(convCoeff(2));
549 result *= h(0)*h(1)*0.5;
561 const std::bitset<Dimensions*2>& boundaryFaceNormals,
568 if ( boundaryFaceNormals[0] ) {
569 stencil(4) += stencil(0);
571 stencil(4) += stencil(3);
573 stencil(4) += stencil(6);
577 if ( boundaryFaceNormals[1] ) {
578 stencil(4) += stencil(0);
580 stencil(4) += stencil(1);
582 stencil(4) += stencil(2);
586 if ( boundaryFaceNormals[2] ) {
587 stencil(4) += stencil(2);
589 stencil(4) += stencil(5);
591 stencil(4) += stencil(8);
605 if ( boundaryFaceNormals[3] ) {
606 stencil(4) += stencil(6);
608 stencil(4) += stencil(7);
610 stencil(4) += stencil(8);
631 ) * h(0) * h(1) * h(2);
654 )) * h(0) * h(1) * h(2);
665 const std::complex<double>& phi ) {
688 const std::complex<double>& phi
741 * h(0) * h(1) * h(2);
765 * h(0) * h(1) * h(2);
773 * h(0) * h(1) * h(2) * h(3);
782 * h(0) * h(1) * h(2) * h(3) * h(4);
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;
851 int integrationPointsPerAxis,
861 const double termScaling = 1.0 / std::pow(integrationPointsPerAxis,Dimensions);
863 dfor( subVolume, integrationPointsPerAxis ) {
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) );
877 result(targetVertexScalar,testVertexScalar) += cellScaling * referenceStiffnessMatrix(targetVertexScalar,testVertexScalar);
#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.
#define dfor2(counter)
Shortcut For dfor(counter,2)
#define dfore(counter, max, dim, value)
This is an exclusive d-dimensional for loop.
#define dfor(counter, max)
d-dimensional Loop
#define enddforx
I prefer to use this macro for dforx instead of a closing bracket as many syntax parser fail otherwis...
const float const float const float struct part *restrict struct part *restrict const float a
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
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)