Peano 4
Loading...
Searching...
No Matches
Tracer.cpp
Go to the documentation of this file.
1// This file is part of the ExaHyPE2 project. For conditions of distribution and
2// use, please see the copyright notice at www.peano-framework.org
3#include "Tracer.h"
4
5#include "PatchUtils.h"
6#include "peano4/utils/Loop.h"
7
10 int voxelsPerAxis,
12) {
14
15 assertion(voxelsPerAxis > 0);
16 double cartesianMeshH = marker.h()(0) / voxelsPerAxis;
17
18 for (int d = 0; d < Dimensions; d++) {
19 double positionWithinPatch = particleX(d) - marker.getOffset()(d);
20 result(d) = static_cast<int>(std::floor(positionWithinPatch / cartesianMeshH));
21 result(d) = std::max(0, result(d));
22 result(d) = std::min(voxelsPerAxis - 1, result(d));
23 }
24
25 return result;
26}
27
30 int voxelsPerAxis,
32) {
34
35 assertion3(voxelsPerAxis>=2, marker, voxelsPerAxis, particleX);
36 double cartesianMeshH = marker.h()(0) / voxelsPerAxis;
37
38 for (int d = 0; d < Dimensions; d++) {
39 double biasedPositionWithinPatch = particleX(d) - marker.getOffset()(d) - cartesianMeshH / 2.0;
40 result(d) = static_cast<int>(std::floor(biasedPositionWithinPatch / cartesianMeshH));
41 result(d) = std::max(0, result(d));
42 result(d) = std::min(voxelsPerAxis - 2, result(d));
43 }
44
45 return result;
46}
47
50 int voxelsPerAxis,
51 int unknownsPerVoxel,
52 const double* __restrict__ Q,
54 int unknown
55) {
57 marker, voxelsPerAxis, particleX
58 );
59
60 int voxelIndex = peano4::utils::dLinearised(voxel, voxelsPerAxis);
61
62 return Q[voxelIndex * unknownsPerVoxel + unknown];
63}
64
67 int voxelsPerAxis,
68 int unknownsPerVoxel,
69 const double* __restrict__ Q,
71 int unknown
72) {
73 tarch::la::Vector<Dimensions, int> biasedVoxel = mapBiasedParticleOntoVoxel(marker, voxelsPerAxis, particleX);
74
75 double h = getVolumeLength(marker.h(), voxelsPerAxis);
76
77 double result = 0.0;
78 dfor2(k)
79 tarch::la::Vector<Dimensions, int> currentVoxel = biasedVoxel + k;
80 tarch::la::Vector<Dimensions, double> currentVoxelX = getVolumeCentre(
81 marker.x(), marker.h(), voxelsPerAxis, currentVoxel
82 );
83 int currentVoxelIndex = peano4::utils::dLinearised(currentVoxel, voxelsPerAxis);
84 double weight = 1.0;
85 for (int d = 0; d < Dimensions; d++) {
86 double relativeDistance = std::abs(currentVoxelX(d) - particleX(d)) / h;
87 if (particleX(d) < currentVoxelX(d) and currentVoxel(d) == 0) {
88 weight = 1.0;
89 }
90 else if (particleX(d) > currentVoxelX(d) and currentVoxel(d) == voxelsPerAxis - 1) {
91 weight = 1.0;
92 } else if (relativeDistance < 1.0) {
93 weight *= (1.0 - relativeDistance);
94 } else {
95 weight = 0.0;
96 }
97 }
98 if (weight > 0.0) {
99 result += weight * Q[currentVoxelIndex * unknownsPerVoxel + unknown];
100 }
102 return result;
103}
#define assertion3(expr, param0, param1, param2)
#define assertion(expr)
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 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
double projectValueOntoParticle_piecewiseConstant(const peano4::datamanagement::CellMarker &marker, int voxelsPerAxis, int unknownsPerVoxel, const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &particleX, int unknown)
Project one quantity from the patch data onto the particle.
Definition Tracer.cpp:48
double projectValueOntoParticle_piecewiseLinear(const peano4::datamanagement::CellMarker &marker, int voxelsPerAxis, int unknownsPerVoxel, const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &particleX, int unknown)
Definition Tracer.cpp:65
tarch::la::Vector< Dimensions, int > mapParticleOntoVoxel(const peano4::datamanagement::CellMarker &marker, int voxelsPerAxis, const tarch::la::Vector< Dimensions, double > &particleX)
Assume that we have a particle suspended in a cell.
Definition Tracer.cpp:8
tarch::la::Vector< Dimensions, int > mapBiasedParticleOntoVoxel(const peano4::datamanagement::CellMarker &marker, int voxelsPerAxis, const tarch::la::Vector< Dimensions, double > &particleX)
Similar to mapParticleOntoVoxel() but this time, we always take the biased voxel to the left.
Definition Tracer.cpp:28
double getVolumeLength(const tarch::la::Vector< 2, double > &h, int numberOfVolumesPerAxisInPatch)
With GCC 10, it was impossible to return/copy the vector class.
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
Definition Loop.cpp:106
Simple vector class.
Definition Vector.h:134