Changeset 18928


Ignore:
Timestamp:
12/04/14 09:39:27 (10 years ago)
Author:
seroussi
Message:

CHG: ordering StressbalanceAnalysis

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r18843 r18928  
    1111
    1212/*Model processing*/
     13void StressbalanceAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
     14
     15        /*Intermediary*/
     16        int        i,j;
     17        int        count,finiteelement;
     18        IssmDouble g;
     19        IssmDouble rho_ice;
     20        IssmDouble FSreconditioning;
     21        bool       isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
     22        bool       spcpresent = false;
     23        int        Mx,Nx;
     24        int        My,Ny;
     25        int        Mz,Nz;
     26        IssmDouble *spcvx          = NULL;
     27        IssmDouble *spcvy          = NULL;
     28        IssmDouble *spcvz          = NULL;
     29        IssmDouble *nodeonSSA = NULL;
     30        IssmDouble *nodeonHO   = NULL;
     31        IssmDouble *nodeonFS   = NULL;
     32        IssmDouble *nodeonbase      = NULL;
     33        IssmDouble *groundedice_ls = NULL;
     34        IssmDouble *vertices_type  = NULL;
     35        IssmDouble *surface        = NULL;
     36        IssmDouble *z              = NULL;
     37        IssmDouble *timesx=NULL;
     38        IssmDouble *timesy=NULL;
     39        IssmDouble *timesz=NULL;
     40   IssmDouble* values=NULL;
     41
     42        /*Fetch parameters: */
     43        iomodel->Constant(&g,ConstantsGEnum);
     44        iomodel->Constant(&rho_ice,MaterialsRhoIceEnum);
     45        iomodel->Constant(&FSreconditioning,StressbalanceFSreconditioningEnum);
     46        iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
     47        iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
     48        iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
     49        iomodel->Constant(&isHO,FlowequationIsHOEnum);
     50        iomodel->Constant(&isFS,FlowequationIsFSEnum);
     51
     52        /*Now, is the flag macayaealHO on? otherwise, do nothing: */
     53        if(!isSSA && !isHO && !isFS && !isL1L2) return;
     54
     55        /*Do we have coupling*/
     56        if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     57         iscoupling = true;
     58        else
     59         iscoupling = false;
     60
     61        /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
     62        if(!iscoupling){
     63
     64                /*Get finite element type*/
     65                if(isSSA)       iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
     66                else if(isL1L2) finiteelement = P1Enum;
     67                else if(isHO)   iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
     68                else if(isFS){  iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
     69                        /*Deduce velocity interpolation from finite element*/
     70                        switch(finiteelement){
     71                                case P1P1Enum              : finiteelement = P1Enum;       break;
     72                                case P1P1GLSEnum           : finiteelement = P1Enum;       break;
     73                                case MINIcondensedEnum     : finiteelement = P1bubbleEnum; break;
     74                                case MINIEnum              : finiteelement = P1bubbleEnum; break;
     75                                case TaylorHoodEnum        : finiteelement = P2Enum;       break;
     76                                case XTaylorHoodEnum       : finiteelement = P2Enum;       break;
     77                                case LATaylorHoodEnum      : finiteelement = P2Enum;       break;
     78                                case LACrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break;
     79                                case OneLayerP4zEnum       : finiteelement = P2xP4Enum;    break;
     80                                case CrouzeixRaviartEnum   : finiteelement = P2bubbleEnum; break;
     81                                default: _error_("finite element "<<EnumToStringx(finiteelement)<<" not supported");
     82                        }
     83                }
     84                else{
     85                        _error_("model not supported yet");
     86                }
     87
     88                if(isFS){
     89
     90                        /*Constraint at the bedrock interface (v.n = vz = 0) (Coordinates will be updated according to the bed slope)*/
     91                        iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
     92                        iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum);
     93                        iomodel->FetchData(&nodeonbase,NULL,NULL,MeshVertexonbaseEnum);
     94                        iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum);
     95                        if(iomodel->domaintype==Domain3DEnum){
     96                                iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
     97                        }
     98                        else if (iomodel->domaintype==Domain2DverticalEnum){
     99                                iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvyEnum);
     100                        }
     101                        else{
     102                                _error_("not supported yet");
     103                        }
     104                        if(iomodel->domaintype==Domain3DEnum){
     105                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0);
     106                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,1);
     107                                IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,2);
     108                                iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
     109                        }
     110                        else if (iomodel->domaintype==Domain2DverticalEnum){
     111                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0);
     112                                IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,1);
     113                                iomodel->DeleteData(spcvz,StressbalanceSpcvyEnum);
     114                        }
     115                        else{
     116                                _error_("not supported yet");
     117                        }
     118                        iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
     119                        iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum);
     120                        iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum);
     121                        iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum);
     122
     123                        /*Pressure spc*/
     124                        count = constraints->Size();
     125                        iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
     126                        iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum);
     127                        iomodel->FetchData(&z,NULL,NULL,MeshZEnum);
     128                        switch(finiteelement){
     129                                case P1Enum:
     130                                        for(i=0;i<iomodel->numberofvertices;i++){
     131                                                if(iomodel->my_vertices[i]){
     132                                                        if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
     133                                                                constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
     134                                                                count++;
     135                                                        }
     136                                                }
     137                                        }
     138                                        break;
     139                                case P1bubbleEnum:
     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,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
     144                                                                count++;
     145                                                        }
     146                                                }
     147                                        }
     148                                        break;
     149                                case P2Enum:
     150                                        for(i=0;i<iomodel->numberofvertices;i++){
     151                                                if(iomodel->my_vertices[i]){
     152                                                        if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
     153                                                                constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
     154                                                                count++;
     155                                                        }
     156                                                }
     157                                        }
     158                                        break;
     159                                case P2bubbleEnum:
     160                                        for(i=0;i<iomodel->numberofvertices;i++){
     161                                                if(iomodel->my_vertices[i]){
     162                                                        if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
     163                                                                constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+iomodel->numberoffaces+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
     164                                                                count++;
     165                                                        }
     166                                                }
     167                                        }
     168                                        break;
     169                                case P2xP4Enum:
     170                                        //Nothing yet
     171                                        break;
     172                                default:
     173                                        _error_("not implemented yet");
     174                        }
     175                        iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
     176                        iomodel->DeleteData(surface,SurfaceEnum);
     177                        iomodel->DeleteData(z,MeshZEnum);
     178                }
     179                else{
     180                        IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0);
     181                        if(iomodel->domaintype!=Domain2DverticalEnum){
     182                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,1);
     183                        }
     184                }
     185
     186                return;
     187        }
     188
     189        /*Constraints: fetch data: */
     190        iomodel->FetchData(&spcvx,&Mx,&Nx,StressbalanceSpcvxEnum);
     191        iomodel->FetchData(&spcvy,&My,&Ny,StressbalanceSpcvyEnum);
     192        iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
     193        iomodel->FetchData(&nodeonSSA,NULL,NULL,FlowequationBorderSSAEnum);
     194        if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonHO,NULL,NULL,FlowequationBorderHOEnum);
     195        if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum);
     196        if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonbase,NULL,NULL,MeshVertexonbaseEnum);
     197        if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum);
     198        iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
     199        iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum);
     200        iomodel->FetchData(&z,NULL,NULL,MeshZEnum);
     201
     202        /*Initialize counter: */
     203        count=0;
     204
     205        /*figure out times: */
     206        timesx=xNew<IssmDouble>(Nx);
     207        for(j=0;j<Nx;j++){
     208                timesx[j]=spcvx[(Mx-1)*Nx+j];
     209        }
     210        /*figure out times: */
     211        timesy=xNew<IssmDouble>(Ny);
     212        for(j=0;j<Ny;j++){
     213                timesy[j]=spcvy[(My-1)*Ny+j];
     214        }
     215        /*figure out times: */
     216        timesz=xNew<IssmDouble>(Nz);
     217        for(j=0;j<Nz;j++){
     218                timesz[j]=spcvz[(Mz-1)*Nz+j];
     219        }
     220
     221        /*Create spcs from x,y,z, as well as the spc values on those spcs: */
     222        for(i=0;i<iomodel->numberofvertices;i++){
     223                if(iomodel->my_vertices[i]){
     224
     225                        /*Start with adding spcs of coupling: zero at the border SSA/HO for the appropriate dofs*/
     226                        if(reCast<int,IssmDouble>(vertices_type[i]==SSAHOApproximationEnum)){
     227                                /*If grionSSA, spc HO dofs: 3 & 4*/
     228                                        if (reCast<int,IssmDouble>(nodeonHO[i])){
     229                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     230                                                count++;
     231                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     232                                                count++;
     233                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     234                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     235                                                        count++;
     236                                                }
     237                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     238                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     239                                                        count++;
     240                                                }
     241
     242                                        }
     243                                        else if (reCast<int,IssmDouble>(nodeonSSA[i])){
     244                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     245                                                count++;
     246                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     247                                                count++;
     248                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     249                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     250                                                        count++;
     251                                                }
     252                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     253                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     254                                                        count++;
     255                                                }
     256
     257                                        }
     258                                        else _error_("if vertices_type is SSAHO, you shoud have nodeonHO or nodeonSSA");
     259                        }
     260                        /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/
     261                        else if (reCast<int,IssmDouble>(vertices_type[i])==HOFSApproximationEnum){
     262                                /*If grion,HO spc FS dofs: 3 4 & 5*/
     263                                        if (reCast<int,IssmDouble>(nodeonHO[i])){
     264                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     265                                                count++;
     266                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     267                                                count++;
     268                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     269                                                count++;
     270                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     271                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     272                                                        count++;
     273                                                }
     274                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     275                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     276                                                        count++;
     277                                                }
     278
     279                                        }
     280                                        else if (reCast<int,IssmDouble>(nodeonFS[i])){ //spc HO nodes: 1 & 2
     281                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     282                                                count++;
     283                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     284                                                count++;
     285                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     286                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     287                                                        count++;
     288                                                }
     289                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     290                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     291                                                        count++;
     292                                                }
     293                                                if (!xIsNan<IssmDouble>(spcvz[i])){
     294                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     295                                                        count++;
     296                                                }
     297                                        }
     298                                        else _error_("if vertices_type is HOFS, you shoud have nodeonHO or nodeonFS");
     299                        }
     300                        /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/
     301                        else if (reCast<int,IssmDouble>(vertices_type[i])==SSAFSApproximationEnum){
     302                                /*If grion,HO spc FS dofs: 3 4 & 5*/
     303                                        if (reCast<int,IssmDouble>(nodeonSSA[i])){
     304                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     305                                                count++;
     306                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     307                                                count++;
     308                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     309                                                count++;
     310                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     311                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     312                                                        count++;
     313                                                }
     314                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     315                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     316                                                        count++;
     317                                                }
     318
     319                                        }
     320                                        else if (reCast<int,IssmDouble>(nodeonFS[i])){ //spc SSA nodes: 1 & 2
     321                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     322                                                count++;
     323                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     324                                                count++;
     325                                                if (!xIsNan<IssmDouble>(spcvx[i])){
     326                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     327                                                        count++;
     328                                                }
     329                                                if (!xIsNan<IssmDouble>(spcvy[i])){
     330                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     331                                                        count++;
     332                                                }
     333                                                if (!xIsNan<IssmDouble>(spcvz[i])){
     334                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     335                                                        count++;
     336                                                }
     337                                        }
     338                                        else _error_("if vertices_type is SSAFS, you shoud have nodeonSSA or nodeonFS");
     339                        }
     340                        /*Now add the regular spcs*/
     341                        else{
     342                                if (Mx==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvx[i])){
     343                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     344                                        count++;
     345
     346                                }
     347                                else if (Mx==iomodel->numberofvertices+1) {
     348                                        /*figure out times and values: */
     349                                        values=xNew<IssmDouble>(Nx);
     350                                        spcpresent=false;
     351                                        for(j=0;j<Nx;j++){
     352                                                values[j]=spcvx[i*Nx+j];
     353                                                if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
     354                                        }
     355
     356                                        if(spcpresent){
     357                                                constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,Nx,timesx,values,StressbalanceAnalysisEnum));
     358                                                count++;
     359                                        }
     360                                        xDelete<IssmDouble>(values);
     361                                }
     362                                else if (vertices_type[i]==SIAApproximationEnum){
     363                                        constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,StressbalanceAnalysisEnum));
     364                                        count++;
     365                                }
     366
     367                                if (My==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvy[i])){
     368                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy.
     369                                        count++;
     370                                }
     371                                else if (My==iomodel->numberofvertices+1){
     372                                        /*figure out times and values: */
     373                                        values=xNew<IssmDouble>(Ny);
     374                                        spcpresent=false;
     375                                        for(j=0;j<Ny;j++){
     376                                                values[j]=spcvy[i*Ny+j];
     377                                                if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
     378                                        }
     379                                        if(spcpresent){
     380                                                constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,Ny,timesy,values,StressbalanceAnalysisEnum));
     381                                                count++;
     382                                        }
     383                                        xDelete<IssmDouble>(values);
     384                                }
     385                                else if (vertices_type[i]==SIAApproximationEnum){
     386                                        constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,StressbalanceAnalysisEnum));
     387                                        count++;
     388                                }
     389
     390                                if (reCast<int,IssmDouble>(vertices_type[i])==FSApproximationEnum ||  (reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum)){
     391                                        if (Mz==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvz[i])){
     392                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
     393                                                count++;
     394                                        }
     395                                        else if (Mz==iomodel->numberofvertices+1){
     396                                                /*figure out times and values: */
     397                                                values=xNew<IssmDouble>(Nz);
     398                                                spcpresent=false;
     399                                                for(j=0;j<Nz;j++){
     400                                                        values[j]=spcvz[i*Nz+j];
     401                                                        if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
     402                                                }
     403                                                if(spcpresent){
     404                                                        constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,Nz,timesz,values,StressbalanceAnalysisEnum));
     405                                                        count++;
     406                                                }
     407                                                xDelete<IssmDouble>(values);
     408                                        }
     409
     410                                }
     411                                if (reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
     412                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
     413                                        count++;
     414                                }
     415                        }
     416                }
     417        }
     418
     419        /*Free data: */
     420        iomodel->DeleteData(spcvx,StressbalanceSpcvxEnum);
     421        iomodel->DeleteData(spcvy,StressbalanceSpcvyEnum);
     422        iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
     423        iomodel->DeleteData(nodeonSSA,FlowequationBorderSSAEnum);
     424        if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonHO,FlowequationBorderHOEnum);
     425        if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum);
     426        if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum);
     427        if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum);
     428        iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
     429        iomodel->DeleteData(surface,SurfaceEnum);
     430        iomodel->DeleteData(z,MeshZEnum);
     431
     432        /*Free resources:*/
     433        xDelete<IssmDouble>(timesx);
     434        xDelete<IssmDouble>(timesy);
     435        xDelete<IssmDouble>(timesz);
     436        xDelete<IssmDouble>(values);
     437
     438}/*}}}*/
     439void StressbalanceAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
     440
     441        /*Intermediary*/
     442        const int   RIFTINFOSIZE = 12;
     443        int         i;
     444        int         count;
     445        int         penpair_ids[2];
     446        bool        isSSA,isL1L2,isHO,isFS;
     447        int         numpenalties,numrifts,numriftsegments;
     448        IssmDouble *riftinfo       = NULL;
     449        IssmDouble *penalties      = NULL;
     450        int         assert_int;
     451
     452        /*Fetch parameters: */
     453        iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
     454        iomodel->Constant(&isFS,FlowequationIsFSEnum);
     455        iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
     456        iomodel->Constant(&isHO,FlowequationIsHOEnum);
     457        iomodel->Constant(&numrifts,RiftsNumriftsEnum);
     458
     459        /*Now, is the flag macayaealHO on? otherwise, do nothing: */
     460        if(!isSSA && !isHO && !isFS && !isL1L2) return;
     461
     462        /*Initialize counter: */
     463        count=0;
     464
     465        /*Create Penpair for penalties: */
     466        iomodel->FetchData(&penalties,&numpenalties,NULL,StressbalanceVertexPairingEnum);
     467
     468        for(i=0;i<numpenalties;i++){
     469
     470                if(iomodel->my_vertices[reCast<int,IssmDouble>(penalties[2*i+0]-1)]){
     471
     472                        /*In debugging mode, check that the second node is in the same cpu*/
     473                        assert_int=iomodel->my_vertices[reCast<int,IssmDouble>(penalties[2*i+1]-1)]; _assert_(assert_int);
     474
     475                        /*Get node ids*/
     476                        penpair_ids[0]=iomodel->nodecounter+reCast<int,IssmDouble>(penalties[2*i+0]);
     477                        penpair_ids[1]=iomodel->nodecounter+reCast<int,IssmDouble>(penalties[2*i+1]);
     478
     479                        /*Create Load*/
     480                        loads->AddObject(new Penpair(iomodel->loadcounter+count+1,&penpair_ids[0],StressbalanceAnalysisEnum));
     481                        count++;
     482                }
     483        }
     484
     485        /*free ressources: */
     486        iomodel->DeleteData(penalties,StressbalanceVertexPairingEnum);
     487
     488        /*Create Riffront loads for rifts: */
     489        if(numrifts){
     490                iomodel->FetchData(&riftinfo,&numriftsegments,NULL,RiftsRiftstructEnum);
     491                iomodel->FetchData(5,RiftsRiftstructEnum,ThicknessEnum,BaseEnum,SurfaceEnum,MaskGroundediceLevelsetEnum);
     492                for(i=0;i<numriftsegments;i++){
     493                        if(iomodel->my_elements[reCast<int,IssmDouble>(*(riftinfo+RIFTINFOSIZE*i+2))-1]){
     494                                loads->AddObject(new Riftfront(iomodel->loadcounter+count+1,i,iomodel,StressbalanceAnalysisEnum));
     495                                count++;
     496                        }
     497                }
     498                iomodel->DeleteData(5,RiftsRiftstructEnum,ThicknessEnum,BaseEnum,SurfaceEnum,MaskGroundediceLevelsetEnum);
     499                xDelete<IssmDouble>(riftinfo);
     500        }
     501}/*}}}*/
     502void StressbalanceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
     503
     504        /*Intermediary*/
     505        bool isSSA,isL1L2,isHO,isFS,iscoupling;
     506        int  finiteelement=-1,approximation=-1;
     507
     508        /*Fetch parameters: */
     509        iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
     510        iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
     511        iomodel->Constant(&isHO,FlowequationIsHOEnum);
     512        iomodel->Constant(&isFS,FlowequationIsFSEnum);
     513
     514        /*Now, check that we have non SIA elements */
     515        if(!isSSA & !isL1L2 & !isHO & !isFS) return;
     516
     517        /*Do we have coupling*/
     518        if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     519         iscoupling = true;
     520        else
     521         iscoupling = false;
     522
     523        /*If no coupling, call Regular CreateNodes, else, use P1 elements only*/
     524        if(!iscoupling){
     525
     526                /*Get finite element type*/
     527                if(isSSA){
     528                        approximation=SSAApproximationEnum;
     529                        iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
     530                }
     531                else if(isL1L2){
     532                        approximation = L1L2ApproximationEnum;
     533                        finiteelement = P1Enum;
     534                }
     535                else if(isHO){
     536                        approximation = HOApproximationEnum;
     537                        iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
     538                }
     539                else if(isFS){
     540                        approximation = FSApproximationEnum;
     541                        iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
     542                }
     543                iomodel->FetchData(3,FlowequationBorderSSAEnum,FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
     544                if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(3,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderFSEnum);
     545                ::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,finiteelement,approximation);
     546                iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
     547                                        FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
     548        }
     549        else{
     550                /*Coupling: we are going to create P1 Elements only*/
     551
     552                Node*  node  = NULL;
     553                int    lid=0;
     554                if(!nodes) nodes = new Nodes();
     555
     556                iomodel->FetchData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
     557                                        FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
     558                if(isFS){
     559                        /*P1+ velocity*/
     560                        for(int i=0;i<iomodel->numberofvertices;i++){
     561                                if(iomodel->my_vertices[i]){
     562                                        approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
     563                                        if(approximation==FSApproximationEnum)  approximation=FSvelocityEnum;
     564                                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,approximation));
     565                                }
     566                        }
     567                        for(int i=0;i<iomodel->numberofelements;i++){
     568                                if(iomodel->my_elements[i]){
     569                                        node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum);
     570                                        node->Deactivate();
     571                                        nodes->AddObject(node);
     572                                }
     573                        }
     574                        /*P1 pressure*/
     575                        for(int i=0;i<iomodel->numberofvertices;i++){
     576                                if(iomodel->my_vertices[i]){
     577                                        approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
     578                                        node = new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,i,iomodel,StressbalanceAnalysisEnum,FSpressureEnum);
     579                                        if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){
     580                                                node->Deactivate();
     581                                        }
     582                                        nodes->AddObject(node);
     583                                }
     584                        }
     585                }
     586                else{
     587                        for(int i=0;i<iomodel->numberofvertices;i++){
     588                                if(iomodel->my_vertices[i]){
     589                                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i])));
     590                                }
     591                        }
     592                }
     593                iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
     594                                        FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
     595        }
     596}/*}}}*/
    13597int  StressbalanceAnalysis::DofsPerNode(int** pdoftype,int domaintype,int approximation){/*{{{*/
    14598
    … …  
    84668        return numdofs;
    85669}/*}}}*/
    86 void StressbalanceAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    87 
    88         /*Intermediaries*/
    89         int     fe_FS;
    90         int     numoutputs;
    91         char**  requestedoutputs = NULL;
    92         int     materials_type;
    93 
    94         parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum));
    95         parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSSAEnum));
    96         parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsL1L2Enum));
    97         parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsHOEnum));
    98         parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
    99         parameters->AddObject(iomodel->CopyConstantObject(FlowequationFeFSEnum));
    100         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
    101         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
    102         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
    103         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum));
    104         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
    105         parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum));
    106         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
    107         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceFSreconditioningEnum));
    108         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum));
    109         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum));
    110         parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum));
    111 
    112         /*XTH LATH parameters*/
    113         iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
    114         if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
    115                 parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum));
    116                 parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum));
    117                 parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum));
    118         }
    119 
    120         iomodel->Constant(&materials_type,MaterialsEnum);
    121         if(materials_type==MatdamageiceEnum){
    122                 parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum));
    123                 parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum));
    124         }
    125 
    126         /*Requested outputs*/
    127         iomodel->FetchData(&requestedoutputs,&numoutputs,StressbalanceRequestedOutputsEnum);
    128         parameters->AddObject(new IntParam(StressbalanceNumRequestedOutputsEnum,numoutputs));
    129         if(numoutputs)parameters->AddObject(new StringArrayParam(StressbalanceRequestedOutputsEnum,requestedoutputs,numoutputs));
    130         iomodel->DeleteData(&requestedoutputs,numoutputs,StressbalanceRequestedOutputsEnum);
    131 
    132         /*Deal with friction parameters*/
    133         int frictionlaw;
    134         iomodel->Constant(&frictionlaw,FrictionLawEnum);
    135         if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject(FrictionGammaEnum));
    136 
    137 }/*}}}*/
    138670void StressbalanceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    139671
    … …  
    304836        xDelete<int>(finiteelement_list);
    305837}/*}}}*/
    306 void StressbalanceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    307 
    308         /*Intermediary*/
    309         bool isSSA,isL1L2,isHO,isFS,iscoupling;
    310         int  finiteelement=-1,approximation=-1;
    311 
    312         /*Fetch parameters: */
    313         iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
    314         iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
    315         iomodel->Constant(&isHO,FlowequationIsHOEnum);
    316         iomodel->Constant(&isFS,FlowequationIsFSEnum);
    317 
    318         /*Now, check that we have non SIA elements */
    319         if(!isSSA & !isL1L2 & !isHO & !isFS) return;
    320 
    321         /*Do we have coupling*/
    322         if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
    323          iscoupling = true;
    324         else
    325          iscoupling = false;
    326 
    327         /*If no coupling, call Regular CreateNodes, else, use P1 elements only*/
    328         if(!iscoupling){
    329 
    330                 /*Get finite element type*/
    331                 if(isSSA){
    332                         approximation=SSAApproximationEnum;
    333                         iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
    334                 }
    335                 else if(isL1L2){
    336                         approximation = L1L2ApproximationEnum;
    337                         finiteelement = P1Enum;
    338                 }
    339                 else if(isHO){
    340                         approximation = HOApproximationEnum;
    341                         iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
    342                 }
    343                 else if(isFS){
    344                         approximation = FSApproximationEnum;
    345                         iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
    346                 }
    347                 iomodel->FetchData(3,FlowequationBorderSSAEnum,FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
    348                 if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(3,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderFSEnum);
    349                 ::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,finiteelement,approximation);
    350                 iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
    351                                         FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
    352         }
    353         else{
    354                 /*Coupling: we are going to create P1 Elements only*/
    355 
    356                 Node*  node  = NULL;
    357                 int    lid=0;
    358                 if(!nodes) nodes = new Nodes();
    359 
    360                 iomodel->FetchData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
    361                                         FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
    362                 if(isFS){
    363                         /*P1+ velocity*/
    364                         for(int i=0;i<iomodel->numberofvertices;i++){
    365                                 if(iomodel->my_vertices[i]){
    366                                         approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
    367                                         if(approximation==FSApproximationEnum)  approximation=FSvelocityEnum;
    368                                         nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,approximation));
    369                                 }
    370                         }
    371                         for(int i=0;i<iomodel->numberofelements;i++){
    372                                 if(iomodel->my_elements[i]){
    373                                         node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum);
    374                                         node->Deactivate();
    375                                         nodes->AddObject(node);
    376                                 }
    377                         }
    378                         /*P1 pressure*/
    379                         for(int i=0;i<iomodel->numberofvertices;i++){
    380                                 if(iomodel->my_vertices[i]){
    381                                         approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
    382                                         node = new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,i,iomodel,StressbalanceAnalysisEnum,FSpressureEnum);
    383                                         if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){
    384                                                 node->Deactivate();
    385                                         }
    386                                         nodes->AddObject(node);
    387                                 }
    388                         }
    389                 }
    390                 else{
    391                         for(int i=0;i<iomodel->numberofvertices;i++){
    392                                 if(iomodel->my_vertices[i]){
    393                                         nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i])));
    394                                 }
    395                         }
    396                 }
    397                 iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
    398                                         FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
    399         }
    400 }/*}}}*/
    401 void StressbalanceAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    402 
    403         /*Intermediary*/
    404         int        i,j;
    405         int        count,finiteelement;
    406         IssmDouble g;
    407         IssmDouble rho_ice;
    408         IssmDouble FSreconditioning;
    409         bool       isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
    410         bool       spcpresent = false;
    411         int        Mx,Nx;
    412         int        My,Ny;
    413         int        Mz,Nz;
    414         IssmDouble *spcvx          = NULL;
    415         IssmDouble *spcvy          = NULL;
    416         IssmDouble *spcvz          = NULL;
    417         IssmDouble *nodeonSSA = NULL;
    418         IssmDouble *nodeonHO   = NULL;
    419         IssmDouble *nodeonFS   = NULL;
    420         IssmDouble *nodeonbase      = NULL;
    421         IssmDouble *groundedice_ls = NULL;
    422         IssmDouble *vertices_type  = NULL;
    423         IssmDouble *surface        = NULL;
    424         IssmDouble *z              = NULL;
    425         IssmDouble *timesx=NULL;
    426         IssmDouble *timesy=NULL;
    427         IssmDouble *timesz=NULL;
    428    IssmDouble* values=NULL;
    429 
    430         /*Fetch parameters: */
    431         iomodel->Constant(&g,ConstantsGEnum);
    432         iomodel->Constant(&rho_ice,MaterialsRhoIceEnum);
    433         iomodel->Constant(&FSreconditioning,StressbalanceFSreconditioningEnum);
    434         iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
    435         iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
    436         iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
    437         iomodel->Constant(&isHO,FlowequationIsHOEnum);
    438         iomodel->Constant(&isFS,FlowequationIsFSEnum);
    439 
    440         /*Now, is the flag macayaealHO on? otherwise, do nothing: */
    441         if(!isSSA && !isHO && !isFS && !isL1L2) return;
    442 
    443         /*Do we have coupling*/
    444         if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
    445          iscoupling = true;
    446         else
    447          iscoupling = false;
    448 
    449         /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
    450         if(!iscoupling){
    451 
    452                 /*Get finite element type*/
    453                 if(isSSA)       iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
    454                 else if(isL1L2) finiteelement = P1Enum;
    455                 else if(isHO)   iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
    456                 else if(isFS){  iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
    457                         /*Deduce velocity interpolation from finite element*/
    458                         switch(finiteelement){
    459                                 case P1P1Enum              : finiteelement = P1Enum;       break;
    460                                 case P1P1GLSEnum           : finiteelement = P1Enum;       break;
    461                                 case MINIcondensedEnum     : finiteelement = P1bubbleEnum; break;
    462                                 case MINIEnum              : finiteelement = P1bubbleEnum; break;
    463                                 case TaylorHoodEnum        : finiteelement = P2Enum;       break;
    464                                 case XTaylorHoodEnum       : finiteelement = P2Enum;       break;
    465                                 case LATaylorHoodEnum      : finiteelement = P2Enum;       break;
    466                                 case LACrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break;
    467                                 case OneLayerP4zEnum       : finiteelement = P2xP4Enum;    break;
    468                                 case CrouzeixRaviartEnum   : finiteelement = P2bubbleEnum; break;
    469                                 default: _error_("finite element "<<EnumToStringx(finiteelement)<<" not supported");
    470                         }
    471                 }
    472                 else{
    473                         _error_("model not supported yet");
    474                 }
    475 
    476                 if(isFS){
    477 
    478                         /*Constraint at the bedrock interface (v.n = vz = 0) (Coordinates will be updated according to the bed slope)*/
    479                         iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
    480                         iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum);
    481                         iomodel->FetchData(&nodeonbase,NULL,NULL,MeshVertexonbaseEnum);
    482                         iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum);
    483                         if(iomodel->domaintype==Domain3DEnum){
    484                                 iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
    485                         }
    486                         else if (iomodel->domaintype==Domain2DverticalEnum){
    487                                 iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvyEnum);
    488                         }
    489                         else{
    490                                 _error_("not supported yet");
    491                         }
    492                         if(iomodel->domaintype==Domain3DEnum){
    493                                 IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0);
    494                                 IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,1);
    495                                 IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,2);
    496                                 iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
    497                         }
    498                         else if (iomodel->domaintype==Domain2DverticalEnum){
    499                                 IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0);
    500                                 IoModelToConstraintsx(constraints,iomodel,spcvz,Mz,Nz,StressbalanceAnalysisEnum,finiteelement,1);
    501                                 iomodel->DeleteData(spcvz,StressbalanceSpcvyEnum);
    502                         }
    503                         else{
    504                                 _error_("not supported yet");
    505                         }
    506                         iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
    507                         iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum);
    508                         iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum);
    509                         iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum);
    510 
    511                         /*Pressure spc*/
    512                         count = constraints->Size();
    513                         iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
    514                         iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum);
    515                         iomodel->FetchData(&z,NULL,NULL,MeshZEnum);
    516                         switch(finiteelement){
    517                                 case P1Enum:
    518                                         for(i=0;i<iomodel->numberofvertices;i++){
    519                                                 if(iomodel->my_vertices[i]){
    520                                                         if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
    521                                                                 constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
    522                                                                 count++;
    523                                                         }
    524                                                 }
    525                                         }
    526                                         break;
    527                                 case P1bubbleEnum:
    528                                         for(i=0;i<iomodel->numberofvertices;i++){
    529                                                 if(iomodel->my_vertices[i]){
    530                                                         if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
    531                                                                 constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
    532                                                                 count++;
    533                                                         }
    534                                                 }
    535                                         }
    536                                         break;
    537                                 case P2Enum:
    538                                         for(i=0;i<iomodel->numberofvertices;i++){
    539                                                 if(iomodel->my_vertices[i]){
    540                                                         if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
    541                                                                 constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
    542                                                                 count++;
    543                                                         }
    544                                                 }
    545                                         }
    546                                         break;
    547                                 case P2bubbleEnum:
    548                                         for(i=0;i<iomodel->numberofvertices;i++){
    549                                                 if(iomodel->my_vertices[i]){
    550                                                         if(reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
    551                                                                 constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+iomodel->numberoffaces+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum));
    552                                                                 count++;
    553                                                         }
    554                                                 }
    555                                         }
    556                                         break;
    557                                 case P2xP4Enum:
    558                                         //Nothing yet
    559                                         break;
    560                                 default:
    561                                         _error_("not implemented yet");
    562                         }
    563                         iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
    564                         iomodel->DeleteData(surface,SurfaceEnum);
    565                         iomodel->DeleteData(z,MeshZEnum);
    566                 }
    567                 else{
    568                         IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,0);
    569                         if(iomodel->domaintype!=Domain2DverticalEnum){
    570                                 IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,1);
    571                         }
    572                 }
    573 
    574                 return;
    575         }
    576 
    577         /*Constraints: fetch data: */
    578         iomodel->FetchData(&spcvx,&Mx,&Nx,StressbalanceSpcvxEnum);
    579         iomodel->FetchData(&spcvy,&My,&Ny,StressbalanceSpcvyEnum);
    580         iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
    581         iomodel->FetchData(&nodeonSSA,NULL,NULL,FlowequationBorderSSAEnum);
    582         if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonHO,NULL,NULL,FlowequationBorderHOEnum);
    583         if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonFS,NULL,NULL,FlowequationBorderFSEnum);
    584         if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&nodeonbase,NULL,NULL,MeshVertexonbaseEnum);
    585         if(iomodel->domaintype==Domain3DEnum)iomodel->FetchData(&groundedice_ls,NULL,NULL,MaskGroundediceLevelsetEnum);
    586         iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
    587         iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum);
    588         iomodel->FetchData(&z,NULL,NULL,MeshZEnum);
    589 
    590         /*Initialize counter: */
    591         count=0;
    592 
    593         /*figure out times: */
    594         timesx=xNew<IssmDouble>(Nx);
    595         for(j=0;j<Nx;j++){
    596                 timesx[j]=spcvx[(Mx-1)*Nx+j];
    597         }
    598         /*figure out times: */
    599         timesy=xNew<IssmDouble>(Ny);
    600         for(j=0;j<Ny;j++){
    601                 timesy[j]=spcvy[(My-1)*Ny+j];
    602         }
    603         /*figure out times: */
    604         timesz=xNew<IssmDouble>(Nz);
    605         for(j=0;j<Nz;j++){
    606                 timesz[j]=spcvz[(Mz-1)*Nz+j];
    607         }
    608 
    609         /*Create spcs from x,y,z, as well as the spc values on those spcs: */
    610         for(i=0;i<iomodel->numberofvertices;i++){
    611                 if(iomodel->my_vertices[i]){
    612 
    613                         /*Start with adding spcs of coupling: zero at the border SSA/HO for the appropriate dofs*/
    614                         if(reCast<int,IssmDouble>(vertices_type[i]==SSAHOApproximationEnum)){
    615                                 /*If grionSSA, spc HO dofs: 3 & 4*/
    616                                         if (reCast<int,IssmDouble>(nodeonHO[i])){
    617                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    618                                                 count++;
    619                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    620                                                 count++;
    621                                                 if (!xIsNan<IssmDouble>(spcvx[i])){
    622                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    623                                                         count++;
    624                                                 }
    625                                                 if (!xIsNan<IssmDouble>(spcvy[i])){
    626                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    627                                                         count++;
    628                                                 }
    629 
    630                                         }
    631                                         else if (reCast<int,IssmDouble>(nodeonSSA[i])){
    632                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    633                                                 count++;
    634                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    635                                                 count++;
    636                                                 if (!xIsNan<IssmDouble>(spcvx[i])){
    637                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    638                                                         count++;
    639                                                 }
    640                                                 if (!xIsNan<IssmDouble>(spcvy[i])){
    641                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    642                                                         count++;
    643                                                 }
    644 
    645                                         }
    646                                         else _error_("if vertices_type is SSAHO, you shoud have nodeonHO or nodeonSSA");
    647                         }
    648                         /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/
    649                         else if (reCast<int,IssmDouble>(vertices_type[i])==HOFSApproximationEnum){
    650                                 /*If grion,HO spc FS dofs: 3 4 & 5*/
    651                                         if (reCast<int,IssmDouble>(nodeonHO[i])){
    652                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    653                                                 count++;
    654                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    655                                                 count++;
    656                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    657                                                 count++;
    658                                                 if (!xIsNan<IssmDouble>(spcvx[i])){
    659                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    660                                                         count++;
    661                                                 }
    662                                                 if (!xIsNan<IssmDouble>(spcvy[i])){
    663                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    664                                                         count++;
    665                                                 }
    666 
    667                                         }
    668                                         else if (reCast<int,IssmDouble>(nodeonFS[i])){ //spc HO nodes: 1 & 2
    669                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    670                                                 count++;
    671                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    672                                                 count++;
    673                                                 if (!xIsNan<IssmDouble>(spcvx[i])){
    674                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    675                                                         count++;
    676                                                 }
    677                                                 if (!xIsNan<IssmDouble>(spcvy[i])){
    678                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    679                                                         count++;
    680                                                 }
    681                                                 if (!xIsNan<IssmDouble>(spcvz[i])){
    682                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    683                                                         count++;
    684                                                 }
    685                                         }
    686                                         else _error_("if vertices_type is HOFS, you shoud have nodeonHO or nodeonFS");
    687                         }
    688                         /*Also add spcs of coupling: zero at the border HO/FS for the appropriate dofs*/
    689                         else if (reCast<int,IssmDouble>(vertices_type[i])==SSAFSApproximationEnum){
    690                                 /*If grion,HO spc FS dofs: 3 4 & 5*/
    691                                         if (reCast<int,IssmDouble>(nodeonSSA[i])){
    692                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    693                                                 count++;
    694                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    695                                                 count++;
    696                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    697                                                 count++;
    698                                                 if (!xIsNan<IssmDouble>(spcvx[i])){
    699                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    700                                                         count++;
    701                                                 }
    702                                                 if (!xIsNan<IssmDouble>(spcvy[i])){
    703                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    704                                                         count++;
    705                                                 }
    706 
    707                                         }
    708                                         else if (reCast<int,IssmDouble>(nodeonFS[i])){ //spc SSA nodes: 1 & 2
    709                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    710                                                 count++;
    711                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    712                                                 count++;
    713                                                 if (!xIsNan<IssmDouble>(spcvx[i])){
    714                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    715                                                         count++;
    716                                                 }
    717                                                 if (!xIsNan<IssmDouble>(spcvy[i])){
    718                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    719                                                         count++;
    720                                                 }
    721                                                 if (!xIsNan<IssmDouble>(spcvz[i])){
    722                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    723                                                         count++;
    724                                                 }
    725                                         }
    726                                         else _error_("if vertices_type is SSAFS, you shoud have nodeonSSA or nodeonFS");
    727                         }
    728                         /*Now add the regular spcs*/
    729                         else{
    730                                 if (Mx==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvx[i])){
    731                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    732                                         count++;
    733 
    734                                 }
    735                                 else if (Mx==iomodel->numberofvertices+1) {
    736                                         /*figure out times and values: */
    737                                         values=xNew<IssmDouble>(Nx);
    738                                         spcpresent=false;
    739                                         for(j=0;j<Nx;j++){
    740                                                 values[j]=spcvx[i*Nx+j];
    741                                                 if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
    742                                         }
    743 
    744                                         if(spcpresent){
    745                                                 constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,Nx,timesx,values,StressbalanceAnalysisEnum));
    746                                                 count++;
    747                                         }
    748                                         xDelete<IssmDouble>(values);
    749                                 }
    750                                 else if (vertices_type[i]==SIAApproximationEnum){
    751                                         constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,StressbalanceAnalysisEnum));
    752                                         count++;
    753                                 }
    754 
    755                                 if (My==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvy[i])){
    756                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy.
    757                                         count++;
    758                                 }
    759                                 else if (My==iomodel->numberofvertices+1){
    760                                         /*figure out times and values: */
    761                                         values=xNew<IssmDouble>(Ny);
    762                                         spcpresent=false;
    763                                         for(j=0;j<Ny;j++){
    764                                                 values[j]=spcvy[i*Ny+j];
    765                                                 if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
    766                                         }
    767                                         if(spcpresent){
    768                                                 constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,Ny,timesy,values,StressbalanceAnalysisEnum));
    769                                                 count++;
    770                                         }
    771                                         xDelete<IssmDouble>(values);
    772                                 }
    773                                 else if (vertices_type[i]==SIAApproximationEnum){
    774                                         constraints->AddObject(new SpcDynamic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,StressbalanceAnalysisEnum));
    775                                         count++;
    776                                 }
    777 
    778                                 if (reCast<int,IssmDouble>(vertices_type[i])==FSApproximationEnum ||  (reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum)){
    779                                         if (Mz==iomodel->numberofvertices && !xIsNan<IssmDouble>(spcvz[i])){
    780                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvz[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
    781                                                 count++;
    782                                         }
    783                                         else if (Mz==iomodel->numberofvertices+1){
    784                                                 /*figure out times and values: */
    785                                                 values=xNew<IssmDouble>(Nz);
    786                                                 spcpresent=false;
    787                                                 for(j=0;j<Nz;j++){
    788                                                         values[j]=spcvz[i*Nz+j];
    789                                                         if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
    790                                                 }
    791                                                 if(spcpresent){
    792                                                         constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,Nz,timesz,values,StressbalanceAnalysisEnum));
    793                                                         count++;
    794                                                 }
    795                                                 xDelete<IssmDouble>(values);
    796                                         }
    797 
    798                                 }
    799                                 if (reCast<int,IssmDouble>(vertices_type[i])==NoneApproximationEnum){
    800                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
    801                                         count++;
    802                                 }
    803                         }
    804                 }
    805         }
    806 
    807         /*Free data: */
    808         iomodel->DeleteData(spcvx,StressbalanceSpcvxEnum);
    809         iomodel->DeleteData(spcvy,StressbalanceSpcvyEnum);
    810         iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
    811         iomodel->DeleteData(nodeonSSA,FlowequationBorderSSAEnum);
    812         if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonHO,FlowequationBorderHOEnum);
    813         if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonFS,FlowequationBorderFSEnum);
    814         if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum);
    815         if(iomodel->domaintype==Domain3DEnum)iomodel->DeleteData(groundedice_ls,MaskGroundediceLevelsetEnum);
    816         iomodel->DeleteData(vertices_type,FlowequationVertexEquationEnum);
    817         iomodel->DeleteData(surface,SurfaceEnum);
    818         iomodel->DeleteData(z,MeshZEnum);
    819 
    820         /*Free resources:*/
    821         xDelete<IssmDouble>(timesx);
    822         xDelete<IssmDouble>(timesy);
    823         xDelete<IssmDouble>(timesz);
    824         xDelete<IssmDouble>(values);
    825 
    826 }/*}}}*/
    827 void StressbalanceAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
    828 
    829         /*Intermediary*/
    830         const int   RIFTINFOSIZE = 12;
    831         int         i;
    832         int         count;
    833         int         penpair_ids[2];
    834         bool        isSSA,isL1L2,isHO,isFS;
    835         int         numpenalties,numrifts,numriftsegments;
    836         IssmDouble *riftinfo       = NULL;
    837         IssmDouble *penalties      = NULL;
    838         int         assert_int;
    839 
    840         /*Fetch parameters: */
    841         iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
    842         iomodel->Constant(&isFS,FlowequationIsFSEnum);
    843         iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
    844         iomodel->Constant(&isHO,FlowequationIsHOEnum);
    845         iomodel->Constant(&numrifts,RiftsNumriftsEnum);
    846 
    847         /*Now, is the flag macayaealHO on? otherwise, do nothing: */
    848         if(!isSSA && !isHO && !isFS && !isL1L2) return;
    849 
    850         /*Initialize counter: */
    851         count=0;
    852 
    853         /*Create Penpair for penalties: */
    854         iomodel->FetchData(&penalties,&numpenalties,NULL,StressbalanceVertexPairingEnum);
    855 
    856         for(i=0;i<numpenalties;i++){
    857 
    858                 if(iomodel->my_vertices[reCast<int,IssmDouble>(penalties[2*i+0]-1)]){
    859 
    860                         /*In debugging mode, check that the second node is in the same cpu*/
    861                         assert_int=iomodel->my_vertices[reCast<int,IssmDouble>(penalties[2*i+1]-1)]; _assert_(assert_int);
    862 
    863                         /*Get node ids*/
    864                         penpair_ids[0]=iomodel->nodecounter+reCast<int,IssmDouble>(penalties[2*i+0]);
    865                         penpair_ids[1]=iomodel->nodecounter+reCast<int,IssmDouble>(penalties[2*i+1]);
    866 
    867                         /*Create Load*/
    868                         loads->AddObject(new Penpair(iomodel->loadcounter+count+1,&penpair_ids[0],StressbalanceAnalysisEnum));
    869                         count++;
    870                 }
    871         }
    872 
    873         /*free ressources: */
    874         iomodel->DeleteData(penalties,StressbalanceVertexPairingEnum);
    875 
    876         /*Create Riffront loads for rifts: */
    877         if(numrifts){
    878                 iomodel->FetchData(&riftinfo,&numriftsegments,NULL,RiftsRiftstructEnum);
    879                 iomodel->FetchData(5,RiftsRiftstructEnum,ThicknessEnum,BaseEnum,SurfaceEnum,MaskGroundediceLevelsetEnum);
    880                 for(i=0;i<numriftsegments;i++){
    881                         if(iomodel->my_elements[reCast<int,IssmDouble>(*(riftinfo+RIFTINFOSIZE*i+2))-1]){
    882                                 loads->AddObject(new Riftfront(iomodel->loadcounter+count+1,i,iomodel,StressbalanceAnalysisEnum));
    883                                 count++;
    884                         }
    885                 }
    886                 iomodel->DeleteData(5,RiftsRiftstructEnum,ThicknessEnum,BaseEnum,SurfaceEnum,MaskGroundediceLevelsetEnum);
    887                 xDelete<IssmDouble>(riftinfo);
    888         }
     838void StressbalanceAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     839
     840        /*Intermediaries*/
     841        int     fe_FS;
     842        int     numoutputs;
     843        char**  requestedoutputs = NULL;
     844        int     materials_type;
     845
     846        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum));
     847        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSSAEnum));
     848        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsL1L2Enum));
     849        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsHOEnum));
     850        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
     851        parameters->AddObject(iomodel->CopyConstantObject(FlowequationFeFSEnum));
     852        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
     853        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
     854        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
     855        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum));
     856        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
     857        parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum));
     858        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
     859        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceFSreconditioningEnum));
     860        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum));
     861        parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum));
     862        parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum));
     863
     864        /*XTH LATH parameters*/
     865        iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
     866        if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
     867                parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum));
     868                parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum));
     869                parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum));
     870        }
     871
     872        iomodel->Constant(&materials_type,MaterialsEnum);
     873        if(materials_type==MatdamageiceEnum){
     874                parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum));
     875                parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum));
     876        }
     877
     878        /*Requested outputs*/
     879        iomodel->FetchData(&requestedoutputs,&numoutputs,StressbalanceRequestedOutputsEnum);
     880        parameters->AddObject(new IntParam(StressbalanceNumRequestedOutputsEnum,numoutputs));
     881        if(numoutputs)parameters->AddObject(new StringArrayParam(StressbalanceRequestedOutputsEnum,requestedoutputs,numoutputs));
     882        iomodel->DeleteData(&requestedoutputs,numoutputs,StressbalanceRequestedOutputsEnum);
     883
     884        /*Deal with friction parameters*/
     885        int frictionlaw;
     886        iomodel->Constant(&frictionlaw,FrictionLawEnum);
     887        if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject(FrictionGammaEnum));
     888
    889889}/*}}}*/
    890890
    … …  
    10271027        }
    10281028}/*}}}*/
    1029 void StressbalanceAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
     1029void           StressbalanceAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    10301030
    10311031        int approximation;
    … …  
    10481048        }
    10491049}/*}}}*/
    1050 void StressbalanceAnalysis::GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element){/*{{{*/
     1050void           StressbalanceAnalysis::GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    10511051
    10521052        IssmDouble   vx,vy;
    … …  
    10981098        xDelete<int>(doflist);
    10991099}/*}}}*/
    1100 void StressbalanceAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     1100void           StressbalanceAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
    11011101        _error_("Not implemented yet");
    11021102}/*}}}*/
    1103 void StressbalanceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
     1103void           StressbalanceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    11041104
    11051105        int approximation;
    … …  
    11331133        }
    11341134}/*}}}*/
    1135 void StressbalanceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
     1135void           StressbalanceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    11361136        /*Default, do nothing*/
    11371137        bool islevelset;
    … …  
    16601660        return pe;
    16611661}/*}}}*/
    1662 void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1662void           StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    16631663        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    16641664         * For node i, Bi can be expressed in the actual coordinate system
    … …  
    17001700        xDelete<IssmDouble>(dbasis);
    17011701}/*}}}*/
    1702 void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1702void           StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    17031703        /*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2.
    17041704         * For node i, Bi can be expressed in the actual coordinate system
    … …  
    17371737        xDelete<IssmDouble>(basis);
    17381738}/*}}}*/
    1739 void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1739void           StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    17401740        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    17411741         * For node i, Bi' can be expressed in the actual coordinate system
    … …  
    17771777        xDelete<IssmDouble>(dbasis);
    17781778}/*}}}*/
    1779 void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
     1779void           StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
    17801780
    17811781        int         i,dim,domaintype;
    … …  
    21112111        return pe;
    21122112}/*}}}*/
    2113 void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
     2113void           StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
    21142114
    21152115        int         i,dim,domaintype;
    … …  
    22862286        return Ke;
    22872287}/*}}}*/
     2288ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFriction(Element* element){/*{{{*/
     2289
     2290        /* Check if ice in element */
     2291        if(!element->IsIceInElement()) return NULL;
     2292
     2293        if(element->IsFloating() || !element->IsOnBase()) return NULL;
     2294
     2295        /*Intermediaries*/
     2296        int         dim;
     2297        bool        mainlyfloating;
     2298        int         migration_style,point1;
     2299        IssmDouble  alpha2,Jdet,fraction1,fraction2;
     2300        IssmDouble  gllevelset,phi=1.;
     2301        IssmDouble *xyz_list_base = NULL;
     2302        Gauss*      gauss         = NULL;
     2303
     2304        /*Get problem dimension*/
     2305        element->FindParam(&dim,DomainDimensionEnum);
     2306
     2307        /*Fetch number of nodes and dof for this finite element*/
     2308        int numnodes = element->GetNumberOfNodes();
     2309        int numdof   = numnodes*(dim-1);
     2310
     2311        /*Initialize Element matrix and vectors*/
     2312        ElementMatrix* Ke = element->NewElementMatrix(HOApproximationEnum);
     2313        IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
     2314        IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
     2315
     2316        /*Retrieve all inputs and parameters*/
     2317        element->GetVerticesCoordinatesBase(&xyz_list_base);
     2318        element->FindParam(&migration_style,GroundinglineMigrationEnum);
     2319        Input* gllevelset_input = NULL;
     2320
     2321        /*build friction object, used later on: */
     2322        Friction* friction=new Friction(element,2);
     2323
     2324        /*Recover portion of element that is grounded*/
     2325        if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base);
     2326        if(migration_style==SubelementMigration2Enum){
     2327                gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
     2328                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     2329                //gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
     2330                gauss=element->NewGaussBase(2);
     2331        }
     2332        else{
     2333                gauss=element->NewGaussBase(2);
     2334        }
     2335
     2336        /* Start  looping on the number of gaussian points: */
     2337        for(int ig=gauss->begin();ig<gauss->end();ig++){
     2338                gauss->GaussPoint(ig);
     2339
     2340                friction->GetAlpha2(&alpha2,gauss);
     2341                if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
     2342                if(migration_style==SubelementMigration2Enum){
     2343                        gllevelset_input->GetInputValue(&gllevelset, gauss);
     2344                        if(gllevelset<0.) alpha2=0.;
     2345                }
     2346
     2347                this->GetBHOFriction(B,element,dim,xyz_list_base,gauss);
     2348                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     2349                for(int i=0;i<dim-1;i++) D[i*(dim-1)+i]=alpha2*gauss->weight*Jdet;
     2350
     2351                TripleMultiply(B,dim-1,numdof,1,
     2352                                        D,dim-1,dim-1,0,
     2353                                        B,dim-1,numdof,0,
     2354                                        &Ke->values[0],1);
     2355        }
     2356
     2357        /*Transform Coordinate System*/
     2358        if(dim==3) element->TransformStiffnessMatrixCoord(Ke,XYEnum);
     2359
     2360        /*Clean up and return*/
     2361        delete gauss;
     2362        delete friction;
     2363        xDelete<IssmDouble>(xyz_list_base);
     2364        xDelete<IssmDouble>(B);
     2365        xDelete<IssmDouble>(D);
     2366        return Ke;
     2367}/*}}}*/
    22882368ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOViscous(Element* element){/*{{{*/
    22892369
    … …  
    23582438        return Ke;
    23592439}/*}}}*/
    2360 ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFriction(Element* element){/*{{{*/
    2361 
    2362         /* Check if ice in element */
    2363         if(!element->IsIceInElement()) return NULL;
    2364 
    2365         if(element->IsFloating() || !element->IsOnBase()) return NULL;
    2366 
    2367         /*Intermediaries*/
    2368         int         dim;
    2369         bool        mainlyfloating;
    2370         int         migration_style,point1;
    2371         IssmDouble  alpha2,Jdet,fraction1,fraction2;
    2372         IssmDouble  gllevelset,phi=1.;
    2373         IssmDouble *xyz_list_base = NULL;
    2374         Gauss*      gauss         = NULL;
    2375 
    2376         /*Get problem dimension*/
    2377         element->FindParam(&dim,DomainDimensionEnum);
    2378 
    2379         /*Fetch number of nodes and dof for this finite element*/
    2380         int numnodes = element->GetNumberOfNodes();
    2381         int numdof   = numnodes*(dim-1);
    2382 
    2383         /*Initialize Element matrix and vectors*/
    2384         ElementMatrix* Ke = element->NewElementMatrix(HOApproximationEnum);
    2385         IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
    2386         IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
    2387 
    2388         /*Retrieve all inputs and parameters*/
    2389         element->GetVerticesCoordinatesBase(&xyz_list_base);
    2390         element->FindParam(&migration_style,GroundinglineMigrationEnum);
    2391         Input* gllevelset_input = NULL;
    2392 
    2393         /*build friction object, used later on: */
    2394         Friction* friction=new Friction(element,2);
    2395 
    2396         /*Recover portion of element that is grounded*/
    2397         if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base);
    2398         if(migration_style==SubelementMigration2Enum){
    2399                 gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
    2400                 element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    2401                 //gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
    2402                 gauss=element->NewGaussBase(2);
    2403         }
    2404         else{
    2405                 gauss=element->NewGaussBase(2);
    2406         }
    2407 
    2408         /* Start  looping on the number of gaussian points: */
    2409         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2410                 gauss->GaussPoint(ig);
    2411 
    2412                 friction->GetAlpha2(&alpha2,gauss);
    2413                 if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
    2414                 if(migration_style==SubelementMigration2Enum){
    2415                         gllevelset_input->GetInputValue(&gllevelset, gauss);
    2416                         if(gllevelset<0.) alpha2=0.;
    2417                 }
    2418 
    2419                 this->GetBHOFriction(B,element,dim,xyz_list_base,gauss);
    2420                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    2421                 for(int i=0;i<dim-1;i++) D[i*(dim-1)+i]=alpha2*gauss->weight*Jdet;
    2422 
    2423                 TripleMultiply(B,dim-1,numdof,1,
    2424                                         D,dim-1,dim-1,0,
    2425                                         B,dim-1,numdof,0,
    2426                                         &Ke->values[0],1);
    2427         }
    2428 
    2429         /*Transform Coordinate System*/
    2430         if(dim==3) element->TransformStiffnessMatrixCoord(Ke,XYEnum);
    2431 
    2432         /*Clean up and return*/
    2433         delete gauss;
    2434         delete friction;
    2435         xDelete<IssmDouble>(xyz_list_base);
    2436         xDelete<IssmDouble>(B);
    2437         xDelete<IssmDouble>(D);
    2438         return Ke;
    2439 }/*}}}*/
    24402440#ifdef FSANALYTICAL
    24412441ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
    … …  
    26332633        return pe;
    26342634}/*}}}*/
    2635 void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2635void           StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    26362636        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    26372637         * For node i, Bi can be expressed in the actual coordinate system
    … …  
    26822682        xDelete<IssmDouble>(dbasis);
    26832683}/*}}}*/
    2684 void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2684void           StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2685        /*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2.
     2686         * For node i, Bi can be expressed in the actual coordinate system
     2687         * by:
     2688         *                       3D           2D
     2689         *                 Bi=[ N   0 ]    Bi=N
     2690         *                    [ 0   N ]
     2691         * where N is the finiteelement function for node i.
     2692         *
     2693         * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
     2694         */
     2695
     2696        /*Fetch number of nodes for this finite element*/
     2697        int numnodes = element->GetNumberOfNodes();
     2698
     2699        /*Get nodal functions derivatives*/
     2700        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     2701        element->NodalFunctions(basis,gauss);
     2702
     2703        /*Build L: */
     2704        if(dim==3){
     2705                for(int i=0;i<numnodes;i++){
     2706                        B[2*numnodes*0+2*i+0] = basis[i];
     2707                        B[2*numnodes*0+2*i+1] = 0.;
     2708                        B[2*numnodes*1+2*i+0] = 0.;
     2709                        B[2*numnodes*1+2*i+1] = basis[i];
     2710                }
     2711        }
     2712        else{
     2713                for(int i=0;i<numnodes;i++){
     2714                        B[i] = basis[i];
     2715                }
     2716        }
     2717
     2718        /*Clean-up*/
     2719        xDelete<IssmDouble>(basis);
     2720}/*}}}*/
     2721void           StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    26852722        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    26862723         * For node i, Bi' can be expressed in the actual coordinate system
    … …  
    27272764        xDelete<IssmDouble>(dbasis);
    27282765}/*}}}*/
    2729 void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    2730         /*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2.
    2731          * For node i, Bi can be expressed in the actual coordinate system
    2732          * by:
    2733          *                       3D           2D
    2734          *                 Bi=[ N   0 ]    Bi=N
    2735          *                    [ 0   N ]
    2736          * where N is the finiteelement function for node i.
    2737          *
    2738          * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
    2739          */
    2740 
    2741         /*Fetch number of nodes for this finite element*/
    2742         int numnodes = element->GetNumberOfNodes();
    2743 
    2744         /*Get nodal functions derivatives*/
    2745         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    2746         element->NodalFunctions(basis,gauss);
    2747 
    2748         /*Build L: */
    2749         if(dim==3){
    2750                 for(int i=0;i<numnodes;i++){
    2751                         B[2*numnodes*0+2*i+0] = basis[i];
    2752                         B[2*numnodes*0+2*i+1] = 0.;
    2753                         B[2*numnodes*1+2*i+0] = 0.;
    2754                         B[2*numnodes*1+2*i+1] = basis[i];
    2755                 }
    2756         }
    2757         else{
    2758                 for(int i=0;i<numnodes;i++){
    2759                         B[i] = basis[i];
    2760                 }
    2761         }
    2762 
    2763         /*Clean-up*/
    2764         xDelete<IssmDouble>(basis);
    2765 }/*}}}*/
    2766 void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
     2766void           StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
    27672767
    27682768        int         i,dim;
    … …  
    29822982        return Ke;
    29832983}/*}}}*/
     2984ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSShelf(Element* element){/*{{{*/
     2985
     2986        if(!element->IsFloating() || !element->IsOnBase()) return NULL;
     2987
     2988        /*If on not water or not FS, skip stiffness: */
     2989        int approximation,shelf_dampening;
     2990        element->GetInputValue(&approximation,ApproximationEnum);
     2991        if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     2992        element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
     2993        if(shelf_dampening==0) return NULL;
     2994
     2995        /*Intermediaries*/
     2996        bool        mainlyfloating;
     2997        int         j,i,dim;
     2998        IssmDouble  Jdet,slope2,scalar,dt;
     2999        IssmDouble  slope[3];
     3000        IssmDouble *xyz_list_base = NULL;
     3001        IssmDouble *xyz_list      = NULL;
     3002        Gauss*      gauss         = NULL;
     3003
     3004        /*Get problem dimension*/
     3005        element->FindParam(&dim,DomainDimensionEnum);
     3006
     3007        /*Fetch number of nodes and dof for this finite element*/
     3008        int vnumnodes = element->NumberofNodesVelocity();
     3009        int pnumnodes = element->NumberofNodesPressure();
     3010        int numdof    = vnumnodes*dim + pnumnodes;
     3011
     3012        /*Initialize Element matrix and vectors*/
     3013        ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
     3014        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3015
     3016        /*Retrieve all inputs and parameters*/
     3017        element->GetVerticesCoordinatesBase(&xyz_list_base);
     3018        element->GetVerticesCoordinates(&xyz_list);
     3019        element->FindParam(&dt,TimesteppingTimeStepEnum);
     3020        if(dt==0)   dt=1.e+5;
     3021        IssmDouble  rho_water     = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     3022        IssmDouble  gravity       = element->GetMaterialParameter(ConstantsGEnum);
     3023        Input*      surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
     3024
     3025        /* Start  looping on the number of gaussian points: */
     3026        gauss=element->NewGaussBase(3);
     3027        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3028                gauss->GaussPoint(ig);
     3029
     3030                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
     3031                element->NodalFunctionsVelocity(vbasis,gauss);
     3032                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     3033                if(dim==2) slope2=slope[0]*slope[0];
     3034                else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1];
     3035                scalar  = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt;
     3036                for(i=0;i<vnumnodes;i++){
     3037                        for(j=0;j<vnumnodes;j++){
     3038                                Ke->values[numdof*((i+1)*dim-1)+(j+1)*dim-1] += scalar*vbasis[i]*vbasis[j];
     3039                        }
     3040                }
     3041        }
     3042
     3043        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
     3044
     3045        /*Clean up and return*/
     3046        delete gauss;
     3047        xDelete<IssmDouble>(xyz_list_base);
     3048        xDelete<IssmDouble>(xyz_list);
     3049        xDelete<IssmDouble>(vbasis);
     3050        return Ke;
     3051}/*}}}*/
     3052ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscous(Element* element){/*{{{*/
     3053
     3054        /*Intermediaries*/
     3055        int         i,dim,epssize;
     3056        IssmDouble  viscosity,FSreconditioning,Jdet;
     3057        IssmDouble *xyz_list = NULL;
     3058
     3059        /*Get problem dimension*/
     3060        element->FindParam(&dim,DomainDimensionEnum);
     3061        if(dim==2) epssize = 3;
     3062        else       epssize = 6;
     3063
     3064        /*Fetch number of nodes and dof for this finite element*/
     3065        int vnumnodes = element->NumberofNodesVelocity();
     3066        int pnumnodes = element->NumberofNodesPressure();
     3067        int numdof    = vnumnodes*dim + pnumnodes;
     3068        int bsize     = epssize + 2;
     3069
     3070        /*Prepare coordinate system list*/
     3071        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     3072        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     3073        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     3074        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     3075
     3076        /*Initialize Element matrix and vectors*/
     3077        ElementMatrix* Ke     = element->NewElementMatrix(FSvelocityEnum);
     3078        IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
     3079        IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
     3080        IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
     3081
     3082        /*Retrieve all inputs and parameters*/
     3083        element->GetVerticesCoordinates(&xyz_list);
     3084        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     3085        Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
     3086        Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
     3087        Input* vz_input;
     3088        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
     3089
     3090        /* Start  looping on the number of gaussian points: */
     3091        Gauss* gauss = element->NewGauss(5);
     3092        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3093                gauss->GaussPoint(ig);
     3094
     3095                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     3096                this->GetBFS(B,element,dim,xyz_list,gauss);
     3097                this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
     3098
     3099                element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     3100                for(i=0;i<epssize;i++)     D[i*bsize+i] = + 2.*viscosity*gauss->weight*Jdet;
     3101                for(i=epssize;i<bsize;i++) D[i*bsize+i] = - FSreconditioning*gauss->weight*Jdet;
     3102
     3103                TripleMultiply(B,bsize,numdof,1,
     3104                                        D,bsize,bsize,0,
     3105                                        Bprime,bsize,numdof,0,
     3106                                        &Ke->values[0],1);
     3107        }
     3108
     3109        /*Transform Coordinate System*/
     3110        element->TransformStiffnessMatrixCoord(Ke,cs_list);
     3111
     3112        /*Clean up and return*/
     3113        delete gauss;
     3114        xDelete<IssmDouble>(xyz_list);
     3115        xDelete<IssmDouble>(D);
     3116        xDelete<IssmDouble>(Bprime);
     3117        xDelete<IssmDouble>(B);
     3118        xDelete<int>(cs_list);
     3119        return Ke;
     3120}/*}}}*/
    29843121ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousLA(Element* element){/*{{{*/
    29853122
    … …  
    31813318        return Ke;
    31823319}/*}}}*/
    3183 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscous(Element* element){/*{{{*/
    3184 
    3185         /*Intermediaries*/
    3186         int         i,dim,epssize;
    3187         IssmDouble  viscosity,FSreconditioning,Jdet;
    3188         IssmDouble *xyz_list = NULL;
    3189 
    3190         /*Get problem dimension*/
    3191         element->FindParam(&dim,DomainDimensionEnum);
    3192         if(dim==2) epssize = 3;
    3193         else       epssize = 6;
    3194 
    3195         /*Fetch number of nodes and dof for this finite element*/
    3196         int vnumnodes = element->NumberofNodesVelocity();
    3197         int pnumnodes = element->NumberofNodesPressure();
    3198         int numdof    = vnumnodes*dim + pnumnodes;
    3199         int bsize     = epssize + 2;
    3200 
    3201         /*Prepare coordinate system list*/
    3202         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    3203         if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
    3204         else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
    3205         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    3206 
    3207         /*Initialize Element matrix and vectors*/
    3208         ElementMatrix* Ke     = element->NewElementMatrix(FSvelocityEnum);
    3209         IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
    3210         IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
    3211         IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
    3212 
    3213         /*Retrieve all inputs and parameters*/
    3214         element->GetVerticesCoordinates(&xyz_list);
    3215         element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    3216         Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
    3217         Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
    3218         Input* vz_input;
    3219         if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
    3220 
    3221         /* Start  looping on the number of gaussian points: */
    3222         Gauss* gauss = element->NewGauss(5);
    3223         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3224                 gauss->GaussPoint(ig);
    3225 
    3226                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    3227                 this->GetBFS(B,element,dim,xyz_list,gauss);
    3228                 this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
    3229 
    3230                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    3231                 for(i=0;i<epssize;i++)     D[i*bsize+i] = + 2.*viscosity*gauss->weight*Jdet;
    3232                 for(i=epssize;i<bsize;i++) D[i*bsize+i] = - FSreconditioning*gauss->weight*Jdet;
    3233 
    3234                 TripleMultiply(B,bsize,numdof,1,
    3235                                         D,bsize,bsize,0,
    3236                                         Bprime,bsize,numdof,0,
    3237                                         &Ke->values[0],1);
    3238         }
    3239 
    3240         /*Transform Coordinate System*/
    3241         element->TransformStiffnessMatrixCoord(Ke,cs_list);
    3242 
    3243         /*Clean up and return*/
    3244         delete gauss;
    3245         xDelete<IssmDouble>(xyz_list);
    3246         xDelete<IssmDouble>(D);
    3247         xDelete<IssmDouble>(Bprime);
    3248         xDelete<IssmDouble>(B);
    3249         xDelete<int>(cs_list);
    3250         return Ke;
    3251 }/*}}}*/
    3252 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSShelf(Element* element){/*{{{*/
    3253 
    3254         if(!element->IsFloating() || !element->IsOnBase()) return NULL;
    3255 
    3256         /*If on not water or not FS, skip stiffness: */
    3257         int approximation,shelf_dampening;
    3258         element->GetInputValue(&approximation,ApproximationEnum);
    3259         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    3260         element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
    3261         if(shelf_dampening==0) return NULL;
    3262 
    3263         /*Intermediaries*/
    3264         bool        mainlyfloating;
    3265         int         j,i,dim;
    3266         IssmDouble  Jdet,slope2,scalar,dt;
    3267         IssmDouble  slope[3];
    3268         IssmDouble *xyz_list_base = NULL;
    3269         IssmDouble *xyz_list      = NULL;
    3270         Gauss*      gauss         = NULL;
    3271 
    3272         /*Get problem dimension*/
    3273         element->FindParam(&dim,DomainDimensionEnum);
    3274 
    3275         /*Fetch number of nodes and dof for this finite element*/
    3276         int vnumnodes = element->NumberofNodesVelocity();
    3277         int pnumnodes = element->NumberofNodesPressure();
    3278         int numdof    = vnumnodes*dim + pnumnodes;
    3279 
    3280         /*Initialize Element matrix and vectors*/
    3281         ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
    3282         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    3283 
    3284         /*Retrieve all inputs and parameters*/
    3285         element->GetVerticesCoordinatesBase(&xyz_list_base);
    3286         element->GetVerticesCoordinates(&xyz_list);
    3287         element->FindParam(&dt,TimesteppingTimeStepEnum);
    3288         if(dt==0)   dt=1.e+5;
    3289         IssmDouble  rho_water     = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    3290         IssmDouble  gravity       = element->GetMaterialParameter(ConstantsGEnum);
    3291         Input*      surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
    3292 
    3293         /* Start  looping on the number of gaussian points: */
    3294         gauss=element->NewGaussBase(3);
    3295         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3296                 gauss->GaussPoint(ig);
    3297 
    3298                 surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
    3299                 element->NodalFunctionsVelocity(vbasis,gauss);
    3300                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    3301                 if(dim==2) slope2=slope[0]*slope[0];
    3302                 else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1];
    3303                 scalar  = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt;
    3304                 for(i=0;i<vnumnodes;i++){
    3305                         for(j=0;j<vnumnodes;j++){
    3306                                 Ke->values[numdof*((i+1)*dim-1)+(j+1)*dim-1] += scalar*vbasis[i]*vbasis[j];
    3307                         }
    3308                 }
    3309         }
    3310 
    3311         /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
    3312 
    3313         /*Clean up and return*/
    3314         delete gauss;
    3315         xDelete<IssmDouble>(xyz_list_base);
    3316         xDelete<IssmDouble>(xyz_list);
    3317         xDelete<IssmDouble>(vbasis);
    3318         return Ke;
    3319 }/*}}}*/
    33203320#ifdef FSANALYTICAL
    33213321ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
    … …  
    34023402
    34033403        /*clean-up and return*/
     3404        return pe;
     3405}/*}}}*/
     3406ElementVector* StressbalanceAnalysis::CreatePVectorFSFriction(Element* element){/*{{{*/
     3407
     3408        if(!element->IsOnBase()) return NULL;
     3409
     3410        /*Intermediaries*/
     3411        int         dim;
     3412        IssmDouble  alpha2,Jdet;
     3413        IssmDouble  bed_normal[3];
     3414        IssmDouble *xyz_list_base = NULL;
     3415        Gauss*      gauss         = NULL;
     3416
     3417        /*Get problem dimension*/
     3418        element->FindParam(&dim,DomainDimensionEnum);
     3419
     3420        /*Fetch number of nodes and dof for this finite element*/
     3421        int vnumnodes = element->NumberofNodesVelocity();
     3422
     3423        /*Initialize Element matrix and vectors*/
     3424        ElementVector* pe = element->NewElementVector(FSvelocityEnum);
     3425        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3426
     3427        /*Retrieve all inputs and parameters*/
     3428        element->GetVerticesCoordinatesBase(&xyz_list_base);
     3429        Input*  alpha2_input=element->GetInput(FrictionCoefficientEnum); _assert_(alpha2_input);
     3430
     3431        /* Start  looping on the number of gaussian points: */
     3432        gauss=element->NewGaussBase(3);
     3433        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3434                gauss->GaussPoint(ig);
     3435
     3436                alpha2_input->GetInputValue(&alpha2, gauss);
     3437                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     3438                element->NodalFunctionsVelocity(vbasis,gauss);
     3439                element->NormalBase(&bed_normal[0],xyz_list_base);
     3440
     3441                for(int i=0;i<vnumnodes;i++){
     3442                        pe->values[i*dim+0] += - alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[1];
     3443                        pe->values[i*dim+1] += alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[0];
     3444                        if(dim==3){
     3445                                pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i];
     3446                        }
     3447                }
     3448
     3449        }
     3450
     3451        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
     3452
     3453        /*Clean up and return*/
     3454        delete gauss;
     3455        xDelete<IssmDouble>(xyz_list_base);
     3456        xDelete<IssmDouble>(vbasis);
     3457        return pe;
     3458}/*}}}*/
     3459ElementVector* StressbalanceAnalysis::CreatePVectorFSStress(Element* element){/*{{{*/
     3460
     3461        if(!element->IsOnBase()) return NULL;
     3462
     3463        /*Intermediaries*/
     3464        int         dim;
     3465        IssmDouble  sigmann,sigmant,Jdet,bedslope,beta;
     3466        IssmDouble *xyz_list_base = NULL;
     3467        Gauss*      gauss         = NULL;
     3468
     3469        /*Get problem dimension*/
     3470        element->FindParam(&dim,DomainDimensionEnum);
     3471
     3472        /*Fetch number of nodes and dof for this finite element*/
     3473        int vnumnodes = element->NumberofNodesVelocity();
     3474
     3475        /*Initialize Element matrix and vectors*/
     3476        ElementVector* pe = element->NewElementVector(FSvelocityEnum);
     3477        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3478
     3479        /*Retrieve all inputs and parameters*/
     3480        element->GetVerticesCoordinatesBase(&xyz_list_base);
     3481        Input*  sigmann_input=element->GetInput(VzEnum); _assert_(sigmann_input);
     3482        Input*  sigmant_input=element->GetInput(TemperatureEnum); _assert_(sigmant_input);
     3483        Input*  bedslope_input=element->GetInput(BedSlopeXEnum);     _assert_(bedslope_input);
     3484
     3485        /* Start  looping on the number of gaussian points: */
     3486        gauss=element->NewGaussBase(3);
     3487        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3488                gauss->GaussPoint(ig);
     3489
     3490                sigmann_input->GetInputValue(&sigmann, gauss);
     3491                sigmant_input->GetInputValue(&sigmant, gauss);
     3492                bedslope_input->GetInputValue(&bedslope, gauss);
     3493                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     3494                element->NodalFunctionsVelocity(vbasis,gauss);
     3495
     3496                beta=sqrt(1+bedslope*bedslope);
     3497                for(int i=0;i<vnumnodes;i++){
     3498                        pe->values[i*dim+0] += - (1./beta)*(-bedslope*sigmann + sigmant)*gauss->weight*Jdet*vbasis[i];
     3499                        pe->values[i*dim+1] += - (1./beta)*(sigmann + bedslope*sigmant)*gauss->weight*Jdet*vbasis[i];
     3500                        if(dim==3){
     3501                                //pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i];
     3502                                _error_("3d not supported yet");
     3503                        }
     3504                }
     3505
     3506        }
     3507
     3508        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
     3509
     3510        /*Clean up and return*/
     3511        delete gauss;
     3512        xDelete<IssmDouble>(xyz_list_base);
     3513        xDelete<IssmDouble>(vbasis);
    34043514        return pe;
    34053515}/*}}}*/
    … …  
    34773587                return pe3;
    34783588        }
    3479         return pe;
    3480 }/*}}}*/
    3481 ElementVector* StressbalanceAnalysis::CreatePVectorFSFriction(Element* element){/*{{{*/
    3482 
    3483         if(!element->IsOnBase()) return NULL;
    3484 
    3485         /*Intermediaries*/
    3486         int         dim;
    3487         IssmDouble  alpha2,Jdet;
    3488         IssmDouble  bed_normal[3];
    3489         IssmDouble *xyz_list_base = NULL;
    3490         Gauss*      gauss         = NULL;
    3491 
    3492         /*Get problem dimension*/
    3493         element->FindParam(&dim,DomainDimensionEnum);
    3494 
    3495         /*Fetch number of nodes and dof for this finite element*/
    3496         int vnumnodes = element->NumberofNodesVelocity();
    3497 
    3498         /*Initialize Element matrix and vectors*/
    3499         ElementVector* pe = element->NewElementVector(FSvelocityEnum);
    3500         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    3501 
    3502         /*Retrieve all inputs and parameters*/
    3503         element->GetVerticesCoordinatesBase(&xyz_list_base);
    3504         Input*  alpha2_input=element->GetInput(FrictionCoefficientEnum); _assert_(alpha2_input);
    3505 
    3506         /* Start  looping on the number of gaussian points: */
    3507         gauss=element->NewGaussBase(3);
    3508         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3509                 gauss->GaussPoint(ig);
    3510 
    3511                 alpha2_input->GetInputValue(&alpha2, gauss);
    3512                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    3513                 element->NodalFunctionsVelocity(vbasis,gauss);
    3514                 element->NormalBase(&bed_normal[0],xyz_list_base);
    3515 
    3516                 for(int i=0;i<vnumnodes;i++){
    3517                         pe->values[i*dim+0] += - alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[1];
    3518                         pe->values[i*dim+1] += alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[0];
    3519                         if(dim==3){
    3520                                 pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i];
    3521                         }
    3522                 }
    3523 
    3524         }
    3525 
    3526         /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
    3527 
    3528         /*Clean up and return*/
    3529         delete gauss;
    3530         xDelete<IssmDouble>(xyz_list_base);
    3531         xDelete<IssmDouble>(vbasis);
    3532         return pe;
    3533 }/*}}}*/
    3534 ElementVector* StressbalanceAnalysis::CreatePVectorFSStress(Element* element){/*{{{*/
    3535 
    3536         if(!element->IsOnBase()) return NULL;
    3537 
    3538         /*Intermediaries*/
    3539         int         dim;
    3540         IssmDouble  sigmann,sigmant,Jdet,bedslope,beta;
    3541         IssmDouble *xyz_list_base = NULL;
    3542         Gauss*      gauss         = NULL;
    3543 
    3544         /*Get problem dimension*/
    3545         element->FindParam(&dim,DomainDimensionEnum);
    3546 
    3547         /*Fetch number of nodes and dof for this finite element*/
    3548         int vnumnodes = element->NumberofNodesVelocity();
    3549 
    3550         /*Initialize Element matrix and vectors*/
    3551         ElementVector* pe = element->NewElementVector(FSvelocityEnum);
    3552         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    3553 
    3554         /*Retrieve all inputs and parameters*/
    3555         element->GetVerticesCoordinatesBase(&xyz_list_base);
    3556         Input*  sigmann_input=element->GetInput(VzEnum); _assert_(sigmann_input);
    3557         Input*  sigmant_input=element->GetInput(TemperatureEnum); _assert_(sigmant_input);
    3558         Input*  bedslope_input=element->GetInput(BedSlopeXEnum);     _assert_(bedslope_input);
    3559 
    3560         /* Start  looping on the number of gaussian points: */
    3561         gauss=element->NewGaussBase(3);
    3562         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3563                 gauss->GaussPoint(ig);
    3564 
    3565                 sigmann_input->GetInputValue(&sigmann, gauss);
    3566                 sigmant_input->GetInputValue(&sigmant, gauss);
    3567                 bedslope_input->GetInputValue(&bedslope, gauss);
    3568                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    3569                 element->NodalFunctionsVelocity(vbasis,gauss);
    3570 
    3571                 beta=sqrt(1+bedslope*bedslope);
    3572                 for(int i=0;i<vnumnodes;i++){
    3573                         pe->values[i*dim+0] += - (1./beta)*(-bedslope*sigmann + sigmant)*gauss->weight*Jdet*vbasis[i];
    3574                         pe->values[i*dim+1] += - (1./beta)*(sigmann + bedslope*sigmant)*gauss->weight*Jdet*vbasis[i];
    3575                         if(dim==3){
    3576                                 //pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i];
    3577                                 _error_("3d not supported yet");
    3578                         }
    3579                 }
    3580 
    3581         }
    3582 
    3583         /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
    3584 
    3585         /*Clean up and return*/
    3586         delete gauss;
    3587         xDelete<IssmDouble>(xyz_list_base);
    3588         xDelete<IssmDouble>(vbasis);
    35893589        return pe;
    35903590}/*}}}*/
    … …  
    37903790}/*}}}*/
    37913791#endif
     3792ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/
     3793
     3794        /*If no front, return NULL*/
     3795        if(!element->IsIcefront()) return NULL;
     3796
     3797        /*Intermediaries*/
     3798        int         i,dim;
     3799        IssmDouble  Jdet,pressure,surface,z;
     3800        IssmDouble      normal[3];
     3801        IssmDouble *xyz_list       = NULL;
     3802        IssmDouble *xyz_list_front = NULL;
     3803        Gauss      *gauss          = NULL;
     3804
     3805        /*Make sure current element is floating*/
     3806        if(!element->IsFloating()) return NULL;
     3807
     3808        /*Get problem dimension*/
     3809        element->FindParam(&dim,DomainDimensionEnum);
     3810
     3811        /*Fetch number of nodes and dof for this finite element*/
     3812        int vnumnodes   = element->NumberofNodesVelocity();
     3813        int pnumnodes   = element->NumberofNodesPressure();
     3814        int numvertices = element->GetNumberOfVertices();
     3815
     3816        /*Prepare coordinate system list*/
     3817        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     3818        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     3819        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     3820        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     3821
     3822        /*Initialize vectors*/
     3823        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
     3824        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3825
     3826        /*Retrieve all inputs and parameters*/
     3827        element->GetVerticesCoordinates(&xyz_list);
     3828        element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
     3829        element->NormalSection(&normal[0],xyz_list_front);
     3830        Input* surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
     3831        IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     3832        IssmDouble  gravity   = element->GetMaterialParameter(ConstantsGEnum);
     3833
     3834        /*Initialize gauss points*/
     3835        IssmDouble zmax=xyz_list[0*3+(dim-1)]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+(dim-1)]>zmax) zmax=xyz_list[i*3+(dim-1)];
     3836        IssmDouble zmin=xyz_list[0*3+(dim-1)]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+(dim-1)]<zmin) zmin=xyz_list[i*3+(dim-1)];
     3837        if(zmax>0. && zmin<0.) gauss=element->NewGauss(xyz_list,xyz_list_front,3,30);//refined in vertical because of the sea level discontinuity
     3838        else                   gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
     3839
     3840        /* Start  looping on the number of gaussian points: */
     3841        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3842                gauss->GaussPoint(ig);
     3843
     3844                element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
     3845                element->NodalFunctionsVelocity(vbasis,gauss);
     3846                surface_input->GetInputValue(&surface,gauss);
     3847                if(dim==3) z=element->GetZcoord(xyz_list,gauss);
     3848                else       z=element->GetYcoord(xyz_list,gauss);
     3849                pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
     3850
     3851                for (int i=0;i<vnumnodes;i++){
     3852                        pe->values[dim*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
     3853                        pe->values[dim*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
     3854                        if(dim==3) pe->values[dim*i+2]+= pressure*Jdet*gauss->weight*normal[2]*vbasis[i];
     3855                }
     3856        }
     3857
     3858        /*Transform coordinate system*/
     3859        element->TransformLoadVectorCoord(pe,cs_list);
     3860
     3861        /*Clean up and return*/
     3862        delete gauss;
     3863        xDelete<int>(cs_list);
     3864        xDelete<IssmDouble>(vbasis);
     3865        xDelete<IssmDouble>(xyz_list);
     3866        xDelete<IssmDouble>(xyz_list_front);
     3867        return pe;
     3868}/*}}}*/
     3869ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/
     3870
     3871        int         i,dim;
     3872        IssmDouble  Jdet,water_pressure,bed;
     3873        IssmDouble      normal[3];
     3874        IssmDouble *xyz_list_base = NULL;
     3875
     3876        /*Get basal element*/
     3877        if(!element->IsOnBase() || !element->IsFloating()) return NULL;
     3878
     3879        /*Get problem dimension*/
     3880        element->FindParam(&dim,DomainDimensionEnum);
     3881
     3882        /*Fetch number of nodes and dof for this finite element*/
     3883        int vnumnodes = element->NumberofNodesVelocity();
     3884        int pnumnodes = element->NumberofNodesPressure();
     3885
     3886        /*Prepare coordinate system list*/
     3887        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     3888        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     3889        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     3890        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     3891
     3892        /*Initialize vectors*/
     3893        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
     3894        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3895
     3896        /*Retrieve all inputs and parameters*/
     3897        element->GetVerticesCoordinatesBase(&xyz_list_base);
     3898        Input*      base_input=element->GetInput(BaseEnum); _assert_(base_input);
     3899        IssmDouble  rho_water=element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     3900        IssmDouble  gravity  =element->GetMaterialParameter(ConstantsGEnum);
     3901
     3902        /* Start  looping on the number of gaussian points: */
     3903        Gauss* gauss=element->NewGaussBase(5);
     3904        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3905                gauss->GaussPoint(ig);
     3906
     3907                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     3908                element->NodalFunctionsVelocity(vbasis,gauss);
     3909
     3910                element->NormalBase(&normal[0],xyz_list_base);
     3911                _assert_(normal[dim-1]<0.);
     3912                base_input->GetInputValue(&bed, gauss);
     3913                water_pressure=gravity*rho_water*bed;
     3914
     3915                for(i=0;i<vnumnodes;i++){
     3916                        pe->values[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0];
     3917                        pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1];
     3918                        if(dim==3){
     3919                                pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2];
     3920                        }
     3921                }
     3922        }
     3923
     3924        /*Transform coordinate system*/
     3925        element->TransformLoadVectorCoord(pe,cs_list);
     3926
     3927        /* shelf dampening*/
     3928        int shelf_dampening;
     3929        element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
     3930        if(shelf_dampening) {
     3931                Input*      mb_input=element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(mb_input);
     3932                IssmDouble dt,mb,normal_b;
     3933                element->FindParam(&dt,TimesteppingTimeStepEnum);
     3934                for(int ig=gauss->begin();ig<gauss->end();ig++){
     3935                        gauss->GaussPoint(ig);
     3936                        element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     3937                        element->NodalFunctionsVelocity(vbasis,gauss);
     3938                        element->NormalBase(&normal[0],xyz_list_base);
     3939                        if (dim==2) normal_b=normal[1];
     3940                        else if (dim==3) normal_b=sqrt(normal[0]*normal[0]+normal[1]*normal[1]);
     3941                        mb_input->GetInputValue(&mb, gauss);
     3942                        for(i=0;i<vnumnodes;i++){
     3943                                pe->values[i*dim+1] += dt*rho_water*gravity*mb*gauss->weight*Jdet*vbasis[i]*normal_b;
     3944                        }
     3945                }
     3946        }
     3947
     3948        /*Clean up and return*/
     3949        delete gauss;
     3950        xDelete<int>(cs_list);
     3951        xDelete<IssmDouble>(vbasis);
     3952        xDelete<IssmDouble>(xyz_list_base);
     3953        return pe;
     3954}/*}}}*/
    37923955ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLA(Element* element){/*{{{*/
    37933956
    … …  
    40444207        return pe;
    40454208}/*}}}*/
    4046 ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/
    4047 
    4048         int         i,dim;
    4049         IssmDouble  Jdet,water_pressure,bed;
    4050         IssmDouble      normal[3];
    4051         IssmDouble *xyz_list_base = NULL;
    4052 
    4053         /*Get basal element*/
    4054         if(!element->IsOnBase() || !element->IsFloating()) return NULL;
    4055 
    4056         /*Get problem dimension*/
    4057         element->FindParam(&dim,DomainDimensionEnum);
    4058 
    4059         /*Fetch number of nodes and dof for this finite element*/
    4060         int vnumnodes = element->NumberofNodesVelocity();
    4061         int pnumnodes = element->NumberofNodesPressure();
    4062 
    4063         /*Prepare coordinate system list*/
    4064         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    4065         if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
    4066         else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
    4067         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    4068 
    4069         /*Initialize vectors*/
    4070         ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
    4071         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    4072 
    4073         /*Retrieve all inputs and parameters*/
    4074         element->GetVerticesCoordinatesBase(&xyz_list_base);
    4075         Input*      base_input=element->GetInput(BaseEnum); _assert_(base_input);
    4076         IssmDouble  rho_water=element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    4077         IssmDouble  gravity  =element->GetMaterialParameter(ConstantsGEnum);
    4078 
    4079         /* Start  looping on the number of gaussian points: */
    4080         Gauss* gauss=element->NewGaussBase(5);
    4081         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4082                 gauss->GaussPoint(ig);
    4083 
    4084                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    4085                 element->NodalFunctionsVelocity(vbasis,gauss);
    4086 
    4087                 element->NormalBase(&normal[0],xyz_list_base);
    4088                 _assert_(normal[dim-1]<0.);
    4089                 base_input->GetInputValue(&bed, gauss);
    4090                 water_pressure=gravity*rho_water*bed;
    4091 
    4092                 for(i=0;i<vnumnodes;i++){
    4093                         pe->values[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0];
    4094                         pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1];
    4095                         if(dim==3){
    4096                                 pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2];
    4097                         }
    4098                 }
    4099         }
    4100 
    4101         /*Transform coordinate system*/
    4102         element->TransformLoadVectorCoord(pe,cs_list);
    4103 
    4104         /* shelf dampening*/
    4105         int shelf_dampening;
    4106         element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
    4107         if(shelf_dampening) {
    4108                 Input*      mb_input=element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(mb_input);
    4109                 IssmDouble dt,mb,normal_b;
    4110                 element->FindParam(&dt,TimesteppingTimeStepEnum);
    4111                 for(int ig=gauss->begin();ig<gauss->end();ig++){
    4112                         gauss->GaussPoint(ig);
    4113                         element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    4114                         element->NodalFunctionsVelocity(vbasis,gauss);
    4115                         element->NormalBase(&normal[0],xyz_list_base);
    4116                         if (dim==2) normal_b=normal[1];
    4117                         else if (dim==3) normal_b=sqrt(normal[0]*normal[0]+normal[1]*normal[1]);
    4118                         mb_input->GetInputValue(&mb, gauss);
    4119                         for(i=0;i<vnumnodes;i++){
    4120                                 pe->values[i*dim+1] += dt*rho_water*gravity*mb*gauss->weight*Jdet*vbasis[i]*normal_b;
    4121                         }
    4122                 }
    4123         }
    4124 
    4125         /*Clean up and return*/
    4126         delete gauss;
    4127         xDelete<int>(cs_list);
    4128         xDelete<IssmDouble>(vbasis);
    4129         xDelete<IssmDouble>(xyz_list_base);
    4130         return pe;
    4131 }/*}}}*/
    4132 ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/
    4133 
    4134         /*If no front, return NULL*/
    4135         if(!element->IsIcefront()) return NULL;
    4136 
    4137         /*Intermediaries*/
    4138         int         i,dim;
    4139         IssmDouble  Jdet,pressure,surface,z;
    4140         IssmDouble      normal[3];
    4141         IssmDouble *xyz_list       = NULL;
    4142         IssmDouble *xyz_list_front = NULL;
    4143         Gauss      *gauss          = NULL;
    4144 
    4145         /*Make sure current element is floating*/
    4146         if(!element->IsFloating()) return NULL;
    4147 
    4148         /*Get problem dimension*/
    4149         element->FindParam(&dim,DomainDimensionEnum);
    4150 
    4151         /*Fetch number of nodes and dof for this finite element*/
    4152         int vnumnodes   = element->NumberofNodesVelocity();
    4153         int pnumnodes   = element->NumberofNodesPressure();
    4154         int numvertices = element->GetNumberOfVertices();
    4155 
    4156         /*Prepare coordinate system list*/
    4157         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    4158         if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
    4159         else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
    4160         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    4161 
    4162         /*Initialize vectors*/
    4163         ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
    4164         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    4165 
    4166         /*Retrieve all inputs and parameters*/
    4167         element->GetVerticesCoordinates(&xyz_list);
    4168         element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
    4169         element->NormalSection(&normal[0],xyz_list_front);
    4170         Input* surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
    4171         IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    4172         IssmDouble  gravity   = element->GetMaterialParameter(ConstantsGEnum);
    4173 
    4174         /*Initialize gauss points*/
    4175         IssmDouble zmax=xyz_list[0*3+(dim-1)]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+(dim-1)]>zmax) zmax=xyz_list[i*3+(dim-1)];
    4176         IssmDouble zmin=xyz_list[0*3+(dim-1)]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+(dim-1)]<zmin) zmin=xyz_list[i*3+(dim-1)];
    4177         if(zmax>0. && zmin<0.) gauss=element->NewGauss(xyz_list,xyz_list_front,3,30);//refined in vertical because of the sea level discontinuity
    4178         else                   gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
    4179 
    4180         /* Start  looping on the number of gaussian points: */
    4181         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4182                 gauss->GaussPoint(ig);
    4183 
    4184                 element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
    4185                 element->NodalFunctionsVelocity(vbasis,gauss);
    4186                 surface_input->GetInputValue(&surface,gauss);
    4187                 if(dim==3) z=element->GetZcoord(xyz_list,gauss);
    4188                 else       z=element->GetYcoord(xyz_list,gauss);
    4189                 pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
    4190 
    4191                 for (int i=0;i<vnumnodes;i++){
    4192                         pe->values[dim*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
    4193                         pe->values[dim*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
    4194                         if(dim==3) pe->values[dim*i+2]+= pressure*Jdet*gauss->weight*normal[2]*vbasis[i];
    4195                 }
    4196         }
    4197 
    4198         /*Transform coordinate system*/
    4199         element->TransformLoadVectorCoord(pe,cs_list);
    4200 
    4201         /*Clean up and return*/
    4202         delete gauss;
    4203         xDelete<int>(cs_list);
    4204         xDelete<IssmDouble>(vbasis);
    4205         xDelete<IssmDouble>(xyz_list);
    4206         xDelete<IssmDouble>(xyz_list_front);
    4207         return pe;
    4208 }/*}}}*/
    4209 void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4209void           StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    42104210        /*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
    42114211         * For node i, Bvi can be expressed in the actual coordinate system
    … …  
    43194319        xDelete<IssmDouble>(pbasis);
    43204320}/*}}}*/
    4321 void StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4321void           StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4322        /* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     4323         * For node i, Li can be expressed in the actual coordinate system
     4324         * by in 3d
     4325         *       Li=[ h 0 0 0 ]
     4326         *                    [ 0 h 0 0 ]
     4327         *      in 2d:
     4328         *       Li=[ h 0 0 ]
     4329         * where h is the interpolation function for node i.
     4330         */
     4331
     4332        /*Fetch number of nodes for this finite element*/
     4333        int pnumnodes = element->NumberofNodesPressure();
     4334        int vnumnodes = element->NumberofNodesVelocity();
     4335        int pnumdof   = pnumnodes;
     4336        int vnumdof   = vnumnodes*dim;
     4337
     4338        /*Get nodal functions derivatives*/
     4339        IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
     4340        element->NodalFunctionsVelocity(vbasis,gauss);
     4341
     4342        /*Build B: */
     4343        if(dim==3){
     4344                for(int i=0;i<vnumnodes;i++){
     4345                        B[(vnumdof+pnumdof)*0+3*i+0] = vbasis[i];
     4346                        B[(vnumdof+pnumdof)*0+3*i+1] = 0.;
     4347                        B[(vnumdof+pnumdof)*0+3*i+2] = 0.;
     4348
     4349                        B[(vnumdof+pnumdof)*1+3*i+0] = 0.;
     4350                        B[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];
     4351                        B[(vnumdof+pnumdof)*1+3*i+2] = 0.;
     4352
     4353                        B[(vnumdof+pnumdof)*2+3*i+0] = 0.;
     4354                        B[(vnumdof+pnumdof)*2+3*i+1] = 0.;
     4355                        B[(vnumdof+pnumdof)*2+3*i+2] = vbasis[i];
     4356                }
     4357                for(int i=0;i<pnumnodes;i++){
     4358                        B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
     4359                        B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
     4360                        B[(vnumdof+pnumdof)*2+i+vnumdof+0] = 0.;
     4361                }
     4362        }
     4363        else{
     4364                for(int i=0;i<vnumnodes;i++){
     4365                        B[(vnumdof+pnumdof)*0+2*i+0] = vbasis[i];
     4366                        B[(vnumdof+pnumdof)*0+2*i+1] = 0.;
     4367
     4368                        B[(vnumdof+pnumdof)*1+2*i+0] = 0.;
     4369                        B[(vnumdof+pnumdof)*1+2*i+1] = vbasis[i];
     4370                }
     4371
     4372                for(int i=0;i<pnumnodes;i++){
     4373                        B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
     4374                        B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
     4375                }
     4376        }
     4377
     4378        /*Clean-up*/
     4379        xDelete<IssmDouble>(vbasis);
     4380}/*}}}*/
     4381void           StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    43224382        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
    43234383         *      For node i, Bi' can be expressed in the actual coordinate system
    … …  
    44324492        xDelete<IssmDouble>(pbasis);
    44334493}/*}}}*/
    4434 void StressbalanceAnalysis::GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4494void           StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4495        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2.
     4496         *      For node i, Bi' can be expressed in the actual coordinate system
     4497         *      by:
     4498         *                      Bvi' = [  dphi/dx   dphi/dy ]
     4499         *
     4500         *      In 3d
     4501         *         Bvi=[ dh/dx   dh/dy    dh/dz  ]
     4502         *      where phi is the finiteelement function for node i.
     4503         */
     4504
     4505        /*Fetch number of nodes for this finite element*/
     4506        int vnumnodes = element->NumberofNodesVelocity();
     4507
     4508        /*Get nodal functions derivatives*/
     4509        IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
     4510        element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     4511
     4512        /*Build B_prime: */
     4513        if(dim==2){
     4514                for(int i=0;i<vnumnodes;i++){
     4515                        Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
     4516                        Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
     4517                }
     4518        }
     4519        else{
     4520                for(int i=0;i<vnumnodes;i++){
     4521                        Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
     4522                        Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
     4523                        Bprime[dim*i+2] = vdbasis[2*vnumnodes+i];
     4524                }
     4525        }
     4526
     4527        /*Clean up*/
     4528        xDelete<IssmDouble>(vdbasis);
     4529}/*}}}*/
     4530void           StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4531        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
     4532         *      For node i, Bi' can be expressed in the actual coordinate system
     4533         *      by:
     4534         *                      Bvi' = [  dphi/dx     0     ]
     4535         *                                       [     0      dphi/dy ]
     4536         *                                       [  dphi/dy   dphi/dx ]
     4537         *
     4538         *      In 3d
     4539         *         Bvi=[ dh/dx     0        0    ]
     4540         *                                      [   0      dh/dy      0    ]
     4541         *                                      [   0        0      dh/dz  ]
     4542         *                                      [ dh/dy    dh/dx      0    ]
     4543         *                                      [ dh/dz      0      dh/dx  ]
     4544         *                                      [   0      dh/dz    dh/dy  ]
     4545         *      where phi is the finiteelement function for node i.
     4546         *      In 3d:
     4547         */
     4548
     4549        /*Fetch number of nodes for this finite element*/
     4550        int vnumnodes = element->NumberofNodesVelocity();
     4551
     4552        /*Get nodal functions derivatives*/
     4553        IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
     4554        element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     4555
     4556        /*Build B_prime: */
     4557        if(dim==2){
     4558                for(int i=0;i<vnumnodes;i++){
     4559                        Bprime[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
     4560                        Bprime[(dim*vnumnodes)*0+dim*i+1] = 0.;
     4561                        Bprime[(dim*vnumnodes)*1+dim*i+0] = 0.;
     4562                        Bprime[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
     4563                        Bprime[(dim*vnumnodes)*2+dim*i+0] = vdbasis[1*vnumnodes+i];
     4564                        Bprime[(dim*vnumnodes)*2+dim*i+1] = vdbasis[0*vnumnodes+i];
     4565                }
     4566        }
     4567        else{
     4568                for(int i=0;i<vnumnodes;i++){
     4569                        Bprime[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
     4570                        Bprime[(dim*vnumnodes)*0+dim*i+1] = 0.;
     4571                        Bprime[(dim*vnumnodes)*0+dim*i+2] = 0.;
     4572                        Bprime[(dim*vnumnodes)*1+dim*i+0] = 0.;
     4573                        Bprime[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
     4574                        Bprime[(dim*vnumnodes)*1+dim*i+2] = 0.;
     4575                        Bprime[(dim*vnumnodes)*2+dim*i+0] = 0.;
     4576                        Bprime[(dim*vnumnodes)*2+dim*i+1] = 0.;
     4577                        Bprime[(dim*vnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
     4578                        Bprime[(dim*vnumnodes)*3+dim*i+0] = vdbasis[1*vnumnodes+i];
     4579                        Bprime[(dim*vnumnodes)*3+dim*i+1] = vdbasis[0*vnumnodes+i];
     4580                        Bprime[(dim*vnumnodes)*3+dim*i+2] = 0.;
     4581                        Bprime[(dim*vnumnodes)*4+dim*i+0] = vdbasis[2*vnumnodes+i];
     4582                        Bprime[(dim*vnumnodes)*4+dim*i+1] = 0.;
     4583                        Bprime[(dim*vnumnodes)*4+dim*i+2] = vdbasis[0*vnumnodes+i];
     4584                        Bprime[(dim*vnumnodes)*5+dim*i+0] = 0.;
     4585                        Bprime[(dim*vnumnodes)*5+dim*i+1] = vdbasis[2*vnumnodes+i];
     4586                        Bprime[(dim*vnumnodes)*5+dim*i+2] = vdbasis[1*vnumnodes+i];
     4587                }
     4588        }
     4589
     4590        /*Clean up*/
     4591        xDelete<IssmDouble>(vdbasis);
     4592}/*}}}*/
     4593void           StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4594        /*Compute B  matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi.
     4595         */
     4596
     4597        /*Fetch number of nodes for this finite element*/
     4598        int pnumnodes;
     4599        if(dim==2) pnumnodes=3;
     4600        else pnumnodes=6;
     4601        //int pnumnodes = element->NumberofNodes(P1Enum);
     4602
     4603        /*Get nodal functions derivatives*/
     4604        IssmDouble* basis =xNew<IssmDouble>(pnumnodes);
     4605        element->NodalFunctionsP1(basis,gauss);
     4606
     4607        /*Build B: */
     4608        for(int i=0;i<pnumnodes;i++){
     4609                B[i] = basis[i];
     4610        }
     4611
     4612        /*Clean up*/
     4613        xDelete<IssmDouble>(basis);
     4614}/*}}}*/
     4615void           StressbalanceAnalysis::GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    44354616        /*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
    44364617         * For node i, Bvi can be expressed in the actual coordinate system
    … …  
    44964677        xDelete<IssmDouble>(vdbasis);
    44974678}/*}}}*/
    4498 void StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    4499         /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
    4500          *      For node i, Bi' can be expressed in the actual coordinate system
    4501          *      by:
    4502          *                      Bvi' = [  dphi/dx     0     ]
    4503          *                                       [     0      dphi/dy ]
    4504          *                                       [  dphi/dy   dphi/dx ]
    4505          *
    4506          *      In 3d
    4507          *         Bvi=[ dh/dx     0        0    ]
    4508          *                                      [   0      dh/dy      0    ]
    4509          *                                      [   0        0      dh/dz  ]
    4510          *                                      [ dh/dy    dh/dx      0    ]
    4511          *                                      [ dh/dz      0      dh/dx  ]
    4512          *                                      [   0      dh/dz    dh/dy  ]
    4513          *      where phi is the finiteelement function for node i.
    4514          *      In 3d:
    4515          */
    4516 
    4517         /*Fetch number of nodes for this finite element*/
    4518         int vnumnodes = element->NumberofNodesVelocity();
    4519 
    4520         /*Get nodal functions derivatives*/
    4521         IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
    4522         element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
    4523 
    4524         /*Build B_prime: */
    4525         if(dim==2){
    4526                 for(int i=0;i<vnumnodes;i++){
    4527                         Bprime[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
    4528                         Bprime[(dim*vnumnodes)*0+dim*i+1] = 0.;
    4529                         Bprime[(dim*vnumnodes)*1+dim*i+0] = 0.;
    4530                         Bprime[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
    4531                         Bprime[(dim*vnumnodes)*2+dim*i+0] = vdbasis[1*vnumnodes+i];
    4532                         Bprime[(dim*vnumnodes)*2+dim*i+1] = vdbasis[0*vnumnodes+i];
    4533                 }
    4534         }
    4535         else{
    4536                 for(int i=0;i<vnumnodes;i++){
    4537                         Bprime[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
    4538                         Bprime[(dim*vnumnodes)*0+dim*i+1] = 0.;
    4539                         Bprime[(dim*vnumnodes)*0+dim*i+2] = 0.;
    4540                         Bprime[(dim*vnumnodes)*1+dim*i+0] = 0.;
    4541                         Bprime[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
    4542                         Bprime[(dim*vnumnodes)*1+dim*i+2] = 0.;
    4543                         Bprime[(dim*vnumnodes)*2+dim*i+0] = 0.;
    4544                         Bprime[(dim*vnumnodes)*2+dim*i+1] = 0.;
    4545                         Bprime[(dim*vnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
    4546                         Bprime[(dim*vnumnodes)*3+dim*i+0] = vdbasis[1*vnumnodes+i];
    4547                         Bprime[(dim*vnumnodes)*3+dim*i+1] = vdbasis[0*vnumnodes+i];
    4548                         Bprime[(dim*vnumnodes)*3+dim*i+2] = 0.;
    4549                         Bprime[(dim*vnumnodes)*4+dim*i+0] = vdbasis[2*vnumnodes+i];
    4550                         Bprime[(dim*vnumnodes)*4+dim*i+1] = 0.;
    4551                         Bprime[(dim*vnumnodes)*4+dim*i+2] = vdbasis[0*vnumnodes+i];
    4552                         Bprime[(dim*vnumnodes)*5+dim*i+0] = 0.;
    4553                         Bprime[(dim*vnumnodes)*5+dim*i+1] = vdbasis[2*vnumnodes+i];
    4554                         Bprime[(dim*vnumnodes)*5+dim*i+2] = vdbasis[1*vnumnodes+i];
    4555                 }
    4556         }
    4557 
    4558         /*Clean up*/
    4559         xDelete<IssmDouble>(vdbasis);
    4560 }/*}}}*/
    4561 void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    4562         /*Compute B  matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi.
    4563          */
    4564 
    4565         /*Fetch number of nodes for this finite element*/
    4566         int pnumnodes;
    4567         if(dim==2) pnumnodes=3;
    4568         else pnumnodes=6;
    4569         //int pnumnodes = element->NumberofNodes(P1Enum);
    4570 
    4571         /*Get nodal functions derivatives*/
    4572         IssmDouble* basis =xNew<IssmDouble>(pnumnodes);
    4573         element->NodalFunctionsP1(basis,gauss);
    4574 
    4575         /*Build B: */
    4576         for(int i=0;i<pnumnodes;i++){
    4577                 B[i] = basis[i];
    4578         }
    4579 
    4580         /*Clean up*/
    4581         xDelete<IssmDouble>(basis);
    4582 }/*}}}*/
    4583 void StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    4584         /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2.
    4585          *      For node i, Bi' can be expressed in the actual coordinate system
    4586          *      by:
    4587          *                      Bvi' = [  dphi/dx   dphi/dy ]
    4588          *
    4589          *      In 3d
    4590          *         Bvi=[ dh/dx   dh/dy    dh/dz  ]
    4591          *      where phi is the finiteelement function for node i.
    4592          */
    4593 
    4594         /*Fetch number of nodes for this finite element*/
    4595         int vnumnodes = element->NumberofNodesVelocity();
    4596 
    4597         /*Get nodal functions derivatives*/
    4598         IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
    4599         element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
    4600 
    4601         /*Build B_prime: */
    4602         if(dim==2){
    4603                 for(int i=0;i<vnumnodes;i++){
    4604                         Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
    4605                         Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
    4606                 }
    4607         }
    4608         else{
    4609                 for(int i=0;i<vnumnodes;i++){
    4610                         Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
    4611                         Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
    4612                         Bprime[dim*i+2] = vdbasis[2*vnumnodes+i];
    4613                 }
    4614         }
    4615 
    4616         /*Clean up*/
    4617         xDelete<IssmDouble>(vdbasis);
    4618 }/*}}}*/
    4619 void StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    4620         /* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    4621          * For node i, Li can be expressed in the actual coordinate system
    4622          * by in 3d
    4623          *       Li=[ h 0 0 0 ]
    4624          *                    [ 0 h 0 0 ]
    4625          *      in 2d:
    4626          *       Li=[ h 0 0 ]
    4627          * where h is the interpolation function for node i.
    4628          */
    4629 
    4630         /*Fetch number of nodes for this finite element*/
    4631         int pnumnodes = element->NumberofNodesPressure();
    4632         int vnumnodes = element->NumberofNodesVelocity();
    4633         int pnumdof   = pnumnodes;
    4634         int vnumdof   = vnumnodes*dim;
    4635 
    4636         /*Get nodal functions derivatives*/
    4637         IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
    4638         element->NodalFunctionsVelocity(vbasis,gauss);
    4639 
    4640         /*Build B: */
    4641         if(dim==3){
    4642                 for(int i=0;i<vnumnodes;i++){
    4643                         B[(vnumdof+pnumdof)*0+3*i+0] = vbasis[i];
    4644                         B[(vnumdof+pnumdof)*0+3*i+1] = 0.;
    4645                         B[(vnumdof+pnumdof)*0+3*i+2] = 0.;
    4646 
    4647                         B[(vnumdof+pnumdof)*1+3*i+0] = 0.;
    4648                         B[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];
    4649                         B[(vnumdof+pnumdof)*1+3*i+2] = 0.;
    4650 
    4651                         B[(vnumdof+pnumdof)*2+3*i+0] = 0.;
    4652                         B[(vnumdof+pnumdof)*2+3*i+1] = 0.;
    4653                         B[(vnumdof+pnumdof)*2+3*i+2] = vbasis[i];
    4654                 }
    4655                 for(int i=0;i<pnumnodes;i++){
    4656                         B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
    4657                         B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
    4658                         B[(vnumdof+pnumdof)*2+i+vnumdof+0] = 0.;
    4659                 }
    4660         }
    4661         else{
    4662                 for(int i=0;i<vnumnodes;i++){
    4663                         B[(vnumdof+pnumdof)*0+2*i+0] = vbasis[i];
    4664                         B[(vnumdof+pnumdof)*0+2*i+1] = 0.;
    4665 
    4666                         B[(vnumdof+pnumdof)*1+2*i+0] = 0.;
    4667                         B[(vnumdof+pnumdof)*1+2*i+1] = vbasis[i];
    4668                 }
    4669 
    4670                 for(int i=0;i<pnumnodes;i++){
    4671                         B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
    4672                         B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
    4673                 }
    4674         }
    4675 
    4676         /*Clean-up*/
    4677         xDelete<IssmDouble>(vbasis);
    4678 }/*}}}*/
    4679 void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4679void           StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    46804680        /*Compute C  matrix. C=[Cp1 Cp2 ...] where:
    46814681         *     Cpi=[phi phi].
    … …  
    46984698        xDelete<IssmDouble>(basis);
    46994699}/*}}}*/
    4700 void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4700void           StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    47014701        /*      Compute C'  matrix. C'=[C1' C2' ...]
    47024702         *                      Ci' = [  phi  0  ]
    … …  
    47464746        xDelete<IssmDouble>(vbasis);
    47474747}/*}}}*/
    4748 void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
     4748void           StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    47494749
    47504750        int*         vdoflist=NULL;
    … …  
    48094809        xDelete<IssmDouble>(vvalues);
    48104810}/*}}}*/
    4811 void StressbalanceAnalysis::InitializeXTH(Elements* elements,Parameters* parameters){/*{{{*/
     4811void           StressbalanceAnalysis::InitializeXTH(Elements* elements,Parameters* parameters){/*{{{*/
    48124812
    48134813        /*Intermediaries*/
    … …  
    48884888
    48894889}/*}}}*/
    4890 void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
     4890void           StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
    48914891
    48924892        bool         results_on_nodes;
    … …  
    49784978        xDelete<int>(cs_list);
    49794979}/*}}}*/
    4980 void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters){/*{{{*/
     4980void           StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters){/*{{{*/
    49814981
    49824982        /*Intermediaries*/
    … …  
    51915191        }
    51925192}/*}}}*/
    5193 void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters){/*{{{*/
     5193void           StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters){/*{{{*/
    51945194
    51955195        /*Intermediaries*/
    … …  
    53275327
    53285328/*Coupling (Tiling)*/
    5329 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3d(Element* element){/*{{{*/
    5330 
    5331         /*compute all stiffness matrices for this element*/
    5332         ElementMatrix* Ke1=CreateKMatrixSSA3dViscous(element);
    5333         ElementMatrix* Ke2=CreateKMatrixSSA3dFriction(element);
    5334         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    5335 
    5336         /*clean-up and return*/
    5337         delete Ke1;
    5338         delete Ke2;
    5339         return Ke;
    5340 }/*}}}*/
    5341 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dFriction(Element* element){/*{{{*/
    5342 
    5343         /*Initialize Element matrix and return if necessary*/
    5344         if(element->IsFloating() || !element->IsOnBase()) return NULL;
    5345 
    5346         /*Build a tria element using the 3 nodes of the base of the penta. Then use
    5347          * the tria functionality to build a friction stiffness matrix on these 3
    5348          * nodes: */
    5349         Element* basalelement = element->SpawnBasalElement();
    5350         ElementMatrix* Ke=CreateKMatrixSSAFriction(basalelement);
    5351         basalelement->DeleteMaterials(); delete basalelement;
    5352 
    5353         /*clean-up and return*/
    5354         return Ke;
    5355 }/*}}}*/
    5356 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dViscous(Element* element){/*{{{*/
    5357 
    5358         /*Constants*/
    5359         const int    numdof2d=2*3;
    5360 
    5361         /*Intermediaries */
    5362         int         i,j,approximation;
    5363         int         dim=3;
    5364         IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot;
    5365         IssmDouble  epsilon[5],oldepsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
    5366         IssmDouble  epsilons[6];                    //6 for FS
    5367         IssmDouble  B[3][numdof2d];
    5368         IssmDouble  Bprime[3][numdof2d];
    5369         IssmDouble  D[3][3]= {0.0};                 // material matrix, simple scalar matrix.
    5370         IssmDouble  D_scalar;
    5371         IssmDouble  Ke_gg[numdof2d][numdof2d]={0.0};
    5372         IssmDouble  *xyz_list  = NULL;
    5373 
    5374         /*Find penta on bed as this is a SSA elements: */
    5375         Element* pentabase=element->GetBasalElement();
    5376         Element* basaltria=pentabase->SpawnBasalElement();
    5377 
    5378         /*Initialize Element matrix*/
    5379         ElementMatrix* Ke=basaltria->NewElementMatrix(SSAApproximationEnum);
    5380         element->GetInputValue(&approximation,ApproximationEnum);
    5381 
    5382         /*Retrieve all inputs and parameters*/
    5383         element->GetVerticesCoordinates(&xyz_list);
    5384         element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
    5385         Input* vx_input   =element->GetInput(VxEnum);       _assert_(vx_input);
    5386         Input* vy_input   =element->GetInput(VyEnum);       _assert_(vy_input);
    5387         Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);
    5388         Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
    5389         Input* vz_input   =element->GetInput(VzEnum);       _assert_(vz_input);
    5390 
    5391         /* Start  looping on the number of gaussian points: */
    5392         Gauss* gauss=element->NewGauss(5);
    5393         Gauss* gauss_tria=new GaussTria();
    5394         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5395 
    5396                 gauss->GaussPoint(ig);
    5397                 gauss->SynchronizeGaussBase(gauss_tria);
    5398 
    5399                 element->JacobianDeterminant(&Jdet, xyz_list,gauss);
    5400                 this->GetBSSA(&B[0][0],basaltria,2,xyz_list, gauss_tria);
    5401                 this->GetBSSAprime(&Bprime[0][0],basaltria,2,xyz_list, gauss_tria);
    5402 
    5403                 if(approximation==SSAHOApproximationEnum){
    5404                         element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
    5405                         element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
    5406                         newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    5407                 }
    5408                 else if (approximation==SSAFSApproximationEnum){
    5409                         element->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    5410                 }
    5411                 else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet");
    5412 
    5413                 D_scalar=2*newviscosity*gauss->weight*Jdet;
    5414                 for (i=0;i<3;i++) D[i][i]=D_scalar;
    5415 
    5416                 TripleMultiply( &B[0][0],3,numdof2d,1,
    5417                                         &D[0][0],3,3,0,
    5418                                         &Bprime[0][0],3,numdof2d,0,
    5419                                         &Ke_gg[0][0],1);
    5420 
    5421         }
    5422         for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg[i][j];
    5423 
    5424         /*Transform Coordinate System*/
    5425         basaltria->TransformStiffnessMatrixCoord(Ke,XYEnum);
    5426 
    5427         /*Clean up and return*/
    5428         xDelete<IssmDouble>(xyz_list);
    5429         delete basaltria->material;
    5430         delete basaltria;
    5431         delete gauss_tria;
    5432         delete gauss;
    5433         return Ke;
    5434 
    5435 }/*}}}*/
    5436 ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFS(Element* element){/*{{{*/
    5437 
    5438         /*compute all stiffness matrices for this element*/
    5439         ElementMatrix* Ke1=CreateKMatrixFS(element);
    5440         int indices[3]={18,19,20};
    5441         Ke1->StaticCondensation(3,&indices[0]);
    5442         int init = element->FiniteElement();
    5443         element->SetTemporaryElementType(P1Enum); // P1 needed for HO
    5444         ElementMatrix* Ke2=CreateKMatrixHO(element);
    5445         element->SetTemporaryElementType(init); // P1 needed for HO
    5446         ElementMatrix* Ke3=CreateKMatrixCouplingHOFS(element);
    5447         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    5448 
    5449         /*clean-up and return*/
    5450         delete Ke1;
    5451         delete Ke2;
    5452         delete Ke3;
    5453         return Ke;
    5454 }/*}}}*/
    5455 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAHO(Element* element){/*{{{*/
    5456 
    5457         /*compute all stiffness matrices for this element*/
    5458         ElementMatrix* Ke1=CreateKMatrixSSA3d(element);
    5459         ElementMatrix* Ke2=CreateKMatrixHO(element);
    5460         ElementMatrix* Ke3=CreateKMatrixCouplingSSAHO(element);
    5461         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    5462 
    5463         /*clean-up and return*/
    5464         delete Ke1;
    5465         delete Ke2;
    5466         delete Ke3;
    5467         return Ke;
    5468 }/*}}}*/
    5469 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFS(Element* element){/*{{{*/
    5470 
    5471         /*compute all stiffness matrices for this element*/
    5472         ElementMatrix* Ke1=CreateKMatrixFS(element);
    5473         int indices[3]={18,19,20};
    5474         Ke1->StaticCondensation(3,&indices[0]);
    5475         int init = element->FiniteElement();
    5476         element->SetTemporaryElementType(P1Enum);
    5477         ElementMatrix* Ke2=CreateKMatrixSSA3d(element);
    5478         element->SetTemporaryElementType(init);
    5479         ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element);
    5480         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    5481 
    5482         /*clean-up and return*/
    5483         delete Ke1;
    5484         delete Ke2;
    5485         delete Ke3;
    5486         return Ke;
    5487 }/*}}}*/
    54885329ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingHOFS(Element* element){/*{{{*/
    54895330
    … …  
    55515392        delete Ke2;
    55525393        return Ke;
    5553 }/*}}}*/
    5554 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHO(Element* element){/*{{{*/
    5555 
    5556         /*compute all stiffness matrices for this element*/
    5557         ElementMatrix* Ke1=CreateKMatrixCouplingSSAHOViscous(element);
    5558         ElementMatrix* Ke2=CreateKMatrixCouplingSSAHOFriction(element);
    5559         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    5560 
    5561         /*clean-up and return*/
    5562         delete Ke1;
    5563         delete Ke2;
    5564         return Ke;
    5565 }/*}}}*/
    5566 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOFriction(Element* element){/*{{{*/
    5567 
    5568         if(element->IsFloating() || !element->IsOnBase()) return NULL;
    5569 
    5570         /*Constants*/
    5571         int numnodes    = element->GetNumberOfNodes();
    5572         int numdof      = 2*numnodes;
    5573         int numdoftotal = 4*numnodes;
    5574 
    5575         /*Intermediaries */
    5576         int         i,j;
    5577         IssmDouble  Jdet2d,alpha2;
    5578         IssmDouble *xyz_list_tria = NULL;
    5579         IssmDouble* L             = xNewZeroInit<IssmDouble>(2*numdof);
    5580         IssmDouble  DL[2][2]      = {{ 0,0 },{0,0}}; //for basal drag
    5581         IssmDouble  DL_scalar;
    5582         IssmDouble* Ke_gg         = xNewZeroInit<IssmDouble>(numdof*numdof);
    5583         Node      **node_list     = xNew<Node*>(2*numnodes);
    5584         int*        cs_list       = xNew<int>(2*numnodes);
    5585 
    5586         /*Initialize Element matrix and return if necessary*/
    5587         ElementMatrix* Ke1=element->NewElementMatrix(SSAApproximationEnum);
    5588         ElementMatrix* Ke2=element->NewElementMatrix(HOApproximationEnum);
    5589         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    5590         delete Ke1; delete Ke2;
    5591 
    5592         /*Prepare node list*/
    5593         for(i=0;i<numnodes;i++){
    5594                 node_list[i+0*numnodes] = element->GetNode(i);
    5595                 node_list[i+1*numnodes] = element->GetNode(i);
    5596                 cs_list[i+0*numnodes] = XYEnum;
    5597                 cs_list[i+1*numnodes] = XYEnum;
    5598         }
    5599 
    5600         /*retrieve inputs :*/
    5601         element->GetVerticesCoordinatesBase(&xyz_list_tria);
    5602 
    5603         /*build friction object, used later on: */
    5604         Friction* friction=new Friction(element,2);
    5605 
    5606         /* Start  looping on the number of gaussian points: */
    5607         Gauss* gauss=element->NewGaussBase(2);
    5608         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5609 
    5610                 gauss->GaussPoint(ig);
    5611 
    5612                 /*Friction: */
    5613                 friction->GetAlpha2(&alpha2,gauss);
    5614                 element->JacobianDeterminantBase(&Jdet2d, xyz_list_tria,gauss);
    5615                 this->GetBHOFriction(L,element,3,xyz_list_tria,gauss);
    5616 
    5617                 DL_scalar=alpha2*gauss->weight*Jdet2d;
    5618                 for (i=0;i<2;i++) DL[i][i]=DL_scalar;
    5619 
    5620                 /*  Do the triple producte tL*D*L: */
    5621                 TripleMultiply( L,2,numdof,1,
    5622                                         &DL[0][0],2,2,0,
    5623                                         L,2,numdof,0,
    5624                                         Ke_gg,1);
    5625         }
    5626 
    5627         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i*numdof+j];
    5628         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i*numdof+j];
    5629 
    5630         /*Transform Coordinate System*/
    5631         element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
    5632 
    5633         /*Clean up and return*/
    5634         delete gauss;
    5635         delete friction;
    5636         xDelete<int>(cs_list);
    5637         xDelete<Node*>(node_list);
    5638         xDelete<IssmDouble>(xyz_list_tria);
    5639         xDelete<IssmDouble>(Ke_gg);
    5640         xDelete<IssmDouble>(L);
    5641         return Ke;
    5642 }/*}}}*/
    5643 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOViscous(Element* element){/*{{{*/
    5644 
    5645         /*Constants*/
    5646         int numnodes    = element->GetNumberOfNodes();
    5647         int numdofm     = 1 *numnodes; //*2/2
    5648         int numdofp     = 2 *numnodes;
    5649         int numdoftotal = 2 *2 *numnodes;//2 dof per nodes and 2 sets of nodes for HO and SSA
    5650 
    5651         /*Intermediaries */
    5652         int         i,j;
    5653         IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
    5654         IssmDouble  *xyz_list      = NULL;
    5655         IssmDouble* B              = xNew<IssmDouble>(3*numdofp);
    5656         IssmDouble* Bprime         = xNew<IssmDouble>(3*numdofm);
    5657         IssmDouble  D[3][3]={0.0}; // material matrix, simple scalar matrix.
    5658         IssmDouble  D_scalar;
    5659         IssmDouble* Ke_gg          = xNewZeroInit<IssmDouble>(numdofp*numdofm);
    5660         Node       **node_list     = xNew<Node*>(2*numnodes);
    5661         int*         cs_list= xNew<int>(2*numnodes);
    5662 
    5663         /*Find penta on bed as HO must be coupled to the dofs on the bed: */
    5664         Element* pentabase=element->GetBasalElement();
    5665         Element* basaltria=pentabase->SpawnBasalElement();
    5666 
    5667         /*prepare node list*/
    5668         for(i=0;i<numnodes;i++){
    5669                 node_list[i+0*numnodes] = pentabase->GetNode(i);
    5670                 node_list[i+1*numnodes] = element  ->GetNode(i);
    5671                 cs_list[i+0*numnodes] = XYEnum;
    5672                 cs_list[i+1*numnodes] = XYEnum;
    5673         }
    5674 
    5675         /*Initialize Element matrix*/
    5676         ElementMatrix* Ke1= pentabase->NewElementMatrix(SSAApproximationEnum);
    5677         ElementMatrix* Ke2= element  ->NewElementMatrix(HOApproximationEnum);
    5678         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    5679         delete Ke1; delete Ke2;
    5680 
    5681         /* Get node coordinates and dof list: */
    5682         element->GetVerticesCoordinates(&xyz_list);
    5683         element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
    5684         Input* vx_input   =element->GetInput(VxEnum);       _assert_(vx_input);
    5685         Input* vy_input   =element->GetInput(VyEnum);       _assert_(vy_input);
    5686         Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);
    5687         Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
    5688 
    5689         /* Start  looping on the number of gaussian points: */
    5690         Gauss* gauss=element->NewGauss(5);
    5691         Gauss* gauss_tria=new GaussTria();
    5692         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5693 
    5694                 gauss->GaussPoint(ig);
    5695                 gauss->SynchronizeGaussBase(gauss_tria);
    5696 
    5697                 element->JacobianDeterminant(&Jdet, xyz_list,gauss);
    5698                 this->GetBSSAHO(B, element,xyz_list, gauss);
    5699                 this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria);
    5700                 element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input);
    5701                 element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input);
    5702 
    5703                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    5704                 D_scalar=2*newviscosity*gauss->weight*Jdet;
    5705                 for (i=0;i<3;i++) D[i][i]=D_scalar;
    5706 
    5707                 TripleMultiply( B,3,numdofp,1,
    5708                                         &D[0][0],3,3,0,
    5709                                         Bprime,3,numdofm,0,
    5710                                         Ke_gg,1);
    5711         }
    5712         for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j];
    5713         for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i];
    5714 
    5715         /*Transform Coordinate System*/
    5716         element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
    5717 
    5718         /*Clean-up and return*/
    5719         basaltria->DeleteMaterials(); delete basaltria;
    5720        
    5721         delete gauss;
    5722         delete gauss_tria;
    5723         xDelete<IssmDouble>(B);
    5724         xDelete<IssmDouble>(Bprime);
    5725         xDelete<IssmDouble>(Ke_gg);
    5726         xDelete<IssmDouble>(xyz_list);
    5727         xDelete<Node*>(node_list);
    5728         xDelete<int>(cs_list);
    5729         return Ke;
    5730 
    57315394}/*}}}*/
    57325395ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFS(Element* element){/*{{{*/
    … …  
    59705633
    59715634}/*}}}*/
    5972 ElementVector* StressbalanceAnalysis::CreatePVectorSSAHO(Element* element){/*{{{*/
     5635ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHO(Element* element){/*{{{*/
     5636
     5637        /*compute all stiffness matrices for this element*/
     5638        ElementMatrix* Ke1=CreateKMatrixCouplingSSAHOViscous(element);
     5639        ElementMatrix* Ke2=CreateKMatrixCouplingSSAHOFriction(element);
     5640        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
     5641
     5642        /*clean-up and return*/
     5643        delete Ke1;
     5644        delete Ke2;
     5645        return Ke;
     5646}/*}}}*/
     5647ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOFriction(Element* element){/*{{{*/
     5648
     5649        if(element->IsFloating() || !element->IsOnBase()) return NULL;
     5650
     5651        /*Constants*/
     5652        int numnodes    = element->GetNumberOfNodes();
     5653        int numdof      = 2*numnodes;
     5654        int numdoftotal = 4*numnodes;
     5655
     5656        /*Intermediaries */
     5657        int         i,j;
     5658        IssmDouble  Jdet2d,alpha2;
     5659        IssmDouble *xyz_list_tria = NULL;
     5660        IssmDouble* L             = xNewZeroInit<IssmDouble>(2*numdof);
     5661        IssmDouble  DL[2][2]      = {{ 0,0 },{0,0}}; //for basal drag
     5662        IssmDouble  DL_scalar;
     5663        IssmDouble* Ke_gg         = xNewZeroInit<IssmDouble>(numdof*numdof);
     5664        Node      **node_list     = xNew<Node*>(2*numnodes);
     5665        int*        cs_list       = xNew<int>(2*numnodes);
     5666
     5667        /*Initialize Element matrix and return if necessary*/
     5668        ElementMatrix* Ke1=element->NewElementMatrix(SSAApproximationEnum);
     5669        ElementMatrix* Ke2=element->NewElementMatrix(HOApproximationEnum);
     5670        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
     5671        delete Ke1; delete Ke2;
     5672
     5673        /*Prepare node list*/
     5674        for(i=0;i<numnodes;i++){
     5675                node_list[i+0*numnodes] = element->GetNode(i);
     5676                node_list[i+1*numnodes] = element->GetNode(i);
     5677                cs_list[i+0*numnodes] = XYEnum;
     5678                cs_list[i+1*numnodes] = XYEnum;
     5679        }
     5680
     5681        /*retrieve inputs :*/
     5682        element->GetVerticesCoordinatesBase(&xyz_list_tria);
     5683
     5684        /*build friction object, used later on: */
     5685        Friction* friction=new Friction(element,2);
     5686
     5687        /* Start  looping on the number of gaussian points: */
     5688        Gauss* gauss=element->NewGaussBase(2);
     5689        for(int ig=gauss->begin();ig<gauss->end();ig++){
     5690
     5691                gauss->GaussPoint(ig);
     5692
     5693                /*Friction: */
     5694                friction->GetAlpha2(&alpha2,gauss);
     5695                element->JacobianDeterminantBase(&Jdet2d, xyz_list_tria,gauss);
     5696                this->GetBHOFriction(L,element,3,xyz_list_tria,gauss);
     5697
     5698                DL_scalar=alpha2*gauss->weight*Jdet2d;
     5699                for (i=0;i<2;i++) DL[i][i]=DL_scalar;
     5700
     5701                /*  Do the triple producte tL*D*L: */
     5702                TripleMultiply( L,2,numdof,1,
     5703                                        &DL[0][0],2,2,0,
     5704                                        L,2,numdof,0,
     5705                                        Ke_gg,1);
     5706        }
     5707
     5708        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i*numdof+j];
     5709        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i*numdof+j];
     5710
     5711        /*Transform Coordinate System*/
     5712        element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
     5713
     5714        /*Clean up and return*/
     5715        delete gauss;
     5716        delete friction;
     5717        xDelete<int>(cs_list);
     5718        xDelete<Node*>(node_list);
     5719        xDelete<IssmDouble>(xyz_list_tria);
     5720        xDelete<IssmDouble>(Ke_gg);
     5721        xDelete<IssmDouble>(L);
     5722        return Ke;
     5723}/*}}}*/
     5724ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOViscous(Element* element){/*{{{*/
     5725
     5726        /*Constants*/
     5727        int numnodes    = element->GetNumberOfNodes();
     5728        int numdofm     = 1 *numnodes; //*2/2
     5729        int numdofp     = 2 *numnodes;
     5730        int numdoftotal = 2 *2 *numnodes;//2 dof per nodes and 2 sets of nodes for HO and SSA
     5731
     5732        /*Intermediaries */
     5733        int         i,j;
     5734        IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
     5735        IssmDouble  *xyz_list      = NULL;
     5736        IssmDouble* B              = xNew<IssmDouble>(3*numdofp);
     5737        IssmDouble* Bprime         = xNew<IssmDouble>(3*numdofm);
     5738        IssmDouble  D[3][3]={0.0}; // material matrix, simple scalar matrix.
     5739        IssmDouble  D_scalar;
     5740        IssmDouble* Ke_gg          = xNewZeroInit<IssmDouble>(numdofp*numdofm);
     5741        Node       **node_list     = xNew<Node*>(2*numnodes);
     5742        int*         cs_list= xNew<int>(2*numnodes);
     5743
     5744        /*Find penta on bed as HO must be coupled to the dofs on the bed: */
     5745        Element* pentabase=element->GetBasalElement();
     5746        Element* basaltria=pentabase->SpawnBasalElement();
     5747
     5748        /*prepare node list*/
     5749        for(i=0;i<numnodes;i++){
     5750                node_list[i+0*numnodes] = pentabase->GetNode(i);
     5751                node_list[i+1*numnodes] = element  ->GetNode(i);
     5752                cs_list[i+0*numnodes] = XYEnum;
     5753                cs_list[i+1*numnodes] = XYEnum;
     5754        }
     5755
     5756        /*Initialize Element matrix*/
     5757        ElementMatrix* Ke1= pentabase->NewElementMatrix(SSAApproximationEnum);
     5758        ElementMatrix* Ke2= element  ->NewElementMatrix(HOApproximationEnum);
     5759        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     5760        delete Ke1; delete Ke2;
     5761
     5762        /* Get node coordinates and dof list: */
     5763        element->GetVerticesCoordinates(&xyz_list);
     5764        element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
     5765        Input* vx_input   =element->GetInput(VxEnum);       _assert_(vx_input);
     5766        Input* vy_input   =element->GetInput(VyEnum);       _assert_(vy_input);
     5767        Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);
     5768        Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
     5769
     5770        /* Start  looping on the number of gaussian points: */
     5771        Gauss* gauss=element->NewGauss(5);
     5772        Gauss* gauss_tria=new GaussTria();
     5773        for(int ig=gauss->begin();ig<gauss->end();ig++){
     5774
     5775                gauss->GaussPoint(ig);
     5776                gauss->SynchronizeGaussBase(gauss_tria);
     5777
     5778                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
     5779                this->GetBSSAHO(B, element,xyz_list, gauss);
     5780                this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria);
     5781                element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input);
     5782                element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input);
     5783
     5784                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
     5785                D_scalar=2*newviscosity*gauss->weight*Jdet;
     5786                for (i=0;i<3;i++) D[i][i]=D_scalar;
     5787
     5788                TripleMultiply( B,3,numdofp,1,
     5789                                        &D[0][0],3,3,0,
     5790                                        Bprime,3,numdofm,0,
     5791                                        Ke_gg,1);
     5792        }
     5793        for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j];
     5794        for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i];
     5795
     5796        /*Transform Coordinate System*/
     5797        element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
     5798
     5799        /*Clean-up and return*/
     5800        basaltria->DeleteMaterials(); delete basaltria;
     5801       
     5802        delete gauss;
     5803        delete gauss_tria;
     5804        xDelete<IssmDouble>(B);
     5805        xDelete<IssmDouble>(Bprime);
     5806        xDelete<IssmDouble>(Ke_gg);
     5807        xDelete<IssmDouble>(xyz_list);
     5808        xDelete<Node*>(node_list);
     5809        xDelete<int>(cs_list);
     5810        return Ke;
     5811
     5812}/*}}}*/
     5813ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFS(Element* element){/*{{{*/
     5814
     5815        /*compute all stiffness matrices for this element*/
     5816        ElementMatrix* Ke1=CreateKMatrixFS(element);
     5817        int indices[3]={18,19,20};
     5818        Ke1->StaticCondensation(3,&indices[0]);
     5819        int init = element->FiniteElement();
     5820        element->SetTemporaryElementType(P1Enum); // P1 needed for HO
     5821        ElementMatrix* Ke2=CreateKMatrixHO(element);
     5822        element->SetTemporaryElementType(init); // P1 needed for HO
     5823        ElementMatrix* Ke3=CreateKMatrixCouplingHOFS(element);
     5824        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
     5825
     5826        /*clean-up and return*/
     5827        delete Ke1;
     5828        delete Ke2;
     5829        delete Ke3;
     5830        return Ke;
     5831}/*}}}*/
     5832ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFS(Element* element){/*{{{*/
     5833
     5834        /*compute all stiffness matrices for this element*/
     5835        ElementMatrix* Ke1=CreateKMatrixFS(element);
     5836        int indices[3]={18,19,20};
     5837        Ke1->StaticCondensation(3,&indices[0]);
     5838        int init = element->FiniteElement();
     5839        element->SetTemporaryElementType(P1Enum);
     5840        ElementMatrix* Ke2=CreateKMatrixSSA3d(element);
     5841        element->SetTemporaryElementType(init);
     5842        ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element);
     5843        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
     5844
     5845        /*clean-up and return*/
     5846        delete Ke1;
     5847        delete Ke2;
     5848        delete Ke3;
     5849        return Ke;
     5850}/*}}}*/
     5851ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAHO(Element* element){/*{{{*/
     5852
     5853        /*compute all stiffness matrices for this element*/
     5854        ElementMatrix* Ke1=CreateKMatrixSSA3d(element);
     5855        ElementMatrix* Ke2=CreateKMatrixHO(element);
     5856        ElementMatrix* Ke3=CreateKMatrixCouplingSSAHO(element);
     5857        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
     5858
     5859        /*clean-up and return*/
     5860        delete Ke1;
     5861        delete Ke2;
     5862        delete Ke3;
     5863        return Ke;
     5864}/*}}}*/
     5865ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3d(Element* element){/*{{{*/
     5866
     5867        /*compute all stiffness matrices for this element*/
     5868        ElementMatrix* Ke1=CreateKMatrixSSA3dViscous(element);
     5869        ElementMatrix* Ke2=CreateKMatrixSSA3dFriction(element);
     5870        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     5871
     5872        /*clean-up and return*/
     5873        delete Ke1;
     5874        delete Ke2;
     5875        return Ke;
     5876}/*}}}*/
     5877ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dFriction(Element* element){/*{{{*/
     5878
     5879        /*Initialize Element matrix and return if necessary*/
     5880        if(element->IsFloating() || !element->IsOnBase()) return NULL;
     5881
     5882        /*Build a tria element using the 3 nodes of the base of the penta. Then use
     5883         * the tria functionality to build a friction stiffness matrix on these 3
     5884         * nodes: */
     5885        Element* basalelement = element->SpawnBasalElement();
     5886        ElementMatrix* Ke=CreateKMatrixSSAFriction(basalelement);
     5887        basalelement->DeleteMaterials(); delete basalelement;
     5888
     5889        /*clean-up and return*/
     5890        return Ke;
     5891}/*}}}*/
     5892ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dViscous(Element* element){/*{{{*/
     5893
     5894        /*Constants*/
     5895        const int    numdof2d=2*3;
     5896
     5897        /*Intermediaries */
     5898        int         i,j,approximation;
     5899        int         dim=3;
     5900        IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot;
     5901        IssmDouble  epsilon[5],oldepsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
     5902        IssmDouble  epsilons[6];                    //6 for FS
     5903        IssmDouble  B[3][numdof2d];
     5904        IssmDouble  Bprime[3][numdof2d];
     5905        IssmDouble  D[3][3]= {0.0};                 // material matrix, simple scalar matrix.
     5906        IssmDouble  D_scalar;
     5907        IssmDouble  Ke_gg[numdof2d][numdof2d]={0.0};
     5908        IssmDouble  *xyz_list  = NULL;
     5909
     5910        /*Find penta on bed as this is a SSA elements: */
     5911        Element* pentabase=element->GetBasalElement();
     5912        Element* basaltria=pentabase->SpawnBasalElement();
     5913
     5914        /*Initialize Element matrix*/
     5915        ElementMatrix* Ke=basaltria->NewElementMatrix(SSAApproximationEnum);
     5916        element->GetInputValue(&approximation,ApproximationEnum);
     5917
     5918        /*Retrieve all inputs and parameters*/
     5919        element->GetVerticesCoordinates(&xyz_list);
     5920        element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
     5921        Input* vx_input   =element->GetInput(VxEnum);       _assert_(vx_input);
     5922        Input* vy_input   =element->GetInput(VyEnum);       _assert_(vy_input);
     5923        Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);
     5924        Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
     5925        Input* vz_input   =element->GetInput(VzEnum);       _assert_(vz_input);
     5926
     5927        /* Start  looping on the number of gaussian points: */
     5928        Gauss* gauss=element->NewGauss(5);
     5929        Gauss* gauss_tria=new GaussTria();
     5930        for(int ig=gauss->begin();ig<gauss->end();ig++){
     5931
     5932                gauss->GaussPoint(ig);
     5933                gauss->SynchronizeGaussBase(gauss_tria);
     5934
     5935                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
     5936                this->GetBSSA(&B[0][0],basaltria,2,xyz_list, gauss_tria);
     5937                this->GetBSSAprime(&Bprime[0][0],basaltria,2,xyz_list, gauss_tria);
     5938
     5939                if(approximation==SSAHOApproximationEnum){
     5940                        element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
     5941                        element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
     5942                        newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
     5943                }
     5944                else if (approximation==SSAFSApproximationEnum){
     5945                        element->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     5946                }
     5947                else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet");
     5948
     5949                D_scalar=2*newviscosity*gauss->weight*Jdet;
     5950                for (i=0;i<3;i++) D[i][i]=D_scalar;
     5951
     5952                TripleMultiply( &B[0][0],3,numdof2d,1,
     5953                                        &D[0][0],3,3,0,
     5954                                        &Bprime[0][0],3,numdof2d,0,
     5955                                        &Ke_gg[0][0],1);
     5956
     5957        }
     5958        for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg[i][j];
     5959
     5960        /*Transform Coordinate System*/
     5961        basaltria->TransformStiffnessMatrixCoord(Ke,XYEnum);
     5962
     5963        /*Clean up and return*/
     5964        xDelete<IssmDouble>(xyz_list);
     5965        delete basaltria->material;
     5966        delete basaltria;
     5967        delete gauss_tria;
     5968        delete gauss;
     5969        return Ke;
     5970
     5971}/*}}}*/
     5972ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFS(Element* element){/*{{{*/
    59735973
    59745974        /*compute all load vectors for this element*/
    5975         ElementVector* pe1=CreatePVectorSSA(element);
    5976         ElementVector* pe2=CreatePVectorHO(element);
     5975        ElementVector* pe1=CreatePVectorCouplingHOFSViscous(element);
     5976        ElementVector* pe2=CreatePVectorCouplingHOFSFriction(element);
    59775977        ElementVector* pe =new ElementVector(pe1,pe2);
    59785978
    … …  
    59805980        delete pe1;
    59815981        delete pe2;
     5982        return pe;
     5983}/*}}}*/
     5984ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSFriction(Element* element){/*{{{*/
     5985
     5986        /*Intermediaries*/
     5987        int         i,approximation;
     5988        int         dim=3;
     5989        IssmDouble  Jdet,Jdet2d,FSreconditioning;
     5990        IssmDouble      bed_normal[3];
     5991        IssmDouble  viscosity, w, alpha2_gauss;
     5992        IssmDouble  dw[3];
     5993        IssmDouble      *xyz_list_tria = NULL;
     5994        IssmDouble  *xyz_list      = NULL;
     5995        IssmDouble  basis[6]; //for the six nodes of the penta
     5996
     5997        /*Initialize Element vector and return if necessary*/
     5998        if(!element->IsOnBase() || element->IsFloating()) return NULL;
     5999        element->GetInputValue(&approximation,ApproximationEnum);
     6000        if(approximation!=HOFSApproximationEnum) return NULL;
     6001
     6002        int vnumnodes = element->NumberofNodesVelocity();
     6003        int pnumnodes = element->NumberofNodesPressure();
     6004        int numnodes  = vnumnodes+pnumnodes;
     6005
     6006        /*Prepare coordinate system list*/
     6007        int*   cs_list   = xNew<int>(vnumnodes+pnumnodes);
     6008        Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
     6009        for(i=0;i<vnumnodes;i++){
     6010                cs_list[i]           = XYZEnum;
     6011                node_list[i]           = element->GetNode(i);
     6012        }
     6013        for(i=0;i<pnumnodes;i++){
     6014                cs_list[vnumnodes+i] = PressureEnum;
     6015                node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
     6016        }
     6017
     6018        ElementVector* pe=element->NewElementVector(FSvelocityEnum);
     6019
     6020        /*Retrieve all inputs and parameters*/
     6021        element->GetVerticesCoordinates(&xyz_list);
     6022        element->GetVerticesCoordinatesBase(&xyz_list_tria);
     6023        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     6024        Input* vx_input=  element->GetInput(VxEnum);   _assert_(vx_input);
     6025        Input* vy_input=  element->GetInput(VyEnum);   _assert_(vy_input);
     6026        Input* vz_input=  element->GetInput(VzEnum);   _assert_(vz_input);
     6027        Input* vzHO_input=element->GetInput(VzHOEnum); _assert_(vzHO_input);
     6028
     6029        /*build friction object, used later on: */
     6030        Friction* friction=new Friction(element,3);
     6031
     6032        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     6033        Gauss* gauss=element->NewGaussBase(2);
     6034        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6035
     6036                gauss->GaussPoint(ig);
     6037
     6038                element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss);
     6039                element->NodalFunctionsP1(basis, gauss);
     6040
     6041                vzHO_input->GetInputValue(&w, gauss);
     6042                vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
     6043
     6044                element->NormalBase(&bed_normal[0],xyz_list_tria);
     6045                element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     6046                friction->GetAlpha2(&alpha2_gauss,gauss);
     6047
     6048                for(i=0;i<3;i++){
     6049                        pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
     6050                        pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
     6051                        pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
     6052                }
     6053        }
     6054
     6055        /*Transform coordinate system*/
     6056        element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
     6057
     6058        /*Clean up and return*/
     6059        xDelete<int>(cs_list);
     6060        xDelete<Node*>(node_list);
     6061        xDelete<IssmDouble>(xyz_list);
     6062        xDelete<IssmDouble>(xyz_list_tria);
     6063        delete gauss;
     6064        delete friction;
     6065        return pe;
     6066}/*}}}*/
     6067ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSViscous(Element* element){/*{{{*/
     6068
     6069        /*Intermediaries */
     6070        int         i,approximation;
     6071        int         dim=3;
     6072        IssmDouble  viscosity,Jdet,FSreconditioning;
     6073        IssmDouble  dw[3];
     6074        IssmDouble  *xyz_list = NULL;
     6075        IssmDouble  basis[6]; //for the six nodes of the penta
     6076        IssmDouble  dbasis[3][6]; //for the six nodes of the penta
     6077
     6078        /*Initialize Element vector and return if necessary*/
     6079        element->GetInputValue(&approximation,ApproximationEnum);
     6080        if(approximation!=HOFSApproximationEnum) return NULL;
     6081        int   vnumnodes = element->NumberofNodesVelocity();
     6082        int   pnumnodes = element->NumberofNodesPressure();
     6083
     6084        /*Prepare coordinate system list*/
     6085        int*   cs_list   = xNew<int>(vnumnodes+pnumnodes);
     6086        Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
     6087        for(i=0;i<vnumnodes;i++){
     6088                cs_list[i]             = XYZEnum;
     6089                node_list[i]           = element->GetNode(i);
     6090        }
     6091        for(i=0;i<pnumnodes;i++){
     6092                cs_list[vnumnodes+i]   = PressureEnum;
     6093                node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
     6094        }
     6095        ElementVector* pe = element->NewElementVector(FSvelocityEnum);
     6096
     6097        /*Retrieve all inputs and parameters*/
     6098        element->GetVerticesCoordinates(&xyz_list);
     6099        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     6100        Input* vx_input   =element->GetInput(VxEnum);   _assert_(vx_input);
     6101        Input* vy_input   =element->GetInput(VyEnum);   _assert_(vy_input);
     6102        Input* vz_input   =element->GetInput(VzEnum);   _assert_(vz_input);
     6103        Input* vzHO_input=element->GetInput(VzHOEnum);  _assert_(vzHO_input);
     6104
     6105        /* Start  looping on the number of gaussian points: */
     6106        Gauss* gauss=element->NewGauss(5);
     6107        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6108
     6109                gauss->GaussPoint(ig);
     6110
     6111                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
     6112                element->NodalFunctionsP1(&basis[0],gauss);
     6113                element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
     6114               
     6115                element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     6116                vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
     6117
     6118                for(i=0;i<6;i++){
     6119                        pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
     6120                        pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
     6121                        pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
     6122                        pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
     6123                }
     6124        }
     6125
     6126        /*Transform coordinate system*/
     6127        element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
     6128
     6129        /*Clean up and return*/
     6130        xDelete<int>(cs_list);
     6131        xDelete<Node*>(node_list);
     6132        xDelete<IssmDouble>(xyz_list);
     6133        delete gauss;
     6134        return pe;
     6135}/*}}}*/
     6136ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFS(Element* element){/*{{{*/
     6137
     6138        /*compute all load vectors for this element*/
     6139        ElementVector* pe1=CreatePVectorCouplingSSAFSViscous(element);
     6140        ElementVector* pe2=CreatePVectorCouplingSSAFSFriction(element);
     6141        ElementVector* pe =new ElementVector(pe1,pe2);
     6142
     6143        /*clean-up and return*/
     6144        delete pe1;
     6145        delete pe2;
     6146        return pe;
     6147}/*}}}*/
     6148ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSFriction(Element* element){/*{{{*/
     6149
     6150        /*Intermediaries*/
     6151        int         i,j,approximation;
     6152        int         dim=3;
     6153        IssmDouble  Jdet,Jdet2d,FSreconditioning;
     6154        IssmDouble      bed_normal[3];
     6155        IssmDouble  viscosity, w, alpha2_gauss;
     6156        IssmDouble  dw[3];
     6157        IssmDouble  basis[6]; //for the six nodes of the penta
     6158        IssmDouble      *xyz_list_tria = NULL;
     6159        IssmDouble  *xyz_list      = NULL;
     6160
     6161        /*Initialize Element vector and return if necessary*/
     6162        if(!element->IsOnBase() || element->IsFloating()) return NULL;
     6163        element->GetInputValue(&approximation,ApproximationEnum);
     6164        if(approximation!=SSAFSApproximationEnum) return NULL;
     6165        int vnumnodes = element->NumberofNodesVelocity();
     6166        int pnumnodes = element->NumberofNodesPressure();
     6167
     6168        /*Prepare coordinate system list*/
     6169        int* cs_list     = xNew<int>(vnumnodes+pnumnodes);
     6170        Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
     6171        for(i=0;i<vnumnodes;i++){
     6172                cs_list[i]             = XYZEnum;
     6173                node_list[i]           = element->GetNode(i);
     6174        }
     6175        for(i=0;i<pnumnodes;i++){
     6176                cs_list[vnumnodes+i]   = PressureEnum;
     6177                node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
     6178        }
     6179        ElementVector* pe=element->NewElementVector(FSvelocityEnum);
     6180
     6181        /*Retrieve all inputs and parameters*/
     6182        element->GetVerticesCoordinates(&xyz_list);
     6183        element->GetVerticesCoordinatesBase(&xyz_list_tria);
     6184        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     6185        Input* vx_input=   element->GetInput(VxEnum);    _assert_(vx_input);
     6186        Input* vy_input=   element->GetInput(VyEnum);    _assert_(vy_input);
     6187        Input* vz_input=   element->GetInput(VzEnum);    _assert_(vz_input);
     6188        Input* vzSSA_input=element->GetInput(VzSSAEnum); _assert_(vzSSA_input);
     6189
     6190        /*build friction object, used later on: */
     6191        Friction* friction=new Friction(element,3);
     6192
     6193        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     6194        Gauss* gauss=element->NewGaussBase(2);
     6195        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6196
     6197                gauss->GaussPoint(ig);
     6198
     6199                element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss);
     6200                element->NodalFunctionsP1(basis, gauss);
     6201
     6202                vzSSA_input->GetInputValue(&w, gauss);
     6203                vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
     6204
     6205                element->NormalBase(&bed_normal[0],xyz_list_tria);
     6206                element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     6207                friction->GetAlpha2(&alpha2_gauss,gauss);
     6208
     6209                for(i=0;i<3;i++){
     6210                        pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
     6211                        pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
     6212                        pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
     6213                }
     6214        }
     6215
     6216        /*Transform coordinate system*/
     6217        element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
     6218
     6219        /*Clean up and return*/
     6220        xDelete<int>(cs_list);
     6221        xDelete<IssmDouble>(xyz_list);
     6222        xDelete<IssmDouble>(xyz_list_tria);
     6223        xDelete<Node*>(node_list);
     6224        delete gauss;
     6225        delete friction;
     6226        return pe;
     6227}/*}}}*/
     6228ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSViscous(Element* element){/*{{{*/
     6229
     6230        /*Intermediaries */
     6231        int         i,approximation;
     6232        IssmDouble  viscosity,Jdet,FSreconditioning;
     6233        IssmDouble  dw[3];
     6234        IssmDouble  *xyz_list = NULL;
     6235        IssmDouble  basis[6]; //for the six nodes of the penta
     6236        IssmDouble  dbasis[3][6]; //for the six nodes of the penta
     6237
     6238        /*Initialize Element vector and return if necessary*/
     6239        element->GetInputValue(&approximation,ApproximationEnum);
     6240        if(approximation!=SSAFSApproximationEnum) return NULL;
     6241        int vnumnodes = element->NumberofNodesVelocity();
     6242        int pnumnodes = element->NumberofNodesPressure();
     6243
     6244        /*Prepare coordinate system list*/
     6245        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     6246        Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
     6247        for(i=0;i<vnumnodes;i++){
     6248                cs_list[i]             = XYZEnum;
     6249                node_list[i]           = element->GetNode(i);
     6250        }
     6251        for(i=0;i<pnumnodes;i++){
     6252                cs_list[vnumnodes+i]   = PressureEnum;
     6253                node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
     6254        }
     6255        ElementVector* pe=element->NewElementVector(FSvelocityEnum);
     6256
     6257        /*Retrieve all inputs and parameters*/
     6258        element->GetVerticesCoordinates(&xyz_list);
     6259        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     6260        Input* vx_input   =element->GetInput(VxEnum);      _assert_(vx_input);
     6261        Input* vy_input   =element->GetInput(VyEnum);      _assert_(vy_input);
     6262        Input* vz_input   =element->GetInput(VzEnum);      _assert_(vz_input);
     6263        Input* vzSSA_input=element->GetInput(VzSSAEnum);   _assert_(vzSSA_input);
     6264
     6265        /* Start  looping on the number of gaussian points: */
     6266        Gauss* gauss=element->NewGauss(5);
     6267        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6268
     6269                gauss->GaussPoint(ig);
     6270                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     6271                element->NodalFunctionsP1(&basis[0], gauss);
     6272                element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
     6273
     6274                vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
     6275                element->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
     6276
     6277                for(i=0;i<6;i++){
     6278                        pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
     6279                        pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
     6280                        pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
     6281                        pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
     6282                }
     6283        }
     6284
     6285        /*Transform coordinate system*/
     6286        element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
     6287
     6288        /*Clean up and return*/
     6289        xDelete<int>(cs_list);
     6290        xDelete<Node*>(node_list);
     6291        xDelete<IssmDouble>(xyz_list);
     6292        delete gauss;
     6293        return pe;
     6294}/*}}}*/
     6295ElementVector* StressbalanceAnalysis::CreatePVectorHOFS(Element* element){/*{{{*/
     6296
     6297        /*compute all load vectors for this element*/
     6298        int init = element->FiniteElement();
     6299        element->SetTemporaryElementType(P1Enum);
     6300        ElementVector* pe1=CreatePVectorHO(element);
     6301        element->SetTemporaryElementType(init);
     6302        ElementVector* pe2=CreatePVectorFS(element);
     6303        int indices[3]={18,19,20};
     6304        element->SetTemporaryElementType(MINIcondensedEnum);
     6305        ElementMatrix* Ke = CreateKMatrixFS(element);
     6306        element->SetTemporaryElementType(init);
     6307        pe2->StaticCondensation(Ke,3,&indices[0]);
     6308        delete Ke;
     6309        ElementVector* pe3=CreatePVectorCouplingHOFS(element);
     6310        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     6311
     6312        /*clean-up and return*/
     6313        delete pe1;
     6314        delete pe2;
     6315        delete pe3;
    59826316        return pe;
    59836317}/*}}}*/
    … …  
    60056339        return pe;
    60066340}/*}}}*/
    6007 ElementVector* StressbalanceAnalysis::CreatePVectorHOFS(Element* element){/*{{{*/
     6341ElementVector* StressbalanceAnalysis::CreatePVectorSSAHO(Element* element){/*{{{*/
    60086342
    60096343        /*compute all load vectors for this element*/
    6010         int init = element->FiniteElement();
    6011         element->SetTemporaryElementType(P1Enum);
    6012         ElementVector* pe1=CreatePVectorHO(element);
    6013         element->SetTemporaryElementType(init);
    6014         ElementVector* pe2=CreatePVectorFS(element);
    6015         int indices[3]={18,19,20};
    6016         element->SetTemporaryElementType(MINIcondensedEnum);
    6017         ElementMatrix* Ke = CreateKMatrixFS(element);
    6018         element->SetTemporaryElementType(init);
    6019         pe2->StaticCondensation(Ke,3,&indices[0]);
    6020         delete Ke;
    6021         ElementVector* pe3=CreatePVectorCouplingHOFS(element);
    6022         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    6023 
    6024         /*clean-up and return*/
    6025         delete pe1;
    6026         delete pe2;
    6027         delete pe3;
    6028         return pe;
    6029 }/*}}}*/
    6030 ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFS(Element* element){/*{{{*/
    6031 
    6032         /*compute all load vectors for this element*/
    6033         ElementVector* pe1=CreatePVectorCouplingHOFSViscous(element);
    6034         ElementVector* pe2=CreatePVectorCouplingHOFSFriction(element);
     6344        ElementVector* pe1=CreatePVectorSSA(element);
     6345        ElementVector* pe2=CreatePVectorHO(element);
    60356346        ElementVector* pe =new ElementVector(pe1,pe2);
    60366347
    … …  
    60406351        return pe;
    60416352}/*}}}*/
    6042 ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSFriction(Element* element){/*{{{*/
    6043 
    6044         /*Intermediaries*/
    6045         int         i,approximation;
    6046         int         dim=3;
    6047         IssmDouble  Jdet,Jdet2d,FSreconditioning;
    6048         IssmDouble      bed_normal[3];
    6049         IssmDouble  viscosity, w, alpha2_gauss;
    6050         IssmDouble  dw[3];
    6051         IssmDouble      *xyz_list_tria = NULL;
    6052         IssmDouble  *xyz_list      = NULL;
    6053         IssmDouble  basis[6]; //for the six nodes of the penta
    6054 
    6055         /*Initialize Element vector and return if necessary*/
    6056         if(!element->IsOnBase() || element->IsFloating()) return NULL;
    6057         element->GetInputValue(&approximation,ApproximationEnum);
    6058         if(approximation!=HOFSApproximationEnum) return NULL;
    6059 
    6060         int vnumnodes = element->NumberofNodesVelocity();
    6061         int pnumnodes = element->NumberofNodesPressure();
    6062         int numnodes  = vnumnodes+pnumnodes;
    6063 
    6064         /*Prepare coordinate system list*/
    6065         int*   cs_list   = xNew<int>(vnumnodes+pnumnodes);
    6066         Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
    6067         for(i=0;i<vnumnodes;i++){
    6068                 cs_list[i]           = XYZEnum;
    6069                 node_list[i]           = element->GetNode(i);
    6070         }
    6071         for(i=0;i<pnumnodes;i++){
    6072                 cs_list[vnumnodes+i] = PressureEnum;
    6073                 node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
    6074         }
    6075 
    6076         ElementVector* pe=element->NewElementVector(FSvelocityEnum);
    6077 
    6078         /*Retrieve all inputs and parameters*/
    6079         element->GetVerticesCoordinates(&xyz_list);
    6080         element->GetVerticesCoordinatesBase(&xyz_list_tria);
    6081         element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    6082         Input* vx_input=  element->GetInput(VxEnum);   _assert_(vx_input);
    6083         Input* vy_input=  element->GetInput(VyEnum);   _assert_(vy_input);
    6084         Input* vz_input=  element->GetInput(VzEnum);   _assert_(vz_input);
    6085         Input* vzHO_input=element->GetInput(VzHOEnum); _assert_(vzHO_input);
    6086 
    6087         /*build friction object, used later on: */
    6088         Friction* friction=new Friction(element,3);
    6089 
    6090         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    6091         Gauss* gauss=element->NewGaussBase(2);
    6092         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6093 
    6094                 gauss->GaussPoint(ig);
    6095 
    6096                 element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss);
    6097                 element->NodalFunctionsP1(basis, gauss);
    6098 
    6099                 vzHO_input->GetInputValue(&w, gauss);
    6100                 vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
    6101 
    6102                 element->NormalBase(&bed_normal[0],xyz_list_tria);
    6103                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    6104                 friction->GetAlpha2(&alpha2_gauss,gauss);
    6105 
    6106                 for(i=0;i<3;i++){
    6107                         pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
    6108                         pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
    6109                         pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
    6110                 }
    6111         }
    6112 
    6113         /*Transform coordinate system*/
    6114         element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
    6115 
    6116         /*Clean up and return*/
    6117         xDelete<int>(cs_list);
    6118         xDelete<Node*>(node_list);
    6119         xDelete<IssmDouble>(xyz_list);
    6120         xDelete<IssmDouble>(xyz_list_tria);
    6121         delete gauss;
    6122         delete friction;
    6123         return pe;
    6124 }/*}}}*/
    6125 ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSViscous(Element* element){/*{{{*/
    6126 
    6127         /*Intermediaries */
    6128         int         i,approximation;
    6129         int         dim=3;
    6130         IssmDouble  viscosity,Jdet,FSreconditioning;
    6131         IssmDouble  dw[3];
    6132         IssmDouble  *xyz_list = NULL;
    6133         IssmDouble  basis[6]; //for the six nodes of the penta
    6134         IssmDouble  dbasis[3][6]; //for the six nodes of the penta
    6135 
    6136         /*Initialize Element vector and return if necessary*/
    6137         element->GetInputValue(&approximation,ApproximationEnum);
    6138         if(approximation!=HOFSApproximationEnum) return NULL;
    6139         int   vnumnodes = element->NumberofNodesVelocity();
    6140         int   pnumnodes = element->NumberofNodesPressure();
    6141 
    6142         /*Prepare coordinate system list*/
    6143         int*   cs_list   = xNew<int>(vnumnodes+pnumnodes);
    6144         Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
    6145         for(i=0;i<vnumnodes;i++){
    6146                 cs_list[i]             = XYZEnum;
    6147                 node_list[i]           = element->GetNode(i);
    6148         }
    6149         for(i=0;i<pnumnodes;i++){
    6150                 cs_list[vnumnodes+i]   = PressureEnum;
    6151                 node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
    6152         }
    6153         ElementVector* pe = element->NewElementVector(FSvelocityEnum);
    6154 
    6155         /*Retrieve all inputs and parameters*/
    6156         element->GetVerticesCoordinates(&xyz_list);
    6157         element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    6158         Input* vx_input   =element->GetInput(VxEnum);   _assert_(vx_input);
    6159         Input* vy_input   =element->GetInput(VyEnum);   _assert_(vy_input);
    6160         Input* vz_input   =element->GetInput(VzEnum);   _assert_(vz_input);
    6161         Input* vzHO_input=element->GetInput(VzHOEnum);  _assert_(vzHO_input);
    6162 
    6163         /* Start  looping on the number of gaussian points: */
    6164         Gauss* gauss=element->NewGauss(5);
    6165         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6166 
    6167                 gauss->GaussPoint(ig);
    6168 
    6169                 element->JacobianDeterminant(&Jdet, xyz_list,gauss);
    6170                 element->NodalFunctionsP1(&basis[0],gauss);
    6171                 element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    6172                
    6173                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    6174                 vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
    6175 
    6176                 for(i=0;i<6;i++){
    6177                         pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
    6178                         pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
    6179                         pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
    6180                         pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
    6181                 }
    6182         }
    6183 
    6184         /*Transform coordinate system*/
    6185         element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
    6186 
    6187         /*Clean up and return*/
    6188         xDelete<int>(cs_list);
    6189         xDelete<Node*>(node_list);
    6190         xDelete<IssmDouble>(xyz_list);
    6191         delete gauss;
    6192         return pe;
    6193 }/*}}}*/
    6194 ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFS(Element* element){/*{{{*/
    6195 
    6196         /*compute all load vectors for this element*/
    6197         ElementVector* pe1=CreatePVectorCouplingSSAFSViscous(element);
    6198         ElementVector* pe2=CreatePVectorCouplingSSAFSFriction(element);
    6199         ElementVector* pe =new ElementVector(pe1,pe2);
    6200 
    6201         /*clean-up and return*/
    6202         delete pe1;
    6203         delete pe2;
    6204         return pe;
    6205 }/*}}}*/
    6206 ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSFriction(Element* element){/*{{{*/
    6207 
    6208         /*Intermediaries*/
    6209         int         i,j,approximation;
    6210         int         dim=3;
    6211         IssmDouble  Jdet,Jdet2d,FSreconditioning;
    6212         IssmDouble      bed_normal[3];
    6213         IssmDouble  viscosity, w, alpha2_gauss;
    6214         IssmDouble  dw[3];
    6215         IssmDouble  basis[6]; //for the six nodes of the penta
    6216         IssmDouble      *xyz_list_tria = NULL;
    6217         IssmDouble  *xyz_list      = NULL;
    6218 
    6219         /*Initialize Element vector and return if necessary*/
    6220         if(!element->IsOnBase() || element->IsFloating()) return NULL;
    6221         element->GetInputValue(&approximation,ApproximationEnum);
    6222         if(approximation!=SSAFSApproximationEnum) return NULL;
    6223         int vnumnodes = element->NumberofNodesVelocity();
    6224         int pnumnodes = element->NumberofNodesPressure();
    6225 
    6226         /*Prepare coordinate system list*/
    6227         int* cs_list     = xNew<int>(vnumnodes+pnumnodes);
    6228         Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
    6229         for(i=0;i<vnumnodes;i++){
    6230                 cs_list[i]             = XYZEnum;
    6231                 node_list[i]           = element->GetNode(i);
    6232         }
    6233         for(i=0;i<pnumnodes;i++){
    6234                 cs_list[vnumnodes+i]   = PressureEnum;
    6235                 node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
    6236         }
    6237         ElementVector* pe=element->NewElementVector(FSvelocityEnum);
    6238 
    6239         /*Retrieve all inputs and parameters*/
    6240         element->GetVerticesCoordinates(&xyz_list);
    6241         element->GetVerticesCoordinatesBase(&xyz_list_tria);
    6242         element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    6243         Input* vx_input=   element->GetInput(VxEnum);    _assert_(vx_input);
    6244         Input* vy_input=   element->GetInput(VyEnum);    _assert_(vy_input);
    6245         Input* vz_input=   element->GetInput(VzEnum);    _assert_(vz_input);
    6246         Input* vzSSA_input=element->GetInput(VzSSAEnum); _assert_(vzSSA_input);
    6247 
    6248         /*build friction object, used later on: */
    6249         Friction* friction=new Friction(element,3);
    6250 
    6251         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    6252         Gauss* gauss=element->NewGaussBase(2);
    6253         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6254 
    6255                 gauss->GaussPoint(ig);
    6256 
    6257                 element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss);
    6258                 element->NodalFunctionsP1(basis, gauss);
    6259 
    6260                 vzSSA_input->GetInputValue(&w, gauss);
    6261                 vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
    6262 
    6263                 element->NormalBase(&bed_normal[0],xyz_list_tria);
    6264                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    6265                 friction->GetAlpha2(&alpha2_gauss,gauss);
    6266 
    6267                 for(i=0;i<3;i++){
    6268                         pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
    6269                         pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
    6270                         pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
    6271                 }
    6272         }
    6273 
    6274         /*Transform coordinate system*/
    6275         element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
    6276 
    6277         /*Clean up and return*/
    6278         xDelete<int>(cs_list);
    6279         xDelete<IssmDouble>(xyz_list);
    6280         xDelete<IssmDouble>(xyz_list_tria);
    6281         xDelete<Node*>(node_list);
    6282         delete gauss;
    6283         delete friction;
    6284         return pe;
    6285 }/*}}}*/
    6286 ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSViscous(Element* element){/*{{{*/
    6287 
    6288         /*Intermediaries */
    6289         int         i,approximation;
    6290         IssmDouble  viscosity,Jdet,FSreconditioning;
    6291         IssmDouble  dw[3];
    6292         IssmDouble  *xyz_list = NULL;
    6293         IssmDouble  basis[6]; //for the six nodes of the penta
    6294         IssmDouble  dbasis[3][6]; //for the six nodes of the penta
    6295 
    6296         /*Initialize Element vector and return if necessary*/
    6297         element->GetInputValue(&approximation,ApproximationEnum);
    6298         if(approximation!=SSAFSApproximationEnum) return NULL;
    6299         int vnumnodes = element->NumberofNodesVelocity();
    6300         int pnumnodes = element->NumberofNodesPressure();
    6301 
    6302         /*Prepare coordinate system list*/
    6303         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    6304         Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
    6305         for(i=0;i<vnumnodes;i++){
    6306                 cs_list[i]             = XYZEnum;
    6307                 node_list[i]           = element->GetNode(i);
    6308         }
    6309         for(i=0;i<pnumnodes;i++){
    6310                 cs_list[vnumnodes+i]   = PressureEnum;
    6311                 node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
    6312         }
    6313         ElementVector* pe=element->NewElementVector(FSvelocityEnum);
    6314 
    6315         /*Retrieve all inputs and parameters*/
    6316         element->GetVerticesCoordinates(&xyz_list);
    6317         element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    6318         Input* vx_input   =element->GetInput(VxEnum);      _assert_(vx_input);
    6319         Input* vy_input   =element->GetInput(VyEnum);      _assert_(vy_input);
    6320         Input* vz_input   =element->GetInput(VzEnum);      _assert_(vz_input);
    6321         Input* vzSSA_input=element->GetInput(VzSSAEnum);   _assert_(vzSSA_input);
    6322 
    6323         /* Start  looping on the number of gaussian points: */
    6324         Gauss* gauss=element->NewGauss(5);
    6325         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6326 
    6327                 gauss->GaussPoint(ig);
    6328                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    6329                 element->NodalFunctionsP1(&basis[0], gauss);
    6330                 element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
    6331 
    6332                 vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
    6333                 element->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
    6334 
    6335                 for(i=0;i<6;i++){
    6336                         pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
    6337                         pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
    6338                         pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
    6339                         pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
    6340                 }
    6341         }
    6342 
    6343         /*Transform coordinate system*/
    6344         element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
    6345 
    6346         /*Clean up and return*/
    6347         xDelete<int>(cs_list);
    6348         xDelete<Node*>(node_list);
    6349         xDelete<IssmDouble>(xyz_list);
    6350         delete gauss;
    6351         return pe;
    6352 }/*}}}*/
    6353 void StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    6354         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2.
    6355          * For node i, Bi can be expressed in the actual coordinate system
     6353void           StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     6354        /*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2.
     6355         * For node i, Bprimei can be expressed in the actual coordinate system
    63566356         * by:
    6357          *       Bi=[ dh/dx          0      ]
    6358          *          [   0           dh/dy  ]
    6359          *          [ 1/2*dh/dy  1/2*dh/dx ]
     6357         *       Bprimei=[ 2*dh/dx    dh/dy   0   0 ]
     6358         *               [  dh/dx    2*dh/dy  0   0 ]
     6359         *               [  dh/dy     dh/dx   0   0 ]
    63606360         * where h is the interpolation function for node i.
    63616361         *
    6362          * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
     6362         * We assume Bprime has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
    63636363         */
    63646364
     6365        int    i;
     6366        IssmDouble dbasismini[3][7];
     6367
     6368        /*Get dbasis in actual coordinate system: */
     6369        element->NodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
     6370
     6371        /*Build Bprime: */
     6372        for(i=0;i<6;i++){
     6373                Bprime[(3*7+6)*0+3*i+0] = 2.*dbasismini[0][i];
     6374                Bprime[(3*7+6)*0+3*i+1] = dbasismini[1][i];
     6375                Bprime[(3*7+6)*0+3*i+2] = 0.;
     6376                Bprime[(3*7+6)*1+3*i+0] = dbasismini[0][i];
     6377                Bprime[(3*7+6)*1+3*i+1] = 2.*dbasismini[1][i];
     6378                Bprime[(3*7+6)*1+3*i+2] = 0.;
     6379                Bprime[(3*7+6)*2+3*i+0] = dbasismini[1][i];
     6380                Bprime[(3*7+6)*2+3*i+1] = dbasismini[0][i];
     6381                Bprime[(3*7+6)*2+3*i+2] = 0.;
     6382        }
     6383
     6384        for(i=0;i<1;i++){ //Add zeros for the bubble function
     6385                Bprime[(3*7+6)*0+3*(6+i)+0] = 0.;
     6386                Bprime[(3*7+6)*0+3*(6+i)+1] = 0.;
     6387                Bprime[(3*7+6)*0+3*(6+i)+2] = 0.;
     6388                Bprime[(3*7+6)*1+3*(6+i)+0] = 0.;
     6389                Bprime[(3*7+6)*1+3*(6+i)+1] = 0.;
     6390                Bprime[(3*7+6)*1+3*(6+i)+2] = 0.;
     6391                Bprime[(3*7+6)*2+3*(6+i)+0] = 0.;
     6392                Bprime[(3*7+6)*2+3*(6+i)+1] = 0.;
     6393                Bprime[(3*7+6)*2+3*(6+i)+2] = 0.;
     6394        }
     6395
     6396        for(i=0;i<6;i++){ //last column not for the bubble function
     6397                Bprime[(3*7+6)*0+7*3+i] = 0.;
     6398                Bprime[(3*7+6)*1+7*3+i] = 0.;
     6399                Bprime[(3*7+6)*2+7*3+i] = 0.;
     6400        }
     6401}/*}}}*/
     6402void           StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     6403        /*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2.
     6404         * For node i, Bprimei can be expressed in the actual coordinate system
     6405         * by:
     6406         *       Bprimei=[  dN/dx    0   ]
     6407         *               [    0    dN/dy ]
     6408         *               [  dN/dy  dN/dx ]
     6409         N               [  dN/dx  dN/dy ]
     6410         * where N is the finiteelement function for node i.
     6411         *
     6412         * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)
     6413         */
     6414
     6415        /*Fetch number of nodes for this finite element*/
    63656416        int numnodes = element->GetNumberOfNodes();
    6366         IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
    6367 
    6368         /*Get dbasis in actual coordinate system: */
     6417
     6418        /*Get nodal functions*/
     6419        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    63696420        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    63706421
    6371         /*Build B: */
     6422        /*Build Bprime: */
    63726423        for(int i=0;i<numnodes;i++){
    6373                 B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
    6374                 B[2*numnodes*0+2*i+1] = 0.;
    6375                 B[2*numnodes*1+2*i+0] = 0.;
    6376                 B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
    6377                 B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
    6378                 B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
     6424                Bprime[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
     6425                Bprime[2*numnodes*0+2*i+1] = 0.;
     6426                Bprime[2*numnodes*1+2*i+0] = 0.;
     6427                Bprime[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
     6428                Bprime[2*numnodes*2+2*i+0] = dbasis[1*numnodes+i];
     6429                Bprime[2*numnodes*2+2*i+1] = dbasis[0*numnodes+i];
     6430                Bprime[2*numnodes*3+2*i+0] = dbasis[0*numnodes+i];
     6431                Bprime[2*numnodes*3+2*i+1] = dbasis[1*numnodes+i];
    63796432        }
    63806433
    … …  
    63826435        xDelete<IssmDouble>(dbasis);
    63836436}/*}}}*/
    6384 void StressbalanceAnalysis::GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     6437void           StressbalanceAnalysis::GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    63856438        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    63866439         * For node i, Bi can be expressed in the actual coordinate system
    … …  
    64406493        }
    64416494}/*}}}*/
    6442 void StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    6443         /*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2.
    6444          * For node i, Bprimei can be expressed in the actual coordinate system
    6445          * by:
    6446          *       Bprimei=[ 2*dh/dx    dh/dy   0   0 ]
    6447          *               [  dh/dx    2*dh/dy  0   0 ]
    6448          *               [  dh/dy     dh/dx   0   0 ]
    6449          * where h is the interpolation function for node i.
    6450          *
    6451          * We assume Bprime has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
    6452          */
    6453 
    6454         int    i;
    6455         IssmDouble dbasismini[3][7];
    6456 
    6457         /*Get dbasis in actual coordinate system: */
    6458         element->NodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
    6459 
    6460         /*Build Bprime: */
    6461         for(i=0;i<6;i++){
    6462                 Bprime[(3*7+6)*0+3*i+0] = 2.*dbasismini[0][i];
    6463                 Bprime[(3*7+6)*0+3*i+1] = dbasismini[1][i];
    6464                 Bprime[(3*7+6)*0+3*i+2] = 0.;
    6465                 Bprime[(3*7+6)*1+3*i+0] = dbasismini[0][i];
    6466                 Bprime[(3*7+6)*1+3*i+1] = 2.*dbasismini[1][i];
    6467                 Bprime[(3*7+6)*1+3*i+2] = 0.;
    6468                 Bprime[(3*7+6)*2+3*i+0] = dbasismini[1][i];
    6469                 Bprime[(3*7+6)*2+3*i+1] = dbasismini[0][i];
    6470                 Bprime[(3*7+6)*2+3*i+2] = 0.;
    6471         }
    6472 
    6473         for(i=0;i<1;i++){ //Add zeros for the bubble function
    6474                 Bprime[(3*7+6)*0+3*(6+i)+0] = 0.;
    6475                 Bprime[(3*7+6)*0+3*(6+i)+1] = 0.;
    6476                 Bprime[(3*7+6)*0+3*(6+i)+2] = 0.;
    6477                 Bprime[(3*7+6)*1+3*(6+i)+0] = 0.;
    6478                 Bprime[(3*7+6)*1+3*(6+i)+1] = 0.;
    6479                 Bprime[(3*7+6)*1+3*(6+i)+2] = 0.;
    6480                 Bprime[(3*7+6)*2+3*(6+i)+0] = 0.;
    6481                 Bprime[(3*7+6)*2+3*(6+i)+1] = 0.;
    6482                 Bprime[(3*7+6)*2+3*(6+i)+2] = 0.;
    6483         }
    6484 
    6485         for(i=0;i<6;i++){ //last column not for the bubble function
    6486                 Bprime[(3*7+6)*0+7*3+i] = 0.;
    6487                 Bprime[(3*7+6)*1+7*3+i] = 0.;
    6488                 Bprime[(3*7+6)*2+7*3+i] = 0.;
    6489         }
    6490 }/*}}}*/
    6491 void StressbalanceAnalysis::GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     6495void           StressbalanceAnalysis::GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    64926496        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    64936497         * For node i, Bi can be expressed in the actual coordinate system
    … …  
    65216525        xDelete<IssmDouble>(dbasis);
    65226526}/*}}}*/
    6523 void StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    6524         /*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2.
    6525          * For node i, Bprimei can be expressed in the actual coordinate system
     6527void           StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     6528        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2.
     6529         * For node i, Bi can be expressed in the actual coordinate system
    65266530         * by:
    6527          *       Bprimei=[  dN/dx    0   ]
    6528          *               [    0    dN/dy ]
    6529          *               [  dN/dy  dN/dx ]
    6530          N               [  dN/dx  dN/dy ]
    6531          * where N is the finiteelement function for node i.
     6531         *       Bi=[ dh/dx          0      ]
     6532         *          [   0           dh/dy   ]
     6533         *          [ 1/2*dh/dy  1/2*dh/dx  ]
     6534         * where h is the interpolation function for node i.
    65326535         *
    6533          * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)
     6536         * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
    65346537         */
    65356538
    6536         /*Fetch number of nodes for this finite element*/
    65376539        int numnodes = element->GetNumberOfNodes();
    6538 
    6539         /*Get nodal functions*/
    6540         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     6540        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
     6541
     6542        /*Get dbasis in actual coordinate system: */
    65416543        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    65426544
    6543         /*Build Bprime: */
     6545        /*Build B: */
    65446546        for(int i=0;i<numnodes;i++){
    6545                 Bprime[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
    6546                 Bprime[2*numnodes*0+2*i+1] = 0.;
    6547                 Bprime[2*numnodes*1+2*i+0] = 0.;
    6548                 Bprime[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
    6549                 Bprime[2*numnodes*2+2*i+0] = dbasis[1*numnodes+i];
    6550                 Bprime[2*numnodes*2+2*i+1] = dbasis[0*numnodes+i];
    6551                 Bprime[2*numnodes*3+2*i+0] = dbasis[0*numnodes+i];
    6552                 Bprime[2*numnodes*3+2*i+1] = dbasis[1*numnodes+i];
     6547                B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
     6548                B[2*numnodes*0+2*i+1] = 0.;
     6549                B[2*numnodes*1+2*i+0] = 0.;
     6550                B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
     6551                B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
     6552                B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
    65536553        }
    65546554
    … …  
    65566556        xDelete<IssmDouble>(dbasis);
    65576557}/*}}}*/
    6558 void StressbalanceAnalysis::GetLSSAFS(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
    6559         /*
    6560          * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    6561          * For node i, Li can be expressed in the actual coordinate system
     6558void           StressbalanceAnalysis::GetLprimeFSSSA(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
     6559        /* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
     6560         * For node i, Lpi can be expressed in the actual coordinate system
    65626561         * by:
    6563          *       Li=[ h    0 ]
    6564          *                    [ 0    h ]
    6565          *                    [ h    0 ]
    6566          *                    [ 0    h ]
    6567          *                    [ h    0 ]
    6568          *                    [ 0    h ]
    6569          *                    [ h    0 ]
    6570          *                    [ 0    h ]
     6562         *       Lpi=[ h    0 ]
     6563         *                     [ 0    h ]
     6564         *                     [ h    0 ]
     6565         *                     [ 0    h ]
    65716566         * where h is the interpolation function for node i.
    65726567         */
    6573 
    65746568        int num_dof=2;
    65756569        IssmDouble basis[3];
    … …  
    65846578        basis[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    65856579
    6586         /*Build LFS: */
     6580        /*Build LprimeFS: */
    65876581        for(int i=0;i<3;i++){
    6588                 LFS[num_dof*3*0+num_dof*i+0] = basis[i];
    6589                 LFS[num_dof*3*0+num_dof*i+1] = 0;
    6590                 LFS[num_dof*3*1+num_dof*i+0] = 0;
    6591                 LFS[num_dof*3*1+num_dof*i+1] = basis[i];
    6592                 LFS[num_dof*3*2+num_dof*i+0] = basis[i];
    6593                 LFS[num_dof*3*2+num_dof*i+1] = 0;
    6594                 LFS[num_dof*3*3+num_dof*i+0] = 0;
    6595                 LFS[num_dof*3*3+num_dof*i+1] = basis[i];
    6596                 LFS[num_dof*3*4+num_dof*i+0] = basis[i];
    6597                 LFS[num_dof*3*4+num_dof*i+1] = 0;
    6598                 LFS[num_dof*3*5+num_dof*i+0] = 0;
    6599                 LFS[num_dof*3*5+num_dof*i+1] = basis[i];
    6600                 LFS[num_dof*3*6+num_dof*i+0] = basis[i];
    6601                 LFS[num_dof*3*6+num_dof*i+1] = 0;
    6602                 LFS[num_dof*3*7+num_dof*i+0] = 0;
    6603                 LFS[num_dof*3*7+num_dof*i+1] = basis[i];
    6604         }
    6605 }/*}}}*/
    6606 void StressbalanceAnalysis::GetLprimeSSAFS(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
     6582                LprimeFS[num_dof*3*0+num_dof*i+0] = basis[i];
     6583                LprimeFS[num_dof*3*0+num_dof*i+1] = 0.;
     6584                LprimeFS[num_dof*3*1+num_dof*i+0] = 0.;
     6585                LprimeFS[num_dof*3*1+num_dof*i+1] = basis[i];
     6586                LprimeFS[num_dof*3*2+num_dof*i+0] = basis[i];
     6587                LprimeFS[num_dof*3*2+num_dof*i+1] = 0.;
     6588                LprimeFS[num_dof*3*3+num_dof*i+0] = 0.;
     6589                LprimeFS[num_dof*3*3+num_dof*i+1] = basis[i];
     6590        }
     6591}/*}}}*/
     6592void           StressbalanceAnalysis::GetLprimeSSAFS(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
    66076593        /* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    66086594         * For node i, Lpi can be expressed in the actual coordinate system
    … …  
    67096695        }
    67106696}/*}}}*/
    6711 void StressbalanceAnalysis::GetLFSSSA(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
     6697void           StressbalanceAnalysis::GetLFSSSA(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
    67126698        /* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    67136699         * For node i, Li can be expressed in the actual coordinate system
    … …  
    67486734        }
    67496735}/*}}}*/
    6750 void StressbalanceAnalysis::GetLprimeFSSSA(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
    6751         /* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    6752          * For node i, Lpi can be expressed in the actual coordinate system
     6736void           StressbalanceAnalysis::GetLSSAFS(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
     6737        /*
     6738         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     6739         * For node i, Li can be expressed in the actual coordinate system
    67536740         * by:
    6754          *       Lpi=[ h    0 ]
    6755          *                     [ 0    h ]
    6756          *                     [ h    0 ]
    6757          *                     [ 0    h ]
     6741         *       Li=[ h    0 ]
     6742         *                    [ 0    h ]
     6743         *                    [ h    0 ]
     6744         *                    [ 0    h ]
     6745         *                    [ h    0 ]
     6746         *                    [ 0    h ]
     6747         *                    [ h    0 ]
     6748         *                    [ 0    h ]
    67586749         * where h is the interpolation function for node i.
    67596750         */
     6751
    67606752        int num_dof=2;
    67616753        IssmDouble basis[3];
    … …  
    67706762        basis[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    67716763
    6772         /*Build LprimeFS: */
     6764        /*Build LFS: */
    67736765        for(int i=0;i<3;i++){
    6774                 LprimeFS[num_dof*3*0+num_dof*i+0] = basis[i];
    6775                 LprimeFS[num_dof*3*0+num_dof*i+1] = 0.;
    6776                 LprimeFS[num_dof*3*1+num_dof*i+0] = 0.;
    6777                 LprimeFS[num_dof*3*1+num_dof*i+1] = basis[i];
    6778                 LprimeFS[num_dof*3*2+num_dof*i+0] = basis[i];
    6779                 LprimeFS[num_dof*3*2+num_dof*i+1] = 0.;
    6780                 LprimeFS[num_dof*3*3+num_dof*i+0] = 0.;
    6781                 LprimeFS[num_dof*3*3+num_dof*i+1] = basis[i];
    6782         }
    6783 }/*}}}*/
    6784 void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
     6766                LFS[num_dof*3*0+num_dof*i+0] = basis[i];
     6767                LFS[num_dof*3*0+num_dof*i+1] = 0;
     6768                LFS[num_dof*3*1+num_dof*i+0] = 0;
     6769                LFS[num_dof*3*1+num_dof*i+1] = basis[i];
     6770                LFS[num_dof*3*2+num_dof*i+0] = basis[i];
     6771                LFS[num_dof*3*2+num_dof*i+1] = 0;
     6772                LFS[num_dof*3*3+num_dof*i+0] = 0;
     6773                LFS[num_dof*3*3+num_dof*i+1] = basis[i];
     6774                LFS[num_dof*3*4+num_dof*i+0] = basis[i];
     6775                LFS[num_dof*3*4+num_dof*i+1] = 0;
     6776                LFS[num_dof*3*5+num_dof*i+0] = 0;
     6777                LFS[num_dof*3*5+num_dof*i+1] = basis[i];
     6778                LFS[num_dof*3*6+num_dof*i+0] = basis[i];
     6779                LFS[num_dof*3*6+num_dof*i+1] = 0;
     6780                LFS[num_dof*3*7+num_dof*i+0] = 0;
     6781                LFS[num_dof*3*7+num_dof*i+1] = basis[i];
     6782        }
     6783}/*}}}*/
     6784void           StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
    67856785
    67866786        int         i;
    … …  
    68796879        xDelete<int>(cs_list);
    68806880}/*}}}*/
    6881 void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/
     6881void           StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/
    68826882
    68836883        int         i;
    … …  
    69856985        xDelete<int>(cs_list);
    69866986}/*}}}*/
    6987 void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
     6987void           StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
    69886988
    69896989        int         i,domaintype;

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

    r18560 r18928  
    1313  public:
    1414                /*Model processing*/
    15                 int  DofsPerNode(int** doflist,int domaintype,int approximation);
    16                 void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
    17                 void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    18                 void CreateNodes(Nodes* nodes,IoModel* iomodel);
    1915                void CreateConstraints(Constraints* constraints,IoModel* iomodel);
    2016                void CreateLoads(Loads* loads, IoModel* iomodel);
     17                void CreateNodes(Nodes* nodes,IoModel* iomodel);
     18                int  DofsPerNode(int** doflist,int domaintype,int approximation);
     19                void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
     20                void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
    2121
    2222                /*Finite element Analysis*/
    … …  
    2626                ElementMatrix* CreateKMatrix(Element* element);
    2727                ElementVector* CreatePVector(Element* element);
    28                 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
    30                 void GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
    31                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    32                 void UpdateConstraints(FemModel* femmodel);
     28                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
     29                void           GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
     30                void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     31                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
     32                void           UpdateConstraints(FemModel* femmodel);
    3333
    3434                /*SSA*/
    3535                ElementMatrix* CreateJacobianMatrixSSA(Element* element);
    3636                ElementMatrix* CreateKMatrixSSA(Element* element);
    37                 ElementMatrix* CreateKMatrixSSAViscous(Element* element);
    3837                ElementMatrix* CreateKMatrixSSAFriction(Element* element);
    3938                ElementMatrix* CreateKMatrixSSALateralFriction(Element* element);
     39                ElementMatrix* CreateKMatrixSSAViscous(Element* element);
    4040                ElementVector* CreatePVectorSSA(Element* element);
     41                ElementVector* CreatePVectorSSAFront(Element* element);
    4142                ElementVector* CreatePVectorSSADrivingStress(Element* element);
    42                 ElementVector* CreatePVectorSSAFront(Element* element);
    43                 void GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    44                 void GetBSSAprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    45                 void GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    46                 void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
     43                void           GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     44                void           GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     45                void           GetBSSAprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     46                void           InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
    4747                /*L1L2*/
    4848                ElementMatrix* CreateKMatrixL1L2(Element* element);
     49                ElementMatrix* CreateKMatrixL1L2Friction(Element* element);
    4950                ElementMatrix* CreateKMatrixL1L2Viscous(Element* element);
    50                 ElementMatrix* CreateKMatrixL1L2Friction(Element* element);
    5151                ElementVector* CreatePVectorL1L2(Element* element);
     52                ElementVector* CreatePVectorL1L2Front(Element* element);
    5253                ElementVector* CreatePVectorL1L2DrivingStress(Element* element);
    53                 ElementVector* CreatePVectorL1L2Front(Element* element);
    54                 void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
     54                void           InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
    5555                /*HO*/
    5656                ElementMatrix* CreateJacobianMatrixHO(Element* element);
    5757                ElementMatrix* CreateKMatrixHO(Element* element);
     58                ElementMatrix* CreateKMatrixHOFriction(Element* element);
    5859                ElementMatrix* CreateKMatrixHOViscous(Element* element);
    59                 ElementMatrix* CreateKMatrixHOFriction(Element* element);
    6060                ElementVector* CreatePVectorHO(Element* element);
     61                ElementVector* CreatePVectorHOFront(Element* element);
    6162                ElementVector* CreatePVectorHODrivingStress(Element* element);
    62                 ElementVector* CreatePVectorHOFront(Element* element);
    63                 void GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    64                 void GetBHOprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    65                 void GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    66                 void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
     63                void           GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     64                void           GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     65                void           GetBHOprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     66                void           InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
    6767                /*FS*/
    6868                ElementVector* CreateDVectorFS(Element* element);
    6969                ElementMatrix* CreateJacobianMatrixFS(Element* element);
    7070                ElementMatrix* CreateKMatrixFS(Element* element);
     71                ElementMatrix* CreateKMatrixFSFriction(Element* element);
     72                ElementMatrix* CreateKMatrixFSShelf(Element* element);
     73                ElementMatrix* CreateKMatrixFSViscous(Element* element);
    7174                ElementMatrix* CreateKMatrixFSViscousLA(Element* element);
    7275                ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
    73                 ElementMatrix* CreateKMatrixFSViscous(Element* element);
    74                 ElementMatrix* CreateKMatrixFSFriction(Element* element);
    75                 ElementMatrix* CreateKMatrixFSShelf(Element* element);
    7676                ElementVector* CreatePVectorFS(Element* element);
     77                ElementVector* CreatePVectorFSFriction(Element* element);
     78                ElementVector* CreatePVectorFSFront(Element* element);
     79                ElementVector* CreatePVectorFSShelf(Element* element);
     80                ElementVector* CreatePVectorFSStress(Element* element);
    7781                ElementVector* CreatePVectorFSViscous(Element* element);
    7882                ElementVector* CreatePVectorFSViscousLA(Element* element);
    7983                ElementVector* CreatePVectorFSViscousXTH(Element* element);
    80                 ElementVector* CreatePVectorFSShelf(Element* element);
    81                 ElementVector* CreatePVectorFSFront(Element* element);
    82                 ElementVector* CreatePVectorFSFriction(Element* element);
    83                 ElementVector* CreatePVectorFSStress(Element* element);
    84                 void GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    85                 void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    86                 void GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    87                 void GetBFSprimevel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    88                 void GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    89                 void GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    90                 void GetBFSprimeUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    91                 void GetCFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    92                 void GetCFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    93                 void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
    94                 void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
    95                 void InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters);
    96                 void InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters);
    97                 void InitializeXTH(Elements* elements,Parameters* parameters);
     84                void           GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     85                void           GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     86                void           GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     87                void           GetBFSprimeUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     88                void           GetBFSprimevel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     89                void           GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     90                void           GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     91                void           GetCFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     92                void           GetCFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     93                void           GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
     94                void           InitializeXTH(Elements* elements,Parameters* parameters);
     95                void           InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
     96                void           InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters);
     97                void           InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters);
    9898                /*Coupling*/
     99                ElementMatrix* CreateKMatrixCouplingHOFS(Element* element);
     100                ElementMatrix* CreateKMatrixCouplingSSAFS(Element* element);
     101                ElementMatrix* CreateKMatrixCouplingSSAFSFriction(Element* element);
     102                ElementMatrix* CreateKMatrixCouplingSSAFSViscous(Element* element);
     103                ElementMatrix* CreateKMatrixCouplingSSAHO(Element* element);
     104                ElementMatrix* CreateKMatrixCouplingSSAHOFriction(Element* element);
     105                ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element);
     106                ElementMatrix* CreateKMatrixHOFS(Element* element);
     107                ElementMatrix* CreateKMatrixSSAFS(Element* element);
     108                ElementMatrix* CreateKMatrixSSAHO(Element* element);
    99109                ElementMatrix* CreateKMatrixSSA3d(Element* element);
    100110                ElementMatrix* CreateKMatrixSSA3dFriction(Element* element);
    101111                ElementMatrix* CreateKMatrixSSA3dViscous(Element* element);
    102                 ElementMatrix* CreateKMatrixHOFS(Element* element);
    103                 ElementMatrix* CreateKMatrixSSAHO(Element* element);
    104                 ElementMatrix* CreateKMatrixSSAFS(Element* element);
    105                 ElementMatrix* CreateKMatrixCouplingHOFS(Element* element);
    106                 ElementMatrix* CreateKMatrixCouplingSSAHO(Element* element);
    107                 ElementMatrix* CreateKMatrixCouplingSSAHOFriction(Element* element);
    108                 ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element);
    109                 ElementMatrix* CreateKMatrixCouplingSSAFS(Element* element);
    110                 ElementMatrix* CreateKMatrixCouplingSSAFSFriction(Element* element);
    111                 ElementMatrix* CreateKMatrixCouplingSSAFSViscous(Element* element);
     112                ElementVector* CreatePVectorSSAFS(Element* element);
    112113                ElementVector* CreatePVectorSSAHO(Element* element);
    113                 ElementVector* CreatePVectorSSAFS(Element* element);
    114114                ElementVector* CreatePVectorCouplingSSAFS(Element* element);
    115115                ElementVector* CreatePVectorCouplingSSAFSFriction(Element* element);
    … …  
    119119                ElementVector* CreatePVectorCouplingHOFSFriction(Element* element);
    120120                ElementVector* CreatePVectorCouplingHOFSViscous(Element* element);
    121                 void GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    122                 void GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    123                 void GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    124                 void GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    125                 void GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    126                 void GetLFSSSA(IssmDouble* L,Element* element,Gauss* gauss);
    127                 void GetLSSAFS(IssmDouble* L,Element* element,Gauss* gauss);
    128                 void GetLprimeFSSSA(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    129                 void GetLprimeSSAFS(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    130                 void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
    131                 void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);
    132                 void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);
     121                void           GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     122                void           GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     123                void           GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     124                void           GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     125                void           GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     126                void           GetLFSSSA(IssmDouble* L,Element* element,Gauss* gauss);
     127                void           GetLprimeFSSSA(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     128                void           GetLprimeSSAFS(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     129                void           GetLSSAFS(IssmDouble* L,Element* element,Gauss* gauss);
     130                void           InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
     131                void           InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);
     132                void           InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);
    133133};
    134134#endif
Note: See TracChangeset for help on using the changeset viewer.