Changeset 8117


Ignore:
Timestamp:
05/03/11 09:08:29 (14 years ago)
Author:
Mathieu Morlighem
Message:

Added mesh velocity for thermal model. Still need to initialize it

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

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r8073 r8117  
    436436        RheologyLawEnum,
    437437        PatersonEnum,
    438         ArrheniusEnum
    439         /*}}}*/
     438        ArrheniusEnum,
     439        /*}}}*/
     440        VxMeshEnum,
     441        VyMeshEnum,
     442        VzMeshEnum
    440443};
    441444
  • issm/trunk/src/c/EnumDefinitions/EnumToString.cpp

    r8073 r8117  
    384384                case PatersonEnum : return "Paterson";
    385385                case ArrheniusEnum : return "Arrhenius";
     386                case VxMeshEnum : return "VxMesh";
     387                case VyMeshEnum : return "VyMesh";
     388                case VzMeshEnum : return "VzMesh";
    386389                default : return "unknown";
    387390
  • issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp

    r8073 r8117  
    382382        else if (strcmp(name,"Paterson")==0) return PatersonEnum;
    383383        else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
     384        else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
     385        else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
     386        else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
    384387        else _error_("Enum %s not found",name);
    385388
  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r8073 r8117  
    17031703        int        artdiff;
    17041704        int        i,j,ig,found=0;
    1705         double     Jdet,u,v,w,epsvel;
     1705        double     Jdet,u,v,w,um,vm,wm,epsvel;
    17061706        double     gravity,rho_ice,rho_water;
    17071707        double     heatcapacity,thermalconductivity,dt;
     
    17401740        this->parameters->FindParam(&artdiff,ArtDiffEnum);
    17411741        this->parameters->FindParam(&epsvel,EpsVelEnum);
    1742         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    1743         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    1744         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     1742        Input* vx_input=inputs->GetInput(VxEnum);      _assert_(vx_input);
     1743        Input* vy_input=inputs->GetInput(VyEnum);      _assert_(vy_input);
     1744        Input* vz_input=inputs->GetInput(VzEnum);      _assert_(vz_input);
     1745        Input* vxm_input=inputs->GetInput(VxMeshEnum); _assert_(vxm_input);
     1746        Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input);
     1747        Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input);
    17451748        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
    17461749
     
    17771780                vy_input->GetParameterValue(&v, gauss);
    17781781                vz_input->GetParameterValue(&w, gauss);
     1782                vxm_input->GetParameterValue(&um,gauss);
     1783                vym_input->GetParameterValue(&vm,gauss);
     1784                vzm_input->GetParameterValue(&wm,gauss);
    17791785
    17801786                D_scalar_advec=gauss->weight*Jdet;
    17811787                if(dt) D_scalar_advec=D_scalar_advec*dt;
    17821788
    1783                 D[0][0]=D_scalar_advec*u;D[0][1]=0;         D[0][2]=0;
    1784                 D[1][0]=0;         D[1][1]=D_scalar_advec*v;D[1][2]=0;
    1785                 D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar_advec*w;
     1789                D[0][0]=D_scalar_advec*(um-u);D[0][1]=0;                    D[0][2]=0;
     1790                D[1][0]=0;                    D[1][1]=D_scalar_advec*(v-vm);D[1][2]=0;
     1791                D[2][0]=0;                    D[2][1]=0;                    D[2][2]=D_scalar_advec*(w-wm);
    17861792
    17871793                TripleMultiply(&B_advec[0][0],3,numdof,1,
     
    18101816                if(artdiff==1){
    18111817                        /*Build K: */
    1812                         D_scalar_artdiff=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
     1818                        D_scalar_artdiff=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
    18131819                        if(dt) D_scalar_artdiff=D_scalar_artdiff*dt;
    18141820                        K[0][0]=D_scalar_artdiff*pow(u,2);       K[0][1]=D_scalar_artdiff*fabs(u)*fabs(v);
     
    18261832                        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
    18271833
    1828                         tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
     1834                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
    18291835
    18301836                        for(i=0;i<numdof;i++){
    18311837                                for(j=0;j<numdof;j++){
    1832                                         Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i])*(u*dh1dh6[0][j]+v*dh1dh6[1][j]+w*dh1dh6[2][j]);
     1838                                        Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*
     1839                                          ((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i])*((u-um)*dh1dh6[0][j]+(v-vm)*dh1dh6[1][j]+(w-wm)*dh1dh6[2][j]);
    18331840                                }
    18341841                        }
     
    18361843                                for(i=0;i<numdof;i++){
    18371844                                        for(j=0;j<numdof;j++){
    1838                                                 Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
     1845                                                Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i]);
    18391846                                        }
    18401847                                }
     
    62386245                        break;
    62396246
     6247                case ThermalAnalysisEnum:
     6248                        /*Initialize mesh velocity*/
     6249                        for(i=0;i<6;i++)nodeinputs[i]=0;
     6250                        this->inputs->AddInput(new PentaVertexInput(VxMeshEnum,nodeinputs));
     6251                        this->inputs->AddInput(new PentaVertexInput(VyMeshEnum,nodeinputs));
     6252                        this->inputs->AddInput(new PentaVertexInput(VzMeshEnum,nodeinputs));
     6253                        break;
     6254
    62406255                default:
    62416256                        /*No update for other solution types*/
Note: See TracChangeset for help on using the changeset viewer.