Ogl3

EulerIntegrator.cpp

@@ -0,0 +1,203 @@
+#include "EulerIntegrator.h"
+
+//##ModelId=45F4D79800AE
+EulerIntegrator::EulerIntegrator(Object& object) : Integrator(object)// constructor
+{
+}
+
+//##ModelId=45F4D79800B0
+EulerIntegrator::~EulerIntegrator() // destructor
+{
+}
+
+
+// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+/***********************
+ * Euler Integrator integrate all the particles on the object *
+ ***********************/
+//void EulerIntegrator::Derivatives(Particle[] points, int NUMP, float DT)
+//void EulerIntegrator::integrate(Object object, int NUMP)
+
+//##ModelId=45F4D79800B2
+void EulerIntegrator::Derivatives(float deltaT, float k)
+{
+
+    int NUMP = object->GetNumberOfParticles();
+
+//    cout<<"derivative,,,,,,,,,,,,,,,,,,object->GetNumberOfParticles======"<<object->GetNumberOfParticles()<<endl;
+    for(int i=0; i<NUMP; i++)
+    {
+        ynew(i, k, deltaT);
+        CollisionDetection(i);
+    }
+}
+
+
+//##ModelId=45F4D79800BE
+void EulerIntegrator::k1(int i, float k, float deltaT)
+{
+
+    for(int j=0; j<object->outer_points.size(); j++)
+    {
+
+        temp_outer_points0[j]->mass = object->outer_points[j]->mass;
+        temp_outer_points0[j]->r = object->outer_points[j]->r;
+        temp_outer_points0[j]->v = object->outer_points[j]->v;
+        temp_outer_points0[j]->f = object->outer_points[j]->f;
+        temp_outer_points0[j]->dr = object->outer_points[j]->dr;
+        temp_outer_points0[j]->dv = object->outer_points[j]->dv;
+        temp_outer_points0[j]->OneOverMass = object->outer_points[j]->OneOverMass;
+
+        temp_inner_points0[j]->mass = object->inner_points[j]->mass;
+        temp_inner_points0[j]->r = object->inner_points[j]->r;
+        temp_inner_points0[j]->v = object->inner_points[j]->v;
+        temp_inner_points0[j]->f = object->inner_points[j]->f;
+        temp_inner_points0[j]->dr = object->inner_points[j]->dr;
+        temp_inner_points0[j]->dv = object->inner_points[j]->dv;
+        temp_inner_points0[j]->OneOverMass = object->inner_points[j]->OneOverMass;
+    }
+
+
+    temp_outer_points1[i]->dv->x = k * (temp_outer_points0[i]->f->x/temp_outer_points0[i]->mass) * deltaT;
+    temp_outer_points1[i]->dv->y = k * (temp_outer_points0[i]->f->y/temp_outer_points0[i]->mass) * deltaT;
+    temp_outer_points1[i]->dv->z = k * (temp_outer_points0[i]->f->z/temp_outer_points0[i]->mass) * deltaT;
+
+    temp_outer_points1[i]->dr->x = k * temp_outer_points0[i]->v->x * deltaT;
+    temp_outer_points1[i]->dr->y = k * temp_outer_points0[i]->v->y * deltaT;
+    temp_outer_points1[i]->dr->z = k * temp_outer_points0[i]->v->z * deltaT;
+
+
+    temp_inner_points1[i]->dv->x = k * (temp_inner_points0[i]->f->x/temp_inner_points0[i]->mass) * deltaT;
+    temp_inner_points1[i]->dv->y = k * (temp_inner_points0[i]->f->y/temp_inner_points0[i]->mass) * deltaT;
+    temp_inner_points1[i]->dv->z = k * (temp_inner_points0[i]->f->z/temp_inner_points0[i]->mass) * deltaT;
+
+    temp_inner_points1[i]->dr->x = k * temp_inner_points0[i]->v->x * deltaT;
+    temp_inner_points1[i]->dr->y = k * temp_inner_points0[i]->v->y * deltaT;
+    temp_inner_points1[i]->dr->z = k * temp_inner_points0[i]->v->z * deltaT;
+
+
+
+
+}
+
+//##ModelId=45F4D79800C2
+void EulerIntegrator::ynew(int i, float k, float deltaT)
+{
+    k1(i, k, deltaT);
+
+    /*if(object->dim==1){
+        if (i==0){
+    object->outer_points[i]->dv->x = 0;//temp_outer_points1[i]->dv->x;
+    object->outer_points[i]->dv->y = 0;//temp_outer_points1[i]->dv->y;
+    object->outer_points[i]->dv->z = 0;//temp_outer_points1[i]->dv->z;
+
+    object->outer_points[i]->dr->x = 0;//temp_outer_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->outer_points[i]->dr->y = 0;//temp_outer_points1[i]->dr->y;
+    object->outer_points[i]->dr->z = 0;//temp_outer_points1[i]->dr->z;
+
+    object->outer_points[i]->v->x = 0;//temp_outer_points0[i]->v->x + temp_outer_points1[i]->dv->x;                // Change in velocity is added to the velocity->
+    object->outer_points[i]->v->y = 0;//temp_outer_points0[i]->v->y + temp_outer_points1[i]->dv->y;
+    object->outer_points[i]->v->z = 0;//temp_outer_points0[i]->v->z + temp_outer_points1[i]->dv->z;
+
+    object->outer_points[i]->r->x = 0;//temp_outer_points0[i]->r->x + temp_outer_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->outer_points[i]->r->y = 0;//temp_outer_points0[i]->r->y + temp_outer_points1[i]->dr->y;
+    object->outer_points[i]->r->z = 0;//temp_outer_points0[i]->r->z + temp_outer_points1[i]->dr->z;
+
+
+
+    object->inner_points[i]->dv->x = 0;//temp_inner_points1[i]->dv->x;                                             // Change in position is velocity times the change in time
+    object->inner_points[i]->dv->y = 0;//temp_inner_points1[i]->dv->y;
+    object->inner_points[i]->dv->z = 0;//temp_inner_points1[i]->dv->z;
+
+    object->inner_points[i]->dr->x = 0;//temp_inner_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->inner_points[i]->dr->y = 0;//temp_inner_points1[i]->dr->y;
+    object->inner_points[i]->dr->z = 0;//temp_inner_points1[i]->dr->z;
+
+
+
+    object->inner_points[i]->v->x = 0;//temp_inner_points0[i]->v->x + temp_inner_points1[i]->dv->x;                // Change in velocity is added to the velocity->
+    object->inner_points[i]->v->y = 0;//temp_inner_points0[i]->v->y + temp_inner_points1[i]->dv->y;
+    object->inner_points[i]->v->z = 0;//temp_inner_points0[i]->v->z + temp_inner_points1[i]->dv->z;
+
+    object->inner_points[i]->r->x = 0;//temp_inner_points0[i]->r->x + temp_inner_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->inner_points[i]->r->y = 0;//temp_inner_points0[i]->r->y + temp_inner_points1[i]->dr->y;
+    object->inner_points[i]->r->z = 0;//temp_inner_points0[i]->r->z + temp_inner_points1[i]->dr->z;
+    }
+        else{
+        object->outer_points[i]->dv->x = temp_outer_points1[i]->dv->x;
+    object->outer_points[i]->dv->y = temp_outer_points1[i]->dv->y;
+    object->outer_points[i]->dv->z = temp_outer_points1[i]->dv->z;
+
+    object->outer_points[i]->dr->x = temp_outer_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->outer_points[i]->dr->y = temp_outer_points1[i]->dr->y;
+    object->outer_points[i]->dr->z = temp_outer_points1[i]->dr->z;
+
+    object->outer_points[i]->v->x = temp_outer_points0[i]->v->x + temp_outer_points1[i]->dv->x;                // Change in velocity is added to the velocity->
+    object->outer_points[i]->v->y = temp_outer_points0[i]->v->y + temp_outer_points1[i]->dv->y;
+    object->outer_points[i]->v->z = temp_outer_points0[i]->v->z + temp_outer_points1[i]->dv->z;
+
+    object->outer_points[i]->r->x = temp_outer_points0[i]->r->x + temp_outer_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->outer_points[i]->r->y = temp_outer_points0[i]->r->y + temp_outer_points1[i]->dr->y;
+    object->outer_points[i]->r->z = temp_outer_points0[i]->r->z + temp_outer_points1[i]->dr->z;
+
+
+
+    object->inner_points[i]->dv->x = temp_inner_points1[i]->dv->x;                                             // Change in position is velocity times the change in time
+    object->inner_points[i]->dv->y = temp_inner_points1[i]->dv->y;
+    object->inner_points[i]->dv->z = temp_inner_points1[i]->dv->z;
+
+    object->inner_points[i]->dr->x = temp_inner_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->inner_points[i]->dr->y = temp_inner_points1[i]->dr->y;
+    object->inner_points[i]->dr->z = temp_inner_points1[i]->dr->z;
+
+
+
+    object->inner_points[i]->v->x = temp_inner_points0[i]->v->x + temp_inner_points1[i]->dv->x;                // Change in velocity is added to the velocity->
+    object->inner_points[i]->v->y = temp_inner_points0[i]->v->y + temp_inner_points1[i]->dv->y;
+    object->inner_points[i]->v->z = temp_inner_points0[i]->v->z + temp_inner_points1[i]->dv->z;
+
+    object->inner_points[i]->r->x = temp_inner_points0[i]->r->x + temp_inner_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->inner_points[i]->r->y = temp_inner_points0[i]->r->y + temp_inner_points1[i]->dr->y;
+    object->inner_points[i]->r->z = temp_inner_points0[i]->r->z + temp_inner_points1[i]->dr->z;
+        }
+    }
+
+    else{*/
+    object->outer_points[i]->dv->x = temp_outer_points1[i]->dv->x;
+    object->outer_points[i]->dv->y = temp_outer_points1[i]->dv->y;
+    object->outer_points[i]->dv->z = temp_outer_points1[i]->dv->z;
+
+    object->outer_points[i]->dr->x = temp_outer_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->outer_points[i]->dr->y = temp_outer_points1[i]->dr->y;
+    object->outer_points[i]->dr->z = temp_outer_points1[i]->dr->z;
+
+    object->outer_points[i]->v->x = temp_outer_points0[i]->v->x + temp_outer_points1[i]->dv->x;                // Change in velocity is added to the velocity->
+    object->outer_points[i]->v->y = temp_outer_points0[i]->v->y + temp_outer_points1[i]->dv->y;
+    object->outer_points[i]->v->z = temp_outer_points0[i]->v->z + temp_outer_points1[i]->dv->z;
+
+    object->outer_points[i]->r->x = temp_outer_points0[i]->r->x + temp_outer_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->outer_points[i]->r->y = temp_outer_points0[i]->r->y + temp_outer_points1[i]->dr->y;
+    object->outer_points[i]->r->z = temp_outer_points0[i]->r->z + temp_outer_points1[i]->dr->z;
+
+
+
+    object->inner_points[i]->dv->x = temp_inner_points1[i]->dv->x;                                             // Change in position is velocity times the change in time
+    object->inner_points[i]->dv->y = temp_inner_points1[i]->dv->y;
+    object->inner_points[i]->dv->z = temp_inner_points1[i]->dv->z;
+
+    object->inner_points[i]->dr->x = temp_inner_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->inner_points[i]->dr->y = temp_inner_points1[i]->dr->y;
+    object->inner_points[i]->dr->z = temp_inner_points1[i]->dr->z;
+
+
+
+    object->inner_points[i]->v->x = temp_inner_points0[i]->v->x + temp_inner_points1[i]->dv->x;                // Change in velocity is added to the velocity->
+    object->inner_points[i]->v->y = temp_inner_points0[i]->v->y + temp_inner_points1[i]->dv->y;
+    object->inner_points[i]->v->z = temp_inner_points0[i]->v->z + temp_inner_points1[i]->dv->z;
+
+    object->inner_points[i]->r->x = temp_inner_points0[i]->r->x + temp_inner_points1[i]->dr->x;                                             // Change in position is velocity times the change in time
+    object->inner_points[i]->r->y = temp_inner_points0[i]->r->y + temp_inner_points1[i]->dr->y;
+    object->inner_points[i]->r->z = temp_inner_points0[i]->r->z + temp_inner_points1[i]->dr->z;
+//    }
+}
+

EulerIntegrator.h

@@ -0,0 +1,26 @@
+#ifndef EULERINTEGRATOR_H
+#define EULERINTEGRATOR_H
+
+#include "Integrator.h"
+
+
+class EulerIntegrator : public Integrator
+{
+public:
+
+    EulerIntegrator(Object&);        // constructor
+
+    virtual ~EulerIntegrator();      // destructor
+
+
+    virtual void Derivatives(float, float); // integrate all the particles on the object
+
+protected:
+
+    void k1(int, float, float);     // get the derivatives of velocity and position
+
+    void ynew(int, float, float);   // get the new velocity and position
+
+};
+
+#endif /* EULERINTEGRATOR_H */

Face.cpp

@@ -0,0 +1,29 @@
+#include "Face.h"
+#include<assert.h>
+
+// Class triface function definition
+// Calculate normal vector to a triangle and its field
+//##ModelId=45F4D7980066
+void Face::CalNormalNField(void)
+{
+    // the normal can be calculated by vector cross multiply
+
+/*    cout<<"fp1->r->x=="<<fp1->r->x<<"==fp1->r->y=="<<fp1->r->y<<"==fp1->r->z=="<<fp1->r->z<<endl;
+    cout<<"fp2->r->x=="<<fp2->r->x<<"==fp2->r->y=="<<fp2->r->y<<"==fp2->r->z=="<<fp2->r->z<<endl;
+    cout<<"fp3->r->x=="<<fp3->r->x<<"==fp3->r->y=="<<fp3->r->y<<"==fp3->r->z=="<<fp3->r->z<<endl;
+*/
+
+
+    Vector N = (*fp1->r - *fp2->r) ^ (*fp1->r - *fp3->r);
+//    Vector N = (*fs1->sp1->r - *fs2->sp1->r) ^ (*fs1->sp1->r - *fs3->sp1->r);
+
+
+//    cout<<"N.x=="<<N.x<<"==N.y=="<<N.y<<"==N.z=="<<N.z<<endl;
+    //field = N.getLength();      // triangle field, normal vector scalar
+    //assert(false);
+
+    N.Normalize();                                   // normalize the vector
+
+//    cout<<"after normalize==="<<"N.x=="<<N.x<<"==N.y=="<<N.y<<"==N.z=="<<N.z<<endl;
+    *normal = N;                                     // as the triangle face normal vector
+}

Face.h

@@ -0,0 +1,115 @@
+#ifndef FACE_H
+#define FACE_H
+
+#include "Particle.h"
+#include "Spring.h"
+
+// Class Face is an element cell of 3D ball entire surface
+//##ModelId=45F4D7980050
+class Face
+{
+public:
+
+    //face constructor with three particles point to the three points on each face, for 2D
+    //##ModelId=45F4D7980051
+    Face(Particle *Ap1, Particle *Ap2, Particle *Ap3) :
+           fp1(Ap1), fp2(Ap2), fp3(Ap3)
+    {
+        normal = new Vector(0,0,0);
+    }
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    // face constructor with three particles point to the three points on each face,
+    // and store the spring information into the spring vector for 3D
+    //##ModelId=45F4D7980060
+    Face(Particle *Ap1, Particle *Ap2, Particle *Ap3, vector<Spring*> &springs) :
+        fp1(Ap1), fp2(Ap2), fp3(Ap3)
+    {
+        normal = new Vector(0,0,0);                    // normal vecotr for face
+//        force  = new Vector(0,0,0);
+
+        fs1 = new Spring(Ap1, Ap2);                    // first spring points to particle p1&p2
+        fs2 = new Spring(Ap2, Ap3);                    // second spring points to particle p2&p3
+        fs3 = new Spring(Ap3, Ap1);                    // third spring points to particle p3&p1
+
+
+        bool a = false, b = false, c = false;            // the three new particles have not existed
+
+        for(int o = 0; o < springs.size(); o++)            // for all the spring in the requested spring vecotr
+        {
+            if                                // if the spring exist with particle p1&p2
+            (
+                   springs[o]->sp1 == Ap1
+                && springs[o]->sp2 == Ap2
+            )
+            {
+                   delete fs1;                    // delete the builded spring fs1
+                   fs1 = springs[o];                // the new spring will point to the existed spring
+                   a = true;                    // set the spring's particles have existed
+            }
+
+            if                                // if the spring exist with particle p2&p3
+            (
+                   springs[o]->sp1 == Ap2
+                && springs[o]->sp2 == Ap3
+            )
+            {
+                   delete fs2;                    // delete the builded spring fs2
+                   fs2 = springs[o];                // the new spring will point to the existed spring
+                    b = true;                    // set the spring's particles have existed
+            }
+
+            if                                // if the spring exist with particle p3&p1
+            (
+                   springs[o]->sp1 == Ap3
+                && springs[o]->sp2 == Ap1
+            )
+            {
+                   delete fs3;                    // delete the builded spring fs3
+                   fs3 = springs[o];                // the new spring will point to the existed spring
+                   c = true;                    // set the spring's particles have existed
+            }
+        }
+
+        // The new points are added to the general collection of points
+        if(!a) springs.push_back(fs1);
+        if(!b) springs.push_back(fs2);
+        if(!c) springs.push_back(fs3);
+    }
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+    //##ModelId=45F4D7980065
+    ~Face()                                    // face destructor
+    {
+//        delete fp1;
+//        delete fp2;
+//        delete fp3;
+        delete normal;
+//        delete force;
+    }
+///////////////////////////////////////////////////////////////////////////////////////////////////
+    //##ModelId=45F4D7980066
+    void CalNormalNField(void);                        // Calculate normal vector to a triangle and its field
+public:
+
+    //##ModelId=45F4D7980070
+    Particle     *fp1;         // first vertex
+    //##ModelId=45F4D7980075
+    Particle     *fp2;       // second vertex
+    //##ModelId=45F4D798007F
+    Particle     *fp3;       // third vertex*/
+
+    //##ModelId=45F4D7980084
+    Spring     *fs1;                             // first spring on the face
+    //##ModelId=45F4D7980089
+    Spring     *fs2;                               // second spring on the face
+    //##ModelId=45F4D798008F
+    Spring     *fs3;                               // third spring on the face
+//    float        field;                            // normal field
+    //##ModelId=45F4D7980094
+    Vector*     normal;                            // triangle normal Vector
+//    Vector*     force;                            // force on the triangle
+};
+
+#endif