Changeset 13081


Ignore:
Timestamp:
08/17/12 13:42:31 (13 years ago)
Author:
Mathieu Morlighem
Message:

CHG: cleaning up Tria (pow can use integers as exponents)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp

    r13073 r13081  
    370370                if(stabilization==2){
    371371                        /*Streamline upwinding*/
    372                         vel=sqrt(pow(vx,2.)+pow(vy,2.))+1.e-8;
     372                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
    373373                        K[0][0]=h/(2*vel)*vx*vx;
    374374                        K[1][0]=h/(2*vel)*vy*vx;
    … …  
    23282328IssmDouble Tria::SurfaceArea(void){
    23292329
    2330         int    i;
    23312330        IssmDouble S;
    23322331        IssmDouble normal[3];
    … …  
    23392338        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    23402339
    2341         for (i=0;i<3;i++){
     2340        for(int i=0;i<3;i++){
    23422341                v13[i]=xyz_list[0][i]-xyz_list[2][i];
    23432342                v23[i]=xyz_list[1][i]-xyz_list[2][i];
    … …  
    23482347        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
    23492348
    2350         S = 0.5 * sqrt(pow(normal[0],(IssmDouble)2)+pow(normal[1],(IssmDouble)2)+pow(normal[2],(IssmDouble)2));
     2349        S = 0.5 * sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
    23512350
    23522351        /*Return: */
    … …  
    23572356void Tria::SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]){
    23582357
    2359         int i;
    23602358        IssmDouble v13[3],v23[3];
    23612359        IssmDouble normal[3];
    23622360        IssmDouble normal_norm;
    23632361
    2364         for (i=0;i<3;i++){
     2362        for(int i=0;i<3;i++){
    23652363                v13[i]=xyz_list[0][i]-xyz_list[2][i];
    23662364                v23[i]=xyz_list[1][i]-xyz_list[2][i];
    … …  
    23712369        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
    23722370
    2373         normal_norm=sqrt( pow(normal[0],(IssmDouble)2)+pow(normal[1],(IssmDouble)2)+pow(normal[2],(IssmDouble)2) );
    2374 
    2375         *(surface_normal)=normal[0]/normal_norm;
    2376         *(surface_normal+1)=normal[1]/normal_norm;
    2377         *(surface_normal+2)=normal[2]/normal_norm;
     2371        normal_norm=sqrt( normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
     2372
     2373        *(surface_normal+0) = normal[0]/normal_norm;
     2374        *(surface_normal+1) = normal[1]/normal_norm;
     2375        *(surface_normal+2) = normal[2]/normal_norm;
    23782376}
    23792377/*}}}*/
    … …  
    26082606        normal[1]=sin(atan2(x1-x2,y2-y1));
    26092607
    2610         length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
     2608        length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
    26112609
    26122610        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    … …  
    29742972                //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
    29752973                surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    2976                 slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
     2974                slope_magnitude=sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
    29772975                if(slope_magnitude>MAXSLOPE) alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT);
    29782976                else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
    … …  
    33043302        /*Get Vz and compute vel*/
    33053303        GetInputListOnVertices(&vz[0],VzEnum,0);
    3306         for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     3304        for(i=0;i<NUMVERTICES;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    33073305
    33083306        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
    … …  
    33643362        /*Now Compute vel*/
    33653363        GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
    3366         for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     3364        for(i=0;i<NUMVERTICES;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    33673365
    33683366        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
    … …  
    39833981                 *      S                obs            obs
    39843982                 */
    3985                 misfit=1/S*pow( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) ,0.5);
     3983                misfit=1/S*sqrt( pow(vx-vxobs,2) + pow(vy-vyobs,2));
    39863984
    39873985                if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum);
    … …  
    40464044                 *                            obs
    40474045                 */
    4048                 velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
    4049                 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
    4050                 misfit=4*pow(meanvel,2.)*pow(log(velocity_mag/obs_velocity_mag),2.);
     4046                velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
     4047                obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
     4048                misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
    40514049
    40524050                if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVelMisfitEnum);
    … …  
    41124110                 *                              obs                       obs
    41134111                 */
    4114                 misfit=0.5*pow(meanvel,2.)*(
    4115                                         pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2.) +
    4116                                         pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2.) );
     4112                misfit=0.5*pow(meanvel,2)*(
     4113                                        pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2) +
     4114                                        pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2) );
    41174115
    41184116                if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVxVyMisfitEnum);
    … …  
    41754173                 *
    41764174                 */
    4177                 misfit=0.5*( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) );
     4175                misfit=0.5*( pow(vx-vxobs,2) + pow(vy-vyobs,2) );
    41784176
    41794177                if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum);
    … …  
    42384236                 *              obs                        obs                     
    42394237                 */
    4240                 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
    4241                 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
    4242                 misfit=0.5*(scalex*pow((vx-vxobs),2.)+scaley*pow((vy-vyobs),2.));
     4238                scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
     4239                scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
     4240                misfit=0.5*(scalex*pow((vx-vxobs),2)+scaley*pow((vy-vyobs),2));
    42434241                if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceRelVelMisfitEnum);
    42444242
    … …  
    42884286
    42894287                /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
    4290                 Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
     4288                Jelem+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
    42914289        }
    42924290
    … …  
    44374435
    44384436                /*compute ThicknessAbsMisfit*/
    4439                 Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight;
     4437                Jelem+=0.5*(thickness-thicknessobs)*(thickness-thicknessobs)*weight*Jdet*gauss->weight;
    44404438        }
    44414439
    … …  
    46254623                                         */
    46264624                                        for (i=0;i<NUMVERTICES;i++){
    4627                                                 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
    4628                                                 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
     4625                                                scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
     4626                                                scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
    46294627                                                dux=scalex*(vxobs-vx);
    46304628                                                duy=scaley*(vyobs-vy);
    … …  
    46464644                                         */
    46474645                                        for (i=0;i<NUMVERTICES;i++){
    4648                                                 velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
    4649                                                 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
    4650                                                 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
     4646                                                velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
     4647                                                obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
     4648                                                scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
    46514649                                                dux=scale*vx;
    46524650                                                duy=scale*vy;
    … …  
    46664664                                         */
    46674665                                        for (i=0;i<NUMVERTICES;i++){
    4668                                                 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
     4666                                                scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
    46694667                                                dux=scale*(vxobs-vx);
    46704668                                                duy=scale*(vyobs-vy);
    … …  
    46844682                                         */
    46854683                                        for (i=0;i<NUMVERTICES;i++){
    4686                                                 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
    4687                                                 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
     4684                                                dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
     4685                                                duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
    46884686                                                pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
    46894687                                                pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    … …  
    48084806                                         */
    48094807                                        for (i=0;i<NUMVERTICES;i++){
    4810                                                 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
    4811                                                 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
     4808                                                scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
     4809                                                scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
    48124810                                                dux=scalex*(vxobs-vx);
    48134811                                                duy=scaley*(vyobs-vy);
    … …  
    48294827                                         */
    48304828                                        for (i=0;i<NUMVERTICES;i++){
    4831                                                 velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
    4832                                                 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
    4833                                                 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
     4829                                                velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
     4830                                                obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
     4831                                                scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
    48344832                                                dux=scale*vx;
    48354833                                                duy=scale*vy;
    … …  
    48494847                                         */
    48504848                                        for (i=0;i<NUMVERTICES;i++){
    4851                                                 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
     4849                                                scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
    48524850                                                dux=scale*(vxobs-vx);
    48534851                                                duy=scale*(vyobs-vy);
    … …  
    48674865                                         */
    48684866                                        for (i=0;i<NUMVERTICES;i++){
    4869                                                 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
    4870                                                 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
     4867                                                dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
     4868                                                duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
    48714869                                                pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
    48724870                                                pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    … …  
    49364934
    49374935                /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
    4938                 Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
     4936                Jelem+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
    49394937        }
    49404938
    … …  
    51935191
    51945192        /* compute VelocityFactor */
    5195         VelocityFactor= n_man*pow(CR,2)*rho_water*g/mu_water;
     5193        VelocityFactor= n_man*CR*CR*rho_water*g/mu_water;
    51965194       
    51975195        gauss=new GaussTria();
    … …  
    52055203
    52065204                /* Water velocity x and y components */
    5207         //      vx[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
    5208         //      vy[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
    5209        
    5210                 vx[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
    5211                 vy[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
     5205        //      vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
     5206        //      vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
     5207                vx[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
     5208                vy[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
    52125209        }
    52135210
    … …  
    53045301
    53055302                /*Artificial diffusivity*/
    5306                 vel=sqrt(pow(vx,2.)+pow(vy,2.));
     5303                vel=sqrt(vx*vx+vy*vy);
    53075304                K[0][0]=diffusivity*h/(2*vel)*vx*vx;
    53085305                K[1][0]=diffusivity*h/(2*vel)*vy*vx;
    … …  
    54145411
    54155412        /*Intermediaries*/
    5416         const int numdof = NDOF1*NUMVERTICES;
    5417 
    5418         int       i;
    5419         int*      doflist=NULL;
    5420         IssmDouble    values[numdof];
     5413        const int   numdof         = NDOF1 *NUMVERTICES;
     5414        int        *doflist        = NULL;
     5415        IssmDouble  values[numdof];
    54215416
    54225417        /*Get dof list: */
    … …  
    54245419
    54255420        /*Use the dof list to index into the solution vector: */
    5426         for(i=0;i<numdof;i++){
     5421        for(int i=0;i<numdof;i++){
    54275422                values[i]=solution[doflist[i]];
    54285423                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    5429                 if (values[i]<pow((IssmDouble)10,(IssmDouble)-10))values[i]=pow((IssmDouble)10,(IssmDouble)-10); //correcting the water column to positive values
    5430  
     5424                if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values
    54315425        }
    54325426
    … …  
    56795673                if(stabilization==1){
    56805674                        /*Streamline upwinding*/
    5681                         vel=sqrt(pow(vx,2.)+pow(vy,2.));
     5675                        vel=sqrt(vx*vx+vy*vy);
    56825676                        K[0][0]=h/(2*vel)*vx*vx;
    56835677                        K[1][0]=h/(2*vel)*vy*vx;
Note: See TracChangeset for help on using the changeset viewer.