Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5909) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5910) @@ -2405,4 +2405,5 @@ ElementMatrix* Tria::CreateKMatrixBalancedthickness_CG(void){ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; @@ -3399,21 +3400,18 @@ ElementMatrix* Tria::CreateKMatrixThermal(void){ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - int i,j; - int ig; - double xyz_list[NUMVERTICES][3]; - double mixed_layer_capacity; - double thermal_exchange_velocity; - double rho_water; - double rho_ice; - double heatcapacity; + + /*Intermediaries */ + int i,j,ig; + double mixed_layer_capacity,thermal_exchange_velocity; + double rho_ice,rho_water,heatcapacity; + double Jdet,dt; + double xyz_list[NUMVERTICES][3]; + double l1l2l3[NUMVERTICES]; + double tl1l2l3D[3]; + double D_scalar; + double Ke_gaussian[numdof][numdof]={0.0}; GaussTria *gauss=NULL; - double Jdet; - double K_terms[numdof][numdof]={0.0}; - double Ke_gaussian[numdof][numdof]={0.0}; - double l1l2l3[NUMVERTICES]; - double tl1l2l3D[3]; - double D_scalar; - double dt; /*Initialize Element matrix and return if necessary*/ @@ -3436,7 +3434,7 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss); GetNodalFunctions(&l1l2l3[0], gauss); - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss); D_scalar=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice); if(dt) D_scalar=dt*D_scalar; @@ -3664,27 +3662,23 @@ void Tria::CreatePVectorDiagnosticMacAyeal( Vec pg, Vec pf){ - int i,j,ig; - - /* node data: */ + /*Constants*/ const int numdof=NDOF2*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - - /* parameters: */ - double plastic_stress; - double slope[2]; - double driving_stress_baseline; - GaussTria* gauss=NULL; - - double Jdet; - double l1l2l3[3]; - double pe_g[numdof]={0.0}; - double pe_g_gaussian[numdof]; - double thickness; - int drag_type; - - /*intermediary: */ + + /*Intermediaries */ + int i,j,ig,drag_type; + double plastic_stress,driving_stress_baseline,thickness; + double Jdet; + double xyz_list[NUMVERTICES][3]; + double slope[2]; + double l1l2l3[3]; + double pe_g[numdof]={0.0}; + double pe_g_gaussian[numdof]; + GaussTria* gauss=NULL; ElementVector* pe=NULL; - /*retrieve inputs :*/ + /*First, if we are on water, return empty vector: */ + if(IsOnWater()) return; + + /*Retrieve all Inputs and parameters: */ inputs->GetParameterValue(&drag_type,DragTypeEnum); Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); @@ -3692,8 +3686,5 @@ Input* drag_input=inputs->GetInput(DragCoefficientEnum);ISSMASSERT(drag_input); - /*First, if we are on water, return empty vector: */ - if(IsOnWater()) return; - - /* Get node coordinates and dof list: */ + /*Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); @@ -3703,4 +3694,7 @@ gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(l1l2l3, gauss); thickness_input->GetParameterValue(&thickness,gauss); @@ -3709,15 +3703,6 @@ /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the * element itself: */ - if(drag_type==1){ - drag_input->GetParameterValue(&plastic_stress,gauss); - } - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /*Get nodal functions: */ - GetNodalFunctions(l1l2l3, gauss); - - /*Compute driving stress: */ + if(drag_type==1) drag_input->GetParameterValue(&plastic_stress,gauss); + driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness; @@ -3738,22 +3723,13 @@ } - /*Add pe_g_gaussian vector to pe_g: */ for( i=0; iparameters,MacAyealApproximationEnum); - - /*Add pe_g values to pe element stifness load: */ pe->AddValues(&pe_g[0]); - - /*Add pe element load vector to global load vector: */ pe->AddToGlobal(pg,pf); /*Free ressources:*/ delete pe; - - /*Free ressources:*/ delete gauss; @@ -3764,27 +3740,17 @@ void Tria::CreatePVectorAdjointBalancedthickness(Vec p_g){ - int i; - - /* node data: */ + /*Constants*/ const int numdof=1*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - - /* gaussian points: */ - int ig; + + /*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]; + double pe_g[numdof] ={0.0}; GaussTria* gauss=NULL; - - /* parameters: */ - double thickness,thicknessobs,weight; - - /*element vector : */ - double pe_g[numdof]={0.0}; - double pe_g_gaussian[numdof]; - - /* Jacobian: */ - double Jdet; - - /*nodal functions: */ - double l1l2l3[3]; /* Get node coordinates and dof list: */ @@ -3792,5 +3758,5 @@ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - /*Retrieve all inputs we will be needing: */ + /*Retrieve all Inputs and parameters: */ Input* thickness_input =inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input); @@ -3803,11 +3769,7 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /* Get nodal functions value at gaussian point:*/ GetNodalFunctions(l1l2l3, gauss); - /*Get parameters at gauss point*/ thickness_input->GetParameterValue(&thickness, gauss); thicknessobs_input->GetParameterValue(&thicknessobs, gauss); @@ -3819,10 +3781,7 @@ /*Add pe_g_gaussian vector to pe_g: */ - for( i=0; i