Changeset 5316

Timestamp:

08/17/10 13:32:50 (15 years ago)

Author:

Eric.Larour

Message:

New routine to update dakota inputs in Tria.cpp

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (2 diffs)
• Tria.cpp (modified) (3 diffs)

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

@@ -414 +414 @@
         if (!IsInput(name)) return;
 
+        /*Penta update B in InputUpdateFromSolutionThermal, so don't look for B update here.*/
+
         switch(type){
 
@@ -448 +450 @@
 /*FUNCTION Penta::InputUpdateFromVectorDakota(double* vector, int name, int type);{{{1*/
 void  Penta::InputUpdateFromVectorDakota(double* vector, int name, int type){
-
-        /*Check that name is an element input*/
-        if (!IsInput(name)) return;
-
-        switch(type){
-
-                case VertexEnum:
-
-                        /*New PentaVertexInpu*/
-                        double values[6];
-
-                        /*Get values on the 6 vertices*/
-                        for (int i=0;i<6;i++){
-                                values[i]=vector[this->nodes[i]->GetVertexDof()];
-                        }
-
-                        /*update input*/
-                        this->inputs->AddInput(new PentaVertexInput(name,values));
-                        return;
-
-                default:
-
-                        ISSMERROR("type %i (%s) not implemented yet",type,EnumToString(type));
-        }
+        ISSMERROR(" not supported yet!");
 }
 /*}}}*/

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

@@ -413 +413 @@
                         double values[3];
 
-                        /*Get values on the 6 vertices*/
+                        /*Get values on the 3 vertices*/
                         for (int i=0;i<3;i++){
                                 values[i]=vector[this->nodes[i]->GetVertexDof()];
@@ -445 +445 @@
 /*FUNCTION Tria::InputUpdateFromVectorDakota(double* vector, int name, int type);{{{1*/
 void  Tria::InputUpdateFromVectorDakota(double* vector, int name, int type){
+        
+        const int numgrids=3;
+        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
+        int i,j;
 
         /*Check that name is an element input*/
@@ -456 +460 @@
                         double values[3];
 
-                        /*Get values on the 6 vertices*/
-                        for (int i=0;i<3;i++){
+                        /*Get values on the 3 vertices*/
+                        for (i=0;i<3;i++){
                                 values[i]=vector[this->nodes[i]->GetVertexDof()];
                         }
 
-                        /*update input*/
-                        if (name==RheologyB2dEnum || name==RheologyBEnum){
-                                matice->inputs->AddInput(new TriaVertexInput(name,values));
+                        /*Branch on the specified type of update: */
+                        switch(name){
+                                case RheologyB2dEnum:
+                                        matice->inputs->AddInput(new TriaVertexInput(name,values));
+                                        break;
+                                case RheologyBEnum:
+                                        matice->inputs->AddInput(new TriaVertexInput(name,values));
+                                        break;
+                                case ThicknessEnum:
+                                        /*Update thickness + surface: assume bed is constant. On ice shelves, takes hydrostatic equilibrium {{{2*/
+                                        double  thickness[3];
+                                        double  surface[3];
+                                        double  bed[3];
+                                        
+                                        /*retrieve inputs: */
+                                        inputs->GetParameterValues(&bed[0],&gaussgrids[0][0],3,BedEnum);
+                                        inputs->GetParameterValues(&surface[0],&gaussgrids[0][0],3,SurfaceEnum);
+
+                                        /*build new thickness: */
+                                        for(j=0;j<3;j++)thickness[j]=values[j];
+
+                                        /*build new bed and surface: */
+                                        if (this->GetShelf()){
+                                                /*hydrostatic equilibrium: */
+                                                double rho_ice,rho_water,di;
+                                                rho_ice=this->matpar->GetRhoIce();
+                                                rho_water=this->matpar->GetRhoWater();
+
+                                                di=rho_ice/rho_water;
+
+                                                for(j=0;j<3;j++){
+                                                        surface[j]=(1-di)*thickness[j];
+                                                        bed[j]=-di*thickness[j];
+                                                }
+                                        }
+                                        else{
+                                                /*update bed on grounded ice: */
+                                                for(j=0;j<3;j++)surface[j]=bed[j]+thickness[j];
+                                        }
+
+                                        /*Add new inputs: */
+                                        this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,thickness));
+                                        this->inputs->AddInput(new TriaVertexInput(BedEnum,bed));
+                                        this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,surface));
+
+                                        /*}}}*/
+                                        break;
+                                default:
+                                        this->inputs->AddInput(new TriaVertexInput(name,values));
                         }
-                        else{
-                                this->inputs->AddInput(new TriaVertexInput(name,values));
-                        }
-                        return;
+                        break;
 
                 default:
-
                         ISSMERROR("type %i (%s) not implemented yet",type,EnumToString(type));
         }
+
 }
 /*}}}*/