Changeset 15740


Ignore:
Timestamp:
08/07/13 10:57:10 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: no need to cast numbers to IssmDouble (needed for AD, these variables are NOT active)

Location:
issm/trunk-jpl/src/c/classes
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp

    r15654 r15740  
    221221
    222222        switch(analysis_type){
    223                 #ifdef _HAVE_DIAGNOSTIC_
    224                 case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
    225                         Ke=PenaltyCreateKMatrixDiagnosticFS(kmax);
    226                         break;
    227                 #endif
    228223                #ifdef _HAVE_THERMAL_
    229224                case ThermalAnalysisEnum:
    … …  
    448443}
    449444/*}}}*/
    450 #ifdef _HAVE_DIAGNOSTIC_
    451 /*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticFS {{{*/
    452 ElementMatrix* Pengrid::PenaltyCreateKMatrixDiagnosticFS(IssmDouble kmax){
    453 
    454         const int numdof = NUMVERTICES *NDOF4;
    455         IssmDouble    slope[2];
    456         IssmDouble    penalty_offset;
    457         int       approximation;
    458 
    459         Penta* penta=(Penta*)element;
    460 
    461         /*Initialize Element vector and return if necessary*/
    462         penta->inputs->GetInputValue(&approximation,ApproximationEnum);
    463         if(approximation!=FSApproximationEnum &&  approximation!=HOFSApproximationEnum) return NULL;
    464         ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters,FSApproximationEnum);
    465 
    466         /*Retrieve all inputs and parameters*/
    467         parameters->FindParam(&penalty_offset,DiagnosticPenaltyFactorEnum);
    468         penta->GetInputValue(&slope[0],node,BedSlopeXEnum);
    469         penta->GetInputValue(&slope[1],node,BedSlopeYEnum);
    470 
    471         /*Create elementary matrix: add penalty to constrain wb (wb=ub*db/dx+vb*db/dy)*/
    472         Ke->values[2*NDOF4+0]=-slope[0]*kmax*pow((IssmDouble)10.0,penalty_offset);
    473         Ke->values[2*NDOF4+1]=-slope[1]*kmax*pow((IssmDouble)10.0,penalty_offset);
    474         Ke->values[2*NDOF4+2]= kmax*pow((IssmDouble)10,penalty_offset);
    475 
    476         /*Transform Coordinate System*/
    477         TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,XYZEnum);
    478 
    479         /*Clean up and return*/
    480         return Ke;
    481 }
    482 /*}}}*/
    483 #endif
    484445#ifdef _HAVE_THERMAL_
    485446/*FUNCTION Pengrid::ConstraintActivateThermal {{{*/
    … …  
    574535        /*Add penalty load*/
    575536        if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp, instead of using spcs, use penalties
    576                 Ke->values[0]=kmax*pow((IssmDouble)10,penalty_factor);
     537                Ke->values[0]=kmax*pow(10.,penalty_factor);
    577538        }
    578539
    … …  
    594555        parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
    595556
    596         Ke->values[0]=kmax*pow((IssmDouble)10,penalty_factor);
     557        Ke->values[0]=kmax*pow(10.,penalty_factor);
    597558
    598559        /*Clean up and return*/
    … …  
    635596        }
    636597        else{
    637                 if (reCast<bool>(dt)) pe->values[0]=melting_offset*pow((IssmDouble)10,penalty_factor)*(temperature-t_pmp)/dt;
    638                 else    pe->values[0]=melting_offset*pow((IssmDouble)10,penalty_factor)*(temperature-t_pmp);
     598                if (reCast<bool>(dt)) pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp)/dt;
     599                else    pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp);
    639600        }
    640601
    … …  
    664625        t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
    665626
    666         pe->values[0]=kmax*pow((IssmDouble)10,penalty_factor)*t_pmp;
     627        pe->values[0]=kmax*pow(10.,penalty_factor)*t_pmp;
    667628
    668629        /*Clean up and return*/

  • issm/trunk-jpl/src/c/classes/Loads/Pengrid.h

    r15564 r15740  
    8787                /*}}}*/
    8888                /*Pengrid management {{{*/
    89                 #ifdef _HAVE_DIAGNOSTIC_
    90                 ElementMatrix* PenaltyCreateKMatrixDiagnosticFS(IssmDouble kmax);
    91                 #endif
    9289                #ifdef _HAVE_THERMAL_
    9390                ElementMatrix* PenaltyCreateKMatrixThermal(IssmDouble kmax);

  • issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp

    r15567 r15740  
    340340
    341341        //Create elementary matrix: add penalty to
    342         Ke->values[0*numdof+0]=+kmax*pow((IssmDouble)10.0,penalty_offset);
    343         Ke->values[0*numdof+2]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    344         Ke->values[2*numdof+0]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    345         Ke->values[2*numdof+2]=+kmax*pow((IssmDouble)10.0,penalty_offset);
    346 
    347         Ke->values[1*numdof+1]=+kmax*pow((IssmDouble)10.0,penalty_offset);
    348         Ke->values[1*numdof+3]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    349         Ke->values[3*numdof+1]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    350         Ke->values[3*numdof+3]=+kmax*pow((IssmDouble)10.0,penalty_offset);
     342        Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset);
     343        Ke->values[0*numdof+2]=-kmax*pow(10.,penalty_offset);
     344        Ke->values[2*numdof+0]=-kmax*pow(10.,penalty_offset);
     345        Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset);
     346
     347        Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset);
     348        Ke->values[1*numdof+3]=-kmax*pow(10.,penalty_offset);
     349        Ke->values[3*numdof+1]=-kmax*pow(10.,penalty_offset);
     350        Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset);
    351351
    352352        /*Clean up and return*/
    … …  
    367367
    368368        //Create elementary matrix: add penalty to
    369         Ke->values[0*numdof+0]=+kmax*pow((IssmDouble)10.0,penalty_offset);
    370         Ke->values[0*numdof+4]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    371         Ke->values[4*numdof+0]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    372         Ke->values[4*numdof+4]=+kmax*pow((IssmDouble)10.0,penalty_offset);
    373 
    374         Ke->values[1*numdof+1]=+kmax*pow((IssmDouble)10.0,penalty_offset);
    375         Ke->values[1*numdof+5]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    376         Ke->values[5*numdof+1]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    377         Ke->values[5*numdof+5]=+kmax*pow((IssmDouble)10.0,penalty_offset);
    378 
    379         Ke->values[2*numdof+2]=+kmax*pow((IssmDouble)10.0,penalty_offset);
    380         Ke->values[2*numdof+6]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    381         Ke->values[6*numdof+2]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    382         Ke->values[6*numdof+6]=+kmax*pow((IssmDouble)10.0,penalty_offset);
    383 
    384         Ke->values[3*numdof+3]=+kmax*pow((IssmDouble)10.0,penalty_offset);
    385         Ke->values[3*numdof+7]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    386         Ke->values[7*numdof+3]=-kmax*pow((IssmDouble)10.0,penalty_offset);
    387         Ke->values[7*numdof+7]=+kmax*pow((IssmDouble)10.0,penalty_offset);
     369        Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset);
     370        Ke->values[0*numdof+4]=-kmax*pow(10.,penalty_offset);
     371        Ke->values[4*numdof+0]=-kmax*pow(10.,penalty_offset);
     372        Ke->values[4*numdof+4]=+kmax*pow(10.,penalty_offset);
     373
     374        Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset);
     375        Ke->values[1*numdof+5]=-kmax*pow(10.,penalty_offset);
     376        Ke->values[5*numdof+1]=-kmax*pow(10.,penalty_offset);
     377        Ke->values[5*numdof+5]=+kmax*pow(10.,penalty_offset);
     378
     379        Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset);
     380        Ke->values[2*numdof+6]=-kmax*pow(10.,penalty_offset);
     381        Ke->values[6*numdof+2]=-kmax*pow(10.,penalty_offset);
     382        Ke->values[6*numdof+6]=+kmax*pow(10.,penalty_offset);
     383
     384        Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset);
     385        Ke->values[3*numdof+7]=-kmax*pow(10.,penalty_offset);
     386        Ke->values[7*numdof+3]=-kmax*pow(10.,penalty_offset);
     387        Ke->values[7*numdof+7]=+kmax*pow(10.,penalty_offset);
    388388
    389389        /*Clean up and return*/
    … …  
    404404
    405405        //Create elementary matrix: add penalty to
    406         Ke->values[0*numdof+0]=+kmax*pow((IssmDouble)10.0,penalty_factor);
    407         Ke->values[0*numdof+1]=-kmax*pow((IssmDouble)10.0,penalty_factor);
    408         Ke->values[1*numdof+0]=-kmax*pow((IssmDouble)10.0,penalty_factor);
    409         Ke->values[1*numdof+1]=+kmax*pow((IssmDouble)10.0,penalty_factor);
     406        Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_factor);
     407        Ke->values[0*numdof+1]=-kmax*pow(10.,penalty_factor);
     408        Ke->values[1*numdof+0]=-kmax*pow(10.,penalty_factor);
     409        Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_factor);
    410410
    411411        /*Clean up and return*/

  • issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp

    r15450 r15740  
    576576                if(shelf){
    577577                        /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */
    578                         pressure=rho_ice*gravity*pow(thickness,(IssmDouble)2)/(IssmDouble)2  - rho_water*gravity*pow(bed,(IssmDouble)2)/(IssmDouble)2;
     578                        pressure=rho_ice*gravity*pow(thickness,2)/2.- rho_water*gravity*pow(bed,2)/2.;
    579579                }
    580580                else{
    581581                        //We are on an icesheet, we assume the water column fills the entire front: */
    582                         pressure=rho_ice*gravity*pow(thickness,(IssmDouble)2)/(IssmDouble)2  - rho_water*gravity*pow(thickness,(IssmDouble)2)/(IssmDouble)2;
     582                        pressure=rho_ice*gravity*pow(thickness,2)/2.- rho_water*gravity*pow(thickness,2)/2.;
    583583                }
    584584        }
    585585        else if(fill==AirEnum){
    586                 pressure=rho_ice*gravity*pow(thickness,(IssmDouble)2)/(IssmDouble)2;   //icefront on an ice sheet, pressure imbalance ice vs air.
     586                pressure=rho_ice*gravity*pow(thickness,2)/2.;   //icefront on an ice sheet, pressure imbalance ice vs air.
    587587        }
    588588        else if(fill==IceEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
    … …  
    593593                if(!shelf) _error_("fill type " << fill << " not supported on ice sheets yet.");
    594594
    595                 pressure_litho=rho_ice*gravity*pow(thickness,(IssmDouble)2)/(IssmDouble)2;
     595                pressure_litho=rho_ice*gravity*pow(thickness,2)/2.;
    596596                pressure_air=0;
    597                 pressure_melange=rho_ice*gravity*pow(fraction*thickness,(IssmDouble)2)/(IssmDouble)2;
     597                pressure_melange=rho_ice*gravity*pow(fraction*thickness,2)/2.;
    598598                pressure_water=1.0/2.0*rho_water*gravity*  ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) );
    599599
    … …  
    684684                /*increase melange fraction: */
    685685                this->fraction+=fractionincrement;
    686                 if (this->fraction>1)this->fraction=(IssmDouble)1.0;
     686                if (this->fraction>1)this->fraction=1.;
    687687                //_printf_("riftfront " << this->Id() << " fraction: " << this->fraction << "\n");
    688688        }

  • issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp

    r15654 r15740  
    249249        else{
    250250                if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
    251                         viscosity=0.5*pow((IssmDouble)10,(IssmDouble)14);
     251                        viscosity=0.5*pow(10.,14);
    252252                }
    253253                else{
    … …  
    261261                        if(A==0){
    262262                                /*Maxiviscositym viscosity for 0 shear areas: */
    263                                 viscosity=2.5*pow(10.,17.);
     263                                viscosity=2.5*pow(10.,17);
    264264                        }
    265265                        else{
    … …  
    317317                if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
    318318                                (epsilon[3]==0) && (epsilon[4]==0)){
    319                         viscosity3d=0.5*pow((IssmDouble)10,(IssmDouble)14);
     319                        viscosity3d=0.5*pow(10.,14);
    320320                }
    321321                else{
    … …  
    332332                        if(A==0){
    333333                                /*Maxiviscosity3dm viscosity for 0 shear areas: */
    334                                 viscosity3d=2.25*pow((IssmDouble)10,(IssmDouble)17);
     334                                viscosity3d=2.25*pow(10.,17);
    335335                        }
    336336                        else{
    … …  
    378378
    379379        /*Get B and n*/
    380         eps0=pow((IssmDouble)10,(IssmDouble)-27);
     380        eps0=pow(10.,(IssmDouble)-27);
    381381        n=GetN();
    382382        Z=GetZ();
    … …  
    390390                if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
    391391                                (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){
    392                         viscosity3d=0.5*pow((IssmDouble)10,(IssmDouble)14);
     392                        viscosity3d=0.5*pow(10.,14);
    393393                }
    394394                else{
    … …  
    406406                        if(A==0){
    407407                                /*Maxiviscosity3dm viscosity for 0 shear areas: */
    408                                 viscosity3d=2.25*pow((IssmDouble)10,(IssmDouble)17);
     408                                viscosity3d=2.25*pow(10.,17);
    409409                        }
    410410                        else{
    … …  
    459459                if(A==0){
    460460                        /*Maximum viscosity_complement for 0 shear areas: */
    461                         viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17);
     461                        viscosity_complement=2.25*pow(10.,17);
    462462                }
    463463                else{
    … …  
    468468        }
    469469        else{
    470                 viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17);
     470                viscosity_complement=4.5*pow(10.,17);
    471471        }
    472472
    … …  
    515515                if(A==0){
    516516                        /*Maximum viscosity_complement for 0 shear areas: */
    517                         viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17);
     517                        viscosity_complement=2.25*pow(10.,17);
    518518                }
    519519                else{
    … …  
    524524        }
    525525        else{
    526                 viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17);
     526                viscosity_complement=4.5*pow(10.,17);
    527527        }
    528528
    … …  
    552552        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
    553553                                (epsilon[3]==0) && (epsilon[4]==0)){
    554                 mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14);
     554                mu_prime=0.5*pow(10.,14);
    555555        }
    556556        else{
    … …  
    586586        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
    587587                                (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){
    588                 mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14);
     588                mu_prime=0.5*pow(10.,14);
    589589        }
    590590        else{
    … …  
    620620
    621621        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
    622                 mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14);
     622                mu_prime=0.5*pow(10.,14);
    623623        }
    624624        else{
Note: See TracChangeset for help on using the changeset viewer.