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
00026
00027 #ifndef SOFA_CORE_COMPONENTMODEL_BEHAVIOR_MIXEDINTERACTIONFORCEFIELD_INL
00028 #define SOFA_CORE_COMPONENTMODEL_BEHAVIOR_MIXEDINTERACTIONFORCEFIELD_INL
00029
00030 #include <sofa/core/objectmodel/DataPtr.h>
00031 #include "MixedInteractionForceField.h"
00032
00033 namespace sofa
00034 {
00035
00036 namespace core
00037 {
00038
00039 namespace componentmodel
00040 {
00041
00042 namespace behavior
00043 {
00044
00045 template<class DataTypes1, class DataTypes2>
00046 MixedInteractionForceField<DataTypes1, DataTypes2>::MixedInteractionForceField(MechanicalState<DataTypes1> *mm1, MechanicalState<DataTypes2> *mm2)
00047 : mstate1(mm1), mstate2(mm2)
00048 {
00049 }
00050
00051 template<class DataTypes1, class DataTypes2>
00052 MixedInteractionForceField<DataTypes1, DataTypes2>::~MixedInteractionForceField()
00053 {
00054 }
00055
00056 template<class DataTypes1, class DataTypes2>
00057 void MixedInteractionForceField<DataTypes1, DataTypes2>::init()
00058 {
00059 InteractionForceField::init();
00060 }
00061
00062 template<class DataTypes1, class DataTypes2>
00063 void MixedInteractionForceField<DataTypes1, DataTypes2>::addForce()
00064 {
00065 if (mstate1 && mstate2)
00066 addForce(*mstate1->getF(), *mstate2->getF(),
00067 *mstate1->getX(), *mstate2->getX(),
00068 *mstate1->getV(), *mstate2->getV());
00069 }
00070
00071 template<class DataTypes1, class DataTypes2>
00072 void MixedInteractionForceField<DataTypes1, DataTypes2>::addDForce(double kFactor, double bFactor)
00073 {
00074 if (mstate1 && mstate2)
00075 addDForce(*mstate1->getF(), *mstate2->getF(),
00076 *mstate1->getDx(), *mstate2->getDx(),
00077 kFactor, bFactor);
00078 }
00079
00080 template<class DataTypes1, class DataTypes2>
00081 void MixedInteractionForceField<DataTypes1, DataTypes2>::addDForceV(double kFactor, double bFactor)
00082 {
00083 if (mstate1 && mstate2)
00084 addDForce(*mstate1->getF(), *mstate2->getF(),
00085 *mstate1->getV(), *mstate2->getV(),
00086 kFactor, bFactor);
00087 }
00088
00089 template<class DataTypes1, class DataTypes2>
00090 void MixedInteractionForceField<DataTypes1, DataTypes2>::addDForce(VecDeriv1& , VecDeriv2& , const VecDeriv1& , const VecDeriv2& )
00091 {
00092 std::cerr << "ERROR("<<getClassName()<<"): addDForce not implemented.\n";
00093 }
00094
00095 template<class DataTypes1, class DataTypes2>
00096 void MixedInteractionForceField<DataTypes1, DataTypes2>::addDForce(VecDeriv1& df1, VecDeriv2& df2, const VecDeriv1& dx1, const VecDeriv2& dx2, double kFactor, double )
00097 {
00098 if (kFactor == 1.0)
00099 addDForce(df1, df2, dx1, dx2);
00100 else if (kFactor != 0.0)
00101 {
00102 BaseMechanicalState::VecId vtmp1(BaseMechanicalState::VecId::V_DERIV,BaseMechanicalState::VecId::V_FIRST_DYNAMIC_INDEX);
00103 mstate1->vAvail(vtmp1);
00104 mstate1->vAlloc(vtmp1);
00105 BaseMechanicalState::VecId vdx1(BaseMechanicalState::VecId::V_DERIV,0);
00107 for (vdx1.index=0;vdx1.index<vtmp1.index;++vdx1.index)
00108 if (mstate1->getVecDeriv(vdx1.index) == &dx1)
00109 break;
00110 mstate1->vOp(vtmp1,BaseMechanicalState::VecId::null(),vdx1,kFactor);
00111
00112 BaseMechanicalState::VecId vtmp2(BaseMechanicalState::VecId::V_DERIV,BaseMechanicalState::VecId::V_FIRST_DYNAMIC_INDEX);
00113 mstate2->vAvail(vtmp2);
00114 mstate2->vAlloc(vtmp2);
00115 BaseMechanicalState::VecId vdx2(BaseMechanicalState::VecId::V_DERIV,0);
00117 for (vdx2.index=0;vdx2.index<vtmp2.index;++vdx2.index)
00118 if (mstate2->getVecDeriv(vdx2.index) == &dx2)
00119 break;
00120 mstate2->vOp(vtmp2,BaseMechanicalState::VecId::null(),vdx2,kFactor);
00121
00122 addDForce(df1, df2, *mstate1->getVecDeriv(vtmp1.index), *mstate2->getVecDeriv(vtmp2.index));
00123
00124 mstate1->vFree(vtmp1);
00125 mstate2->vFree(vtmp2);
00126 }
00127 }
00128
00129 template<class DataTypes1, class DataTypes2>
00130 double MixedInteractionForceField<DataTypes1, DataTypes2>::getPotentialEnergy()
00131 {
00132 if (mstate1 && mstate2)
00133 return getPotentialEnergy(*mstate1->getX(), *mstate2->getX());
00134 else return 0;
00135 }
00136
00137 }
00138
00139 }
00140
00141 }
00142
00143 }
00144
00145 #endif