Changeset 4280

Timestamp:

06/28/10 16:17:30 (15 years ago)

Author:

Mathieu Morlighem

Message:

fixed diagnosticvert

Location:

issm/trunk/src/c

Files:

• modules/ModelProcessorx/DiagnosticVert/UpdateElementsDiagnosticVert.cpp (modified) (2 diffs)
• objects/Elements/Penta.cpp (modified) (3 diffs)
• objects/Elements/Penta.h (modified) (1 diff)

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/UpdateElementsDiagnosticVert.cpp

@@ -32 +32 @@
         IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
         IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
-        IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
-        IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
+        IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
+        IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
         IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
         IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
@@ -59 +59 @@
         xfree((void**)&iomodel->elementonsurface);
         xfree((void**)&iomodel->elementonwater);
+        xfree((void**)&iomodel->melting_rate);
+        xfree((void**)&iomodel->accumulation_rate);
         xfree((void**)&iomodel->vx);
         xfree((void**)&iomodel->vy);

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

@@ -318 +318 @@
                 InputUpdateFromSolutionDiagnosticHoriz( solution);
         }
+        else if (analysis_type==DiagnosticVertAnalysisEnum){
+                InputUpdateFromSolutionDiagnosticVert( solution);
+        }
         else if (analysis_type==DiagnosticStokesAnalysisEnum){
                 InputUpdateFromSolutionDiagnosticStokes( solution);
@@ -1575 +1578 @@
 
         /*Recover nodes ids needed to initialize the node hook.*/
-        for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook.
-                penta_node_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
+        for(i=0;i<6;i++){ 
+                //go recover node ids, needed to initialize the node hook.
+                //WARNING: We assume P1 elements here!!!!!
+                penta_node_ids[i]=iomodel->nodecounter+(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
         }
 
@@ -1943 +1948 @@
         inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
         
+        for(i=0;i<numvertices;i++){
+                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+        }
+        /*Now, we have to move the previous Vx and Vy inputs  to old 
+         * status, otherwise, we'll wipe them off: */
+        this->inputs->ChangeEnum(VxEnum,VxOldEnum);
+        this->inputs->ChangeEnum(VyEnum,VyOldEnum);
+        this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
+
+        /*Add vx and vy as inputs to the tria element: */
+        this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
+        this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
+        this->inputs->AddInput(new TriaVertexInput(VelEnum,vel));
+        this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
+}
+
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticVert {{{1*/
+void  Penta::InputUpdateFromSolutionDiagnosticVert(double* solution){
+
+        int i;
+
+        const int    numvertices=6;
+        const int    numdofpervertex=1;
+        const int    numdof=numdofpervertex*numvertices;
+
+        int          doflist[numdof];
+        double       values[numdof];
+        double       vx[numvertices];
+        double       vy[numvertices];
+        double       vz[numvertices];
+        double       vel[numvertices];
+        int          dummy;
+        double       pressure[numvertices];
+        double       surface[numvertices];
+        double       rho_ice,g;
+        double       xyz_list[numvertices][3];
+        double       gauss[numvertices][numvertices]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+
+        Input*       VxInput=NULL;
+        double*      VxPtr=NULL;
+        Input*       VyInput=NULL;
+        double*      VyPtr=NULL;
+
+        /*Get dof list: */
+        GetDofList(&doflist[0],&dummy);
+
+        /*Get node data: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+
+        /*Use the dof list to index into the solution vector: */
+        for(i=0;i<numdof;i++){
+                values[i]=solution[doflist[i]];
+        }
+
+        /*Ok, we have vz in values, fill in vz array: */
+        for(i=0;i<numvertices;i++){
+                vz[i]=values[i*numdofpervertex+0];
+        }
+
+        /*Get Vx and Vy*/
+        VxInput=inputs->GetInput(VxEnum);
+        if (VxInput){
+                if (VxInput->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vx is of type %s",EnumAsString(VxInput->Enum()));
+                VxInput->GetValuesPtr(&VxPtr,&dummy);
+                for(i=0;i<numvertices;i++) vx[i]=VxPtr[i];
+        }
+        else for(i=0;i<numvertices;i++) vx[i]=0.0;
+
+        VyInput=inputs->GetInput(VyEnum);
+        if (VyInput){
+                if (VyInput->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vy is of type %s",EnumAsString(VyInput->Enum()));
+                VyInput->GetValuesPtr(&VyPtr,&dummy);
+                for(i=0;i<numvertices;i++) vy[i]=VyPtr[i];
+        }
+        else for(i=0;i<numvertices;i++) vy[i]=0.0;
+
+        /*Now Compute vel*/
+        for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+
+        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
+         *so the pressure is just the pressure at the z elevation: */
+        rho_ice=matpar->GetRhoIce();
+        g=matpar->GetG();
+        inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
+
         for(i=0;i<numvertices;i++){
                 pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);

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

@@ -164 +164 @@
                 void      GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
                 void      GetStrainRateStokes(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
-                Penta*    GetUpperElement(void);
+                Penta*  GetUpperElement(void);
                 void      InputExtrude(int enum_type);
-                void      InputUpdateFromSolutionBalancedthickness( double* solutiong);
-                void      InputUpdateFromSolutionBalancedthickness2( double* solutiong);
-                void      InputUpdateFromSolutionBalancedvelocities( double* solutiong);
-                void      InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
-                void      InputUpdateFromSolutionDiagnosticStokes( double* solutiong);
-                void      InputUpdateFromSolutionPrognostic( double* solutiong);
-                void      InputUpdateFromSolutionPrognostic2(double* solutiong);
-                void      InputUpdateFromSolutionSlopeCompute( double* solutiong);
+                void    InputUpdateFromSolutionBalancedthickness( double* solutiong);
+                void    InputUpdateFromSolutionBalancedthickness2( double* solutiong);
+                void    InputUpdateFromSolutionBalancedvelocities( double* solutiong);
+                void    InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
+                void    InputUpdateFromSolutionDiagnosticVert( double* solutiong);
+                void    InputUpdateFromSolutionDiagnosticStokes( double* solutiong);
+                void    InputUpdateFromSolutionPrognostic( double* solutiong);
+                void    InputUpdateFromSolutionPrognostic2(double* solutiong);
+                void    InputUpdateFromSolutionSlopeCompute( double* solutiong);
                 bool      IsInput(int name);
                 bool      IsOnSurface(void);