Commit 3ce09f30 authored by Matthias Messner's avatar Matthias Messner

now forces with correct sign

parent 2a8f6a22
......@@ -158,7 +158,7 @@ private:
const FReal ws = Source.getPhysicalValue();
// potential
Target.incPotential(one_over_r * ws);
F3DPosition force(Target.getPosition() - Source.getPosition());
F3DPosition force(Source.getPosition() - Target.getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
// force
Target.incForces(force.getX(), force.getY(), force.getZ());
......
......@@ -391,7 +391,7 @@ inline void FChebInterpolator<ORDER>::applyL2PGradient(const F3DPosition& center
P[0] *= FReal(2.);
P[1] *= FReal(2.); P[1] += FReal(1.);
P[2] *= FReal(2.); P[2] += FReal(1.);
forces[0] -= P[0] * P[1] * P[2] * localExpansion[n];
forces[0] += P[0] * P[1] * P[2] * localExpansion[n];
// f1 component //////////////////////////////////////
P[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
......@@ -405,7 +405,7 @@ inline void FChebInterpolator<ORDER>::applyL2PGradient(const F3DPosition& center
P[0] *= FReal(2.); P[0] += FReal(1.);
P[1] *= FReal(2.);
P[2] *= FReal(2.); P[2] += FReal(1.);
forces[1] -= P[0] * P[1] * P[2] * localExpansion[n];
forces[1] += P[0] * P[1] * P[2] * localExpansion[n];
// f2 component //////////////////////////////////////
P[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
......@@ -419,7 +419,7 @@ inline void FChebInterpolator<ORDER>::applyL2PGradient(const F3DPosition& center
P[0] *= FReal(2.); P[0] += FReal(1.);
P[1] *= FReal(2.); P[1] += FReal(1.);
P[2] *= FReal(2.);
forces[2] -= P[0] * P[1] * P[2] * localExpansion[n];
forces[2] += P[0] * P[1] * P[2] * localExpansion[n];
}
// scale forces
......
......@@ -172,7 +172,7 @@ int main(int argc, char* argv[])
// potential
Potential[counter] += one_over_r * ws;
// force
F3DPosition force(iTarget.data().getPosition() - iSource.data().getPosition());
F3DPosition force(iSource.data().getPosition() - iTarget.data().getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
Force[counter*3 + 0] += force.getX();
Force[counter*3 + 1] += force.getY();
......
......@@ -140,7 +140,7 @@ public:
// potential
p[counter] += one_over_r * ws;
// force
F3DPosition force(iTarget.data().getPosition() - iSource.data().getPosition());
F3DPosition force(iSource.data().getPosition() - iTarget.data().getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
f[counter*3 + 0] += force.getX();
f[counter*3 + 1] += force.getY();
......@@ -181,16 +181,7 @@ int main(int argc, char* argv[])
// init timer
FTic time;
// // potentials
// FReal* p1; p1 = NULL;
// FReal* p2; p2 = NULL;
// FReal* p10; p10 = NULL; unsigned int nt1 = 0;
// FReal* p20; p20 = NULL; unsigned int nt2 = 0;
// // number of particles
// unsigned int NumParticles = 0;
// only for direct computation of nt1 target particles
unsigned int nt1 = 0;
FReal* p10; p10 = NULL;
FReal* f10; p10 = NULL;
......@@ -199,14 +190,11 @@ int main(int argc, char* argv[])
{ // begin Chebyshef kernel
// accuracy
const unsigned int ORDER = 7;
const FReal epsilon = FReal(1e-7);
const unsigned int ORDER = 5;
const FReal epsilon = FReal(1e-5);
//unsigned int nt1 = 0;
FReal* p1; p1 = NULL;
//FReal* p10; p10 = NULL;
FReal* f1; p1 = NULL;
//FReal* f10; p10 = NULL;
// typedefs
typedef FChebParticle ParticleClass;
......@@ -226,9 +214,6 @@ int main(int argc, char* argv[])
FFmaScanfLoader<ParticleClass> loader(filename);
if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");
// // loader set number of particles
// NumParticles = loader.getNumberOfParticles();
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
......@@ -374,8 +359,6 @@ int main(int argc, char* argv[])
} while(leavesIterator.moveRight());
} // -----------------------------------------------------
FBlas::scal(nt1*3, FReal(-1.), f10);
std::cout << "\nPotential error:" << std::endl;
std::cout << "Relative L2 error = "
<< computeL2norm( nt1, p10, fmmParticlesPotential) << std::endl;
......@@ -396,15 +379,6 @@ int main(int argc, char* argv[])
} // end FFmaBlas kernel
// analyse error
//std::cout << "\nRelative L2 error = " << computeL2norm( NumParticles, p1, p2) << std::endl;
//std::cout << "Relative Lmax error = " << computeINFnorm(NumParticles, p1, p2) << "\n" << std::endl;
// free memory
if (p10!=NULL) delete [] p10;
if (f10!=NULL) delete [] f10;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment