Changeset 7672

Timestamp:

03/22/11 14:14:39 (14 years ago)

Author:

jschierm

Message:

Tria.cpp: fix scaled thicknesses for qmu baseline solution to be the same as a non-qmu solution (1/06/11).

File:

• issm/trunk/src/c/objects/Elements/Tria.cpp (modified) (3 diffs)

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -3854 +3854 @@
                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_vertex_ids[i]-1];
                 this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,nodeinputs));
+        }
+        if (iomodel->thickness_coeff) {
+                for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness_coeff[tria_vertex_ids[i]-1];
+                this->inputs->AddInput(new TriaVertexInput(ThicknessCoeffEnum,nodeinputs));
         }
         if (iomodel->surface) {
@@ -4432 +4436 @@
                                         /*Update thickness + surface: assume bed is constant. On ice shelves, takes hydrostatic equilibrium {{{2*/
                                         double  thickness[3];
+                                        double  thickness_init[3];
+                                        double  thickness_coeff[3];
                                         double  surface[3];
                                         double  bed[3];
                                         
                                         /*retrieve inputs: */
+                                        GetParameterListOnVertices(&thickness_init[0],ThicknessEnum);
+                                        GetParameterListOnVertices(&thickness_coeff[0],ThicknessCoeffEnum);
                                         GetParameterListOnVertices(&bed[0],BedEnum);
                                         GetParameterListOnVertices(&surface[0],SurfaceEnum);
 
                                         /*build new thickness: */
-                                        for(j=0;j<3;j++)thickness[j]=values[j];
+//                                      for(j=0;j<3;j++)thickness[j]=values[j];
+
+                    printf("Tria::InputUpdateFromVectorDakota,id=%d\n",this->Id());
+                    printf("nodes:            %d,%d,%d\n",this->nodes[0]->Sid()+1,this->nodes[1]->Sid()+1,this->nodes[2]->Sid())+1;
+                    printf("init thickness:   %e,%e,%e\n",thickness_init[0],thickness_init[1],thickness_init[2]);
+                    printf("thickness_coeff:  %e,%e,%e\n",thickness_coeff[0],thickness_coeff[1],thickness_coeff[2]);
+                    printf("values:           %e,%e,%e\n",values[0],values[1],values[2]);
+                    printf("old bed:          %e,%e,%e\n",bed[0],bed[1],bed[2]);
+                    printf("old surface:      %e,%e,%e\n",surface[0],surface[1],surface[2]);
 
                                         /*build new bed and surface: */
                                         if (this->IsOnShelf()){
+                        printf("shelf\n");
                                                 /*hydrostatic equilibrium: */
                                                 double rho_ice,rho_water,di;
@@ -4451 +4468 @@
                                                 di=rho_ice/rho_water;
 
-                                                for(j=0;j<3;j++){
-                                                        surface[j]=(1-di)*thickness[j];
-                                                        bed[j]=-di*thickness[j];
+                                                /*build new thickness: */
+                                                for (j=0; j<3; j++) {
+                                                /*  for observed/interpolated/hydrostatic thickness, remove scaling from any hydrostatic thickness  */
+                                                        if     (thickness_coeff[j] >= 0.)
+                                                                thickness[j]=values[j]-(values[j]/thickness_init[j]-1.)*thickness_coeff[j]*surface[j]/(1.-di);
+                                                /*  for minimum thickness, don't scale  */
+                                                        else
+                                                                thickness[j]=thickness_init[j];
+
+                                                /*  check the computed thickness and update bed  */
+                                                        if (thickness[j] < 0.)
+                                                                thickness[j]=1./(1.-di);
+                                                        bed[j]=surface[j]-thickness[j];
                                                 }
+
+//                                              for(j=0;j<3;j++){
+//                                                      surface[j]=(1-di)*thickness[j];
+//                                                      bed[j]=-di*thickness[j];
+//                                              }
                                         }
                                         else{
+                        printf("sheet\n");
+                                                /*build new thickness: */
+                                                for (j=0; j<3; j++) {
+                                                /*  for observed thickness, use scaled value  */
+                                                        if     (thickness_coeff[j] >= 0.)
+                                                                thickness[j]=values[j];
+                                                /*  for minimum thickness, don't scale  */
+                                                        else
+                                                                thickness[j]=thickness_init[j];
+                                                }
+
                                                 /*update bed on grounded ice: */
-                                                for(j=0;j<3;j++)surface[j]=bed[j]+thickness[j];
+//                                              for(j=0;j<3;j++)surface[j]=bed[j]+thickness[j];
+                                                for(j=0;j<3;j++)bed[j]=surface[j]-thickness[j];
                                         }
+
+                    printf("new thickness:    %e,%e,%e\n",thickness[0],thickness[1],thickness[2]);
+                    printf("new bed:          %e,%e,%e\n",bed[0],bed[1],bed[2]);
+                    printf("new surface:      %e,%e,%e\n",surface[0],surface[1],surface[2]);
 
                                         /*Add new inputs: */