Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 15739) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 15740) @@ -221,9 +221,4 @@ switch(analysis_type){ - #ifdef _HAVE_DIAGNOSTIC_ - case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: - Ke=PenaltyCreateKMatrixDiagnosticFS(kmax); - break; - #endif #ifdef _HAVE_THERMAL_ case ThermalAnalysisEnum: @@ -448,38 +443,4 @@ } /*}}}*/ -#ifdef _HAVE_DIAGNOSTIC_ -/*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticFS {{{*/ -ElementMatrix* Pengrid::PenaltyCreateKMatrixDiagnosticFS(IssmDouble kmax){ - - const int numdof = NUMVERTICES *NDOF4; - IssmDouble slope[2]; - IssmDouble penalty_offset; - int approximation; - - Penta* penta=(Penta*)element; - - /*Initialize Element vector and return if necessary*/ - penta->inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL; - ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters,FSApproximationEnum); - - /*Retrieve all inputs and parameters*/ - parameters->FindParam(&penalty_offset,DiagnosticPenaltyFactorEnum); - penta->GetInputValue(&slope[0],node,BedSlopeXEnum); - penta->GetInputValue(&slope[1],node,BedSlopeYEnum); - - /*Create elementary matrix: add penalty to constrain wb (wb=ub*db/dx+vb*db/dy)*/ - Ke->values[2*NDOF4+0]=-slope[0]*kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[2*NDOF4+1]=-slope[1]*kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[2*NDOF4+2]= kmax*pow((IssmDouble)10,penalty_offset); - - /*Transform Coordinate System*/ - TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,XYZEnum); - - /*Clean up and return*/ - return Ke; -} -/*}}}*/ -#endif #ifdef _HAVE_THERMAL_ /*FUNCTION Pengrid::ConstraintActivateThermal {{{*/ @@ -574,5 +535,5 @@ /*Add penalty load*/ if (temperaturevalues[0]=kmax*pow((IssmDouble)10,penalty_factor); + Ke->values[0]=kmax*pow(10.,penalty_factor); } @@ -594,5 +555,5 @@ parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); - Ke->values[0]=kmax*pow((IssmDouble)10,penalty_factor); + Ke->values[0]=kmax*pow(10.,penalty_factor); /*Clean up and return*/ @@ -635,6 +596,6 @@ } else{ - if (reCast(dt)) pe->values[0]=melting_offset*pow((IssmDouble)10,penalty_factor)*(temperature-t_pmp)/dt; - else pe->values[0]=melting_offset*pow((IssmDouble)10,penalty_factor)*(temperature-t_pmp); + if (reCast(dt)) pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp)/dt; + else pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp); } @@ -664,5 +625,5 @@ t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; - pe->values[0]=kmax*pow((IssmDouble)10,penalty_factor)*t_pmp; + pe->values[0]=kmax*pow(10.,penalty_factor)*t_pmp; /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 15739) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 15740) @@ -87,7 +87,4 @@ /*}}}*/ /*Pengrid management {{{*/ - #ifdef _HAVE_DIAGNOSTIC_ - ElementMatrix* PenaltyCreateKMatrixDiagnosticFS(IssmDouble kmax); - #endif #ifdef _HAVE_THERMAL_ ElementMatrix* PenaltyCreateKMatrixThermal(IssmDouble kmax); Index: /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp (revision 15739) +++ /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp (revision 15740) @@ -340,13 +340,13 @@ //Create elementary matrix: add penalty to - Ke->values[0*numdof+0]=+kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[0*numdof+2]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[2*numdof+0]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[2*numdof+2]=+kmax*pow((IssmDouble)10.0,penalty_offset); - - Ke->values[1*numdof+1]=+kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[1*numdof+3]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[3*numdof+1]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[3*numdof+3]=+kmax*pow((IssmDouble)10.0,penalty_offset); + Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset); + Ke->values[0*numdof+2]=-kmax*pow(10.,penalty_offset); + Ke->values[2*numdof+0]=-kmax*pow(10.,penalty_offset); + Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset); + + Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset); + Ke->values[1*numdof+3]=-kmax*pow(10.,penalty_offset); + Ke->values[3*numdof+1]=-kmax*pow(10.,penalty_offset); + Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset); /*Clean up and return*/ @@ -367,23 +367,23 @@ //Create elementary matrix: add penalty to - Ke->values[0*numdof+0]=+kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[0*numdof+4]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[4*numdof+0]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[4*numdof+4]=+kmax*pow((IssmDouble)10.0,penalty_offset); - - Ke->values[1*numdof+1]=+kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[1*numdof+5]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[5*numdof+1]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[5*numdof+5]=+kmax*pow((IssmDouble)10.0,penalty_offset); - - Ke->values[2*numdof+2]=+kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[2*numdof+6]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[6*numdof+2]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[6*numdof+6]=+kmax*pow((IssmDouble)10.0,penalty_offset); - - Ke->values[3*numdof+3]=+kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[3*numdof+7]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[7*numdof+3]=-kmax*pow((IssmDouble)10.0,penalty_offset); - Ke->values[7*numdof+7]=+kmax*pow((IssmDouble)10.0,penalty_offset); + Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset); + Ke->values[0*numdof+4]=-kmax*pow(10.,penalty_offset); + Ke->values[4*numdof+0]=-kmax*pow(10.,penalty_offset); + Ke->values[4*numdof+4]=+kmax*pow(10.,penalty_offset); + + Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset); + Ke->values[1*numdof+5]=-kmax*pow(10.,penalty_offset); + Ke->values[5*numdof+1]=-kmax*pow(10.,penalty_offset); + Ke->values[5*numdof+5]=+kmax*pow(10.,penalty_offset); + + Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset); + Ke->values[2*numdof+6]=-kmax*pow(10.,penalty_offset); + Ke->values[6*numdof+2]=-kmax*pow(10.,penalty_offset); + Ke->values[6*numdof+6]=+kmax*pow(10.,penalty_offset); + + Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset); + Ke->values[3*numdof+7]=-kmax*pow(10.,penalty_offset); + Ke->values[7*numdof+3]=-kmax*pow(10.,penalty_offset); + Ke->values[7*numdof+7]=+kmax*pow(10.,penalty_offset); /*Clean up and return*/ @@ -404,8 +404,8 @@ //Create elementary matrix: add penalty to - Ke->values[0*numdof+0]=+kmax*pow((IssmDouble)10.0,penalty_factor); - Ke->values[0*numdof+1]=-kmax*pow((IssmDouble)10.0,penalty_factor); - Ke->values[1*numdof+0]=-kmax*pow((IssmDouble)10.0,penalty_factor); - Ke->values[1*numdof+1]=+kmax*pow((IssmDouble)10.0,penalty_factor); + Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_factor); + Ke->values[0*numdof+1]=-kmax*pow(10.,penalty_factor); + Ke->values[1*numdof+0]=-kmax*pow(10.,penalty_factor); + Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_factor); /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 15739) +++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 15740) @@ -576,13 +576,13 @@ if(shelf){ /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */ - pressure=rho_ice*gravity*pow(thickness,(IssmDouble)2)/(IssmDouble)2 - rho_water*gravity*pow(bed,(IssmDouble)2)/(IssmDouble)2; + pressure=rho_ice*gravity*pow(thickness,2)/2.- rho_water*gravity*pow(bed,2)/2.; } else{ //We are on an icesheet, we assume the water column fills the entire front: */ - pressure=rho_ice*gravity*pow(thickness,(IssmDouble)2)/(IssmDouble)2 - rho_water*gravity*pow(thickness,(IssmDouble)2)/(IssmDouble)2; + pressure=rho_ice*gravity*pow(thickness,2)/2.- rho_water*gravity*pow(thickness,2)/2.; } } else if(fill==AirEnum){ - pressure=rho_ice*gravity*pow(thickness,(IssmDouble)2)/(IssmDouble)2; //icefront on an ice sheet, pressure imbalance ice vs air. + pressure=rho_ice*gravity*pow(thickness,2)/2.; //icefront on an ice sheet, pressure imbalance ice vs air. } else if(fill==IceEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice) @@ -593,7 +593,7 @@ if(!shelf) _error_("fill type " << fill << " not supported on ice sheets yet."); - pressure_litho=rho_ice*gravity*pow(thickness,(IssmDouble)2)/(IssmDouble)2; + pressure_litho=rho_ice*gravity*pow(thickness,2)/2.; pressure_air=0; - pressure_melange=rho_ice*gravity*pow(fraction*thickness,(IssmDouble)2)/(IssmDouble)2; + pressure_melange=rho_ice*gravity*pow(fraction*thickness,2)/2.; pressure_water=1.0/2.0*rho_water*gravity* ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) ); @@ -684,5 +684,5 @@ /*increase melange fraction: */ this->fraction+=fractionincrement; - if (this->fraction>1)this->fraction=(IssmDouble)1.0; + if (this->fraction>1)this->fraction=1.; //_printf_("riftfront " << this->Id() << " fraction: " << this->fraction << "\n"); } Index: /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp (revision 15739) +++ /issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp (revision 15740) @@ -249,5 +249,5 @@ else{ if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ - viscosity=0.5*pow((IssmDouble)10,(IssmDouble)14); + viscosity=0.5*pow(10.,14); } else{ @@ -261,5 +261,5 @@ if(A==0){ /*Maxiviscositym viscosity for 0 shear areas: */ - viscosity=2.5*pow(10.,17.); + viscosity=2.5*pow(10.,17); } else{ @@ -317,5 +317,5 @@ if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && (epsilon[3]==0) && (epsilon[4]==0)){ - viscosity3d=0.5*pow((IssmDouble)10,(IssmDouble)14); + viscosity3d=0.5*pow(10.,14); } else{ @@ -332,5 +332,5 @@ if(A==0){ /*Maxiviscosity3dm viscosity for 0 shear areas: */ - viscosity3d=2.25*pow((IssmDouble)10,(IssmDouble)17); + viscosity3d=2.25*pow(10.,17); } else{ @@ -378,5 +378,5 @@ /*Get B and n*/ - eps0=pow((IssmDouble)10,(IssmDouble)-27); + eps0=pow(10.,(IssmDouble)-27); n=GetN(); Z=GetZ(); @@ -390,5 +390,5 @@ if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){ - viscosity3d=0.5*pow((IssmDouble)10,(IssmDouble)14); + viscosity3d=0.5*pow(10.,14); } else{ @@ -406,5 +406,5 @@ if(A==0){ /*Maxiviscosity3dm viscosity for 0 shear areas: */ - viscosity3d=2.25*pow((IssmDouble)10,(IssmDouble)17); + viscosity3d=2.25*pow(10.,17); } else{ @@ -459,5 +459,5 @@ if(A==0){ /*Maximum viscosity_complement for 0 shear areas: */ - viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17); + viscosity_complement=2.25*pow(10.,17); } else{ @@ -468,5 +468,5 @@ } else{ - viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17); + viscosity_complement=4.5*pow(10.,17); } @@ -515,5 +515,5 @@ if(A==0){ /*Maximum viscosity_complement for 0 shear areas: */ - viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17); + viscosity_complement=2.25*pow(10.,17); } else{ @@ -524,5 +524,5 @@ } else{ - viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17); + viscosity_complement=4.5*pow(10.,17); } @@ -552,5 +552,5 @@ if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && (epsilon[3]==0) && (epsilon[4]==0)){ - mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14); + mu_prime=0.5*pow(10.,14); } else{ @@ -586,5 +586,5 @@ if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){ - mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14); + mu_prime=0.5*pow(10.,14); } else{ @@ -620,5 +620,5 @@ if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ - mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14); + mu_prime=0.5*pow(10.,14); } else{