Changeset 15654


Ignore:
Timestamp:
07/31/13 15:18:17 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: FS can now use different types of interpolation for velocity and pressure

Location:
issm/trunk-jpl/src
Files:
5 added
1 deleted
31 edited

Legend:

Unmodified
Added
Removed

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

    r15643 r15654  
    865865}
    866866/*}}}*/
     867/*FUNCTION Penta::GetDofListVelocity{{{*/
     868void  Penta::GetDofListVelocity(int** pdoflist,int setenum){
     869
     870        /*Fetch number of nodes and dof for this finite element*/
     871        int numnodes = this->NumberofNodesVelocity();
     872
     873        /*First, figure out size of doflist and create it: */
     874        int numberofdofs=0;
     875        for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
     876
     877        /*Allocate output*/
     878        int* doflist=xNew<int>(numberofdofs);
     879
     880        /*Populate: */
     881        int count=0;
     882        for(int i=0;i<numnodes;i++){
     883                nodes[i]->GetDofList(doflist+count,FSApproximationEnum,setenum);
     884                count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
     885        }
     886
     887        /*Assign output pointers:*/
     888        *pdoflist=doflist;
     889}
     890/*}}}*/
     891/*FUNCTION Penta::GetDofListPressure{{{*/
     892void  Penta::GetDofListPressure(int** pdoflist,int setenum){
     893
     894        /*Fetch number of nodes and dof for this finite element*/
     895        int vnumnodes = this->NumberofNodesVelocity();
     896        int pnumnodes = this->NumberofNodesPressure();
     897
     898        /*First, figure out size of doflist and create it: */
     899        int numberofdofs=0;
     900        for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
     901
     902        /*Allocate output*/
     903        int* doflist=xNew<int>(numberofdofs);
     904
     905        /*Populate: */
     906        int count=0;
     907        for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++){
     908                nodes[i]->GetDofList(doflist+count,FSApproximationEnum,setenum);
     909                count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
     910        }
     911
     912        /*Assign output pointers:*/
     913        *pdoflist=doflist;
     914}
     915/*}}}*/
    867916/*FUNCTION Penta::GetGroundedPortion{{{*/
    868917IssmDouble Penta::GetGroundedPortion(IssmDouble* xyz_list){
    … …  
    10341083
    10351084        _assert_(nodes);
    1036         for(int i=0;i<NUMVERTICES;i++){
    1037                 if(node==nodes[i])
    1038                  return i;
     1085        int numnodes = this->NumberofNodes();
     1086
     1087        for(int i=0;i<numnodes;i++){
     1088                if(node==nodes[i]) return i;
    10391089        }
    10401090        _error_("Node provided not found among element nodes");
    … …  
    13251375         */
    13261376
    1327         int i;
    13281377        IssmDouble epsilonvx[6];
    13291378        IssmDouble epsilonvy[6];
    … …  
    13411390
    13421391        /*Sum all contributions*/
    1343         for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
     1392        for(int i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
    13441393}
    13451394/*}}}*/
    … …  
    26792728     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], pdds, pds,
    26802729                                  signorm, yts, h[iv], s[iv], rho_ice, rho_water, desfac, s0p);
    2681      //printf("mass balance %f \n",agd[iv]);
    26822730   }
    26832731
    … …  
    32253273                        penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
    32263274                        break;
     3275                case P1P1Enum: case P1P1GLSEnum: case MINIcondensedEnum:
     3276                        numnodes         = 12;
     3277                        penta_node_ids   = xNew<int>(numnodes);
     3278                        penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
     3279                        penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1];
     3280                        penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2];
     3281                        penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3];
     3282                        penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4];
     3283                        penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5];
     3284                        penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+0];
     3285                        penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+1];
     3286                        penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+2];
     3287                        penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+3];
     3288                        penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+4];
     3289                        penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+5];
     3290                        break;
     3291                case MINIEnum:
     3292                        numnodes         = 13;
     3293                        penta_node_ids   = xNew<int>(numnodes);
     3294                        penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
     3295                        penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1];
     3296                        penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2];
     3297                        penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3];
     3298                        penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4];
     3299                        penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5];
     3300                        penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+index+1;
     3301                        penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+0];
     3302                        penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+1];
     3303                        penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+2];
     3304                        penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+3];
     3305                        penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+4];
     3306                        penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*index+5];
     3307                        break;
    32273308                default:
    32283309                        _error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");
    … …  
    49775058ElementMatrix* Penta::CreateKMatrixAdjointFS(void){
    49785059
    4979         /*Constants*/
    4980         const int    numdof=NDOF4*NUMVERTICES;
    4981 
    49825060        /*Intermediaries */
    49835061        int        i,j;
    … …  
    49995077        if(incomplete_adjoint) return Ke;
    50005078
     5079        /*Fetch number of nodes and dof for this finite element*/
     5080        int vnumnodes = this->NumberofNodesVelocity();
     5081        int pnumnodes = this->NumberofNodesPressure();
     5082        int numdof    = vnumnodes*NDOF3 + pnumnodes*NDOF1;
     5083
    50015084        /*Retrieve all inputs and parameters*/
    50025085        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    … …  
    50455128
    50465129        /*Transform Coordinate System*/
    5047         TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
     5130        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZEnum);
    50485131
    50495132        /*Clean up and return*/
    … …  
    51035186ElementVector* Penta::CreatePVectorAdjointFS(void){
    51045187
    5105         if (!IsOnSurface()) return NULL;
    5106 
    5107         /*Call Tria function*/
    5108         Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).
    5109         ElementVector* pe=tria->CreatePVectorAdjointFS();
    5110         delete tria->material; delete tria;
    5111 
    5112         /*clean up and return*/
     5188        if(!IsOnSurface()) return NULL;
     5189
     5190        /*Intermediaries */
     5191        int         i,j,resp;
     5192        int        *responses= NULL;
     5193        int         num_responses;
     5194        IssmDouble  Jdet2d;
     5195        IssmDouble  obs_velocity_mag,velocity_mag;
     5196        IssmDouble  dux,duy;
     5197        IssmDouble  epsvel  = 2.220446049250313e-16;
     5198        IssmDouble  meanvel = 3.170979198376458e-05;  /*1000 m/yr */
     5199        IssmDouble  scalex  = 0,scaley=0,scale=0,S=0;
     5200        IssmDouble  vx,vy,vxobs,vyobs,weight;
     5201        IssmDouble  xyz_list[NUMVERTICES][3];
     5202        IssmDouble      xyz_list_tria[NUMVERTICES2D][3];
     5203
     5204        /*Fetch number of nodes and dof for this finite element*/
     5205        int vnumnodes = this->NumberofNodesVelocity();
     5206        int pnumnodes = this->NumberofNodesPressure();
     5207        int vnumdof   = vnumnodes*NDOF3;
     5208
     5209        /*Prepare coordinate system list*/
     5210        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     5211        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     5212        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     5213
     5214        /*Initialize Element matrix and vectors*/
     5215        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
     5216        IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
     5217
     5218        /*Retrieve all inputs and parameters*/
     5219        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     5220        this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     5221        this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
     5222        Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     5223        Input* vx_input      = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
     5224        Input* vy_input      = inputs->GetInput(VyEnum);                                 _assert_(vy_input);
     5225        Input* vxobs_input   = inputs->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
     5226        Input* vyobs_input   = inputs->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
     5227
     5228        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     5229
     5230        /*Get Surface if required by one response*/
     5231        for(resp=0;resp<num_responses;resp++){
     5232                if(responses[resp]==SurfaceAverageVelMisfitEnum){
     5233                        inputs->GetInputValue(&S,SurfaceAreaEnum); break;
     5234                }
     5235        }
     5236
     5237        /* Start  looping on the number of gaussian points: */
     5238        GaussPenta* gauss=new GaussPenta(3,4,5,4);
     5239        for(int ig=gauss->begin();ig<gauss->end();ig++){
     5240
     5241                gauss->GaussPoint(ig);
     5242
     5243                /* Get Jacobian determinant: */
     5244                GetTriaJacobianDeterminant(&Jdet2d,&xyz_list_tria[0][0],gauss);
     5245
     5246                /*Get all parameters at gaussian point*/
     5247                vx_input->GetInputValue(&vx,gauss);
     5248                vy_input->GetInputValue(&vy,gauss);
     5249                vxobs_input->GetInputValue(&vxobs,gauss);
     5250                vyobs_input->GetInputValue(&vyobs,gauss);
     5251                GetNodalFunctionsVelocity(vbasis,gauss);
     5252
     5253                /*Loop over all requested responses*/
     5254                for(resp=0;resp<num_responses;resp++){
     5255
     5256                        weights_input->GetInputValue(&weight,gauss,resp);
     5257
     5258                        switch(responses[resp]){
     5259
     5260                                case SurfaceAbsVelMisfitEnum:
     5261                                        /*
     5262                                         *      1  [           2              2 ]
     5263                                         * J = --- | (u - u   )  +  (v - v   )  |
     5264                                         *      2  [       obs            obs   ]
     5265                                         *
     5266                                         *        dJ
     5267                                         * DU = - -- = (u   - u )
     5268                                         *        du     obs
     5269                                         */
     5270                                        for(i=0;i<vnumnodes;i++){
     5271                                                dux=vxobs-vx;
     5272                                                duy=vyobs-vy;
     5273                                                pe->values[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i];
     5274                                                pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i];
     5275                                        }
     5276                                        break;
     5277                                case SurfaceRelVelMisfitEnum:
     5278                                        /*
     5279                                         *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
     5280                                         * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
     5281                                         *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
     5282                                         *              obs                        obs                     
     5283                                         *
     5284                                         *        dJ     \bar{v}^2
     5285                                         * DU = - -- = ------------- (u   - u )
     5286                                         *        du   (u   + eps)^2    obs
     5287                                         *               obs
     5288                                         */
     5289                                        for(i=0;i<vnumnodes;i++){
     5290                                                scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
     5291                                                scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
     5292                                                dux=scalex*(vxobs-vx);
     5293                                                duy=scaley*(vyobs-vy);
     5294                                                pe->values[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i];
     5295                                                pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i];
     5296                                        }
     5297                                        break;
     5298                                case SurfaceLogVelMisfitEnum:
     5299                                        /*
     5300                                         *                 [        vel + eps     ] 2
     5301                                         * J = 4 \bar{v}^2 | log ( -----------  ) | 
     5302                                         *                 [       vel   + eps    ]
     5303                                         *                            obs
     5304                                         *
     5305                                         *        dJ                 2 * log(...)
     5306                                         * DU = - -- = - 4 \bar{v}^2 -------------  u
     5307                                         *        du                 vel^2 + eps
     5308                                         *           
     5309                                         */
     5310                                        for(i=0;i<vnumnodes;i++){
     5311                                                velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
     5312                                                obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
     5313                                                scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
     5314                                                dux=scale*vx;
     5315                                                duy=scale*vy;
     5316                                                pe->values[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i];
     5317                                                pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i];
     5318                                        }
     5319                                        break;
     5320                                case SurfaceAverageVelMisfitEnum:
     5321                                        /*
     5322                                         *      1                    2              2
     5323                                         * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
     5324                                         *      S                obs            obs
     5325                                         *
     5326                                         *        dJ      1       1
     5327                                         * DU = - -- = - --- ----------- * 2 (u - u   )
     5328                                         *        du      S  2 sqrt(...)           obs
     5329                                         */
     5330                                        for(i=0;i<vnumnodes;i++){
     5331                                                scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
     5332                                                dux=scale*(vxobs-vx);
     5333                                                duy=scale*(vyobs-vy);
     5334                                                pe->values[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i];
     5335                                                pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i];
     5336                                        }
     5337                                        break;
     5338                                case SurfaceLogVxVyMisfitEnum:
     5339                                        /*
     5340                                         *      1            [        |u| + eps     2          |v| + eps     2  ]
     5341                                         * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
     5342                                         *      2            [       |u    |+ eps              |v    |+ eps     ]
     5343                                         *                              obs                       obs
     5344                                         *        dJ                              1      u                             1
     5345                                         * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
     5346                                         *        du                         |u| + eps  |u|                           u + eps
     5347                                         */
     5348                                        for(i=0;i<vnumnodes;i++){
     5349                                                dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
     5350                                                duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
     5351                                                pe->values[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i];
     5352                                                pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i];
     5353                                        }
     5354                                        break;
     5355                                case DragCoefficientAbsGradientEnum:
     5356                                        /*Nothing in P vector*/
     5357                                        break;
     5358                                case ThicknessAbsGradientEnum:
     5359                                        /*Nothing in P vector*/
     5360                                        break;
     5361                                case ThicknessAcrossGradientEnum:
     5362                                        /*Nothing in P vector*/
     5363                                        break;
     5364                                case ThicknessAlongGradientEnum:
     5365                                        /*Nothing in P vector*/
     5366                                        break;
     5367                                case RheologyBbarAbsGradientEnum:
     5368                                        /*Nothing in P vector*/
     5369                                        break;
     5370                                default:
     5371                                        _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
     5372                        }
     5373                }
     5374        }
     5375
     5376        /*Clean up and return*/
     5377        xDelete<int>(responses);
     5378        xDelete<int>(cs_list);
     5379        xDelete<IssmDouble>(vbasis);
     5380        delete gauss;
    51135381        return pe;
    51145382}
    … …  
    54945762void  Penta::InputUpdateFromSolutionAdjointFS(IssmDouble* solution){
    54955763
    5496         const int    numdof=NDOF4*NUMVERTICES;
    5497 
    5498         int    i;
    5499         IssmDouble values[numdof];
    5500         IssmDouble lambdax[NUMVERTICES];
    5501         IssmDouble lambday[NUMVERTICES];
    5502         IssmDouble lambdaz[NUMVERTICES];
    5503         IssmDouble lambdap[NUMVERTICES];
    5504         int*   doflist=NULL;
     5764        int          i;
     5765        int*         vdoflist=NULL;
     5766        int*         pdoflist=NULL;
     5767        IssmDouble   FSreconditioning;
     5768        GaussPenta  *gauss;
     5769
     5770        /*Fetch number of nodes and dof for this finite element*/
     5771        int vnumnodes = this->NumberofNodesVelocityFinal();
     5772        int pnumnodes = this->NumberofNodesPressure();
     5773        int vnumdof   = vnumnodes*NDOF3;
     5774        int pnumdof   = pnumnodes*NDOF1;
     5775
     5776        /*Initialize values*/
     5777        IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
     5778        IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
     5779        IssmDouble* lambdax = xNew<IssmDouble>(vnumnodes);
     5780        IssmDouble* lambday = xNew<IssmDouble>(vnumnodes);
     5781        IssmDouble* lambdaz = xNew<IssmDouble>(vnumnodes);
     5782        IssmDouble* lambdap = xNew<IssmDouble>(pnumnodes);
    55055783
    55065784        /*Get dof list: */
    5507         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     5785        GetDofListVelocity(&vdoflist,GsetEnum);
     5786        GetDofListPressure(&pdoflist,GsetEnum);
    55085787
    55095788        /*Use the dof list to index into the solution vector: */
    5510         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    5511 
    5512         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    5513         for(i=0;i<NUMVERTICES;i++){
    5514                 lambdax[i]=values[i*NDOF4+0];
    5515                 lambday[i]=values[i*NDOF4+1];
    5516                 lambdaz[i]=values[i*NDOF4+2];
    5517                 lambdap[i]=values[i*NDOF4+3];
    5518 
    5519                 /*Check solution*/
    5520                 if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
    5521                 if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
    5522                 if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector");
    5523                 if(xIsNan<IssmDouble>(lambdap[i])) _error_("NaN found in solution vector");
    5524         }
     5789        for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
     5790        for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
     5791
     5792        /*Transform solution in Cartesian Space*/
     5793        TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYZEnum);
     5794
     5795        /*fill in all arrays: */
     5796        for(i=0;i<vnumnodes;i++){
     5797                lambdax[i] = vvalues[i*NDOF3+0]; if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
     5798                lambday[i] = vvalues[i*NDOF3+1]; if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
     5799                lambdaz[i] = vvalues[i*NDOF3+2]; if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector");
     5800        }
     5801        for(i=0;i<pnumnodes;i++){
     5802                lambdap[i] = pvalues[i]; if(xIsNan<IssmDouble>(lambdap[i])) _error_("NaN found in solution vector");
     5803        }
     5804
     5805        /*Recondition pressure and compute vel: */
     5806        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
     5807        for(i=0;i<pnumnodes;i++) lambdap[i]=lambdap[i]*FSreconditioning;
    55255808
    55265809        /*Add vx and vy as inputs to the tria element: */
    … …  
    55315814
    55325815        /*Free ressources:*/
    5533         xDelete<int>(doflist);
     5816        xDelete<int>(vdoflist);
     5817        xDelete<int>(pdoflist);
     5818        xDelete<IssmDouble>(lambdap);
     5819        xDelete<IssmDouble>(lambdaz);
     5820        xDelete<IssmDouble>(lambday);
     5821        xDelete<IssmDouble>(lambdax);
     5822        xDelete<IssmDouble>(vvalues);
     5823        xDelete<IssmDouble>(pvalues);
    55345824}
    55355825/*}}}*/
    … …  
    60516341        /*output: */
    60526342        ElementVector* De=NULL;
    6053         /*intermediary: */
     6343
     6344        /*Initialize Element vector and return if necessary*/
    60546345        int approximation;
    6055         int i;
    6056 
    6057         /*Initialize Element vector and return if necessary*/
    60586346        inputs->GetInputValue(&approximation,ApproximationEnum);
    60596347        if(approximation!=FSApproximationEnum) return NULL;
    60606348
    6061         De=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    6062 
    6063         for (i=0;i<NUMVERTICES;i++){
    6064                 De->values[i*4+0]=VelocityEnum;
    6065                 De->values[i*4+1]=VelocityEnum;
    6066                 De->values[i*4+2]=VelocityEnum;
    6067                 De->values[i*4+3]=PressureEnum;
     6349        /*Fetch number of nodes and dof for this finite element*/
     6350        int vnumnodes = this->NumberofNodesVelocity();
     6351        int pnumnodes = this->NumberofNodesPressure();
     6352
     6353        De=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
     6354
     6355        for(int i=0;i<vnumnodes;i++){
     6356                De->values[i*3+0]=VelocityEnum;
     6357                De->values[i*3+1]=VelocityEnum;
     6358                De->values[i*3+2]=VelocityEnum;
     6359        }
     6360        for(int i=0;i<pnumnodes;i++){
     6361                De->values[vnumnodes*3+i]=PressureEnum;
    60686362        }
    60696363
    … …  
    63176611                node_list[i+1*NUMVERTICES] = this->nodes[i];
    63186612                cs_list[i+0*NUMVERTICES] = XYEnum;
    6319                 cs_list[i+1*NUMVERTICES] = XYZPEnum;
     6613                cs_list[i+1*NUMVERTICES] = XYZEnum;
    63206614        }
    63216615
    … …  
    64276721                node_list[i+1*NUMVERTICES] = this->nodes[i];
    64286722                cs_list[i+0*NUMVERTICES] = XYEnum;
    6429                 cs_list[i+1*NUMVERTICES] = XYZPEnum;
     6723                cs_list[i+1*NUMVERTICES] = XYZEnum;
    64306724        }
    64316725
    … …  
    65166810                node_list[i+1*NUMVERTICES] = this->nodes[i];
    65176811                cs_list[i+0*NUMVERTICES] = XYEnum;
    6518                 cs_list[i+1*NUMVERTICES] = XYZPEnum;
     6812                cs_list[i+1*NUMVERTICES] = XYZEnum;
    65196813        }
    65206814
    … …  
    65266820        delete Ke2;
    65276821        Ke1=CreateKMatrixDiagnosticHO(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum);
    6528         Ke2=CreateKMatrixDiagnosticFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,NUMVERTICES,XYZPEnum);
     6822        Ke2=CreateKMatrixDiagnosticFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,NUMVERTICES,XYZEnum);
    65296823
    65306824        for(i=0;i<numdofs;i++) for(j=0;j<NUMVERTICES;j++){
    … …  
    70997393ElementMatrix* Penta::CreateKMatrixDiagnosticFS(void){
    71007394
    7101         int fe_FS;
    7102         ElementMatrix* Ke1;
    7103         ElementMatrix* Ke2;
    7104         ElementMatrix* Ke;
    7105         parameters->FindParam(&fe_FS,FlowequationFeFSEnum);
    7106 
    7107         switch(fe_FS){
    7108                 case 0:
    7109                         /*compute all stiffness matrices for this element*/
    7110                         Ke1=CreateKMatrixDiagnosticFSViscous();
    7111                         Ke2=CreateKMatrixDiagnosticFSFriction();
    7112                         Ke =new ElementMatrix(Ke1,Ke2);
    7113                         break;
    7114                 case 1:
    7115                         /*compute all stiffness matrices for this element*/
    7116                         Ke1=CreateKMatrixDiagnosticFSGLSViscous();
    7117                         Ke2=CreateKMatrixDiagnosticFSFriction();
    7118                         Ke =new ElementMatrix(Ke1,Ke2);
    7119                         break;
    7120                 default:
    7121                         _error_("Finite element" << fe_FS << " not supported yet");
    7122         }
     7395        ElementMatrix* Ke1 = NULL;
     7396        ElementMatrix* Ke2 = NULL;
     7397        ElementMatrix* Ke  = NULL;
     7398
     7399        /*compute all stiffness matrices for this element*/
     7400        Ke1=CreateKMatrixDiagnosticFSViscous();
     7401        Ke2=CreateKMatrixDiagnosticFSFriction();
     7402        Ke =new ElementMatrix(Ke1,Ke2);
    71237403
    71247404        /*clean-up and return*/
    … …  
    71297409}
    71307410/*}}}*/
    7131 /*FUNCTION Penta::CreateKMatrixDiagnosticFSViscous {{{*/
    7132 ElementMatrix* Penta::CreateKMatrixDiagnosticFSViscous(void){
    7133 
    7134         /*Intermediaries */
    7135         int        i,approximation;
    7136         IssmDouble Jdet,viscosity,FSreconditioning;
    7137         IssmDouble xyz_list[NUMVERTICES][3];
    7138         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    7139         IssmDouble B[8][27];
    7140         IssmDouble B_prime[8][27];
    7141         IssmDouble D_scalar;
    7142         IssmDouble D[8][8]={0.0};
    7143         IssmDouble Ke_temp[27][27]={0.0}; //for the six nodes and the bubble
    7144         GaussPenta *gauss=NULL;
    7145 
    7146         /*If on water or not FS, skip stiffness: */
    7147         inputs->GetInputValue(&approximation,ApproximationEnum);
    7148         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    7149         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    7150 
    7151         /*Retrieve all inputs and parameters*/
    7152         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7153         parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    7154         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    7155         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    7156         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    7157 
    7158         /* Start  looping on the number of gaussian points: */
    7159         gauss=new GaussPenta(5,5);
    7160         for(int ig=gauss->begin();ig<gauss->end();ig++){
    7161 
    7162                 gauss->GaussPoint(ig);
    7163 
    7164                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7165                 GetBFS(&B[0][0],&xyz_list[0][0],gauss);
    7166                 GetBprimeFS(&B_prime[0][0],&xyz_list[0][0],gauss);
    7167 
    7168                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    7169                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    7170 
    7171                 D_scalar=gauss->weight*Jdet;
    7172                 for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
    7173                 for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning;
    7174 
    7175                 TripleMultiply( &B[0][0],8,27,1,
    7176                                         &D[0][0],8,8,0,
    7177                                         &B_prime[0][0],8,27,0,
    7178                                         &Ke_temp[0][0],1);
    7179         }
    7180 
    7181         /*Condensation*/
    7182         ReduceMatrixFS(Ke->values, &Ke_temp[0][0]);
    7183         //Ke->Echo();
    7184         //_error_("stop");
    7185 
    7186         /*Transform Coordinate System*/
    7187         TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
    7188 
    7189         /*Clean up and return*/
    7190         delete gauss;
    7191         return Ke;
    7192 }
    7193 /*}}}*/
    7194 /*FUNCTION Penta::CreateKMatrixDiagnosticFSGLSViscous {{{*/
    7195 ElementMatrix* Penta::CreateKMatrixDiagnosticFSGLSViscous(void){
     7411/*FUNCTION Penta::KMatrixGLSstabilization{{{*/
     7412void Penta::KMatrixGLSstabilization(ElementMatrix* Ke){
    71967413
    71977414        int        numdof  = NUMVERTICES*NDOF4;
    … …  
    72237440        int c=3; //index of pressure
    72247441
    7225         /*If on water or not FS, skip stiffness: */
    7226         inputs->GetInputValue(&approximation,ApproximationEnum);
    7227         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    7228         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    7229 
    72307442        /*Retrieve all inputs and parameters*/
    72317443        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    … …  
    72397451        diameter=MinEdgeLength(xyz_list);
    72407452
    7241         if(stabilization){
    72427453                gauss=new GaussPenta();
    72437454                for(int iv=0;iv<6;iv++){
    … …  
    72517462                }
    72527463                delete gauss;
    7253         }
    72547464
    72557465        /* Start  looping on the number of gaussian points: */
    … …  
    72607470
    72617471                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7262                 GetBFSGLS(&B[0][0],&xyz_list[0][0],gauss);
    7263                 GetBprimeFSGLS(&B_prime[0][0],&xyz_list[0][0],gauss);
    72647472
    72657473                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    … …  
    72677475
    72687476                D_scalar=gauss->weight*Jdet;
    7269                 for (i=0;i<6;i++) D[i][i]=D_scalar*2.*2.*viscosity;
    7270                 for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning;
    7271 
    7272                 TripleMultiply( &B[0][0],8,24,1,
    7273                                         &D[0][0],8,8,0,
    7274                                         &B_prime[0][0],8,24,0,
    7275                                         &Ke_temp[0][0],0);
    7276 
    7277                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_temp[i][j];
     7477
    72787478
    72797479                /*Add stabilization*/
    7280                 if(stabilization){
    7281                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    7282                         dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
    7283                         mu = 2.*viscosity;
    7284                         for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
    7285                                 SU[p][i][j]=0.;
    7286                                 SW[p][i][j]=0.;
    7287                         }
    7288                         for(p=0;p<6;p++){
    7289                                 for(i=0;i<3;i++){
    7290                                         SU[p][i][c] += dbasis[i][p];
    7291                                         SW[p][c][i] += dbasis[i][p];
    7292                                         for(j=0;j<3;j++){
    7293                                                 SU[p][i][i] += -dmu[j]*dbasis[j][p];
    7294                                                 SU[p][i][j] += -dmu[i]*dbasis[i][p];
    7295                                                 for(ii=0;ii<6;ii++){
    7296                                                         SU[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
    7297                                                         SU[p][i][j] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
    7298                                                 }
    7299                                                 SW[p][i][i] += -dmu[j]*dbasis[j][p];
    7300                                                 SW[p][j][i] += -dmu[j]*dbasis[i][p];
    7301                                                 for(ii=0;ii<6;ii++){
    7302                                                         SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
    7303                                                         SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
    7304                                                 }
     7480                GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
     7481                dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
     7482                mu = 2.*viscosity;
     7483                for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
     7484                        SU[p][i][j]=0.;
     7485                        SW[p][i][j]=0.;
     7486                }
     7487                for(p=0;p<6;p++){
     7488                        for(i=0;i<3;i++){
     7489                                SU[p][i][c] += dbasis[i][p];
     7490                                SW[p][c][i] += dbasis[i][p];
     7491                                for(j=0;j<3;j++){
     7492                                        SU[p][i][i] += -dmu[j]*dbasis[j][p];
     7493                                        SU[p][i][j] += -dmu[i]*dbasis[i][p];
     7494                                        for(ii=0;ii<6;ii++){
     7495                                                SU[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
     7496                                                SU[p][i][j] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
     7497                                        }
     7498                                        SW[p][i][i] += -dmu[j]*dbasis[j][p];
     7499                                        SW[p][j][i] += -dmu[j]*dbasis[i][p];
     7500                                        for(ii=0;ii<6;ii++){
     7501                                                SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
     7502                                                SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
    73057503                                        }
    73067504                                }
    73077505                        }
    7308                         IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity);
    7309                         for(p=0;p<6;p++){
    7310                                 for(q=0;q<6;q++){
    7311                                         for(i=0;i<4;i++){
    7312                                                 for(j=0;j<4;j++){
    7313                                                         Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][0]*SU[q][0][j];
    7314                                                         Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][1]*SU[q][1][j];
    7315                                                         Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][2]*SU[q][2][j];
    7316                                                 }
     7506                }
     7507                IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity);
     7508                for(p=0;p<6;p++){
     7509                        for(q=0;q<6;q++){
     7510                                for(i=0;i<4;i++){
     7511                                        for(j=0;j<4;j++){
     7512                                                Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][0]*SU[q][0][j];
     7513                                                Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][1]*SU[q][1][j];
     7514                                                Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][2]*SU[q][2][j];
    73177515                                        }
    73187516                                }
    73197517                        }
    7320 
    7321                 }
     7518                }
     7519        }
     7520
     7521        /*Clean up*/
     7522        delete gauss;
     7523}
     7524/*}}}*/
     7525/*FUNCTION Penta::CreateKMatrixDiagnosticFSViscous {{{*/
     7526ElementMatrix* Penta::CreateKMatrixDiagnosticFSViscous(void){
     7527
     7528        /*Intermediaries */
     7529        int        i,j,approximation;
     7530        IssmDouble Jdet,viscosity,FSreconditioning,D_scalar;
     7531        IssmDouble xyz_list[NUMVERTICES][3];
     7532        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     7533        GaussPenta *gauss=NULL;
     7534
     7535        /*If on water or not FS, skip stiffness: */
     7536        inputs->GetInputValue(&approximation,ApproximationEnum);
     7537        if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     7538
     7539        /*Fetch number of nodes and dof for this finite element*/
     7540        int vnumnodes = this->NumberofNodesVelocity();
     7541        int pnumnodes = this->NumberofNodesPressure();
     7542        int numdof    = vnumnodes*NDOF3 + pnumnodes*NDOF1;
     7543
     7544        /*Prepare coordinate system list*/
     7545        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     7546        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     7547        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     7548
     7549        /*Initialize Element matrix and vectors*/
     7550        ElementMatrix* Ke     = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
     7551        IssmDouble*    B      = xNew<IssmDouble>(8*numdof);
     7552        IssmDouble*    Bprime = xNew<IssmDouble>(8*numdof);
     7553        IssmDouble*    D      = xNewZeroInit<IssmDouble>(8*8);
     7554
     7555        /*Retrieve all inputs and parameters*/
     7556        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     7557        parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
     7558        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     7559        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     7560        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     7561
     7562        /* Start  looping on the number of gaussian points: */
     7563        gauss=new GaussPenta(5,5);
     7564        for(int ig=gauss->begin();ig<gauss->end();ig++){
     7565
     7566                gauss->GaussPoint(ig);
     7567
     7568                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     7569                GetBFS(B,&xyz_list[0][0],gauss);
     7570                GetBprimeFS(Bprime,&xyz_list[0][0],gauss);
     7571
     7572                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     7573                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
     7574
     7575                D_scalar=gauss->weight*Jdet;
     7576                for(i=0;i<6;i++) D[i*8+i] = +D_scalar*2.*viscosity;
     7577                for(i=6;i<8;i++) D[i*8+i] = -D_scalar*FSreconditioning;
     7578
     7579                TripleMultiply(B,8,numdof,1,
     7580                                        D,8,8,0,
     7581                                        Bprime,8,numdof,0,
     7582                                        Ke->values,1);
    73227583        }
    73237584
    73247585        /*Transform Coordinate System*/
    7325         TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
     7586        TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);
    73267587
    73277588        /*Clean up and return*/
    73287589        delete gauss;
     7590        xDelete<int>(cs_list);
     7591        xDelete<IssmDouble>(B);
     7592        xDelete<IssmDouble>(Bprime);
     7593        xDelete<IssmDouble>(D);
    73297594        return Ke;
    73307595}
    … …  
    73337598ElementMatrix* Penta::CreateKMatrixDiagnosticFSFriction(void){
    73347599
    7335         /*Constants*/
    7336         const int numdof=NUMVERTICES*NDOF4;
    7337         const int numdof2d=NUMVERTICES2D*NDOF4;
    7338 
    73397600        /*Intermediaries */
    7340         int        i,j;
    7341         int        analysis_type,approximation;
    7342         IssmDouble alpha2,Jdet2d;
    7343         IssmDouble FSreconditioning,viscosity;
    7344         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    7345         IssmDouble xyz_list[NUMVERTICES][3];
    7346         IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    7347         IssmDouble LFS[2][numdof2d];
    7348         IssmDouble DLFS[2][2]={0.0};
    7349         IssmDouble Ke_drag_gaussian[numdof2d][numdof2d];
    7350         Friction*  friction=NULL;
    7351         GaussPenta *gauss=NULL;
     7601        int         i,j;
     7602        int         analysis_type,approximation;
     7603        IssmDouble  alpha2,Jdet2d;
     7604        IssmDouble  FSreconditioning,viscosity;
     7605        IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     7606        IssmDouble  xyz_list[NUMVERTICES][3];
     7607        IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
     7608        Friction   *friction = NULL;
     7609        GaussPenta *gauss    = NULL;
    73527610
    73537611        /*If on water or not FS, skip stiffness: */
    73547612        inputs->GetInputValue(&approximation,ApproximationEnum);
    7355         if(IsFloating() || !IsOnBed() || (approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum &&  approximation!=HOFSApproximationEnum)) return NULL;
    7356         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
     7613        if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     7614
     7615        /*Fetch number of nodes and dof for this finite element*/
     7616        int vnumnodes = this->NumberofNodesVelocity();
     7617        int pnumnodes = this->NumberofNodesPressure();
     7618        int  numdof   = vnumnodes*NDOF3 + pnumnodes*NDOF1;
     7619
     7620        /*Initialize Element matrix and vectors*/
     7621        ElementMatrix* Ke        = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
     7622        IssmDouble*    BFriction = xNew<IssmDouble>(2*numdof);
     7623        IssmDouble*    D         = xNewZeroInit<IssmDouble>(2*2);
    73577624
    73587625        /*Retrieve all inputs and parameters*/
    … …  
    73757642
    73767643                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    7377                 GetLFS(&LFS[0][0], gauss);
     7644                GetLFS(BFriction,gauss);
    73787645
    73797646                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    73807647                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    73817648
    7382                 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
    7383 
    7384                 DLFS[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
    7385                 DLFS[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy
    7386 
    7387                 TripleMultiply( &LFS[0][0],2,numdof2d,1,
    7388                                         &DLFS[0][0],2,2,0,
    7389                                         &LFS[0][0],2,numdof2d,0,
    7390                                         &Ke_drag_gaussian[0][0],0);
    7391 
    7392                 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof+j]+=Ke_drag_gaussian[i][j];
     7649                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     7650
     7651                D[0*2+0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
     7652                D[1*2+1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy
     7653
     7654                TripleMultiply(BFriction,2,numdof,1,
     7655                                        D,2,2,0,
     7656                                        BFriction,2,numdof,0,
     7657                                        Ke->values,1);
    73937658        }
    73947659
    73957660        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
    7396         //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
     7661        //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZEnum);
    73977662
    73987663        /*Clean up and return*/
    73997664        delete gauss;
    74007665        delete friction;
     7666        xDelete<IssmDouble>(BFriction);
     7667        xDelete<IssmDouble>(D);
    74017668        return Ke;
    74027669}
    … …  
    75767843
    75777844        /*Transform coordinate system*/
    7578         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
     7845        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum);
    75797846
    75807847        /*Clean up and return*/
    … …  
    76507917
    76517918        /*Transform coordinate system*/
    7652         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
     7919        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum);
    76537920
    76547921        /*Clean up and return*/
    … …  
    77277994
    77287995        /*Transform coordinate system*/
    7729         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
     7996        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum);
    77307997
    77317998        /*Clean up and return*/
    … …  
    78018068
    78028069        /*Transform coordinate system*/
    7803         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
     8070        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum);
    78048071
    78058072        /*Clean up and return*/
    … …  
    81478414ElementVector* Penta::CreatePVectorDiagnosticFS(void){
    81488415
    8149         int fe_FS;
    81508416        ElementVector* pe1;
    81518417        ElementVector* pe2;
    81528418        ElementVector* pe3;
    81538419        ElementVector* pe;
    8154         parameters->FindParam(&fe_FS,FlowequationFeFSEnum);
    8155 
    8156         switch(fe_FS){
    8157                 case 0:
    8158                         /*compute all stiffness matrices for this element*/
    8159                         pe1=CreatePVectorDiagnosticFSViscous();
    8160                         pe2=CreatePVectorDiagnosticFSShelf();
    8161                         pe3=CreatePVectorDiagnosticFSFront();
    8162                         pe =new ElementVector(pe1,pe2,pe3);
    8163                         break;
    8164                 case 1:
    8165                         /*compute all stiffness matrices for this element*/
    8166                         pe1=CreatePVectorDiagnosticFSGLSViscous();
    8167                         pe2=CreatePVectorDiagnosticFSShelf();
    8168                         pe3=CreatePVectorDiagnosticFSFront();
    8169                         pe =new ElementVector(pe1,pe2,pe3);
    8170                         break;
    8171                 default:
    8172                         _error_("Finite element" << fe_FS << " not supported yet");
    8173         }
     8420
     8421        /*compute all stiffness matrices for this element*/
     8422        pe1=CreatePVectorDiagnosticFSViscous();
     8423        pe2=CreatePVectorDiagnosticFSShelf();
     8424        pe3=CreatePVectorDiagnosticFSFront();
     8425        pe =new ElementVector(pe1,pe2,pe3);
    81748426
    81758427        /*clean-up and return*/
    … …  
    81808432}
    81818433/*}}}*/
    8182 /*FUNCTION Penta::CreatePVectorDiagnosticFSViscous {{{*/
    8183 ElementVector* Penta::CreatePVectorDiagnosticFSViscous(void){
    8184 
    8185         /*Constants*/
    8186         const int numdofbubble=NDOF4*NUMVERTICES+NDOF3*1;
    8187 
    8188         /*Intermediaries*/
    8189         int        i,j;
    8190         int        approximation;
    8191         IssmDouble Jdet,viscosity;
    8192         IssmDouble gravity,rho_ice,FSreconditioning;
    8193         IssmDouble forcex,forcey,forcez;
    8194         IssmDouble xyz_list[NUMVERTICES][3];
    8195         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    8196         IssmDouble l1l7[7]; //for the six nodes and the bubble
    8197         IssmDouble B[8][numdofbubble];
    8198         IssmDouble B_prime[8][numdofbubble];
    8199         IssmDouble B_prime_bubble[8][3];
    8200         IssmDouble D[8][8]={0.0};
    8201         IssmDouble D_scalar;
    8202         IssmDouble Pe_gaussian[numdofbubble]={0.0}; //for the six nodes and the bubble
    8203         IssmDouble Ke_temp[numdofbubble][3]={0.0}; //for the six nodes and the bubble
    8204         GaussPenta *gauss=NULL;
    8205 
    8206         /*Initialize Element vector and return if necessary*/
    8207         inputs->GetInputValue(&approximation,ApproximationEnum);
    8208         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    8209         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    8210 
    8211         /*Retrieve all inputs and parameters*/
    8212         this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    8213         rho_ice=matpar->GetRhoIce();
    8214         gravity=matpar->GetG();
    8215         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8216         Input* vx_input=inputs->GetInput(VxEnum);   _assert_(vx_input);
    8217         Input* vy_input=inputs->GetInput(VyEnum);   _assert_(vy_input);
    8218         Input* vz_input=inputs->GetInput(VzEnum);   _assert_(vz_input);
    8219         Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
    8220         Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
    8221         Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
    8222 
    8223         /* Start  looping on the number of gaussian points: */
    8224         gauss=new GaussPenta(5,5);
    8225         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8226 
    8227                 gauss->GaussPoint(ig);
    8228 
    8229                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8230                 GetBFS(&B[0][0],&xyz_list[0][0],gauss);
    8231                 GetBprimeFS(&B_prime[0][0],&xyz_list[0][0], gauss);
    8232                 GetNodalFunctionsMINI(&l1l7[0], gauss);
    8233 
    8234                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    8235                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    8236 
    8237                 loadingforcex_input->GetInputValue(&forcex, gauss);
    8238                 loadingforcey_input->GetInputValue(&forcey, gauss);
    8239                 loadingforcez_input->GetInputValue(&forcez, gauss);
    8240 
    8241                 for(i=0;i<NUMVERTICES+1;i++){
    8242                         Pe_gaussian[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l7[i];
    8243                         Pe_gaussian[i*NDOF4+0]+=rho_ice*forcex*Jdet*gauss->weight*l1l7[i];
    8244                         Pe_gaussian[i*NDOF4+1]+=rho_ice*forcey*Jdet*gauss->weight*l1l7[i];
    8245                         Pe_gaussian[i*NDOF4+2]+=rho_ice*forcez*Jdet*gauss->weight*l1l7[i];
    8246                 }
    8247 
    8248                 /*Get bubble part of Bprime */
    8249                 for(i=0;i<8;i++) for(j=0;j<3;j++) B_prime_bubble[i][j]=B_prime[i][j+24];
    8250 
    8251                 D_scalar=gauss->weight*Jdet;
    8252                 for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
    8253                 for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning;
    8254 
    8255                 TripleMultiply(&B[0][0],8,numdofbubble,1,
    8256                                         &D[0][0],8,8,0,
    8257                                         &B_prime_bubble[0][0],8,3,0,
    8258                                         &Ke_temp[0][0],1);
    8259         }
    8260 
    8261         /*Condensation*/
    8262         ReduceVectorFS(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]);
    8263 
    8264         /*Transform coordinate system*/
    8265         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
    8266 
    8267         /*Clean up and return*/
    8268         delete gauss;
    8269         return pe;
    8270 }
    8271 /*}}}*/
    82728434/*FUNCTION Penta::CreatePVectorDiagnosticFSFront{{{*/
    82738435ElementVector* Penta::CreatePVectorDiagnosticFSFront(void){
    82748436
    82758437        /*Intermediaries */
     8438        int         i;
    82768439        IssmDouble  ls[NUMVERTICES];
    82778440        IssmDouble  xyz_list[NUMVERTICES][3];
    … …  
    82938456
    82948457        /*Fetch number of nodes and dof for this finite element*/
    8295         int         numnodes = this->NumberofNodes();
    8296         int         numdof   = numnodes*NDOF4;
    82978458        IssmDouble  rho_ice,rho_water,gravity;
    82988459        IssmDouble  Jdet,z_g,water_pressure,air_pressure;
    82998460        IssmDouble  surface_under_water,base_under_water,pressure;
    8300         GaussPenta* gauss;
    8301         IssmDouble* basis = xNew<IssmDouble>(numnodes);
    83028461        IssmDouble  xyz_list_front[4][3];
    83038462        IssmDouble  area_coordinates[4][3];
    83048463        IssmDouble  normal[3];
    8305 
     8464        GaussPenta* gauss;
     8465
     8466        /*Fetch number of nodes and dof for this finite element*/
     8467        int vnumnodes = this->NumberofNodesVelocity();
     8468        int pnumnodes = this->NumberofNodesPressure();
     8469        int vnumdof   = vnumnodes*NDOF3;
     8470
     8471        /*Prepare coordinate system list*/
     8472        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     8473        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     8474        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     8475
     8476        /*Initialize Element matrix and vectors*/
     8477        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
     8478        IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
     8479
     8480        /*Retrieve all inputs and parameters*/
    83068481        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8307         ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,FSApproximationEnum);
    83088482        rho_water=matpar->GetRhoWater();
    83098483        rho_ice  =matpar->GetRhoIce();
    … …  
    83238497                gauss->GaussPoint(ig);
    83248498                z_g=GetZcoord(gauss);
    8325                 GetNodalFunctions(basis,gauss);
     8499                GetNodalFunctionsVelocity(vbasis,gauss);
    83268500                GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss);
    83278501
    … …  
    83308504                pressure = water_pressure + air_pressure;
    83318505
    8332                 for (int i=0;i<numnodes;i++){
    8333                         pe->values[4*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
    8334                         pe->values[4*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
    8335                         pe->values[4*i+2]+= pressure*Jdet*gauss->weight*normal[2]*basis[i];
    8336                         pe->values[4*i+3]+= 0;
     8506                for(i=0;i<vnumnodes;i++){
     8507                        pe->values[3*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
     8508                        pe->values[3*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
     8509                        pe->values[3*i+2]+= pressure*Jdet*gauss->weight*normal[2]*vbasis[i];
    83378510                }
    83388511        }
    83398512
    83408513        /*Transform coordinate system*/
    8341         TransformLoadVectorCoord(pe,nodes,numnodes,XYZPEnum);
     8514        TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
    83428515
    83438516        /*Clean up and return*/
    8344         xDelete<IssmDouble>(basis);
     8517        xDelete<int>(cs_list);
     8518        xDelete<IssmDouble>(vbasis);
    83458519        delete gauss;
    83468520        return pe;
    83478521}
    83488522/*}}}*/
    8349 /*FUNCTION Penta::CreatePVectorDiagnosticFSGLSViscous {{{*/
    8350 ElementVector* Penta::CreatePVectorDiagnosticFSGLSViscous(void){
     8523/*FUNCTION Penta::PVectorGLSstabilization{{{*/
     8524void Penta::PVectorGLSstabilization(ElementVector* pe){
    83518525
    83528526        /*Constants*/
    … …  
    83748548        int c=3; //index of pressure
    83758549
    8376         /*Initialize Element vector and return if necessary*/
    8377         inputs->GetInputValue(&approximation,ApproximationEnum);
    8378         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    8379         parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    8380         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    8381 
    83828550        /*Retrieve all inputs and parameters*/
    83838551        rho_ice=matpar->GetRhoIce();
    … …  
    83858553        B=material->GetB();
    83868554        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8387         Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
    8388         Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
    8389         Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
    83908555        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    83918556        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    … …  
    83958560        diameter=MinEdgeLength(xyz_list);
    83968561
    8397         if(stabilization){
    83988562                gauss=new GaussPenta();
    83998563                for(int iv=0;iv<6;iv++){
    … …  
    84078571                }
    84088572                delete gauss;
    8409         }
    84108573
    84118574        /* Start  looping on the number of gaussian points: */
    … …  
    84208583                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    84218584
    8422                 loadingforcex_input->GetInputValue(&forcex, gauss);
    8423                 loadingforcey_input->GetInputValue(&forcey, gauss);
    8424                 loadingforcez_input->GetInputValue(&forcez, gauss);
    8425 
    8426                 for(i=0;i<NUMVERTICES;i++){
    8427                         pe->values[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l6[i];
    8428                         pe->values[i*NDOF4+0]+=rho_ice*forcex*Jdet*gauss->weight*l1l6[i];
    8429                         pe->values[i*NDOF4+1]+=rho_ice*forcey*Jdet*gauss->weight*l1l6[i];
    8430                         pe->values[i*NDOF4+2]+=rho_ice*forcez*Jdet*gauss->weight*l1l6[i];
    8431                 }
    8432 
    8433                 /*Add stabilization*/
    8434                 if(stabilization){
    8435                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    8436                         dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
    8437                         mu = 2.*viscosity;
    8438                         for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
    8439                                 SW[p][i][j]=0.;
    8440                         }
    8441                         for(p=0;p<6;p++){
    8442                                 for(i=0;i<3;i++){
    8443                                         SW[p][c][i] += dbasis[i][p];
    8444                                         for(j=0;j<3;j++){
    8445                                                 SW[p][i][i] += -dmu[j]*dbasis[j][p];
    8446                                                 SW[p][j][i] += -dmu[j]*dbasis[i][p];
    8447                                                 for(ii=0;ii<6;ii++){
    8448                                                         SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
    8449                                                         SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
    8450                                                 }
     8585                GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
     8586                dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
     8587                mu = 2.*viscosity;
     8588                for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
     8589                        SW[p][i][j]=0.;
     8590                }
     8591                for(p=0;p<6;p++){
     8592                        for(i=0;i<3;i++){
     8593                                SW[p][c][i] += dbasis[i][p];
     8594                                for(j=0;j<3;j++){
     8595                                        SW[p][i][i] += -dmu[j]*dbasis[j][p];
     8596                                        SW[p][j][i] += -dmu[j]*dbasis[i][p];
     8597                                        for(ii=0;ii<6;ii++){
     8598                                                SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
     8599                                                SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
    84518600                                        }
    84528601                                }
    84538602                        }
    8454                         IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity);
    8455                         for(p=0;p<6;p++){
    8456                                 for(j=0;j<4;j++){
    8457                                         pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcex*SW[p][j][0];
    8458                                         pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcey*SW[p][j][1];
    8459                                         pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*(forcez-gravity)*SW[p][j][2];
    8460                                 }
     8603                }
     8604                IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity);
     8605                for(p=0;p<6;p++){
     8606                        for(j=0;j<4;j++){
     8607                                pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcex*SW[p][j][0];
     8608                                pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcey*SW[p][j][1];
     8609                                pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*(forcez-gravity)*SW[p][j][2];
    84618610                        }
    8462 
     8611                }
     8612        }
     8613
     8614        /*Clean up*/
     8615        delete gauss;
     8616}
     8617/*}}}*/
     8618/*FUNCTION Penta::CreatePVectorDiagnosticFSViscous {{{*/
     8619ElementVector* Penta::CreatePVectorDiagnosticFSViscous(void){
     8620
     8621        /*Intermediaries*/
     8622        int        i,j;
     8623        int        approximation;
     8624        IssmDouble Jdet,gravity,rho_ice;
     8625        IssmDouble forcex,forcey,forcez;
     8626        IssmDouble xyz_list[NUMVERTICES][3];
     8627        GaussPenta *gauss=NULL;
     8628
     8629        /*Initialize Element vector and return if necessary*/
     8630        inputs->GetInputValue(&approximation,ApproximationEnum);
     8631        if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     8632
     8633        /*Fetch number of nodes and dof for this finite element*/
     8634        int vnumnodes = this->NumberofNodesVelocity();
     8635        int pnumnodes = this->NumberofNodesPressure();
     8636        int vnumdof   = vnumnodes*NDOF3;
     8637
     8638        /*Prepare coordinate system list*/
     8639        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     8640        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     8641        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     8642
     8643        /*Initialize Element matrix and vectors*/
     8644        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
     8645        IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
     8646
     8647        /*Retrieve all inputs and parameters*/
     8648        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     8649        Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
     8650        Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
     8651        Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
     8652        rho_ice = matpar->GetRhoIce();
     8653        gravity = matpar->GetG();
     8654
     8655        /* Start  looping on the number of gaussian points: */
     8656        gauss=new GaussPenta(5,5);
     8657        for(int ig=gauss->begin();ig<gauss->end();ig++){
     8658
     8659                gauss->GaussPoint(ig);
     8660
     8661                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     8662                GetNodalFunctionsVelocity(vbasis, gauss);
     8663
     8664                loadingforcex_input->GetInputValue(&forcex,gauss);
     8665                loadingforcey_input->GetInputValue(&forcey,gauss);
     8666                loadingforcez_input->GetInputValue(&forcez,gauss);
     8667
     8668                for(i=0;i<vnumnodes;i++){
     8669                        pe->values[i*NDOF3+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
     8670                        pe->values[i*NDOF3+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
     8671                        pe->values[i*NDOF3+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
     8672                        pe->values[i*NDOF3+2] += +rho_ice*forcez *Jdet*gauss->weight*vbasis[i];
    84638673                }
    84648674        }
    84658675
    84668676        /*Transform coordinate system*/
    8467         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
     8677        TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
    84688678
    84698679        /*Clean up and return*/
    84708680        delete gauss;
     8681        xDelete<int>(cs_list);
     8682        xDelete<IssmDouble>(vbasis);
    84718683        return pe;
    84728684}
    … …  
    84848696        IssmDouble      bed_normal[3];
    84858697        IssmDouble  dz[3];
    8486         IssmDouble  basis[6]; //for the six nodes of the penta
    84878698        IssmDouble  Jdet2d;
    84888699        GaussPenta  *gauss=NULL;
    … …  
    84918702        if(!IsOnBed() || !IsFloating()) return NULL;
    84928703        inputs->GetInputValue(&approximation,ApproximationEnum);
     8704        if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     8705
     8706        /*Fetch number of nodes and dof for this finite element*/
     8707        int vnumnodes = this->NumberofNodesVelocity();
     8708        int pnumnodes = this->NumberofNodesPressure();
     8709        int vnumdof   = vnumnodes*NDOF3;
     8710
     8711        /*Prepare coordinate system list*/
     8712        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     8713        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     8714        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     8715
     8716        /*Initialize Element matrix and vectors*/
     8717        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
     8718        IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
     8719
     8720        /*Retrieve all inputs and parameters*/
    84938721        this->parameters->FindParam(&shelf_dampening,DiagnosticShelfDampeningEnum);
    8494         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    8495         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    8496 
    8497         /*Retrieve all inputs and parameters*/
    84988722        rho_water=matpar->GetRhoWater();
    84998723        gravity=matpar->GetG();
    … …  
    85138737
    85148738                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    8515                 GetNodalFunctionsP1(basis, gauss);
     8739                GetNodalFunctionsVelocity(vbasis, gauss);
    85168740
    85178741                BedNormal(&bed_normal[0],xyz_list_tria);
    … …  
    85228746                        vy_input->GetInputValue(&vy, gauss);
    85238747                        vz_input->GetInputValue(&vz, gauss);
    8524                         dt=0;
     8748                        dt=0.;
    85258749                        normal_vel=bed_normal[0]*vx+bed_normal[1]*vy+bed_normal[2]*vz;
    85268750                        damper=gravity*rho_water*pow(1+pow(dz[0],2)+pow(dz[1],2),0.5)*normal_vel*dt;
    85278751                }
    8528                 else damper=0;
     8752                else damper=0.;
    85298753                water_pressure=gravity*rho_water*bed;
    85308754
    8531                 for(i=0;i<NUMVERTICES;i++) for(j=0;j<3;j++) pe->values[i*NDOF4+j]+=(water_pressure+damper)*gauss->weight*Jdet2d*basis[i]*bed_normal[j];
     8755                for(i=0;i<vnumnodes;i++){
     8756                        pe->values[3*i+0]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[0];
     8757                        pe->values[3*i+1]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[1];
     8758                        pe->values[3*i+2]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[2];
     8759                }
    85328760        }
    85338761
    85348762        /*Transform coordinate system*/
    8535         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
     8763        TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
    85368764
    85378765        /*Clean up and return*/
    85388766        delete gauss;
     8767        xDelete<int>(cs_list);
     8768        xDelete<IssmDouble>(vbasis);
    85398769        return pe;
    85408770}
    … …  
    87949024ElementMatrix* Penta::CreateJacobianDiagnosticFS(void){
    87959025
    8796         /*Constants*/
    8797         const int    numdof=NDOF4*NUMVERTICES;
    8798 
    87999026        /*Intermediaries */
    8800         int        i,j;
     9027        int        i,j,approximation;
    88019028        IssmDouble xyz_list[NUMVERTICES][3];
    88029029        IssmDouble Jdet;
    … …  
    88059032        IssmDouble eps3dotdphii,eps3dotdphij;
    88069033        IssmDouble mu_prime;
    8807         IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     9034        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    88089035        IssmDouble eps1[3],eps2[3],eps3[3];
    8809         IssmDouble dphi[3][NUMVERTICES];
    88109036        GaussPenta *gauss=NULL;
     9037
     9038        /*If on water or not FS, skip stiffness: */
     9039        //inputs->GetInputValue(&approximation,ApproximationEnum);
     9040        //if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     9041
     9042        /*Fetch number of nodes and dof for this finite element*/
     9043        int vnumnodes = this->NumberofNodesVelocity();
     9044        int pnumnodes = this->NumberofNodesPressure();
     9045        int numdof    = vnumnodes*NDOF3 + pnumnodes*NDOF1;
     9046
     9047        /*Prepare coordinate system list*/
     9048        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     9049        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     9050        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    88119051
    88129052        /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
    88139053        ElementMatrix* Ke=CreateKMatrixDiagnosticFS();
     9054        IssmDouble*    dbasis = xNew<IssmDouble>(3*vnumnodes);
    88149055
    88159056        /*Retrieve all inputs and parameters*/
    88169057        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8817         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    8818         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    8819         Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
     9058        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     9059        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     9060        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    88209061
    88219062        /* Start  looping on the number of gaussian points: */
    … …  
    88259066                gauss->GaussPoint(ig);
    88269067
    8827                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8828                 GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    8829 
     9068                GetJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     9069                GetNodalFunctionsDerivativesVelocity(dbasis,&xyz_list[0][0],gauss);
     9070
     9071                //this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     9072                //material->GetViscosityDerivativeEpsSquareFS(&mu_prime,&epsilon[0]);
     9073                //eps1[0]=epsilon[0];   eps2[0]=epsilon[3];   eps3[0]=epsilon[4];
     9074                //eps1[1]=epsilon[3];   eps2[1]=epsilon[1];   eps3[1]=epsilon[5];
     9075                //eps1[2]=epsilon[4];   eps2[2]=epsilon[5];   eps3[2]=epsilon[2];
    88309076                this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    88319077                material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    … …  
    88349080                eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
    88359081
    8836                 for(i=0;i<6;i++){
    8837                         for(j=0;j<6;j++){
    8838                                 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
    8839                                 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
    8840                                 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
    8841                                 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
    8842                                 eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i];
    8843                                 eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j];
    8844 
    8845                                 Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
    8846                                 Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
    8847                                 Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
    8848 
    8849                                 Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
    8850                                 Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
    8851                                 Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
    8852 
    8853                                 Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
    8854                                 Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
    8855                                 Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
     9082                for(i=0;i<vnumnodes-1;i++){
     9083                        for(j=0;j<vnumnodes-1;j++){
     9084                                eps1dotdphii=eps1[0]*dbasis[0*vnumnodes+i]+eps1[1]*dbasis[1*vnumnodes+i]+eps1[2]*dbasis[2*vnumnodes+i];
     9085                                eps1dotdphij=eps1[0]*dbasis[0*vnumnodes+j]+eps1[1]*dbasis[1*vnumnodes+j]+eps1[2]*dbasis[2*vnumnodes+j];
     9086                                eps2dotdphii=eps2[0]*dbasis[0*vnumnodes+i]+eps2[1]*dbasis[1*vnumnodes+i]+eps2[2]*dbasis[2*vnumnodes+i];
     9087                                eps2dotdphij=eps2[0]*dbasis[0*vnumnodes+j]+eps2[1]*dbasis[1*vnumnodes+j]+eps2[2]*dbasis[2*vnumnodes+j];
     9088                                eps3dotdphii=eps3[0]*dbasis[0*vnumnodes+i]+eps3[1]*dbasis[1*vnumnodes+i]+eps3[2]*dbasis[2*vnumnodes+i];
     9089                                eps3dotdphij=eps3[0]*dbasis[0*vnumnodes+j]+eps3[1]*dbasis[1*vnumnodes+j]+eps3[2]*dbasis[2*vnumnodes+j];
     9090
     9091                                Ke->values[numdof*(3*i+0)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
     9092                                Ke->values[numdof*(3*i+0)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
     9093                                Ke->values[numdof*(3*i+0)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
     9094
     9095                                Ke->values[numdof*(3*i+1)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
     9096                                Ke->values[numdof*(3*i+1)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
     9097                                Ke->values[numdof*(3*i+1)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
     9098
     9099                                Ke->values[numdof*(3*i+2)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
     9100                                Ke->values[numdof*(3*i+2)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
     9101                                Ke->values[numdof*(3*i+2)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
    88569102                        }
    88579103                }
    … …  
    88599105
    88609106        /*Transform Coordinate System*/
    8861         TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
     9107        TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);
    88629108
    88639109        /*Clean up and return*/
     9110        xDelete<int>(cs_list);
     9111        xDelete<IssmDouble>(dbasis);
    88649112        delete gauss;
    88659113        return Ke;
    … …  
    88699117void  Penta::GetSolutionFromInputsDiagnosticHoriz(Vector<IssmDouble>* solution){
    88709118
    8871         const int    numdof=NDOF2*NUMVERTICES;
    8872 
    8873         int         i;
    88749119        int         approximation;
    8875         int        *doflist        = NULL;
     9120        int        *doflist = NULL;
    88769121        IssmDouble  vx,vy;
    8877         IssmDouble  values[numdof];
    8878         GaussPenta *gauss;
     9122
     9123        /*Fetch number of nodes and dof for this finite element*/
     9124        int numnodes = this->NumberofNodes();
     9125        int numdof   = numnodes*NDOF2;
    88799126
    88809127        /*Get approximation enum and dof list: */
    … …  
    88839130        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    88849131
    8885         /*If the element is a coupling, do nothing: every node is also on an other elements
    8886          * (as coupling is between SSA and HO) so the other element will take care of it*/
     9132        /*Fetch dof list and allocate solution vectors*/
    88879133        GetDofList(&doflist,approximation,GsetEnum);
     9134        IssmDouble* values = xNew<IssmDouble>(numdof);
    88889135
    88899136        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    8890         /*P1 element only for now*/
    8891         gauss=new GaussPenta();
    8892         for(i=0;i<NUMVERTICES;i++){
     9137        GaussPenta* gauss=new GaussPenta();
     9138        for(int i=0;i<numnodes;i++){
     9139                gauss->GaussNode(numnodes,i);
    88939140
    88949141                /*Recover vx and vy*/
    8895                 gauss->GaussVertex(i);
    88969142                vx_input->GetInputValue(&vx,gauss);
    88979143                vy_input->GetInputValue(&vy,gauss);
    … …  
    89059151        /*Free ressources:*/
    89069152        delete gauss;
     9153        xDelete<IssmDouble>(values);
    89079154        xDelete<int>(doflist);
    89089155}
    … …  
    89809227void  Penta::GetSolutionFromInputsDiagnosticFS(Vector<IssmDouble>* solution){
    89819228
    8982         const int    numdof=NDOF4*NUMVERTICES;
    8983 
    8984         int          i;
    8985         int*         doflist=NULL;
    8986         IssmDouble       vx,vy,vz,p;
    8987         IssmDouble       FSreconditioning;
    8988         IssmDouble       values[numdof];
    8989         GaussPenta   *gauss;
     9229        int*         vdoflist=NULL;
     9230        int*         pdoflist=NULL;
     9231        IssmDouble   vx,vy,vz,p;
     9232        IssmDouble   FSreconditioning;
     9233        GaussPenta  *gauss;
     9234
     9235        /*Fetch number of nodes and dof for this finite element*/
     9236        int vnumnodes = this->NumberofNodesVelocityFinal();
     9237        int pnumnodes = this->NumberofNodesPressure();
     9238        int vnumdof   = vnumnodes*NDOF3;
     9239        int pnumdof   = pnumnodes*NDOF1;
     9240
     9241        /*Initialize values*/
     9242        IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
     9243        IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
    89909244
    89919245        /*Get dof list: */
    8992         GetDofList(&doflist,FSApproximationEnum,GsetEnum);
     9246        GetDofListVelocity(&vdoflist,GsetEnum);
     9247        GetDofListPressure(&pdoflist,GsetEnum);
    89939248        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    89949249        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    … …  
    89969251        Input* p_input =inputs->GetInput(PressureEnum); _assert_(p_input);
    89979252
    8998         /*Recondition pressure: */
    89999253        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    90009254
    9001         /*Ok, we have vx vy vz and P in values, fill in vx vy vz P arrays: */
    9002         /*P1 element only for now*/
    9003         gauss=new GaussPenta();
    9004         for(i=0;i<NUMVERTICES;i++){
    9005                 gauss->GaussVertex(i);
     9255        /*Ok, we have vx vy vz in values, fill in vx vy vz arrays: */
     9256        gauss = new GaussPenta();
     9257        for(int i=0;i<vnumnodes;i++){
     9258                gauss->GaussNode(vnumnodes,i);
     9259
    90069260                vx_input->GetInputValue(&vx,gauss);
    90079261                vy_input->GetInputValue(&vy,gauss);
    90089262                vz_input->GetInputValue(&vz,gauss);
     9263                vvalues[i*NDOF3+0]=vx;
     9264                vvalues[i*NDOF3+1]=vy;
     9265                vvalues[i*NDOF3+2]=vz;
     9266        }
     9267        for(int i=0;i<pnumnodes;i++){
     9268                gauss->GaussNode(pnumnodes,i);
     9269
    90099270                p_input ->GetInputValue(&p ,gauss);
    9010                 values[i*NDOF4+0]=vx;
    9011                 values[i*NDOF4+1]=vy;
    9012                 values[i*NDOF4+2]=vz;
    9013                 values[i*NDOF4+3]=p/FSreconditioning;
     9271                pvalues[i]=p/FSreconditioning;
    90149272        }
    90159273
    90169274        /*Add value to global vector*/
    9017         solution->SetValues(numdof,doflist,values,INS_VAL);
     9275        solution->SetValues(vnumdof,vdoflist,vvalues,INS_VAL);
     9276        solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL);
    90189277
    90199278        /*Free ressources:*/
    90209279        delete gauss;
    9021         xDelete<int>(doflist);
     9280        xDelete<int>(pdoflist);
     9281        xDelete<int>(vdoflist);
     9282        xDelete<IssmDouble>(pvalues);
     9283        xDelete<IssmDouble>(vvalues);
    90229284}
    90239285/*}}}*/
    … …  
    93439605        /*Transform solution in Cartesian Space*/
    93449606        TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
    9345         TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZPEnum);
     9607        TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZEnum);
    93469608
    93479609        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    95919853        /*Transform solution in Cartesian Space*/
    95929854        TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
    9593         TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZPEnum);
     9855        TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZEnum);
    95949856
    95959857        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    981410076void  Penta::InputUpdateFromSolutionDiagnosticFS(IssmDouble* solution){
    981510077
    9816         const int numdof=NDOF4*NUMVERTICES;
    9817 
    9818         int     i;
    9819         IssmDouble  values[numdof];
    9820         IssmDouble  vx[NUMVERTICES];
    9821         IssmDouble  vy[NUMVERTICES];
    9822         IssmDouble  vz[NUMVERTICES];
    9823         IssmDouble  vel[NUMVERTICES];
    9824         IssmDouble  pressure[NUMVERTICES];
    9825         IssmDouble  FSreconditioning;
    9826         int*    doflist=NULL;
     10078        int          i;
     10079        int*         vdoflist=NULL;
     10080        int*         pdoflist=NULL;
     10081        IssmDouble   FSreconditioning;
     10082        GaussPenta  *gauss;
     10083
     10084        /*Fetch number of nodes and dof for this finite element*/
     10085        int vnumnodes = this->NumberofNodesVelocityFinal();
     10086        int pnumnodes = this->NumberofNodesPressure();
     10087        int vnumdof   = vnumnodes*NDOF3;
     10088        int pnumdof   = pnumnodes*NDOF1;
     10089
     10090        /*Initialize values*/
     10091        IssmDouble* vvalues  = xNew<IssmDouble>(vnumdof);
     10092        IssmDouble* pvalues  = xNew<IssmDouble>(pnumdof);
     10093        IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
     10094        IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
     10095        IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
     10096        IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
     10097        IssmDouble* pressure = xNew<IssmDouble>(vnumnodes);
    982710098
    982810099        /*Get dof list: */
    9829         GetDofList(&doflist,FSApproximationEnum,GsetEnum);
     10100        GetDofListVelocity(&vdoflist,GsetEnum);
     10101        GetDofListPressure(&pdoflist,GsetEnum);
    983010102
    983110103        /*Use the dof list to index into the solution vector: */
    9832         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
     10104        for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
     10105        for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
    983310106
    983410107        /*Transform solution in Cartesian Space*/
    9835         TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYZPEnum);
     10108        TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYZEnum);
    983610109
    983710110        /*Ok, we have vx and vy in values, fill in all arrays: */
    9838         for(i=0;i<NUMVERTICES;i++){
    9839                 vx[i]=values[i*NDOF4+0];
    9840                 vy[i]=values[i*NDOF4+1];
    9841                 vz[i]=values[i*NDOF4+2];
    9842                 pressure[i]=values[i*NDOF4+3];
    9843 
    9844                 /*Check solution*/
    9845                 if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
    9846                 if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
    9847                 if(xIsNan<IssmDouble>(vz[i]))       _error_("NaN found in solution vector");
     10111        for(i=0;i<vnumnodes;i++){
     10112                vx[i] = vvalues[i*NDOF3+0];
     10113                vy[i] = vvalues[i*NDOF3+1];
     10114                vz[i] = vvalues[i*NDOF3+2];
     10115                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
     10116                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     10117                if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
     10118        }
     10119        for(i=0;i<pnumnodes;i++){
     10120                pressure[i] = pvalues[i];
    984810121                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    984910122        }
    … …  
    985110124        /*Recondition pressure and compute vel: */
    985210125        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    9853         for(i=0;i<NUMVERTICES;i++) pressure[i]=pressure[i]*FSreconditioning;
    9854         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);
     10126        for(i=0;i<pnumnodes;i++) pressure[i]=pressure[i]*FSreconditioning;
     10127        for(i=0;i<vnumnodes;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    985510128
    985610129        /*Now, we have to move the previous inputs  to old
    … …  
    986910142
    987010143        /*Free ressources:*/
    9871         xDelete<int>(doflist);
     10144        xDelete<IssmDouble>(pressure);
     10145        xDelete<IssmDouble>(vel);
     10146        xDelete<IssmDouble>(vz);
     10147        xDelete<IssmDouble>(vy);
     10148        xDelete<IssmDouble>(vx);
     10149        xDelete<IssmDouble>(vvalues);
     10150        xDelete<IssmDouble>(pvalues);
     10151        xDelete<int>(vdoflist);
     10152        xDelete<int>(pdoflist);
    987210153}
    987310154/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r15643 r15654  
    184184                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints);
    185185                void             GetDofList(int** pdoflist,int approximation_enum,int setenum);
     186                void             GetDofListVelocity(int** pdoflist,int setenum);
     187                void             GetDofListPressure(int** pdoflist,int setenum);
    186188                void             GetVertexPidList(int* doflist);
    187189                void           GetVertexSidList(int* sidlist);
    … …  
    252254                ElementMatrix* CreateKMatrixDiagnosticFS(void);
    253255                ElementMatrix* CreateKMatrixDiagnosticFSViscous(void);
    254                 ElementMatrix* CreateKMatrixDiagnosticFSGLSViscous(void);
     256                void           KMatrixGLSstabilization(ElementMatrix* Ke);
    255257                ElementMatrix* CreateKMatrixDiagnosticFSFriction(void);
    256                 ElementMatrix* CreateKMatrixDiagnosticFSGLSFriction(void);
    257258                ElementMatrix* CreateKMatrixDiagnosticVert(void);
    258259                ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
    … …  
    295296                ElementVector* CreatePVectorDiagnosticFSFront(void);
    296297                ElementVector* CreatePVectorDiagnosticFSViscous(void);
    297                 ElementVector* CreatePVectorDiagnosticFSGLSViscous(void);
     298                void           PVectorGLSstabilization(ElementVector* pe);
    298299                ElementVector* CreatePVectorDiagnosticFSShelf(void);
    299300                ElementVector* CreatePVectorDiagnosticVert(void);

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

    r15643 r15654  
    250250}
    251251/*}}}*/
     252/*FUNCTION PentaRef::GetBFSstrainrate {{{*/
     253void PentaRef::GetBFSstrainrate(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
     254
     255        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4.
     256         * For node i, Bi can be expressed in the actual coordinate system
     257         * by:          Bi=[ dh/dx          0              0     ]
     258         *                                      [   0           dh/dy           0     ]
     259         *                                      [   0             0           dh/dy   ]
     260         *                                      [ 1/2*dh/dy    1/2*dh/dx        0     ]
     261         *                                      [ 1/2*dh/dz       0         1/2*dh/dx ]
     262         *                                      [   0          1/2*dh/dz    1/2*dh/dy ]
     263         *      where h is the interpolation function for node i.
     264         *      Same thing for Bb except the last column that does not exist.
     265         */
     266
     267        /*Fetch number of nodes for this finite element*/
     268        int numnodes = this->NumberofNodes();
     269
     270        /*Get nodal functions derivatives*/
     271        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
     272        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     273
     274        /*Build B: */
     275        for(int i=0;i<numnodes;i++){
     276                B[3*numnodes*0+3*i+0] = dbasis[0*numnodes+i+0];
     277                B[3*numnodes*0+3*i+1] = 0.;
     278                B[3*numnodes*0+3*i+2] = 0.;
     279
     280                B[3*numnodes*1+3*i+0] = 0.;
     281                B[3*numnodes*1+3*i+1] = dbasis[1*numnodes+i+0];
     282                B[3*numnodes*1+3*i+2] = 0.;
     283
     284                B[3*numnodes*2+3*i+0] = 0.;
     285                B[3*numnodes*2+3*i+1] = 0.;
     286                B[3*numnodes*2+3*i+2] = dbasis[2*numnodes+i+0];
     287
     288                B[3*numnodes*3+3*i+0] = .5*dbasis[1*numnodes+i+0];
     289                B[3*numnodes*3+3*i+1] = .5*dbasis[0*numnodes+i+0];
     290                B[3*numnodes*3+3*i+2] = 0.;
     291
     292                B[3*numnodes*4+3*i+0] = .5*dbasis[2*numnodes+i+0];
     293                B[3*numnodes*4+3*i+1] = 0.;
     294                B[3*numnodes*4+3*i+2] = .5*dbasis[0*numnodes+i+0];
     295
     296                B[3*numnodes*5+3*i+0] = 0.;
     297                B[3*numnodes*5+3*i+1] = .5*dbasis[2*numnodes+i+0];
     298                B[3*numnodes*5+3*i+2] = .5*dbasis[1*numnodes+i+0];
     299        }
     300
     301        /*Clean up*/
     302        xDelete<IssmDouble>(dbasis);
     303}
     304/*}}}*/
    252305/*FUNCTION PentaRef::GetBFS {{{*/
    253306void PentaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
    254307
    255         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4.
    256          * For node i, Bi can be expressed in the actual coordinate system
    257          * by:          Bi=[ dh/dx          0              0       0  ]
    258          *                                      [   0           dh/dy           0       0  ]
    259          *                                      [   0             0           dh/dy     0  ]
    260          *                                      [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
    261          *                                      [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
    262          *                                      [   0          1/2*dh/dz    1/2*dh/dy   0  ]
    263          *                                      [   0             0             0       h  ]
    264          *                                      [ dh/dx         dh/dy         dh/dz     0  ]
     308        /*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
     309         * For node i, Bvi can be expressed in the actual coordinate system
     310         * by:     Bvi=[ dh/dx          0              0      ]
     311         *                                      [   0           dh/dy           0      ]
     312         *                                      [   0             0           dh/dy    ]
     313         *                                      [ 1/2*dh/dy    1/2*dh/dx        0      ]
     314         *                                      [ 1/2*dh/dz       0         1/2*dh/dx  ]
     315         *                                      [   0          1/2*dh/dz    1/2*dh/dy  ]
     316         *                                      [   0             0             0      ]
     317         *                                      [ dh/dx         dh/dy         dh/dz    ]
     318         *
     319         * by:    Bpi=[ 0 ]
     320         *                                      [ 0 ]
     321         *                                      [ 0 ]
     322         *                                      [ 0 ]
     323         *                                      [ 0 ]
     324         *                                      [ 0 ]
     325         *                                      [ h ]
     326         *                                      [ 0 ]
    265327         *      where h is the interpolation function for node i.
    266328         *      Same thing for Bb except the last column that does not exist.
    267329         */
    268330
    269         int i;
    270 
    271         IssmDouble dbasismini[3][NUMNODESP1b];
    272         IssmDouble basis[NUMNODESP1];
    273 
    274         /*Get dbasismini in actual coordinate system: */
    275         GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
    276         GetNodalFunctionsP1(basis, gauss);
     331        /*Fetch number of nodes for this finite element*/
     332        int pnumnodes = this->NumberofNodesPressure();
     333        int vnumnodes = this->NumberofNodesVelocity();
     334
     335        /*Get nodal functions derivatives*/
     336        IssmDouble* vdbasis=xNew<IssmDouble>(3*vnumnodes);
     337        IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
     338        GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     339        GetNodalFunctionsPressure(pbasis,gauss);
    277340
    278341        /*Build B: */
    279         for(i=0;i<NUMNODESP1b;i++){
    280                 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i+0]; //B[0][NDOF4*i+0] = dbasis[0][i+0];
    281                 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
    282                 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
    283                 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
    284                 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i+0];
    285                 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
    286                 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.;
    287                 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.;
    288                 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dbasismini[2][i+0];
    289                 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = .5*dbasismini[1][i+0];
    290                 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = .5*dbasismini[0][i+0];
    291                 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
    292                 B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = .5*dbasismini[2][i+0];
    293                 B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1] = 0.;
    294                 B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = .5*dbasismini[0][i+0];
    295                 B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+0] = 0.;
    296                 B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = .5*dbasismini[2][i+0];
    297                 B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = .5*dbasismini[1][i+0];
    298                 B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+0] = 0.;
    299                 B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1] = 0.;
    300                 B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2] = 0.;
    301                 B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = dbasismini[0][i+0];
    302                 B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = dbasismini[1][i+0];
    303                 B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = dbasismini[2][i+0];
    304         }
    305 
    306         for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
    307                 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
    308                 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
    309                 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
    310                 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = 0.;
    311                 B[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3] = 0.;
    312                 B[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3] = 0.;
    313                 B[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3] = basis[i];
    314                 B[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3] = 0.;
    315         }
    316 
     342        for(int i=0;i<vnumnodes;i++){
     343                B[(3*vnumnodes+pnumnodes)*0+3*i+0] = vdbasis[0*vnumnodes+i];
     344                B[(3*vnumnodes+pnumnodes)*0+3*i+1] = 0.;
     345                B[(3*vnumnodes+pnumnodes)*0+3*i+2] = 0.;
     346                B[(3*vnumnodes+pnumnodes)*1+3*i+0] = 0.;
     347                B[(3*vnumnodes+pnumnodes)*1+3*i+1] = vdbasis[1*vnumnodes+i];
     348                B[(3*vnumnodes+pnumnodes)*1+3*i+2] = 0.;
     349                B[(3*vnumnodes+pnumnodes)*2+3*i+0] = 0.;
     350                B[(3*vnumnodes+pnumnodes)*2+3*i+1] = 0.;
     351                B[(3*vnumnodes+pnumnodes)*2+3*i+2] = vdbasis[2*vnumnodes+i];
     352                B[(3*vnumnodes+pnumnodes)*3+3*i+0] = .5*vdbasis[1*vnumnodes+i];
     353                B[(3*vnumnodes+pnumnodes)*3+3*i+1] = .5*vdbasis[0*vnumnodes+i];
     354                B[(3*vnumnodes+pnumnodes)*3+3*i+2] = 0.;
     355                B[(3*vnumnodes+pnumnodes)*4+3*i+0] = .5*vdbasis[2*vnumnodes+i];
     356                B[(3*vnumnodes+pnumnodes)*4+3*i+1] = 0.;
     357                B[(3*vnumnodes+pnumnodes)*4+3*i+2] = .5*vdbasis[0*vnumnodes+i];
     358                B[(3*vnumnodes+pnumnodes)*5+3*i+0] = 0.;
     359                B[(3*vnumnodes+pnumnodes)*5+3*i+1] = .5*vdbasis[2*vnumnodes+i];
     360                B[(3*vnumnodes+pnumnodes)*5+3*i+2] = .5*vdbasis[1*vnumnodes+i];
     361                B[(3*vnumnodes+pnumnodes)*6+3*i+0] = 0.;
     362                B[(3*vnumnodes+pnumnodes)*6+3*i+1] = 0.;
     363                B[(3*vnumnodes+pnumnodes)*6+3*i+2] = 0.;
     364                B[(3*vnumnodes+pnumnodes)*7+3*i+0] = vdbasis[0*vnumnodes+i];
     365                B[(3*vnumnodes+pnumnodes)*7+3*i+1] = vdbasis[1*vnumnodes+i];
     366                B[(3*vnumnodes+pnumnodes)*7+3*i+2] = vdbasis[2*vnumnodes+i];
     367        }
     368        for(int i=0;i<pnumnodes;i++){
     369                B[(3*vnumnodes+pnumnodes)*0+(3*vnumnodes)+i] = 0.;
     370                B[(3*vnumnodes+pnumnodes)*1+(3*vnumnodes)+i] = 0.;
     371                B[(3*vnumnodes+pnumnodes)*2+(3*vnumnodes)+i] = 0.;
     372                B[(3*vnumnodes+pnumnodes)*3+(3*vnumnodes)+i] = 0.;
     373                B[(3*vnumnodes+pnumnodes)*4+(3*vnumnodes)+i] = 0.;
     374                B[(3*vnumnodes+pnumnodes)*5+(3*vnumnodes)+i] = 0.;
     375                B[(3*vnumnodes+pnumnodes)*6+(3*vnumnodes)+i] = pbasis[i];
     376                B[(3*vnumnodes+pnumnodes)*7+(3*vnumnodes)+i] = 0.;
     377        }
     378
     379        /*Clean up*/
     380        xDelete<IssmDouble>(vdbasis);
     381        xDelete<IssmDouble>(pbasis);
    317382}
    318383/*}}}*/
    … …  
    388453         *      For node i, Bi' can be expressed in the actual coordinate system
    389454         *      by:
    390          *                              Bi'=[  dh/dx   0          0       0]
    391          *                                       [   0      dh/dy      0       0]
    392          *                                       [   0      0         dh/dz    0]
    393          *                                       [  dh/dy   dh/dx      0       0]
    394          *                                       [  dh/dz   0        dh/dx     0]
    395          *                                       [   0      dh/dz    dh/dy     0]
    396          *                                       [  dh/dx   dh/dy    dh/dz     0]
    397          *                                       [   0      0          0       h]
     455         *                      Bvi' = [  dh/dx   0          0    ]
     456         *                                       [   0      dh/dy      0    ]
     457         *                                       [   0      0         dh/dz ]
     458         *                                       [  dh/dy   dh/dx      0    ]
     459         *                                       [  dh/dz   0        dh/dx  ]
     460         *                                       [   0      dh/dz    dh/dy  ]
     461         *                                       [  dh/dx   dh/dy    dh/dz  ]
     462         *                                       [   0      0          0    ]
     463         *
     464         * by:    Bpi=[ 0 ]
     465         *                                      [ 0 ]
     466         *                                      [ 0 ]
     467         *                                      [ 0 ]
     468         *                                      [ 0 ]
     469         *                                      [ 0 ]
     470         *                                      [ 0 ]
     471         *                                      [ h ]
    398472         *      where h is the interpolation function for node i.
    399          *
    400          *      Same thing for the bubble fonction except that there is no fourth column
    401          */
    402 
    403         int i;
    404         IssmDouble dbasismini[3][NUMNODESP1b];
    405         IssmDouble basis[NUMNODESP1];
    406 
    407         /*Get dbasismini in actual coordinate system: */
    408         GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
    409         GetNodalFunctionsP1(basis, gauss);
    410 
    411         /*B_primeuild B_prime: */
    412         for(i=0;i<NUMNODESP1b;i++){
    413                 B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i];
    414                 B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
    415                 B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
    416                 B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
    417                 B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i];
    418                 B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
    419                 B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.;
    420                 B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.;
    421                 B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = dbasismini[2][i];
    422                 B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = dbasismini[1][i];
    423                 B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = dbasismini[0][i];
    424                 B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
    425                 B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+0] = dbasismini[2][i];
    426                 B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+1] = 0.;
    427                 B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+2] = dbasismini[0][i];
    428                 B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+0] = 0.;
    429                 B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+1] = dbasismini[2][i];
    430                 B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+2] = dbasismini[1][i];
    431                 B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+0] = dbasismini[0][i];
    432                 B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+1] = dbasismini[1][i];
    433                 B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+2] = dbasismini[2][i];
    434                 B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+0] = 0.;
    435                 B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+1] = 0.;
    436                 B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+2] = 0.;
    437         }
    438 
    439         for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
    440                 B_prime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
    441                 B_prime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
    442                 B_prime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
    443                 B_prime[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = 0.;
    444                 B_prime[(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3] = 0.;
    445                 B_prime[(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3] = 0.;
    446                 B_prime[(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3] = 0.;
    447                 B_prime[(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3] = basis[i];
    448         }
    449 
     473         */
     474
     475        /*Fetch number of nodes for this finite element*/
     476        int pnumnodes = this->NumberofNodesPressure();
     477        int vnumnodes = this->NumberofNodesVelocity();
     478
     479        /*Get nodal functions derivatives*/
     480        IssmDouble* vdbasis=xNew<IssmDouble>(3*vnumnodes);
     481        IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
     482        GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     483        GetNodalFunctionsPressure(pbasis,gauss);
     484
     485        /*Build B_prime: */
     486        for(int i=0;i<vnumnodes;i++){
     487                B_prime[(3*vnumnodes+pnumnodes)*0+3*i+0] = vdbasis[0*vnumnodes+i];
     488                B_prime[(3*vnumnodes+pnumnodes)*0+3*i+1] = 0.;
     489                B_prime[(3*vnumnodes+pnumnodes)*0+3*i+2] = 0.;
     490                B_prime[(3*vnumnodes+pnumnodes)*1+3*i+0] = 0.;
     491                B_prime[(3*vnumnodes+pnumnodes)*1+3*i+1] = vdbasis[1*vnumnodes+i];
     492                B_prime[(3*vnumnodes+pnumnodes)*1+3*i+2] = 0.;
     493                B_prime[(3*vnumnodes+pnumnodes)*2+3*i+0] = 0.;
     494                B_prime[(3*vnumnodes+pnumnodes)*2+3*i+1] = 0.;
     495                B_prime[(3*vnumnodes+pnumnodes)*2+3*i+2] = vdbasis[2*vnumnodes+i];
     496                B_prime[(3*vnumnodes+pnumnodes)*3+3*i+0] = vdbasis[1*vnumnodes+i];
     497                B_prime[(3*vnumnodes+pnumnodes)*3+3*i+1] = vdbasis[0*vnumnodes+i];
     498                B_prime[(3*vnumnodes+pnumnodes)*3+3*i+2] = 0.;
     499                B_prime[(3*vnumnodes+pnumnodes)*4+3*i+0] = vdbasis[2*vnumnodes+i];
     500                B_prime[(3*vnumnodes+pnumnodes)*4+3*i+1] = 0.;
     501                B_prime[(3*vnumnodes+pnumnodes)*4+3*i+2] = vdbasis[0*vnumnodes+i];
     502                B_prime[(3*vnumnodes+pnumnodes)*5+3*i+0] = 0.;
     503                B_prime[(3*vnumnodes+pnumnodes)*5+3*i+1] = vdbasis[2*vnumnodes+i];
     504                B_prime[(3*vnumnodes+pnumnodes)*5+3*i+2] = vdbasis[1*vnumnodes+i];
     505                B_prime[(3*vnumnodes+pnumnodes)*6+3*i+0] = vdbasis[0*vnumnodes+i];
     506                B_prime[(3*vnumnodes+pnumnodes)*6+3*i+1] = vdbasis[1*vnumnodes+i];
     507                B_prime[(3*vnumnodes+pnumnodes)*6+3*i+2] = vdbasis[2*vnumnodes+i];
     508                B_prime[(3*vnumnodes+pnumnodes)*7+3*i+0] = 0.;
     509                B_prime[(3*vnumnodes+pnumnodes)*7+3*i+1] = 0.;
     510                B_prime[(3*vnumnodes+pnumnodes)*7+3*i+2] = 0.;
     511        }
     512        for(int i=0;i<pnumnodes;i++){
     513                B_prime[(3*vnumnodes+pnumnodes)*0+(3*vnumnodes)+i] = 0.;
     514                B_prime[(3*vnumnodes+pnumnodes)*1+(3*vnumnodes)+i] = 0.;
     515                B_prime[(3*vnumnodes+pnumnodes)*2+(3*vnumnodes)+i] = 0.;
     516                B_prime[(3*vnumnodes+pnumnodes)*3+(3*vnumnodes)+i] = 0.;
     517                B_prime[(3*vnumnodes+pnumnodes)*4+(3*vnumnodes)+i] = 0.;
     518                B_prime[(3*vnumnodes+pnumnodes)*5+(3*vnumnodes)+i] = 0.;
     519                B_prime[(3*vnumnodes+pnumnodes)*6+(3*vnumnodes)+i] = 0.;
     520                B_prime[(3*vnumnodes+pnumnodes)*7+(3*vnumnodes)+i] = pbasis[i];
     521        }
     522
     523        /*Clean up*/
     524        xDelete<IssmDouble>(vdbasis);
     525        xDelete<IssmDouble>(pbasis);
    450526}
    451527/*}}}*/
    … …  
    673749         * For node i, Li can be expressed in the actual coordinate system
    674750         * by:
    675          *       Li=[ h 0 ]
    676          *                    [ 0 h ]
    677          *                    [ 0 0 ]
    678          *                    [ 0 0 ]
     751         *       Li=[ h 0 0 0 ]
     752         *                    [ 0 h 0 0 ]
    679753         * where h is the interpolation function for node i.
    680754         */
    681755
    682         const int num_dof=4;
    683         IssmDouble L1L2l3[NUMNODESP1_2d];
    684 
    685         /*Get L1L2l3 in actual coordinate system: */
    686         L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    687         L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    688         L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     756        /*Fetch number of nodes for this finite element*/
     757        int pnumnodes = this->NumberofNodesPressure();
     758        int vnumnodes = this->NumberofNodesVelocity();
     759        int pnumdof   = pnumnodes;
     760        int vnumdof   = vnumnodes*NDOF3;
     761
     762        /*Get nodal functions derivatives*/
     763        IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
     764        GetNodalFunctionsVelocity(vbasis,gauss);
    689765
    690766        /*Build LFS: */
    691         for(int i=0;i<3;i++){
    692                 LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
    693                 LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
    694                 LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
    695                 LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
    696 
    697                 LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
    698                 LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
    699                 LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
    700                 LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
    701         }
     767        for(int i=0;i<vnumnodes;i++){
     768                LFS[(vnumdof+pnumdof)*0+3*i+0] = vbasis[i];
     769                LFS[(vnumdof+pnumdof)*0+3*i+1] = 0.;
     770                LFS[(vnumdof+pnumdof)*0+3*i+2] = 0.;
     771
     772                LFS[(vnumdof+pnumdof)*1+3*i+0] = 0.;
     773                LFS[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];
     774                LFS[(vnumdof+pnumdof)*1+3*i+2] = 0.;
     775        }
     776
     777        for(int i=0;i<pnumnodes;i++){
     778                LFS[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
     779                LFS[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
     780        }
     781
     782        /*Clean-up*/
     783        xDelete<IssmDouble>(vbasis);
    702784}
    703785/*}}}*/
    … …  
    11961278}
    11971279/*}}}*/
     1280/*FUNCTION PentaRef::GetNodalFunctionsVelocity{{{*/
     1281void PentaRef::GetNodalFunctionsVelocity(IssmDouble* basis,GaussPenta* gauss){
     1282        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     1283
     1284        switch(this->element_type){
     1285                case P1P1Enum:
     1286                        this->element_type = P1Enum;
     1287                        this->GetNodalFunctions(basis,gauss);
     1288                        this->element_type = P1P1Enum;
     1289                        return;
     1290                case P1P1GLSEnum:
     1291                        this->element_type = P1Enum;
     1292                        this->GetNodalFunctions(basis,gauss);
     1293                        this->element_type = P1P1GLSEnum;
     1294                        return;
     1295                case MINIcondensedEnum:
     1296                        this->element_type = P1bubbleEnum;
     1297                        this->GetNodalFunctions(basis,gauss);
     1298                        this->element_type = MINIcondensedEnum;
     1299                        return;
     1300                case MINIEnum:
     1301                        this->element_type = P1bubbleEnum;
     1302                        this->GetNodalFunctions(basis,gauss);
     1303                        this->element_type = MINIEnum;
     1304                        return;
     1305                default:
     1306                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1307        }
     1308}
     1309/*}}}*/
     1310/*FUNCTION PentaRef::GetNodalFunctionsPressure{{{*/
     1311void PentaRef::GetNodalFunctionsPressure(IssmDouble* basis,GaussPenta* gauss){
     1312        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     1313
     1314        switch(this->element_type){
     1315                case P1P1Enum:
     1316                        this->element_type = P1Enum;
     1317                        this->GetNodalFunctions(basis,gauss);
     1318                        this->element_type = P1P1Enum;
     1319                        return;
     1320                case P1P1GLSEnum:
     1321                        this->element_type = P1Enum;
     1322                        this->GetNodalFunctions(basis,gauss);
     1323                        this->element_type = P1P1GLSEnum;
     1324                        return;
     1325                case MINIcondensedEnum:
     1326                        this->element_type = P1Enum;
     1327                        this->GetNodalFunctions(basis,gauss);
     1328                        this->element_type = MINIcondensedEnum;
     1329                        return;
     1330                case MINIEnum:
     1331                        this->element_type = P1Enum;
     1332                        this->GetNodalFunctions(basis,gauss);
     1333                        this->element_type = MINIEnum;
     1334                        return;
     1335                default:
     1336                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1337        }
     1338}
     1339/*}}}*/
    11981340/*FUNCTION PentaRef::GetNodalFunctionsDerivatives{{{*/
    11991341void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){
    … …  
    12291371        xDelete<IssmDouble>(dbasis_ref);
    12301372
     1373}
     1374/*}}}*/
     1375/*FUNCTION PentaRef::GetNodalFunctionsDerivativesVelocity{{{*/
     1376void PentaRef::GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){
     1377        switch(this->element_type){
     1378                case P1P1Enum:
     1379                        this->element_type = P1Enum;
     1380                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1381                        this->element_type = P1P1Enum;
     1382                        return;
     1383                case P1P1GLSEnum:
     1384                        this->element_type = P1Enum;
     1385                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1386                        this->element_type = P1P1GLSEnum;
     1387                        return;
     1388                case MINIcondensedEnum:
     1389                        this->element_type = P1bubbleEnum;
     1390                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1391                        this->element_type = MINIcondensedEnum;
     1392                        return;
     1393                case MINIEnum:
     1394                        this->element_type = P1bubbleEnum;
     1395                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1396                        this->element_type = MINIEnum;
     1397                        return;
     1398                default:
     1399                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1400        }
     1401}
     1402/*}}}*/
     1403/*FUNCTION PentaRef::GetNodalFunctionsDerivativesPressure{{{*/
     1404void PentaRef::GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){
     1405        switch(this->element_type){
     1406                case P1P1Enum:
     1407                        this->element_type = P1Enum;
     1408                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1409                        this->element_type = P1P1Enum;
     1410                        return;
     1411                case P1P1GLSEnum:
     1412                        this->element_type = P1Enum;
     1413                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1414                        this->element_type = P1P1GLSEnum;
     1415                        return;
     1416                case MINIcondensedEnum:
     1417                        this->element_type = P1Enum;
     1418                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1419                        this->element_type = MINIcondensedEnum;
     1420                        return;
     1421                case MINIEnum:
     1422                        this->element_type = P1Enum;
     1423                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1424                        this->element_type = MINIEnum;
     1425                        return;
     1426                default:
     1427                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1428        }
    12311429}
    12321430/*}}}*/
    … …  
    12711469                case P1bubbleEnum:
    12721470                        /*Nodal function 1*/
    1273                         dbasis[NUMNODESP1*0+0]   = (zeta-1.)/4.;
    1274                         dbasis[NUMNODESP1*1+0]   = SQRT3/12.*(zeta-1.);
    1275                         dbasis[NUMNODESP1*2+0]   = -.5*gauss->coord1;
     1471                        dbasis[NUMNODESP1b*0+0]   = (zeta-1.)/4.;
     1472                        dbasis[NUMNODESP1b*1+0]   = SQRT3/12.*(zeta-1.);
     1473                        dbasis[NUMNODESP1b*2+0]   = -.5*gauss->coord1;
    12761474                        /*Nodal function 2*/
    1277                         dbasis[NUMNODESP1*0+1]   = (1.-zeta)/4.;
    1278                         dbasis[NUMNODESP1*1+1]   = SQRT3/12.*(zeta-1);
    1279                         dbasis[NUMNODESP1*2+1]   = -.5*gauss->coord2;
     1475                        dbasis[NUMNODESP1b*0+1]   = (1.-zeta)/4.;
     1476                        dbasis[NUMNODESP1b*1+1]   = SQRT3/12.*(zeta-1);
     1477                        dbasis[NUMNODESP1b*2+1]   = -.5*gauss->coord2;
    12801478                        /*Nodal function 3*/
    1281                         dbasis[NUMNODESP1*0+2]   = 0.;
    1282                         dbasis[NUMNODESP1*1+2]   = SQRT3/6.*(1.-zeta);
    1283                         dbasis[NUMNODESP1*2+2]   = -.5*gauss->coord3;
     1479                        dbasis[NUMNODESP1b*0+2]   = 0.;
     1480                        dbasis[NUMNODESP1b*1+2]   = SQRT3/6.*(1.-zeta);
     1481                        dbasis[NUMNODESP1b*2+2]   = -.5*gauss->coord3;
    12841482                        /*Nodal function 4*/
    1285                         dbasis[NUMNODESP1*0+3]   = -(1.+zeta)/4.;
    1286                         dbasis[NUMNODESP1*1+3]   = -SQRT3/12.*(1.+zeta);
    1287                         dbasis[NUMNODESP1*2+3]   = .5*gauss->coord1;
     1483                        dbasis[NUMNODESP1b*0+3]   = -(1.+zeta)/4.;
     1484                        dbasis[NUMNODESP1b*1+3]   = -SQRT3/12.*(1.+zeta);
     1485                        dbasis[NUMNODESP1b*2+3]   = .5*gauss->coord1;
    12881486                        /*Nodal function 5*/
    1289                         dbasis[NUMNODESP1*0+4]   = (1.+zeta)/4.;
    1290                         dbasis[NUMNODESP1*1+4]   = -SQRT3/12.*(1.+zeta);
    1291                         dbasis[NUMNODESP1*2+4]   = .5*gauss->coord2;
     1487                        dbasis[NUMNODESP1b*0+4]   = (1.+zeta)/4.;
     1488                        dbasis[NUMNODESP1b*1+4]   = -SQRT3/12.*(1.+zeta);
     1489                        dbasis[NUMNODESP1b*2+4]   = .5*gauss->coord2;
    12921490                        /*Nodal function 6*/
    1293                         dbasis[NUMNODESP1*0+5]   = 0.;
    1294                         dbasis[NUMNODESP1*1+5]   = SQRT3/6.*(1.+zeta);
    1295                         dbasis[NUMNODESP1*2+5]   = .5*gauss->coord3;
     1491                        dbasis[NUMNODESP1b*0+5]   = 0.;
     1492                        dbasis[NUMNODESP1b*1+5]   = SQRT3/6.*(1.+zeta);
     1493                        dbasis[NUMNODESP1b*2+5]   = .5*gauss->coord3;
    12961494                        /*Nodal function 7*/
    12971495                        dbasis[NUMNODESP1b*0+6] = 27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
    … …  
    17041902
    17051903        switch(this->element_type){
    1706                 case P1Enum:    return NUMNODESP1;
    1707                 case P2Enum:    return NUMNODESP2;
    1708                 case P2xP1Enum: return NUMNODESP2xP1;
    1709                 case P1xP2Enum: return NUMNODESP1xP2;
    1710                 case MINIEnum:  return NUMNODESP1b;
     1904                case P1Enum:            return NUMNODESP1;
     1905                case P1bubbleEnum:      return NUMNODESP1b;
     1906                case P2Enum:            return NUMNODESP2;
     1907                case P2xP1Enum:         return NUMNODESP2xP1;
     1908                case P1xP2Enum:         return NUMNODESP1xP2;
     1909                case P1P1Enum:          return NUMNODESP1*2;
     1910                case P1P1GLSEnum:       return NUMNODESP1*2;
     1911                case MINIcondensedEnum: return NUMNODESP1*2;
     1912                case MINIEnum:          return NUMNODESP1b+NUMNODESP1;
    17111913                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    17121914        }
    … …  
    17151917}
    17161918/*}}}*/
     1919/*FUNCTION PentaRef::NumberofNodesPressure{{{*/
     1920int PentaRef::NumberofNodesPressure(void){
     1921
     1922        switch(this->element_type){
     1923                case P1P1Enum:          return NUMNODESP1;
     1924                case P1P1GLSEnum:       return NUMNODESP1;
     1925                case MINIcondensedEnum: return NUMNODESP1;
     1926                case MINIEnum:          return NUMNODESP1;
     1927                default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1928        }
     1929
     1930        return -1;
     1931}
     1932/*}}}*/
     1933/*FUNCTION PentaRef::NumberofNodesVelocity{{{*/
     1934int PentaRef::NumberofNodesVelocity(void){
     1935
     1936        switch(this->element_type){
     1937                case P1P1Enum:          return NUMNODESP1;
     1938                case P1P1GLSEnum:       return NUMNODESP1;
     1939                case MINIcondensedEnum: return NUMNODESP1b;
     1940                case MINIEnum:          return NUMNODESP1b;
     1941                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1942        }
     1943
     1944        return -1;
     1945}
     1946/*}}}*/
     1947/*FUNCTION PentaRef::NumberofNodesVelocityFinal{{{*/
     1948int PentaRef::NumberofNodesVelocityFinal(void){
     1949
     1950        /*When static condensation is applied, the final number of nodes might be different*/
     1951
     1952        switch(this->element_type){
     1953                case P1P1Enum:          return NUMNODESP1;
     1954                case P1P1GLSEnum:       return NUMNODESP1;
     1955                case MINIcondensedEnum: return NUMNODESP1;
     1956                case MINIEnum:          return NUMNODESP1b;
     1957                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1958        }
     1959
     1960        return -1;
     1961}
     1962/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/PentaRef.h

    r15567 r15654  
    2323                /*Numerics*/
    2424                void GetNodalFunctions(IssmDouble* basis, GaussPenta* gauss);
     25                void GetNodalFunctionsVelocity(IssmDouble* basis, GaussPenta* gauss);
     26                void GetNodalFunctionsPressure(IssmDouble* basis, GaussPenta* gauss);
    2527                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,GaussPenta* gauss);
     28                void GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,GaussPenta* gauss);
     29                void GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list,GaussPenta* gauss);
    2630                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussPenta* gauss);
    2731                void GetNodalFunctionsP1(IssmDouble* l1l6, GaussPenta* gauss);
    … …  
    4145                void GetBSSAFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    4246                void GetBHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
     47                void GetBFSstrainrate(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    4348                void GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    4449                void GetBFSGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    … …  
    6570
    6671                int  NumberofNodes(void);
     72                int  NumberofNodesVelocity(void);
     73                int  NumberofNodesVelocityFinal(void);
     74                int  NumberofNodesPressure(void);
    6775};
    6876#endif

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

    r15645 r15654  
    33533353
    33543354        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    3355         /*P1 element only for now*/
    33563355        gauss=new GaussTria();
    33573356        for(int i=0;i<numnodes;i++){
    3358 
    33593357                gauss->GaussNode(numnodes,i);
    33603358

  • issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.cpp

    r15104 r15654  
    110110}
    111111/*}}}*/
     112/*FUNCTION DatasetInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index){{{*/
     113void DatasetInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index){
     114
     115        /*Get requested input within dataset*/
     116        if(index<0 || index > inputs->Size()-1) _error_("index requested (" << index << ") exceeds dataset size (" << inputs->Size() << ")");
     117        Input* input=dynamic_cast<Input*>(this->inputs->GetObjectByOffset(index));
     118
     119        input->GetInputValue(pvalue,gauss);
     120}
     121/*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h

    r15564 r15654  
    5050                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error_("not implemented yet");};
    5151                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index);
    52                 void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error_("not implemented yet");};
     52                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index);
    5353                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");};
    5454                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp

    r15643 r15654  
    134134/*FUNCTION PentaInput::GetVxStrainRate3d{{{*/
    135135void PentaInput::GetVxStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){
    136         int i,j;
    137 
    138         const int numnodes=6;
    139         const int DOFVELOCITY=3;
    140         IssmDouble B[8][27];
    141         IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
    142         IssmDouble velocity[numnodes][DOFVELOCITY];
    143         _assert_(this->NumberofNodes()==6); //Check Tria too
     136
     137        /*Intermediary*/
     138        int         numnodes=this->NumberofNodes();
     139        IssmDouble* B=xNew<IssmDouble>(numnodes*(NDOF3*numnodes));
     140        IssmDouble* velocity=xNew<IssmDouble>(numnodes*NDOF3);
    144141
    145142        /*Get B matrix: */
    146         GetBFS(&B[0][0], xyz_list, gauss);
    147 
    148         /*Create a reduced matrix of B to get rid of pressure */
    149         for (i=0;i<6;i++){
    150                 for (j=0;j<3;j++){
    151                         B_reduced[i][j]=B[i][j];
    152                 }
    153                 for (j=4;j<7;j++){
    154                         B_reduced[i][j-1]=B[i][j];
    155                 }
    156                 for (j=8;j<11;j++){
    157                         B_reduced[i][j-2]=B[i][j];
    158                 }
    159                 for (j=12;j<15;j++){
    160                         B_reduced[i][j-3]=B[i][j];
    161                 }
    162                 for (j=16;j<19;j++){
    163                         B_reduced[i][j-4]=B[i][j];
    164                 }
    165                 for (j=20;j<23;j++){
    166                         B_reduced[i][j-5]=B[i][j];
    167                 }
    168         }
     143        GetBFSstrainrate(B,xyz_list,gauss);
    169144
    170145        /*Here, we are computing the strain rate of (vx,0,0)*/
    171         for(i=0;i<numnodes;i++){
    172                 velocity[i][0]=this->values[i];
    173                 velocity[i][1]=0.0;
    174                 velocity[i][2]=0.0;
     146        for(int i=0;i<numnodes;i++){
     147                velocity[NDOF3*i+0]=this->values[i];
     148                velocity[NDOF3*i+1]=0.;
     149                velocity[NDOF3*i+2]=0.;
    175150        }
    176151        /*Multiply B by velocity, to get strain rate: */
    177         MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvx,0);
     152        MatrixMultiply(B,numnodes,NDOF3*numnodes,0,velocity,NDOF3*numnodes,1,0,epsilonvx,0);
     153
     154        /*Clean-up*/
     155        xDelete<IssmDouble>(B);
     156        xDelete<IssmDouble>(velocity);
    178157
    179158}
    180159/*}}}*/
    181160/*FUNCTION PentaInput::GetVyStrainRate3d{{{*/
    182 void PentaInput::GetVyStrainRate3d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){
    183         int i,j;
    184 
    185         const int numnodes=6;
    186         const int DOFVELOCITY=3;
    187         IssmDouble B[8][27];
    188         IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
    189         IssmDouble velocity[numnodes][DOFVELOCITY];
    190 
    191         _assert_(this->NumberofNodes()==6); //Check Tria too
     161void PentaInput::GetVyStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){
     162
     163        /*Intermediary*/
     164        int         numnodes=this->NumberofNodes();
     165        IssmDouble* B=xNew<IssmDouble>(numnodes*(NDOF3*numnodes));
     166        IssmDouble* velocity=xNew<IssmDouble>(numnodes*NDOF3);
    192167
    193168        /*Get B matrix: */
    194         GetBFS(&B[0][0], xyz_list, gauss);
    195         /*Create a reduced matrix of B to get rid of pressure */
    196         for (i=0;i<6;i++){
    197                 for (j=0;j<3;j++){
    198                         B_reduced[i][j]=B[i][j];
    199                 }
    200                 for (j=4;j<7;j++){
    201                         B_reduced[i][j-1]=B[i][j];
    202                 }
    203                 for (j=8;j<11;j++){
    204                         B_reduced[i][j-2]=B[i][j];
    205                 }
    206                 for (j=12;j<15;j++){
    207                         B_reduced[i][j-3]=B[i][j];
    208                 }
    209                 for (j=16;j<19;j++){
    210                         B_reduced[i][j-4]=B[i][j];
    211                 }
    212                 for (j=20;j<23;j++){
    213                         B_reduced[i][j-5]=B[i][j];
    214                 }
    215         }
     169        GetBFSstrainrate(B,xyz_list,gauss);
    216170
    217171        /*Here, we are computing the strain rate of (0,vy,0)*/
    218         for(i=0;i<numnodes;i++){
    219                 velocity[i][0]=0.0;
    220                 velocity[i][1]=this->values[i];
    221                 velocity[i][2]=0.0;
     172        for(int i=0;i<numnodes;i++){
     173                velocity[NDOF3*i+0]=0.;
     174                velocity[NDOF3*i+1]=this->values[i];
     175                velocity[NDOF3*i+2]=0.;
    222176        }
    223177        /*Multiply B by velocity, to get strain rate: */
    224         MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvy,0);
    225 
     178        MatrixMultiply(B,numnodes,NDOF3*numnodes,0,velocity,NDOF3*numnodes,1,0,epsilonvx,0);
     179
     180        /*Clean-up*/
     181        xDelete<IssmDouble>(B);
     182        xDelete<IssmDouble>(velocity);
    226183}
    227184/*}}}*/
    228185/*FUNCTION PentaInput::GetVzStrainRate3d{{{*/
    229 void PentaInput::GetVzStrainRate3d(IssmDouble* epsilonvz,IssmDouble* xyz_list, GaussPenta* gauss){
    230         int i,j;
    231 
    232         const int numnodes=6;
    233         const int DOFVELOCITY=3;
    234         IssmDouble B[8][27];
    235         IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
    236         IssmDouble velocity[numnodes][DOFVELOCITY];
     186void PentaInput::GetVzStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){
     187
     188        /*Intermediary*/
     189        int         numnodes=this->NumberofNodes();
     190        IssmDouble* B=xNew<IssmDouble>(numnodes*(NDOF3*numnodes));
     191        IssmDouble* velocity=xNew<IssmDouble>(numnodes*NDOF3);
    237192
    238193        /*Get B matrix: */
    239         GetBFS(&B[0][0], xyz_list, gauss);
    240 
    241         _assert_(this->NumberofNodes()==6); //Check Tria too
    242 
    243         /*Create a reduced matrix of B to get rid of pressure */
    244         for (i=0;i<6;i++){
    245                 for (j=0;j<3;j++){
    246                         B_reduced[i][j]=B[i][j];
    247                 }
    248                 for (j=4;j<7;j++){
    249                         B_reduced[i][j-1]=B[i][j];
    250                 }
    251                 for (j=8;j<11;j++){
    252                         B_reduced[i][j-2]=B[i][j];
    253                 }
    254                 for (j=12;j<15;j++){
    255                         B_reduced[i][j-3]=B[i][j];
    256                 }
    257                 for (j=16;j<19;j++){
    258                         B_reduced[i][j-4]=B[i][j];
    259                 }
    260                 for (j=20;j<23;j++){
    261                         B_reduced[i][j-5]=B[i][j];
    262                 }
    263         }
     194        GetBFSstrainrate(B,xyz_list,gauss);
    264195
    265196        /*Here, we are computing the strain rate of (0,0,vz)*/
    266         for(i=0;i<numnodes;i++){
    267                 velocity[i][0]=0.0;
    268                 velocity[i][1]=0.0;
    269                 velocity[i][2]=this->values[i];
    270         }
    271 
     197        for(int i=0;i<numnodes;i++){
     198                velocity[NDOF3*i+0]=0.;
     199                velocity[NDOF3*i+1]=0.;
     200                velocity[NDOF3*i+2]=this->values[i];
     201        }
    272202        /*Multiply B by velocity, to get strain rate: */
    273         MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvz,0);
    274 
     203        MatrixMultiply(B,numnodes,NDOF3*numnodes,0,velocity,NDOF3*numnodes,1,0,epsilonvx,0);
     204
     205        /*Clean-up*/
     206        xDelete<IssmDouble>(B);
     207        xDelete<IssmDouble>(velocity);
    275208}
    276209/*}}}*/

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

    r15564 r15654  
    475475
    476476        /*Transform Coordinate System*/
    477         TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,XYZPEnum);
     477        TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,XYZEnum);
    478478
    479479        /*Clean up and return*/

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

    r15567 r15654  
    561561                exz=epsilon[3];
    562562                eyz=epsilon[4];
     563                eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;
     564
     565                mu_prime=(1-n)/(2*n) * mu/eff2;
     566        }
     567
     568        /*Assign output pointers:*/
     569        *pmu_prime=mu_prime;
     570}
     571/*}}}*/
     572/*FUNCTION Matdamageice::GetViscosityDerivativeEpsSquareFS{{{*/
     573void  Matdamageice::GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* epsilon){
     574
     575        /*output: */
     576        IssmDouble mu_prime;
     577        IssmDouble mu,n,eff2;
     578
     579        /*input strain rate: */
     580        IssmDouble exx,eyy,ezz,exy,exz,eyz;
     581
     582        /*Get visocisty and n*/
     583        GetViscosity3d(&mu,epsilon);
     584        n=GetN();
     585
     586        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
     587                                (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){
     588                mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14);
     589        }
     590        else{
     591                /*Retrive strain rate components: */
     592                exx=epsilon[0];
     593                eyy=epsilon[1];
     594                ezz=epsilon[2];
     595                exy=epsilon[3];
     596                exz=epsilon[4];
     597                eyz=epsilon[5];
    563598                eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;
    564599

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

    r15564 r15654  
    5656                void   GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
    5757                void   GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
     58                void   GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    5859                void   GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    5960                IssmDouble GetA();

  • issm/trunk-jpl/src/c/classes/Materials/Material.h

    r15564 r15654  
    3232                virtual void       GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;
    3333                virtual void       GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0;
     34                virtual void       GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0;
    3435                virtual void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0;
    3536                virtual IssmDouble GetA()=0;

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

    r15567 r15654  
    495495                exz=epsilon[3];
    496496                eyz=epsilon[4];
     497                eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;
     498
     499                mu_prime=(1-n)/(2*n) * mu/eff2;
     500        }
     501
     502        /*Assign output pointers:*/
     503        *pmu_prime=mu_prime;
     504}
     505/*}}}*/
     506/*FUNCTION Matice::GetViscosityDerivativeEpsSquareFS{{{*/
     507void  Matice::GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* epsilon){
     508
     509        /*output: */
     510        IssmDouble mu_prime;
     511        IssmDouble mu,n,eff2;
     512
     513        /*input strain rate: */
     514        IssmDouble exx,eyy,exy,exz,eyz,ezz;
     515
     516        /*Get visocisty and n*/
     517        GetViscosity3d(&mu,epsilon);
     518        n=GetN();
     519
     520        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
     521                                (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){
     522                mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)14);
     523        }
     524        else{
     525                /*Retrive strain rate components: */
     526                exx=epsilon[0];
     527                eyy=epsilon[1];
     528                ezz=epsilon[2];
     529                exy=epsilon[3];
     530                exz=epsilon[4];
     531                eyz=epsilon[5];
    497532                eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;
    498533

  • issm/trunk-jpl/src/c/classes/Materials/Matice.h

    r15564 r15654  
    6363                void GetViscosityZComplement(IssmDouble*, IssmDouble*){_error_("not supported");};
    6464                void       GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
     65                void       GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    6566                void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    6667                IssmDouble GetA();

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r15564 r15654  
    9393                void       GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};
    9494                void       GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");};
     95                void       GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");};
    9596                void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");};
    9697                IssmDouble GetA(){_error_("not supported");};

  • issm/trunk-jpl/src/c/classes/Node.cpp

    r15634 r15654  
    6767                /*Coordinate system provided, convert to coord_system matrix*/
    6868                _assert_(iomodel->Data(DiagnosticReferentialEnum));
    69                 XZvectorsToCoordinateSystem(&this->coord_system[0][0],iomodel->Data(DiagnosticReferentialEnum)+io_index*6);
     69                XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(DiagnosticReferentialEnum)[io_index*6]);
    7070
    7171                if(iomodel->dim==3){
    … …  
    438438        /*Put dof for this node into the s set (ie, this dof will be constrained
    439439         * to a fixed value during computations. */
    440 
    441         if(dof>=this->indexing.gsize){
    442                 printf("dof spc = %i\n",dof);
    443                 this->Echo();
    444         }
    445440        _assert_(dof<this->indexing.gsize);
    446441
    … …  
    10101005        for(i=0;i<numnodes;i++){
    10111006                switch(cs_array[i]){
    1012                         case XYEnum:   numdofs+=2; break;
    1013                         case XYZPEnum: numdofs+=4; break;
     1007                        case PressureEnum: numdofs+=1; break;
     1008                        case XYEnum:       numdofs+=2; break;
     1009                        case XYZEnum:      numdofs+=3; break;
    10141010                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    10151011                }
    … …  
    10591055        for(i=0;i<numnodes;i++){
    10601056                switch(cs_array[i]){
    1061                         case XYEnum:   numdofs+=2; break;
    1062                         case XYZPEnum: numdofs+=4; break;
     1057                        case PressureEnum: numdofs+=1; break;
     1058                        case XYEnum:       numdofs+=2; break;
     1059                        case XYZEnum:      numdofs+=3; break;
    10631060                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    10641061                }
    … …  
    11071104        for(i=0;i<numnodes;i++){
    11081105                switch(cs_array[i]){
    1109                         case XYEnum:   numdofs+=2; break;
    1110                         case XYZPEnum: numdofs+=4; break;
     1106                        case PressureEnum: numdofs+=1; break;
     1107                        case XYEnum:       numdofs+=2; break;
     1108                        case XYZEnum:      numdofs+=3; break;
    11111109                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    11121110                }
    … …  
    11541152        for(int i=0;i<numnodes;i++){
    11551153                switch(cs_array[i]){
    1156                         case XYEnum:   numdofs+=2; break;
    1157                         case XYZPEnum: numdofs+=4; break;
     1154                        case PressureEnum: numdofs+=1; break;
     1155                        case XYEnum:       numdofs+=2; break;
     1156                        case XYZEnum:      numdofs+=3; break;
    11581157                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    11591158                }
    … …  
    11931192        for(i=0;i<numnodes;i++){
    11941193                switch(cs_array[i]){
    1195                         case XYEnum:   numdofs+=2; break;
    1196                         case XYZPEnum: numdofs+=4; break;
     1194                        case PressureEnum: numdofs+=1; break;
     1195                        case XYEnum:       numdofs+=2; break;
     1196                        case XYZEnum:      numdofs+=3; break;
    11971197                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    11981198                }
    … …  
    12161216                nodes[i]->GetCoordinateSystem(&coord_system[0][0]);
    12171217                switch(cs_array[i]){
     1218                        case PressureEnum:
     1219                                /*DO NOT change anything*/
     1220                                transform[(numdofs)*(counter) + counter] = 1.;
     1221                                counter+=1;
     1222                                break;
    12181223                        case XYEnum:
    12191224                                /*We remove the z component, we need to renormalize x and y: x=[x1 x2 0] y=[-x2 x1 0]*/
    … …  
    12251230                                counter+=2;
    12261231                                break;
    1227                         case XYZPEnum:
    1228                                 /*Only the first 3 coordinates are changed (x,y,z), leave the others (P) unchanged*/
     1232                        case XYZEnum:
     1233                                /*The 3 coordinates are changed (x,y,z)*/
    12291234                                transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0];
    12301235                                transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1];
    … …  
    12361241                                transform[(numdofs)*(counter+2) + counter+1] = coord_system[2][1];
    12371242                                transform[(numdofs)*(counter+2) + counter+2] = coord_system[2][2];
    1238                                 transform[(numdofs)*(counter+3) + counter+3] = 1.0;
    1239                                 counter+=4;
     1243                                counter+=3;
    12401244                                break;
    12411245                        default:

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

    r15636 r15654  
    448448                        }
    449449                        break;
     450                case 7: //P1+ Lagrange element
     451                        switch(iv){
     452                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
     453                                case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
     454                                case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
     455                                case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
     456                                case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
     457                                case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
     458
     459                                case 6: coord1=0.; coord2=0.; coord3=0.; coord4=0.;break;
     460                                default: _error_("node index should be in [0 5]");
     461                        }
     462                        break;
    450463                case 15: //P2 Lagrange element
    451464                        switch(iv){
    … …  
    470483                        }
    471484                        break;
    472                 default: _error_("supported number of nodes are 6 and 15");
     485                default: _error_("Number of nodes "<<numnodes<<" not supported");
    473486        }
    474487

  • issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp

    r15104 r15654  
    223223        this->col_slocaldoflist=NULL;
    224224        this->col_sglobaldoflist=NULL;
    225 
    226225}
    227226/*}}}*/
    … …  
    423422
    424423        _assert_(Ke);
     424        _assert_(this);
    425425
    426426        this->nrows =Ke->nrows;

  • issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp

    r15636 r15654  
    1111
    1212        /*intermediary: */
    13         int         i,j,count;
    14         IssmDouble  yts;
    15         FILE       *fid              = NULL;
     13        FILE       *fid = NULL;
    1614        int         code,vector_layout;
    17         IssmDouble *times            = NULL;
    18         IssmDouble *values           = NULL;
    19         bool        spcpresent       = false;
    20 
    21         /*P2 finite elements*/
    22         int   v1,v2;
    23         bool *my_edges = NULL;
    24 
    25         /*variables being fetched: */
    2615        IssmDouble *spcdata = NULL;
    2716        int         M,N;
    28 
    29         /*Fetch parameters: */
    30         iomodel->Constant(&yts,ConstantsYtsEnum);
    3117
    3218        /*First of, find the record for the enum, and get code  of data type: */
    … …  
    3824        iomodel->FetchData(&spcdata,&M,&N,vector_enum);
    3925
     26        /*Call IoModelToConstraintsx*/
     27        IoModelToConstraintsx(constraints,iomodel,spcdata,M,N,analysis_type,finite_element,dof);
     28
     29        /*Clean up*/
     30        xDelete<IssmDouble>(spcdata);
     31}
     32
     33void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,IssmDouble* spcdata,int M,int N,int analysis_type,int finite_element,int dof){
     34
     35        /*intermediary: */
     36        int         i,j,count,elementnbv;
     37        IssmDouble  value;
     38        IssmDouble *times            = NULL;
     39        IssmDouble *values           = NULL;
     40        bool        spcpresent       = false;
     41
     42        /*P2 finite elements*/
     43        int   v1,v2;
     44        bool *my_edges = NULL;
     45
    4046        switch(finite_element){
    4147                case P1Enum:
    4248                        /*Nothing else to do*/
     49                        break;
     50                case P1bubbleEnum:
     51                        if(iomodel->dim!=3) _error_("3d is the only supported dimension");
     52                        elementnbv = 6;
    4353                        break;
    4454                case P1xP2Enum:
    … …  
    8494                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
    8595                                                                                        dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
     96                                                        count++;
     97                                                }
     98                                        }
     99                                }
     100                                break;
     101                        case P1bubbleEnum:
     102                                for(i=0;i<iomodel->numberofvertices;i++){
     103                                        if((iomodel->my_vertices[i])){
     104                                                if (!xIsNan<IssmDouble>(spcdata[i])){
     105                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,dof,spcdata[i],analysis_type));
     106                                                        count++;
     107                                                }
     108                                        }
     109                                }
     110                                for(i=0;i<iomodel->numberofelements;i++){
     111                                        if(iomodel->my_elements[i]){
     112                                                value = spcdata[iomodel->elements[i*elementnbv+0]-1];
     113                                                for(j=1;j<elementnbv;j++) value += spcdata[iomodel->elements[i*elementnbv+j]-1];
     114                                                value = value/reCast<IssmDouble,int>(elementnbv+0);
     115                                                if(!xIsNan<IssmDouble>(value)){
     116                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
     117                                                                                        dof,value,analysis_type));
    86118                                                        count++;
    87119                                                }
    … …  
    145177                /*figure out times: */
    146178                times=xNew<IssmDouble>(N);
    147                 for(j=0;j<N;j++) times[j]=spcdata[(M-1)*N+j]*yts;
     179                for(j=0;j<N;j++) times[j]=spcdata[(M-1)*N+j];
    148180
    149181                switch(finite_element){
    … …  
    291323        }
    292324        else{
    293                 _error_("Size of field " << EnumToStringx(vector_enum) << " not supported");
     325                _error_("Size of spc field not supported");
    294326        }
    295327
    296328        /*Free ressources:*/
    297         xDelete<IssmDouble>(spcdata);
    298329        xDelete<IssmDouble>(times);
    299330        xDelete<IssmDouble>(values);

  • issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h

    r15621 r15654  
    99/* local prototypes: */
    1010void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element,int dof=1);
     11void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,IssmDouble* spcdata,int M,int N,int analysis_type,int finite_element,int dof=1);
    1112
    1213#endif  /* _IOMODELTOELEMENTINPUTX_H */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp

    r15636 r15654  
    9797                        break;
    9898
    99                 case MINIcondensedEnum:
     99                /*Stokes elements*/
     100                case P1P1Enum:
    100101                        _assert_(approximation==FSApproximationEnum);
    101102                        /*P1 velocity*/
    102103                        for(i=0;i<iomodel->numberofvertices;i++){
    103104                                if(iomodel->my_vertices[i]){
    104                                         nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,approximation));
     105                                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum));
     106                                }
     107                        }
     108                        /*P1 pressure*/
     109                        for(i=0;i<iomodel->numberofvertices;i++){
     110                                if(iomodel->my_vertices[i]){
     111                                        nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,i,iomodel,analysis,FSpressureEnum));
     112                                }
     113                        }
     114                        break;
     115                case P1P1GLSEnum:
     116                        _assert_(approximation==FSApproximationEnum);
     117                        /*P1 velocity*/
     118                        for(i=0;i<iomodel->numberofvertices;i++){
     119                                if(iomodel->my_vertices[i]){
     120                                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum));
     121                                }
     122                        }
     123                        /*P1 pressure*/
     124                        for(i=0;i<iomodel->numberofvertices;i++){
     125                                if(iomodel->my_vertices[i]){
     126                                        nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,i,iomodel,analysis,FSpressureEnum));
     127                                }
     128                        }
     129                        break;
     130                case MINIcondensedEnum:
     131                        _assert_(approximation==FSApproximationEnum);
     132                        /*P1 velocity (bubble statically condensed)*/
     133                        for(i=0;i<iomodel->numberofvertices;i++){
     134                                if(iomodel->my_vertices[i]){
     135                                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum));
     136                                }
     137                        }
     138                        /*P1 pressure*/
     139                        for(i=0;i<iomodel->numberofvertices;i++){
     140                                if(iomodel->my_vertices[i]){
     141                                        nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,i,iomodel,analysis,FSpressureEnum));
     142                                }
     143                        }
     144                        break;
     145                case MINIEnum:
     146                        _assert_(approximation==FSApproximationEnum);
     147                        /*P1+ velocity*/
     148                        for(i=0;i<iomodel->numberofvertices;i++){
     149                                if(iomodel->my_vertices[i]){
     150                                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum));
     151                                }
     152                        }
     153                        for(i=0;i<iomodel->numberofelements;i++){
     154                                if(iomodel->my_elements[i]){
     155                                        nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,analysis,FSvelocityEnum));
    105156                                }
    106157                        }
    … …  
    108159                        for(i=0;i<iomodel->numberofvertices;i++){
    109160                                if(iomodel->my_vertices[i]){
    110                                         nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,i,i,iomodel,analysis,approximation));
     161                                        nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,i,iomodel,analysis,FSpressureEnum));
    111162                                }
    112163                        }

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp

    r15644 r15654  
    6868
    6969        /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
    70         if(!iscoupling && !isFS){
     70        if(!iscoupling){
    7171
    7272                /*Get finite element type*/
    … …  
    9494                else if(isFS){
    9595                        finiteelement = P1Enum;
     96                        iomodel->Constant(&temp,FlowequationFeFSEnum);
     97                        switch(temp){
     98                                case 0 : finiteelement = P1Enum; break;//P1P1
     99                                case 1 : finiteelement = P1Enum; break;//P1P1GSL
     100                                case 2 : finiteelement = P1Enum; break;//MINIcondensedEnum
     101                                case 3 : finiteelement = P1bubbleEnum; break;//MINIEnum
     102                                default: _error_("finite element "<<temp<<" not supported");
     103                        }
    96104                }
    97105                IoModelToConstraintsx(constraints,iomodel,DiagnosticSpcvxEnum,DiagnosticHorizAnalysisEnum,finiteelement,1);
    98106                IoModelToConstraintsx(constraints,iomodel,DiagnosticSpcvyEnum,DiagnosticHorizAnalysisEnum,finiteelement,2);
     107
     108                if(isFS){
     109
     110                        /*Constraint at the bedrock interface (v.n = vz = 0) (Coordinates will be updated according to the bed slope)*/
     111                        iomodel->FetchData(&spcvz,&Mz,&Nz,DiagnosticSpcvzEnum);
     112                        iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
     113                        iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum);
     114                        iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
     115                        iomodel->FetchData(&nodeonicesheet,NULL,NULL,MaskVertexongroundediceEnum);
     116                        for(i=0;i<iomodel->numberofvertices;i++){
     117                                if(iomodel->my_vertices[i]){
     118                                        if(reCast<int,IssmDouble>(nodeonbed[i]) && reCast<int,IssmDouble>(nodeonicesheet[i]) && reCast<int,IssmDouble>(nodeonFS[i])){
     119                                                if(vertices_type[i] == FSApproximationEnum){
     120                                                        for(j=0;j<Nz;j++) spcvz[i*Nz+j] = 0.;
     121                                                }
     122                                                else{
     123                                                        _error_("not supported");
     124                                                }
     125                                        }
     126                                }
     127                        }
     128                        IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,DiagnosticHorizAnalysisEnum,finiteelement,3);
     129                        iomodel->DeleteData(spcvz,DiagnosticSpcvzEnum);
     130                        iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
     131                        iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum);
     132                        iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum);
     133                        iomodel->DeleteData(nodeonicesheet,MaskVertexongroundediceEnum);
     134
     135                        /*Pressure spc*/
     136                        count = constraints->Size();
     137                        iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
     138                        iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum);
     139                        iomodel->FetchData(&z,NULL,NULL,MeshZEnum);
     140                        for(i=0;i<iomodel->numberofvertices;i++){
     141                                if(iomodel->my_vertices[i]){
     142                                        if (reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
     143                                                constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,DiagnosticHorizAnalysisEnum));
     144                                                count++;
     145                                        }
     146                                }
     147                        }
     148                        iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
     149                        iomodel->DeleteData(surface,SurfaceEnum);
     150                        iomodel->DeleteData(z,MeshZEnum);
     151                }
    99152
    100153                *pconstraints=constraints;
    … …  
    325378                                }
    326379                                if (reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
    327                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,g*rho_ice*(surface[i]-z[i])/FSreconditioning,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
     380                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
    328381                                        count++;
    329382                                }

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp

    r15644 r15654  
    6060                else if(isFS){
    6161                        approximation = FSApproximationEnum;
    62                         finiteelement = P1Enum;
     62                        iomodel->Constant(&temp,FlowequationFeFSEnum);
     63                        switch(temp){
     64                                case 0 : finiteelement = P1P1Enum;          break;
     65                                case 1 : finiteelement = P1P1GLSEnum;       break;
     66                                case 2 : finiteelement = MINIcondensedEnum; break;
     67                                case 3 : finiteelement = MINIEnum;          break;
     68                                default: _error_("finite element "<<temp<<" not supported");
     69                        }
    6370                }
    6471

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp

    r15644 r15654  
    6464                }
    6565                else if(isFS){
    66                         finiteelement = P1Enum;
     66                        iomodel->Constant(&temp,FlowequationFeFSEnum);
     67                        switch(temp){
     68                                case 0 : finiteelement = P1P1Enum;          break;
     69                                case 1 : finiteelement = P1P1GLSEnum;       break;
     70                                case 2 : finiteelement = MINIcondensedEnum; break;
     71                                case 3 : finiteelement = MINIEnum;          break;
     72                                default: _error_("finite element "<<temp<<" not supported");
     73                        }
    6774                }
    6875        }

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp

    r15567 r15654  
    3232                                case FSApproximationEnum:
    3333                                        numdofs=4;
     34                                        break;
     35                                case FSvelocityEnum:
     36                                        numdofs=3;
     37                                        break;
     38                                case FSpressureEnum:
     39                                        numdofs=1;
    3440                                        break;
    3541                                case NoneApproximationEnum:

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r15636 r15654  
    319319        HOFSApproximationEnum,
    320320        FSApproximationEnum,
     321        FSvelocityEnum,
     322        FSpressureEnum,
    321323        /*}}}*/
    322324        /*Datasets {{{*/
    … …  
    504506        P2xP1Enum,
    505507        P1xP2Enum,
     508        P1P1Enum,
     509        P1P1GLSEnum,
    506510        MINIEnum,
    507511        MINIcondensedEnum,
    … …  
    594598        /*Coordinate Systems{{{*/
    595599        XYEnum,
    596         XYZPEnum,
     600        XYZEnum,
    597601        /*}}}*/
    598602        /*Toolkits{{{*/

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r15636 r15654  
    323323                case HOFSApproximationEnum : return "HOFSApproximation";
    324324                case FSApproximationEnum : return "FSApproximation";
     325                case FSvelocityEnum : return "FSvelocity";
     326                case FSpressureEnum : return "FSpressure";
    325327                case ConstraintsEnum : return "Constraints";
    326328                case LoadsEnum : return "Loads";
    … …  
    496498                case P2xP1Enum : return "P2xP1";
    497499                case P1xP2Enum : return "P1xP2";
     500                case P1P1Enum : return "P1P1";
     501                case P1P1GLSEnum : return "P1P1GLS";
    498502                case MINIEnum : return "MINI";
    499503                case MINIcondensedEnum : return "MINIcondensed";
    … …  
    570574                case NearestInterpEnum : return "NearestInterp";
    571575                case XYEnum : return "XY";
    572                 case XYZPEnum : return "XYZP";
     576                case XYZEnum : return "XYZ";
    573577                case DenseEnum : return "Dense";
    574578                case MpiDenseEnum : return "MpiDense";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r15636 r15654  
    329329              else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum;
    330330              else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
     331              else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
     332              else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
    331333              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    332334              else if (strcmp(name,"Loads")==0) return LoadsEnum;
    … …  
    381383              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    382384              else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    383               else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
    384               else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum;
     388              if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
     389              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
     390              else if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum;
    389391              else if (strcmp(name,"Segment")==0) return SegmentEnum;
    390392              else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
    … …  
    504506              else if (strcmp(name,"P2")==0) return P2Enum;
    505507              else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
    506               else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
    507               else if (strcmp(name,"MINI")==0) return MINIEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
     511              if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
     512              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
     513              else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
     514              else if (strcmp(name,"MINI")==0) return MINIEnum;
     515              else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
    512516              else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
    513517              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
    … …  
    582586              else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
    583587              else if (strcmp(name,"XY")==0) return XYEnum;
    584               else if (strcmp(name,"XYZP")==0) return XYZPEnum;
     588              else if (strcmp(name,"XYZ")==0) return XYZEnum;
    585589              else if (strcmp(name,"Dense")==0) return DenseEnum;
    586590              else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;

  • issm/trunk-jpl/src/m/classes/flowequation.m

    r15644 r15654  
    6969                function obj = setdefaultparameters(obj) % {{{
    7070
     71                        %MINI element for FS by default
     72                        obj.fe_FS = 3;
    7173                end % }}}
    7274                function md = checkconsistency(obj,md,solution,analyses) % {{{
    … …  
    8183                                md = checkfield(md,'flowequation.fe_SSA','numel',[1],'values',[0 1]);
    8284                                md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',[0:3]);
    83                                 md = checkfield(md,'flowequation.fe_FS','numel',[1],'values',[0]);
     85                                md = checkfield(md,'flowequation.fe_FS','numel',[1],'values',[0:3]);
    8486                                md = checkfield(md,'flowequation.borderSSA','size',[md.mesh.numberofvertices 1],'values',[0 1]);
    8587                                md = checkfield(md,'flowequation.borderHO','size',[md.mesh.numberofvertices 1],'values',[0 1]);

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r15636 r15654  
    42934293        return StringToEnum('FSApproximation')[0]
    42944294
     4295def FSvelocityEnum():
     4296        """
     4297        FSVELOCITYENUM - Enum of FSvelocity
     4298
     4299        WARNING: DO NOT MODIFY THIS FILE
     4300                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     4301                                Please read src/c/shared/Enum/README for more information
     4302
     4303           Usage:
     4304              macro=FSvelocityEnum()
     4305        """
     4306
     4307        return StringToEnum('FSvelocity')[0]
     4308
     4309def FSpressureEnum():
     4310        """
     4311        FSPRESSUREENUM - Enum of FSpressure
     4312
     4313        WARNING: DO NOT MODIFY THIS FILE
     4314                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     4315                                Please read src/c/shared/Enum/README for more information
     4316
     4317           Usage:
     4318              macro=FSpressureEnum()
     4319        """
     4320
     4321        return StringToEnum('FSpressure')[0]
     4322
    42954323def ConstraintsEnum():
    42964324        """
    … …  
    67156743        return StringToEnum('P1xP2')[0]
    67166744
     6745def P1P1Enum():
     6746        """
     6747        P1P1ENUM - Enum of P1P1
     6748
     6749        WARNING: DO NOT MODIFY THIS FILE
     6750                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     6751                                Please read src/c/shared/Enum/README for more information
     6752
     6753           Usage:
     6754              macro=P1P1Enum()
     6755        """
     6756
     6757        return StringToEnum('P1P1')[0]
     6758
     6759def P1P1GLSEnum():
     6760        """
     6761        P1P1GLSENUM - Enum of P1P1GLS
     6762
     6763        WARNING: DO NOT MODIFY THIS FILE
     6764                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     6765                                Please read src/c/shared/Enum/README for more information
     6766
     6767           Usage:
     6768              macro=P1P1GLSEnum()
     6769        """
     6770
     6771        return StringToEnum('P1P1GLS')[0]
     6772
    67176773def MINIEnum():
    67186774        """
    … …  
    77517807        return StringToEnum('XY')[0]
    77527808
    7753 def XYZPEnum():
    7754         """
    7755         XYZPENUM - Enum of XYZP
    7756 
    7757         WARNING: DO NOT MODIFY THIS FILE
    7758                                 this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
    7759                                 Please read src/c/shared/Enum/README for more information
    7760 
    7761            Usage:
    7762               macro=XYZPEnum()
    7763         """
    7764 
    7765         return StringToEnum('XYZP')[0]
     7809def XYZEnum():
     7810        """
     7811        XYZENUM - Enum of XYZ
     7812
     7813        WARNING: DO NOT MODIFY THIS FILE
     7814                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     7815                                Please read src/c/shared/Enum/README for more information
     7816
     7817           Usage:
     7818              macro=XYZEnum()
     7819        """
     7820
     7821        return StringToEnum('XYZ')[0]
    77667822
    77677823def DenseEnum():
    … …  
    79598015        """
    79608016
    7961         return 567
    7962 
     8017        return 571
     8018

  • issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m

    r15636 r15654  
    99%      macro=MaximumNumberOfEnums()
    1010
    11 macro=567;
     11macro=571;
Note: See TracChangeset for help on using the changeset viewer.