← Back to branch summary

~ubuntu-branches/debian/jessie/yade/jessie

« back to all changes in this revision

Viewing changes to pkg/dem/Tetra.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2014-08-04 19:34:58 UTC
  • mfrom: (1.1.11)
  • Revision ID: package-import@ubuntu.com-20140804193458-cw8qhnujxe9wzi15
Tags: 1.11.0-1
* [a0600ae] Imported Upstream version 1.11.0
* [a3055e0] Do not use parallel build on kfreebsd-amd64 and s390x.
* [f86b405] Remove applied patches.

Show diffs side-by-side

added added

removed removed

Lines of Context:
pkg/dem/Tetra.cpp
659
659
                        #define v1 iABinfo[i].V[1]
660
660
                        #define v2 iABinfo[i].V[2]
661
661
                        #define v3 iABinfo[i].V[3]
662
 
                        Real dV=fabs(Vector3r(v1-v0).Dot((v2-v0).Cross(v3-v0)))/6.;
 
662
                        Real dV=std::abs(Vector3r(v1-v0).Dot((v2-v0).Cross(v3-v0)))/6.;
663
663
                        V+=dV;
664
664
                        Sg+=dV*(v0+v1+v2+v3)*.25;
665
665
                        vector<Vector3r> t; t.push_back(v0); t.push_back(v1); t.push_back(v2); t.push_back(v3);
956
956
        FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
957
957
                // normally, we would test isReal(), but TetraVolumetricLaw doesn't use phys at all
958
958
                if (!I->geom) continue; // Ig2_Tetra_Tetra_TTetraGeom::go returned false for this interaction, skip it
959
 
                const shared_ptr<TTetraGeom>& contactGeom(boost::dynamic_pointer_cast<TTetraGeom>(I->geom));
 
959
                const shared_ptr<TTetraGeom>& contactGeom(YADE_PTR_DYN_CAST<TTetraGeom>(I->geom));
960
960
                if(!contactGeom) continue;
961
961
 
962
962
                const Body::id_t idA=I->getId1(), idB=I->getId2();
963
963
                const shared_ptr<Body>& A=Body::byId(idA), B=Body::byId(idB);
964
964
                        
965
 
                const shared_ptr<ElastMat>& physA(boost::dynamic_pointer_cast<ElastMat>(A->material));
966
 
                const shared_ptr<ElastMat>& physB(boost::dynamic_pointer_cast<ElastMat>(B->material));
 
965
                const shared_ptr<ElastMat>& physA(YADE_PTR_DYN_CAST<ElastMat>(A->material));
 
966
                const shared_ptr<ElastMat>& physB(YADE_PTR_DYN_CAST<ElastMat>(B->material));
967
967
                
968
968
 
969
969
                /* Cross-section is volumetrically equivalent to the penetration configuration */
1073
1073
        // Jacobian of transformation to the reference 4hedron
1074
1074
        double detJ=(x2-x1)*(y3-y1)*(z4-z1)+(x3-x1)*(y4-y1)*(z2-z1)+(x4-x1)*(y2-y1)*(z3-z1)
1075
1075
                -(x2-x1)*(y4-y1)*(z3-z1)-(x3-x1)*(y2-y1)*(z4-z1)-(x4-x1)*(y3-y1)*(z2-z1);
1076
 
        detJ=fabs(detJ);
 
1076
        detJ=std::abs(detJ);
1077
1077
        double a=detJ*(y1*y1+y1*y2+y2*y2+y1*y3+y2*y3+
1078
1078
                y3*y3+y1*y4+y2*y4+y3*y4+y4*y4+z1*z1+z1*z2+
1079
1079
                z2*z2+z1*z3+z2*z3+z3*z3+z1*z4+z2*z4+z3*z4+z4*z4)/60.;
1143
1143
Quaternionr TetrahedronWithLocalAxesPrincipal(shared_ptr<Body>& tetraBody){
1144
1144
        //const shared_ptr<RigidBodyParameters>& rbp(YADE_PTR_CAST<RigidBodyParameters>(tetraBody->physicalParameters));
1145
1145
        State* rbp=tetraBody->state.get();
1146
 
        const shared_ptr<Tetra>& tMold(boost::dynamic_pointer_cast<Tetra>(tetraBody->shape));
 
1146
        const shared_ptr<Tetra>& tMold(YADE_PTR_DYN_CAST<Tetra>(tetraBody->shape));
1147
1147
 
1148
1148
        #define v0 tMold->v[0]
1149
1149
        #define v1 tMold->v[1]
1181
1181
 
1182
1182
 
1183
1183
Real TetrahedronSignedVolume(const Vector3r v[4]) { return (Vector3r(v[3])-Vector3r(v[0])).dot((Vector3r(v[3])-Vector3r(v[1])).cross(Vector3r(v[3])-Vector3r(v[2])))/6.; }
1184
 
Real TetrahedronVolume(const Vector3r v[4]) { return fabs(TetrahedronSignedVolume(v)); }
 
1184
Real TetrahedronVolume(const Vector3r v[4]) { return std::abs(TetrahedronSignedVolume(v)); }
1185
1185
Real TetrahedronSignedVolume(const vector<Vector3r>& v) { return Vector3r(v[1]-v[0]).dot(Vector3r(v[2]-v[0]).cross(v[3]-v[0]))/6.; }
1186
 
Real TetrahedronVolume(const vector<Vector3r>& v) { return fabs(TetrahedronSignedVolume(v)); }
 
1186
Real TetrahedronVolume(const vector<Vector3r>& v) { return std::abs(TetrahedronSignedVolume(v)); }
1187
1187
#ifdef YADE_CGAL
1188
1188
Real TetrahedronVolume(const CGAL::Point_3<CGAL::Cartesian<Real> >* v[4]) {
1189
1189
        Vector3r vv[4];
1208
1208
 
1209
1209
 
1210
1210
#ifdef YADE_CGAL
1211
 
void Law2_TTetraSimpleGeom_NormPhys_Simple::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
 
1211
bool Law2_TTetraSimpleGeom_NormPhys_Simple::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
1212
1212
        int id1 = contact->getId1(), id2 = contact->getId2();
1213
1213
        TTetraSimpleGeom* geom= static_cast<TTetraSimpleGeom*>(ig.get());
1214
1214
        NormPhys* phys = static_cast<NormPhys*>(ip.get());
1215
1215
        if ( geom->flag == 0 || geom->penetrationVolume <= 0. ) {
1216
 
                scene->interactions->requestErase(contact);
1217
 
                return;
 
1216
                return false;
1218
1217
        }
1219
1218
        Real& un=geom->penetrationVolume;
1220
1219
        phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
1223
1222
        State* de2 = Body::byId(id2,scene)->state.get();
1224
1223
        applyForceAtContactPoint(-phys->normalForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
1225
1224
        // TODO periodic
 
1225
        return true;
1226
1226
}
1227
1227
#endif
1228
1228