Changeset 14363

Timestamp:

02/21/13 14:51:31 (12 years ago)

Author:

seroussi

Message:

FIX: fixed initialization of GLlevelset

Location:

issm/trunk-jpl/src/c

Files:

• classes/objects/Elements/Tria.cpp (modified) (2 diffs)
• modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp

@@ -1129 +1129 @@
         
         /*Check that not all nodes are grounded or floating*/
-        if(gl[0]>0 && gl[1]>0 && gl[2]>0) _error_("Error. Three nodes of the element are grounded, should not compute the fraction of grounded element");
-        if(gl[0]<0 && gl[1]<0 && gl[2]<0) _error_("Error. Three nodes of the element are floating, should not compute the fraction of grounded element");
-
-        /*Figure out if two nodes are floating or grounded*/
-        if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=true;
-
-        if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
-                /*Coordinates of point 2: same as initial point 2*/
-                xyz_bis[2][0]=*(xyz_list+3*2+0);
-                xyz_bis[2][1]=*(xyz_list+3*2+1);
-                xyz_bis[2][2]=*(xyz_list+3*2+2);
-
-                /*Portion of the segments*/
-                s1=gl[2]/(gl[2]-gl[1]);
-                s2=gl[2]/(gl[2]-gl[0]);
-                
-                /*New point 1*/
-                xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0));
-                xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1));
-                xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2));
-
-                /*New point 0*/
-                xyz_bis[0][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0));
-                xyz_bis[0][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1));
-                xyz_bis[0][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2));
-        }
-        else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
-                /*Coordinates of point 0: same as initial point 2*/
-                xyz_bis[0][0]=*(xyz_list+3*0+0);
-                xyz_bis[0][1]=*(xyz_list+3*0+1);
-                xyz_bis[0][2]=*(xyz_list+3*0+2);
-
-                /*Portion of the segments*/
-                s1=gl[0]/(gl[0]-gl[1]);
-                s2=gl[0]/(gl[0]-gl[2]);
-                
-                /*New point 1*/
-                xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0));
-                xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1));
-                xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2));
-
-                /*New point 2*/
-                xyz_bis[2][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0));
-                xyz_bis[2][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1));
-                xyz_bis[2][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*2+2)-*(xyz_list+3*0+2));
-        }
-        else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
-                /*Coordinates of point 1: same as initial point 2*/
-                xyz_bis[1][0]=*(xyz_list+3*1+0);
-                xyz_bis[1][1]=*(xyz_list+3*1+1);
-                xyz_bis[1][2]=*(xyz_list+3*1+2);
-
-                /*Portion of the segments*/
-                s1=gl[1]/(gl[1]-gl[0]);
-                s2=gl[1]/(gl[1]-gl[2]);
-                
-                /*New point 0*/
-                xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0));
-                xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1));
-                xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2));
-
-                /*New point 2*/
-                xyz_bis[2][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0));
-                xyz_bis[2][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1));
-                xyz_bis[2][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2));
-        }
-
-        /*Compute fraction of grounded element*/
-        GetJacobianDeterminant2d(&area_init, xyz_list,gauss);
-        GetJacobianDeterminant2d(&area_grounded, &xyz_bis[0][0],gauss);
-        if(mainlyfloating==true) area_grounded=area_init-area_grounded;
-        phi=area_grounded/area_init;
+        if(gl[0]>0 && gl[1]>0 && gl[2]>0){
+                phi=1;
+        }
+        else if(gl[0]<0 && gl[1]<0 && gl[2]<0){
+                _error_("Error. Three nodes of the element are floating, should not compute the fraction of grounded element");
+        }
+        else{
+
+                /*Figure out if two nodes are floating or grounded*/
+                if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=true;
+
+                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+                        /*Coordinates of point 2: same as initial point 2*/
+                        xyz_bis[2][0]=*(xyz_list+3*2+0);
+                        xyz_bis[2][1]=*(xyz_list+3*2+1);
+                        xyz_bis[2][2]=*(xyz_list+3*2+2);
+
+                        /*Portion of the segments*/
+                        s1=gl[2]/(gl[2]-gl[1]);
+                        s2=gl[2]/(gl[2]-gl[0]);
+
+                        /*New point 1*/
+                        xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0));
+                        xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1));
+                        xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2));
+
+                        /*New point 0*/
+                        xyz_bis[0][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0));
+                        xyz_bis[0][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1));
+                        xyz_bis[0][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2));
+                }
+                else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+                        /*Coordinates of point 0: same as initial point 2*/
+                        xyz_bis[0][0]=*(xyz_list+3*0+0);
+                        xyz_bis[0][1]=*(xyz_list+3*0+1);
+                        xyz_bis[0][2]=*(xyz_list+3*0+2);
+
+                        /*Portion of the segments*/
+                        s1=gl[0]/(gl[0]-gl[1]);
+                        s2=gl[0]/(gl[0]-gl[2]);
+
+                        /*New point 1*/
+                        xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0));
+                        xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1));
+                        xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2));
+
+                        /*New point 2*/
+                        xyz_bis[2][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0));
+                        xyz_bis[2][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1));
+                        xyz_bis[2][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*2+2)-*(xyz_list+3*0+2));
+                }
+                else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+                        /*Coordinates of point 1: same as initial point 2*/
+                        xyz_bis[1][0]=*(xyz_list+3*1+0);
+                        xyz_bis[1][1]=*(xyz_list+3*1+1);
+                        xyz_bis[1][2]=*(xyz_list+3*1+2);
+
+                        /*Portion of the segments*/
+                        s1=gl[1]/(gl[1]-gl[0]);
+                        s2=gl[1]/(gl[1]-gl[2]);
+
+                        /*New point 0*/
+                        xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0));
+                        xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1));
+                        xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2));
+
+                        /*New point 2*/
+                        xyz_bis[2][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0));
+                        xyz_bis[2][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1));
+                        xyz_bis[2][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2));
+                }
+
+                /*Compute fraction of grounded element*/
+                GetJacobianDeterminant2d(&area_init, xyz_list,gauss);
+                GetJacobianDeterminant2d(&area_grounded, &xyz_bis[0][0],gauss);
+                if(mainlyfloating==true) area_grounded=area_init-area_grounded;
+                phi=area_grounded/area_init;
+        }
+
         if(phi>1 || phi<0) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
+
         return phi;
 }
@@ -1985 +1993 @@
                                 name==ThicknessEnum ||
                                 name==SurfaceEnum ||
+                                name==BathymetryEnum ||
                                 name==BedEnum ||
+                                name==GLlevelsetEnum ||
                                 name==SurfaceSlopeXEnum ||
                                 name==SurfaceSlopeYEnum ||

issm/trunk-jpl/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp

@@ -29 +29 @@
 
                 /*get parameters and constants: */
-                parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
-                parameters->FindParam(&rho_water,MaterialsRhoWaterEnum);
+                iomodel->Constant(&rho_ice,MaterialsRhoIceEnum);
+                iomodel->Constant(&rho_water,MaterialsRhoWaterEnum);
                 iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
-                iomodel->FetchData(2,ThicknessEnum,BathymetryEnum);
+                iomodel->FetchData(3,MeshElementsEnum,ThicknessEnum,BathymetryEnum);
 
                 /*Create phi vector */
@@ -45 +45 @@
                 for(i=0;i<elements->Size();i++){
                         Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-                        element->InputUpdateFromVector(phi,GLlevelsetEnum,VertexEnum);
+                        element->InputCreate(phi,i,iomodel,numberofvertices,1,1,GLlevelsetEnum,1);
                 }
 
                 /*Free ressources:*/
+                iomodel->DeleteData(3,MeshElementsEnum,ThicknessEnum,BathymetryEnum);
                 xDelete<IssmDouble>(phi);
         }