Changeset 5211


Ignore:
Timestamp:
08/12/10 15:57:40 (15 years ago)
Author:
seroussi
Message:

added InputUpdateFromSolutionMacAyealPattyn

Location:
issm/trunk/src/c/objects/Elements
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r5207 r5211  
    33213321                                        D_scalar=D_scalar*dt;
    33223322                                }
    3323                                 K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
    3324                                 K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
     3323                                K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*(u)*(v);
     3324                                K[1][0]=D_scalar*(u)*(v);        K[1][1]=D_scalar*pow(v,2);
     3325                                //K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
     3326                                //K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
    33253327
    33263328                                /*Get B_artdiff: */
     
    51915193                }
    51925194        }
    5193         else{
     5195        else if (approximation==PattynApproximationEnum){
    51945196                InputUpdateFromSolutionDiagnosticPattyn(solution);
     5197        }
     5198        else if (approximation==MacAyealPattynApproximationEnum){
     5199                InputUpdateFromSolutionDiagnosticMacAyealPattyn(solution);
    51955200        }
    51965201}
     
    52885293        /*Free ressources:*/
    52895294        xfree((void**)&doflist);
     5295}
     5296/*}}}*/
     5297/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn {{{1*/
     5298void  Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn(double* solution){
     5299
     5300        int i;
     5301
     5302        const int    numvertices=6;
     5303        const int    numdofpervertex=4;
     5304        const int    solutionnumdof=2;
     5305        const int    numdof=solutionnumdof*numvertices;
     5306        int*         doflist=NULL;
     5307        int*         doflistbase=NULL;
     5308        double       macayeal_values[numdof];
     5309        double       pattyn_values[numdof];
     5310        double       vx[numvertices];
     5311        double       vy[numvertices];
     5312        double       vz[numvertices];
     5313        double       vel[numvertices];
     5314        int          dummy;
     5315        double       pressure[numvertices];
     5316        double       surface[numvertices];
     5317        double       rho_ice,g;
     5318        double       xyz_list[numvertices][3];
     5319        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}};
     5320
     5321        Input  *vz_input       = NULL;
     5322        double *vz_ptr         = NULL;
     5323        Penta  *penta          = NULL;
     5324
     5325        /*OK, we have to add results of this element for pattyn
     5326         * and results from the penta at base for macayeal. Now recover results*/
     5327        penta=GetBasalElement();
     5328
     5329        /*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
     5330        GetDofList(&doflist);
     5331        penta->GetDofList(&doflistbase);
     5332
     5333        /*Get node data: */
     5334        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
     5335
     5336        /*Use the dof list to index into the solution vector: */
     5337        for(i=0;i<numdof;i++){
     5338                pattyn_values[i]=solution[doflist[i]];
     5339                macayeal_values[i]=solution[doflistbase[i]];
     5340        }
     5341
     5342        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     5343        for(i=0;i<numvertices;i++){
     5344                vx[i]=macayeal_values[i*numdofpervertex+0]+pattyn_values[i*numdofpervertex+0];
     5345                vy[i]=macayeal_values[i*numdofpervertex+1]+pattyn_values[i*numdofpervertex+1];
     5346        }
     5347
     5348        /*Get Vz*/
     5349        vz_input=inputs->GetInput(VzEnum);
     5350        if (vz_input){
     5351                if (vz_input->Enum()!=PentaVertexInputEnum){
     5352                        ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumToString(vz_input->Enum()));
     5353                }
     5354                vz_input->GetValuesPtr(&vz_ptr,&dummy);
     5355                for(i=0;i<numvertices;i++) vz[i]=vz_ptr[i];
     5356        }
     5357        else{
     5358                for(i=0;i<numvertices;i++) vz[i]=0.0;
     5359        }
     5360
     5361        /*Now Compute vel*/
     5362        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);
     5363
     5364        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
     5365         *so the pressure is just the pressure at the z elevation: */
     5366        rho_ice=matpar->GetRhoIce();
     5367        g=matpar->GetG();
     5368        inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
     5369
     5370        for(i=0;i<numvertices;i++){
     5371                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
     5372        }
     5373        /*Now, we have to move the previous Vx and Vy inputs  to old
     5374         * status, otherwise, we'll wipe them off: */
     5375        this->inputs->ChangeEnum(VxEnum,VxOldEnum);
     5376        this->inputs->ChangeEnum(VyEnum,VyOldEnum);
     5377        this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
     5378
     5379        /*Add vx and vy as inputs to the tria element: */
     5380        this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
     5381        this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
     5382        this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
     5383        this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
     5384
     5385        /*Free ressources:*/
     5386        xfree((void**)&doflist);
     5387        xfree((void**)&doflistbase);
    52905388}
    52915389/*}}}*/
  • issm/trunk/src/c/objects/Elements/Penta.h

    r5203 r5211  
    167167                void    InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
    168168                void    InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);
     169                void    InputUpdateFromSolutionDiagnosticMacAyealPattyn( double* solutiong);
    169170                void    InputUpdateFromSolutionDiagnosticPattyn( double* solutiong);
    170171                void    InputUpdateFromSolutionDiagnosticHutter( double* solutiong);
Note: See TracChangeset for help on using the changeset viewer.