Peano 4
Loading...
Searching...
No Matches
DGUtils.cpp
Go to the documentation of this file.
2
3
5 int unknownsPlusAuxiliaryVariables,
6 int order,
7 int numberOfProjectedQuantities,
8 int normal,
9 int isRightFaceHalf,
10 const double* __restrict__ srcQ,
11 double* __restrict__ destQ
12) {
14 unknownsPlusAuxiliaryVariables,
15 order+1,
16 numberOfProjectedQuantities,
17 normal,
18 isRightFaceHalf,
19 srcQ,
20 destQ
21 );
22}
23
24
26 const tarch::la::Vector<3,double>& cellSize,
27 const tarch::la::Vector<3,int>& index,
28 const double* __restrict__ quadratureWeights
29) {
30 double result = 1.0;
31 for (int d=0; d<Dimensions; d++) {
32 result *= cellSize(d) * quadratureWeights[ index(d) ];
33 }
34 return result;
35}
36
37
39 const tarch::la::Vector<2,double>& cellCentre,
40 const tarch::la::Vector<2,double>& cellSize,
41 const tarch::la::Vector<2,int>& index,
42 int polynomialOrder,
43 const double* __restrict__ quadraturePoints
44) {
46
47 result(0) = quadraturePoints[index(0)];
48 result(1) = quadraturePoints[index(1)];
49
50 return tarch::la::multiplyComponents(result,cellSize) + cellCentre - 0.5*cellSize;
51}
52
53
55 const tarch::la::Vector<3,double>& cellCentre,
56 const tarch::la::Vector<3,double>& cellSize,
57 const tarch::la::Vector<3,int>& index,
58 int polynomialOrder,
59 const double* __restrict__ quadraturePoints
60) {
62
63 result(0) = quadraturePoints[index(0)];
64 result(1) = quadraturePoints[index(1)];
65 result(2) = quadraturePoints[index(2)];
66
67 return tarch::la::multiplyComponents(result,cellSize) + cellCentre - 0.5*cellSize;
68}
69
70
71int exahype2::dg::getNodesPerCell(int nodesPerAxis){
72 int nodesPerCell=nodesPerAxis;
73 for(int d=1; d<Dimensions; d++) {
74 nodesPerCell *= nodesPerAxis;
75 }
76 return nodesPerCell;
77}
78
81 for(int i=1; i<Dimensions; i++){
82 strides[i] = strides[i-1]*nodesPerAxis;
83 }
84 return strides;
85}
86
88 int node,
90){
91
93 int tmp = 0;
94 for(int d=Dimensions-1; d>=0; d--){
95 index[d] = (node - tmp) / strides[d];
96 tmp += strides[d]*index[d];
97 }
98 return index;
99
100}
101
102
104 const tarch::la::Vector<Dimensions,int>& indexCell,
105 const int direction,
106 const int orientation,
107 const int nodesPerAxis
108 ){
109
110 /*
111 * Imagine a 3d cell of order 2, so there are 9 nodes
112 * for a node [i,j,k], we want face element [j,k] if in
113 * direction i, [i,k] in direction j and [i,j] if in
114 * direction k.
115 * Element [a,b] of a face is element a*strides[0]+b
116 * (concept is same in other dimensions but with
117 * additional terms and greater strides)
118 *
119 * However since the Hull is one array, we need to skip
120 * the previous faces, so ignore nodesPerFace*dimensions
121 * indexes (nodesPerFace*dimensions+1) if right.
122 */
123
124 int HullIndex = 0;
125 int stride = 1;
126
127 for(int dim=0; dim<Dimensions; dim++){
128 if(dim!=direction){
129 HullIndex += indexCell[dim]*stride;
130 stride *= nodesPerAxis;
131 }
132 }
133
134 /* if left face need to skip left and right faces for each previous dimension
135 * if right face, need to skip those plus the left face
136 */
137 return HullIndex + stride*(2*direction+orientation); //stride is nodesPerAxis^(d)
138
139}
140
141
143 const double* __restrict__ const QCell,
144 const double* __restrict__ const derivativeOperator,
145 const double invDx,
146 const int nodesPerAxis,
147 const int strideQ,
148 const int scalarIndex,
149 double* __restrict__ gradQ
150 ){
151
152 tarch::la::Vector<Dimensions,int> strides = getStrides(nodesPerAxis);
153 tarch::la::Vector<Dimensions,int> index = getIndex(scalarIndex, strides);
154
155 for(int dim=0; dim<Dimensions; dim++){
156
157 for(int var=0; var<strideQ; var++){
158 gradQ[dim*strideQ+var] = 0.0;
159 }//var
160
161 for(int node=0; node<nodesPerAxis; node++){
162// const double coeff = invDx * derivativeOperator[node+nodesPerAxis*index[dim]];
163 const double coeff = invDx * derivativeOperator[node+nodesPerAxis*index[dim]];
164 for(int var=0; var<strideQ; var++){
165 gradQ[dim*strideQ+var] += coeff*QCell[(scalarIndex + (node-index[dim])*strides[dim])*strideQ + var];
166 }//var
167 }//node
168 }//dim
169
170}
171
172
174 double* __restrict__ QOut,
175 const double* __restrict__ Qsubstract,
176 const int order,
177 const int unknowns,
178 const int auxiliaryVariables
179) {
180 const int nodesPerAxis = order+1;
181 const int nodesPerCell = getNodesPerCell(nodesPerAxis);
182 const int strideQ = unknowns+auxiliaryVariables;
183
184 for(int node=0; node<nodesPerCell; node++){
185 for(int var=0; var<strideQ; var++){
186 QOut[node*strideQ+var] -= Qsubstract[node*strideQ+var];
187 }
188 }
189}
190
191
193 const double* __restrict__ Q,
194 const int order,
195 const int unknowns,
196 const int auxiliaryVariables
197) {
198 return ::exahype2::fv::plotPatch(
199 Q,
200 order+1,
201 unknowns,
202 auxiliaryVariables,
203 Dimensions==2
204 );
205}
206
207
209 const double* __restrict__ Q,
210 const int order,
211 const int unknowns,
212 const int auxiliaryVariables,
213 int normal,
214 int numberOfQuantitiesProjectedOntoFace
215) {
216 return ::exahype2::fv::plotPatchOverlap(
217 Q,
218 unknowns,
219 auxiliaryVariables,
220 order+1,
221 numberOfQuantitiesProjectedOntoFace,
222 normal,
223 Dimensions==2
224 );
225}
std::string plotFace(const double *__restrict__ Q, const int order, const int unknowns, const int auxiliaryVariables, int normal, int numberOfQuantitiesProjectedOntoFace)
Definition DGUtils.cpp:208
tarch::la::Vector< Dimensions, int > getStrides(int nodesPerAxis)
Definition DGUtils.cpp:79
void computeGradient(const double *__restrict__ const QCell, const double *__restrict__ const derivativeOperator, const double invDx, const int nodesPerAxis, const int strideQ, const int scalarIndex, double *__restrict__ gradQ)
Definition DGUtils.cpp:142
double getQuadratureWeight(const tarch::la::Vector< 3, double > &cellSize, const tarch::la::Vector< 3, int > &index, const double *__restrict__ quadratureWeights)
Compute integral over shape function over cell defined by index.
Definition DGUtils.cpp:25
tarch::la::Vector< Dimensions, int > getIndex(int node, tarch::la::Vector< Dimensions, int > strides)
Definition DGUtils.cpp:87
tarch::la::Vector< 2, double > getQuadraturePoint(const tarch::la::Vector< 2, double > &cellCentre, const tarch::la::Vector< 2, double > &cellSize, const tarch::la::Vector< 2, int > &index, int polynomialOrder, const double *__restrict__ quadraturePoints)
Construct location of a quadrature point.
Definition DGUtils.cpp:38
int getNodesPerCell(int nodesPerAxis)
The number of nodes in a cell is basically the input to the power of d.
Definition DGUtils.cpp:71
void subtractCell(double *__restrict__ QOut, const double *__restrict__ Qsubstract, const int order, const int unknowns, const int auxiliaryVariables)
Definition DGUtils.cpp:173
void copyOneSideOfFaceProjection(int unknownsPlusAuxiliaryVariables, int order, int numberOfProjectedQuantities, int normal, int isRightFaceHalf, const double *__restrict__ srcQ, double *__restrict__ destQ)
Delegate to PatchUtils of the Finite Volume scheme.
Definition DGUtils.cpp:4
std::string plotCell(const double *__restrict__ Q, const int order, const int unknowns, const int auxiliaryVariables)
Definition DGUtils.cpp:192
int cellIndexToHullIndex(const tarch::la::Vector< Dimensions, int > &indexCell, const int direction, const int orientation, const int nodesPerAxis)
Definition DGUtils.cpp:103
void copyHalfOfHalo(int unknownsPlusAuxiliaryVariables, int numberOfGridCellsPerPatchPerAxis, int haloSize, int normal, bool isRightLayer, const double *__restrict__ srcQ, double *__restrict__ destQ)
A face always holds a left and a right overlap.
Matrix< Rows, Cols, Scalar > multiplyComponents(const Matrix< Rows, X, Scalar > &lhs, const Matrix< X, Cols, Scalar > &rhs)
Simple vector class.
Definition Vector.h:134