Peano 4
Loading...
Searching...
No Matches
ElementMatrix.cpp
Go to the documentation of this file.
1#include "ElementMatrix.h"
2#include "StencilFactory.h"
3
5#include "peano4/utils/Loop.h"
6
7
8
17
18 dfor2(j)
19 dfor2(i)
21 double commonFacesPowerTwo = 1.0;
22 for (int d=0; d<Dimensions; d++) {
23 stencilEntry(d) = i(d)-j(d)+1;
24 if (i(d)==j(d)) commonFacesPowerTwo *= 2.0;
25 }
26 result(jScalar,iScalar) = stencil(peano4::utils::dLinearised(stencilEntry,3)) / commonFacesPowerTwo;
29
30 return result;
31}
32
33
42
43 dfor2(j)
44 dfor2(i)
46 double commonFacesPowerTwo = 1.0;
47 for (int d=0; d<Dimensions; d++) {
48 stencilEntry(d) = i(d)-j(d)+1;
49 if (i(d)==j(d)) commonFacesPowerTwo *= 2.0;
50 }
51 result(jScalar,iScalar) = complexStencil(peano4::utils::dLinearised(stencilEntry,3)) / commonFacesPowerTwo;
54
55 return result;
56}
57
58
62
63
65 return stencil(ThreePowerD / 2);
66}
67
68
72
73 dfor2(j)
74 dfor2(k)
75 const int stencilEntry = mapElementMatrixEntryOntoStencilEntry(jScalar,kScalar);
76 result(stencilEntry) += matrix(jScalar,kScalar);
79
80 return result;
81}
82
83
86
87 dfor2(j)
88 dfor2(k)
89 const int stencilEntry = mapElementMatrixEntryOntoStencilEntry(jScalar,kScalar);
90 assertionEquals1(result(jScalar*ThreePowerD + stencilEntry),0.0,result);
91 result(jScalar*ThreePowerD + stencilEntry) = matrix(jScalar,kScalar);
94
95 return result;
96}
97
98
102
103 dfor2(j)
104 dfor2(i)
105 const int stencilOffset = jScalar * ThreePowerD;
107 double commonFacesPowerTwo = 1.0;
108 for (int d=0; d<Dimensions; d++) {
109 if (i(d)==j(d)) {
110 stencilEntry(d) = 1;
111 commonFacesPowerTwo *= 2.0;
112 }
113 else if (i(d)<j(d)) {
114 stencilEntry(d) = 0;
115 }
116 else {
117 stencilEntry(d) = 2;
118 }
119 }
120 const int vectorOfStencilIndex = stencilOffset+peano4::utils::dLinearised(stencilEntry,3);
121 result(jScalar,iScalar) = vectorOfStencils( vectorOfStencilIndex ) / commonFacesPowerTwo;
124
125 #ifdef Dim2
126 assertionNumericalEquals2( result(0,0), vectorOfStencils(0*ThreePowerD+4)/4.0, result, vectorOfStencils );
127 assertionNumericalEquals2( result(0,1), vectorOfStencils(0*ThreePowerD+5)/2.0, result, vectorOfStencils );
128 assertionNumericalEquals2( result(0,2), vectorOfStencils(0*ThreePowerD+7)/2.0, result, vectorOfStencils );
129 assertionNumericalEquals2( result(0,3), vectorOfStencils(0*ThreePowerD+8)/1.0, result, vectorOfStencils );
130
131 assertionNumericalEquals2( result(1,0), vectorOfStencils(1*ThreePowerD+3)/2.0, result, vectorOfStencils );
132 assertionNumericalEquals2( result(1,1), vectorOfStencils(1*ThreePowerD+4)/4.0, result, vectorOfStencils );
133 assertionNumericalEquals2( result(1,2), vectorOfStencils(1*ThreePowerD+6)/1.0, result, vectorOfStencils );
134 assertionNumericalEquals2( result(1,3), vectorOfStencils(1*ThreePowerD+7)/2.0, result, vectorOfStencils );
135
136 assertionNumericalEquals2( result(2,0), vectorOfStencils(2*ThreePowerD+1)/2.0, result, vectorOfStencils );
137 assertionNumericalEquals2( result(2,1), vectorOfStencils(2*ThreePowerD+2)/1.0, result, vectorOfStencils );
138 assertionNumericalEquals2( result(2,2), vectorOfStencils(2*ThreePowerD+4)/4.0, result, vectorOfStencils );
139 assertionNumericalEquals2( result(2,3), vectorOfStencils(2*ThreePowerD+5)/2.0, result, vectorOfStencils );
140
141 assertionNumericalEquals2( result(3,0), vectorOfStencils(3*ThreePowerD+0)/1.0, result, vectorOfStencils );
142 assertionNumericalEquals2( result(3,1), vectorOfStencils(3*ThreePowerD+1)/2.0, result, vectorOfStencils );
143 assertionNumericalEquals2( result(3,2), vectorOfStencils(3*ThreePowerD+3)/2.0, result, vectorOfStencils );
144 assertionNumericalEquals2( result(3,3), vectorOfStencils(3*ThreePowerD+4)/4.0, result, vectorOfStencils );
145 #endif
146
147 return result;
148}
149
150
152 assertion2( cellIsInside[ThreePowerD/2], vectorOfStencils, cellIsInside );
153
155
156 dfor2(j)
157 dfor2(i)
158 const int stencilOffset = jScalar * ThreePowerD;
160 int adjacentInnerCells = 0;
161
162 dfor2(jj)
163 dfor2(ii)
164 if (jj+j==ii+i) {
165 adjacentInnerCells += ( cellIsInside[peano4::utils::dLinearised(jj+j,3)] ) ? 1 : 0;
166 }
169
170 for (int d=0; d<Dimensions; d++) {
171 if (i(d)==j(d)) {
172 stencilEntry(d) = 1;
173 }
174 else if (i(d)<j(d)) {
175 stencilEntry(d) = 0;
176 }
177 else {
178 stencilEntry(d) = 2;
179 }
180 }
181 const int vectorOfStencilIndex = stencilOffset+peano4::utils::dLinearised(stencilEntry,3);
182 result(jScalar,iScalar) = vectorOfStencils( vectorOfStencilIndex ) / static_cast<double>(adjacentInnerCells);
185
186 assertion4( !cellIsInside.all() || tarch::la::equals(result,getElementWiseAssemblyMatrix(vectorOfStencils) ), vectorOfStencils, cellIsInside, result, getElementWiseAssemblyMatrix(vectorOfStencils) );
187
188 return result;
189}
190
191
196
197
201) {
203
204 for (int d=0; d<Dimensions; d++) {
205 if (col(d)>row(d)) {
206 stencilEntry(d)=2;
207 }
208 else if (col(d)<row(d)) {
209 stencilEntry(d)=0;
210 }
211 else {
212 stencilEntry(d)=1;
213 }
214 }
215
216 return getStencilEntryLinearisedIndex(stencilEntry);
217}
218
219
221 assertion( row>=0 );
222 assertion( row<TwoPowerD );
223 assertion( col>=0 );
224 assertion( col<TwoPowerD );
225 #if defined(Dim2)
226 static const int result[] = {4,5,7,8,3,4,6,7,1,2,4,5,0,1,3,4};
227 int thisResult = result[row*4+col];
228
229 #ifdef Asserts
230 tarch::la::Vector<Dimensions,int> toIndex = peano::utils::dDelinearised(row,2);
231 tarch::la::Vector<Dimensions,int> fromIndex = peano::utils::dDelinearised(col,2);
232 assertionEquals2( thisResult, mapElementMatrixEntryOntoStencilEntry(toIndex,fromIndex), row, col );
233 #endif
234
235 return thisResult;
236 #else
239 return mapElementMatrixEntryOntoStencilEntry(toIndex,fromIndex);
240 #endif
241}
242
243
245toolbox::finiteelements::getElementWiseAssemblyMatrix( const VectorOfComplexStencils& vectorOfComplexStencils ) {
246 ComplexElementWiseAssemblyMatrix result;
247
248 dfor2(j)
249 dfor2(i)
250 const int stencilOffset = jScalar * ThreePowerD;
252 double commonFacesPowerTwo = 1.0;
253 for (int d=0; d<Dimensions; d++) {
254 if (i(d)==j(d)) {
255 stencilEntry(d) = 1;
256 commonFacesPowerTwo *= 2.0;
257 }
258 else if (i(d)<j(d)) {
259 stencilEntry(d) = 0;
260 }
261 else {
262 stencilEntry(d) = 2;
263 }
264 }
265 const int vectorOfComplexStencilIndex = stencilOffset+peano4::utils::dLinearised(stencilEntry,3);
266 result(jScalar,iScalar) = vectorOfComplexStencils( vectorOfComplexStencilIndex ) / commonFacesPowerTwo;
269
270 #ifdef Dim2
271 assertionNumericalEquals2( result(0,0), vectorOfComplexStencils(0*ThreePowerD+4)/4.0, result, vectorOfComplexStencils );
272 assertionNumericalEquals2( result(0,1), vectorOfComplexStencils(0*ThreePowerD+5)/2.0, result, vectorOfComplexStencils );
273 assertionNumericalEquals2( result(0,2), vectorOfComplexStencils(0*ThreePowerD+7)/2.0, result, vectorOfComplexStencils );
274 assertionNumericalEquals2( result(0,3), vectorOfComplexStencils(0*ThreePowerD+8)/1.0, result, vectorOfComplexStencils );
275
276 assertionNumericalEquals2( result(1,0), vectorOfComplexStencils(1*ThreePowerD+3)/2.0, result, vectorOfComplexStencils );
277 assertionNumericalEquals2( result(1,1), vectorOfComplexStencils(1*ThreePowerD+4)/4.0, result, vectorOfComplexStencils );
278 assertionNumericalEquals2( result(1,2), vectorOfComplexStencils(1*ThreePowerD+6)/1.0, result, vectorOfComplexStencils );
279 assertionNumericalEquals2( result(1,3), vectorOfComplexStencils(1*ThreePowerD+7)/2.0, result, vectorOfComplexStencils );
280
281 assertionNumericalEquals2( result(2,0), vectorOfComplexStencils(2*ThreePowerD+1)/2.0, result, vectorOfComplexStencils );
282 assertionNumericalEquals2( result(2,1), vectorOfComplexStencils(2*ThreePowerD+2)/1.0, result, vectorOfComplexStencils );
283 assertionNumericalEquals2( result(2,2), vectorOfComplexStencils(2*ThreePowerD+4)/4.0, result, vectorOfComplexStencils );
284 assertionNumericalEquals2( result(2,3), vectorOfComplexStencils(2*ThreePowerD+5)/2.0, result, vectorOfComplexStencils );
285
286 assertionNumericalEquals2( result(3,0), vectorOfComplexStencils(3*ThreePowerD+0)/1.0, result, vectorOfComplexStencils );
287 assertionNumericalEquals2( result(3,1), vectorOfComplexStencils(3*ThreePowerD+1)/2.0, result, vectorOfComplexStencils );
288 assertionNumericalEquals2( result(3,2), vectorOfComplexStencils(3*ThreePowerD+3)/2.0, result, vectorOfComplexStencils );
289 assertionNumericalEquals2( result(3,3), vectorOfComplexStencils(3*ThreePowerD+4)/4.0, result, vectorOfComplexStencils );
290/*
291 result = stencil0[4]/4.0, stencil0[5]/2.0, stencil0[7]/2.0, stencil0[8]/1.0,
292 stencil1[3]/2.0, stencil1[4]/4.0, stencil1[6]/1.0, stencil1[7]/2.0,
293 stencil2[1]/2.0, stencil2[2]/1.0, stencil2[4]/4.0, stencil2[5]/2.0,
294 stencil3[0]/1.0, stencil3[1]/2.0, stencil3[3]/2.0, stencil3[4]/4.0;
295*/
296 #endif
297
298 return result;
299}
300
301
305 double weight
306) {
307 static toolbox::finiteelements::ElementWiseAssemblyMatrix referenceStiffnessMatrix =
310 );
311
313 for (int row=0; row<TwoPowerD; row++)
314 for (int col=0; col<TwoPowerD; col++) {
315 result(row,col) = matrix(row,col) - weight * referenceStiffnessMatrix(row,col);
316 }
317
318 return result;
319}
320
321
#define assertion2(expr, param0, param1)
#define assertion4(expr, param0, param1, param2, param3)
#define assertionEquals1(lhs, rhs, a)
#define assertionNumericalEquals2(lhs, rhs, a, b)
#define assertionMsg(expr, message)
#define assertion(expr)
#define assertionEquals2(lhs, rhs, a, b)
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ our matrix elements will nabla phi_i dx f By this will be a sparse matrix
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
#define ThreePowerD
Definition Globals.h:24
#define TwoPowerD
Definition Globals.h:19
#define dfor2(counter)
Shortcut For dfor(counter,2)
Definition Loop.h:906
#define enddforx
I prefer to use this macro for dforx instead of a closing bracket as many syntax parser fail otherwis...
Definition Loop.h:928
Static (i.e.
Definition Matrix.h:31
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
Definition Loop.cpp:106
tarch::la::Vector< Dimensions, int > dDelinearised(int value, int max)
Counterpart of dLinearised().
Definition Loop.cpp:146
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.
int getStencilEntryLinearisedIndex(const tarch::la::Vector< Dimensions, int > stencilEntry)
Definition Stencil.cpp:4
double getDiagonalElement(const ElementWiseAssemblyMatrix &matrix)
ElementWiseAssemblyMatrix inverseHierarchicalTransform(const ElementWiseAssemblyMatrix &matrix, const tarch::la::Vector< Dimensions, double > &h, double weight=1.0)
int mapElementMatrixEntryOntoStencilEntry(int row, int col)
An element-wise assembly matrix is distilled from a sequence of stencils which are in turn double a...
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.
ElementWiseAssemblyMatrix hierarchicalTransform(const ElementWiseAssemblyMatrix &matrix, const tarch::la::Vector< Dimensions, double > &h, double weight=1.0)
Stencil reconstructUniformStencilFragments(const ElementWiseAssemblyMatrix &matrix)
Reconstruct stencils from element matrices.
tarch::la::Matrix< TwoPowerD, TwoPowerD, std::complex< double > > ComplexElementWiseAssemblyMatrix
Definition Stencil.h:20
ElementWiseAssemblyMatrix getElementWiseAssemblyMatrix(const Stencil &stencil)
Derive element-wise matrix from stencil.
tarch::la::Vector< TwoPowerD *ThreePowerD, double > reconstructStencilFragments(const ElementWiseAssemblyMatrix &matrix)
Reconstruct stencils from element matrices.
Simple vector class.
Definition Vector.h:134