Want to take part in these discussions? Sign in if you have an account, or apply for one below
Vanilla 1.1.10 is a product of Lussumo. More Information: Documentation, Community Support.
//get neighborhood
unsigned int neighborhood = hashPhyxel(*i);
//calculate moment matrix A and inverse A^-1
Mat3 moment, invMoment;
//if there are too few neighbors (or if they're coplanar/colinear), the
//moment matrix will be singular (non-invertable)
//if(hashtable[neighborhood].size() < 4) {
// cout << "This is not a good moment: " << hashtable[neighborhood].size()
// << endl;
//}
for(list<Phyxel *>::iterator j = hashtable[neighborhood].begin();
j != hashtable[neighborhood].end(); j++) {
// for(list<Phyxel *>::iterator j = allPhyxels.begin();
// j != allPhyxels.end(); j++) {
if(*(*j) != *i) {
Vec3 diff = (*j)->getCurrentPosition() - i->getCurrentPosition();
Mat3 momentJ = ops::squareVector(diff);
momentJ *= weight(norm(diff),i->getSupportRadius());
moment += momentJ;
}
}
double aDet = invert(invMoment, moment);
if(aDet != 0 ) {
if(i->id == 499) {
//cout << "Phyxel #" << i->id << endl;
//cout << "Moment matrix: " << endl << moment << endl;
//cout << "Determinant: " << aDet << endl;
//cout << "Inverted moment: " << endl << invMoment << endl;
}
Vec3 sumU, sumV, sumW, sumF;
Vec3 dU, dV, dW;
//Equation 15: calculate derivatives
for(list<Phyxel *>::iterator j = hashtable[neighborhood].begin();
j != hashtable[neighborhood].end(); j++) {
//for(list<Phyxel *>::iterator j = allPhyxels.begin();
// j != allPhyxels.end(); j++) {
if(*(*j) != *i) {
Vec3 diff = (*j)->getCurrentPosition() - i->getCurrentPosition();
Vec3 spaceFactor = diff * weight(norm(diff),i->getSupportRadius());
double u = (*j)->getDisplacement()[0] - i->getDisplacement()[0];
sumU += u * spaceFactor;
double v = (*j)->getDisplacement()[1] - i->getDisplacement()[1];
sumV += v * spaceFactor;
double w = (*j)->getDisplacement()[2] - i->getDisplacement()[2];
sumW += w * spaceFactor;
}
}
//get x displacement derivative
dU = invMoment * sumU;
//get y displacement derivative
dV = invMoment * sumV;
//get z displacement derivative
dW = invMoment * sumW;
//assemble the Jacobian matrix J
Mat3 jacobian(dU, dV, dW);
jacobian += ops::IDENTITY;
//if(i->id == 499) {
// cout << "Jacobian: " << endl <<jacobian << endl;
//}
//calculate strain & stress tensors
Mat3 strain = transpose(jacobian) * jacobian - ops::IDENTITY;
Mat3 stress = strain * i->getModulusOfElasticity();
//calculate forces
//Eq 24-25: fixed term
Mat3 externalForce = -2 * i->getVolume() * jacobian * stress;
Mat3 conservation(cross(dV, dW), cross(dW, dU), cross(dU, dV));
Mat3 volumeForce = -(i->getVolume()) * k * (det(jacobian)-1) * conservation;
Mat3 fixedTerm = (externalForce + volumeForce) * invMoment;
//cout << i->id << " volume " << i->getVolume() << endl;
//cout << "External F: " << externalForce << endl;
//cout << "Volume F: " << volumeForce << endl;
//cout << "Force term: " << fixedTerm << endl;
for(list<Phyxel *>::iterator j = hashtable[neighborhood].begin();
j != hashtable[neighborhood].end(); j++) {
//for(list<Phyxel *>::iterator j = allPhyxels.begin();
// j != allPhyxels.end(); j++) {
if((*(*j)) != *i) {
Vec3 diff = (*j)->getCurrentPosition() - i->getCurrentPosition();
Vec3 spaceFactor = diff * weight(norm(diff),i->getSupportRadius());
sumF -= spaceFactor;
(*j)->addExternalForce( fixedTerm * spaceFactor);
}
}
i->addExternalForce(fixedTerm * -sumF);
//reset displacement
i->updatePlasticState();
} else {
//moment matrix couldn't be inverted
//Usually, this means the phyxel was kicked out of range of any
//neighbors. For now, I'm just trusting gravity to bring it back.
i->clearExternalForces();
i->updatePlasticState();
}
Posted By: gedavidshttp://dl.getdropbox.com/u/1932888/motivatordf4d70b3e402cc41ad05064cf381410bb054225a.jpgWow, my dropbox account!