Changeset 5946


Ignore:
Timestamp:
09/22/10 10:40:57 (14 years ago)
Author:
seroussi
Message:

cleaned tria and penta

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

Legend:

Unmodified
Added
Removed

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

    r5943 r5946  
    19791979ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynViscous(void){
    19801980
    1981         const int    NUMVERTICESm=3;  //MacAyealNUMVERTICES
    1982         const int    numdofm=2*NUMVERTICESm;
    1983         const int    NUMVERTICESp=6; //Pattyn NUMVERTICES
    1984         const int    numdofp=2*NUMVERTICESp;
     1981        /*Constants*/
     1982        const int    numdofm=NDOF2*NUMVERTICES2D;
     1983        const int    numdofp=NDOF2*NUMVERTICES;
    19851984        const int    numdoftotal=2*NDOF2*NUMVERTICES;
    1986         int             i,j,ig;
    1987         double       xyz_list[NUMVERTICESp][3];
     1985
     1986        /*Intermediaries */
     1987        int         i,j,ig;
     1988        double      Jdet;
     1989        double      viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
     1990        double      epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     1991        double      xyz_list[NUMVERTICES][3];
     1992        double      B[3][numdofp];
     1993        double      Bprime[3][numdofm];
     1994        double      D[3][3]={0.0};            // material matrix, simple scalar matrix.
     1995        double      D_scalar;
     1996        double      Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix
     1997        double      Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point.
    19881998        GaussPenta *gauss=NULL;
    19891999        GaussTria  *gauss_tria=NULL;
    1990         double viscosity; //viscosity
    1991         double oldviscosity; //viscosity
    1992         double newviscosity; //viscosity
    1993         double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    1994         double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    1995         double B[3][numdofp];
    1996         double Bprime[3][numdofm];
    1997         double L[2][numdofp];
    1998         double D[3][3]={0.0};            // material matrix, simple scalar matrix.
    1999         double D_scalar;
    2000         double DL[2][2]={0.0}; //for basal drag
    2001         double DL_scalar;
    2002         double Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix
    2003         double Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point.
    2004         double Jdet;
    2005         double  alpha2_list[3];
    2006         double  alpha2;
    2007         double viscosity_overshoot;
    20082000
    20092001        /*Find penta on bed as pattyn must be coupled to the dofs on the bed: */
    … …  
    20192011
    20202012        /* Get node coordinates and dof list: */
    2021         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESp);
     2013        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    20222014        this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
    20232015        Input* vx_input=inputs->GetInput(VxEnum);       ISSMASSERT(vx_input);
    … …  
    20342026                gauss->SynchronizeGaussTria(gauss_tria);
    20352027
     2028                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     2029                GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss);
     2030                tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
     2031
    20362032                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    20372033                this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    … …  
    20392035                matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    20402036
    2041                 GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss);
    2042                 tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    2043 
    2044                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    20452037                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    20462038                D_scalar=2*newviscosity*gauss->weight*Jdet;
    … …  
    21392131ElementMatrix* Penta::CreateKMatrixDiagnosticHutter(void){
    21402132
     2133        /*Constants*/
     2134        const int numdof=NDOF2*NUMVERTICES;
     2135
    21412136        /*Intermediaries*/
    2142         const int numdof=NDOF2*NUMVERTICES;
    21432137        int       connectivity[2];
     2138        int       i,i0,i1,j0,j1;
    21442139        double    one0,one1;
    2145         int       i,i0,i1,j0,j1;
    21462140
    21472141        /*Initialize Element matrix and return if necessary*/
    … …  
    22382232ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyeal3dViscous(void){
    22392233
     2234        /*Constants*/
    22402235        const int    numdof2d=2*NUMVERTICES2D;
    2241         int             i,j,ig;
    2242         double       xyz_list[NUMVERTICES][3];
     2236
     2237        /*Intermediaries */
     2238        int         i,j,ig;
     2239        double      Jdet;
     2240        double      viscosity, oldviscosity, newviscosity, viscosity_overshoot;
     2241        double      epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     2242        double      xyz_list[NUMVERTICES][3];
     2243        double      B[3][numdof2d];
     2244        double      Bprime[3][numdof2d];
     2245        double      D[3][3]={0.0};            // material matrix, simple scalar matrix.
     2246        double      D_scalar;
     2247        double      Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point.
     2248        Tria*       tria=NULL;
     2249        Penta*      pentabase=NULL;
    22432250        GaussPenta *gauss=NULL;
    22442251        GaussTria  *gauss_tria=NULL;
    2245         double viscosity, oldviscosity, newviscosity, viscosity_overshoot;
    2246         double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    2247         double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    2248         double B[3][numdof2d];
    2249         double Bprime[3][numdof2d];
    2250         double L[2][numdof2d];
    2251         double D[3][3]={0.0};            // material matrix, simple scalar matrix.
    2252         double D_scalar;
    2253         double DL[2][2]={0.0}; //for basal drag
    2254         double DL_scalar;
    2255         double Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point.
    2256         double Jdet;
    2257         double  slope[2]={0.0};
    2258         double  slope_magnitude;
    2259         double  alpha2_list[3];
    2260         double  alpha2;
    2261         double MAXSLOPE=.06; // 6 %
    2262         double MOUNTAINKEXPONENT=10;
    2263         Tria*  tria=NULL;
    2264         Penta* pentabase=NULL;
    22652252
    22662253        /*Find penta on bed as this is a macayeal elements: */
    … …  
    22882275                gauss->SynchronizeGaussTria(gauss_tria);
    22892276
     2277                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    22902278                tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria);
    22912279                tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    2292 
    22932280
    22942281                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    … …  
    22962283                matice->GetViscosity3d(&viscosity, &epsilon[0]);
    22972284                matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
     2285
    22982286                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    2299                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    23002287                D_scalar=2*newviscosity*gauss->weight*Jdet;
    23012288                for (i=0;i<3;i++) D[i][i]=D_scalar;
    … …  
    25552542        double     LprimeStokes[14][numdof2d];
    25562543        double     DLStokes[14][14]={0.0};
    2557         double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble
    25582544        double     Ke_drag_gaussian[numdof2d][numdof2d];
    25592545        Friction*  friction=NULL;
    … …  
    26132599                                        &Ke_drag_gaussian[0][0],0);
    26142600
    2615                 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_temp[i][j]+=Ke_drag_gaussian[i][j];
    2616         }
    2617 
    2618         for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof+j]+=Ke_temp[i][j];
     2601                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof+j]+=Ke_drag_gaussian[i][j];
     2602        }
    26192603
    26202604        /*Clean up and return*/
    … …  
    26412625ElementMatrix* Penta::CreateKMatrixDiagnosticVertVolume(void){
    26422626
    2643         /* local declarations */
     2627        /*Constants*/
    26442628        const int    numdof=NDOF1*NUMVERTICES;
    2645         int             i,j,ig;
    2646         double       xyz_list[NUMVERTICES][3];
    2647         GaussPenta *gauss=NULL;
    2648         double  Ke_gg[numdof][numdof]={0.0};
    2649         double  Jdet;
    2650         double  B[NDOF1][NUMVERTICES];
    2651         double  Bprime[NDOF1][NUMVERTICES];
    2652         double  DL_scalar;
     2629
     2630        /*Intermediaries */
     2631        int         i,j,ig;
     2632        double      Jdet;
     2633        double      xyz_list[NUMVERTICES][3];
     2634        double      B[NDOF1][NUMVERTICES];
     2635        double      Bprime[NDOF1][NUMVERTICES];
     2636        double      DL_scalar;
     2637        double      Ke_gg[numdof][numdof]={0.0};
     2638        GaussPenta  *gauss=NULL;
    26532639
    26542640        /*Initialize Element matrix and return if necessary*/
    … …  
    26652651                gauss->GaussPoint(ig);
    26662652
     2653                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    26672654                GetBVert(&B[0][0], &xyz_list[0][0], gauss);
    26682655                GetBprimeVert(&Bprime[0][0], &xyz_list[0][0], gauss);
    26692656
    2670                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    26712657                DL_scalar=gauss->weight*Jdet;
    26722658
    … …  
    27612747ElementMatrix* Penta::CreateKMatrixThermalVolume(void){
    27622748
    2763         int i,j;
    2764         int     ig;
    2765         int found=0;
     2749        /*Constants*/
    27662750        const int    numdof=NDOF1*NUMVERTICES;
    2767         double       xyz_list[NUMVERTICES][3];
    2768         int*         doflist=NULL;
    2769         GaussPenta *gauss=NULL;
    2770         double  K[2][2]={0.0};
    2771         double  u,v,w;
    2772         double     Ke_gaussian_conduct[numdof][numdof];
    2773         double     Ke_gaussian_advec[numdof][numdof];
    2774         double     Ke_gaussian_artdiff[numdof][numdof];
    2775         double     Ke_gaussian_transient[numdof][numdof];
     2751
     2752        /*Intermediaries */
     2753        bool       artdiff;
     2754        int        i,j,ig,found=0;
     2755        double     Jdet,u,v,w,epsvel;
     2756        double     gravity,rho_ice,rho_water;
     2757        double     heatcapacity,thermalconductivity,dt;
     2758        double     xyz_list[NUMVERTICES][3];
    27762759        double     B[3][numdof];
    27772760        double     Bprime[3][numdof];
    … …  
    27832766        double     D_scalar;
    27842767        double     D[3][3];
    2785         double     l1l2l3[3];
    2786         double     tl1l2l3D[3];
     2768        double     K[2][2]={0.0};
    27872769        double     tBD[3][numdof];
    27882770        double     tBD_conduct[3][numdof];
    … …  
    27902772        double     tBD_artdiff[3][numdof];
    27912773        double     tLD[numdof];
    2792         double     Jdet;
    2793         double     gravity,rho_ice,rho_water;
    2794         double     heatcapacity,thermalconductivity;
    2795         double     mixed_layer_capacity,thermal_exchange_velocity;
    2796         double dt,epsvel;
    2797         bool   artdiff;
    2798         Tria*  tria=NULL;
     2774        double     Ke_gaussian_conduct[numdof][numdof];
     2775        double     Ke_gaussian_advec[numdof][numdof];
     2776        double     Ke_gaussian_artdiff[numdof][numdof];
     2777        double     Ke_gaussian_transient[numdof][numdof];
     2778        Tria*      tria=NULL;
     2779        GaussPenta *gauss=NULL;
    27992780
    28002781        /*Initialize Element matrix and return if necessary*/
    … …  
    30193000ElementVector* Penta::CreatePVectorCouplingPattynStokesViscous(void){
    30203001
    3021         const int numdof=NUMVERTICES*NDOF4;
    3022         int i,j;
    3023         int     ig;
    3024         double             xyz_list_tria[NUMVERTICES2D][3];
    3025         double         xyz_list[NUMVERTICES][3];
    3026         double             bed_normal[3];
     3002        /*Constants*/
     3003        const int   numdof=NUMVERTICES*NDOF4;
     3004
     3005        /*Intermediaries */
     3006        int         i,j,ig;
     3007        int         approximation;
     3008        double      viscosity,Jdet;
     3009        double      stokesreconditioning;
     3010        double      epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     3011        double      dw[3];
     3012        double      xyz_list[NUMVERTICES][3];
     3013        double      l1l6[6]; //for the six nodes of the penta
     3014        double      dh1dh6[3][6]; //for the six nodes of the penta
    30273015        GaussPenta *gauss=NULL;
    3028         double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    3029         double  viscosity, w, alpha2_gauss;
    3030         double  dw[3];
    3031         double     Pe_gaussian[24]={0.0}; //for the six nodes
    3032         double     l1l6[6]; //for the six nodes of the penta
    3033         double     dh1dh6[3][6]; //for the six nodes of the penta
    3034         double     Jdet;
    3035         double     Jdet2d;
    3036         Tria*     tria=NULL;
    3037         Friction* friction=NULL;
    3038         double stokesreconditioning;
    3039         int    approximation;
    3040         int    analysis_type;
    30413016
    30423017        /*Initialize Element vector and return if necessary*/
    30433018        if(IsOnWater()) return NULL;
     3019        inputs->GetParameterValue(&approximation,ApproximationEnum);
     3020        if(approximation!=PattynStokesApproximationEnum) return NULL;
     3021        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     3022
     3023        /*Retrieve all inputs and parameters*/
     3024        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3025        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
     3026        Input* vx_input=inputs->GetInput(VxEnum);               ISSMASSERT(vx_input);
     3027        Input* vy_input=inputs->GetInput(VyEnum);               ISSMASSERT(vy_input);
     3028        Input* vz_input=inputs->GetInput(VzEnum);               ISSMASSERT(vz_input);
     3029        Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   ISSMASSERT(vzpattyn_input);
     3030
     3031        /* Start  looping on the number of gaussian points: */
     3032        gauss=new GaussPenta(5,5);
     3033        for (ig=gauss->begin();ig<gauss->end();ig++){
     3034
     3035                gauss->GaussPoint(ig);
     3036
     3037                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     3038                GetNodalFunctionsP1(&l1l6[0], gauss);
     3039                GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     3040
     3041                vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
     3042
     3043                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     3044                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     3045
     3046                for(i=0;i<NUMVERTICES;i++){
     3047                        pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]/2;
     3048                        pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2;
     3049                        pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+dw[2]*dh1dh6[2][i])/2;
     3050                        pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i];
     3051                }
     3052        }
     3053
     3054        /*Clean up and return*/
     3055        delete gauss;
     3056        return pe;
     3057}
     3058/*}}}*/
     3059/*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{1*/
     3060ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){
     3061
     3062        /*Constants*/
     3063        const int numdof=NUMVERTICES*NDOF4;
     3064
     3065        /*Intermediaries*/
     3066        int         i,j,ig;
     3067        int         approximation,analysis_type;
     3068        double      Jdet,Jdet2d;
     3069        double      stokesreconditioning;
     3070        double     bed_normal[3];
     3071        double      epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     3072        double      viscosity, w, alpha2_gauss;
     3073        double      dw[3];
     3074        double     xyz_list_tria[NUMVERTICES2D][3];
     3075        double      xyz_list[NUMVERTICES][3];
     3076        double      l1l6[6]; //for the six nodes of the penta
     3077        Tria*       tria=NULL;
     3078        Friction*   friction=NULL;
     3079        GaussPenta  *gauss=NULL;
     3080
     3081        /*Initialize Element vector and return if necessary*/
     3082        if(IsOnWater() || !IsOnBed() || IsOnShelf()) return NULL;
    30443083        inputs->GetParameterValue(&approximation,ApproximationEnum);
    30453084        if(approximation!=PattynStokesApproximationEnum) return NULL;
    … …  
    30553094        Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   ISSMASSERT(vzpattyn_input);
    30563095
    3057         /* Start  looping on the number of gaussian points: */
    3058         gauss=new GaussPenta(5,5);
    3059         for (ig=gauss->begin();ig<gauss->end();ig++){
    3060 
    3061                 gauss->GaussPoint(ig);
    3062 
    3063                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    3064                 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    3065 
    3066                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3067 
    3068                 GetNodalFunctionsP1(&l1l6[0], gauss);
    3069                 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
    3070 
    3071                 vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    3072 
    3073                 for(i=0;i<NUMVERTICES;i++){
    3074                         pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]/2;
    3075                         pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2;
    3076                         pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+dw[2]*dh1dh6[2][i])/2;
    3077                         pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i];
    3078                 }
    3079         }
    3080 
    3081         /*Clean up and return*/
    3082         delete gauss;
    3083         return pe;
    3084 }
    3085 /*}}}*/
    3086 /*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{1*/
    3087 ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){
    3088 
    3089         const int numdof=NUMVERTICES*NDOF4;
    3090         int i,j;
    3091         int     ig;
    3092         double             xyz_list_tria[NUMVERTICES2D][3];
    3093         double         xyz_list[NUMVERTICES][3];
    3094         double             bed_normal[3];
    3095         GaussPenta *gauss=NULL;
    3096         double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    3097         double  viscosity, w, alpha2_gauss;
    3098         double  dw[3];
    3099         double     Pe_gaussian[24]={0.0}; //for the six nodes
    3100         double     l1l6[6]; //for the six nodes of the penta
    3101         double     dh1dh6[3][6]; //for the six nodes of the penta
    3102         double     Jdet;
    3103         double     Jdet2d;
    3104         Tria*     tria=NULL;
    3105         Friction* friction=NULL;
    3106         double stokesreconditioning;
    3107         int    approximation;
    3108         int    analysis_type;
    3109 
    3110         /*Initialize Element vector and return if necessary*/
    3111         if(IsOnWater() || !IsOnBed() || IsOnShelf()) return NULL;
    3112         inputs->GetParameterValue(&approximation,ApproximationEnum);
    3113         if(approximation!=PattynStokesApproximationEnum) return NULL;
    3114         ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
    3115 
    3116         /*Retrieve all inputs and parameters*/
    3117         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3118         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3119         this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    3120         Input* vx_input=inputs->GetInput(VxEnum);               ISSMASSERT(vx_input);
    3121         Input* vy_input=inputs->GetInput(VyEnum);               ISSMASSERT(vy_input);
    3122         Input* vz_input=inputs->GetInput(VzEnum);               ISSMASSERT(vz_input);
    3123         Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   ISSMASSERT(vzpattyn_input);
    3124 
    3125 
    3126         for(i=0;i<NUMVERTICES2D;i++){
    3127                 for(j=0;j<3;j++){
    3128                         xyz_list_tria[i][j]=xyz_list[i][j];
    3129                 }
    3130         }
     3096        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    31313097
    31323098        /*build friction object, used later on: */
    … …  
    31393105                gauss->GaussPoint(ig);
    31403106
    3141                 /*Get the Jacobian determinant */
    31423107                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3143 
    3144                 /*Get L matrix : */
    31453108                GetNodalFunctionsP1(l1l6, gauss);
    31463109
    3147                 /*Get normal vecyor to the bed */
     3110                vzpattyn_input->GetParameterValue(&w, gauss);
     3111                vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
     3112
    31483113                BedNormal(&bed_normal[0],xyz_list_tria);
    3149 
    3150                 /*Get Viscosity*/
    31513114                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    31523115                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    3153 
    3154                 /*Get friction*/
    31553116                friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
    3156 
    3157                 /*Get w and its derivatives*/
    3158                 vzpattyn_input->GetParameterValue(&w, gauss);
    3159                 vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    31603117
    31613118                for(i=0;i<NUMVERTICES2D;i++){
    … …  
    32303187ElementVector* Penta::CreatePVectorDiagnosticHutter(void){
    32313188
    3232         int i,j,k;
    3233         int      ig;
     3189        /*Constants*/
    32343190        const int numdofs=NDOF2*NUMVERTICES;
    3235         double    xyz_list[NUMVERTICES][3];
    3236         double    xyz_list_segment[2][3];
    3237         double    z_list[NUMVERTICES];
    3238         double    z_segment[2];
    3239         double    slope2,constant_part;
    3240         int       node0,node1;
    3241         GaussPenta* gauss=NULL;
    3242         double   slope[2];
    3243         double   z_g;
    3244         double   rho_ice,gravity,n,B;
    3245         double   Jdet;
    3246         double   ub,vb;
    3247         double   surface,thickness;
    3248         int  connectivity[2];
     3191
     3192        /*Intermediaries*/
     3193        int          i,j,k,ig;
     3194        int          node0,node1;
     3195        int          connectivity[2];
     3196        double       Jdet;
     3197        double       xyz_list[NUMVERTICES][3];
     3198        double       xyz_list_segment[2][3];
     3199        double       z_list[NUMVERTICES];
     3200        double       z_segment[2],slope[2];
     3201        double       slope2,constant_part;
     3202        double       rho_ice,gravity,n,B;
     3203        double       ub,vb,z_g,surface,thickness;
     3204        GaussPenta*  gauss=NULL;
    32493205
    32503206        /*Initialize Element vector and return if necessary*/
    … …  
    32663222        /*Loop on the three segments*/
    32673223        for(i=0;i<3;i++){
    3268 
    32693224                node0=i;
    32703225                node1=i+3;
    … …  
    33353290ElementVector* Penta::CreatePVectorDiagnosticPattyn(void){
    33363291
    3337         /* node data: */
     3292        /*Constants*/
    33383293        const int    numdof=NDOF2*NUMVERTICES;
    3339         int i,j;
    3340         int     ig;
    3341         double       xyz_list[NUMVERTICES][3];
    3342         double  slope[3]; //do not put 2! this goes into GetParameterDerivativeValue, which addresses slope[3] also!
    3343         double  driving_stress_baseline;
    3344         double  thickness;
    3345         GaussPenta *gauss=NULL;
    3346         double Jdet;
    3347         double l1l6[6];
    3348         double  pe_g_gaussian[numdof];
     3294
     3295        /*Intermediaries*/
     3296        int         i,j,ig;
     3297        double      Jdet;
     3298        double      slope[3]; //do not put 2! this goes into GetParameterDerivativeValue, which addresses slope[3] also!
     3299        double      driving_stress_baseline,thickness;
     3300        double      xyz_list[NUMVERTICES][3];
     3301        double      l1l6[6];
     3302        GaussPenta  *gauss=NULL;
    33493303
    33503304        /*Initialize Element vector and return if necessary*/
    … …  
    33633317                gauss->GaussPoint(ig);
    33643318
     3319                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     3320                GetNodalFunctionsP1(l1l6, gauss);
     3321
    33653322                thickness_input->GetParameterValue(&thickness, gauss);
    33663323                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    3367 
    3368                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3369                 GetNodalFunctionsP1(l1l6, gauss);
    33703324
    33713325                driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
    … …  
    33963350ElementVector* Penta::CreatePVectorDiagnosticStokesViscous(void){
    33973351
    3398         int i,j;
    3399         int     ig;
    3400         const int numdof=NUMVERTICES*NDOF4;
    3401         int       numdof2d=NUMVERTICES2D*NDOF4;
    3402         double         gravity,rho_ice,rho_water;
    3403         double             xyz_list_tria[NUMVERTICES2D][3];
    3404         double         xyz_list[NUMVERTICES][3];
    3405         double             bed_normal[3];
    3406         double         bed;
    3407         GaussPenta *gauss=NULL;
    3408         double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    3409         double  viscosity;
    3410         double  water_pressure;
    3411         double     Pe_temp[27]={0.0}; //for the six nodes and the bubble
    3412         double     Pe_gaussian[27]={0.0}; //for the six nodes and the bubble
    3413         double     Ke_temp[27][3]={0.0}; //for the six nodes and the bubble
    3414         double     Pe_reduced[numdof]; //for the six nodes only
    3415         double     Ke_gaussian[27][3];
    3416         double     l1l6[6]; //for the six nodes of the penta
     3352        /*Constants*/
     3353        const int numdofbubble=NDOF4*NUMVERTICES+NDOF3*1;
     3354
     3355        /*Intermediaries*/
     3356        int        i,j,ig;
     3357        int        approximation;
     3358        double     Jdet,viscosity;
     3359        double     gravity,rho_ice,stokesreconditioning;
     3360        double     xyz_list[NUMVERTICES][3];
     3361        double     epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    34173362        double     l1l7[7]; //for the six nodes and the bubble
    3418         double     B[8][27];
    3419         double     B_prime[8][27];
     3363        double     B[8][numdofbubble];
     3364        double     B_prime[8][numdofbubble];
    34203365        double     B_prime_bubble[8][3];
    3421         double     Jdet;
    3422         double     Jdet2d;
    34233366        double     D[8][8]={0.0};
    34243367        double     D_scalar;
    3425         double     tBD[27][8];
    3426         Tria*            tria=NULL;
    3427         double stokesreconditioning;
    3428         int  approximation;
     3368        double     tBD[numdofbubble][8];
     3369        double     Pe_gaussian[numdofbubble]={0.0}; //for the six nodes and the bubble
     3370        double     Ke_temp[numdofbubble][3]={0.0}; //for the six nodes and the bubble
     3371        double     Ke_gaussian[numdofbubble][3];
     3372        GaussPenta *gauss=NULL;
    34293373
    34303374        /*Initialize Element vector and return if necessary*/
    … …  
    34363380        /*Retrieve all inputs and parameters*/
    34373381        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    3438         rho_water=matpar->GetRhoWater();
    34393382        rho_ice=matpar->GetRhoIce();
    34403383        gravity=matpar->GetG();
    … …  
    34433386        Input* vy_input=inputs->GetInput(VyEnum);   ISSMASSERT(vy_input);
    34443387        Input* vz_input=inputs->GetInput(VzEnum);   ISSMASSERT(vz_input);
    3445         Input* bed_input=inputs->GetInput(BedEnum); ISSMASSERT(bed_input);
    34463388
    34473389        /* Start  looping on the number of gaussian points: */
    … …  
    34513393                gauss->GaussPoint(ig);
    34523394
     3395                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     3396                GetBStokes(&B[0][0],&xyz_list[0][0],gauss);
     3397                GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss);
     3398                GetNodalFunctionsMINI(&l1l7[0], gauss);
     3399
    34533400                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    34543401                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    3455                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3456                 GetNodalFunctionsMINI(&l1l7[0], gauss);
    3457 
    3458                 /* Build gaussian vector */
     3402
    34593403                for(i=0;i<NUMVERTICES+1;i++){
    3460                         Pe_gaussian[i*NDOF4+2]=-rho_ice*gravity*Jdet*gauss->weight*l1l7[i];
    3461                 }
    3462 
    3463                 /*Add Pe_gaussian to terms in Pe_temp. Watch out for column orientation from matlab: */
    3464                 for(i=0;i<27;i++) Pe_temp[i]+=Pe_gaussian[i];
    3465 
    3466                 /*Get B and Bprime matrices: */
    3467                 GetBStokes(&B[0][0],&xyz_list[0][0],gauss);
    3468                 GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss);
     3404                        Pe_gaussian[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l7[i];
     3405                }
    34693406
    34703407                /*Get bubble part of Bprime */
    34713408                for(i=0;i<8;i++) for(j=0;j<3;j++) B_prime_bubble[i][j]=B_prime[i][j+24];
    34723409
    3473                 /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    3474                  * onto this scalar matrix, so that we win some computational time: */
    34753410                D_scalar=gauss->weight*Jdet;
    34763411                for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
    34773412                for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
    34783413
    3479                 /*  Do the triple product tB*D*Bprime: */
    3480                 MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0);
    3481                 MatrixMultiply(&tBD[0][0],27,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0);
    3482 
    3483                 /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */
    3484                 for(i=0;i<27;i++) for(j=0;j<3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
     3414                MatrixMultiply(&B[0][0],8,numdofbubble,1,&D[0][0],8,8,0,&tBD[0][0],0);
     3415                MatrixMultiply(&tBD[0][0],numdofbubble,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0);
     3416
     3417                for(i=0;i<numdofbubble;i++) for(j=0;j<NDOF3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
    34853418        }
    34863419
    34873420        /*Condensation*/
    3488         ReduceVectorStokes(pe->values, &Ke_temp[0][0], &Pe_temp[0]);
     3421        ReduceVectorStokes(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]);
    34893422
    34903423        /*Clean up and return*/
    … …  
    34963429ElementVector* Penta::CreatePVectorDiagnosticStokesShelf(void){
    34973430
    3498         int i,j;
    3499         int     ig;
    3500         const int numdof=NUMVERTICES*NDOF4;
    3501         int       numdof2d=NUMVERTICES2D*NDOF4;
    3502         double         gravity,rho_ice,rho_water;
    3503         double             xyz_list_tria[NUMVERTICES2D][3];
    3504         double         xyz_list[NUMVERTICES][3];
    3505         double             bed_normal[3];
    3506         double         bed;
    3507         GaussPenta *gauss=NULL;
    3508         double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    3509         double  viscosity;
    3510         double  water_pressure;
    3511         double     Pe_temp[27]={0.0}; //for the six nodes and the bubble
    3512         double     Pe_gaussian[27]={0.0}; //for the six nodes and the bubble
    3513         double     Ke_temp[27][3]={0.0}; //for the six nodes and the bubble
    3514         double     Pe_reduced[numdof]; //for the six nodes only
    3515         double     Ke_gaussian[27][3];
    3516         double     l1l6[6]; //for the six nodes of the penta
    3517         double     l1l7[7]; //for the six nodes and the bubble
    3518         double     B[8][27];
    3519         double     B_prime[8][27];
    3520         double     B_prime_bubble[8][3];
    3521         double     Jdet;
    3522         double     Jdet2d;
    3523         double     D[8][8]={0.0};
    3524         double     D_scalar;
    3525         double     tBD[27][8];
    3526         Tria*            tria=NULL;
    3527         double stokesreconditioning;
    3528         int  approximation;
     3431        /*Intermediaries*/
     3432        int         i,j,ig;
     3433        int         approximation;
     3434        double      gravity,rho_water,bed,water_pressure;
     3435        double          xyz_list_tria[NUMVERTICES2D][3];
     3436        double      xyz_list[NUMVERTICES][3];
     3437        double          bed_normal[3];
     3438        double      l1l6[6]; //for the six nodes of the penta
     3439        double      Jdet2d;
     3440        Tria*       tria=NULL;
     3441        GaussPenta  *gauss=NULL;
    35293442
    35303443        /*Initialize Element vector and return if necessary*/
    … …  
    35353448
    35363449        /*Retrieve all inputs and parameters*/
    3537         this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    35383450        rho_water=matpar->GetRhoWater();
    3539         rho_ice=matpar->GetRhoIce();
    35403451        gravity=matpar->GetG();
    35413452        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3542         Input* vx_input=inputs->GetInput(VxEnum);   ISSMASSERT(vx_input);
    3543         Input* vy_input=inputs->GetInput(VyEnum);   ISSMASSERT(vy_input);
    3544         Input* vz_input=inputs->GetInput(VzEnum);   ISSMASSERT(vz_input);
    35453453        Input* bed_input=inputs->GetInput(BedEnum); ISSMASSERT(bed_input);
    35463454
    3547 
    3548         for(i=0;i<NUMVERTICES2D;i++){
    3549                 for(j=0;j<3;j++){
    3550                         xyz_list_tria[i][j]=xyz_list[i][j];
    3551                 }
    3552         }
     3455        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    35533456
    35543457        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    … …  
    35593462
    35603463                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     3464                GetNodalFunctionsP1(l1l6, gauss);
     3465
    35613466                bed_input->GetParameterValue(&bed, gauss);
    3562                 GetNodalFunctionsP1(l1l6, gauss);
    3563 
     3467                BedNormal(&bed_normal[0],xyz_list_tria);
    35643468                water_pressure=gravity*rho_water*bed;
    3565 
    3566                 BedNormal(&bed_normal[0],xyz_list_tria);
    35673469
    35683470                for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) pe->values[i*NDOF4+j]+=water_pressure*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j];
    … …  
    36053507ElementVector* Penta::CreatePVectorDiagnosticVertVolume(void){
    36063508
    3607         const int    numdof=NDOF1*NUMVERTICES;
    3608         double       xyz_list[NUMVERTICES][3];
    3609         int*         doflist=NULL;
    3610         int i;
    3611         int     ig;
     3509        /*Constants*/
     3510        const int  numdof=NDOF1*NUMVERTICES;
     3511
     3512        /*Intermediaries*/
     3513        int        i,ig;
     3514        int        approximation;
     3515        double     Jdet;
     3516        double     xyz_list[NUMVERTICES][3];
     3517        double     dudx,dvdy,dwdz;
     3518        double     du[3],dv[3],dw[3];
     3519        double     l1l6[6];
    36123520        GaussPenta *gauss=NULL;
    3613         double Jdet;
    3614         double l1l6[6];
    3615         Tria* tria=NULL;
    3616         double du[3];
    3617         double dv[3];
    3618         double dw[3];
    3619         double dudx,dvdy,dwdz;
    3620         int approximation;
    36213521
    36223522        /*Initialize Element vector and return if necessary*/
    … …  
    36403540                gauss->GaussPoint(ig);
    36413541
     3542                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     3543                GetNodalFunctionsP1(l1l6, gauss);
     3544
    36423545                vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss);
    36433546                vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss);
    … …  
    36503553                dvdy=dv[1];
    36513554
    3652                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3653                 GetNodalFunctionsP1(l1l6, gauss);
    3654 
    36553555                for (i=0;i<numdof;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*l1l6[i];
    36563556        }
    … …  
    37353635ElementVector* Penta::CreatePVectorThermalVolume(void){
    37363636
     3637        /*Constants*/
    37373638        const int  numdof=NUMVERTICES*NDOF1;
    3738         int i,j;
    3739         int     ig;
    3740         int found=0;
    3741         double        xyz_list[NUMVERTICES][3];
     3639
     3640        /*Intermediaries*/
     3641        int    i,j,ig,found=0;
     3642        int    friction_type;
     3643        double Jdet,phi,dt;
     3644        double rho_ice,heatcapacity;
     3645        double viscosity,temperature;
     3646        double scalar_def,scalar_transient;
     3647        double temperature_list[NUMVERTICES];
     3648        double xyz_list[NUMVERTICES][3];
     3649        double L[numdof];
     3650        double epsilon[6];
    37423651        GaussPenta *gauss=NULL;
    3743         double temperature_list[NUMVERTICES];
    3744         double temperature;
    3745         double gravity,rho_ice,rho_water;
    3746         double mixed_layer_capacity,heatcapacity;
    3747         double beta,meltingpoint,thermal_exchange_velocity;
    3748         int    friction_type;
    3749         double P_terms[numdof]={0.0};
    3750         double L[numdof];
    3751         double l1l2l3[3];
    3752         double basalfriction;
    3753         double epsilon[6];
    3754         double epsilon_sqr[3][3];
    3755         double epsilon_matrix[3][3];
    3756         double Jdet;
    3757         double viscosity;
    3758         double epsilon_eff;
    3759         double phi;
    3760         double t_pmp;
    3761         double scalar;
    3762         double scalar_def;
    3763         double scalar_ocean;
    3764         double scalar_transient;
    3765         Tria*  tria=NULL;
    3766         double dt;
    37673652
    37683653        /*Initialize Element vector and return if necessary*/
    … …  
    37733658        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    37743659        this->parameters->FindParam(&dt,DtEnum);
    3775         rho_water=matpar->GetRhoWater();
    37763660        rho_ice=matpar->GetRhoIce();
    3777         gravity=matpar->GetG();
    37783661        heatcapacity=matpar->GetHeatCapacity();
    3779         beta=matpar->GetBeta();
    3780         meltingpoint=matpar->GetMeltingPoint();
    37813662        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
    37823663        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
    37833664        Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input);
    37843665        Input* temperature_input=NULL;
    3785         if (dt){
    3786                 temperature_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(inputs);
    3787         }
     3666        if (dt) temperature_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(inputs);
    37883667
    37893668        /* Start  looping on the number of gaussian points: */
    … …  
    37933672                gauss->GaussPoint(ig);
    37943673
     3674                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     3675                GetNodalFunctionsP1(&L[0], gauss);
     3676
    37953677                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    37963678                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    3797                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3798                 GetNodalFunctionsP1(&L[0], gauss);
    37993679                GetPhi(&phi, &epsilon[0], viscosity);
    38003680

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

    r5934 r5946  
    30383038        /*Intermediaries */
    30393039        int       i,j,ig;
    3040         int*      doflist=NULL;
    30413040        double    xyz_list[NUMVERTICES][3];
    30423041        double    surface_normal[3];
    … …  
    32743273        /*Intermediaries */
    32753274        int        i,j,ig,dim;
    3276         int*       doflist=NULL;
    32773275        double     xyz_list[NUMVERTICES][3];
    32783276        double     Jdettria,dt,vx,vy;
    … …  
    34683466        /*Intermediaries */
    34693467        int        i,j,ig;
    3470         int*       doflist=NULL;
    34713468        double     xyz_list[NUMVERTICES][3];
    34723469        double     dhdt_g,melting_g,accumulation_g,Jdettria;
    … …  
    35133510        /*Intermediaries */
    35143511        int        i,j,ig;
    3515         int*       doflist=NULL;
    35163512        double     xyz_list[NUMVERTICES][3];
    35173513        double     melting_g,accumulation_g,dhdt_g,Jdettria;
    … …  
    35583554        /*Intermediaries */
    35593555        int        i,j,ig;
    3560         int*       doflist=NULL;
    35613556        double     xyz_list[NUMVERTICES][3];
    35623557        double     Jdettria,accumulation_g,melting_g;
    … …  
    37033698                        for (i=0;i<NUMVERTICES;i++){
    37043699                                for (j=0;j<NDOF2;j++){
    3705                                         pe_g_gaussian[i*NDOF2+j]=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss->weight*l1l2l3[i];
     3700                                        pe->values[i*NDOF2+j]+=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss->weight*l1l2l3[i];
    37063701                                }
    37073702                        }
    … …  
    37103705                        for (i=0;i<NUMVERTICES;i++){
    37113706                                for (j=0;j<NDOF2;j++){
    3712                                         pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l2l3[i];
     3707                                        pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l2l3[i];
    37133708                                }
    37143709                        }
    37153710                }
    3716                 for(i=0;i<numdof;i++) pe->values[i]+=pe_g_gaussian[i];
    37173711        }
    37183712
    … …  
    37303724        /*Intermediaries */
    37313725        int         i,ig;
    3732         int*        doflist=NULL;
    37333726        double      Jdet;
    37343727        double      thickness,thicknessobs,weight;
    37353728        double      xyz_list[NUMVERTICES][3];
    37363729        double      l1l2l3[3];
    3737         double      pe_g_gaussian[numdof];
    3738         GaussTria* gauss=NULL;
     3730        GaussTria*  gauss=NULL;
    37393731
    37403732        /*Initialize Element vector and return if necessary*/
    … …  
    37743766        /*Constants*/
    37753767        const int    numdof=NDOF2*NUMVERTICES;
    3776         double       xyz_list[NUMVERTICES][3];
    3777         double vx_list[NUMVERTICES];
    3778         double vy_list[NUMVERTICES];
    3779         double obs_vx_list[NUMVERTICES];
    3780         double obs_vy_list[NUMVERTICES];
    3781         double dux_list[NUMVERTICES];
    3782         double duy_list[NUMVERTICES];
    3783         double weights_list[NUMVERTICES];
    3784         int i;
    3785         int     ig;
     3768
     3769        /*Intermediaries */
     3770        int        i,ig,response;
     3771        double     Jdet;
     3772        double     obs_velocity_mag,velocity_mag;
     3773        double     dux,duy,meanvel,epsvel;
     3774        double     scalex=0,scaley=0,scale=0,S=0;
     3775        double     xyz_list[NUMVERTICES][3];
     3776        double     vx_list[NUMVERTICES];
     3777        double     vy_list[NUMVERTICES];
     3778        double     obs_vx_list[NUMVERTICES];
     3779        double     obs_vy_list[NUMVERTICES];
     3780        double     dux_list[NUMVERTICES];
     3781        double     duy_list[NUMVERTICES];
     3782        double     weights_list[NUMVERTICES];
     3783        double     l1l2l3[3];
    37863784        GaussTria* gauss=NULL;
    3787         double  obs_velocity_mag,velocity_mag;
    3788         double  dux,duy;
    3789         double  meanvel, epsvel;
    3790         double  pe_g_gaussian[numdof];
    3791         double Jdet;
    3792         double l1l2l3[3];
    3793         double scalex=0;
    3794         double scaley=0;
    3795         double scale=0;
    3796         double S=0;
    3797         int    response;
    37983785
    37993786        /*Initialize Element vector and return if necessary*/
    … …  
    39463933
    39473934                for (i=0;i<NUMVERTICES;i++){
    3948                         pe_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss->weight*l1l2l3[i];
    3949                         pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss->weight*l1l2l3[i];
     3935                        pe->values[i*NDOF2+0]+=dux*Jdet*gauss->weight*l1l2l3[i];
     3936                        pe->values[i*NDOF2+1]+=duy*Jdet*gauss->weight*l1l2l3[i];
    39503937                }
    3951 
    3952                 for(i=0;i<numdof;i++) pe->values[i]+=pe_g_gaussian[i];
    39533938        }
    39543939
    … …  
    39613946ElementVector* Tria::CreatePVectorAdjointStokes(void){
    39623947
    3963         const int    numdof=NDOF4*NUMVERTICES;
    3964         int i;
    3965         int     ig;
    3966         double       xyz_list[NUMVERTICES][3];
    3967         double vx_list[NUMVERTICES];
    3968         double vy_list[NUMVERTICES];
    3969         double obs_vx_list[NUMVERTICES];
    3970         double obs_vy_list[NUMVERTICES];
    3971         double dux_list[NUMVERTICES];
    3972         double duy_list[NUMVERTICES];
    3973         double weights_list[NUMVERTICES];
     3948        /*Intermediaries */
     3949        int        i,ig;
     3950        int        fit=-1;
     3951        int        response;
     3952        double     Jdet;
     3953        double     obs_velocity_mag,velocity_mag;
     3954        double     dux,duy,meanvel,epsvel;
     3955        double     scalex=0,scaley=0,scale=0,S=0;
     3956        double     xyz_list[NUMVERTICES][3];
     3957        double     vx_list[NUMVERTICES];
     3958        double     vy_list[NUMVERTICES];
     3959        double     obs_vx_list[NUMVERTICES];
     3960        double     obs_vy_list[NUMVERTICES];
     3961        double     dux_list[NUMVERTICES];
     3962        double     duy_list[NUMVERTICES];
     3963        double     weights_list[NUMVERTICES];
     3964        double     l1l2l3[3];
    39743965        GaussTria* gauss=NULL;
    3975         double  obs_velocity_mag,velocity_mag;
    3976         double  dux,duy;
    3977         double  meanvel, epsvel;
    3978         double  pe_g[numdof]={0.0};
    3979         double Jdet;
    3980         double l1l2l3[3];
    3981         double scalex=0;
    3982         double scaley=0;
    3983         double scale=0;
    3984         double S=0;
    3985         int    fit=-1;
    3986         int    response;
    39873966
    39883967        /*Initialize Element vector and return if necessary*/
    … …  
    41484127ElementVector* Tria::CreatePVectorDiagnosticHutter(void){
    41494128
    4150         /*Collapsed formulation: */
    4151         int       i;
    4152         const int numdofs=NDOF2*NUMVERTICES;
    4153         int*         doflist=NULL;
    4154         double    constant_part,ub,vb;
    4155         double    rho_ice,gravity,n,B;
    4156         double    slope[2];
    4157         double    thickness;
    4158         double    slope2;
    4159         int       connectivity;
     4129        /*Intermediaries */
     4130        int        i,connectivity;
     4131        double     constant_part,ub,vb;
     4132        double     rho_ice,gravity,n,B;
     4133        double     slope2,thickness;
     4134        double     slope[2];
    41604135        GaussTria* gauss=NULL;
    41614136
    … …  
    42164191ElementVector* Tria::CreatePVectorPrognostic_CG(void){
    42174192
    4218         /* local declarations */
    4219         int             i,j;
    4220 
    4221         /* node data: */
     4193        /*Constants*/
    42224194        const int    numdof=NDOF1*NUMVERTICES;
    4223         double       xyz_list[NUMVERTICES][3];
    4224         int*         doflist=NULL;
    4225         int          numberofdofspernode=1;
    4226         int     ig;
     4195
     4196        /*Intermediaries */
     4197        int        i,j,ig;
     4198        double     Jdettria,dt;
     4199        double     accumulation_g,melting_g,thickness_g;
     4200        double     xyz_list[NUMVERTICES][3];
     4201        double     L[NUMVERTICES];
    42274202        GaussTria* gauss=NULL;
    4228         double L[NUMVERTICES];
    4229         double Jdettria;
    4230         double  accumulation_g;
    4231         double  melting_g;
    4232         double  thickness_g;
    4233         double  dt;
    42344203
    42354204        /*Initialize Element vector and return if necessary*/
    … …  
    42514220
    42524221                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    4253                 GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode);
     4222                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    42544223
    42554224                accumulation_input->GetParameterValue(&accumulation_g,gauss);
    … …  
    42684237ElementVector* Tria::CreatePVectorPrognostic_DG(void){
    42694238
    4270         /* local declarations */
     4239        /*Constants*/
    42714240        const int    numdof=NDOF1*NUMVERTICES;
    4272         int             i,j;
    4273         int     ig;
    4274         double       xyz_list[NUMVERTICES][3];
    4275         int*         doflist=NULL;
    4276         int          numberofdofspernode=1;
     4241
     4242        /*Intermediaries */
     4243        int        i,j,ig;
     4244        double     Jdettria,dt;
     4245        double     accumulation_g,melting_g,thickness_g;
     4246        double     xyz_list[NUMVERTICES][3];
     4247        double     L[NUMVERTICES];
    42774248        GaussTria* gauss=NULL;
    4278         double L[NUMVERTICES];
    4279         double Jdettria;
    4280         double  accumulation_g;
    4281         double  melting_g;
    4282         double  thickness_g;
    4283         double  dt;
    42844249
    42854250        /*Initialize Element vector and return if necessary*/
    … …  
    43014266
    43024267                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    4303                 GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode);
     4268                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    43044269
    43054270                accumulation_input->GetParameterValue(&accumulation_g,gauss);
    … …  
    43184283ElementVector* Tria::CreatePVectorSlope(void){
    43194284
    4320         /* node data: */
     4285        /*Constants*/
    43214286        const int    numdof=NDOF1*NUMVERTICES;
    4322         int             i,j;
    4323         int     ig;
    4324         double       xyz_list[NUMVERTICES][3];
    4325         int*         doflist=NULL;
     4287       
     4288        /*Intermediaries */
     4289        int        i,j,ig;
     4290        int        analysis_type;
     4291        double     Jdet;
     4292        double     xyz_list[NUMVERTICES][3];
     4293        double     slope[2];
     4294        double     l1l2l3[3];
    43264295        GaussTria* gauss=NULL;
    4327         double Jdet;
    4328         double l1l2l3[3];
    4329         double  pe_g_gaussian[numdof];
    4330         double  slope[2];
    4331         int     analysis_type;
    43324296
    43334297        /*Initialize Element vector and return if necessary*/
    … …  
    43524316                gauss->GaussPoint(ig);
    43534317
     4318                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     4319                GetNodalFunctions(l1l2l3, gauss);
     4320
    43544321                slope_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    43554322
    4356                 /* Get Jacobian determinant: */
    4357                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    4358                
    4359                  /*Get nodal functions: */
    4360                 GetNodalFunctions(l1l2l3, gauss);
    4361 
    4362                 /*Build pe_g_gaussian vector: */
    43634323                if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){
    4364                         for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss->weight*slope[0]*l1l2l3[i];
     4324                        for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[0]*l1l2l3[i];
    43654325                }
    43664326                if ( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){
    4367                         for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss->weight*slope[1]*l1l2l3[i];
     4327                        for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[1]*l1l2l3[i];
    43684328                }
    4369 
    4370                 /*Add pe_g_gaussian vector to pe_g: */
    4371                 for(i=0;i<numdof;i++) pe->values[i]+=pe_g_gaussian[i];
    4372 
    43734329        }
    43744330
    … …  
    43814337ElementVector* Tria::CreatePVectorThermalShelf(void){
    43824338       
     4339        /*Constants*/
    43834340        const int  numdof=NUMVERTICES*NDOF1;
    4384         int i;
    4385         double       xyz_list[NUMVERTICES][3];
    4386         double mixed_layer_capacity;
    4387         double thermal_exchange_velocity;
    4388         double rho_water;
    4389         double rho_ice;
    4390         double heatcapacity;
    4391         double beta;
    4392         double meltingpoint;
    4393         double dt;
    4394         double pressure;
    4395         int     ig;
     4341
     4342        /*Intermediaries */
     4343        int        i,ig;
     4344        double     Jdet;
     4345        double     mixed_layer_capacity,thermal_exchange_velocity;
     4346        double     rho_ice,rho_water,pressure,dt,scalar_ocean;
     4347        double     meltingpoint,beta,heatcapacity,t_pmp;
     4348        double     xyz_list[NUMVERTICES][3];
     4349        double     l1l2l3[NUMVERTICES];
    43964350        GaussTria* gauss=NULL;
    4397         double  Jdet;
    4398         double  l1l2l3[NUMVERTICES];
    4399         double  t_pmp;
    4400         double  scalar_ocean;
    44014351
    44024352        /*Initialize Element vector and return if necessary*/
    … …  
    44294379
    44304380                scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
    4431                 if(dt){
    4432                         scalar_ocean=dt*scalar_ocean;
    4433                 }
     4381                if(dt) scalar_ocean=dt*scalar_ocean;
    44344382
    44354383                for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*l1l2l3[i];
    … …  
    44444392ElementVector* Tria::CreatePVectorThermalSheet(void){
    44454393
     4394        /*Constants*/
    44464395        const int  numdof=NUMVERTICES*NDOF1;
    4447         int i;
    4448         int     ig;
     4396
     4397        /*Intermediaries */
     4398        int        i,ig;
     4399        int        analysis_type,drag_type;
    44494400        double     xyz_list[NUMVERTICES][3];
    4450         double rho_ice;
    4451         double heatcapacity;
    4452         double dt;
    4453         double pressure_list[3];
    4454         double pressure;
    4455         int    drag_type;
    4456         double basalfriction;
    4457         Friction* friction=NULL;
    4458         double alpha2,vx,vy;
    4459         double geothermalflux_value;
     4401        double     Jdet,dt;
     4402        double     rho_ice,heatcapacity,geothermalflux_value;
     4403        double     basalfriction,alpha2,vx,vy,pressure;
     4404        double     pressure_list[3];
     4405        double     scalar;
     4406        double     l1l2l3[NUMVERTICES];
     4407        Friction*  friction=NULL;
    44604408        GaussTria* gauss=NULL;
    4461         double  Jdet;
    4462         double  P_terms[numdof]={0.0};
    4463         double  l1l2l3[NUMVERTICES];
    4464         double  scalar;
    4465         int analysis_type;
    44664409
    44674410        /*Initialize Element vector and return if necessary*/
    … …  
    44984441               
    44994442                scalar=gauss->weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
    4500                 if(dt){
    4501                         scalar=dt*scalar;
    4502                 }
     4443                if(dt) scalar=dt*scalar;
    45034444
    45044445                for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l2l3[i];
Note: See TracChangeset for help on using the changeset viewer.