MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
Functions
maia::dg::interpolation Namespace Reference

Holds helper functions for the interpolation. More...

Functions

void calcLegendrePolyAndDeriv (MInt Nmax, MFloat x, MFloat *polynomial, MFloat *derivative)
 Evaluates the Legendre polynomial and its derivative of degree Nmax at point x. More...
 
void calcLegendreGaussNodesAndWeights (MInt Nmax, MFloat *nodes, MFloat *wInt)
 Calculate the Gauss integration nodes and weight for the Legendre polynomials on the interval [-1,1]. More...
 
void calcQandL (MInt Nmax, MFloat x, MFloat &q, MFloat &qDeriv, MFloat &poly)
 Auxiliary function (only used by calcLegendreGaussLobattoNodesAndWeights()) More...
 
void calcLegendreGaussLobattoNodesAndWeights (MInt Nmax, MFloat *nodes, MFloat *wInt)
 Calculate the Gauss-Lobatto integration nodes and weight for the Legendre polynomials on the interval [-1,1]. More...
 
void calcBarycentricWeights (MInt Nmax, const MFloat *nodes, MFloat *weights)
 Calculates the barycentric weights for Lagrange interpolation at thei specified nodes. More...
 
MFloat getLagrangeInterpolation (MFloat x, MInt Nmax, const MFloat *nodes, const MFloat *values, const MFloat *wBary)
 Calculates the interpolated value at point x given a set of nodes, values, and weights. More...
 
void calcLagrangeInterpolatingPolynomials (const MFloat x, const MInt polyDeg, const MFloat *nodes, const MFloat *wBary, MFloat *polynomials)
 Calculates the values of the Lagrangian polynomials l_j for a given point x in [-1,1]. More...
 
void calcLhat (const MFloat x, const MInt polyDeg, const MFloat *nodes, const MFloat *wInt, const MFloat *wBary, MFloat *Lhat)
 Calculates the Lagrange polynomials evaluated at point x in [-1,1] and pre-divides them by the integration weights. More...
 
MBool ascendingAbsVal (MFloat i, MFloat j)
 
void calcPolynomialDerivativeMatrix (MInt Nmax, const MFloat *nodes, const MFloat *wBary, MFloat *derivMatrix)
 Calculates the first derivative approximation matrix. More...
 
void calcDhat (MInt Nmax, const MFloat *nodes, const MFloat *wInt, const MFloat *wBary, MFloat *dhatMatrix)
 Calculates the polynomial derivative matrix and normalizes it using the integration weights. More...
 
template<MInt nDim, MInt noVariables>
void prolongToFaceGauss (const MFloat *const source, const MInt faceId, const MInt noNodes1D, const MFloat *const LFaceN, const MFloat *const LFaceP, MFloat *const destination)
 Extrapolates ("prolongs") the Gauss node values of an element to a given face. More...
 
template<MInt nDim, MInt noVariables, class T , class U >
void prolongToFaceGaussLobatto (const T source, const MInt faceId, const MInt noNodes1D, U destination)
 Extrapolates ("prolongs") the Gauss-Lobatto node values of an element to a given face. More...
 
template<class T , class U , class V , class W >
void calcPolynomialInterpolationMatrix (MInt noNodesIn, const T nodesIn, MInt noNodesOut, const U nodesOut, const V wBary, W vandermonde)
 Calculate the polynomial interpolation matrix (Vandermonde) to interpolate from one set of nodes to another. More...
 
template<MInt nDim, class T , class U , class V >
void interpolateNodes (const T source, U vandermonde, MInt noNodesIn, MInt noNodesOut, MInt noVariables, V destination)
 
void calcSBPNodes (const MInt noNodes, MFloat *nodes)
 Generates an equidistant node distribution. More...
 
void calcSBPWeights (const MInt noNodes, MFloatVector sbpP, MFloat *wInt)
 Calulates the diagonal weight matrix (H) entries. More...
 
MBool checkSBPProp (const MFloatMatrix &Q)
 Helper function for checking the SBP property of a given operator. More...
 
void readCSV (std::string path, std::vector< std::vector< MFloat > > &data)
 Reads .csv-file on rank 0 and broadcasts it to all ranks. More...
 
void readDDRP (MFloat *nodesVector, MFloat *wIntVector, MFloat *dhatMatrix)
 Reads DDRP operator and nodal distribution from file. More...
 
void calcDhatSBP (const MInt noNodes, MFloatVector sbpA, MFloatVector sbpQ, const MFloat *wInt, MFloat *dhatMatrix)
 Calculates SBP operator from coefficients. More...
 
void readSbpOperator (MInt noNodes1D, MString opName, MFloatVector &sbpA, MFloatVector &sbpP, MFloatVector &sbpQ)
 Reads all coefficient files to contruct SBP operator. More...
 
MBool getSBPOperator (const MInt noNodes1D, const MString opName, MFloatVector &sbpA, MFloatVector &sbpP, MFloatVector &sbpQ)
 Gets SBP Coefficients from corresponding header if existent. More...
 
void findNodeIndicesAtPoint (const MFloat point, const MFloat *nodes, const MInt noNodes, MInt &idx1, MInt &idx2)
 Finds the neighboring indices to point in one dimension. More...
 
void calcLinearInterpolationBase (const MFloat x, const MInt noNodes, const MFloat *nodes, MFloat *const polynomials)
 Calculates linear interpolation base for point (1D). More...
 
template<class T , class U , class V >
void calcLinearInterpolationMatrix (const MInt noNodesIn, const T nodesIn, const MInt noNodesOut, const U nodesOut, V vandermonde)
 Calculates the linear interpolation matrix (Vandermonde) to interpolate from one set of nodes to another. More...
 
void calcBilinearInterpolation (MFloat *point, const MFloat *nodes, const MInt noNodes, const MInt noVars, const MFloat *u, MFloat *const out)
 Calculates bilinear interpolation at given point (2D) More...
 
void calcTrilinearInterpolation (MFloat *point, const MFloat *nodes, const MInt noNodes, const MInt noVars, const MFloat *u, MFloat *const out)
 Calculates trilinear interpolation at given point (3D) More...
 

Detailed Description

Holds helper functions for the construction of operators and interpolation in SBP mode.

Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-05

The following two sources may be cited:

Kopriva09: David A. Kopriva, Implementing Spectral Methods for Partial Differential Equations, Springer, 2009
HesthavenWarburton08: J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods, Springer, 2008

Author
Julian Vorspohlj.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-11-05

Function Documentation

◆ ascendingAbsVal()

MBool maia::dg::interpolation::ascendingAbsVal ( MFloat  i,
MFloat  j 
)
inline

Definition at line 409 of file dgcartesianinterpolation.h.

409{ return (std::fabs(i) < std::fabs(j)); }

◆ calcBarycentricWeights()

void maia::dg::interpolation::calcBarycentricWeights ( MInt  Nmax,
const MFloat *  nodes,
MFloat *  weights 
)
inline
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-05
Parameters
[in]NmaxThe polynomial degree.
[in]nodesNodes of the interpolating polynomial.
[out]weightsThe barycentric weights.

Taken from Kopriva09, p. 75, algorithm 30

Definition at line 277 of file dgcartesianinterpolation.h.

277 {
278 TRACE();
279
280 const MInt noNodes = Nmax + 1;
281 std::fill_n(&weights[0], noNodes, F1);
282
283 // Main iteration loop to calculate inverse weights
284 for(MInt j = 1; j < noNodes; j++) {
285 for(MInt k = 0; k < j; k++) {
286 weights[k] = weights[k] * (nodes[k] - nodes[j]);
287 weights[j] = weights[j] * (nodes[j] - nodes[k]);
288 }
289 }
290
291 // Invert previous results to get final weights
292 for(MInt j = 0; j < noNodes; j++) {
293 weights[j] = F1 / weights[j];
294 }
295}
MInt
int32_t MInt
Definition: maiatypes.h:62

◆ calcBilinearInterpolation()

void maia::dg::interpolation::calcBilinearInterpolation ( MFloat *  point,
const MFloat *  nodes,
const MInt  noNodes,
const MInt  noVars,
const MFloat *  u,
MFloat *const  out 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-05-14
Parameters
[in]pointPoint at wich interpolation is evaluated
[in]nodesNodal distribution
[in]noNodesNumber of nodes 1D
[in]noVarsNumber of variables per node
[in]uSolution field at all nodes
[out]outInterpolated solution at point

Definition at line 532 of file sbpcartesianinterpolation.h.

533 {
534 // Framing indices of point on grid
535 // nDim x 2
536 MInt idx[2][2] = {{0, 0}, {0, 0}};
537 findNodeIndicesAtPoint(point[0], nodes, noNodes, idx[0][0], idx[0][1]);
538 findNodeIndicesAtPoint(point[1], nodes, noNodes, idx[1][0], idx[1][1]);
539
540 MFloatTensor U(const_cast<MFloat*>(u), noNodes, noNodes, 1, noVars);
541 const MFloat dx = 2.0 / (noNodes - 1);
542 const MFloat dx2 = pow(dx, 2);
543 // Not sure if necessary
544 std::fill_n(out, noVars, 0.0);
545
546 // aij is the corresponding weight to the value Uij
547 // aij is the fraction of the diagonally opposite subarea
548
549 for(MInt i = 0; i < 2; i++) {
550 for(MInt j = 0; j < 2; j++) {
551 const MFloat a = fabs((nodes[idx[0][1 - i]] - point[0]) * (nodes[idx[1][1 - j]] - point[1])) / dx2;
552 for(MInt v = 0; v < noVars; v++) {
553 out[v] += a * U(idx[0][i], idx[1][j], 0, v);
554 }
555 }
556 }
557}
MFloat
double MFloat
Definition: maiatypes.h:52
a
Definition: contexttypes.h:19

◆ calcDhat()

void maia::dg::interpolation::calcDhat ( MInt  Nmax,
const MFloat *  nodes,
const MFloat *  wInt,
const MFloat *  wBary,
MFloat *  dhatMatrix 
)
inline
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-08
Parameters
[in]NmaxMaximum polynomial degree.
[in]nodesThe integration nodes.
[in]wIntThe integration weights.
[in]wBaryThe barycentric weights.
Returns
The normalized derivative matrix.

This method returns Dhat(j,n) = -D(n,j) * w(n) / w(j) as defined in Kopriva 09, p. 137, Eq. (4.139).

Definition at line 488 of file dgcartesianinterpolation.h.

488 {
489 TRACE();
490
491 // Get number of nodes and standard derivative matrix
492 const MInt noNodes = Nmax + 1;
493 MFloatMatrix Dhat(&dhatMatrix[0], noNodes, noNodes);
494 calcPolynomialDerivativeMatrix(Nmax, &nodes[0], &wBary[0], &Dhat[0]);
495
496 Dhat.transpose();
497
498 // Add normalization using integration weights
499 for(MInt j = 0; j < noNodes; j++) {
500 for(MInt n = 0; n < noNodes; n++) {
501 Dhat(j, n) = -Dhat(j, n) * wInt[n] / wInt[j];
502 }
503 }
504}
maia::dg::interpolation::calcPolynomialDerivativeMatrix
void calcPolynomialDerivativeMatrix(MInt Nmax, const MFloat *nodes, const MFloat *wBary, MFloat *derivMatrix)
Calculates the first derivative approximation matrix.
Definition: dgcartesianinterpolation.h:426

◆ calcDhatSBP()

void maia::dg::interpolation::calcDhatSBP ( const MInt  noNodes,
MFloatVector  sbpA,
MFloatVector  sbpQ,
const MFloat *  wInt,
MFloat *  dhatMatrix 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-11-05
Parameters
[in]noNodesNumber of nodes
[in]sbpACoefficients for the inner stencil
[in]sbpQCoefficients for the boundary solver
[in]wIntIntegration weights
[out]dHatMatrixSBP operator (already scaled by wInt)

Definition at line 246 of file sbpcartesianinterpolation.h.

247 {
248 TRACE();
249
250 if(noNodes <= 0) {
251 TERM(-1);
252 }
253
254 MInt lenA = sbpA.dim0();
255 MInt lenQ = sbpQ.dim0();
256 MInt lenR = (1 + sqrt(1 + 8 * lenQ)) / 2;
257
258 // Create Q Matrix from 'q' and 'a' Vector
259 MFloatMatrix Q = MFloatMatrix(noNodes, noNodes);
260 Q.set(F0);
261
262 // Write upper rxr solver
263 Q(0, 0) = -0.5;
264 MInt k = 0;
265 for(MInt i = 0; i < lenR; i++) {
266 for(MInt j = i; j < lenR - 1; j++) {
267 Q(i, j + 1) = sbpQ(k);
268 Q(j + 1, i) = -sbpQ(k);
269 k++;
270 }
271 }
272
273 // Write interior points
274 for(MInt i = lenR; i < noNodes - lenR; i++) {
275 k = lenA - 1;
276 for(MInt j = i - lenA; j < i; j++) {
277 Q(i, j) = -sbpA(k);
278 k--;
279 }
280 k = 0;
281 for(MInt j = i + 1; j < i + lenA + 1; j++) {
282 Q(i, j) = sbpA(k);
283 k++;
284 }
285 }
286
287 // Write 'a' points around the upper rxr solver
288 for(MInt i = 0; i < lenA; i++) {
289 for(MInt j = i; j < lenA; j++) {
290 Q(lenR - i - 1, lenR - i + j) = sbpA(j);
291 }
292 }
293
294 // Mirror upper rxr solver onto lower half
295 for(MInt i = 0; i < lenR; i++) {
296 for(MInt j = 0; j < lenR + lenA; j++) {
297 Q(noNodes - i - 1, noNodes - j - 1) = -Q(i, j);
298 }
299 }
300
301 Q.transpose();
302
303 // checkSBPProp(Q);
304
305 MFloatMatrix Dhat(&dhatMatrix[0], noNodes, noNodes);
306 // Add normalization using integration weights
307 for(MInt i = 0; i < noNodes; i++) {
308 for(MInt j = 0; j < noNodes; j++) {
309 Dhat(i, j) = -Q(i, j) / wInt[i];
310 }
311 }
312}
maia::tensor::Tensor::transpose
void transpose()
Transposes the first two dimensions (i.e. acting as if the tensor is a matrix).
Definition: tensor.h:306
maia::tensor::Tensor::dim0
size_type dim0() const
Return the size of dimension 0.
Definition: tensor.h:366
maia::tensor::Tensor::set
void set(const T &value)
Initializes tensor to constant value.
Definition: tensor.h:271
MFloatMatrix
maia::tensor::Tensor< MFloat > MFloatMatrix
Definition: tensor.h:1098

◆ calcLagrangeInterpolatingPolynomials()

void maia::dg::interpolation::calcLagrangeInterpolatingPolynomials ( const MFloat  x,
const MInt  polyDeg,
const MFloat *  nodes,
const MFloat *  wBary,
MFloat *  polynomials 
)
inline
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-05
Parameters
[in]xThe point at which the polynomials should be evaluated.
[in]noNodesNumber of nodes..
[in]nodesThe nodes of the Lagrange polynomials.
[in]wBaryThe barycentric weights of the Lagrange polynomials.
[out]polynomialsThe Nmax+1 Lagrange polynomials l_j, evaluated at x.

Taken from Kopriva09, p. 77, algorithm 34.

Definition at line 351 of file dgcartesianinterpolation.h.

352 {
353 // TRACE();
354 const MInt noNodes = polyDeg + 1;
355 std::fill_n(&polynomials[0], noNodes, F0);
356
357 // Catch situation where evaluation point is on node, which saves considerable
358 // computational time
359 for(MInt j = 0; j < noNodes; j++) {
360 if(approx(x, nodes[j], MFloatEps)) {
361 polynomials[j] = F1;
362 return;
363 }
364 }
365
366 // Otherwise calculate Lagrange interpolating polynomials according to
367 // Eq. (3.39) in Kopriva09
368 MFloat s = F0;
369 MFloat t = F0;
370 for(MInt j = 0; j < noNodes; j++) {
371 t = wBary[j] / (x - nodes[j]);
372 polynomials[j] = t;
373 s += t;
374 }
375 for(MInt j = 0; j < noNodes; j++) {
376 polynomials[j] = polynomials[j] / s;
377 }
378}
approx
MBool approx(const T &, const U &, const T)
Definition: functions.h:272

◆ calcLegendreGaussLobattoNodesAndWeights()

void maia::dg::interpolation::calcLegendreGaussLobattoNodesAndWeights ( MInt  Nmax,
MFloat *  nodes,
MFloat *  wInt 
)
inline
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-07
Parameters
[in]NmaxMaximum degree of the polynomials.
[out]nodesThe resulting integration nodes.
[out]wIntThe resulting integration weights.

Taken from Kopriva09, p. 66, algorithm 25.

Definition at line 198 of file dgcartesianinterpolation.h.

198 {
199 TRACE();
200 // Gauss-Lobatto only valid for degree >= 1
201 ASSERT(Nmax > 0, "Number of nodes must be greater than zero!");
202
203 // Reset nodes and weights
204 const MInt noNodes = Nmax + 1;
205 std::fill_n(&nodes[0], noNodes, F0);
206 std::fill_n(&wInt[0], noNodes, F0);
207
208 // Return values for polynomial and derivative
209 MFloat q, qDeriv, poly;
210
211 // Set tolerance and number of iterations. According to Kopriva09,
212 // 1.0E-15 and 10 should be more than sufficient.
213 const MFloat tol = F4 * MFloatEps;
214 const MInt noIterations = 10;
215
216 // Catch simple cases before going into the full loop
217 if(Nmax == 1) {
218 nodes[0] = -F1;
219 wInt[0] = F1;
220 nodes[1] = F1;
221 wInt[1] = wInt[0];
222 } else {
223 nodes[0] = -F1;
224 wInt[0] = F2 / (Nmax * (Nmax + F1));
225 nodes[Nmax] = F1;
226 wInt[Nmax] = wInt[0];
227
228 for(MInt j = 1; j < (Nmax + 1) / 2; j++) {
229 // Calculate starting guess for Newton method
230 nodes[j] = -cos((j + F1B4) * PI / Nmax - F3B8 / Nmax / PI * F1 / (j + F1B4));
231
232 // Use Newton method to find root of Legendre polynomial
233 // -> this is also the integration node
234 for(MInt k = 0; k < noIterations; k++) {
235 calcQandL(Nmax, nodes[j], q, qDeriv, poly);
236 const MFloat delta = -q / qDeriv;
237 nodes[j] += delta;
238
239 // Stop iterations if error is small enough
240 if(fabs(delta) <= tol * fabs(nodes[j])) {
241 break;
242 }
243 }
244
245 // Calculate weight
246 calcQandL(Nmax, nodes[j], q, qDeriv, poly);
247 wInt[j] = F2 / (Nmax * (Nmax + F1) * poly * poly);
248
249 // Set nodes and weights according to symmetry properties
250 nodes[Nmax - j] = -nodes[j];
251 wInt[Nmax - j] = wInt[j];
252 }
253 }
254
255 // If odd number of nodes (noNodes = Nmax + 1), set center node toi
256 // origin (0.0) and calculate weight
257 if(Nmax % 2 == 0) {
258 calcQandL(Nmax, F0, q, qDeriv, poly);
259 nodes[Nmax / 2] = F0;
260 wInt[Nmax / 2] = F2 / (Nmax * (Nmax + F1) * poly * poly);
261 }
262}
maia::dg::interpolation::calcQandL
void calcQandL(MInt Nmax, MFloat x, MFloat &q, MFloat &qDeriv, MFloat &poly)
Auxiliary function (only used by calcLegendreGaussLobattoNodesAndWeights())
Definition: dgcartesianinterpolation.h:151

◆ calcLegendreGaussNodesAndWeights()

void maia::dg::interpolation::calcLegendreGaussNodesAndWeights ( MInt  Nmax,
MFloat *  nodes,
MFloat *  wInt 
)
inline
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-07
Parameters
[in]NmaxMaximum degree of the polynomials.
[out]nodesThe resulting integration nodes.
[out]wIntThe resulting integration weights.

Taken from Kopriva09, p. 64, algorithm 23.

Definition at line 131 of file dgcartesianinterpolation.h.

131 {
132 maia::math::calculateLegendreGaussNodesAndWeights(Nmax, nodes, wInt);
133}
maia::math::calculateLegendreGaussNodesAndWeights
void calculateLegendreGaussNodesAndWeights(MInt Nmax, MFloat *nodes, MFloat *wInt)
Calculate the Gauss integration nodes and weight for the Legendre polynomials on the interval [-1,...
Definition: maiamath.h:806

◆ calcLegendrePolyAndDeriv()

void maia::dg::interpolation::calcLegendrePolyAndDeriv ( MInt  Nmax,
MFloat  x,
MFloat *  polynomial,
MFloat *  derivative 
)
inline
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-06
Parameters
[in]NmaxThe polynomial degree.
[in]xThe evaluation point.
[out]polynomialThe resulting value of the Legendre polynomial.
[out]derivativeThe resulting value of the derivative.

Taken from Kopriva09, p. 63, algorithm 22.

Definition at line 114 of file dgcartesianinterpolation.h.

114 {
115 maia::math::calculateLegendrePolyAndDeriv(Nmax, x, polynomial, derivative);
116}
maia::math::calculateLegendrePolyAndDeriv
void calculateLegendrePolyAndDeriv(MInt Nmax, MFloat x, MFloat *polynomial, MFloat *derivative)
Evaluates the Legendre polynomial and its derivative of degree Nmax at point x.
Definition: maiamath.h:754

◆ calcLhat()

void maia::dg::interpolation::calcLhat ( const MFloat  x,
const MInt  polyDeg,
const MFloat *  nodes,
const MFloat *  wInt,
const MFloat *  wBary,
MFloat *  Lhat 
)
inline
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-08
Parameters
[in]xEvaluation point.
[in]polyDegMaximum polynomial degree.
[in]nodesThe integration nodes.
[in]wIntThe integration weights.
[in]wBaryThe barycentric weights.
[out]LhatThe Lagrange polynomials divided by the integration weights, i.e. Lhat_j(x) = L_j(x) / w(j).

Definition at line 395 of file dgcartesianinterpolation.h.

396 {
397 TRACE();
398 const MInt noNodes = polyDeg + 1;
399
400 // Calculate the Lagrange polynomials and evaluate them at x
401 calcLagrangeInterpolatingPolynomials(x, polyDeg, &nodes[0], &wBary[0], Lhat);
402 // Normalize by integration weights
403 for(MInt j = 0; j < noNodes; j++) {
404 Lhat[j] = Lhat[j] / wInt[j];
405 }
406}
maia::dg::interpolation::calcLagrangeInterpolatingPolynomials
void calcLagrangeInterpolatingPolynomials(const MFloat x, const MInt polyDeg, const MFloat *nodes, const MFloat *wBary, MFloat *polynomials)
Calculates the values of the Lagrangian polynomials l_j for a given point x in [-1,...
Definition: dgcartesianinterpolation.h:351

◆ calcLinearInterpolationBase()

void maia::dg::interpolation::calcLinearInterpolationBase ( const MFloat  x,
const MInt  noNodes,
const MFloat *  nodes,
MFloat *const  polynomials 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-11-05
Note
f(x) = nodalValues * polynomials
Parameters
[in]xPoint at which the solution is interpolated
[in]noNodesNumber of nodes
[in]nodesNodal distribution
[out]polynomialsCoefficients for the linear interpolation

Definition at line 467 of file sbpcartesianinterpolation.h.

468 {
469 std::fill_n(&polynomials[0], noNodes, F0);
470
471 // Framing indices of point on grid
472 MInt idx[2] = {0, 0};
473 findNodeIndicesAtPoint(x, nodes, noNodes, idx[0], idx[1]);
474
475 const MFloat dx = 2.0 / (noNodes - 1);
476 polynomials[idx[0]] = fabs(nodes[idx[1]] - x) / dx;
477 polynomials[idx[1]] = fabs(nodes[idx[0]] - x) / dx;
478}

◆ calcLinearInterpolationMatrix()

template<class T , class U , class V >
void maia::dg::interpolation::calcLinearInterpolationMatrix ( const MInt  noNodesIn,
const T  nodesIn,
const MInt  noNodesOut,
const U  nodesOut,
V  vandermonde 
)
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-05-14
Template Parameters
T,U,VAny container type that overloads operator[] for element access (including pointers).
Parameters
[in]noNodesInNumber of incoming nodes
[in]nodesInIncoming nodal distribution
[in]noNodesOutNumber of outgoing nodes
[in]nodesOutOutgoing nodal distribution
[out]vandermondeInterpolation matrix

Definition at line 496 of file sbpcartesianinterpolation.h.

497 {
498 TRACE();
499
500 MFloatMatrix vdm(&vandermonde[0], noNodesOut, noNodesIn);
501
502 for(MInt k = 0; k < noNodesOut; k++) {
503 MBool rowHasMatch = false;
504 for(MInt j = 0; j < noNodesIn; j++) {
505 vdm(k, j) = F0;
506 if(approx(nodesOut[k], nodesIn[j], MFloatEps)) {
507 rowHasMatch = true;
508 vdm(k, j) = F1;
509 }
510 }
511
512 if(rowHasMatch == false) {
513 calcLinearInterpolationBase(nodesOut[k], noNodesIn, nodesIn, &vdm(k, 0));
514 }
515 }
516}
MBool
bool MBool
Definition: maiatypes.h:58

◆ calcPolynomialDerivativeMatrix()

void maia::dg::interpolation::calcPolynomialDerivativeMatrix ( MInt  Nmax,
const MFloat *  nodes,
const MFloat *  wBary,
MFloat *  derivMatrix 
)
inline
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-05
Parameters
[in]NmaxMaximum polynomial degree.
[in]nodesThe nodes of the Lagrange polynomial.
[in]wBaryThe barycentric weights of the Lagrange polynomial.
[out]derivMatrixThe derivative matrix.

Taken from Kopriva09, p. 82, algorithm 37.

Definition at line 426 of file dgcartesianinterpolation.h.

426 {
427 TRACE();
428
429#ifdef CALC_POLY_DERIV_SORTING
430 const MInt noNodes = Nmax + 1;
431 MFloatMatrix D(&derivMatrix[0], noNodes, noNodes);
432 std::vector<MFloat> column(noNodes);
433 D.set(F0);
434
435 // Iterate over all indices and calculate the D_ij according to Eq. (3.48).
436 // For i == j, use the negative sum trick.
437 for(MInt i = 0; i < noNodes; i++) {
438 for(MInt j = 0; j < noNodes; j++) {
439 if(j != i) {
440 D(i, j) = wBary[j] / wBary[i] * F1 / (nodes[i] - nodes[j]);
441 column[j] = D(i, j);
442 } else {
443 column[j] = F0;
444 }
445 }
446 std::sort(column.begin(), column.end(), ascendingAbsVal);
447 for(MInt j = 0; j < noNodes; j++) {
448 D(i, i) -= column[j];
449 }
450 }
451
452
453#else
454 const MInt noNodes = Nmax + 1;
455 MFloatMatrix D(&derivMatrix[0], noNodes, noNodes);
456 D.set(F0);
457
458 // Iterate over all indices and calculate the D_ij according to Eq. (3.48).
459 // For i == j, use the negative sum trick.
460 for(MInt i = 0; i < noNodes; i++) {
461 for(MInt j = 0; j < noNodes; j++) {
462 if(j != i) {
463 D(i, j) = wBary[j] / wBary[i] * F1 / (nodes[i] - nodes[j]);
464 D(i, i) -= D(i, j);
465 }
466 }
467 }
468#endif
469}

◆ calcPolynomialInterpolationMatrix()

template<class T , class U , class V , class W >
void maia::dg::interpolation::calcPolynomialInterpolationMatrix ( MInt  noNodesIn,
const T  nodesIn,
MInt  noNodesOut,
const U  nodesOut,
const V  wBary,
W  vandermonde 
)
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2013-05-14
Template Parameters
T,U,V,WAny container type that overloads operator[] for element access (including pointers).
Parameters
[in]noNodesIn
[in]nodesIn
[in]noNodesOut
[in]nodesOut
[in]wBary
[out]vandermonde

Definition at line 729 of file dgcartesianinterpolation.h.

730 {
731 TRACE();
732
733 MFloatMatrix vdm(&vandermonde[0], noNodesOut, noNodesIn);
734
735 for(MInt k = 0; k < noNodesOut; k++) {
736 MBool rowHasMatch = false;
737 for(MInt j = 0; j < noNodesIn; j++) {
738 vdm(k, j) = F0;
739 if(approx(nodesOut[k], nodesIn[j], MFloatEps)) {
740 rowHasMatch = true;
741 vdm(k, j) = F1;
742 }
743 }
744
745 if(rowHasMatch == false) {
746 MFloat s = 0;
747 for(MInt j = 0; j < noNodesIn; j++) {
748 MFloat t = wBary[j] / (nodesOut[k] - nodesIn[j]);
749 vdm(k, j) = t;
750 s += t;
751 }
752 for(MInt j = 0; j < noNodesIn; j++) {
753 vdm(k, j) = vdm(k, j) / s;
754 }
755 }
756 }
757}

◆ calcQandL()

void maia::dg::interpolation::calcQandL ( MInt  Nmax,
MFloat  x,
MFloat &  q,
MFloat &  qDeriv,
MFloat &  poly 
)
inline
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-07
Parameters
[in]NmaxThe polynomial degree.
[in]xThe evaluation point.
[out]qThe resulting q value.
[out]qDerivThe resulting q' value.
[out]polyThe resulting value of the Legendre polynomial.

Taken from Kopriva09, p. 65, algorithm 24.

Definition at line 151 of file dgcartesianinterpolation.h.

151 {
152 TRACE();
153 // This may only be called for Nmax >= 2 (see Kopriva09)
154 if(Nmax < 2) {
155 mTerm(1, AT_, "Gauss-Lobatto needs at least two nodes!");
156 }
157
158 MFloat polyLast1, polyLast2, derivLast1, derivLast2, deriv;
159
160 polyLast2 = F1;
161 polyLast1 = x;
162 derivLast2 = F0;
163 derivLast1 = F1;
164
165 for(MInt k = 2; k <= Nmax; k++) {
166 poly = (F2 * k - F1) / k * x * polyLast1 - (k - F1) / k * polyLast2;
167 deriv = derivLast2 + (F2 * k - F1) * polyLast1;
168 polyLast2 = polyLast1;
169 polyLast1 = poly;
170 derivLast2 = derivLast1;
171 derivLast1 = deriv;
172 }
173
174 MFloat polyNp1, derivNp1;
175 MInt k = Nmax + 1;
176
177 polyNp1 = (F2 * k - F1) / k * x * poly - (k - F1) / k * polyLast2;
178 derivNp1 = derivLast2 + (F2 * k - F1) * polyLast1;
179
180 q = polyNp1 - polyLast2;
181 qDeriv = derivNp1 - derivLast2;
182}
mTerm
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29

◆ calcSBPNodes()

void maia::dg::interpolation::calcSBPNodes ( const MInt  noNodes,
MFloat *  nodes 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-11-05
Parameters
[in]noNodesNumber of nodes
[out]nodesNode distribution

Definition at line 42 of file sbpcartesianinterpolation.h.

42 {
43 TRACE();
44
45 ASSERT(noNodes > 0, "Number of nodes must be greater than zero!");
46
47 // Create equidistant nodes on [-1,1]
48 const MFloat dx = F2 / (noNodes - 1);
49 std::fill_n(&nodes[0], noNodes, F0);
50 for(MInt i = 0; i < noNodes; i++) {
51 nodes[i] = -F1 + i * dx;
52 }
53}

◆ calcSBPWeights()

void maia::dg::interpolation::calcSBPWeights ( const MInt  noNodes,
MFloatVector  sbpP,
MFloat *  wInt 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-11-05
Parameters
[in]noNodesNumber of nodes
[in]sbpPCoefficients for the boundary nodes
[out]wIntIntegration matrix (1D diagonal)

Definition at line 65 of file sbpcartesianinterpolation.h.

65 {
66 TRACE();
67
68 ASSERT(noNodes > 0, "Number of nodes must be greater than zero!");
69
70 const MInt lenP = sbpP.dim0();
71 const MFloat dx = F2 / (noNodes - 1);
72
73 // Create (diagonal) H Matrix as Vector from 'p' Vector
74 std::fill_n(&wInt[0], noNodes, dx);
75 for(MInt i = 0; i < lenP; i++) {
76 wInt[i] = sbpP(i) * dx;
77 wInt[noNodes - i - 1] = sbpP(i) * dx;
78 }
79}

◆ calcTrilinearInterpolation()

void maia::dg::interpolation::calcTrilinearInterpolation ( MFloat *  point,
const MFloat *  nodes,
const MInt  noNodes,
const MInt  noVars,
const MFloat *  u,
MFloat *const  out 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-05-14
Parameters
[in]pointPoint at wich interpolation is evaluated
[in]nodesNodal distribution
[in]noNodesNumber of nodes 1D
[in]noVarsNumber of variables per node
[in]uSolution field at all nodes
[out]outInterpolated solution at point

Definition at line 572 of file sbpcartesianinterpolation.h.

573 {
574 // Framing indices of point on grid
575 // nDim x 2
576 MInt idx[3][2] = {{0, 0}, {0, 0}, {0, 0}};
577 findNodeIndicesAtPoint(point[0], nodes, noNodes, idx[0][0], idx[0][1]);
578 findNodeIndicesAtPoint(point[1], nodes, noNodes, idx[1][0], idx[1][1]);
579 findNodeIndicesAtPoint(point[2], nodes, noNodes, idx[2][0], idx[2][1]);
580
581 MFloatTensor U(const_cast<MFloat*>(u), noNodes, noNodes, noNodes, noVars);
582 const MFloat dx = 2.0 / (noNodes - 1);
583 const MFloat dx3 = pow(dx, 3);
584 // Not sure if necessary
585 std::fill_n(out, noVars, 0.0);
586
587 // aijk is the corresponding weight to the value Uijk
588 // aijk is the fraction of the diagonally opposite subvolume
589
590 for(MInt i = 0; i < 2; i++) {
591 for(MInt j = 0; j < 2; j++) {
592 for(MInt k = 0; k < 2; k++) {
593 const MFloat a = fabs((nodes[idx[0][1 - i]] - point[0]) * (nodes[idx[1][1 - j]] - point[1])
594 * (nodes[idx[2][1 - k]] - point[2]))
595 / dx3;
596
597 for(MInt v = 0; v < noVars; v++) {
598 out[v] += a * U(idx[0][i], idx[1][j], idx[2][k], v);
599 }
600 }
601 }
602 }
603}
maia::dg::interpolation::findNodeIndicesAtPoint
void findNodeIndicesAtPoint(const MFloat point, const MFloat *nodes, const MInt noNodes, MInt &idx1, MInt &idx2)
Finds the neighboring indices to point in one dimension.
Definition: sbpcartesianinterpolation.h:438

◆ checkSBPProp()

MBool maia::dg::interpolation::checkSBPProp ( const MFloatMatrix &  Q)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-11-05
Parameters
[in]Q(hopefully) SBP operator

Definition at line 89 of file sbpcartesianinterpolation.h.

89 {
90 const MInt noNodes = Q.dim(0);
91
92 MFloatMatrix B(noNodes, noNodes);
93 B.set(F0);
94 for(MInt i = 0; i < noNodes; i++) {
95 for(MInt j = 0; j < noNodes; j++) {
96 B(i, j) = Q(i, j) + Q(j, i);
97 }
98 }
99
100 MBool sbpProp = true;
101 if(!approx(B(0, 0), -F1, MFloatEps) || !approx(B(noNodes - 1, noNodes - 1), F1, MFloatEps)) {
102 sbpProp = false;
103 }
104
105 for(MInt i = 0; i < noNodes; i++) {
106 for(MInt j = 0; j < noNodes; j++) {
107 if(!approx(B(i, j), 0.0, MFloatEps) && i != j) {
108 sbpProp = false;
109 break;
110 }
111 if(!sbpProp) break;
112 }
113 }
114 std::cout << "SBP Propery " << sbpProp << std::endl;
115 return sbpProp;
116}
maia::tensor::Tensor::dim
size_type dim(size_type d) const
Return the size of dimension d.
Definition: tensor.h:349

◆ findNodeIndicesAtPoint()

void maia::dg::interpolation::findNodeIndicesAtPoint ( const MFloat  point,
const MFloat *  nodes,
const MInt  noNodes,
MInt &  idx1,
MInt &  idx2 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-11-05
Note
Even if point is collocating with a node, two different indices are returned.
Parameters
[in]point1D Point
[in]nodesNodal distribution
[in]noNodesNumber of nodes
[out]idx1First index (more negativer)

Definition at line 438 of file sbpcartesianinterpolation.h.

439 {
440 // TODO labels:DG Change to bisection/Catch collocating points
441 for(MInt i = 0; i < noNodes; i++) {
442 if(nodes[i] > point || approx(nodes[i], point, MFloatEps)) {
443 idx1 = i - 1;
444 idx2 = i;
445 break;
446 }
447 }
448 if(idx1 == -1) {
449 idx1++;
450 idx2++;
451 }
452}

◆ getLagrangeInterpolation()

MFloat maia::dg::interpolation::getLagrangeInterpolation ( MFloat  x,
MInt  Nmax,
const MFloat *  nodes,
const MFloat *  values,
const MFloat *  wBary 
)
inline
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-11-07
Parameters
[in]xThe point which should be evaluated.
[in]NmaxThe maximum polynomial degree.
[in]nodesThe set of Nmax+1 node locations.
[in]valuesThe set of values at the nodes.
[in]wBaryThe barycentric weights for the given set of nodes.
Returns
The interpolated value.

Taken from Kopriva09, p. 75, algorithm 31.

Definition at line 315 of file dgcartesianinterpolation.h.

316 {
317 // TRACE();
318
319 const MInt noNodes = Nmax + 1;
320 MFloat numerator = F0;
321 MFloat denominator = F0;
322
323 for(MInt j = 0; j < noNodes; j++) {
324 if(approx(x, nodes[j], MFloatEps)) {
325 return values[j];
326 }
327 MFloat t = wBary[j] / (x - nodes[j]);
328 numerator += t * values[j];
329 denominator += t;
330 }
331
332 return numerator / denominator;
333}

◆ getSBPOperator()

MBool maia::dg::interpolation::getSBPOperator ( const MInt  noNodes1D,
const MString  opName,
MFloatVector &  sbpA,
MFloatVector &  sbpP,
MFloatVector &  sbpQ 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2020-01-01
Parameters
[in]noNodes1DNumber of nodes
[in]opNameName of SBP operator
[out]sbpACoefficients for the inner stencil
[out]sbpPCoefficients for the integration matrix
[out]sbpQCoefficients for the boundary solver
Returns
true if operator exists in header

Definition at line 381 of file sbpcartesianinterpolation.h.

382 {
383 SBPOperator* op = nullptr;
384 MBool exists = false;
385
386 if(opName == "go4/s306") {
387 op = &s306;
388 exists = true;
389 } else {
390 return exists = false;
391 }
392
393 MInt lenA = op->a.size();
394 MInt lenP = op->p.size();
395 MInt lenR = lenP;
396 MInt lenQ = lenP * (lenP - 1) / 2;
397
398 ASSERT(noNodes1D >= 2 * lenR,
399 "Not enough nodes for the " + opName + " operator!" << noNodes1D << " set but " << 2 * lenR << " required.");
400
401 sbpA = MFloatVector(lenA);
402 MInt i = 0;
403 for(auto& a : op->a) {
404 sbpA(i) = a;
405 i++;
406 }
407
408 sbpQ = MFloatVector(lenQ);
409 i = 0;
410 for(auto& q : op->q) {
411 sbpQ(i) = q;
412 i++;
413 }
414
415 sbpP = MFloatVector(lenP);
416 i = 0;
417 for(auto& p : op->p) {
418 sbpP(i) = p;
419 i++;
420 }
421
422 return exists;
423}
maia::fc::cell::p
constexpr std::underlying_type< FcCell >::type p(const FcCell property)
Converts property name to underlying integer value.
Definition: fccellproperties.h:32
SBPOperator::a
std::vector< MFloat > a
Definition: sbpcartesianoperators.h:19
SBPOperator::p
std::vector< MFloat > p
Definition: sbpcartesianoperators.h:20
SBPOperator::q
std::vector< MFloat > q
Definition: sbpcartesianoperators.h:21
MFloatVector
maia::tensor::Tensor< MFloat > MFloatVector
Definition: tensor.h:1095

◆ interpolateNodes()

template<MInt nDim, class T , class U , class V >
void maia::dg::interpolation::interpolateNodes ( const T  source,
U  vandermonde,
MInt  noNodesIn,
MInt  noNodesOut,
MInt  noVariables,
V  destination 
)
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2013-05-14
Template Parameters
nDimNumber of spatial dimensions of the element.
T,U,VAny container type that overloads operator[] for element access (including pointers).
Parameters
[in]source
[in]vandermonde
[in]noNodesIn
[in]noNodesOut
[in]noVariables
[out]destination

Definition at line 777 of file dgcartesianinterpolation.h.

777 {
778 // TRACE();
779 const MInt noNodesIn1D = noNodesIn;
780 const MInt noNodesIn1D2 = (nDim == 2 || nDim == 3) ? noNodesIn1D : 1;
781 const MInt noNodesIn1D3 = (nDim == 3) ? noNodesIn1D : 1;
782 const MInt noNodesOut1D = noNodesOut;
783 const MInt noNodesOut1D2 = (nDim == 2 || nDim == 3) ? noNodesOut1D : 1;
784 const MInt noNodesOut1D3 = (nDim == 3) ? noNodesOut1D : 1;
785 const MFloatTensor src(&source[0], noNodesIn, noNodesIn1D2, noNodesIn1D3, noVariables);
786 MFloatTensor dest(&destination[0], noNodesOut, noNodesOut1D2, noNodesOut1D3, noVariables);
787 const MFloatMatrix vdm(&vandermonde[0], noNodesOut, noNodesIn);
788 dest.set(F0);
789
790 // Allocate buffers (if necessary) for intermediate results
791 MFloat* buf1Ptr = nullptr;
792 MFloat* buf2Ptr = nullptr;
793 MFloatTensor buffer1, buffer2;
794 switch(nDim) {
795 case 1: {
796 buf1Ptr = &destination[0];
797 buf2Ptr = nullptr;
798 } break;
799
800 case 2: {
801 buffer1.resize(noNodesOut, noNodesIn1D2, noNodesIn1D3, noVariables);
802 buf1Ptr = &buffer1[0];
803 buf2Ptr = &destination[0];
804 } break;
805
806 case 3: {
807 buffer1.resize(noNodesOut, noNodesIn1D2, noNodesIn1D3, noVariables);
808 buffer2.resize(noNodesOut, noNodesOut, noNodesIn, noVariables);
809 buf1Ptr = &buffer1[0];
810 buf2Ptr = &buffer2[0];
811 } break;
812
813 default:
814 mTerm(1, AT_, "Bad dimension - must be 1, 2 or 3.");
815 }
816 MFloatTensor buf1(buf1Ptr, noNodesOut, noNodesIn1D2, noNodesIn1D3, noVariables);
817 buf1.set(F0);
818
819 // x-direction
820 for(MInt i = 0; i < noNodesOut1D; i++) {
821 for(MInt j = 0; j < noNodesIn1D2; j++) {
822 for(MInt k = 0; k < noNodesIn1D3; k++) {
823 for(MInt ii = 0; ii < noNodesIn1D; ii++) {
824 for(MInt l = 0; l < noVariables; l++) {
825 buf1(i, j, k, l) += vdm(i, ii) * src(ii, j, k, l);
826 }
827 }
828 }
829 }
830 }
831
832 // Return here for 1D since buf2Ptr is a nullptr
833 IF_CONSTEXPR(nDim == 1) { return; }
834
835 MFloatTensor buf2(buf2Ptr, noNodesOut, noNodesOut1D2, noNodesIn1D3, noVariables);
836 buf2.set(F0);
837
838 // Interpolate in y-direction only for 2D & 3D
839 IF_CONSTEXPR(nDim == 2 || nDim == 3) {
840 // y-direction
841 for(MInt i = 0; i < noNodesOut1D; i++) {
842 for(MInt j = 0; j < noNodesOut1D; j++) {
843 for(MInt k = 0; k < noNodesIn1D3; k++) {
844 for(MInt jj = 0; jj < noNodesIn1D; jj++) {
845 for(MInt l = 0; l < noVariables; l++) {
846 buf2(i, j, k, l) += vdm(j, jj) * buf1(i, jj, k, l);
847 }
848 }
849 }
850 }
851 }
852 }
853
854 // Interpolate in z-direction only for 3D
855 IF_CONSTEXPR(nDim == 3) {
856 // z-direction
857 for(MInt i = 0; i < noNodesOut1D; i++) {
858 for(MInt j = 0; j < noNodesOut1D; j++) {
859 for(MInt k = 0; k < noNodesOut1D; k++) {
860 for(MInt kk = 0; kk < noNodesIn1D; kk++) {
861 for(MInt l = 0; l < noVariables; l++) {
862 dest(i, j, k, l) += vdm(k, kk) * buf2(i, j, kk, l);
863 }
864 }
865 }
866 }
867 }
868 }
869}
maia::tensor::Tensor::resize
size_type resize(size_type n0, size_type n1=1, size_type n2=1, size_type n3=1, size_type n4=1)
Deletes the old data structure and creates a new one with the requested dimensions.
Definition: tensor.h:230

◆ prolongToFaceGauss()

template<MInt nDim, MInt noVariables>
void maia::dg::interpolation::prolongToFaceGauss ( const MFloat *const  source,
const MInt  faceId,
const MInt  noNodes1D,
const MFloat *const  LFaceN,
const MFloat *const  LFaceP,
MFloat *const  destination 
)
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-12-07
Template Parameters
nDimNumber of spatial dimensions of the element.
noVariablesNumber of variables per node.
T,U,VAny container type that overloads operator[] for element access (including pointers).
Parameters
[in]sourcePointer to memory where the values at the nodes are located.
[in]faceIdFace to which the values are prolonged ([-x, x, -y, y] = [0, 1, 2, 3])
[in]polyDegPolynomial degree inside the element.
[in]LFaceArray with Lagrange polynomials at x = -1.0 (index 0) and +1.0 (index 1)
[out]destinationPointer to memory where the extrapolated values should be stored.

Definition at line 528 of file dgcartesianinterpolation.h.

533 {
534 // TRACE();
535
536 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
537 const MFloatTensor src(const_cast<MFloat*>(source), noNodes1D, noNodes1D, noNodes1D3, noVariables);
538 MFloatTensor dest(destination, noNodes1D, noNodes1D3, noVariables);
539
540 // Reset destination to zero
541 std::fill_n(destination, noNodes1D * noNodes1D3 * noVariables, F0);
542
543 // Index for x-direction: i
544 // Index for y-direction: j
545 // Index for z-direction: k
546 // Index for variables: l
547 switch(faceId) {
548 case 0: // prolong to negative x-direction
549 for(MInt i = 0; i < noNodes1D; i++) {
550 for(MInt j = 0; j < noNodes1D; j++) {
551 for(MInt k = 0; k < noNodes1D3; k++) {
552 for(MInt l = 0; l < noVariables; l++) {
553 dest(j, k, l) += src(i, j, k, l) * LFaceN[i];
554 }
555 }
556 }
557 }
558 break;
559 case 1: // prolong to positive x-direction
560 for(MInt i = 0; i < noNodes1D; i++) {
561 for(MInt j = 0; j < noNodes1D; j++) {
562 for(MInt k = 0; k < noNodes1D3; k++) {
563 for(MInt l = 0; l < noVariables; l++) {
564 dest(j, k, l) += src(i, j, k, l) * LFaceP[i];
565 }
566 }
567 }
568 }
569 break;
570 case 2: // prolong to negative y-direction
571 for(MInt i = 0; i < noNodes1D; i++) {
572 for(MInt j = 0; j < noNodes1D; j++) {
573 for(MInt k = 0; k < noNodes1D3; k++) {
574 for(MInt l = 0; l < noVariables; l++) {
575 dest(i, k, l) += src(i, j, k, l) * LFaceN[j];
576 }
577 }
578 }
579 }
580 break;
581 case 3: // prolong to positive y-direction
582 for(MInt i = 0; i < noNodes1D; i++) {
583 for(MInt j = 0; j < noNodes1D; j++) {
584 for(MInt k = 0; k < noNodes1D3; k++) {
585 for(MInt l = 0; l < noVariables; l++) {
586 dest(i, k, l) += src(i, j, k, l) * LFaceP[j];
587 }
588 }
589 }
590 }
591 break;
592 case 4: // prolong to negative z-direction
593 for(MInt i = 0; i < noNodes1D; i++) {
594 for(MInt j = 0; j < noNodes1D; j++) {
595 for(MInt k = 0; k < noNodes1D; k++) {
596 for(MInt l = 0; l < noVariables; l++) {
597 dest(i, j, l) += src(i, j, k, l) * LFaceN[k];
598 }
599 }
600 }
601 }
602 break;
603 case 5: // prolong to positive z-direction
604 for(MInt i = 0; i < noNodes1D; i++) {
605 for(MInt j = 0; j < noNodes1D; j++) {
606 for(MInt k = 0; k < noNodes1D; k++) {
607 for(MInt l = 0; l < noVariables; l++) {
608 dest(i, j, l) += src(i, j, k, l) * LFaceP[k];
609 }
610 }
611 }
612 }
613 break;
614 default:
615 mTerm(1, AT_, "Bad face id.");
616 }
617}

◆ prolongToFaceGaussLobatto()

template<MInt nDim, MInt noVariables, class T , class U >
void maia::dg::interpolation::prolongToFaceGaussLobatto ( const T  source,
const MInt  faceId,
const MInt  noNodes1D,
U  destination 
)
Author
Michael Schlottke (mic) mic@a.nosp@m.ia.r.nosp@m.wth-a.nosp@m.ache.nosp@m.n.de
Date
2012-12-07
Template Parameters
nDimNumber of spatial dimensions of the element.
noVariablesNumber of variables per node.
T,U,V,WAny container type that overloads operator[] for element access (including pointers).
Parameters
[in]sourcePointer to memory where the values at the nodes are located.
[in]faceIdFace to which the values are prolonged ([-x, x, -y, y] = [0, 1, 2, 3])
[in]polyDegPolynomial degree inside the element.
[out]destinationPointer to memory where the extrapolated values should be stored.

Definition at line 640 of file dgcartesianinterpolation.h.

640 {
641 // TRACE();
642
643 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
644 const MFloatTensor src(&source[0], noNodes1D, noNodes1D, noNodes1D3, noVariables);
645 MFloatTensor dest(&destination[0], noNodes1D, noNodes1D3, noVariables);
646
647 // Index for x-direction: i
648 // Index for y-direction: j
649 // Index for z-direction: k
650 // Index for variables: l
651 switch(faceId) {
652 case 0: // prolong to negative x-direction
653 for(MInt j = 0; j < noNodes1D; j++) {
654 for(MInt k = 0; k < noNodes1D3; k++) {
655 for(MInt l = 0; l < noVariables; l++) {
656 dest(j, k, l) = src(0, j, k, l);
657 }
658 }
659 }
660 break;
661 case 1: // prolong to positive x-direction
662 for(MInt j = 0; j < noNodes1D; j++) {
663 for(MInt k = 0; k < noNodes1D3; k++) {
664 for(MInt l = 0; l < noVariables; l++) {
665 dest(j, k, l) = src(noNodes1D - 1, j, k, l);
666 }
667 }
668 }
669 break;
670 case 2: // prolong to negative y-direction
671 for(MInt i = 0; i < noNodes1D; i++) {
672 for(MInt k = 0; k < noNodes1D3; k++) {
673 for(MInt l = 0; l < noVariables; l++) {
674 dest(i, k, l) = src(i, 0, k, l);
675 }
676 }
677 }
678 break;
679 case 3: // prolong to positive y-direction
680 for(MInt i = 0; i < noNodes1D; i++) {
681 for(MInt k = 0; k < noNodes1D3; k++) {
682 for(MInt l = 0; l < noVariables; l++) {
683 dest(i, k, l) = src(i, noNodes1D - 1, k, l);
684 }
685 }
686 }
687 break;
688 case 4: // prolong to negative z-direction
689 for(MInt i = 0; i < noNodes1D; i++) {
690 for(MInt j = 0; j < noNodes1D; j++) {
691 for(MInt l = 0; l < noVariables; l++) {
692 dest(i, j, l) = src(i, j, 0, l);
693 }
694 }
695 }
696 break;
697 case 5: // prolong to positive z-direction
698 for(MInt i = 0; i < noNodes1D; i++) {
699 for(MInt j = 0; j < noNodes1D; j++) {
700 for(MInt l = 0; l < noVariables; l++) {
701 dest(i, j, l) = src(i, j, noNodes1D - 1, l);
702 }
703 }
704 }
705 break;
706 default:
707 mTerm(1, AT_, "Bad face id.");
708 }
709}

◆ readCSV()

void maia::dg::interpolation::readCSV ( std::string  path,
std::vector< std::vector< MFloat > > &  data 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2020-01-01
Parameters
[in]pathAbsolute path to .csv-file
[out]dataRead data

Definition at line 127 of file sbpcartesianinterpolation.h.

127 {
128 std::ifstream stream(path);
129
130 MInt rank;
131 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
132
133 if(!rank) {
134 // Read data on rank 0
135 MInt l = 0;
136 while(stream) {
137 l++;
138 std::string s;
139 if(!getline(stream, s)) break;
140 if(s[0] != '#') {
141 std::istringstream ss(s);
142 std::vector<MFloat> record;
143 while(ss) {
144 std::string line;
145 if(!getline(ss, line, ',')) break;
146 record.push_back(std::stof(line));
147 }
148 data.push_back(record);
149 }
150 }
151
152 if(!stream.eof()) {
153 std::cerr << "Could not read file " << path << "\n";
154 std::__throw_invalid_argument("File not found.");
155 }
156
157 MInt dim0 = data.size();
158 MInt dim1 = data[0].size();
159
160 // transfer read data to sendBuffer
161 MFloatVector sendBuffer;
162 sendBuffer.resize(dim0 * dim1);
163 MInt it = 0;
164 for(auto& line : data) {
165 for(auto& coef : line) {
166 sendBuffer[it] = coef;
167 it++;
168 }
169 }
170
171 // Send data
172 MPI_Bcast(&dim0, 1, MPI_INT, 0, MPI_COMM_WORLD, AT_, "dim0");
173 MPI_Bcast(&dim1, 1, MPI_INT, 0, MPI_COMM_WORLD, AT_, "dim1");
174 MPI_Bcast(&sendBuffer[0], dim0 * dim1, MPI_DOUBLE, 0, MPI_COMM_WORLD, AT_, "sendBuffer");
175
176 } else {
177 // Receive data
178 MInt dim0;
179 MInt dim1;
180 MPI_Bcast(&dim0, 1, MPI_INT, 0, MPI_COMM_WORLD, AT_, "dim0");
181 MPI_Bcast(&dim1, 1, MPI_INT, 0, MPI_COMM_WORLD, AT_, "dim1");
182
183 MFloatVector recvBuffer;
184 recvBuffer.resize(dim0 * dim1);
185 MPI_Bcast(&recvBuffer[0], dim0 * dim1, MPI_DOUBLE, 0, MPI_COMM_WORLD, AT_, "recvBuffer");
186
187 for(MInt i = 0; i < dim0; i++) {
188 std::vector<MFloat> line;
189 for(MInt j = 0; j < dim1; j++) {
190 line.push_back(recvBuffer(dim1 * i + j));
191 }
192 data.push_back(line);
193 }
194 }
195}
MPI_Bcast
int MPI_Bcast(void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm, const MString &name, const MString &varname)
same as MPI_Bcast
Definition: mpioverride.cpp:1114

◆ readDDRP()

void maia::dg::interpolation::readDDRP ( MFloat *  nodesVector,
MFloat *  wIntVector,
MFloat *  dhatMatrix 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-11-05
Parameters
[out]nodesVectorNodal distribution from file
[out]wIntVectorIntegration matrix (1D diagonal)
[out]dHatMatrixDDRP operator (already scaled by wInt)

Definition at line 207 of file sbpcartesianinterpolation.h.

207 {
208 std::string path = "./operatorsSBP/go4/DDRP307_N32/";
209 std::vector<std::vector<MFloat>> dData;
210 std::vector<std::vector<MFloat>> wData;
211 std::vector<std::vector<MFloat>> xData;
212
213 readCSV(path + "diff_matrix.csv", dData);
214 readCSV(path + "weights_int.csv", wData);
215 readCSV(path + "nodes.csv", xData);
216
217 const MInt nNodes = xData.size();
218
219 MFloatMatrix Dhat(&dhatMatrix[0], nNodes, nNodes);
220 MFloatVector wInt(&wIntVector[0], nNodes);
221 MFloatVector nodes(&nodesVector[0], nNodes);
222
223 // Copy the 2D Vector to the corresponding data structures
224 // Transpose and normalize differentiation matrix D
225 for(MInt i = 0; i < nNodes; i++) {
226 for(MInt j = 0; j < nNodes; j++) {
227 Dhat(i, j) = dData[j][i] * wData[j][0] / wData[i][0];
228 }
229 wInt(i) = wData[i][0];
230 nodes(i) = xData[i][0];
231 }
232}
maia::dg::interpolation::readCSV
void readCSV(std::string path, std::vector< std::vector< MFloat > > &data)
Reads .csv-file on rank 0 and broadcasts it to all ranks.
Definition: sbpcartesianinterpolation.h:127

◆ readSbpOperator()

void maia::dg::interpolation::readSbpOperator ( MInt  noNodes1D,
MString  opName,
MFloatVector &  sbpA,
MFloatVector &  sbpP,
MFloatVector &  sbpQ 
)
inline
Author
Julian Vorspohl j.vor.nosp@m.spoh.nosp@m.l@aia.nosp@m..rwt.nosp@m.h-aac.nosp@m.hen..nosp@m.de
Date
2019-11-05
Parameters
[in]noNodes1DNumber of nodes
[in]opNameName of the SBP operator
[out]sbpACoefficients for the inner stencil
[out]sbpPCoefficients for the integration matrix
[out]sbpQCoefficients for the boundary solver

Definition at line 326 of file sbpcartesianinterpolation.h.

327 {
328 MString path = "./operatorsSBP/" + opName;
329 std::vector<std::vector<MFloat>> aData;
330 std::vector<std::vector<MFloat>> pData;
331 std::vector<std::vector<MFloat>> qData;
332
333 readCSV(path + "/a.csv", aData);
334 readCSV(path + "/p.csv", pData);
335 readCSV(path + "/q.csv", qData);
336
337 MInt lenA = aData.size();
338 MInt lenP = pData.size();
339 MInt lenR = lenP;
340 MInt lenQ = lenP * (lenP - 1) / 2;
341
342 ASSERT(noNodes1D >= 2 * lenR,
343 "Not enough nodes for the " + opName + " operator! " << noNodes1D << " set but " << 2 * lenR << " required.");
344
345 sbpA = MFloatVector(lenA);
346 MInt i = 0;
347 for(auto& a : aData) {
348 sbpA(i) = a[0];
349 i++;
350 }
351
352 sbpQ = MFloatVector(lenQ);
353 i = 0;
354 for(auto& q : qData) {
355 sbpQ(i) = q[0];
356 i++;
357 }
358
359 sbpP = MFloatVector(lenP);
360 i = 0;
361 for(auto& p : pData) {
362 sbpP(i) = p[0];
363 i++;
364 }
365}
MString
std::basic_string< char > MString
Definition: maiatypes.h:55