Changeset 4832

Timestamp:

07/27/10 11:41:32 (15 years ago)

Author:

Mathieu Morlighem

Message:

Extrudion of Vx and Vy is now done by the elements

Location:

issm/trunk/src/c

Files:

• objects/Elements/Penta.cpp (modified) (5 diffs)
• objects/Elements/Penta.h (modified) (1 diff)
• solvers/solver_diagnostic_nonlinear.cpp (modified) (1 diff)

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

@@ -2125 +2125 @@
 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{1*/
 void  Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
+
+        /*Inputs*/
+        bool    collapsed,onbed;
+
+        /*Recover inputs*/
+        inputs->GetParameterValue(&collapsed,CollapseEnum);
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
+        /*MacAyeal, everything is done by the element on bed*/
+        if (collapsed){
+                if (!onbed){
+                        return;
+                }
+                else{
+                        InputUpdateFromSolutionDiagnosticMacAyeal(solution);
+                        return;
+                }
+        }
+        else{
+                InputUpdateFromSolutionDiagnosticPattyn(solution);
+        }
+}
+
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyeal {{{1*/
+void  Penta::InputUpdateFromSolutionDiagnosticMacAyeal(double* solution){
+
+        int i;
+
+        const int    numvertices=6;
+        const int    numdofpervertex=2;
+        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][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+
+        Input  *VzInput        = NULL;
+        double *VzPtr          = NULL;
+        Penta  *penta          = NULL;
+
+        /*OK, we are on bed. Now recover results*/
+
+        /*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 vx and vy in values, fill in vx and vy arrays and extrude */
+        for(i=0;i<3;i++){
+                vx[i]  =values[i*numdofpervertex+0];
+                vy[i]  =values[i*numdofpervertex+1];
+                vx[i+3]=vx[i];
+                vy[i+3]=vy[i];
+        }
+
+        /*Get parameters fro pressure computation*/
+        rho_ice=matpar->GetRhoIce();
+        g=matpar->GetG();
+
+        /*Start looping over all elements above current element and update all inputs*/
+        penta=this;
+        for(;;){
+
+                /*Now Compute vel*/
+                VzInput=inputs->GetInput(VzEnum);
+                if (VzInput){
+                        if (VzInput->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumAsString(VzInput->Enum()));
+                        VzInput->GetValuesPtr(&VzPtr,&dummy);
+                        for(i=0;i<numvertices;i++) vz[i]=VzPtr[i];
+                }
+                else{for(i=0;i<numvertices;i++) vz[i]=0.0;}
+                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);
+
+                /*Now compute pressure*/
+                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: */
+                penta->inputs->ChangeEnum(VxEnum,VxOldEnum);
+                penta->inputs->ChangeEnum(VyEnum,VyOldEnum);
+                penta->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
+
+                /*Add vx and vy as inputs to the tria element: */
+                penta->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
+                penta->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
+                penta->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
+                penta->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
+
+                /*Stop if we have reached the surface*/
+                if (penta->IsOnSurface()) break;
+
+                /* get upper Penta*/
+                penta=penta->GetUpperElement(); ISSMASSERT(penta->Id()!=this->id);
+        }
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattyn {{{1*/
+void  Penta::InputUpdateFromSolutionDiagnosticPattyn(double* solution){
         
         int i;
@@ -2145 +2262 @@
         double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
         
-        Input*       VzInput=NULL;
-        double*      VzPtr=NULL;
+        Input  *VzInput        = NULL;
+        double *VzPtr          = NULL;
 
         /*Get dof list: */
@@ -2186 +2303 @@
         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]);
@@ -2199 +2316 @@
         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(VelEnum,vel));
         this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
 }
@@ -4916 +5033 @@
 
         /*Are we on the base, not on the surface?:*/
-        if(onbed==1){
+        if(onbed){
 
                 /*OK, we are on bed. we will follow the steps:

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

@@ -177 +177 @@
                 void    InputUpdateFromSolutionBalancedvelocities( double* solutiong);
                 void    InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
+                void    InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);
+                void    InputUpdateFromSolutionDiagnosticPattyn( double* solutiong);
                 void    InputUpdateFromSolutionDiagnosticHutter( double* solutiong);
                 void    InputUpdateFromSolutionDiagnosticVert( double* solutiong);

issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp

@@ -128 +128 @@
         }
         
-
-        /*extrude if we are in 3D: */
-        if (dim==3){
-                if(verbose)_printf_("%s\n","extruding velocity of collapsed elements...");
-                InputExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,true); //true means that the velocity is only extruded if collapsed
-                InputExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,true);
-        }
-
         //more output might be needed, when running in control
         if(pKff0){