Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5945) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5946) @@ -1979,31 +1979,23 @@ ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynViscous(void){ - const int NUMVERTICESm=3; //MacAyealNUMVERTICES - const int numdofm=2*NUMVERTICESm; - const int NUMVERTICESp=6; //Pattyn NUMVERTICES - const int numdofp=2*NUMVERTICESp; + /*Constants*/ + const int numdofm=NDOF2*NUMVERTICES2D; + const int numdofp=NDOF2*NUMVERTICES; const int numdoftotal=2*NDOF2*NUMVERTICES; - int i,j,ig; - double xyz_list[NUMVERTICESp][3]; + + /*Intermediaries */ + int i,j,ig; + double Jdet; + double viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity + double epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + double xyz_list[NUMVERTICES][3]; + double B[3][numdofp]; + double Bprime[3][numdofm]; + double D[3][3]={0.0}; // material matrix, simple scalar matrix. + double D_scalar; + double Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix + double Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point. GaussPenta *gauss=NULL; GaussTria *gauss_tria=NULL; - double viscosity; //viscosity - double oldviscosity; //viscosity - double newviscosity; //viscosity - double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double B[3][numdofp]; - double Bprime[3][numdofm]; - double L[2][numdofp]; - double D[3][3]={0.0}; // material matrix, simple scalar matrix. - double D_scalar; - double DL[2][2]={0.0}; //for basal drag - double DL_scalar; - double Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix - double Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point. - double Jdet; - double alpha2_list[3]; - double alpha2; - double viscosity_overshoot; /*Find penta on bed as pattyn must be coupled to the dofs on the bed: */ @@ -2019,5 +2011,5 @@ /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESp); + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); @@ -2034,4 +2026,8 @@ gauss->SynchronizeGaussTria(gauss_tria); + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss); + tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria); + this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); @@ -2039,8 +2035,4 @@ matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); - GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss); - tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=2*newviscosity*gauss->weight*Jdet; @@ -2139,9 +2131,11 @@ ElementMatrix* Penta::CreateKMatrixDiagnosticHutter(void){ + /*Constants*/ + const int numdof=NDOF2*NUMVERTICES; + /*Intermediaries*/ - const int numdof=NDOF2*NUMVERTICES; int connectivity[2]; + int i,i0,i1,j0,j1; double one0,one1; - int i,i0,i1,j0,j1; /*Initialize Element matrix and return if necessary*/ @@ -2238,29 +2232,22 @@ ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyeal3dViscous(void){ + /*Constants*/ const int numdof2d=2*NUMVERTICES2D; - int i,j,ig; - double xyz_list[NUMVERTICES][3]; + + /*Intermediaries */ + int i,j,ig; + double Jdet; + double viscosity, oldviscosity, newviscosity, viscosity_overshoot; + double epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + double xyz_list[NUMVERTICES][3]; + double B[3][numdof2d]; + double Bprime[3][numdof2d]; + double D[3][3]={0.0}; // material matrix, simple scalar matrix. + double D_scalar; + double Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point. + Tria* tria=NULL; + Penta* pentabase=NULL; GaussPenta *gauss=NULL; GaussTria *gauss_tria=NULL; - double viscosity, oldviscosity, newviscosity, viscosity_overshoot; - double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - double B[3][numdof2d]; - double Bprime[3][numdof2d]; - double L[2][numdof2d]; - double D[3][3]={0.0}; // material matrix, simple scalar matrix. - double D_scalar; - double DL[2][2]={0.0}; //for basal drag - double DL_scalar; - double Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point. - double Jdet; - double slope[2]={0.0}; - double slope_magnitude; - double alpha2_list[3]; - double alpha2; - double MAXSLOPE=.06; // 6 % - double MOUNTAINKEXPONENT=10; - Tria* tria=NULL; - Penta* pentabase=NULL; /*Find penta on bed as this is a macayeal elements: */ @@ -2288,7 +2275,7 @@ gauss->SynchronizeGaussTria(gauss_tria); + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria); tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria); - this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); @@ -2296,6 +2283,6 @@ matice->GetViscosity3d(&viscosity, &epsilon[0]); matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); D_scalar=2*newviscosity*gauss->weight*Jdet; for (i=0;i<3;i++) D[i][i]=D_scalar; @@ -2555,5 +2542,4 @@ double LprimeStokes[14][numdof2d]; double DLStokes[14][14]={0.0}; - double Ke_temp[27][27]={0.0}; //for the six nodes and the bubble double Ke_drag_gaussian[numdof2d][numdof2d]; Friction* friction=NULL; @@ -2613,8 +2599,6 @@ &Ke_drag_gaussian[0][0],0); - for(i=0;ivalues[i*numdof+j]+=Ke_temp[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_drag_gaussian[i][j]; + } /*Clean up and return*/ @@ -2641,14 +2625,16 @@ ElementMatrix* Penta::CreateKMatrixDiagnosticVertVolume(void){ - /* local declarations */ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - int i,j,ig; - double xyz_list[NUMVERTICES][3]; - GaussPenta *gauss=NULL; - double Ke_gg[numdof][numdof]={0.0}; - double Jdet; - double B[NDOF1][NUMVERTICES]; - double Bprime[NDOF1][NUMVERTICES]; - double DL_scalar; + + /*Intermediaries */ + int i,j,ig; + double Jdet; + double xyz_list[NUMVERTICES][3]; + double B[NDOF1][NUMVERTICES]; + double Bprime[NDOF1][NUMVERTICES]; + double DL_scalar; + double Ke_gg[numdof][numdof]={0.0}; + GaussPenta *gauss=NULL; /*Initialize Element matrix and return if necessary*/ @@ -2665,8 +2651,8 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); GetBVert(&B[0][0], &xyz_list[0][0], gauss); GetBprimeVert(&Bprime[0][0], &xyz_list[0][0], gauss); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); DL_scalar=gauss->weight*Jdet; @@ -2761,17 +2747,14 @@ ElementMatrix* Penta::CreateKMatrixThermalVolume(void){ - int i,j; - int ig; - int found=0; + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - GaussPenta *gauss=NULL; - double K[2][2]={0.0}; - double u,v,w; - double Ke_gaussian_conduct[numdof][numdof]; - double Ke_gaussian_advec[numdof][numdof]; - double Ke_gaussian_artdiff[numdof][numdof]; - double Ke_gaussian_transient[numdof][numdof]; + + /*Intermediaries */ + bool artdiff; + int i,j,ig,found=0; + double Jdet,u,v,w,epsvel; + double gravity,rho_ice,rho_water; + double heatcapacity,thermalconductivity,dt; + double xyz_list[NUMVERTICES][3]; double B[3][numdof]; double Bprime[3][numdof]; @@ -2783,6 +2766,5 @@ double D_scalar; double D[3][3]; - double l1l2l3[3]; - double tl1l2l3D[3]; + double K[2][2]={0.0}; double tBD[3][numdof]; double tBD_conduct[3][numdof]; @@ -2790,11 +2772,10 @@ double tBD_artdiff[3][numdof]; double tLD[numdof]; - double Jdet; - double gravity,rho_ice,rho_water; - double heatcapacity,thermalconductivity; - double mixed_layer_capacity,thermal_exchange_velocity; - double dt,epsvel; - bool artdiff; - Tria* tria=NULL; + double Ke_gaussian_conduct[numdof][numdof]; + double Ke_gaussian_advec[numdof][numdof]; + double Ke_gaussian_artdiff[numdof][numdof]; + double Ke_gaussian_transient[numdof][numdof]; + Tria* tria=NULL; + GaussPenta *gauss=NULL; /*Initialize Element matrix and return if necessary*/ @@ -3019,27 +3000,85 @@ ElementVector* Penta::CreatePVectorCouplingPattynStokesViscous(void){ - const int numdof=NUMVERTICES*NDOF4; - int i,j; - int ig; - double xyz_list_tria[NUMVERTICES2D][3]; - double xyz_list[NUMVERTICES][3]; - double bed_normal[3]; + /*Constants*/ + const int numdof=NUMVERTICES*NDOF4; + + /*Intermediaries */ + int i,j,ig; + int approximation; + double viscosity,Jdet; + double stokesreconditioning; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double dw[3]; + double xyz_list[NUMVERTICES][3]; + double l1l6[6]; //for the six nodes of the penta + double dh1dh6[3][6]; //for the six nodes of the penta GaussPenta *gauss=NULL; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double viscosity, w, alpha2_gauss; - double dw[3]; - double Pe_gaussian[24]={0.0}; //for the six nodes - double l1l6[6]; //for the six nodes of the penta - double dh1dh6[3][6]; //for the six nodes of the penta - double Jdet; - double Jdet2d; - Tria* tria=NULL; - Friction* friction=NULL; - double stokesreconditioning; - int approximation; - int analysis_type; /*Initialize Element vector and return if necessary*/ if(IsOnWater()) return NULL; + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation!=PattynStokesApproximationEnum) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); + Input* vzpattyn_input=inputs->GetInput(VzPattynEnum); ISSMASSERT(vzpattyn_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1(&l1l6[0], gauss); + GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss); + + vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss); + + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + for(i=0;ivalues[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]/2; + pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2; + 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; + pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i]; + } + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{1*/ +ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF4; + + /*Intermediaries*/ + int i,j,ig; + int approximation,analysis_type; + double Jdet,Jdet2d; + double stokesreconditioning; + double bed_normal[3]; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double viscosity, w, alpha2_gauss; + double dw[3]; + double xyz_list_tria[NUMVERTICES2D][3]; + double xyz_list[NUMVERTICES][3]; + double l1l6[6]; //for the six nodes of the penta + Tria* tria=NULL; + Friction* friction=NULL; + GaussPenta *gauss=NULL; + + /*Initialize Element vector and return if necessary*/ + if(IsOnWater() || !IsOnBed() || IsOnShelf()) return NULL; inputs->GetParameterValue(&approximation,ApproximationEnum); if(approximation!=PattynStokesApproximationEnum) return NULL; @@ -3055,78 +3094,5 @@ Input* vzpattyn_input=inputs->GetInput(VzPattynEnum); ISSMASSERT(vzpattyn_input); - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - GetNodalFunctionsP1(&l1l6[0], gauss); - GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss); - - vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - - for(i=0;ivalues[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]/2; - pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2; - 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; - pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i]; - } - } - - /*Clean up and return*/ - delete gauss; - return pe; -} -/*}}}*/ -/*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{1*/ -ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){ - - const int numdof=NUMVERTICES*NDOF4; - int i,j; - int ig; - double xyz_list_tria[NUMVERTICES2D][3]; - double xyz_list[NUMVERTICES][3]; - double bed_normal[3]; - GaussPenta *gauss=NULL; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double viscosity, w, alpha2_gauss; - double dw[3]; - double Pe_gaussian[24]={0.0}; //for the six nodes - double l1l6[6]; //for the six nodes of the penta - double dh1dh6[3][6]; //for the six nodes of the penta - double Jdet; - double Jdet2d; - Tria* tria=NULL; - Friction* friction=NULL; - double stokesreconditioning; - int approximation; - int analysis_type; - - /*Initialize Element vector and return if necessary*/ - if(IsOnWater() || !IsOnBed() || IsOnShelf()) return NULL; - inputs->GetParameterValue(&approximation,ApproximationEnum); - if(approximation!=PattynStokesApproximationEnum) return NULL; - ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); - Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); - Input* vzpattyn_input=inputs->GetInput(VzPattynEnum); ISSMASSERT(vzpattyn_input); - - - for(i=0;iGaussPoint(ig); - /*Get the Jacobian determinant */ GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - - /*Get L matrix : */ GetNodalFunctionsP1(l1l6, gauss); - /*Get normal vecyor to the bed */ + vzpattyn_input->GetParameterValue(&w, gauss); + vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss); + BedNormal(&bed_normal[0],xyz_list_tria); - - /*Get Viscosity*/ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - /*Get friction*/ friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); - - /*Get w and its derivatives*/ - vzpattyn_input->GetParameterValue(&w, gauss); - vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss); for(i=0;iGaussPoint(ig); + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1(l1l6, gauss); + thickness_input->GetParameterValue(&thickness, gauss); surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); - - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(l1l6, gauss); driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG(); @@ -3396,35 +3350,25 @@ ElementVector* Penta::CreatePVectorDiagnosticStokesViscous(void){ - int i,j; - int ig; - const int numdof=NUMVERTICES*NDOF4; - int numdof2d=NUMVERTICES2D*NDOF4; - double gravity,rho_ice,rho_water; - double xyz_list_tria[NUMVERTICES2D][3]; - double xyz_list[NUMVERTICES][3]; - double bed_normal[3]; - double bed; - GaussPenta *gauss=NULL; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double viscosity; - double water_pressure; - double Pe_temp[27]={0.0}; //for the six nodes and the bubble - double Pe_gaussian[27]={0.0}; //for the six nodes and the bubble - double Ke_temp[27][3]={0.0}; //for the six nodes and the bubble - double Pe_reduced[numdof]; //for the six nodes only - double Ke_gaussian[27][3]; - double l1l6[6]; //for the six nodes of the penta + /*Constants*/ + const int numdofbubble=NDOF4*NUMVERTICES+NDOF3*1; + + /*Intermediaries*/ + int i,j,ig; + int approximation; + double Jdet,viscosity; + double gravity,rho_ice,stokesreconditioning; + double xyz_list[NUMVERTICES][3]; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ double l1l7[7]; //for the six nodes and the bubble - double B[8][27]; - double B_prime[8][27]; + double B[8][numdofbubble]; + double B_prime[8][numdofbubble]; double B_prime_bubble[8][3]; - double Jdet; - double Jdet2d; double D[8][8]={0.0}; double D_scalar; - double tBD[27][8]; - Tria* tria=NULL; - double stokesreconditioning; - int approximation; + double tBD[numdofbubble][8]; + double Pe_gaussian[numdofbubble]={0.0}; //for the six nodes and the bubble + double Ke_temp[numdofbubble][3]={0.0}; //for the six nodes and the bubble + double Ke_gaussian[numdofbubble][3]; + GaussPenta *gauss=NULL; /*Initialize Element vector and return if necessary*/ @@ -3436,5 +3380,4 @@ /*Retrieve all inputs and parameters*/ this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); - rho_water=matpar->GetRhoWater(); rho_ice=matpar->GetRhoIce(); gravity=matpar->GetG(); @@ -3443,5 +3386,4 @@ Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); - Input* bed_input=inputs->GetInput(BedEnum); ISSMASSERT(bed_input); /* Start looping on the number of gaussian points: */ @@ -3451,40 +3393,31 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetBStokes(&B[0][0],&xyz_list[0][0],gauss); + GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss); + GetNodalFunctionsMINI(&l1l7[0], gauss); + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsMINI(&l1l7[0], gauss); - - /* Build gaussian vector */ + for(i=0;iweight*l1l7[i]; - } - - /*Add Pe_gaussian to terms in Pe_temp. Watch out for column orientation from matlab: */ - for(i=0;i<27;i++) Pe_temp[i]+=Pe_gaussian[i]; - - /*Get B and Bprime matrices: */ - GetBStokes(&B[0][0],&xyz_list[0][0],gauss); - GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss); + Pe_gaussian[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l7[i]; + } /*Get bubble part of Bprime */ for(i=0;i<8;i++) for(j=0;j<3;j++) B_prime_bubble[i][j]=B_prime[i][j+24]; - /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant - * onto this scalar matrix, so that we win some computational time: */ D_scalar=gauss->weight*Jdet; for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity; for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning; - /* Do the triple product tB*D*Bprime: */ - MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0); - MatrixMultiply(&tBD[0][0],27,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0); - - /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */ - for(i=0;i<27;i++) for(j=0;j<3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j]; + MatrixMultiply(&B[0][0],8,numdofbubble,1,&D[0][0],8,8,0,&tBD[0][0],0); + MatrixMultiply(&tBD[0][0],numdofbubble,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0); + + for(i=0;ivalues, &Ke_temp[0][0], &Pe_temp[0]); + ReduceVectorStokes(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]); /*Clean up and return*/ @@ -3496,35 +3429,15 @@ ElementVector* Penta::CreatePVectorDiagnosticStokesShelf(void){ - int i,j; - int ig; - const int numdof=NUMVERTICES*NDOF4; - int numdof2d=NUMVERTICES2D*NDOF4; - double gravity,rho_ice,rho_water; - double xyz_list_tria[NUMVERTICES2D][3]; - double xyz_list[NUMVERTICES][3]; - double bed_normal[3]; - double bed; - GaussPenta *gauss=NULL; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double viscosity; - double water_pressure; - double Pe_temp[27]={0.0}; //for the six nodes and the bubble - double Pe_gaussian[27]={0.0}; //for the six nodes and the bubble - double Ke_temp[27][3]={0.0}; //for the six nodes and the bubble - double Pe_reduced[numdof]; //for the six nodes only - double Ke_gaussian[27][3]; - double l1l6[6]; //for the six nodes of the penta - double l1l7[7]; //for the six nodes and the bubble - double B[8][27]; - double B_prime[8][27]; - double B_prime_bubble[8][3]; - double Jdet; - double Jdet2d; - double D[8][8]={0.0}; - double D_scalar; - double tBD[27][8]; - Tria* tria=NULL; - double stokesreconditioning; - int approximation; + /*Intermediaries*/ + int i,j,ig; + int approximation; + double gravity,rho_water,bed,water_pressure; + double xyz_list_tria[NUMVERTICES2D][3]; + double xyz_list[NUMVERTICES][3]; + double bed_normal[3]; + double l1l6[6]; //for the six nodes of the penta + double Jdet2d; + Tria* tria=NULL; + GaussPenta *gauss=NULL; /*Initialize Element vector and return if necessary*/ @@ -3535,20 +3448,10 @@ /*Retrieve all inputs and parameters*/ - this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); gravity=matpar->GetG(); GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); Input* bed_input=inputs->GetInput(BedEnum); ISSMASSERT(bed_input); - - for(i=0;iGetParameterValue(&bed, gauss); - GetNodalFunctionsP1(l1l6, gauss); - + BedNormal(&bed_normal[0],xyz_list_tria); water_pressure=gravity*rho_water*bed; - - BedNormal(&bed_normal[0],xyz_list_tria); for(i=0;ivalues[i*NDOF4+j]+=water_pressure*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j]; @@ -3605,18 +3507,16 @@ ElementVector* Penta::CreatePVectorDiagnosticVertVolume(void){ - const int numdof=NDOF1*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - int i; - int ig; + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries*/ + int i,ig; + int approximation; + double Jdet; + double xyz_list[NUMVERTICES][3]; + double dudx,dvdy,dwdz; + double du[3],dv[3],dw[3]; + double l1l6[6]; GaussPenta *gauss=NULL; - double Jdet; - double l1l6[6]; - Tria* tria=NULL; - double du[3]; - double dv[3]; - double dw[3]; - double dudx,dvdy,dwdz; - int approximation; /*Initialize Element vector and return if necessary*/ @@ -3640,4 +3540,7 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1(l1l6, gauss); + vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss); vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss); @@ -3650,7 +3553,4 @@ dvdy=dv[1]; - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(l1l6, gauss); - for (i=0;ivalues[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*l1l6[i]; } @@ -3735,34 +3635,19 @@ ElementVector* Penta::CreatePVectorThermalVolume(void){ + /*Constants*/ const int numdof=NUMVERTICES*NDOF1; - int i,j; - int ig; - int found=0; - double xyz_list[NUMVERTICES][3]; + + /*Intermediaries*/ + int i,j,ig,found=0; + int friction_type; + double Jdet,phi,dt; + double rho_ice,heatcapacity; + double viscosity,temperature; + double scalar_def,scalar_transient; + double temperature_list[NUMVERTICES]; + double xyz_list[NUMVERTICES][3]; + double L[numdof]; + double epsilon[6]; GaussPenta *gauss=NULL; - double temperature_list[NUMVERTICES]; - double temperature; - double gravity,rho_ice,rho_water; - double mixed_layer_capacity,heatcapacity; - double beta,meltingpoint,thermal_exchange_velocity; - int friction_type; - double P_terms[numdof]={0.0}; - double L[numdof]; - double l1l2l3[3]; - double basalfriction; - double epsilon[6]; - double epsilon_sqr[3][3]; - double epsilon_matrix[3][3]; - double Jdet; - double viscosity; - double epsilon_eff; - double phi; - double t_pmp; - double scalar; - double scalar_def; - double scalar_ocean; - double scalar_transient; - Tria* tria=NULL; - double dt; /*Initialize Element vector and return if necessary*/ @@ -3773,17 +3658,11 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); this->parameters->FindParam(&dt,DtEnum); - rho_water=matpar->GetRhoWater(); rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); heatcapacity=matpar->GetHeatCapacity(); - beta=matpar->GetBeta(); - meltingpoint=matpar->GetMeltingPoint(); Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); Input* temperature_input=NULL; - if (dt){ - temperature_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(inputs); - } + if (dt) temperature_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(inputs); /* Start looping on the number of gaussian points: */ @@ -3793,8 +3672,9 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1(&L[0], gauss); + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctionsP1(&L[0], gauss); GetPhi(&phi, &epsilon[0], viscosity); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5945) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5946) @@ -3038,5 +3038,4 @@ /*Intermediaries */ int i,j,ig; - int* doflist=NULL; double xyz_list[NUMVERTICES][3]; double surface_normal[3]; @@ -3274,5 +3273,4 @@ /*Intermediaries */ int i,j,ig,dim; - int* doflist=NULL; double xyz_list[NUMVERTICES][3]; double Jdettria,dt,vx,vy; @@ -3468,5 +3466,4 @@ /*Intermediaries */ int i,j,ig; - int* doflist=NULL; double xyz_list[NUMVERTICES][3]; double dhdt_g,melting_g,accumulation_g,Jdettria; @@ -3513,5 +3510,4 @@ /*Intermediaries */ int i,j,ig; - int* doflist=NULL; double xyz_list[NUMVERTICES][3]; double melting_g,accumulation_g,dhdt_g,Jdettria; @@ -3558,5 +3554,4 @@ /*Intermediaries */ int i,j,ig; - int* doflist=NULL; double xyz_list[NUMVERTICES][3]; double Jdettria,accumulation_g,melting_g; @@ -3703,5 +3698,5 @@ for (i=0;iweight*l1l2l3[i]; + pe->values[i*NDOF2+j]+=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss->weight*l1l2l3[i]; } } @@ -3710,9 +3705,8 @@ for (i=0;iweight*l1l2l3[i]; + pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l2l3[i]; } } } - for(i=0;ivalues[i]+=pe_g_gaussian[i]; } @@ -3730,11 +3724,9 @@ /*Intermediaries */ int i,ig; - int* doflist=NULL; double Jdet; double thickness,thicknessobs,weight; double xyz_list[NUMVERTICES][3]; double l1l2l3[3]; - double pe_g_gaussian[numdof]; - GaussTria* gauss=NULL; + GaussTria* gauss=NULL; /*Initialize Element vector and return if necessary*/ @@ -3774,26 +3766,21 @@ /*Constants*/ const int numdof=NDOF2*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - double vx_list[NUMVERTICES]; - double vy_list[NUMVERTICES]; - double obs_vx_list[NUMVERTICES]; - double obs_vy_list[NUMVERTICES]; - double dux_list[NUMVERTICES]; - double duy_list[NUMVERTICES]; - double weights_list[NUMVERTICES]; - int i; - int ig; + + /*Intermediaries */ + int i,ig,response; + double Jdet; + double obs_velocity_mag,velocity_mag; + double dux,duy,meanvel,epsvel; + double scalex=0,scaley=0,scale=0,S=0; + double xyz_list[NUMVERTICES][3]; + double vx_list[NUMVERTICES]; + double vy_list[NUMVERTICES]; + double obs_vx_list[NUMVERTICES]; + double obs_vy_list[NUMVERTICES]; + double dux_list[NUMVERTICES]; + double duy_list[NUMVERTICES]; + double weights_list[NUMVERTICES]; + double l1l2l3[3]; GaussTria* gauss=NULL; - double obs_velocity_mag,velocity_mag; - double dux,duy; - double meanvel, epsvel; - double pe_g_gaussian[numdof]; - double Jdet; - double l1l2l3[3]; - double scalex=0; - double scaley=0; - double scale=0; - double S=0; - int response; /*Initialize Element vector and return if necessary*/ @@ -3946,9 +3933,7 @@ for (i=0;iweight*l1l2l3[i]; - pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss->weight*l1l2l3[i]; + pe->values[i*NDOF2+0]+=dux*Jdet*gauss->weight*l1l2l3[i]; + pe->values[i*NDOF2+1]+=duy*Jdet*gauss->weight*l1l2l3[i]; } - - for(i=0;ivalues[i]+=pe_g_gaussian[i]; } @@ -3961,28 +3946,22 @@ ElementVector* Tria::CreatePVectorAdjointStokes(void){ - const int numdof=NDOF4*NUMVERTICES; - int i; - int ig; - double xyz_list[NUMVERTICES][3]; - double vx_list[NUMVERTICES]; - double vy_list[NUMVERTICES]; - double obs_vx_list[NUMVERTICES]; - double obs_vy_list[NUMVERTICES]; - double dux_list[NUMVERTICES]; - double duy_list[NUMVERTICES]; - double weights_list[NUMVERTICES]; + /*Intermediaries */ + int i,ig; + int fit=-1; + int response; + double Jdet; + double obs_velocity_mag,velocity_mag; + double dux,duy,meanvel,epsvel; + double scalex=0,scaley=0,scale=0,S=0; + double xyz_list[NUMVERTICES][3]; + double vx_list[NUMVERTICES]; + double vy_list[NUMVERTICES]; + double obs_vx_list[NUMVERTICES]; + double obs_vy_list[NUMVERTICES]; + double dux_list[NUMVERTICES]; + double duy_list[NUMVERTICES]; + double weights_list[NUMVERTICES]; + double l1l2l3[3]; GaussTria* gauss=NULL; - double obs_velocity_mag,velocity_mag; - double dux,duy; - double meanvel, epsvel; - double pe_g[numdof]={0.0}; - double Jdet; - double l1l2l3[3]; - double scalex=0; - double scaley=0; - double scale=0; - double S=0; - int fit=-1; - int response; /*Initialize Element vector and return if necessary*/ @@ -4148,14 +4127,10 @@ ElementVector* Tria::CreatePVectorDiagnosticHutter(void){ - /*Collapsed formulation: */ - int i; - const int numdofs=NDOF2*NUMVERTICES; - int* doflist=NULL; - double constant_part,ub,vb; - double rho_ice,gravity,n,B; - double slope[2]; - double thickness; - double slope2; - int connectivity; + /*Intermediaries */ + int i,connectivity; + double constant_part,ub,vb; + double rho_ice,gravity,n,B; + double slope2,thickness; + double slope[2]; GaussTria* gauss=NULL; @@ -4216,20 +4191,14 @@ ElementVector* Tria::CreatePVectorPrognostic_CG(void){ - /* local declarations */ - int i,j; - - /* node data: */ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - int numberofdofspernode=1; - int ig; + + /*Intermediaries */ + int i,j,ig; + double Jdettria,dt; + double accumulation_g,melting_g,thickness_g; + double xyz_list[NUMVERTICES][3]; + double L[NUMVERTICES]; GaussTria* gauss=NULL; - double L[NUMVERTICES]; - double Jdettria; - double accumulation_g; - double melting_g; - double thickness_g; - double dt; /*Initialize Element vector and return if necessary*/ @@ -4251,5 +4220,5 @@ GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); accumulation_input->GetParameterValue(&accumulation_g,gauss); @@ -4268,18 +4237,14 @@ ElementVector* Tria::CreatePVectorPrognostic_DG(void){ - /* local declarations */ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - int i,j; - int ig; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - int numberofdofspernode=1; + + /*Intermediaries */ + int i,j,ig; + double Jdettria,dt; + double accumulation_g,melting_g,thickness_g; + double xyz_list[NUMVERTICES][3]; + double L[NUMVERTICES]; GaussTria* gauss=NULL; - double L[NUMVERTICES]; - double Jdettria; - double accumulation_g; - double melting_g; - double thickness_g; - double dt; /*Initialize Element vector and return if necessary*/ @@ -4301,5 +4266,5 @@ GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); + GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); accumulation_input->GetParameterValue(&accumulation_g,gauss); @@ -4318,16 +4283,15 @@ ElementVector* Tria::CreatePVectorSlope(void){ - /* node data: */ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - int i,j; - int ig; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; + + /*Intermediaries */ + int i,j,ig; + int analysis_type; + double Jdet; + double xyz_list[NUMVERTICES][3]; + double slope[2]; + double l1l2l3[3]; GaussTria* gauss=NULL; - double Jdet; - double l1l2l3[3]; - double pe_g_gaussian[numdof]; - double slope[2]; - int analysis_type; /*Initialize Element vector and return if necessary*/ @@ -4352,23 +4316,15 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(l1l2l3, gauss); + slope_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /*Get nodal functions: */ - GetNodalFunctions(l1l2l3, gauss); - - /*Build pe_g_gaussian vector: */ if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){ - for(i=0;iweight*slope[0]*l1l2l3[i]; + for(i=0;ivalues[i]+=Jdet*gauss->weight*slope[0]*l1l2l3[i]; } if ( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){ - for(i=0;iweight*slope[1]*l1l2l3[i]; + for(i=0;ivalues[i]+=Jdet*gauss->weight*slope[1]*l1l2l3[i]; } - - /*Add pe_g_gaussian vector to pe_g: */ - for(i=0;ivalues[i]+=pe_g_gaussian[i]; - } @@ -4381,22 +4337,16 @@ ElementVector* Tria::CreatePVectorThermalShelf(void){ + /*Constants*/ const int numdof=NUMVERTICES*NDOF1; - int i; - double xyz_list[NUMVERTICES][3]; - double mixed_layer_capacity; - double thermal_exchange_velocity; - double rho_water; - double rho_ice; - double heatcapacity; - double beta; - double meltingpoint; - double dt; - double pressure; - int ig; + + /*Intermediaries */ + int i,ig; + double Jdet; + double mixed_layer_capacity,thermal_exchange_velocity; + double rho_ice,rho_water,pressure,dt,scalar_ocean; + double meltingpoint,beta,heatcapacity,t_pmp; + double xyz_list[NUMVERTICES][3]; + double l1l2l3[NUMVERTICES]; GaussTria* gauss=NULL; - double Jdet; - double l1l2l3[NUMVERTICES]; - double t_pmp; - double scalar_ocean; /*Initialize Element vector and return if necessary*/ @@ -4429,7 +4379,5 @@ scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice); - if(dt){ - scalar_ocean=dt*scalar_ocean; - } + if(dt) scalar_ocean=dt*scalar_ocean; for(i=0;ivalues[i]+=scalar_ocean*l1l2l3[i]; @@ -4444,24 +4392,19 @@ ElementVector* Tria::CreatePVectorThermalSheet(void){ + /*Constants*/ const int numdof=NUMVERTICES*NDOF1; - int i; - int ig; + + /*Intermediaries */ + int i,ig; + int analysis_type,drag_type; double xyz_list[NUMVERTICES][3]; - double rho_ice; - double heatcapacity; - double dt; - double pressure_list[3]; - double pressure; - int drag_type; - double basalfriction; - Friction* friction=NULL; - double alpha2,vx,vy; - double geothermalflux_value; + double Jdet,dt; + double rho_ice,heatcapacity,geothermalflux_value; + double basalfriction,alpha2,vx,vy,pressure; + double pressure_list[3]; + double scalar; + double l1l2l3[NUMVERTICES]; + Friction* friction=NULL; GaussTria* gauss=NULL; - double Jdet; - double P_terms[numdof]={0.0}; - double l1l2l3[NUMVERTICES]; - double scalar; - int analysis_type; /*Initialize Element vector and return if necessary*/ @@ -4498,7 +4441,5 @@ scalar=gauss->weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice); - if(dt){ - scalar=dt*scalar; - } + if(dt) scalar=dt*scalar; for(i=0;ivalues[i]+=scalar*l1l2l3[i];