RepulsiveSpringForceField.inl
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef SOFA_COMPONENT_INTERACTIONFORCEFIELD_REPULSIVESPRINGFORCEFIELD_INL
00026 #define SOFA_COMPONENT_INTERACTIONFORCEFIELD_REPULSIVESPRINGFORCEFIELD_INL
00027
00028 #include <sofa/component/interactionforcefield/RepulsiveSpringForceField.h>
00029 #include <sofa/component/forcefield/StiffSpringForceField.inl>
00030
00031 namespace sofa
00032 {
00033
00034 namespace component
00035 {
00036
00037 namespace interactionforcefield
00038 {
00039 using std::cerr;
00040 using std::endl;
00041
00042 template<class DataTypes>
00043 void RepulsiveSpringForceField<DataTypes>::addForce(VecDeriv& f1, VecDeriv& f2, const VecCoord& x1, const VecCoord& x2, const VecDeriv& v1, const VecDeriv& v2)
00044 {
00045 const helper::vector<Spring>& springs= this->springs.getValue();
00046 this->dfdx.resize(springs.size());
00047 f1.resize(x1.size());
00048 f2.resize(x2.size());
00049 for (unsigned int i=0; i<springs.size(); i++)
00050 {
00051 #if 1
00052 int a = springs[i].m1;
00053 int b = springs[i].m2;
00054 Coord u = x2[b]-x1[a];
00055 Real d = u.norm();
00056 if (d < springs[i].initpos)
00057 {
00058 Real inverseLength = 1.0f/d;
00059 u *= inverseLength;
00060 Real elongation = (Real)(d - springs[i].initpos);
00061 Deriv relativeVelocity = v2[b]-v1[a];
00062 Real elongationVelocity = dot(u,relativeVelocity);
00063 Real forceIntensity = (Real)(springs[i].ks*elongation+springs[i].kd*elongationVelocity);
00064 Deriv force = u*forceIntensity;
00065 f1[a]+=force;
00066 f2[b]-=force;
00067
00068 Mat& m = this->dfdx[i];
00069 Real tgt = forceIntensity * inverseLength;
00070 for( int j=0; j<N; ++j )
00071 {
00072 for( int k=0; k<N; ++k )
00073 {
00074 m[j][k] = ((Real)springs[i].ks-tgt) * u[j] * u[k];
00075 }
00076 m[j][j] += tgt;
00077 }
00078 }
00079 else
00080 #endif
00081 {
00082 Mat& m = this->dfdx[i];
00083 for( int j=0; j<N; ++j )
00084 for( int k=0; k<N; ++k )
00085 m[j][k] = 0.0;
00086 }
00087 }
00088 }
00089
00090 template <class DataTypes>
00091 double RepulsiveSpringForceField<DataTypes>::getPotentialEnergy()
00092 {
00093 cerr<<"RepulsiveSpringForceField::getPotentialEnergy-not-implemented !!!"<<endl;
00094 return 0;
00095 }
00096
00097
00098 }
00099
00100 }
00101
00102 }
00103
00104 #endif