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_FORCEFIELD_MESHSPRINGFORCEFIELD_INL
00026 #define SOFA_COMPONENT_FORCEFIELD_MESHSPRINGFORCEFIELD_INL
00027
00028 #include <sofa/component/forcefield/MeshSpringForceField.h>
00029 #include <sofa/component/forcefield/StiffSpringForceField.inl>
00030 #include <sofa/core/componentmodel/topology/BaseMeshTopology.h>
00031 #include <iostream>
00032 using std::cerr;
00033 using std::endl;
00034
00035 namespace sofa
00036 {
00037
00038 namespace component
00039 {
00040
00041 namespace forcefield
00042 {
00043
00044 template <class DataTypes>
00045 double MeshSpringForceField<DataTypes>::getPotentialEnergy()
00046 {
00047 cerr<<"MeshSpringForceField::getPotentialEnergy-not-implemented !!!"<<endl;
00048 return 0;
00049 }
00050
00051 template<class DataTypes>
00052 void MeshSpringForceField<DataTypes>::addSpring(std::set<std::pair<int,int> >& sset, int m1, int m2, Real stiffness, Real damping)
00053 {
00054 if (m1<m2)
00055 {
00056 if (sset.count(std::make_pair(m1,m2))>0) return;
00057 sset.insert(std::make_pair(m1,m2));
00058 }
00059 else
00060 {
00061 if (sset.count(std::make_pair(m2,m1))>0) return;
00062 sset.insert(std::make_pair(m2,m1));
00063 }
00064 Real l = ((*this->mstate2->getX0())[m2] - (*this->mstate1->getX0())[m1]).norm();
00065 this->springs.beginEdit()->push_back(typename SpringForceField<DataTypes>::Spring(m1,m2,stiffness/l, damping/l, l));
00066 this->springs.endEdit();
00067 }
00068
00069 template<class DataTypes>
00070 void MeshSpringForceField<DataTypes>::init()
00071 {
00072 this->StiffSpringForceField<DataTypes>::clear();
00073 if(!(this->mstate1) || !(this->mstate2))
00074 this->mstate2 = this->mstate1 = dynamic_cast<sofa::core::componentmodel::behavior::MechanicalState<DataTypes> *>(this->getContext()->getMechanicalState());
00075
00076 if (this->mstate1==this->mstate2)
00077 {
00078 sofa::core::componentmodel::topology::BaseMeshTopology* topology = this->mstate1->getContext()->getMeshTopology();
00079
00080
00081 if (topology != NULL)
00082 {
00083 std::set< std::pair<int,int> > sset;
00084 int n;
00085 Real s, d;
00086 if (this->linesStiffness.getValue() != 0.0 || this->linesDamping.getValue() != 0.0)
00087 {
00088 s = this->linesStiffness.getValue();
00089 d = this->linesDamping.getValue();
00090 n = topology->getNbLines();
00091 for (int i=0;i<n;++i)
00092 {
00093 sofa::core::componentmodel::topology::BaseMeshTopology::Line e = topology->getLine(i);
00094 this->addSpring(sset, e[0], e[1], s, d);
00095 }
00096 }
00097 if (this->trianglesStiffness.getValue() != 0.0 || this->trianglesDamping.getValue() != 0.0)
00098 {
00099 s = this->trianglesStiffness.getValue();
00100 d = this->trianglesDamping.getValue();
00101 n = topology->getNbTriangles();
00102 for (int i=0;i<n;++i)
00103 {
00104 sofa::core::componentmodel::topology::BaseMeshTopology::Triangle e = topology->getTriangle(i);
00105 this->addSpring(sset, e[0], e[1], s, d);
00106 this->addSpring(sset, e[0], e[2], s, d);
00107 this->addSpring(sset, e[1], e[2], s, d);
00108 }
00109 }
00110 if (this->quadsStiffness.getValue() != 0.0 || this->quadsDamping.getValue() != 0.0)
00111 {
00112 s = this->quadsStiffness.getValue();
00113 d = this->quadsDamping.getValue();
00114 n = topology->getNbQuads();
00115 for (int i=0;i<n;++i)
00116 {
00117 sofa::core::componentmodel::topology::BaseMeshTopology::Quad e = topology->getQuad(i);
00118 this->addSpring(sset, e[0], e[1], s, d);
00119 this->addSpring(sset, e[0], e[2], s, d);
00120 this->addSpring(sset, e[0], e[3], s, d);
00121 this->addSpring(sset, e[1], e[2], s, d);
00122 this->addSpring(sset, e[1], e[3], s, d);
00123 this->addSpring(sset, e[2], e[3], s, d);
00124 }
00125 }
00126 if (this->tetrasStiffness.getValue() != 0.0 || this->tetrasDamping.getValue() != 0.0)
00127 {
00128 s = this->tetrasStiffness.getValue();
00129 d = this->tetrasDamping.getValue();
00130 n = topology->getNbTetras();
00131 for (int i=0;i<n;++i)
00132 {
00133 sofa::core::componentmodel::topology::BaseMeshTopology::Tetra e = topology->getTetra(i);
00134 this->addSpring(sset, e[0], e[1], s, d);
00135 this->addSpring(sset, e[0], e[2], s, d);
00136 this->addSpring(sset, e[0], e[3], s, d);
00137 this->addSpring(sset, e[1], e[2], s, d);
00138 this->addSpring(sset, e[1], e[3], s, d);
00139 this->addSpring(sset, e[2], e[3], s, d);
00140 }
00141 }
00142
00143 if (this->cubesStiffness.getValue() != 0.0 || this->cubesDamping.getValue() != 0.0)
00144 {
00145 s = this->cubesStiffness.getValue();
00146 d = this->cubesDamping.getValue();
00147 #ifdef SOFA_NEW_HEXA
00148 n = topology->getNbHexas();
00149 for (int i=0;i<n;++i)
00150 {
00151 if (!topology->isCubeActive(i)) continue;
00152 sofa::core::componentmodel::topology::BaseMeshTopology::Hexa e = topology->getHexa(i);
00153 #else
00154 n = topology->getNbCubes();
00155 for (int i=0;i<n;++i)
00156 {
00157 if (!topology->isCubeActive(i)) continue;
00158 sofa::core::componentmodel::topology::BaseMeshTopology::Cube e = topology->getCube(i);
00159 #endif
00160 for (int i=0;i<8;i++)
00161 for (int j=i+1;j<8;j++)
00162 {
00163 this->addSpring(sset, e[i], e[j], s, d);
00164 }
00165 }
00166 }
00167 }
00168 }
00169 this->StiffSpringForceField<DataTypes>::init();
00170 }
00171
00172 template<class DataTypes>
00173 void MeshSpringForceField<DataTypes>::parse(core::objectmodel::BaseObjectDescription* arg)
00174 {
00175 this->StiffSpringForceField<DataTypes>::parse(arg);
00176 if (arg->getAttribute("stiffness")) this->setStiffness ((Real)atof(arg->getAttribute("stiffness")));
00177 if (arg->getAttribute("damping")) this->setDamping ((Real)atof(arg->getAttribute("damping")));
00178 }
00179
00180 }
00181
00182 }
00183
00184 }
00185
00186 #endif