Changeset 16447

Timestamp:

10/17/13 15:22:40 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: fixing Updatevertex position for 2d vertical meshes

Location:

issm/trunk-jpl/src/c/classes

Files:

• FemModel.cpp (modified) (1 diff)
• Vertex.cpp (modified) (3 diffs)
• Vertex.h (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -369 +369 @@
 int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/
 
-        Vector<IssmDouble> *vz        = NULL;
-        IssmDouble         *thickness = NULL;
-        IssmDouble         *bed       = NULL;
+        IssmDouble         *surface = NULL;
+        IssmDouble         *bed     = NULL;
 
         /*get vertex vectors for bed and thickness: */
-        GetVectorFromInputsx(&thickness,this, ThicknessEnum,VertexEnum);
-        GetVectorFromInputsx(&bed      ,this, BedEnum,      VertexEnum);
+        GetVectorFromInputsx(&surface  ,this, SurfaceEnum,VertexEnum);
+        GetVectorFromInputsx(&bed      ,this, BedEnum,    VertexEnum);
 
         /*Allocate vector*/
-        vz=new Vector<IssmDouble>(vertices->NumberOfVertices());
+        Vector<IssmDouble> *vx=new Vector<IssmDouble>(vertices->NumberOfVertices());
+        Vector<IssmDouble> *vy=new Vector<IssmDouble>(vertices->NumberOfVertices());
+        Vector<IssmDouble> *vz=new Vector<IssmDouble>(vertices->NumberOfVertices());
 
         /*Update verices new geometry: */
         for(int i=0;i<vertices->Size();i++){
                 Vertex* vertex=(Vertex*)vertices->GetObjectByOffset(i);
-                vertex->UpdatePosition(vz,parameters,thickness,bed);
+                vertex->UpdatePosition(vx,vy,vz,parameters,surface,bed);
         }
 
         /*Assemble mesh velocity*/
+        vx->Assemble();
+        vy->Assemble();
         vz->Assemble();
 
         /*Update element inputs*/
+        InputUpdateFromVectorx(this,vx,VxMeshEnum,VertexPIdEnum);
+        InputUpdateFromVectorx(this,vy,VyMeshEnum,VertexPIdEnum);
         InputUpdateFromVectorx(this,vz,VzMeshEnum,VertexPIdEnum);
 
         /*Free ressources:*/
+        delete vx;
+        delete vy;
         delete vz;
         xDelete<IssmDouble>(bed);
-        xDelete<IssmDouble>(thickness);
+        xDelete<IssmDouble>(surface);
         return 1;
 }

issm/trunk-jpl/src/c/classes/Vertex.cpp

@@ -32 +32 @@
         this->y            = iomodel->Data(MeshYEnum)[i];
         this->z            = iomodel->Data(MeshZEnum)[i];
+        this->meshtype     = iomodel->meshtype;
 
         _assert_(iomodel->Data(BedEnum) && iomodel->Data(ThicknessEnum) && iomodel->numbernodetoelementconnectivity);
-        this->sigma        = (iomodel->Data(MeshZEnum)[i]-iomodel->Data(BedEnum)[i])/(iomodel->Data(ThicknessEnum)[i]);
+        switch(iomodel->meshtype){
+                case Mesh3DEnum:
+                case Mesh2DhorizontalEnum:
+                        this->sigma = (iomodel->Data(MeshZEnum)[i]-iomodel->Data(BedEnum)[i])/(iomodel->Data(ThicknessEnum)[i]);
+                        break;
+                case Mesh2DverticalEnum:
+                        this->sigma = (iomodel->Data(MeshYEnum)[i]-iomodel->Data(BedEnum)[i])/(iomodel->Data(ThicknessEnum)[i]);
+                        break;
+        }
+
         this->connectivity = iomodel->numbernodetoelementconnectivity[i];
 
@@ -112 +122 @@
 /*}}}*/
 /*FUNCTION Vertex::UpdatePosition {{{*/
-void  Vertex::UpdatePosition(Vector<IssmDouble>* vz,Parameters* parameters,IssmDouble* thickness,IssmDouble* bed){
-
-        IssmDouble oldz,newz;
-        IssmDouble dt,velz;
+void  Vertex::UpdatePosition(Vector<IssmDouble>* vx,Vector<IssmDouble>* vy,Vector<IssmDouble>* vz,Parameters* parameters,IssmDouble* surface,IssmDouble* bed){
+
+        IssmDouble oldy,newy,vely;
+        IssmDouble oldz,newz,velz;
+        IssmDouble dt;
 
         /*Get time stepping*/
@@ -121 +132 @@
 
         /*sigma remains constant. z=bed+sigma*thickness*/
-        oldz = this->z;
-        newz = bed[this->pid]+sigma*thickness[this->pid];
-        velz = (newz-oldz)/dt;
-        this->z = newz;
-
-        /*put vz in vector*/
-        vz->SetValue(this->pid,velz,INS_VAL);
+        if(this->meshtype==Mesh2DverticalEnum){
+                oldy = this->y;
+                newy = bed[this->pid]+sigma*(surface[this->pid] - bed[this->pid]);
+                vely = (newy-oldy)/dt;
+                this->y = newy;
+                vy->SetValue(this->pid,vely,INS_VAL);
+        }
+        else{
+                oldz = this->z;
+                newz = bed[this->pid]+sigma*(surface[this->pid] - bed[this->pid]);
+                velz = (newz-oldz)/dt;
+                this->z = newz;
+                vz->SetValue(this->pid,velz,INS_VAL);
+        }
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Vertex.h

@@ -21 +21 @@
         public: 
                 bool       clone;
+                int        meshtype;
                 int        id;           // random index
                 int        sid;          // "serial" id (rank of this vertex if the dataset was on 1 cpu)
@@ -49 +50 @@
                 IssmDouble GetY(void); 
                 IssmDouble GetZ(void); 
-                void       UpdatePosition(Vector<IssmDouble>* vz,Parameters* parameters,IssmDouble* thickness,IssmDouble* bed);
+                void       UpdatePosition(Vector<IssmDouble>* vx,Vector<IssmDouble>* vy,Vector<IssmDouble>* vz,Parameters* parameters,IssmDouble* thickness,IssmDouble* bed);
                 void       DistributePids(int* ppidcount);
                 void       OffsetPids(int pidcount);