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

Fixed SPH smoothing. Fixlet for dinterpolate

Browse Source
This commit is contained in:
Guilhem Lavaux 2009-09-29 16:26:56 -05:00
parent afd706c1dc
commit b18f06b606
3 changed files with 23 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/dinterpolate.tcc
View File

@ -98,7 +98,8 @@ namespace CosmoTool {
int signum; int signum;
gsl_linalg_LU_decomp(&M.matrix, p, &signum); gsl_linalg_LU_decomp(&M.matrix, p, &signum);
if (fabs(gsl_linalg_LU_det(&M.matrix, signum)) < 1e-10) double a = fabs(gsl_linalg_LU_det(&M.matrix, signum));
if (a < 1e-10)
throw InvalidArgumentException("Invalid tesselation. One simplex is coplanar."); throw InvalidArgumentException("Invalid tesselation. One simplex is coplanar.");
gsl_linalg_LU_invert(&M.matrix, p, &iM.matrix); gsl_linalg_LU_invert(&M.matrix, p, &iM.matrix);

14
src/sphSmooth.hpp
View File

@ -25,7 +25,7 @@ namespace CosmoTool
typedef void (*runParticleValue)(ValType& t); typedef void (*runParticleValue)(ValType& t);
public: public:
SPHSmooth(SPHTree *tree, uint32_t Nsph, uint32_t deltaN); SPHSmooth(SPHTree *tree, uint32_t Nsph);
virtual ~SPHSmooth(); virtual ~SPHSmooth();
void fetchNeighbours(const typename SPHTree::coords& c); void fetchNeighbours(const typename SPHTree::coords& c);
@ -63,6 +63,18 @@ namespace CosmoTool
SPHTree *getTree() { return tree; } SPHTree *getTree() { return tree; }
void changeNgb(uint32_t newMax) {
delete[] ngb;
delete[] distances;
ngb = new SPHCell *[newMax];
distances = new CoordType[newMax];
maxNgb = newMax;
}
uint32_t getCurrent() const { return currentNgb; }
uint32_t getNgb() const { return maxNgb; }
protected: protected:
SPHCell **ngb; SPHCell **ngb;
CoordType *distances; CoordType *distances;

18
src/sphSmooth.tcc
View File

@ -3,14 +3,13 @@
namespace CosmoTool { namespace CosmoTool {
template<typename ValType, int Ndims> template<typename ValType, int Ndims>
SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph, uint32_t deltaN) SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph)
{ {
this->Nsph = Nsph; this->Nsph = Nsph;
this->deltaNsph = deltaN;
this->tree = tree; this->tree = tree;
this->currentNgb = 0; this->currentNgb = 0;
maxNgb = Nsph; this->maxNgb = Nsph;
ngb = new SPHCell *[maxNgb]; ngb = new SPHCell *[maxNgb];
distances = new CoordType[maxNgb]; distances = new CoordType[maxNgb];
} }
@ -43,7 +42,6 @@ template<typename ValType, int Ndims>
void void
SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c) SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c)
{ {
typename SPHTree::coords radius;
ComputePrecision d2, max_dist = 0; ComputePrecision d2, max_dist = 0;
uint32_t requested = Nsph; uint32_t requested = Nsph;
@ -53,6 +51,7 @@ SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c)
currentNgb = 0; currentNgb = 0;
for (uint32_t i = 0; i < maxNgb && ngb[i] != 0; i++,currentNgb++) for (uint32_t i = 0; i < maxNgb && ngb[i] != 0; i++,currentNgb++)
{ {
distances[i] = sqrt(distances[i]);
d2 = distances[i]; d2 = distances[i];
if (d2 > max_dist) if (d2 > max_dist)
max_dist = d2; max_dist = d2;
@ -66,18 +65,17 @@ void
SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords& c, SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords& c,
ComputePrecision radius) ComputePrecision radius)
{ {
typename SPHTree::coords allRadii = {radius,radius,radius};
uint32_t numPart; uint32_t numPart;
ComputePrecision d2, max_dist = 0; ComputePrecision d2, max_dist = 0;
memcpy(currentCenter, c, sizeof(c)); memcpy(currentCenter, c, sizeof(c));
numPart = tree->getIntersection(c, ngb, distances, currentNgb = tree->getIntersection(c, radius, ngb, distances,
maxNgb); maxNgb);
currentNgb = 0; for (uint32_t i = 0; i < currentNgb; i++)
for (uint32_t i = 0; i < maxNgb && ngb[i] != 0; i++,currentNgb++)
{ {
distances[i] = sqrt(distances[i]);
d2 = distances[i]; d2 = distances[i];
if (d2 > max_dist) if (d2 > max_dist)
max_dist = d2; max_dist = d2;
@ -119,8 +117,8 @@ ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typena
for (uint32_t i = 0; i < currentNgb; i++) for (uint32_t i = 0; i < currentNgb; i++)
{ {
outputValue += computeWValue(c, ngb[i], distances[i], fun); outputValue += computeWValue(c, *ngb[i], distances[i], fun);
weight += computeWValue(c, ngb[i], distances[i], interpolateOne); weight += computeWValue(c, *ngb[i], distances[i], interpolateOne);
} }
return (outputValue == 0) ? 0 : (outputValue / weight); return (outputValue == 0) ? 0 : (outputValue / weight);