Changeset 6150

Timestamp:

10/04/10 13:17:49 (14 years ago)

Author:

seroussi

Message:

added some code for MacAyealStokes

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (8 diffs)
• Penta.h (modified) (1 diff)

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

@@ -830 +830 @@
                         GetSolutionFromInputsDiagnosticHoriz(solution);
                 }
-                else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum){
+                else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
                         return; //the elements around will create the solution
                 }
@@ -4438 +4438 @@
                 InputUpdateFromSolutionDiagnosticPattynStokes(solution);
         }
+        else if (approximation==MacAyealStokesApproximationEnum){
+                InputUpdateFromSolutionDiagnosticMacAyealStokes(solution);
+        }
         else if (approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
                 InputUpdateFromSolutionDiagnosticStokes(solution);
@@ -4627 +4630 @@
         xfree((void**)&doflistp);
         xfree((void**)&doflistm);
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes {{{1*/
+void  Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes(double* solution){
+
+        int i;
+
+        const int    numdofm=NDOF2*NUMVERTICES;
+        const int    numdofs=NDOF4*NUMVERTICES;
+        const int    numdof2d=NDOF2*NUMVERTICES2D;
+        int*         doflistm=NULL;
+        int*         doflists=NULL;
+        double       macayeal_values[numdofm];
+        double       stokes_values[numdofs];
+        double       vx[NUMVERTICES];
+        double       vy[NUMVERTICES];
+        double       vz[NUMVERTICES];
+        double       vzmacayeal[NUMVERTICES];
+        double       vzstokes[NUMVERTICES];
+        double       vel[NUMVERTICES];
+        int          dummy;
+        double       pressure[NUMVERTICES];
+        double       xyz_list[NUMVERTICES][3];
+        double       stokesreconditioning;
+
+        double *vzmacayeal_ptr         = NULL;
+        Penta  *penta          = NULL;
+
+        /*OK, we have to add results of this element for macayeal 
+         * and results from the penta at base for macayeal. Now recover results*/
+        penta=GetBasalElement();
+
+        /*Get dof listof this element (macayeal dofs) and of the penta at base (macayeal dofs): */
+        penta->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum);
+        GetDofList(&doflists,StokesApproximationEnum,GsetEnum);
+        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
+
+        /*Get node data: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+
+        /*Use the dof list to index into the solution vector: */
+        for(i=0;i<numdof2d;i++){
+                macayeal_values[i]=solution[doflistm[i]];
+                macayeal_values[i+numdof2d]=solution[doflistm[i]];
+        }
+        for(i=0;i<numdofs;i++){
+                stokes_values[i]=solution[doflists[i]];
+        }
+
+        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        for(i=0;i<NUMVERTICES;i++){
+                vx[i]=stokes_values[i*NDOF4+0]+macayeal_values[i*NDOF2+0];
+                vy[i]=stokes_values[i*NDOF4+1]+macayeal_values[i*NDOF2+1];
+                vzstokes[i]=stokes_values[i*NDOF4+2];
+                pressure[i]=stokes_values[i*NDOF4+3]*stokesreconditioning;
+        }
+
+        /*Get Vz*/
+        Input* vzmacayeal_input=inputs->GetInput(VzMacAyealEnum);
+        if (vzmacayeal_input){
+                if (vzmacayeal_input->Enum()!=PentaVertexInputEnum){
+                        ISSMERROR("Cannot compute Vel as VzMacAyeal is of type %s",EnumToString(vzmacayeal_input->Enum()));
+                }
+                vzmacayeal_input->GetValuesPtr(&vzmacayeal_ptr,&dummy);
+                for(i=0;i<NUMVERTICES;i++) vzmacayeal[i]=vzmacayeal_ptr[i];
+        }
+        else{
+                ISSMERROR("Cannot update solution as VzMacAyeal is not present");
+        }
+
+        /*Now Compute vel*/
+        for(i=0;i<NUMVERTICES;i++) {
+                vz[i]=vzmacayeal[i]+vzstokes[i];
+                vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+        }
+
+        /*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(VzEnum,VzOldEnum);
+        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 PentaVertexInput(VzEnum,vz));
+        this->inputs->AddInput(new PentaVertexInput(VzStokesEnum,vzstokes));
+        this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
+        this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
+
+        /*Free ressources:*/
+        xfree((void**)&doflistm);
+        xfree((void**)&doflists);
 }
 /*}}}*/
@@ -4889 +4986 @@
         double       vy[NUMVERTICES];
         double       vz[NUMVERTICES];
+        double       vzmacayeal[NUMVERTICES];
         double       vzpattyn[NUMVERTICES];
         double       vzstokes[NUMVERTICES];
@@ -4943 +5041 @@
         else for(i=0;i<NUMVERTICES;i++) vy[i]=0.0;
 
-        /*Do some modifications if we actually have a PattynStokes element*/
+        /*Do some modifications if we actually have a PattynStokes or MacAyealStokes element*/
         if(approximation==PattynStokesApproximationEnum){
                 Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
@@ -4957 +5055 @@
                 }
         }
+        else if(approximation==MacAyealStokesApproximationEnum){
+                Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
+                if (vzstokes_input){
+                        if (vzstokes_input->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as VzStokes is of type %s",EnumToString(vy_input->Enum()));
+                        vzstokes_input->GetValuesPtr(&vzstokes_ptr,&dummy);
+                        for(i=0;i<NUMVERTICES;i++) vzstokes[i]=vzstokes_ptr[i];
+                }
+                else ISSMERROR("Cannot compute Vz as VzStokes in not present in MacAyealStokes element");
+                for(i=0;i<NUMVERTICES;i++){
+                        vzmacayeal[i]=vz[i];
+                        vz[i]=vzmacayeal[i]+vzstokes[i];
+                }
+        }
 
         /*Now Compute vel*/
@@ -4964 +5075 @@
         /*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: except it this is a PattynStokes element */
-        if(approximation!=PattynStokesApproximationEnum){
+        if(approximation!=PattynStokesApproximationEnum &&  approximation!=MacAyealStokesApproximationEnum){
                 rho_ice=matpar->GetRhoIce();
                 g=matpar->GetG();
@@ -4975 +5086 @@
         this->inputs->ChangeEnum(VzEnum,VzOldEnum);
 
-        if(approximation!=PattynStokesApproximationEnum){
+        if(approximation!=PattynStokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum){
                 this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
                 this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
         }
-        else{
+        else if(approximation==PattynStokesApproximationEnum){
                 this->inputs->AddInput(new PentaVertexInput(VzPattynEnum,vzpattyn));
+        }
+        else if(approximation==MacAyealStokesApproximationEnum){
+                this->inputs->AddInput(new PentaVertexInput(VzMacAyealEnum,vzpattyn));
         }
 

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

@@ -202 +202 @@
                 void    InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);
                 void    InputUpdateFromSolutionDiagnosticMacAyealPattyn( double* solutiong);
+                void    InputUpdateFromSolutionDiagnosticMacAyealStokes( double* solutiong);
                 void    InputUpdateFromSolutionDiagnosticPattyn( double* solutiong);
                 void    InputUpdateFromSolutionDiagnosticPattynStokes( double* solutiong);