Changeset 5803


Ignore:
Timestamp:
09/14/10 14:46:19 (15 years ago)
Author:
seroussi
Message:

some modifications for PattynStokes coupling

Location:
issm/trunk/src/c
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r5770 r5803  
    252252        VzObsEnum,
    253253        VzOldEnum,
     254        VzPattynEnum,
     255        VzStokesEnum,
    254256        QmuVzEnum,
    255257        WeightsEnum,

  • issm/trunk/src/c/EnumDefinitions/EnumToString.cpp

    r5770 r5803  
    226226                case VzObsEnum : return "VzObs";
    227227                case VzOldEnum : return "VzOld";
     228                case VzPattynEnum : return "VzPattyn";
     229                case VzStokesEnum : return "VzStokes";
    228230                case QmuVzEnum : return "QmuVz";
    229231                case WeightsEnum : return "Weights";

  • issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp

    r5770 r5803  
    224224        else if (strcmp(name,"VzObs")==0) return VzObsEnum;
    225225        else if (strcmp(name,"VzOld")==0) return VzOldEnum;
     226        else if (strcmp(name,"VzPattyn")==0) return VzPattynEnum;
     227        else if (strcmp(name,"VzStokes")==0) return VzStokesEnum;
    226228        else if (strcmp(name,"QmuVz")==0) return QmuVzEnum;
    227229        else if (strcmp(name,"Weights")==0) return WeightsEnum;

  • issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp

    r4236 r5803  
    3131        /*Fetch data: */
    3232        IoModelFetchData(&iomodel->spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity");
     33        IoModelFetchData(&iomodel->gridonstokes,NULL,NULL,iomodel_handle,"gridonstokes");
    3334
    3435        /*Initialize counter*/
    … …  
    4142                if(iomodel->my_vertices[i]){
    4243
    43                         if ((int)iomodel->spcvelocity[6*i+2]){
     44                        if ((int)iomodel->gridonstokes[i]){
     45                                constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,DiagnosticVertAnalysisEnum)); //spc to zero as vertical velocity is done in Horiz for Stokes
     46                                count++;
     47                        }
     48                        else if ((int)iomodel->spcvelocity[6*i+2]){
    4449                                constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,
    4550                                                                *(iomodel->spcvelocity+6*i+5)/iomodel->yts,DiagnosticVertAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    … …  
    5257        /*Free data: */
    5358        xfree((void**)&iomodel->spcvelocity);
     59        xfree((void**)&iomodel->gridonstokes);
    5460
    5561        cleanup_and_return:

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

    r5785 r5803  
    367367        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
    368368        if (analysis_type==DiagnosticHorizAnalysisEnum){
    369                 int approximation;
    370                 inputs->GetParameterValue(&approximation,ApproximationEnum);
    371                 if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
    372                         InputUpdateFromSolutionDiagnosticStokes( solution);
    373                 }
    374                 else{
    375                         InputUpdateFromSolutionDiagnosticHoriz( solution);
    376                 }
     369                InputUpdateFromSolutionDiagnosticHoriz( solution);
    377370        }
    378371        else if (analysis_type==DiagnosticHutterAnalysisEnum){
    … …  
    901894                        GetSolutionFromInputsDiagnosticStokes(solution);
    902895                }
    903                 else{
     896                else if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum || approximation==HutterApproximationEnum){
    904897                        GetSolutionFromInputsDiagnosticHoriz(solution);
     898                }
     899                else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum){
     900                        return; //the elements around will create the solution
    905901                }
    906902        }
    … …  
    36033599        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
    36043600        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
    3605         Input* vz_input=NULL;
     3601        Input* vzstokes_input=NULL;
    36063602        if(approximation==PattynStokesApproximationEnum){
    3607                 vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input);
     3603                vzstokes_input=inputs->GetInput(VzStokesEnum); ISSMASSERT(vzstokes_input);
    36083604        }
    36093605
    … …  
    36183614                vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss);
    36193615                if(approximation==PattynStokesApproximationEnum){
    3620                         vz_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
     3616                        vzstokes_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    36213617                        dwdz=dw[2];
    36223618                }
    … …  
    40094005        /*If the element is a coupling, do nothing: every grid is also on an other elements
    40104006         * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/
    4011         if(approximation==MacAyealPattynApproximationEnum){
    4012          return;
    4013         }
    4014         else{
    4015                 GetDofList(&doflist,approximation,GsetEnum);
    4016         }
     4007        GetDofList(&doflist,approximation,GsetEnum);
    40174008
    40184009        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    44604451        else if (approximation==PattynApproximationEnum){
    44614452                InputUpdateFromSolutionDiagnosticPattyn(solution);
     4453        }
     4454        else if (approximation==PattynStokesApproximationEnum){
     4455                InputUpdateFromSolutionDiagnosticPattynStokes(solution);
     4456        }
     4457        else if (approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
     4458                InputUpdateFromSolutionDiagnosticStokes(solution);
    44624459        }
    44634460        else if (approximation==MacAyealPattynApproximationEnum){
    … …  
    47264723
    47274724/*}}}*/
     4725/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattynStokes {{{1*/
     4726void  Penta::InputUpdateFromSolutionDiagnosticPattynStokes(double* solution){
     4727
     4728        int i;
     4729
     4730        const int    numdofpervertexp=2;
     4731        const int    numdofpervertexs=4;
     4732        const int    numdofp=numdofpervertexp*NUMVERTICES;
     4733        const int    numdofs=numdofpervertexs*NUMVERTICES;
     4734        int*         doflistp=NULL;
     4735        int*         doflists=NULL;
     4736        double       pattyn_values[numdofp];
     4737        double       stokes_values[numdofs];
     4738        double       vx[NUMVERTICES];
     4739        double       vy[NUMVERTICES];
     4740        double       vz[NUMVERTICES];
     4741        double       vzpattyn[NUMVERTICES];
     4742        double       vzstokes[NUMVERTICES];
     4743        double       vel[NUMVERTICES];
     4744        int          dummy;
     4745        double       pressure[NUMVERTICES];
     4746        double       surface[NUMVERTICES];
     4747        double       rho_ice,g;
     4748        double       xyz_list[NUMVERTICES][3];
     4749
     4750        double *vzpattyn_ptr         = NULL;
     4751        Penta  *penta          = NULL;
     4752
     4753        /*OK, we have to add results of this element for pattyn
     4754         * and results from the penta at base for macayeal. Now recover results*/
     4755        penta=GetBasalElement();
     4756
     4757        /*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
     4758        GetDofList(&doflistp,PattynApproximationEnum,GsetEnum);
     4759        GetDofList(&doflists,StokesApproximationEnum,GsetEnum);
     4760
     4761        /*Get node data: */
     4762        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4763
     4764        /*Use the dof list to index into the solution vector: */
     4765        for(i=0;i<numdofp;i++){
     4766                pattyn_values[i]=solution[doflistp[i]];
     4767        }
     4768        for(i=0;i<numdofs;i++){
     4769                stokes_values[i]=solution[doflists[i]];
     4770        }
     4771
     4772        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     4773        for(i=0;i<NUMVERTICES;i++){
     4774                vx[i]=stokes_values[i*numdofpervertexs+0]+pattyn_values[i*numdofpervertexp+0];
     4775                vy[i]=stokes_values[i*numdofpervertexs+1]+pattyn_values[i*numdofpervertexp+1];
     4776                vzstokes[i]=stokes_values[i*numdofpervertexs+2];
     4777                pressure[i]=stokes_values[i*numdofpervertexs+3];
     4778        }
     4779
     4780        /*Get Vz*/
     4781        Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);
     4782        if (vzpattyn_input){
     4783                if (vzpattyn_input->Enum()!=PentaVertexInputEnum){
     4784                        ISSMERROR("Cannot compute Vel as VzPattyn is of type %s",EnumToString(vzpattyn_input->Enum()));
     4785                }
     4786                vzpattyn_input->GetValuesPtr(&vzpattyn_ptr,&dummy);
     4787                for(i=0;i<NUMVERTICES;i++) vzpattyn[i]=vzpattyn_ptr[i];
     4788        }
     4789        else{
     4790                ISSMERROR("Cannot update solution as VzPattyn is not present");
     4791        }
     4792
     4793        /*Now Compute vel*/
     4794        for(i=0;i<NUMVERTICES;i++) {
     4795                vz[i]=vzpattyn[i]+vzstokes[i];
     4796                vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     4797        }
     4798
     4799        /*Now, we have to move the previous Vx and Vy inputs  to old
     4800         * status, otherwise, we'll wipe them off: */
     4801        this->inputs->ChangeEnum(VxEnum,VxOldEnum);
     4802        this->inputs->ChangeEnum(VyEnum,VyOldEnum);
     4803        this->inputs->ChangeEnum(VzEnum,VzOldEnum);
     4804        this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
     4805
     4806        /*Add vx and vy as inputs to the tria element: */
     4807        this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
     4808        this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
     4809        this->inputs->AddInput(new PentaVertexInput(VzEnum,vz));
     4810        this->inputs->AddInput(new PentaVertexInput(VzStokesEnum,vzstokes));
     4811        this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
     4812        this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
     4813
     4814        /*Free ressources:*/
     4815        xfree((void**)&doflistp);
     4816        xfree((void**)&doflists);
     4817}
     4818/*}}}*/
    47284819/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHutter {{{1*/
    47294820void  Penta::InputUpdateFromSolutionDiagnosticHutter(double* solution){
    … …  
    48154906        double       vy[NUMVERTICES];
    48164907        double       vz[NUMVERTICES];
     4908        double       vzpattyn[NUMVERTICES];
     4909        double       vzstokes[NUMVERTICES];
    48174910        double       vel[NUMVERTICES];
    48184911        int          dummy;
    … …  
    48244917        double*      vx_ptr=NULL;
    48254918        double*      vy_ptr=NULL;
     4919        double*      vzstokes_ptr=NULL;
     4920
     4921        int         approximation;
    48264922
    48274923        /*Get dof list: */
    … …  
    48304926        /*Get node data: */
    48314927        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4928        inputs->GetParameterValue(&approximation,ApproximationEnum);
    48324929
    48334930        /*Use the dof list to index into the solution vector: */
    … …  
    48584955        else for(i=0;i<NUMVERTICES;i++) vy[i]=0.0;
    48594956
     4957        /*Do some modifications if we actually have a PattynStokes element*/
     4958        if(approximation==PattynStokesApproximationEnum){
     4959                Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
     4960                if (vzstokes_input){
     4961                        if (vzstokes_input->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as VzStokes is of type %s",EnumToString(vy_input->Enum()));
     4962                        vzstokes_input->GetValuesPtr(&vzstokes_ptr,&dummy);
     4963                        for(i=0;i<NUMVERTICES;i++) vzstokes[i]=vzstokes_ptr[i];
     4964                }
     4965                else ISSMERROR("Cannot compute Vz as VzStokes in not present in PattynStokes element");
     4966                for(i=0;i<NUMVERTICES;i++){
     4967                        vzpattyn[i]=vz[i];
     4968                        vz[i]=vzpattyn[i]+vzstokes[i];
     4969                }
     4970        }
     4971
    48604972        /*Now Compute vel*/
     4973
    48614974        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);
    48624975
    48634976        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
    4864          *so the pressure is just the pressure at the z elevation: */
    4865         rho_ice=matpar->GetRhoIce();
    4866         g=matpar->GetG();
    4867         GetParameterListOnVertices(&surface[0],SurfaceEnum);
    4868         for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
     4977         *so the pressure is just the pressure at the z elevation: except it this is a PattynStokes element */
     4978        if(!approximation==PattynStokesApproximationEnum){
     4979                rho_ice=matpar->GetRhoIce();
     4980                g=matpar->GetG();
     4981                GetParameterListOnVertices(&surface[0],SurfaceEnum);
     4982                for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
     4983        }
    48694984
    48704985        /*Now, we have to move the previous Vz inputs to old
    4871          * status, otherwise, we'll wipe them off: */
     4986         * status, otherwise, we'll wipe them off and add the new inputs: */
    48724987        this->inputs->ChangeEnum(VzEnum,VzOldEnum);
    4873         this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
    4874 
    4875         /*Add vz and vel as inputs to the penta element: */
     4988
     4989        if(!approximation==PattynStokesApproximationEnum){
     4990                this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
     4991                this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
     4992        }
     4993        else{
     4994                this->inputs->AddInput(new PentaVertexInput(VzPattynEnum,vzpattyn));
     4995        }
     4996
    48764997        this->inputs->AddInput(new PentaVertexInput(VzEnum,vz));
    48774998        this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
    4878         this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
    48794999
    48805000        /*Free ressources:*/

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

    r5785 r5803  
    41524152        Input* vx_input=inputs->GetInput(VxEnum);               ISSMASSERT(vx_input);
    41534153        Input* vy_input=inputs->GetInput(VyEnum);               ISSMASSERT(vy_input);
    4154         Input* vz_input=NULL;
     4154        Input* vzstokes_input=NULL;
    41554155        if(approximation==PattynStokesApproximationEnum){
    4156                 vz_input=inputs->GetInput(VzEnum);            ISSMASSERT(vz_input);
     4156                vzstokes_input=inputs->GetInput(VzStokesEnum);       ISSMASSERT(vzstokes_input);
    41574157        }
    41584158        Input* bed_input=inputs->GetInput(BedEnum);             ISSMASSERT(bed_input);
    … …  
    41724172                vy_input->GetParameterValue(&vy, gauss);
    41734173                if(approximation==PattynStokesApproximationEnum){
    4174                         vz_input->GetParameterValue(&vz, gauss);
     4174                        vzstokes_input->GetParameterValue(&vz, gauss);
    41754175                }
    41764176                else vz=0;
Note: See TracChangeset for help on using the changeset viewer.