Witam. Wzdęło nie na napisanie symulatora ośrodka sprężystego (link: https://github.com/pirogronian/ElasticBodySimulator - proszę wybaczyć obecne tam pliki tymczasowe, robiłem repo pierwszy raz i na szybko). Rzecz wiąże się z pewnymi średnio zaawansowanymi obliczeniami na wektorach 3D. Rzecz tyczy się (najpewniej) symulatora, czyli klasy ElasticBodySimulator, zawartej w plikach https://github.com/pirogronian/ElasticBodySimulator/blob/master/Simulator/v1.0/ElasticBodySimulator.cpp oraz https://github.com/pirogronian/ElasticBodySimulator/blob/master/Simulator/v1.0/ElasticBodySimulator.h. Chodzi o to, że w miarę trwania symulacji energia układu rośnie, prawdopodobnie wykładniczo. Na obecnym etapie (po wyłapaniu kilku ewidentnych błędów) nie mam już żadnej pewności, czy jest to dalej błąd w moim algorytmie, czy też wynik niedoskonałości obliczeń zmiennoprzecinkowych.
Aby samodzielnie dokonać testu symulacji, należy po uruchomieniu programu kliknąć (nie ma tam nic innego) zakładkę Symulacja, następnie Nowa symulacja. Domyślnie układ składa się z sześciennej kostki. Rzutowanie na płaszczyznę odbywa się na stałe w poprzek osi X na zerowym indeksie, więc rezultatem będą cztery kropki na planie kwadratu. Z dostępnych obecnie opcji nowej symulacji istotne są Okres czasu - odcinek czasu na krok - jego zwiększenie bardzo przyspiesza przyrost energii, więc mam wrażenie, jakbym gdzieś w kodzie dał o jedno mnożenie za dużo; reszta powinna być jasna. Masa wiązania jest obecnie nieużywana.
Szczególnie podejrzane dla mnie są funkcje symulatora, odpowiedzialne za wyliczanie wektorów odległości oraz siły i położenia:
QVector3D inline distances(ElasticBodySimulator::Granule a, ElasticBodySimulator::Granule b) {
QVector3D v;
v.setX(b.position.x() - a.position.x());
v.setY(b.position.y() - a.position.y());
v.setZ(b.position.z() - a.position.z());
return v;
}
QVector3D ElasticBodySimulator::forces(ElasticBodySimulator::Granule a, ElasticBodySimulator::Granule b) {
QVector3D dis, fp;
dis = distances(a, b);
qreal f = _couple * (dis.length() - _eqspace);
fp.setX(f * dis.x());
fp.setY(f * dis.y());
fp.setZ(f * dis.z());
return fp;
}
void ElasticBodySimulator::updateForce(int i, int j, int k) {
ElasticBodySimulator::Granule &g = _gMatrix[i][j][k];
g.force.setX(0);
g.force.setY(0);
g.force.setZ(0);
if (i < _xSize - 1) g.force += forces(g, _gMatrix[i+1][j][k]);
if (j < _ySize - 1) g.force += forces(g, _gMatrix[i][j+1][k]);
if (k < _zSize - 1) g.force += forces(g, _gMatrix[i][j][k+1]);
if (i > 0) g.force += forces(g, _gMatrix[i-1][j][k]);
if (j > 0) g.force += forces(g, _gMatrix[i][j-1][k]);
if (k > 0) g.force += forces(g, _gMatrix[i][j][k-1]);
}
void ElasticBodySimulator::updatePosition(ElasticBodySimulator::Granule &g) {
QVector3D accel, newspeed, movement;
accel = g.force / granMass();
newspeed = g.speed + accel * _timeSlice;
movement = (g.speed + newspeed) / 2 * _timeSlice;
g.position += movement;
g.speed = newspeed;
}
Przeglądałem je ju jednak tyle razy, że jeśli faktycznie jest tam błąd, to po prostu jestem na niego ślepy. Będę wdzięczny za wszelkie sugestie.
C++ 