Changeset 25905


Ignore:
Timestamp:
12/20/20 13:01:58 (4 years ago)
Author:
Mathieu Morlighem
Message:

CHG: removing some xNew to speed up code

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r25807 r25905  
    651651        const int NUM_VERTICES = this->GetNumberOfVertices();
    652652
    653         int        i;
    654         IssmDouble accuref, runoffref; //reference values at given altitude
    655         IssmDouble accualti, runoffalti; //reference altitudes
     653        IssmDouble accuref, runoffref;  //reference values at given altitude
    656654        IssmDouble accugrad, runoffgrad; //gradients from reference altitude
    657         IssmDouble rho_water, rho_ice;
    658         IssmDouble time;
    659 
    660         IssmDouble*             smb      = xNew<IssmDouble>(NUM_VERTICES);
    661         IssmDouble*             surf     = xNew<IssmDouble>(NUM_VERTICES);
    662         IssmDouble*             accu     = xNew<IssmDouble>(NUM_VERTICES);
    663         IssmDouble*             runoff = xNew<IssmDouble>(NUM_VERTICES);
     655
     656        IssmDouble smb[MAXVERTICES];
     657        IssmDouble surf[MAXVERTICES];
     658        IssmDouble accu[MAXVERTICES];
     659        IssmDouble runoff[MAXVERTICES];
    664660
    665661        /*Get material parameters :*/
    666         rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
    667         rho_ice=this->FindParam(MaterialsRhoIceEnum);
     662        IssmDouble rho_water = this->FindParam(MaterialsRhoSeawaterEnum);
     663        IssmDouble rho_ice   = this->FindParam(MaterialsRhoIceEnum);
    668664
    669665        /*Recover parameters*/
    670         parameters->FindParam(&time,TimeEnum);
    671         parameters->FindParam(&accualti,SmbAccualtiEnum);
    672         parameters->FindParam(&runoffalti,SmbRunoffaltiEnum);
     666        IssmDouble time       = this->FindParam(TimeEnum);
     667        IssmDouble accualti   = this->FindParam(SmbAccualtiEnum);
     668        IssmDouble runoffalti = this->FindParam(SmbRunoffaltiEnum);
    673669
    674670        /*Recover reference values at current time*/
    … …  
    683679        /*Compute the temperature and precipitation*/
    684680        for(int iv=0;iv<NUM_VERTICES;iv++){
    685                 accu[iv]=max(0.,(accuref+(surf[iv]-accualti)*accugrad));
    686                 runoff[iv]=max(0.,(runoffref+(surf[iv]-runoffalti)*runoffgrad));
    687                 smb[iv]=(accu[iv]-runoff[iv])*rho_ice/rho_water;
     681                accu[iv]   = max(0.,(accuref+(surf[iv]-accualti)*accugrad));
     682                runoff[iv] = max(0.,(runoffref+(surf[iv]-runoffalti)*runoffgrad));
     683                smb[iv]    = (accu[iv]-runoff[iv])*rho_ice/rho_water;
    688684        }
    689685
    … …  
    701697        default: _error_("Not implemented yet");
    702698        }
    703         /*clean-up*/
    704         xDelete<IssmDouble>(surf);
    705         xDelete<IssmDouble>(accu);
    706         xDelete<IssmDouble>(runoff);
    707         xDelete<IssmDouble>(smb);
    708699}
    709700/*}}}*/
    … …  
    12931284}
    12941285/*}}}*/
    1295 /* void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){/\*{{{*\/ */
    1296 
    1297 /*      /\*Fetch number vertices for this element and allocate arrays*\/ */
    1298 /*      const int NUM_VERTICES = this->GetNumberOfVertices(); */
    1299 /*      int*        vertexpidlist = xNew<int>(NUM_VERTICES); */
    1300 /*      IssmDouble* values        = xNew<IssmDouble>(NUM_VERTICES); */
    1301 
    1302 /*      /\*Fill in values*\/ */
    1303 /*      this->GetVerticesPidList(vertexpidlist); */
    1304 /*      this->GetInputListOnVertices(values,input_enum); */
    1305 /*      vector->SetValues(NUM_VERTICES,vertexpidlist,values,INS_VAL); */
    1306 
    1307 /*      /\*Clean up*\/ */
    1308 /*      xDelete<int>(vertexpidlist); */
    1309 /*      xDelete<IssmDouble>(values); */
    1310 
    1311 /* } */
    1312 /* /\*}}}*\/ */
    13131286void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum,int type){/*{{{*/
    13141287
    1315         /*Fetch number vertices for this element and allocate arrays*/
    1316         const int NUM_VERTICES = this->GetNumberOfVertices();
    1317 
    1318         int         numnodes    = this->GetNumberOfNodes();
    1319         int*        doflist     = NULL;
    1320         IssmDouble  value;
    1321         IssmDouble* values      = NULL;
    1322         Input*     input       = NULL;
    13231288
    13241289        switch(type){
    1325                 case ElementSIdEnum:
    1326                         input=this->GetInput(input_enum); _assert_(input);
     1290                case ElementSIdEnum:{
     1291         IssmDouble  value;
     1292                        Input* input=this->GetInput(input_enum); _assert_(input);
    13271293                        input->GetInputAverage(&value);
    13281294                        vector->SetValue(this->sid,value,INS_VAL);
     1295                          }
    13291296                        break;
    1330                 case VertexPIdEnum:
    1331                         doflist = xNew<int>(NUM_VERTICES);
    1332                         values = xNew<IssmDouble>(NUM_VERTICES);
     1297                case VertexPIdEnum:{
     1298         int        doflist[MAXVERTICES];
     1299         IssmDouble values[MAXVERTICES];
     1300                        const int  NUM_VERTICES = this->GetNumberOfVertices();
    13331301                        /*Fill in values*/
    1334                         this->GetVerticesPidList(doflist);
    1335                         this->GetInputListOnVertices(values,input_enum);
     1302                        this->GetVerticesPidList(&doflist[0]);
     1303                        this->GetInputListOnVertices(&values[0],input_enum);
    13361304                        vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
     1305                         }
    13371306                        break;
    1338                 case VertexSIdEnum:
    1339                         doflist = xNew<int>(NUM_VERTICES);
    1340                         values = xNew<IssmDouble>(NUM_VERTICES);
     1307                case VertexSIdEnum:{
     1308         int        doflist[MAXVERTICES];
     1309         IssmDouble values[MAXVERTICES];
     1310                        const int  NUM_VERTICES = this->GetNumberOfVertices();
    13411311                        /*Fill in values*/
    13421312                        this->GetVerticesSidList(doflist);
    13431313                        this->GetInputListOnVertices(values,input_enum);
    13441314                        vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
     1315                         }
    13451316                        break;
    1346                 case NodesEnum:
    1347                         doflist = xNew<int>(numnodes);
    1348                         values = xNew<IssmDouble>(numnodes);
     1317                case NodesEnum:{
     1318         int         numnodes = this->GetNumberOfNodes();
     1319         int        *doflist  = xNew<int>(numnodes);
     1320         IssmDouble *values   = xNew<IssmDouble>(numnodes);
    13491321                        /*Fill in values*/
    13501322                        this->GetInputListOnNodes(values,input_enum);
    13511323                        this->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    13521324                        vector->SetValues(numnodes,doflist,values,INS_VAL);
     1325         xDelete<int>(doflist);
     1326         xDelete<IssmDouble>(values);
     1327                     }
    13531328                        break;
    1354                 case NodeSIdEnum:
    1355                         doflist = xNew<int>(numnodes);
    1356                         values = xNew<IssmDouble>(numnodes);
     1329                case NodeSIdEnum:{
     1330         int         numnodes = this->GetNumberOfNodes();
     1331         int        *doflist  = xNew<int>(numnodes);
     1332         IssmDouble *values   = xNew<IssmDouble>(numnodes);
    13571333                        /*Fill in values*/
    13581334                        this->GetNodesSidList(doflist);
    13591335                        this->GetInputListOnNodes(values,input_enum);
    13601336                        vector->SetValues(numnodes,doflist,values,INS_VAL);
     1337         xDelete<int>(doflist);
     1338         xDelete<IssmDouble>(values);
     1339                       }
    13611340                        break;
    13621341                default:
    … …  
    13641343        }
    13651344
    1366         /*Clean up*/
    1367         xDelete<int>(doflist);
    1368         xDelete<IssmDouble>(values);
    1369 
    1370 }
    1371 /*}}}*/
     1345}/*}}}*/
    13721346void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum,int type, IssmDouble time){/*{{{*/
    13731347
    1374         /*Fetch number vertices for this element and allocate arrays*/
    1375         const int NUM_VERTICES = this->GetNumberOfVertices();
    1376 
    1377         int         numnodes    = this->GetNumberOfNodes();
    1378         int*        doflist     = NULL;
    1379         IssmDouble* values      = NULL;
    1380 
    13811348        switch(type){
    1382                 case VertexPIdEnum:
    1383                         doflist = xNew<int>(NUM_VERTICES);
    1384                         values = xNew<IssmDouble>(NUM_VERTICES);
     1349                case VertexPIdEnum:{
     1350                        int        doflist[MAXVERTICES];
     1351                        IssmDouble values[MAXVERTICES];
     1352                        const int  NUM_VERTICES = this->GetNumberOfVertices();
    13851353                        /*Fill in values*/
    13861354                        this->GetVerticesPidList(doflist);
    13871355                        this->GetInputListOnVerticesAtTime(values,input_enum,time);
    1388                         vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
     1356                        vector->SetValues(NUM_VERTICES,doflist,&values[0],INS_VAL);
     1357                                                                 }
    13891358                        break;
    1390                 case VertexSIdEnum:
    1391                         doflist = xNew<int>(NUM_VERTICES);
    1392                         values = xNew<IssmDouble>(NUM_VERTICES);
     1359                case VertexSIdEnum:{
     1360                        int        doflist[MAXVERTICES];
     1361                        IssmDouble values[MAXVERTICES];
     1362                        const int  NUM_VERTICES = this->GetNumberOfVertices();
    13931363                        /*Fill in values*/
    13941364                        this->GetVerticesSidList(doflist);
    13951365                        this->GetInputListOnVerticesAtTime(values,input_enum,time);
    1396                         vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
     1366                        vector->SetValues(NUM_VERTICES,doflist,&values[0],INS_VAL);
     1367                                                                 }
    13971368                        break;
    13981369                default:
    13991370                        _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
    14001371        }
    1401 
    1402         /*Clean up*/
    1403         xDelete<int>(doflist);
    1404         xDelete<IssmDouble>(values);
    1405 
    14061372}
    14071373/*}}}*/
    … …  
    14491415        /*Create list of x*/
    14501416        const int NUM_VERTICES = this->GetNumberOfVertices();
    1451 
    1452         IssmDouble* x_list = xNew<IssmDouble>(NUM_VERTICES);
    1453 
     1417        IssmDouble x_list[MAXVERTICES];
    14541418        for(int i=0;i<NUM_VERTICES;i++) x_list[i]=xyz_list[i*3+0];
    1455         ValueP1OnGauss(&x,x_list,gauss);
    1456 
    1457         xDelete<IssmDouble>(x_list);
     1419
     1420        /*Get value at gauss point*/
     1421        ValueP1OnGauss(&x,&x_list[0],gauss);
    14581422        return x;
    14591423}/*}}}*/
    … …  
    14651429        /*Create list of y*/
    14661430        const int NUM_VERTICES = this->GetNumberOfVertices();
    1467 
    1468         IssmDouble* y_list      = xNew<IssmDouble>(NUM_VERTICES);
    1469 
     1431        IssmDouble y_list[MAXVERTICES];
    14701432        for(int i=0;i<NUM_VERTICES;i++) y_list[i]=xyz_list[i*3+1];
    1471         ValueP1OnGauss(&y,y_list,gauss);
    1472 
    1473         xDelete<IssmDouble>(y_list);
     1433
     1434        /*Get value at gauss point*/
     1435        ValueP1OnGauss(&y,&y_list[0],gauss);
    14741436        return y;
    14751437}/*}}}*/
    … …  
    14811443        /*Create list of z*/
    14821444        const int NUM_VERTICES = this->GetNumberOfVertices();
    1483 
    1484         IssmDouble* z_list      = xNew<IssmDouble>(NUM_VERTICES);
    1485 
     1445        IssmDouble z_list[MAXVERTICES];
    14861446        for(int i=0;i<NUM_VERTICES;i++) z_list[i]=xyz_list[i*3+2];
    1487         ValueP1OnGauss(&z,z_list,gauss);
    1488 
    1489         xDelete<IssmDouble>(z_list);
     1447
     1448        /*Get value at gauss point*/
     1449        ValueP1OnGauss(&z,&z_list[0],gauss);
    14901450        return z;
    14911451}/*}}}*/
    … …  
    16121572void       Element::InputCreate(IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){/*{{{*/
    16131573
    1614         /*Intermediaries*/
    1615         int i,t;
    1616 
    16171574        /*Branch on type of vector: nodal or elementary: */
    16181575        if(vector_type==1){ //nodal vector
    … …  
    16201577                const int NUM_VERTICES = this->GetNumberOfVertices();
    16211578
    1622                 int        *vertexids  = xNew<int>(NUM_VERTICES);
    1623                 int        *vertexlids = xNew<int>(NUM_VERTICES);
    1624                 IssmDouble *values     = xNew<IssmDouble>(NUM_VERTICES);
     1579                int        vertexids[MAXVERTICES];
     1580                int        vertexlids[MAXVERTICES];
     1581                IssmDouble values[MAXVERTICES];
    16251582
    16261583                /*Recover vertices ids needed to initialize inputs*/
    16271584                _assert_(iomodel->elements);
    1628                 for(i=0;i<NUM_VERTICES;i++){
     1585                for(int i=0;i<NUM_VERTICES;i++){
    16291586                        vertexids[i] =reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
    16301587                        vertexlids[i]=iomodel->my_vertices_lids[vertexids[i]-1];
    … …  
    16371594                }
    16381595                else if(M==iomodel->numberofvertices){
    1639                         for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
     1596                        for(int i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
    16401597                        this->SetElementInput(inputs,NUM_VERTICES,vertexlids,values,vector_enum);
    16411598                }
    … …  
    16431600                        /*create transient input: */
    16441601                        IssmDouble* times = xNew<IssmDouble>(N);
    1645                         for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1602                        for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    16461603                        inputs->SetTransientInput(vector_enum,times,N);
    16471604                        TransientInput* transientinput = inputs->GetTransientInput(vector_enum);
    1648                         for(t=0;t<N;t++){
    1649                                 for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
     1605                        for(int t=0;t<N;t++){
     1606                                for(int i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
    16501607                                switch(this->ObjectEnum()){
    1651                                         case TriaEnum:  transientinput->AddTriaTimeInput( t,NUM_VERTICES,vertexlids,values,P1Enum); break;
    1652                                         case PentaEnum: transientinput->AddPentaTimeInput(t,NUM_VERTICES,vertexlids,values,P1Enum); break;
     1608                                        case TriaEnum:  transientinput->AddTriaTimeInput( t,NUM_VERTICES,vertexlids,&values[0],P1Enum); break;
     1609                                        case PentaEnum: transientinput->AddPentaTimeInput(t,NUM_VERTICES,vertexlids,&values[0],P1Enum); break;
    16531610                                        default: _error_("Not implemented yet");
    16541611                                }
    … …  
    16591616
    16601617                        /*This is a Patch!*/
    1661                         xDelete<IssmDouble>(values);
    1662                         values = xNew<IssmDouble>(N);
    1663                         for(int j=0;j<N;j++) values[j]=vector[this->Sid()*N+j];
     1618                        IssmDouble* evalues = xNew<IssmDouble>(N);
     1619                        for(int j=0;j<N;j++) evalues[j]=vector[this->Sid()*N+j];
    16641620
    16651621                        if (N==this->GetNumberOfNodes(P1Enum)){
    1666                                 this->SetElementInput(inputs,NUM_VERTICES,vertexlids,values,vector_enum);
     1622                                this->SetElementInput(inputs,NUM_VERTICES,vertexlids,evalues,vector_enum);
    16671623                        }
    16681624                        else if(N==this->GetNumberOfNodes(P0Enum)){
    1669                                 this->SetElementInput(inputs,vector_enum,values[0]);
     1625                                this->SetElementInput(inputs,vector_enum,evalues[0]);
    16701626                        }
    16711627                        else if(N==this->GetNumberOfNodes(P1xP2Enum)){ _assert_(this->ObjectEnum()==PentaEnum);
    1672                                 inputs->SetPentaInput(vector_enum,P1xP2Enum,this->lid,N,values);
     1628                                inputs->SetPentaInput(vector_enum,P1xP2Enum,this->lid,N,evalues);
    16731629                        }
    16741630                        else if(N==this->GetNumberOfNodes(P1xP3Enum)){ _assert_(this->ObjectEnum()==PentaEnum);
    1675                                 inputs->SetPentaInput(vector_enum,P1xP3Enum,this->lid,N,values);
     1631                                inputs->SetPentaInput(vector_enum,P1xP3Enum,this->lid,N,evalues);
    16761632                        }
    16771633                        else{
    16781634                                _error_("Patch interpolation not supported yet");
    16791635                        }
     1636                        xDelete<IssmDouble>(evalues);
    16801637
    16811638                }
    … …  
    16831640                        _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
    16841641                }
    1685 
    1686                 xDelete<IssmDouble>(values);
    1687                 xDelete<int>(vertexids);
    1688                 xDelete<int>(vertexlids);
    16891642        }
    16901643        else if(vector_type==2){ //element vector
    … …  
    17081661                        /*create transient input: */
    17091662                        IssmDouble* times = xNew<IssmDouble>(N);
    1710                         for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1663                        for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    17111664                        inputs->SetTransientInput(vector_enum,times,N);
    17121665                        TransientInput* transientinput = inputs->GetTransientInput(vector_enum);
    1713                         for(t=0;t<N;t++){
     1666                        for(int t=0;t<N;t++){
    17141667                                value=vector[N*this->Sid()+t];
    17151668                                switch(this->ObjectEnum()){
    … …  
    17281681                if(M==iomodel->numberofelements){
    17291682                        IssmDouble* layers = xNewZeroInit<IssmDouble>(N);
    1730                         for(t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
     1683                        for(int t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
    17311684                        inputs->SetArrayInput(vector_enum,this->lid,layers,N);
    17321685                        xDelete<IssmDouble>(layers);
    … …  
    18501803         */
    18511804
    1852         /*Intermediaries*/
    1853         int i,t;
    1854 
    18551805        /*Branch on type of vector: nodal or elementary: */
    18561806        if(vector_type==1){ //nodal vector
    … …  
    18581808                const int NUM_VERTICES = this->GetNumberOfVertices();
    18591809
    1860                 int        *vertexids  = xNew<int>(NUM_VERTICES);
    1861                 int        *vertexlids = xNew<int>(NUM_VERTICES);
    1862                 IssmDouble *values     = xNew<IssmDouble>(NUM_VERTICES);
     1810                int        vertexids[MAXVERTICES];
     1811                int        vertexlids[MAXVERTICES];
     1812                IssmDouble values[MAXVERTICES];
    18631813
    18641814                /*Recover vertices ids needed to initialize inputs*/
    18651815                _assert_(iomodel->elements);
    1866                 for(i=0;i<NUM_VERTICES;i++){
     1816                for(int i=0;i<NUM_VERTICES;i++){
    18671817                        vertexids[i] =reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
    18681818                        vertexlids[i]=iomodel->my_vertices_lids[vertexids[i]-1];
    … …  
    18761826                }
    18771827                else if(M==iomodel->numberofvertices){
    1878                         for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
     1828                        for(int i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
    18791829                        switch(this->ObjectEnum()){
    18801830                                case TriaEnum:  inputs->SetTriaDatasetInput(enum_type,input_id,P1Enum,NUM_VERTICES,vertexlids,values); break;
    … …  
    18861836                        /*create transient input: */
    18871837                        IssmDouble* times = xNew<IssmDouble>(N);
    1888                         for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1838                        for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    18891839                        TransientInput* transientinput = inputs->SetDatasetTransientInput(enum_type,input_id,times,N);
    1890                         for(t=0;t<N;t++){
    1891                                 for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
     1840                        for(int t=0;t<N;t++){
     1841                                for(int i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
    18921842                                switch(this->ObjectEnum()){
    18931843                                        case TriaEnum:  transientinput->AddTriaTimeInput( t,NUM_VERTICES,vertexlids,values,P1Enum); break;
    … …  
    19011851                        _error_("not implemented yet (M="<<M<<")");
    19021852                }
    1903 
    1904                 xDelete<IssmDouble>(values);
    1905                 xDelete<int>(vertexids);
    1906                 xDelete<int>(vertexlids);
    19071853        }
    19081854        else if(vector_type==2){ //element vector
    … …  
    20652011        /*Allocate some arrays*/
    20662012        const int numvertices = this->GetNumberOfVertices();
    2067         IssmDouble* basalmeltrate = xNew<IssmDouble>(numvertices);
     2013        IssmDouble basalmeltrate[MAXVERTICES];
    20682014
    20692015        /*Get variables*/
    … …  
    21142060        xDelete<IssmDouble>(mean_tf);
    21152061        xDelete<IssmDouble>(depths);
    2116         xDelete<IssmDouble>(basalmeltrate);
    21172062
    21182063}/*}}}*/
    … …  
    21212066        const int NUM_VERTICES = this->GetNumberOfVertices();
    21222067
    2123         IssmDouble  deepwaterel,upperwaterel,deepwatermelt,upperwatermelt;
    2124         IssmDouble *base         = xNew<IssmDouble>(NUM_VERTICES);
    2125         IssmDouble *values       = xNew<IssmDouble>(NUM_VERTICES);
    2126         IssmDouble *perturbation = xNew<IssmDouble>(NUM_VERTICES);
    2127         IssmDouble      time;
     2068        IssmDouble deepwaterel,upperwaterel,deepwatermelt,upperwatermelt;
     2069        IssmDouble base[MAXVERTICES];
     2070        IssmDouble values[MAXVERTICES];
     2071        IssmDouble perturbation[MAXVERTICES];
     2072        IssmDouble time;
    21282073
    21292074        parameters->FindParam(&time,TimeEnum);
    … …  
    21342079        _assert_(upperwaterel>deepwaterel);
    21352080
    2136         this->GetInputListOnVertices(base,BaseEnum);
    2137    this->GetInputListOnVertices(perturbation,BasalforcingsPerturbationMeltingRateEnum);
     2081        this->GetInputListOnVertices(&base[0],BaseEnum);
     2082   this->GetInputListOnVertices(&perturbation[0],BasalforcingsPerturbationMeltingRateEnum);
    21382083        for(int i=0;i<NUM_VERTICES;i++){
    21392084                if(base[i]>=upperwaterel){
    … …  
    21522097
    21532098        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
    2154         xDelete<IssmDouble>(base);
    2155    xDelete<IssmDouble>(perturbation);
    2156         xDelete<IssmDouble>(values);
    21572099
    21582100}/*}}}*/
    … …  
    21612103        const int NUM_VERTICES = this->GetNumberOfVertices();
    21622104
    2163         IssmDouble *deepwatermelt = xNew<IssmDouble>(NUM_VERTICES);
    2164         IssmDouble *deepwaterel   = xNew<IssmDouble>(NUM_VERTICES);
    2165         IssmDouble *upperwaterel  = xNew<IssmDouble>(NUM_VERTICES);
    2166         IssmDouble *base          = xNew<IssmDouble>(NUM_VERTICES);
    2167         IssmDouble *values        = xNew<IssmDouble>(NUM_VERTICES);
    2168 
    2169         this->GetInputListOnVertices(base,BaseEnum);
    2170         this->GetInputListOnVertices(deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum);
    2171         this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum);
    2172         this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum);
     2105        IssmDouble deepwatermelt[MAXVERTICES];
     2106        IssmDouble deepwaterel[MAXVERTICES];;
     2107        IssmDouble upperwaterel[MAXVERTICES];
     2108        IssmDouble base[MAXVERTICES];
     2109        IssmDouble values[MAXVERTICES];
     2110
     2111        this->GetInputListOnVertices(&base[0],BaseEnum);
     2112        this->GetInputListOnVertices(&deepwatermelt[0],BasalforcingsDeepwaterMeltingRateEnum);
     2113        this->GetInputListOnVertices(&deepwaterel[0],BasalforcingsDeepwaterElevationEnum);
     2114        this->GetInputListOnVertices(&upperwaterel[0],BasalforcingsUpperwaterElevationEnum);
    21732115
    21742116        for(int i=0;i<NUM_VERTICES;i++){
    … …  
    21782120        }
    21792121
    2180         this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
    2181         xDelete<IssmDouble>(base);
    2182         xDelete<IssmDouble>(deepwaterel);
    2183         xDelete<IssmDouble>(deepwatermelt);
    2184         xDelete<IssmDouble>(upperwaterel);
    2185         xDelete<IssmDouble>(values);
    2186 
     2122        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,&values[0],P1Enum);
    21872123}/*}}}*/
    21882124void       Element::MantlePlumeGeothermalFlux(){/*{{{*/
    … …  
    21932129        IssmDouble  crustheat,plumeheat,dt,middleplumedepth,a,e,eprime,A0,lambda,Alambda,dAlambda;
    21942130        IssmDouble  x,y,z,c;
    2195         IssmDouble* values   = xNew<IssmDouble>(NUM_VERTICES);
     2131        IssmDouble  values[MAXVERTICES];
    21962132        IssmDouble *xyz_list = NULL;
    21972133
    … …  
    22282164        }
    22292165
    2230         this->AddInput(BasalforcingsGeothermalfluxEnum,values,P1Enum);
     2166        this->AddInput(BasalforcingsGeothermalfluxEnum,&values[0],P1Enum);
    22312167        xDelete<IssmDouble>(xyz_list);
    2232         xDelete<IssmDouble>(values);
    22332168
    22342169}/*}}}*/
    … …  
    22532188
    22542189        const int  NUM_VERTICES = this->GetNumberOfVertices();
    2255         int        i,migration_style;
     2190        int        migration_style;
    22562191        IssmDouble bed_hydro,yts;
    2257         IssmDouble rho_water,rho_ice,density;
    2258         IssmDouble* melting = xNew<IssmDouble>(NUM_VERTICES);
    2259         IssmDouble* phi     = xNew<IssmDouble>(NUM_VERTICES);
    2260         IssmDouble* h       = xNew<IssmDouble>(NUM_VERTICES);
    2261         IssmDouble* s       = xNew<IssmDouble>(NUM_VERTICES);
    2262         IssmDouble* b       = xNew<IssmDouble>(NUM_VERTICES);
    2263         IssmDouble* r       = xNew<IssmDouble>(NUM_VERTICES);
    2264         IssmDouble* sl      = xNew<IssmDouble>(NUM_VERTICES);
     2192        IssmDouble melting[MAXVERTICES];
     2193        IssmDouble phi[MAXVERTICES];
     2194        IssmDouble h[MAXVERTICES];
     2195        IssmDouble s[MAXVERTICES];
     2196        IssmDouble b[MAXVERTICES];
     2197        IssmDouble r[MAXVERTICES];
     2198        IssmDouble sl[MAXVERTICES];
    22652199
    22662200        /*Recover info at the vertices: */
    … …  
    22732207        GetInputListOnVertices(&sl[0],SealevelEnum);
    22742208        GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum);
    2275         rho_water   = FindParam(MaterialsRhoSeawaterEnum);
    2276         rho_ice     = FindParam(MaterialsRhoIceEnum);
    2277         density     = rho_ice/rho_water;
     2209        IssmDouble rho_water   = FindParam(MaterialsRhoSeawaterEnum);
     2210        IssmDouble rho_ice     = FindParam(MaterialsRhoIceEnum);
     2211        IssmDouble density     = rho_ice/rho_water;
    22782212
    22792213        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
    2280         for(i=0;i<NUM_VERTICES;i++){
     2214        for(int i=0;i<NUM_VERTICES;i++){
    22812215                /* Contact FS*/
    22822216                if(migration_style == ContactEnum){
    … …  
    22842218                        if(phi[i]>=0.) b[i]=r[i];
    22852219                }
    2286                 else if(migration_style == GroundingOnlyEnum && b[i]<r[i]) b[i]=r[i];
    2287                 /*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
     2220                else if(migration_style == GroundingOnlyEnum && b[i]<r[i]){
     2221                        /*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
     2222                        b[i]=r[i];
     2223                }
    22882224                else if(phi[i]<=0.){
    22892225                        if(b[i]<=r[i]){
    … …  
    22992235                                /*Unground only if the element is connected to the ice shelf*/
    23002236                                if(migration_style==AggressiveMigrationEnum || migration_style==SubelementMigrationEnum){
    2301                                         s[i]        = (1-density)*h[i]+sl[i];
    2302                                         b[i]        = -density*h[i]+sl[i];
     2237                                        s[i] = (1.-density)*h[i]+sl[i];
     2238                                        b[i] = -density*h[i]+sl[i];
    23032239                                }
    23042240                                else if(migration_style==SoftMigrationEnum && phi_ungrounding[vertices[i]->Pid()]<0.){
    2305                                         s[i]        = (1-density)*h[i]+sl[i];
    2306                                         b[i]        = -density*h[i]+sl[i];
     2241                                        s[i] = (1.-density)*h[i]+sl[i];
     2242                                        b[i] = -density*h[i]+sl[i];
    23072243                                }
    23082244                                else{
    … …  
    23142250
    23152251        /*Recalculate phi*/
    2316         for(i=0;i<NUM_VERTICES;i++){
     2252        for(int i=0;i<NUM_VERTICES;i++){
    23172253                if(migration_style==SoftMigrationEnum){
    23182254                        bed_hydro=-density*h[i]+sl[i];
    … …  
    23212257                        }
    23222258                }
    2323                 else if(migration_style!=ContactEnum) phi[i]=h[i]+(r[i]-sl[i])/density;
     2259                else if(migration_style!=ContactEnum){
     2260                        phi[i]=h[i]+(r[i]-sl[i])/density;
     2261                }
    23242262                else{
    23252263                        /*do nothing*/
    … …  
    23322270        this->AddInput(BaseEnum,&b[0],P1Enum);
    23332271
    2334         /*Delete*/
    2335         xDelete<IssmDouble>(melting);
    2336         xDelete<IssmDouble>(phi);
    2337         xDelete<IssmDouble>(r);
    2338         xDelete<IssmDouble>(b);
    2339         xDelete<IssmDouble>(s);
    2340         xDelete<IssmDouble>(sl);
    2341         xDelete<IssmDouble>(h);
    2342 
    2343 }
    2344 /*}}}*/
     2272}/*}}}*/
    23452273void       Element::MismipFloatingiceMeltingRate(){/*{{{*/
    23462274
    … …  
    25662494
    25672495        const int NUM_VERTICES = this->GetNumberOfVertices();
    2568 
    25692496        int        basinid, maxbox, num_basins, numnodes, M;
    25702497        IssmDouble gamma_T, overturning_coeff, thickness;
    … …  
    25932520        this->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
    25942521        this->parameters->FindParam(&isplume, BasalforcingsPicoIsplumeEnum);
    2595         Input *thickness_input    = this->GetInput(ThicknessEnum);                         _assert_(thickness_input);
    2596         _assert_(basinid<=num_basins);
     2522        Input *thickness_input = this->GetInput(ThicknessEnum);
     2523   _assert_(basinid<=num_basins);
     2524   _assert_(thickness_input);
    25972525
    25982526        IssmDouble* boxareas = NULL;
    … …  
    26032531        IssmDouble g1               = area_boxi*gamma_T;
    26042532
    2605         IssmDouble* basalmeltrates_shelf                                = xNew<IssmDouble>(NUM_VERTICES);
    2606         IssmDouble* potential_pressure_melting_point    = xNew<IssmDouble>(NUM_VERTICES);
    2607         IssmDouble* Tocs                                                                = xNew<IssmDouble>(NUM_VERTICES);
    2608         IssmDouble* Socs                                                                = xNew<IssmDouble>(NUM_VERTICES);
     2533        IssmDouble basalmeltrates_shelf[MAXVERTICES];
     2534        IssmDouble potential_pressure_melting_point[MAXVERTICES];
     2535        IssmDouble Tocs[MAXVERTICES];
     2536        IssmDouble Socs[MAXVERTICES];
    26092537
    26102538        /* First box calculations */
    … …  
    26122540                /* Get box1 parameters and inputs */
    26132541                IssmDouble time, toc_farocean, soc_farocean;
     2542      IssmDouble overturnings[MAXVERTICES];
    26142543                this->parameters->FindParam(&time,TimeEnum);
    26152544                this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum);
    26162545                this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);
    2617                 IssmDouble      s1                              = soc_farocean/(nu*lambda);
    2618                 IssmDouble* overturnings        = xNew<IssmDouble>(NUM_VERTICES);
     2546                IssmDouble s1 = soc_farocean/(nu*lambda);
    26192547                Input *overturningC_input = this->GetInput(BasalforcingsPicoOverturningCoeffEnum); _assert_(overturningC_input);
    26202548
    … …  
    26402568                }
    26412569
    2642                 if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1DGEnum);
    2643                 this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1DGEnum);
    2644                 this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1DGEnum);
    2645                 this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1DGEnum);
     2570                if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,&basalmeltrates_shelf[0],P1DGEnum);
     2571                this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,&Tocs[0],P1DGEnum);
     2572                this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,&Socs[0],P1DGEnum);
     2573                this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,&overturnings[0],P1DGEnum);
    26462574
    26472575                /*Cleanup and return*/
    … …  
    27162644
    27172645        /*Define arrays*/
    2718         IssmDouble* basalmeltrates_shelf = xNew<IssmDouble>(NUM_VERTICES); //Basal melt-rate
     2646        IssmDouble basalmeltrates_shelf[MAXVERTICES];
    27192647
    27202648        /*Polynomial coefficients*/
    … …  
    27892717
    27902718        /*Save computed melt-rate*/
    2791         this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1DGEnum);
     2719        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,&basalmeltrates_shelf[0],P1DGEnum);
    27922720
    27932721        /*Cleanup and return*/
    27942722        delete gauss;
    2795 
    27962723}/*}}}*/
    27972724void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/
    … …  
    42324159        /*Fetch number vertices and allocate memory*/
    42334160        const int NUM_VERTICES  = this->GetNumberOfVertices();
    4234 
    4235         IssmDouble* maxprincipal = xNew<IssmDouble>(NUM_VERTICES);
     4161        IssmDouble maxprincipal[MAXVERTICES];
    42364162
    42374163        /*Retrieve all inputs and parameters*/
    … …  
    42514177        /*loop over vertices: */
    42524178        Gauss* gauss=this->NewGauss();
    4253         for (int iv=0;iv<NUM_VERTICES;iv++){
     4179        for(int iv=0;iv<NUM_VERTICES;iv++){
    42544180                gauss->GaussVertex(iv);
    42554181
    … …  
    42964222
    42974223        /*Create input*/
    4298         this->AddInput(StressMaxPrincipalEnum,maxprincipal,P1Enum);
     4224        this->AddInput(StressMaxPrincipalEnum,&maxprincipal[0],P1Enum);
    42994225
    43004226        /*Clean up and return*/
    4301         xDelete<IssmDouble>(maxprincipal);
    43024227        xDelete<IssmDouble>(xyz_list);
    43034228        delete gauss;
    … …  
    45754500        /*Fetch number vertices and allocate memory*/
    45764501        const int NUM_VERTICES = this->GetNumberOfVertices();
    4577 
    4578         IssmDouble* viscousheating = xNew<IssmDouble>(NUM_VERTICES);
     4502        IssmDouble viscousheating[MAXVERTICES];
    45794503
    45804504        /*Retrieve all inputs and parameters*/
    … …  
    45954519
    45964520        /*Create PentaVertex input, which will hold the basal friction:*/
    4597         this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum);
     4521        this->AddInput(ViscousHeatingEnum,&viscousheating[0],P1Enum);
    45984522
    45994523        /*Clean up and return*/
    4600         xDelete<IssmDouble>(viscousheating);
    46014524        xDelete<IssmDouble>(xyz_list);
    46024525        delete gauss;
    … …  
    46834606        IssmDouble  kappa,kappa_c,kappa_t;  //enthalpy conductivities
    46844607        IssmDouble  Hc,Ht;
    4685         IssmDouble* PIE   = xNew<IssmDouble>(numvertices);
    4686         IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
     4608        IssmDouble PIE[MAXVERTICES];
     4609        IssmDouble dHpmp[MAXVERTICES];
    46874610
    46884611        for(int iv=0; iv<numvertices; iv++){
    … …  
    47204643
    47214644        /*Clean up and return*/
    4722         xDelete<IssmDouble>(PIE);
    4723         xDelete<IssmDouble>(dHpmp);
    47244645        return kappa;
    47254646}

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

    r25839 r25905  
    22222222
    22232223        /*Intermediaries*/
    2224         int i, numiceverts, numnoiceverts;
     2224        int numiceverts, numnoiceverts;
    22252225        int ind0, ind1, lastindex;
    22262226        int indices_ice[NUMVERTICES],indices_noice[NUMVERTICES];
    … …  
    22342234         * Exploit: ice/no-ice parts are connected, so find starting vertex of segment*/
    22352235        lastindex=0;
    2236         for(i=0;i<NUMVERTICES;i++){ // go backwards along vertices, and check for sign change
     2236        for(int i=0;i<NUMVERTICES;i++){ // go backwards along vertices, and check for sign change
    22372237                ind0=(NUMVERTICES-i)%NUMVERTICES;
    22382238                ind1=(ind0-1+NUMVERTICES)%NUMVERTICES;
    … …  
    22482248        numiceverts=0;
    22492249        numnoiceverts=0;
    2250         for(i=0;i<NUMVERTICES;i++){
     2250        for(int i=0;i<NUMVERTICES;i++){
    22512251                ind0=(lastindex+i)%NUMVERTICES;
    22522252                if(lsf[i]<=level){
    … …  
    22602260        }
    22612261        //merge indices
    2262         for(i=0;i<numiceverts;i++){indices[i]=indices_ice[i];}
    2263         for(i=0;i<numnoiceverts;i++){indices[numiceverts+i]=indices_noice[i];}
    2264 
    2265         switch (numiceverts){
     2262        for(int i=0;i<numiceverts;  i++){
     2263                indices[i]=indices_ice[i];
     2264        }
     2265        for(int i=0;i<numnoiceverts;i++){
     2266                indices[numiceverts+i]=indices_noice[i];
     2267        }
     2268
     2269        switch(numiceverts){
    22662270                case 0: // no vertex has ice: element is ice free, no intersection
    2267                         for(i=0;i<2;i++)
     2271                        for(int i=0;i<2;i++)
    22682272                                fraction[i]=0.;
    22692273                        break;
    22702274                case 1: // one vertex has ice:
    2271                         for(i=0;i<2;i++){
     2275                        for(int i=0;i<2;i++){
    22722276                                fraction[i]=(level-lsf[indices[0]])/(lsf[indices[numiceverts+i]]-lsf[indices[0]]);
    22732277                        }
    22742278                        break;
    22752279                case 2: // two vertices have ice: fraction is computed from first ice vertex to last in CCW fashion
    2276                         for(i=0;i<2;i++){
     2280                        for(int i=0;i<2;i++){
    22772281                                fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
    22782282                        }
    22792283                        break;
    22802284                case NUMVERTICES: // all vertices have ice: return triangle area
    2281                         for(i=0;i<2;i++)
     2285                        for(int i=0;i<2;i++)
    22822286                                fraction[i]=1.;
    22832287                        break;
    … …  
    32573261
    32583262        /*intermediary: */
    3259         IssmDouble* values=NULL;
    3260         Input*     thickness_input=NULL;
    32613263        IssmDouble  thickness;
    3262         IssmDouble  weight;
    32633264        IssmDouble  Jdet;
    3264         IssmDouble  volume;
    3265         IssmDouble  rho_ice;
    32663265        int         point1;
    32673266        IssmDouble  fraction1,fraction2;
    32683267        bool        mainlynegative=true;
    32693268
    3270         /*Output:*/
    3271         volume=0;
    3272 
    32733269        /* Get node coordinates and dof list: */
    32743270        IssmDouble  xyz_list[NUMVERTICES][3];
    … …  
    32763272
    32773273        /*Retrieve inputs required:*/
    3278         thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
     3274        Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
    32793275
    32803276        /*Retrieve material parameters: */
    3281         rho_ice=FindParam(MaterialsRhoIceEnum);
     3277        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
    32823278
    32833279        /*Retrieve values of the levelset defining the masscon: */
    3284         values = xNew<IssmDouble>(NUMVERTICES);
     3280        IssmDouble values[NUMVERTICES];
    32853281        for(int i=0;i<NUMVERTICES;i++){
    32863282                values[i]=levelset[this->vertices[i]->Sid()];
    … …  
    32883284
    32893285        /*Ok, use the level set values to figure out where we put our gaussian points:*/
    3290         this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
     3286        this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,&values[0]);
    32913287        Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
    32923288
    3293         volume=0;
    3294 
     3289        IssmDouble volume=0.;
    32953290        while(gauss->next()){
    3296 
    32973291                this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    32983292                thickness_input->GetInputValue(&thickness, gauss);
    3299 
    33003293                volume+=thickness*gauss->weight*Jdet;
    33013294        }
    33023295
    33033296        /* clean up and Return: */
    3304         xDelete<IssmDouble>(values);
    33053297        delete gauss;
    33063298        return rho_ice*volume;
Note: See TracChangeset for help on using the changeset viewer.