11 const int auxiliaryVariables,
13 NonConservativeProduct nonconservativeProduct,
15 PointSources pointSources,
16 const double* __restrict__ quadratureNodes,
17 const double* __restrict__ massMatrixDiagonal,
18 const double* __restrict__ stiffnessMatrix,
19 const double* __restrict__ derivativeOperator,
21 bool evaluateNonconservativeProduct,
23 bool evaluatePointSources
41 double* fluxValues = evaluateFlux ?
new double[unknowns] :
nullptr;
42 double* sourceValues = evaluateSource ?
new double[unknowns] :
nullptr;
43 double* gradientValues = evaluateNonconservativeProduct ?
new double[(unknowns+auxiliaryVariables)] :
nullptr;
44 double* ncpValues = evaluateNonconservativeProduct ?
new double[unknowns] :
nullptr;
46 for (
int currentCell=0; currentCell<cellData.
numberOfCells; currentCell++) {
47 const double* __restrict__ cellQin = cellData.
QIn[currentCell];
48 double* __restrict__ cellQout = cellData.
QOut[currentCell];
51 const double t = cellData.
t[currentCell];
52 const double dt = cellData.
dt[currentCell];
54 const int nodesPerAxis = order+1;
55 const int strideQ=unknowns+auxiliaryVariables;
58 dfor( dof, order+1 ) {
59 for(
int var=0; var<unknowns; var++){
60 cellQout[enumeratorOut(currentCell,dof,var)] = 0.0;
65 for(
int dim=0; dim<Dimensions; dim++)
66 for(
int i=0; i<nodesPerAxis; i++){
69 contributingNodeIndex[dim] = i;
73 contributingNodeIndex,
77 flux( &cellQin[enumeratorIn(currentCell, contributingNodeIndex, 0)],
89 double coeff = stiffnessMatrix[i+dof(dim)*nodesPerAxis];
90 for (
int dd=0; dd<Dimensions; dd++) {
91 coeff *= (dd==dim) ? 1.0 : massMatrixDiagonal[ dof(dd) ] * h(dd);
94 for(
int var=0; var<unknowns; var++){
95 cellQout[enumeratorOut(currentCell, dof, var)] += coeff * fluxValues[var];
112 &cellQin[enumeratorIn(currentCell, dof, 0)],
120 double coeff = massMatrixDiagonal[ dof(0) ] * massMatrixDiagonal[ dof(1) ] *
tarch::la::volume(h);
122 double coeff = massMatrixDiagonal[ dof(0) ] * massMatrixDiagonal[ dof(1) ] * massMatrixDiagonal[ dof(2) ] *
tarch::la::volume(h);
125 for(
int var=0; var<unknowns; var++){
126 cellQout[enumeratorOut(currentCell, dof, var)] += coeff * sourceValues[var];
132 if (evaluateNonconservativeProduct) {
144 for (
int dd=0; dd<Dimensions; dd++) {
145 coeff *= massMatrixDiagonal[ dof(dd) ] * h(dd);
148 for(
int dim=0; dim<Dimensions; dim++) {
150 std::fill(gradientValues, &gradientValues[strideQ-1], 0.0);
152 const int invDx = 1 / h(dim);
155 for(
int node=0; node<nodesPerAxis; node++){
157 contributingNodeIndex[dim] = node;
158 const double coeffDerivative = invDx * derivativeOperator[node + nodesPerAxis*dof(dim)];
160 for(
int var=0; var<strideQ; var++){
161 gradientValues[var] += coeffDerivative * cellQin[enumeratorIn(currentCell, contributingNodeIndex, var)];
162 assertion(gradientValues[var]==gradientValues[var]);
166 nonconservativeProduct(
167 &cellQin[enumeratorIn(currentCell, dof, 0)],
175 for(
int var=0; var<unknowns; var++){
176 assertion(ncpValues[var]==ncpValues[var]);
177 cellQout[enumeratorOut(currentCell, dof, var)] -= coeff * ncpValues[var];
185 delete [] fluxValues;
187 if (evaluateSource) {
188 delete [] sourceValues;
190 if (evaluateNonconservativeProduct) {
191 delete [] gradientValues;
201 const int auxiliaryVariables,
202 MaxEigenvalue maxEigenvalue,
203 const double* __restrict__ QuadratureNodes1d
213 for (
int cellIndex=0; cellIndex<cellData.
numberOfCells; cellIndex++) {
214 double eigenvalue = std::numeric_limits<double>::min();
217 for(
int direction=0; direction<Dimensions; direction++) {
218 double nodeEigenvalue = maxEigenvalue(
219 cellData.
QOut[cellIndex] + enumerator(0,k,0),
227 cellData.
t[cellIndex],
228 cellData.
dt[cellIndex],
234 eigenvalue = std::max(eigenvalue, nodeEigenvalue);
262 const int auxiliaryVariables,
263 const double* __restrict__ massMatrixDiagonal1d
273 for (
int currentCell=0; currentCell<cellData.
numberOfCells; currentCell++) {
274 dfor( dof, order+1 ) {
275 double* __restrict__ cellQInOut = cellData.
QOut[currentCell];
278 double diagonal = massMatrixDiagonal1d[ dof(0) ] * massMatrixDiagonal1d[ dof(1) ] * (cellData.
cellSize[currentCell])(0) * (cellData.
cellSize[currentCell])(1);
280 double diagonal = massMatrixDiagonal1d[ dof(0) ]
281 * massMatrixDiagonal1d[ dof(1) ]
282 * massMatrixDiagonal1d[ dof(2) ]
283 * (cellData.
cellSize[currentCell])(0)
284 * (cellData.
cellSize[currentCell])(1)
285 * (cellData.
cellSize[currentCell])(2);
290 for(
int var=0; var<unknowns; var++) {
291 cellQInOut[ enumerator(currentCell,dof,var) ] = 1.0 / diagonal * cellQInOut[enumerator(currentCell,dof,var)];
void copySolution(const double *__restrict__ cellQin, const int order, const int unknowns, const int auxiliaryVariables, const tarch::la::Vector< Dimensions, int > &node, double *__restrict__ cellQout)
Copy solution over for one node.
tarch::la::Vector< Dimensions, double > getQuadraturePoint(const tarch::la::Vector< Dimensions, double > &cellCentre, const tarch::la::Vector< Dimensions, double > &cellSize, const tarch::la::Vector< Dimensions, int > &index, int polynomialOrder, const T *__restrict__ quadraturePoints)
Construct location of a quadrature point.
void reduceMaxEigenvalue_patchwise_functors(::exahype2::CellData< double, double > &cellData, const int order, const int unknowns, const int auxiliaryVariables, MaxEigenvalue maxEigenvalue, const double *__restrict__ QuadratureNodes1d)
Compute the maximum eigenvalues over a sequence of cells and store the result in the respective CellD...
void cellIntegral_patchwise_in_situ_GaussLegendre_functors(::exahype2::CellData< double, double > &cellData, const int order, const int unknowns, const int auxiliaryVariables, Flux flux, NonConservativeProduct nonconservativeProduct, Source source, PointSources pointSources, const double *__restrict__ QuadratureNodes1d, const double *__restrict__ MassMatrixDiagonal1d, const double *__restrict__ StiffnessMatrix1d, const double *__restrict__ DerivativeOperator1d, bool evaluateFlux, bool evaluateNonconservativeProduct, bool evaluateSource, bool evaluatePointSources)