Changeset 9009

Timestamp:

07/15/11 17:43:41 (14 years ago)

Author:

Eric.Larour

Message:

Implementing Dakota update in 3D

File:

• issm/trunk/src/c/objects/Elements/Penta.cpp (modified) (2 diffs)

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

@@ -4937 +4937 @@
                         this->inputs->AddInput(transientinput);
                 }
-                else _error_("nodal vector is either numberofnodes (%i), or numberofnodes+1 long. Field provided is %i long",iomodel->numberofvertices,M);
+                else _error_("nodal vector is either numberofnodes (%i), or numberofnodes+1 long. Field provided is %i long. Enum %s",iomodel->numberofvertices,M,EnumToStringx(vector_enum));
         }
         else if(vector_type==2){ //element vector
@@ -6505 +6505 @@
 /*FUNCTION Penta::InputUpdateFromVectorDakota(double* vector, int name, int type);{{{1*/
 void  Penta::InputUpdateFromVectorDakota(double* vector, int name, int type){
-        _error_(" not supported yet!");
+        
+        int i,j;
+
+        /*Check that name is an element input*/
+        if (!IsInput(name)) return;
+
+        switch(type){
+
+                case VertexEnum:
+
+                        /*New PentaVertexInput*/
+                        double values[6];
+
+                        /*Get values on the 6 vertices*/
+                        for (i=0;i<6;i++){
+                                values[i]=vector[this->nodes[i]->GetSidList()]; //careful, vector of values here is not parallel distributed, but serial distributed (from a serial Dakota core!)
+                        }
+
+                        /*Branch on the specified type of update: */
+                        switch(name){
+                                case ThicknessEnum:
+                                        /*Update thickness + surface: assume bed is constant. On ice shelves, takes hydrostatic equilibrium {{{2*/
+                                        double  thickness[6];
+                                        double  thickness_init[6];
+                                        double  thickness_coeff[6];
+                                        double  surface[6];
+                                        double  bed[6];
+                                        
+                                        /*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<6;j++)thickness[j]=values[j];
+
+                                        /*build new bed and surface: */
+                                        if (this->IsOnShelf()){
+                                                /*hydrostatic equilibrium: */
+                                                double rho_ice,rho_water,di;
+                                                rho_ice=this->matpar->GetRhoIce();
+                                                rho_water=this->matpar->GetRhoWater();
+
+                                                di=rho_ice/rho_water;
+
+                                                /*build new thickness: */
+                                                for (j=0; j<6; 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<6;j++){
+//                                                      surface[j]=(1-di)*thickness[j];
+//                                                      bed[j]=-di*thickness[j];
+//                                              }
+                                        }
+                                        else{
+                                                /*build new thickness: */
+                                                for (j=0; j<6; 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<6;j++)surface[j]=bed[j]+thickness[j];
+                                                for(j=0;j<6;j++)bed[j]=surface[j]-thickness[j];
+                                        }
+
+                                        /*Add new inputs: */
+                                        this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,thickness));
+                                        this->inputs->AddInput(new PentaVertexInput(BedEnum,bed));
+                                        this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,surface));
+
+                                        /*}}}*/
+                                        break;
+                                default:
+                                        this->inputs->AddInput(new PentaVertexInput(name,values));
+                        }
+                        break;
+
+                default:
+                        _error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
+        }
+
 }
 /*}}}*/