Changeset 6212


Ignore:
Timestamp:
10/08/10 16:34:31 (14 years ago)
Author:
seroussi
Message:

added Franca2003 artificial diffusivity

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

Legend:

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

    r6200 r6212  
    28582858
    28592859        /*Intermediaries */
    2860         bool       artdiff;
     2860        int        artdiff;
    28612861        int        i,j,ig,found=0;
    28622862        double     Jdet,u,v,w,epsvel;
    28632863        double     gravity,rho_ice,rho_water;
    28642864        double     heatcapacity,thermalconductivity,dt;
     2865        double     scalar_artdiff;
     2866        double     tau_parameter,diameter,normu;
    28652867        double     xyz_list[NUMVERTICES][3];
    28662868        double     B[3][numdof];
     
    28712873        double     Bprime_advec[3][numdof];
    28722874        double     L[numdof];
     2875        double     dh1dh6[3][6];
    28732876        double     D_scalar;
    28742877        double     D[3][3];
     
    29032906        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
    29042907        Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input);
     2908        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
    29052909
    29062910        /* Start  looping on the number of gaussian points: */
     
    29612965                /*Artifficial diffusivity*/
    29622966
    2963                 if(artdiff){
     2967                if(artdiff==1){
    29642968                        /*Build K: */
    29652969                        D_scalar=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
     
    29722976                        MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0);
    29732977                        MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0);
     2978                }
     2979                else if(artdiff==2){
     2980
     2981                        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     2982
     2983                        normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
     2984                        if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
     2985                                tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
     2986                        }
     2987                        else tau_parameter=diameter/(2*normu);
     2988                        scalar_artdiff=tau_parameter*gauss->weight*Jdet;
     2989
     2990                        for(i=0;i<numdof;i++){
     2991                                for(j=0;j<numdof;j++){
     2992                                        Ke_gaussian_artdiff[i][j]=scalar_artdiff*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i])*(u*dh1dh6[0][j]+v*dh1dh6[1][j]+w*dh1dh6[2][j]);
     2993                                }
     2994                        }
    29742995                }
    29752996                else{
     
    37463767        /*Intermediaries*/
    37473768        int    i,j,ig,found=0;
    3748         int    friction_type;
     3769        int    friction_type,artdiff;
    37493770        double Jdet,phi,dt;
    37503771        double rho_ice,heatcapacity;
     3772        double thermalconductivity;
    37513773        double viscosity,temperature;
     3774        double tau_parameter,diameter,normu;
     3775        double u,v,w,scalar_artdiff;
    37523776        double scalar_def,scalar_transient;
    37533777        double temperature_list[NUMVERTICES];
    37543778        double xyz_list[NUMVERTICES][3];
    37553779        double L[numdof];
     3780        double dh1dh6[3][6];
    37563781        double epsilon[6];
    37573782        GaussPenta *gauss=NULL;
     
    37633788        /*Retrieve all inputs and parameters*/
    37643789        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3765         this->parameters->FindParam(&dt,DtEnum);
    37663790        rho_ice=matpar->GetRhoIce();
    37673791        heatcapacity=matpar->GetHeatCapacity();
     3792        thermalconductivity=matpar->GetThermalConductivity();
     3793        this->parameters->FindParam(&dt,DtEnum);
     3794        this->parameters->FindParam(&artdiff,ArtDiffEnum);
    37683795        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
    37693796        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
     
    37713798        Input* temperature_input=NULL;
    37723799        if (dt) temperature_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(inputs);
     3800        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
    37733801
    37743802        /* Start  looping on the number of gaussian points: */
     
    37953823                        scalar_transient=temperature*Jdet*gauss->weight;
    37963824                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
     3825                }
     3826
     3827                if(artdiff==2){
     3828                        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     3829
     3830                        vx_input->GetParameterValue(&u, gauss);
     3831                        vy_input->GetParameterValue(&v, gauss);
     3832                        vz_input->GetParameterValue(&w, gauss);
     3833
     3834                        normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
     3835                        if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
     3836                                tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
     3837                        }
     3838                        else tau_parameter=diameter/(2*normu);
     3839
     3840                        scalar_artdiff=tau_parameter*gauss->weight*Jdet*phi/(rho_ice*heatcapacity);
     3841                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_artdiff*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
    37973842                }
    37983843        }
     
    54415486}
    54425487/*}}}*/
     5488/*FUNCTION Penta::MinEdgeLength{{{1*/
     5489double Penta::MinEdgeLength(double xyz_list[6][3]){
     5490        /*Return the minimum lenght of the nine egdes of the penta*/
     5491
     5492        int    i,node0,node1;
     5493        int    edges[9][2]={{0,1},{0,2},{1,2},{3,4},{3,5},{4,5},{0,3},{1,4},{2,5}}; //list of the nine edges
     5494        double minlength=-1;
     5495        double length;
     5496       
     5497        for(i=0;i<9;i++){
     5498                /*Find the two nodes for this edge*/
     5499                node0=edges[i][0];
     5500                node1=edges[i][1];
     5501
     5502                /*Compute the length of this edge and compare it to the minimal length*/
     5503                length=pow(pow(xyz_list[node0][0]-xyz_list[node1][0],2.0)+pow(xyz_list[node0][1]-xyz_list[node1][1],2.0)+pow(xyz_list[node0][2]-xyz_list[node1][2],2.0),0.5);
     5504                if(length<minlength || minlength<0) minlength=length;
     5505        }
     5506
     5507        return minlength;
     5508}
     5509/*}}}*/
    54435510/*FUNCTION Penta::ReduceMatrixStokes {{{1*/
    54445511void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){
     
    56045671        *(surface_normal+1)=normal[1]/normal_norm;
    56055672        *(surface_normal+2)=normal[2]/normal_norm;
    5606 
    5607 }
    5608 /*}}}*/
     5673}
     5674/*}}}*/
  • TabularUnified issm/trunk/src/c/objects/Elements/Penta.h

    r6200 r6212  
    218218                bool    IsOnShelf(void);
    219219                bool    IsOnWater(void);
    220                 void      ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
     220                double  MinEdgeLength(double xyz_list[6][3]);
     221           void   ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
    221222                void      ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);
    222223                void      SetClone(int* minranks);
  • TabularUnified issm/trunk/src/c/objects/Elements/Tria.cpp

    r6200 r6212  
    24152415
    24162416        /*Intermediaries */
    2417         bool       artdiff;
     2417        int        artdiff;
    24182418        int        i,j,ig,dim;
    24192419        double     Jdettria ,vx,vy,dvxdx,dvydy;
     
    25862586
    25872587        /*Intermediaries */
    2588         bool    artdiff;
     2588        int     artdiff;
    25892589        int     i,j,ig,dim;
    25902590        double  nx,ny,norm,Jdettria;
     
    31573157
    31583158        /*Intermediaries */
    3159         bool       artdiff;
     3159        int        artdiff;
    31603160        int        i,j,ig,dim;
    31613161        double     Jdettria,DL_scalar,dt;
Note: See TracChangeset for help on using the changeset viewer.