guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Reformat

Browse Source
This commit is contained in:
Guilhem Lavaux 2022-11-16 16:48:40 +01:00
parent 4633f6edc9
commit fe06434619
1 changed files with 147 additions and 163 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

310
src/sphSmooth.tcc
View File

@ -3,231 +3,215 @@
namespace CosmoTool { namespace CosmoTool {
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph) SPHSmooth<ValType, Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph) {
{ this->Nsph = Nsph;
this->Nsph = Nsph; this->tree = tree;
this->tree = tree; internal.currentNgb = 0;
internal.currentNgb = 0;
this->maxNgb = Nsph; this->maxNgb = Nsph;
internal.ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[maxNgb]); internal.ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[maxNgb]);
internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]); internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
} }
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
SPHSmooth<ValType,Ndims>::~SPHSmooth() SPHSmooth<ValType, Ndims>::~SPHSmooth() {}
{
}
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
template<typename FuncT> template <typename FuncT>
ComputePrecision SPHSmooth<ValType,Ndims>::computeWValue(const typename SPHTree::coords& c, ComputePrecision SPHSmooth<ValType, Ndims>::computeWValue(
SPHCell& cell, const typename SPHTree::coords &c, SPHCell &cell, CoordType d, FuncT fun,
CoordType d, SPHState *state) {
FuncT fun, SPHState *state) CoordType weight;
{
CoordType weight;
d /= state->smoothRadius; d /= state->smoothRadius;
weight = getKernel(d); weight = getKernel(d);
if (cell.val.weight != 0) if (cell.val.weight != 0)
return weight * fun(cell.val.pValue) / cell.val.weight; return weight * fun(cell.val.pValue) / cell.val.weight;
else else
return 0; return 0;
} }
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
void void SPHSmooth<ValType, Ndims>::fetchNeighbours(
SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNngb) const typename SPHTree::coords &c, uint32_t newNngb) {
{ ComputePrecision d2, max_dist = 0;
ComputePrecision d2, max_dist = 0; uint32_t requested = newNngb;
uint32_t requested = newNngb;
if (requested > maxNgb) if (requested > maxNgb) {
{
maxNgb = requested; maxNgb = requested;
internal.ngb = boost::shared_ptr<P_SPHCell[]>(new P_SPHCell[maxNgb]); internal.ngb = boost::shared_ptr<P_SPHCell[]>(new P_SPHCell[maxNgb]);
internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]); internal.distances =
boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
} }
memcpy(internal.currentCenter, c, sizeof(c)); memcpy(internal.currentCenter, c, sizeof(c));
tree->getNearestNeighbours(c, requested, (SPHCell **)internal.ngb.get(), (CoordType*)internal.distances.get()); tree->getNearestNeighbours(
c, requested, (SPHCell **)internal.ngb.get(),
(CoordType *)internal.distances.get());
internal.currentNgb = 0; internal.currentNgb = 0;
for (uint32_t i = 0; i < requested && (internal.ngb)[i] != 0; i++,internal.currentNgb++) for (uint32_t i = 0; i < requested && (internal.ngb)[i] != 0;
{ i++, internal.currentNgb++) {
internal.distances[i] = sqrt(internal.distances[i]); internal.distances[i] = sqrt(internal.distances[i]);
d2 = internal.distances[i]; d2 = internal.distances[i];
if (d2 > max_dist) if (d2 > max_dist)
max_dist = d2; max_dist = d2;
} }
internal.smoothRadius = max_dist / 2; internal.smoothRadius = max_dist / 2;
} }
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
void SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, SPHState *state) void SPHSmooth<ValType, Ndims>::fetchNeighbours(
{ const typename SPHTree::coords &c, SPHState *state) {
ComputePrecision d2, max_dist = 0; ComputePrecision d2, max_dist = 0;
uint32_t requested = Nsph; uint32_t requested = Nsph;
if (state != 0) { if (state != 0) {
state->distances = boost::shared_ptr<CoordType[]>(new CoordType[Nsph]); state->distances = boost::shared_ptr<CoordType[]>(new CoordType[Nsph]);
state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[Nsph]); state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[Nsph]);
} else } else
state = &internal; state = &internal;
memcpy(state->currentCenter, c, sizeof(c)); memcpy(state->currentCenter, c, sizeof(c));
tree->getNearestNeighbours(c, requested, state->ngb.get(), state->distances.get()); tree->getNearestNeighbours(
c, requested, state->ngb.get(), state->distances.get());
state->currentNgb = 0; state->currentNgb = 0;
for (uint32_t i = 0; i < requested && state->ngb[i] != 0; i++,state->currentNgb++) for (uint32_t i = 0; i < requested && state->ngb[i] != 0;
{ i++, state->currentNgb++) {
d2 = state->distances[i] = sqrt(state->distances[i]); d2 = state->distances[i] = sqrt(state->distances[i]);
if (d2 > max_dist) if (d2 > max_dist)
max_dist = d2; max_dist = d2;
} }
state->smoothRadius = max_dist / 2; state->smoothRadius = max_dist / 2;
} }
template <typename ValType, int Ndims>
void SPHSmooth<ValType, Ndims>::fetchNeighboursOnVolume(
const typename SPHTree::coords &c, ComputePrecision radius) {
uint32_t numPart;
ComputePrecision d2, max_dist = 0;
template<typename ValType, int Ndims> memcpy(internal.currentCenter, c, sizeof(c));
void
SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords& c,
ComputePrecision radius)
{
uint32_t numPart;
ComputePrecision d2, max_dist = 0;
memcpy(internal.currentCenter, c, sizeof(c)); internal.currentNgb = tree->getIntersection(
c, radius, internal.ngb, internal.distances, maxNgb);
internal.currentNgb = tree->getIntersection(c, radius, internal.ngb, internal.distances, for (uint32_t i = 0; i < internal.currentNgb; i++) {
maxNgb);
for (uint32_t i = 0; i < internal.currentNgb; i++)
{
d2 = internal.distances[i] = sqrt(internal.distances[i]); d2 = internal.distances[i] = sqrt(internal.distances[i]);
if (d2 > max_dist) if (d2 > max_dist)
max_dist = d2; max_dist = d2;
} }
internal.smoothRadius = max_dist / 2; internal.smoothRadius = max_dist / 2;
} }
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
template<typename FuncT> template <typename FuncT>
ComputePrecision ComputePrecision SPHSmooth<ValType, Ndims>::computeSmoothedValue(
SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c, const typename SPHTree::coords &c, FuncT fun, SPHState *state) {
FuncT fun, SPHState *state) if (state == 0)
{ state = &internal;
if (state == 0)
state = &internal;
ComputePrecision outputValue = 0; ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0; ComputePrecision max_dist = 0;
ComputePrecision r3 = cube(state->smoothRadius); ComputePrecision r3 = cube(state->smoothRadius);
for (uint32_t i = 0; i < state->currentNgb; i++) for (uint32_t i = 0; i < state->currentNgb; i++) {
{ outputValue +=
outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun, state); computeWValue(c, *state->ngb[i], state->distances[i], fun, state);
} }
return outputValue / r3; return outputValue / r3;
} }
template<typename ValType> template <typename ValType>
ComputePrecision interpolateOne(const ValType& t) ComputePrecision interpolateOne(const ValType &t) {
{ return 1.0;
return 1.0; }
}
// WARNING ! Cell's weight must be 1 !!! // WARNING ! Cell's weight must be 1 !!!
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
template<typename FuncT> template <typename FuncT>
ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typename SPHTree::coords& c, ComputePrecision SPHSmooth<ValType, Ndims>::computeInterpolatedValue(
FuncT fun, SPHState *state) const typename SPHTree::coords &c, FuncT fun, SPHState *state) {
{ if (state == 0)
if (state == 0) state = &internal;
state = &internal;
ComputePrecision outputValue = 0; ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0; ComputePrecision max_dist = 0;
ComputePrecision weight = 0; ComputePrecision weight = 0;
for (uint32_t i = 0; i < state->currentNgb; i++) for (uint32_t i = 0; i < state->currentNgb; i++) {
{
outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun); outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun);
weight += computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne); weight +=
computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
} }
return (outputValue == 0) ? 0 : (outputValue / weight); return (outputValue == 0) ? 0 : (outputValue / weight);
} }
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
template<typename FuncT> template <typename FuncT>
void SPHSmooth<ValType,Ndims>::runForEachNeighbour(FuncT fun, SPHState *state) void
{ SPHSmooth<ValType, Ndims>::runForEachNeighbour(FuncT fun, SPHState *state) {
if (state == 0) if (state == 0)
state = &internal; state = &internal;
for (uint32_t i = 0; i < state->currentNgb; i++) for (uint32_t i = 0; i < state->currentNgb; i++) {
{
fun(state->ngb[i]); fun(state->ngb[i]);
} }
} }
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c, SPHState *state) void SPHSmooth<ValType, Ndims>::addGridSite(
{ const typename SPHTree::coords &c, SPHState *state) {
ComputePrecision outputValue = 0; ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0; ComputePrecision max_dist = 0;
ComputePrecision r3 = cube(state->smoothRadius); ComputePrecision r3 = cube(state->smoothRadius);
for (uint32_t i = 0; i < state->currentNgb; i++) for (uint32_t i = 0; i < state->currentNgb; i++) {
{
ComputePrecision d = state->distances[i]; ComputePrecision d = state->distances[i];
SPHCell& cell = *(state->ngb[i]); SPHCell &cell = *(state->ngb[i]);
double kernel_value = getKernel(d/state->smoothRadius) / r3; double kernel_value = getKernel(d / state->smoothRadius) / r3;
#pragma omp atomic #pragma omp atomic
cell.val.weight += kernel_value; cell.val.weight += kernel_value;
} }
} }
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c) void
{ SPHSmooth<ValType, Ndims>::addGridSite(const typename SPHTree::coords &c) {
addGridSite(c, &internal); addGridSite(c, &internal);
} }
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
ComputePrecision ComputePrecision
SPHSmooth<ValType,Ndims>::getKernel(ComputePrecision x) const SPHSmooth<ValType, Ndims>::getKernel(ComputePrecision x) const {
{ // WARNING !!! This is an unnormalized version of the kernel.
// WARNING !!! This is an unnormalized version of the kernel. if (x < 1)
if (x < 1) return 1 - 1.5 * x * x + 0.75 * x * x * x;
return 1 - 1.5 * x * x + 0.75 * x * x * x; else if (x < 2) {
else if (x < 2)
{
ComputePrecision d = 2 - x; ComputePrecision d = 2 - x;
return 0.25 * d * d * d; return 0.25 * d * d * d;
} } else
else return 0;
return 0; }
}
template<typename ValType, int Ndims> template <typename ValType, int Ndims>
bool SPHSmooth<ValType,Ndims>::hasNeighbours() const bool SPHSmooth<ValType, Ndims>::hasNeighbours() const {
{ return (internal.currentNgb != 0);
return (internal.currentNgb != 0); }
}
template<class ValType1, class ValType2, int Ndims> template <class ValType1, class ValType2, int Ndims>
bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2) bool operator<(
{ const SPHSmooth<ValType1, Ndims> &s1,
return (s1.getSmoothingLen() < s2.getSmoothingLen()); const SPHSmooth<ValType2, Ndims> &s2) {
} return (s1.getSmoothingLen() < s2.getSmoothingLen());
}
}; }; // namespace CosmoTool