Changeset 16829


Ignore:
Timestamp:
11/19/13 10:19:27 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added KMatrix HO

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r16827 r16829  
    818818                case SSAApproximationEnum:
    819819                        return CreateKMatrixSSA(element);
     820                case HOApproximationEnum:
     821                        return CreateKMatrixHO(element);
    820822                case NoneApproximationEnum:
    821823                        return NULL;
    … …  
    823825                        _error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
    824826        }
    825 }/*}}}*/
    826 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA(Element* element){/*{{{*/
    827 
    828         /*Intermediaries*/
    829         int      meshtype;
    830         Element* basalelement;
    831 
    832         /*Get basal element*/
    833         element->FindParam(&meshtype,MeshTypeEnum);
    834         switch(meshtype){
    835                 case Mesh2DhorizontalEnum:
    836                         basalelement = element;
    837                         break;
    838                 case Mesh3DEnum:
    839                         if(!element->IsOnBed()) return NULL;
    840                         basalelement = element->SpawnBasalElement();
    841                         break;
    842                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    843         }
    844 
    845         /*compute all stiffness matrices for this element*/
    846         ElementMatrix* Ke1=CreateKMatrixSSAViscous(basalelement);
    847         ElementMatrix* Ke2=CreateKMatrixSSAFriction(basalelement);
    848         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    849 
    850         /*clean-up and return*/
    851         if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    852         delete Ke1;
    853         delete Ke2;
    854         return Ke;
    855 }/*}}}*/
    856 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAViscous(Element* element){/*{{{*/
    857 
    858         /*Intermediaries*/
    859         IssmDouble  viscosity,newviscosity,oldviscosity;
    860         IssmDouble  viscosity_overshoot,thickness,Jdet;
    861         IssmDouble  D_scalar;
    862         IssmDouble *xyz_list = NULL;
    863 
    864         /*Fetch number of nodes and dof for this finite element*/
    865         int numnodes = element->GetNumberOfNodes();
    866         int numdof   = numnodes*2;
    867 
    868         /*Initialize Element matrix and vectors*/
    869         ElementMatrix* Ke     = element->NewElementMatrix(SSAApproximationEnum);
    870         IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
    871         IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
    872         IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
    873 
    874         /*Retrieve all inputs and parameters*/
    875         element->GetVerticesCoordinates(&xyz_list);
    876         Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
    877         Input* vx_input=element->GetInput(VxEnum);               _assert_(vx_input);
    878         Input* vy_input=element->GetInput(VyEnum);               _assert_(vy_input);
    879         Input* vxold_input=element->GetInput(VxPicardEnum);      _assert_(vxold_input);
    880         Input* vyold_input=element->GetInput(VyPicardEnum);      _assert_(vyold_input);
    881         element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
    882 
    883         /* Start  looping on the number of gaussian points: */
    884         Gauss* gauss = element->NewGauss(2);
    885         for(int ig=gauss->begin();ig<gauss->end();ig++){
    886                 gauss->GaussPoint(ig);
    887 
    888                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    889                 this->GetBSSA(B,element,xyz_list,gauss);
    890                 this->GetBSSAprime(Bprime,element,xyz_list,gauss);
    891 
    892                 element->ViscositySSA(&viscosity,xyz_list,gauss,vx_input,vy_input);
    893                 element->ViscositySSA(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
    894                 thickness_input->GetInputValue(&thickness, gauss);
    895 
    896                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    897                 D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet;
    898                 for(int i=0;i<3;i++) D[i*3+i]=D_scalar;
    899 
    900                 TripleMultiply(B,3,numdof,1,
    901                                         D,3,3,0,
    902                                         Bprime,3,numdof,0,
    903                                         &Ke->values[0],1);
    904         }
    905 
    906         /*Transform Coordinate System*/
    907         element->TransformStiffnessMatrixCoord(Ke,XYEnum);
    908 
    909         /*Clean up and return*/
    910         delete gauss;
    911         xDelete<IssmDouble>(xyz_list);
    912         xDelete<IssmDouble>(D);
    913         xDelete<IssmDouble>(Bprime);
    914         xDelete<IssmDouble>(B);
    915         return Ke;
    916 }/*}}}*/
    917 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
    918         return NULL;
    919827}/*}}}*/
    920828ElementVector* StressbalanceAnalysis::CreatePVector(Element* element){/*{{{*/
    … …  
    934842                        _error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
    935843        }
    936 }/*}}}*/
    937 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
    938 
    939         /*compute all load vectors for this element*/
    940         ElementVector* pe1=CreatePVectorFSViscous(element);
    941         ElementVector* pe2=CreatePVectorFSShelf(element);
    942         ElementVector* pe3=CreatePVectorFSFront(element);
    943         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    944 
    945         /*clean-up and return*/
    946         delete pe1;
    947         delete pe2;
    948         delete pe3;
    949         return pe;
    950 }/*}}}*/
    951 ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
    952 
    953         int         i,meshtype,dim;
    954         IssmDouble  Jdet,forcex,forcey,forcez;
    955         IssmDouble *xyz_list = NULL;
    956 
    957         /*Get problem dimension*/
    958         element->FindParam(&meshtype,MeshTypeEnum);
    959         switch(meshtype){
    960                 case Mesh2DverticalEnum: dim = 2; break;
    961                 case Mesh3DEnum:         dim = 3; break;
    962                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    963         }
    964 
    965         /*Fetch number of nodes and dof for this finite element*/
    966         int vnumnodes = element->GetNumberOfNodesVelocity();
    967         int pnumnodes = element->GetNumberOfNodesPressure();
    968 
    969         /*Prepare coordinate system list*/
    970         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    971         if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
    972         else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
    973         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    974 
    975         /*Initialize vectors*/
    976         ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
    977         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    978 
    979         /*Retrieve all inputs and parameters*/
    980         element->GetVerticesCoordinates(&xyz_list);
    981         Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
    982         Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
    983         Input*      loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
    984         IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
    985         IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
    986 
    987         /* Start  looping on the number of gaussian points: */
    988         Gauss* gauss=element->NewGauss(5);
    989         for(int ig=gauss->begin();ig<gauss->end();ig++){
    990                 gauss->GaussPoint(ig);
    991 
    992                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    993                 element->NodalFunctionsVelocity(vbasis,gauss);
    994 
    995                 loadingforcex_input->GetInputValue(&forcex,gauss);
    996                 loadingforcey_input->GetInputValue(&forcey,gauss);
    997                 if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss);
    998 
    999                 for(i=0;i<vnumnodes;i++){
    1000                         pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
    1001                         pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
    1002                         if(dim==3){
    1003                                 pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
    1004                                 pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
    1005                         }
    1006                         else{
    1007                                 pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
    1008                         }
    1009                 }
    1010         }
    1011 
    1012         /*Transform coordinate system*/
    1013         element->TransformLoadVectorCoord(pe,cs_list);
    1014 
    1015         /*Clean up and return*/
    1016         delete gauss;
    1017         xDelete<int>(cs_list);
    1018         xDelete<IssmDouble>(vbasis);
    1019         xDelete<IssmDouble>(xyz_list);
    1020         return pe;
    1021 }/*}}}*/
    1022 ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/
    1023 
    1024         int         i,meshtype,dim;
    1025         IssmDouble  Jdet,water_pressure,bed;
    1026         IssmDouble      normal[3];
    1027         IssmDouble *xyz_list_base = NULL;
    1028 
    1029         /*Get basal element*/
    1030         if(!element->IsOnBed() || !element->IsFloating()) return NULL;
    1031 
    1032         /*Get problem dimension*/
    1033         element->FindParam(&meshtype,MeshTypeEnum);
    1034         switch(meshtype){
    1035                 case Mesh2DverticalEnum: dim = 2; break;
    1036                 case Mesh3DEnum:         dim = 3; break;
    1037                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    1038         }
    1039 
    1040         /*Fetch number of nodes and dof for this finite element*/
    1041         int vnumnodes = element->GetNumberOfNodesVelocity();
    1042         int pnumnodes = element->GetNumberOfNodesPressure();
    1043 
    1044         /*Prepare coordinate system list*/
    1045         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    1046         if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
    1047         else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
    1048         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    1049 
    1050         /*Initialize vectors*/
    1051         ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
    1052         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    1053 
    1054         /*Retrieve all inputs and parameters*/
    1055         element->GetVerticesCoordinatesBase(&xyz_list_base);
    1056         Input*      bed_input=element->GetInput(BedEnum); _assert_(bed_input);
    1057         IssmDouble  rho_water=element->GetMaterialParameter(MaterialsRhoWaterEnum);
    1058         IssmDouble  gravity  =element->GetMaterialParameter(ConstantsGEnum);
    1059 
    1060         /* Start  looping on the number of gaussian points: */
    1061         Gauss* gauss=element->NewGaussBase(5);
    1062         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1063                 gauss->GaussPoint(ig);
    1064 
    1065                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    1066                 element->NodalFunctionsVelocity(vbasis,gauss);
    1067 
    1068                 element->NormalBase(&normal[0],xyz_list_base);
    1069                 _assert_(normal[dim-1]<0.);
    1070                 bed_input->GetInputValue(&bed, gauss);
    1071                 water_pressure=gravity*rho_water*bed;
    1072 
    1073                 for(i=0;i<vnumnodes;i++){
    1074                         pe->values[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0];
    1075                         pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1];
    1076                         if(dim==3){
    1077                                 pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2];
    1078                         }
    1079                 }
    1080         }
    1081 
    1082         /*Transform coordinate system*/
    1083         element->TransformLoadVectorCoord(pe,cs_list);
    1084 
    1085         /*Clean up and return*/
    1086         delete gauss;
    1087         xDelete<int>(cs_list);
    1088         xDelete<IssmDouble>(vbasis);
    1089         xDelete<IssmDouble>(xyz_list_base);
    1090         return pe;
    1091 }/*}}}*/
    1092 ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/
    1093 
    1094         return NULL;
    1095 }/*}}}*/
    1096 ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
    1097 
    1098         /*compute all load vectors for this element*/
    1099         ElementVector* pe1=CreatePVectorHODrivingStress(element);
    1100         ElementVector* pe2=CreatePVectorHOFront(element);
    1101         ElementVector* pe =new ElementVector(pe1,pe2);
    1102 
    1103         /*clean-up and return*/
    1104         delete pe1;
    1105         delete pe2;
    1106         return pe;
    1107 }/*}}}*/
    1108 ElementVector* StressbalanceAnalysis::CreatePVectorHODrivingStress(Element* element){/*{{{*/
    1109 
    1110         /*Intermediaries */
    1111         IssmDouble  Jdet,slope[3];
    1112         IssmDouble* xyz_list = NULL;
    1113 
    1114         /*Fetch number of nodes and dof for this finite element*/
    1115         int numnodes = element->GetNumberOfNodes();
    1116 
    1117         /*Initialize Element vector and vectors*/
    1118         ElementVector* pe=element->NewElementVector(HOApproximationEnum);
    1119         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    1120 
    1121         /*Retrieve all inputs and parameters*/
    1122         element->GetVerticesCoordinates(&xyz_list);
    1123         Input*     surface_input = element->GetInput(SurfaceEnum);   _assert_(surface_input);
    1124         IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
    1125 
    1126         /* Start  looping on the number of gaussian points: */
    1127         Gauss* gauss=element->NewGauss(3);
    1128         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1129                 gauss->GaussPoint(ig);
    1130 
    1131                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    1132                 element->NodalFunctions(basis, gauss);
    1133                 surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
    1134 
    1135                 for(int i=0;i<numnodes;i++){
    1136                         pe->values[i*2+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i];
    1137                         pe->values[i*2+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i];
    1138                 }
    1139         }
    1140 
    1141         /*Transform coordinate system*/
    1142         element->TransformLoadVectorCoord(pe,XYEnum);
    1143 
    1144         /*Clean up and return*/
    1145         xDelete<IssmDouble>(basis);
    1146         xDelete<IssmDouble>(xyz_list);
    1147         delete gauss;
    1148         return pe;
    1149 }/*}}}*/
    1150 ElementVector* StressbalanceAnalysis::CreatePVectorHOFront(Element* element){/*{{{*/
    1151 
    1152         return NULL;
    1153 }/*}}}*/
    1154 ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
    1155 
    1156         /*Intermediaries*/
    1157         int      meshtype;
    1158         Element* basalelement;
    1159 
    1160         /*Get basal element*/
    1161         element->FindParam(&meshtype,MeshTypeEnum);
    1162         switch(meshtype){
    1163                 case Mesh2DhorizontalEnum:
    1164                         basalelement = element;
    1165                         break;
    1166                 case Mesh3DEnum:
    1167                         if(!element->IsOnBed()) return NULL;
    1168                         basalelement = element->SpawnBasalElement();
    1169                         break;
    1170                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    1171         }
    1172 
    1173         /*compute all load vectors for this element*/
    1174         ElementVector* pe1=CreatePVectorSSADrivingStress(basalelement);
    1175         ElementVector* pe2=CreatePVectorSSAFront(basalelement);
    1176         ElementVector* pe =new ElementVector(pe1,pe2);
    1177 
    1178         /*clean-up and return*/
    1179         if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    1180         delete pe1;
    1181         delete pe2;
    1182         return pe;
    1183 }/*}}}*/
    1184 ElementVector* StressbalanceAnalysis::CreatePVectorSSADrivingStress(Element* element){/*{{{*/
    1185 
    1186         /*Intermediaries */
    1187         IssmDouble  thickness,Jdet,slope[2];
    1188         IssmDouble* xyz_list = NULL;
    1189 
    1190         /*Fetch number of nodes and dof for this finite element*/
    1191         int numnodes = element->GetNumberOfNodes();
    1192 
    1193         /*Initialize Element vector and vectors*/
    1194         ElementVector* pe    = element->NewElementVector(SSAApproximationEnum);
    1195         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    1196 
    1197         /*Retrieve all inputs and parameters*/
    1198         element->GetVerticesCoordinates(&xyz_list);
    1199         Input*     thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
    1200         Input*     surface_input  =element->GetInput(SurfaceEnum);   _assert_(surface_input);
    1201         IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
    1202 
    1203         /* Start  looping on the number of gaussian points: */
    1204         Gauss* gauss=element->NewGauss(2);
    1205         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1206 
    1207                 gauss->GaussPoint(ig);
    1208 
    1209                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    1210                 element->NodalFunctions(basis, gauss);
    1211 
    1212                 thickness_input->GetInputValue(&thickness,gauss);
    1213                 surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
    1214 
    1215                 for(int i=0;i<numnodes;i++){
    1216                         pe->values[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i];
    1217                         pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i];
    1218                 }
    1219         }
    1220 
    1221         /*Transform coordinate system*/
    1222         element->TransformLoadVectorCoord(pe,XYEnum);
    1223 
    1224         /*Clean up and return*/
    1225         xDelete<IssmDouble>(xyz_list);
    1226         xDelete<IssmDouble>(basis);
    1227         delete gauss;
    1228         return pe;
    1229 }/*}}}*/
    1230 ElementVector* StressbalanceAnalysis::CreatePVectorSSAFront(Element* element){/*{{{*/
    1231 
    1232         return NULL;
    1233 }/*}}}*/
    1234 void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    1235         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    1236          * For node i, Bi can be expressed in the actual coordinate system
    1237          * by:
    1238          *       Bi=[ dN/dx           0    ]
    1239          *          [   0           dN/dy  ]
    1240          *          [ 1/2*dN/dy  1/2*dN/dx ]
    1241          * where N is the finiteelement function for node i.
    1242          *
    1243          * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
    1244          */
    1245 
    1246         /*Fetch number of nodes for this finite element*/
    1247         int numnodes = element->GetNumberOfNodes();
    1248 
    1249         /*Get nodal functions derivatives*/
    1250         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    1251         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    1252 
    1253         /*Build B: */
    1254         for(int i=0;i<numnodes;i++){
    1255                 B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
    1256                 B[2*numnodes*0+2*i+1] = 0.;
    1257                 B[2*numnodes*1+2*i+0] = 0.;
    1258                 B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
    1259                 B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
    1260                 B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
    1261         }
    1262 
    1263         /*Clean-up*/
    1264         xDelete<IssmDouble>(dbasis);
    1265 }/*}}}*/
    1266 void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    1267         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    1268          * For node i, Bi' can be expressed in the actual coordinate system
    1269          * by:
    1270          *       Bi_prime=[ 2*dN/dx    dN/dy ]
    1271          *                [   dN/dx  2*dN/dy ]
    1272          *                [   dN/dy    dN/dx ]
    1273          * where hNis the finiteelement function for node i.
    1274          *
    1275          * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
    1276          */
    1277 
    1278         /*Fetch number of nodes for this finite element*/
    1279         int numnodes = element->GetNumberOfNodes();
    1280 
    1281         /*Get nodal functions derivatives*/
    1282         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    1283         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    1284 
    1285         /*Build B': */
    1286         for(int i=0;i<numnodes;i++){
    1287                 Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
    1288                 Bprime[2*numnodes*0+2*i+1] =    dbasis[1*numnodes+i];
    1289                 Bprime[2*numnodes*1+2*i+0] =    dbasis[0*numnodes+i];
    1290                 Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
    1291                 Bprime[2*numnodes*2+2*i+0] =    dbasis[1*numnodes+i];
    1292                 Bprime[2*numnodes*2+2*i+1] =    dbasis[0*numnodes+i];
    1293         }
    1294 
    1295         /*Clean-up*/
    1296         xDelete<IssmDouble>(dbasis);
    1297844}/*}}}*/
    1298845void StressbalanceAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    … …  
    1317864        }
    1318865}/*}}}*/
    1319 void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    1320 
    1321         int*         vdoflist=NULL;
    1322         int*         pdoflist=NULL;
    1323         Input*       vz_input=NULL;
    1324         int          meshtype,dim;
    1325         IssmDouble   vx,vy,vz,p;
    1326         IssmDouble   FSreconditioning;
    1327 
    1328         /*Get some parameters*/
    1329         element->FindParam(&meshtype,MeshTypeEnum);
    1330         element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    1331         switch(meshtype){
    1332                 case Mesh2DverticalEnum: dim = 2; break;
    1333                 case Mesh3DEnum:         dim = 3; break;
    1334                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    1335         }
    1336 
    1337         /*Fetch number of nodes and dof for this finite element*/
    1338         int vnumnodes = element->NumberofNodesVelocity();
    1339         int pnumnodes = element->NumberofNodesPressure();
    1340         int vnumdof   = vnumnodes*dim;
    1341         int pnumdof   = pnumnodes*1;
    1342 
    1343         /*Initialize values*/
    1344         IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
    1345         IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
    1346 
    1347         /*Get dof list: */
    1348         element->GetDofListVelocity(&vdoflist,GsetEnum);
    1349         element->GetDofListPressure(&pdoflist,GsetEnum);
    1350         Input*     vx_input=element->GetInput(VxEnum);       _assert_(vx_input);
    1351         Input*     vy_input=element->GetInput(VyEnum);       _assert_(vy_input);
    1352         if(dim==3){vz_input=element->GetInput(VzEnum);       _assert_(vz_input);}
    1353         Input*     p_input =element->GetInput(PressureEnum); _assert_(p_input);
    1354 
    1355         element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    1356 
    1357         /*Ok, we have the velocities in inputs, fill in solution */
    1358         Gauss* gauss = element->NewGauss();
    1359         for(int i=0;i<vnumnodes;i++){
    1360                 gauss->GaussNode(element->VelocityInterpolation(),i);
    1361                 vx_input->GetInputValue(&vx,gauss);
    1362                 vy_input->GetInputValue(&vy,gauss);
    1363                 vvalues[i*dim+0]=vx;
    1364                 vvalues[i*dim+1]=vy;
    1365                 if(dim==3){
    1366                         vz_input->GetInputValue(&vz,gauss);
    1367                         vvalues[i*dim+2]=vz;
    1368                 }
    1369         }
    1370         for(int i=0;i<pnumnodes;i++){
    1371                 gauss->GaussNode(element->PressureInterpolation(),i);
    1372                 p_input->GetInputValue(&p ,gauss);
    1373                 pvalues[i]=p/FSreconditioning;
    1374         }
    1375 
    1376         /*Add value to global vector*/
    1377         solution->SetValues(vnumdof,vdoflist,vvalues,INS_VAL);
    1378         solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL);
    1379 
    1380         /*Free ressources:*/
    1381         delete gauss;
    1382         xDelete<int>(pdoflist);
    1383         xDelete<int>(vdoflist);
    1384         xDelete<IssmDouble>(pvalues);
    1385         xDelete<IssmDouble>(vvalues);
    1386 }/*}}}*/
    1387866void StressbalanceAnalysis::GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    1388867
    … …  
    1412891                vx_input->GetInputValue(&vx,gauss);
    1413892                vy_input->GetInputValue(&vy,gauss);
    1414                 values[i*NDOF2+0]=vx;
    1415                 values[i*NDOF2+1]=vy;
     893                values[i*2+0]=vx;
     894                values[i*2+1]=vy;
    1416895        }
    1417896
    … …  
    1455934        }
    1456935}/*}}}*/
    1457 void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
    1458 
    1459         int          i,dim,meshtype;
    1460         int*         vdoflist=NULL;
    1461         int*         pdoflist=NULL;
    1462         IssmDouble   FSreconditioning;
    1463 
     936
     937/*SSA*/
     938ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA(Element* element){/*{{{*/
     939
     940        /*Intermediaries*/
     941        int      meshtype;
     942        Element* basalelement;
     943
     944        /*Get basal element*/
    1464945        element->FindParam(&meshtype,MeshTypeEnum);
    1465         element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    1466946        switch(meshtype){
    1467                 case Mesh2DverticalEnum: dim = 2; break;
    1468                 case Mesh3DEnum:         dim = 3; break;
     947                case Mesh2DhorizontalEnum:
     948                        basalelement = element;
     949                        break;
     950                case Mesh3DEnum:
     951                        if(!element->IsOnBed()) return NULL;
     952                        basalelement = element->SpawnBasalElement();
     953                        break;
    1469954                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    1470955        }
    1471956
     957        /*compute all stiffness matrices for this element*/
     958        ElementMatrix* Ke1=CreateKMatrixSSAViscous(basalelement);
     959        ElementMatrix* Ke2=CreateKMatrixSSAFriction(basalelement);
     960        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     961
     962        /*clean-up and return*/
     963        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     964        delete Ke1;
     965        delete Ke2;
     966        return Ke;
     967}/*}}}*/
     968ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAViscous(Element* element){/*{{{*/
     969
     970        /*Intermediaries*/
     971        IssmDouble  viscosity,newviscosity,oldviscosity;
     972        IssmDouble  viscosity_overshoot,thickness,Jdet;
     973        IssmDouble  D_scalar;
     974        IssmDouble *xyz_list = NULL;
     975
    1472976        /*Fetch number of nodes and dof for this finite element*/
    1473         int vnumnodes = element->GetNumberOfNodesVelocity();
    1474         int pnumnodes = element->GetNumberOfNodesPressure();
    1475         int vnumdof   = vnumnodes*dim;
    1476         int pnumdof   = pnumnodes*1;
    1477 
    1478         /*Initialize values*/
    1479         IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
    1480         IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
    1481         IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
    1482         IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
    1483         IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
    1484         IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
    1485 
    1486         /*Prepare coordinate system list*/
    1487         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    1488         if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
    1489         else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
    1490         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    1491 
    1492         /*Get dof list: */
    1493         element->GetDofListVelocity(&vdoflist,GsetEnum);
    1494         element->GetDofListPressure(&pdoflist,GsetEnum);
     977        int numnodes = element->GetNumberOfNodes();
     978        int numdof   = numnodes*2;
     979
     980        /*Initialize Element matrix and vectors*/
     981        ElementMatrix* Ke     = element->NewElementMatrix(SSAApproximationEnum);
     982        IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
     983        IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
     984        IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
     985
     986        /*Retrieve all inputs and parameters*/
     987        element->GetVerticesCoordinates(&xyz_list);
     988        Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
     989        Input* vx_input=element->GetInput(VxEnum);               _assert_(vx_input);
     990        Input* vy_input=element->GetInput(VyEnum);               _assert_(vy_input);
     991        Input* vxold_input=element->GetInput(VxPicardEnum);      _assert_(vxold_input);
     992        Input* vyold_input=element->GetInput(VyPicardEnum);      _assert_(vyold_input);
     993        element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
     994
     995        /* Start  looping on the number of gaussian points: */
     996        Gauss* gauss = element->NewGauss(2);
     997        for(int ig=gauss->begin();ig<gauss->end();ig++){
     998                gauss->GaussPoint(ig);
     999
     1000                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1001                this->GetBSSA(B,element,xyz_list,gauss);
     1002                this->GetBSSAprime(Bprime,element,xyz_list,gauss);
     1003
     1004                element->ViscositySSA(&viscosity,xyz_list,gauss,vx_input,vy_input);
     1005                element->ViscositySSA(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
     1006                thickness_input->GetInputValue(&thickness, gauss);
     1007
     1008                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
     1009                D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet;
     1010                for(int i=0;i<3;i++) D[i*3+i]=D_scalar;
     1011
     1012                TripleMultiply(B,3,numdof,1,
     1013                                        D,3,3,0,
     1014                                        Bprime,3,numdof,0,
     1015                                        &Ke->values[0],1);
     1016        }
     1017
     1018        /*Transform Coordinate System*/
     1019        element->TransformStiffnessMatrixCoord(Ke,XYEnum);
     1020
     1021        /*Clean up and return*/
     1022        delete gauss;
     1023        xDelete<IssmDouble>(xyz_list);
     1024        xDelete<IssmDouble>(D);
     1025        xDelete<IssmDouble>(Bprime);
     1026        xDelete<IssmDouble>(B);
     1027        return Ke;
     1028}/*}}}*/
     1029ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
     1030        return NULL;
     1031}/*}}}*/
     1032ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
     1033
     1034        /*Intermediaries*/
     1035        int      meshtype;
     1036        Element* basalelement;
     1037
     1038        /*Get basal element*/
     1039        element->FindParam(&meshtype,MeshTypeEnum);
     1040        switch(meshtype){
     1041                case Mesh2DhorizontalEnum:
     1042                        basalelement = element;
     1043                        break;
     1044                case Mesh3DEnum:
     1045                        if(!element->IsOnBed()) return NULL;
     1046                        basalelement = element->SpawnBasalElement();
     1047                        break;
     1048                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1049        }
     1050
     1051        /*compute all load vectors for this element*/
     1052        ElementVector* pe1=CreatePVectorSSADrivingStress(basalelement);
     1053        ElementVector* pe2=CreatePVectorSSAFront(basalelement);
     1054        ElementVector* pe =new ElementVector(pe1,pe2);
     1055
     1056        /*clean-up and return*/
     1057        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     1058        delete pe1;
     1059        delete pe2;
     1060        return pe;
     1061}/*}}}*/
     1062ElementVector* StressbalanceAnalysis::CreatePVectorSSADrivingStress(Element* element){/*{{{*/
     1063
     1064        /*Intermediaries */
     1065        IssmDouble  thickness,Jdet,slope[2];
     1066        IssmDouble* xyz_list = NULL;
     1067
     1068        /*Fetch number of nodes and dof for this finite element*/
     1069        int numnodes = element->GetNumberOfNodes();
     1070
     1071        /*Initialize Element vector and vectors*/
     1072        ElementVector* pe    = element->NewElementVector(SSAApproximationEnum);
     1073        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     1074
     1075        /*Retrieve all inputs and parameters*/
     1076        element->GetVerticesCoordinates(&xyz_list);
     1077        Input*     thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
     1078        Input*     surface_input  =element->GetInput(SurfaceEnum);   _assert_(surface_input);
     1079        IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
     1080
     1081        /* Start  looping on the number of gaussian points: */
     1082        Gauss* gauss=element->NewGauss(2);
     1083        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1084
     1085                gauss->GaussPoint(ig);
     1086
     1087                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1088                element->NodalFunctions(basis, gauss);
     1089
     1090                thickness_input->GetInputValue(&thickness,gauss);
     1091                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
     1092
     1093                for(int i=0;i<numnodes;i++){
     1094                        pe->values[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i];
     1095                        pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i];
     1096                }
     1097        }
     1098
     1099        /*Transform coordinate system*/
     1100        element->TransformLoadVectorCoord(pe,XYEnum);
     1101
     1102        /*Clean up and return*/
     1103        xDelete<IssmDouble>(xyz_list);
     1104        xDelete<IssmDouble>(basis);
     1105        delete gauss;
     1106        return pe;
     1107}/*}}}*/
     1108ElementVector* StressbalanceAnalysis::CreatePVectorSSAFront(Element* element){/*{{{*/
     1109
     1110        return NULL;
     1111}/*}}}*/
     1112void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1113        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
     1114         * For node i, Bi can be expressed in the actual coordinate system
     1115         * by:
     1116         *       Bi=[ dN/dx           0    ]
     1117         *          [   0           dN/dy  ]
     1118         *          [ 1/2*dN/dy  1/2*dN/dx ]
     1119         * where N is the finiteelement function for node i.
     1120         *
     1121         * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
     1122         */
     1123
     1124        /*Fetch number of nodes for this finite element*/
     1125        int numnodes = element->GetNumberOfNodes();
     1126
     1127        /*Get nodal functions derivatives*/
     1128        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     1129        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1130
     1131        /*Build B: */
     1132        for(int i=0;i<numnodes;i++){
     1133                B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
     1134                B[2*numnodes*0+2*i+1] = 0.;
     1135                B[2*numnodes*1+2*i+0] = 0.;
     1136                B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
     1137                B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
     1138                B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
     1139        }
     1140
     1141        /*Clean-up*/
     1142        xDelete<IssmDouble>(dbasis);
     1143}/*}}}*/
     1144void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1145        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
     1146         * For node i, Bi' can be expressed in the actual coordinate system
     1147         * by:
     1148         *       Bi_prime=[ 2*dN/dx    dN/dy ]
     1149         *                [   dN/dx  2*dN/dy ]
     1150         *                [   dN/dy    dN/dx ]
     1151         * where hNis the finiteelement function for node i.
     1152         *
     1153         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
     1154         */
     1155
     1156        /*Fetch number of nodes for this finite element*/
     1157        int numnodes = element->GetNumberOfNodes();
     1158
     1159        /*Get nodal functions derivatives*/
     1160        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     1161        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1162
     1163        /*Build B': */
     1164        for(int i=0;i<numnodes;i++){
     1165                Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
     1166                Bprime[2*numnodes*0+2*i+1] =    dbasis[1*numnodes+i];
     1167                Bprime[2*numnodes*1+2*i+0] =    dbasis[0*numnodes+i];
     1168                Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
     1169                Bprime[2*numnodes*2+2*i+0] =    dbasis[1*numnodes+i];
     1170                Bprime[2*numnodes*2+2*i+1] =    dbasis[0*numnodes+i];
     1171        }
     1172
     1173        /*Clean-up*/
     1174        xDelete<IssmDouble>(dbasis);
     1175}/*}}}*/
     1176void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
     1177
     1178        int         i,meshtype;
     1179        IssmDouble  rho_ice,g;
     1180        int*        doflist=NULL;
     1181        IssmDouble* xyz_list=NULL;
     1182        Element*    basalelement=NULL;
     1183
     1184        /*Deal with pressure first*/
     1185        int numvertices = element->GetNumberOfVertices();
     1186        IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
     1187        IssmDouble* thickness = xNew<IssmDouble>(numvertices);
     1188        IssmDouble* surface   = xNew<IssmDouble>(numvertices);
     1189
     1190        element->FindParam(&meshtype,MeshTypeEnum);
     1191        rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
     1192        g       =element->GetMaterialParameter(ConstantsGEnum);
     1193        switch(meshtype){
     1194                case Mesh2DhorizontalEnum:
     1195                        element->GetInputListOnVertices(thickness,ThicknessEnum);
     1196                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
     1197                        break;
     1198                case Mesh3DEnum:
     1199                        element->GetVerticesCoordinates(&xyz_list);
     1200                        element->GetInputListOnVertices(surface,SurfaceEnum);
     1201                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
     1202                        break;
     1203                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1204        }
     1205        element->AddInput(PressureEnum,pressure,P1Enum);
     1206        xDelete<IssmDouble>(pressure);
     1207        xDelete<IssmDouble>(thickness);
     1208        xDelete<IssmDouble>(surface);
     1209
     1210        /*Get basal element*/
     1211        switch(meshtype){
     1212                case Mesh2DhorizontalEnum:
     1213                        basalelement = element;
     1214                        break;
     1215                case Mesh3DEnum:
     1216                        if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
     1217                        basalelement=element->SpawnBasalElement();
     1218                        break;
     1219                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1220        }
     1221
     1222        /*Fetch number of nodes and dof for this finite element*/
     1223        int numnodes = basalelement->GetNumberOfNodes();
     1224        int numdof   = numnodes*2;
     1225
     1226        /*Fetch dof list and allocate solution vectors*/
     1227        basalelement->GetDofList(&doflist,SSAApproximationEnum,GsetEnum);
     1228        IssmDouble* values    = xNew<IssmDouble>(numdof);
     1229        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
     1230        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
     1231        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
     1232        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
    14951233
    14961234        /*Use the dof list to index into the solution vector: */
    1497         for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
    1498         for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
     1235        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    14991236
    15001237        /*Transform solution in Cartesian Space*/
    1501         element->TransformSolutionCoord(values,cs_list);
    1502 
    1503         /*Ok, we have vx and vy in values, fill in all arrays: */
    1504         for(i=0;i<vnumnodes;i++){
    1505                 vx[i] = values[i*dim+0];
    1506                 vy[i] = values[i*dim+1];
     1238        basalelement->TransformSolutionCoord(&values[0],XYEnum);
     1239        basalelement->FindParam(&meshtype,MeshTypeEnum);
     1240
     1241        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     1242        for(i=0;i<numnodes;i++){
     1243                vx[i]=values[i*2+0];
     1244                vy[i]=values[i*2+1];
     1245
     1246                /*Check solution*/
    15071247                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    15081248                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    1509 
    1510                 if(dim==3){
    1511                         vz[i] = values[i*dim+2];
    1512                         if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
    1513                 }
    1514         }
    1515         for(i=0;i<pnumnodes;i++){
    1516                 pressure[i] = values[vnumdof+i];
    1517                 if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    1518         }
    1519 
    1520         /*Recondition pressure and compute vel: */
    1521         for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
    1522         if(dim==3) for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    1523         else       for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
    1524 
    1525         /*Now, we have to move the previous inputs  to old
     1249        }
     1250
     1251        /*Get Vz and compute vel*/
     1252        basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
     1253        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     1254
     1255        /*Now, we have to move the previous Vx and Vy inputs  to old
    15261256         * status, otherwise, we'll wipe them off: */
    15271257        element->InputChangeName(VxEnum,VxPicardEnum);
    15281258        element->InputChangeName(VyEnum,VyPicardEnum);
    1529         element->InputChangeName(PressureEnum,PressurePicardEnum);
    1530         if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);
    15311259
    15321260        /*Add vx and vy as inputs to the tria element: */
    1533         element->AddInput(VxEnum,vx,P1Enum);
    1534         element->AddInput(VyEnum,vy,P1Enum);
    1535         element->AddInput(VelEnum,vel,P1Enum);
    1536         element->AddInput(PressureEnum,pressure,P1Enum);
    1537         if(dim==3) element->AddInput(VzEnum,vz,P1Enum);
     1261        element->AddBasalInput(VxEnum,vx,P1Enum);
     1262        element->AddBasalInput(VyEnum,vy,P1Enum);
     1263        element->AddBasalInput(VelEnum,vel,P1Enum);
    15381264
    15391265        /*Free ressources:*/
    1540         xDelete<IssmDouble>(pressure);
    15411266        xDelete<IssmDouble>(vel);
    15421267        xDelete<IssmDouble>(vz);
    … …  
    15441269        xDelete<IssmDouble>(vx);
    15451270        xDelete<IssmDouble>(values);
    1546         xDelete<int>(vdoflist);
    1547         xDelete<int>(pdoflist);
    1548         xDelete<int>(cs_list);
     1271        xDelete<IssmDouble>(xyz_list);
     1272        xDelete<int>(doflist);
     1273        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     1274}/*}}}*/
     1275
     1276/*L1L2*/
     1277void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
     1278
     1279        int         i,meshtype;
     1280        IssmDouble  rho_ice,g;
     1281        int*        doflist=NULL;
     1282        IssmDouble* xyz_list=NULL;
     1283        Element*    basalelement=NULL;
     1284
     1285        /*Deal with pressure first*/
     1286        int numvertices = element->GetNumberOfVertices();
     1287        IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
     1288        IssmDouble* thickness = xNew<IssmDouble>(numvertices);
     1289        IssmDouble* surface   = xNew<IssmDouble>(numvertices);
     1290
     1291        element->FindParam(&meshtype,MeshTypeEnum);
     1292        rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
     1293        g       =element->GetMaterialParameter(ConstantsGEnum);
     1294        switch(meshtype){
     1295                case Mesh2DhorizontalEnum:
     1296                        element->GetInputListOnVertices(thickness,ThicknessEnum);
     1297                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
     1298                        break;
     1299                case Mesh3DEnum:
     1300                        element->GetVerticesCoordinates(&xyz_list);
     1301                        element->GetInputListOnVertices(surface,SurfaceEnum);
     1302                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
     1303                        break;
     1304                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1305        }
     1306        element->AddInput(PressureEnum,pressure,P1Enum);
     1307        xDelete<IssmDouble>(pressure);
     1308        xDelete<IssmDouble>(thickness);
     1309        xDelete<IssmDouble>(surface);
     1310
     1311        /*Get basal element*/
     1312        switch(meshtype){
     1313                case Mesh2DhorizontalEnum:
     1314                        basalelement = element;
     1315                        break;
     1316                case Mesh3DEnum:
     1317                        if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
     1318                        basalelement=element->SpawnBasalElement();
     1319                        break;
     1320                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1321        }
     1322
     1323        /*Fetch number of nodes and dof for this finite element*/
     1324        int numnodes = basalelement->GetNumberOfNodes();
     1325        int numdof   = numnodes*2;
     1326
     1327        /*Fetch dof list and allocate solution vectors*/
     1328        basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     1329        IssmDouble* values    = xNew<IssmDouble>(numdof);
     1330        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
     1331        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
     1332        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
     1333        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
     1334
     1335        /*Use the dof list to index into the solution vector: */
     1336        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
     1337
     1338        /*Transform solution in Cartesian Space*/
     1339        basalelement->TransformSolutionCoord(&values[0],XYEnum);
     1340        basalelement->FindParam(&meshtype,MeshTypeEnum);
     1341
     1342        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     1343        for(i=0;i<numnodes;i++){
     1344                vx[i]=values[i*2+0];
     1345                vy[i]=values[i*2+1];
     1346
     1347                /*Check solution*/
     1348                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
     1349                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     1350        }
     1351
     1352        /*Get Vz and compute vel*/
     1353        basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
     1354        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     1355
     1356        /*Now, we have to move the previous Vx and Vy inputs  to old
     1357         * status, otherwise, we'll wipe them off: */
     1358        element->InputChangeName(VxEnum,VxPicardEnum);
     1359        element->InputChangeName(VyEnum,VyPicardEnum);
     1360
     1361        /*Add vx and vy as inputs to the tria element: */
     1362        element->AddBasalInput(VxEnum,vx,P1Enum);
     1363        element->AddBasalInput(VyEnum,vy,P1Enum);
     1364        element->AddBasalInput(VelEnum,vel,P1Enum);
     1365
     1366        /*Free ressources:*/
     1367        xDelete<IssmDouble>(vel);
     1368        xDelete<IssmDouble>(vz);
     1369        xDelete<IssmDouble>(vy);
     1370        xDelete<IssmDouble>(vx);
     1371        xDelete<IssmDouble>(values);
     1372        xDelete<IssmDouble>(xyz_list);
     1373        xDelete<int>(doflist);
     1374        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     1375}/*}}}*/
     1376
     1377/*HO*/
     1378ElementMatrix* StressbalanceAnalysis::CreateKMatrixHO(Element* element){/*{{{*/
     1379
     1380        /*compute all stiffness matrices for this element*/
     1381        ElementMatrix* Ke1=CreateKMatrixHOViscous(element);
     1382        ElementMatrix* Ke2=CreateKMatrixHOFriction(element);
     1383        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     1384
     1385        /*clean-up and return*/
     1386        delete Ke1;
     1387        delete Ke2;
     1388        return Ke;
     1389}/*}}}*/
     1390ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOViscous(Element* element){/*{{{*/
     1391
     1392        /*Intermediaries*/
     1393        IssmDouble  viscosity,newviscosity,oldviscosity;
     1394        IssmDouble  viscosity_overshoot,thickness,Jdet;
     1395        IssmDouble  D_scalar;
     1396        IssmDouble *xyz_list = NULL;
     1397
     1398        /*Fetch number of nodes and dof for this finite element*/
     1399        int numnodes = element->GetNumberOfNodes();
     1400        int numdof   = numnodes*2;
     1401
     1402        /*Initialize Element matrix and vectors*/
     1403        ElementMatrix* Ke     = element->NewElementMatrix(HOApproximationEnum);
     1404        IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
     1405        IssmDouble*    Bprime = xNew<IssmDouble>(5*numdof);
     1406        IssmDouble*    D      = xNewZeroInit<IssmDouble>(5*5);
     1407
     1408        /*Retrieve all inputs and parameters*/
     1409        element->GetVerticesCoordinates(&xyz_list);
     1410        Input* vx_input=element->GetInput(VxEnum);               _assert_(vx_input);
     1411        Input* vy_input=element->GetInput(VyEnum);               _assert_(vy_input);
     1412        Input* vxold_input=element->GetInput(VxPicardEnum);      _assert_(vxold_input);
     1413        Input* vyold_input=element->GetInput(VyPicardEnum);      _assert_(vyold_input);
     1414        element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
     1415
     1416        /* Start  looping on the number of gaussian points: */
     1417        Gauss* gauss = element->NewGauss(5);
     1418        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1419                gauss->GaussPoint(ig);
     1420
     1421                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1422                this->GetBHO(B,element,xyz_list,gauss);
     1423                this->GetBHOprime(Bprime,element,xyz_list,gauss);
     1424
     1425                element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input);
     1426                element->ViscosityHO(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
     1427
     1428                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
     1429                D_scalar=2.*newviscosity*gauss->weight*Jdet;
     1430                for(int i=0;i<5;i++) D[i*5+i]=D_scalar;
     1431
     1432                TripleMultiply(B,5,numdof,1,
     1433                                        D,5,5,0,
     1434                                        Bprime,5,numdof,0,
     1435                                        &Ke->values[0],1);
     1436        }
     1437
     1438        /*Transform Coordinate System*/
     1439        element->TransformStiffnessMatrixCoord(Ke,XYEnum);
     1440
     1441        /*Clean up and return*/
     1442        delete gauss;
     1443        xDelete<IssmDouble>(xyz_list);
     1444        xDelete<IssmDouble>(D);
     1445        xDelete<IssmDouble>(Bprime);
     1446        xDelete<IssmDouble>(B);
     1447        return Ke;
     1448}/*}}}*/
     1449ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFriction(Element* element){/*{{{*/
     1450
     1451        if(element->IsFloating() || !element->IsOnBed()) return NULL;
     1452
     1453        return NULL;
     1454}/*}}}*/
     1455ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
     1456
     1457        /*compute all load vectors for this element*/
     1458        ElementVector* pe1=CreatePVectorHODrivingStress(element);
     1459        ElementVector* pe2=CreatePVectorHOFront(element);
     1460        ElementVector* pe =new ElementVector(pe1,pe2);
     1461
     1462        /*clean-up and return*/
     1463        delete pe1;
     1464        delete pe2;
     1465        return pe;
     1466}/*}}}*/
     1467ElementVector* StressbalanceAnalysis::CreatePVectorHODrivingStress(Element* element){/*{{{*/
     1468
     1469        /*Intermediaries */
     1470        IssmDouble  Jdet,slope[3];
     1471        IssmDouble* xyz_list = NULL;
     1472
     1473        /*Fetch number of nodes and dof for this finite element*/
     1474        int numnodes = element->GetNumberOfNodes();
     1475
     1476        /*Initialize Element vector and vectors*/
     1477        ElementVector* pe=element->NewElementVector(HOApproximationEnum);
     1478        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     1479
     1480        /*Retrieve all inputs and parameters*/
     1481        element->GetVerticesCoordinates(&xyz_list);
     1482        Input*     surface_input = element->GetInput(SurfaceEnum);   _assert_(surface_input);
     1483        IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
     1484
     1485        /* Start  looping on the number of gaussian points: */
     1486        Gauss* gauss=element->NewGauss(3);
     1487        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1488                gauss->GaussPoint(ig);
     1489
     1490                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1491                element->NodalFunctions(basis, gauss);
     1492                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
     1493
     1494                for(int i=0;i<numnodes;i++){
     1495                        pe->values[i*2+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i];
     1496                        pe->values[i*2+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i];
     1497                }
     1498        }
     1499
     1500        /*Transform coordinate system*/
     1501        element->TransformLoadVectorCoord(pe,XYEnum);
     1502
     1503        /*Clean up and return*/
     1504        xDelete<IssmDouble>(basis);
     1505        xDelete<IssmDouble>(xyz_list);
     1506        delete gauss;
     1507        return pe;
     1508}/*}}}*/
     1509ElementVector* StressbalanceAnalysis::CreatePVectorHOFront(Element* element){/*{{{*/
     1510
     1511        return NULL;
     1512}/*}}}*/
     1513void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1514        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
     1515         * For node i, Bi can be expressed in the actual coordinate system
     1516         * by:
     1517         *       Bi=[ dh/dx          0      ]
     1518         *          [   0           dh/dy   ]
     1519         *          [ 1/2*dh/dy  1/2*dh/dx  ]
     1520         *          [ 1/2*dh/dz      0      ]
     1521         *          [  0         1/2*dh/dz  ]
     1522         * where h is the interpolation function for node i.
     1523         *
     1524         * We assume B has been allocated already, of size: 5x(NDOF2*numnodes)
     1525         */
     1526
     1527        /*Fetch number of nodes for this finite element*/
     1528        int numnodes = element->GetNumberOfNodes();
     1529
     1530        /*Get nodal functions derivatives*/
     1531        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
     1532        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1533
     1534        /*Build B: */
     1535        for(int i=0;i<numnodes;i++){
     1536                B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
     1537                B[2*numnodes*0+2*i+1] = 0.;
     1538                B[2*numnodes*1+2*i+0] = 0.;
     1539                B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
     1540                B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
     1541                B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
     1542                B[2*numnodes*3+2*i+0] = .5*dbasis[2*numnodes+i];
     1543                B[2*numnodes*3+2*i+1] = 0.;
     1544                B[2*numnodes*4+2*i+0] = 0.;
     1545                B[2*numnodes*4+2*i+1] = .5*dbasis[2*numnodes+i];
     1546        }
     1547
     1548        /*Clean-up*/
     1549        xDelete<IssmDouble>(dbasis);
     1550}/*}}}*/
     1551void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1552        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
     1553         * For node i, Bi' can be expressed in the actual coordinate system
     1554         * by:
     1555         *       Bi_prime=[ 2*dN/dx    dN/dy ]
     1556         *                [   dN/dx  2*dN/dy ]
     1557         *                [   dN/dy    dN/dx ]
     1558         * where hNis the finiteelement function for node i.
     1559         *
     1560         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
     1561         */
     1562
     1563        /*Fetch number of nodes for this finite element*/
     1564        int numnodes = element->GetNumberOfNodes();
     1565
     1566        /*Get nodal functions derivatives*/
     1567        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
     1568        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1569
     1570        /*Build B': */
     1571        for(int i=0;i<numnodes;i++){
     1572                Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
     1573                Bprime[2*numnodes*0+2*i+1] = dbasis[1*numnodes+i];
     1574                Bprime[2*numnodes*1+2*i+0] = dbasis[0*numnodes+i];
     1575                Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
     1576                Bprime[2*numnodes*2+2*i+0] = dbasis[1*numnodes+i];
     1577                Bprime[2*numnodes*2+2*i+1] = dbasis[0*numnodes+i];
     1578                Bprime[2*numnodes*3+2*i+0] = dbasis[2*numnodes+i];
     1579                Bprime[2*numnodes*3+2*i+1] = 0.;
     1580                Bprime[2*numnodes*4+2*i+0] = 0.;
     1581                Bprime[2*numnodes*4+2*i+1] = dbasis[2*numnodes+i];
     1582        }
     1583
     1584        /*Clean-up*/
     1585        xDelete<IssmDouble>(dbasis);
    15491586}/*}}}*/
    15501587void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
    … …  
    16191656        xDelete<int>(doflist);
    16201657}/*}}}*/
     1658
     1659/*FS*/
     1660ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
     1661
     1662        /*compute all load vectors for this element*/
     1663        ElementVector* pe1=CreatePVectorFSViscous(element);
     1664        ElementVector* pe2=CreatePVectorFSShelf(element);
     1665        ElementVector* pe3=CreatePVectorFSFront(element);
     1666        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     1667
     1668        /*clean-up and return*/
     1669        delete pe1;
     1670        delete pe2;
     1671        delete pe3;
     1672        return pe;
     1673}/*}}}*/
     1674ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
     1675
     1676        int         i,meshtype,dim;
     1677        IssmDouble  Jdet,forcex,forcey,forcez;
     1678        IssmDouble *xyz_list = NULL;
     1679
     1680        /*Get problem dimension*/
     1681        element->FindParam(&meshtype,MeshTypeEnum);
     1682        switch(meshtype){
     1683                case Mesh2DverticalEnum: dim = 2; break;
     1684                case Mesh3DEnum:         dim = 3; break;
     1685                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1686        }
     1687
     1688        /*Fetch number of nodes and dof for this finite element*/
     1689        int vnumnodes = element->GetNumberOfNodesVelocity();
     1690        int pnumnodes = element->GetNumberOfNodesPressure();
     1691
     1692        /*Prepare coordinate system list*/
     1693        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     1694        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     1695        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     1696        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     1697
     1698        /*Initialize vectors*/
     1699        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
     1700        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     1701
     1702        /*Retrieve all inputs and parameters*/
     1703        element->GetVerticesCoordinates(&xyz_list);
     1704        Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
     1705        Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
     1706        Input*      loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
     1707        IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
     1708        IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
     1709
     1710        /* Start  looping on the number of gaussian points: */
     1711        Gauss* gauss=element->NewGauss(5);
     1712        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1713                gauss->GaussPoint(ig);
     1714
     1715                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1716                element->NodalFunctionsVelocity(vbasis,gauss);
     1717
     1718                loadingforcex_input->GetInputValue(&forcex,gauss);
     1719                loadingforcey_input->GetInputValue(&forcey,gauss);
     1720                if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss);
     1721
     1722                for(i=0;i<vnumnodes;i++){
     1723                        pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
     1724                        pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
     1725                        if(dim==3){
     1726                                pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
     1727                                pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
     1728                        }
     1729                        else{
     1730                                pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
     1731                        }
     1732                }
     1733        }
     1734
     1735        /*Transform coordinate system*/
     1736        element->TransformLoadVectorCoord(pe,cs_list);
     1737
     1738        /*Clean up and return*/
     1739        delete gauss;
     1740        xDelete<int>(cs_list);
     1741        xDelete<IssmDouble>(vbasis);
     1742        xDelete<IssmDouble>(xyz_list);
     1743        return pe;
     1744}/*}}}*/
     1745ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/
     1746
     1747        int         i,meshtype,dim;
     1748        IssmDouble  Jdet,water_pressure,bed;
     1749        IssmDouble      normal[3];
     1750        IssmDouble *xyz_list_base = NULL;
     1751
     1752        /*Get basal element*/
     1753        if(!element->IsOnBed() || !element->IsFloating()) return NULL;
     1754
     1755        /*Get problem dimension*/
     1756        element->FindParam(&meshtype,MeshTypeEnum);
     1757        switch(meshtype){
     1758                case Mesh2DverticalEnum: dim = 2; break;
     1759                case Mesh3DEnum:         dim = 3; break;
     1760                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1761        }
     1762
     1763        /*Fetch number of nodes and dof for this finite element*/
     1764        int vnumnodes = element->GetNumberOfNodesVelocity();
     1765        int pnumnodes = element->GetNumberOfNodesPressure();
     1766
     1767        /*Prepare coordinate system list*/
     1768        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     1769        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     1770        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     1771        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     1772
     1773        /*Initialize vectors*/
     1774        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
     1775        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     1776
     1777        /*Retrieve all inputs and parameters*/
     1778        element->GetVerticesCoordinatesBase(&xyz_list_base);
     1779        Input*      bed_input=element->GetInput(BedEnum); _assert_(bed_input);
     1780        IssmDouble  rho_water=element->GetMaterialParameter(MaterialsRhoWaterEnum);
     1781        IssmDouble  gravity  =element->GetMaterialParameter(ConstantsGEnum);
     1782
     1783        /* Start  looping on the number of gaussian points: */
     1784        Gauss* gauss=element->NewGaussBase(5);
     1785        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1786                gauss->GaussPoint(ig);
     1787
     1788                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     1789                element->NodalFunctionsVelocity(vbasis,gauss);
     1790
     1791                element->NormalBase(&normal[0],xyz_list_base);
     1792                _assert_(normal[dim-1]<0.);
     1793                bed_input->GetInputValue(&bed, gauss);
     1794                water_pressure=gravity*rho_water*bed;
     1795
     1796                for(i=0;i<vnumnodes;i++){
     1797                        pe->values[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0];
     1798                        pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1];
     1799                        if(dim==3){
     1800                                pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2];
     1801                        }
     1802                }
     1803        }
     1804
     1805        /*Transform coordinate system*/
     1806        element->TransformLoadVectorCoord(pe,cs_list);
     1807
     1808        /*Clean up and return*/
     1809        delete gauss;
     1810        xDelete<int>(cs_list);
     1811        xDelete<IssmDouble>(vbasis);
     1812        xDelete<IssmDouble>(xyz_list_base);
     1813        return pe;
     1814}/*}}}*/
     1815ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/
     1816
     1817        return NULL;
     1818}/*}}}*/
     1819void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
     1820
     1821        int*         vdoflist=NULL;
     1822        int*         pdoflist=NULL;
     1823        Input*       vz_input=NULL;
     1824        int          meshtype,dim;
     1825        IssmDouble   vx,vy,vz,p;
     1826        IssmDouble   FSreconditioning;
     1827
     1828        /*Get some parameters*/
     1829        element->FindParam(&meshtype,MeshTypeEnum);
     1830        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     1831        switch(meshtype){
     1832                case Mesh2DverticalEnum: dim = 2; break;
     1833                case Mesh3DEnum:         dim = 3; break;
     1834                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1835        }
     1836
     1837        /*Fetch number of nodes and dof for this finite element*/
     1838        int vnumnodes = element->NumberofNodesVelocity();
     1839        int pnumnodes = element->NumberofNodesPressure();
     1840        int vnumdof   = vnumnodes*dim;
     1841        int pnumdof   = pnumnodes*1;
     1842
     1843        /*Initialize values*/
     1844        IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
     1845        IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
     1846
     1847        /*Get dof list: */
     1848        element->GetDofListVelocity(&vdoflist,GsetEnum);
     1849        element->GetDofListPressure(&pdoflist,GsetEnum);
     1850        Input*     vx_input=element->GetInput(VxEnum);       _assert_(vx_input);
     1851        Input*     vy_input=element->GetInput(VyEnum);       _assert_(vy_input);
     1852        if(dim==3){vz_input=element->GetInput(VzEnum);       _assert_(vz_input);}
     1853        Input*     p_input =element->GetInput(PressureEnum); _assert_(p_input);
     1854
     1855        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     1856
     1857        /*Ok, we have the velocities in inputs, fill in solution */
     1858        Gauss* gauss = element->NewGauss();
     1859        for(int i=0;i<vnumnodes;i++){
     1860                gauss->GaussNode(element->VelocityInterpolation(),i);
     1861                vx_input->GetInputValue(&vx,gauss);
     1862                vy_input->GetInputValue(&vy,gauss);
     1863                vvalues[i*dim+0]=vx;
     1864                vvalues[i*dim+1]=vy;
     1865                if(dim==3){
     1866                        vz_input->GetInputValue(&vz,gauss);
     1867                        vvalues[i*dim+2]=vz;
     1868                }
     1869        }
     1870        for(int i=0;i<pnumnodes;i++){
     1871                gauss->GaussNode(element->PressureInterpolation(),i);
     1872                p_input->GetInputValue(&p ,gauss);
     1873                pvalues[i]=p/FSreconditioning;
     1874        }
     1875
     1876        /*Add value to global vector*/
     1877        solution->SetValues(vnumdof,vdoflist,vvalues,INS_VAL);
     1878        solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL);
     1879
     1880        /*Free ressources:*/
     1881        delete gauss;
     1882        xDelete<int>(pdoflist);
     1883        xDelete<int>(vdoflist);
     1884        xDelete<IssmDouble>(pvalues);
     1885        xDelete<IssmDouble>(vvalues);
     1886}/*}}}*/
     1887void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
     1888
     1889        int          i,dim,meshtype;
     1890        int*         vdoflist=NULL;
     1891        int*         pdoflist=NULL;
     1892        IssmDouble   FSreconditioning;
     1893
     1894        element->FindParam(&meshtype,MeshTypeEnum);
     1895        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     1896        switch(meshtype){
     1897                case Mesh2DverticalEnum: dim = 2; break;
     1898                case Mesh3DEnum:         dim = 3; break;
     1899                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1900        }
     1901
     1902        /*Fetch number of nodes and dof for this finite element*/
     1903        int vnumnodes = element->GetNumberOfNodesVelocity();
     1904        int pnumnodes = element->GetNumberOfNodesPressure();
     1905        int vnumdof   = vnumnodes*dim;
     1906        int pnumdof   = pnumnodes*1;
     1907
     1908        /*Initialize values*/
     1909        IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
     1910        IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
     1911        IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
     1912        IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
     1913        IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
     1914        IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
     1915
     1916        /*Prepare coordinate system list*/
     1917        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     1918        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     1919        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     1920        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     1921
     1922        /*Get dof list: */
     1923        element->GetDofListVelocity(&vdoflist,GsetEnum);
     1924        element->GetDofListPressure(&pdoflist,GsetEnum);
     1925
     1926        /*Use the dof list to index into the solution vector: */
     1927        for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
     1928        for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
     1929
     1930        /*Transform solution in Cartesian Space*/
     1931        element->TransformSolutionCoord(values,cs_list);
     1932
     1933        /*Ok, we have vx and vy in values, fill in all arrays: */
     1934        for(i=0;i<vnumnodes;i++){
     1935                vx[i] = values[i*dim+0];
     1936                vy[i] = values[i*dim+1];
     1937                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
     1938                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     1939
     1940                if(dim==3){
     1941                        vz[i] = values[i*dim+2];
     1942                        if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
     1943                }
     1944        }
     1945        for(i=0;i<pnumnodes;i++){
     1946                pressure[i] = values[vnumdof+i];
     1947                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
     1948        }
     1949
     1950        /*Recondition pressure and compute vel: */
     1951        for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
     1952        if(dim==3) for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     1953        else       for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
     1954
     1955        /*Now, we have to move the previous inputs  to old
     1956         * status, otherwise, we'll wipe them off: */
     1957        element->InputChangeName(VxEnum,VxPicardEnum);
     1958        element->InputChangeName(VyEnum,VyPicardEnum);
     1959        element->InputChangeName(PressureEnum,PressurePicardEnum);
     1960        if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);
     1961
     1962        /*Add vx and vy as inputs to the tria element: */
     1963        element->AddInput(VxEnum,vx,P1Enum);
     1964        element->AddInput(VyEnum,vy,P1Enum);
     1965        element->AddInput(VelEnum,vel,P1Enum);
     1966        element->AddInput(PressureEnum,pressure,P1Enum);
     1967        if(dim==3) element->AddInput(VzEnum,vz,P1Enum);
     1968
     1969        /*Free ressources:*/
     1970        xDelete<IssmDouble>(pressure);
     1971        xDelete<IssmDouble>(vel);
     1972        xDelete<IssmDouble>(vz);
     1973        xDelete<IssmDouble>(vy);
     1974        xDelete<IssmDouble>(vx);
     1975        xDelete<IssmDouble>(values);
     1976        xDelete<int>(vdoflist);
     1977        xDelete<int>(pdoflist);
     1978        xDelete<int>(cs_list);
     1979}/*}}}*/
     1980
     1981/*Coupling (Tiling)*/
    16211982void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
    16221983
    … …  
    17162077        xDelete<int>(cs_list);
    17172078}/*}}}*/
    1718 void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
    1719 
    1720         int         i,meshtype;
    1721         IssmDouble  rho_ice,g;
    1722         int*        doflist=NULL;
    1723         IssmDouble* xyz_list=NULL;
    1724         Element*    basalelement=NULL;
    1725 
    1726         /*Deal with pressure first*/
    1727         int numvertices = element->GetNumberOfVertices();
    1728         IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
    1729         IssmDouble* thickness = xNew<IssmDouble>(numvertices);
    1730         IssmDouble* surface   = xNew<IssmDouble>(numvertices);
    1731 
    1732         element->FindParam(&meshtype,MeshTypeEnum);
    1733         rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
    1734         g       =element->GetMaterialParameter(ConstantsGEnum);
    1735         switch(meshtype){
    1736                 case Mesh2DhorizontalEnum:
    1737                         element->GetInputListOnVertices(thickness,ThicknessEnum);
    1738                         for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
    1739                         break;
    1740                 case Mesh3DEnum:
    1741                         element->GetVerticesCoordinates(&xyz_list);
    1742                         element->GetInputListOnVertices(surface,SurfaceEnum);
    1743                         for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
    1744                         break;
    1745                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    1746         }
    1747         element->AddInput(PressureEnum,pressure,P1Enum);
    1748         xDelete<IssmDouble>(pressure);
    1749         xDelete<IssmDouble>(thickness);
    1750         xDelete<IssmDouble>(surface);
    1751 
    1752         /*Get basal element*/
    1753         switch(meshtype){
    1754                 case Mesh2DhorizontalEnum:
    1755                         basalelement = element;
    1756                         break;
    1757                 case Mesh3DEnum:
    1758                         if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
    1759                         basalelement=element->SpawnBasalElement();
    1760                         break;
    1761                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    1762         }
    1763 
    1764         /*Fetch number of nodes and dof for this finite element*/
    1765         int numnodes = basalelement->GetNumberOfNodes();
    1766         int numdof   = numnodes*2;
    1767 
    1768         /*Fetch dof list and allocate solution vectors*/
    1769         basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    1770         IssmDouble* values    = xNew<IssmDouble>(numdof);
    1771         IssmDouble* vx        = xNew<IssmDouble>(numnodes);
    1772         IssmDouble* vy        = xNew<IssmDouble>(numnodes);
    1773         IssmDouble* vz        = xNew<IssmDouble>(numnodes);
    1774         IssmDouble* vel       = xNew<IssmDouble>(numnodes);
    1775 
    1776         /*Use the dof list to index into the solution vector: */
    1777         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    1778 
    1779         /*Transform solution in Cartesian Space*/
    1780         basalelement->TransformSolutionCoord(&values[0],XYEnum);
    1781         basalelement->FindParam(&meshtype,MeshTypeEnum);
    1782 
    1783         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    1784         for(i=0;i<numnodes;i++){
    1785                 vx[i]=values[i*2+0];
    1786                 vy[i]=values[i*2+1];
    1787 
    1788                 /*Check solution*/
    1789                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    1790                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    1791         }
    1792 
    1793         /*Get Vz and compute vel*/
    1794         basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
    1795         for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    1796 
    1797         /*Now, we have to move the previous Vx and Vy inputs  to old
    1798          * status, otherwise, we'll wipe them off: */
    1799         element->InputChangeName(VxEnum,VxPicardEnum);
    1800         element->InputChangeName(VyEnum,VyPicardEnum);
    1801 
    1802         /*Add vx and vy as inputs to the tria element: */
    1803         element->AddBasalInput(VxEnum,vx,P1Enum);
    1804         element->AddBasalInput(VyEnum,vy,P1Enum);
    1805         element->AddBasalInput(VelEnum,vel,P1Enum);
    1806 
    1807         /*Free ressources:*/
    1808         xDelete<IssmDouble>(vel);
    1809         xDelete<IssmDouble>(vz);
    1810         xDelete<IssmDouble>(vy);
    1811         xDelete<IssmDouble>(vx);
    1812         xDelete<IssmDouble>(values);
    1813         xDelete<IssmDouble>(xyz_list);
    1814         xDelete<int>(doflist);
    1815         if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    1816 }/*}}}*/
    1817 void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
    1818 
    1819         int         i,meshtype;
    1820         IssmDouble  rho_ice,g;
    1821         int*        doflist=NULL;
    1822         IssmDouble* xyz_list=NULL;
    1823         Element*    basalelement=NULL;
    1824 
    1825         /*Deal with pressure first*/
    1826         int numvertices = element->GetNumberOfVertices();
    1827         IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
    1828         IssmDouble* thickness = xNew<IssmDouble>(numvertices);
    1829         IssmDouble* surface   = xNew<IssmDouble>(numvertices);
    1830 
    1831         element->FindParam(&meshtype,MeshTypeEnum);
    1832         rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
    1833         g       =element->GetMaterialParameter(ConstantsGEnum);
    1834         switch(meshtype){
    1835                 case Mesh2DhorizontalEnum:
    1836                         element->GetInputListOnVertices(thickness,ThicknessEnum);
    1837                         for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
    1838                         break;
    1839                 case Mesh3DEnum:
    1840                         element->GetVerticesCoordinates(&xyz_list);
    1841                         element->GetInputListOnVertices(surface,SurfaceEnum);
    1842                         for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
    1843                         break;
    1844                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    1845         }
    1846         element->AddInput(PressureEnum,pressure,P1Enum);
    1847         xDelete<IssmDouble>(pressure);
    1848         xDelete<IssmDouble>(thickness);
    1849         xDelete<IssmDouble>(surface);
    1850 
    1851         /*Get basal element*/
    1852         switch(meshtype){
    1853                 case Mesh2DhorizontalEnum:
    1854                         basalelement = element;
    1855                         break;
    1856                 case Mesh3DEnum:
    1857                         if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
    1858                         basalelement=element->SpawnBasalElement();
    1859                         break;
    1860                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    1861         }
    1862 
    1863         /*Fetch number of nodes and dof for this finite element*/
    1864         int numnodes = basalelement->GetNumberOfNodes();
    1865         int numdof   = numnodes*2;
    1866 
    1867         /*Fetch dof list and allocate solution vectors*/
    1868         basalelement->GetDofList(&doflist,SSAApproximationEnum,GsetEnum);
    1869         IssmDouble* values    = xNew<IssmDouble>(numdof);
    1870         IssmDouble* vx        = xNew<IssmDouble>(numnodes);
    1871         IssmDouble* vy        = xNew<IssmDouble>(numnodes);
    1872         IssmDouble* vz        = xNew<IssmDouble>(numnodes);
    1873         IssmDouble* vel       = xNew<IssmDouble>(numnodes);
    1874 
    1875         /*Use the dof list to index into the solution vector: */
    1876         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    1877 
    1878         /*Transform solution in Cartesian Space*/
    1879         basalelement->TransformSolutionCoord(&values[0],XYEnum);
    1880         basalelement->FindParam(&meshtype,MeshTypeEnum);
    1881 
    1882         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    1883         for(i=0;i<numnodes;i++){
    1884                 vx[i]=values[i*2+0];
    1885                 vy[i]=values[i*2+1];
    1886 
    1887                 /*Check solution*/
    1888                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    1889                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    1890         }
    1891 
    1892         /*Get Vz and compute vel*/
    1893         basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
    1894         for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    1895 
    1896         /*Now, we have to move the previous Vx and Vy inputs  to old
    1897          * status, otherwise, we'll wipe them off: */
    1898         element->InputChangeName(VxEnum,VxPicardEnum);
    1899         element->InputChangeName(VyEnum,VyPicardEnum);
    1900 
    1901         /*Add vx and vy as inputs to the tria element: */
    1902         element->AddBasalInput(VxEnum,vx,P1Enum);
    1903         element->AddBasalInput(VyEnum,vy,P1Enum);
    1904         element->AddBasalInput(VelEnum,vel,P1Enum);
    1905 
    1906         /*Free ressources:*/
    1907         xDelete<IssmDouble>(vel);
    1908         xDelete<IssmDouble>(vz);
    1909         xDelete<IssmDouble>(vy);
    1910         xDelete<IssmDouble>(vx);
    1911         xDelete<IssmDouble>(values);
    1912         xDelete<IssmDouble>(xyz_list);
    1913         xDelete<int>(doflist);
    1914         if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    1915 }/*}}}*/
    19162079void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/
    19172080

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

    r16818 r16829  
    2222                /*Finite element Analysis*/
    2323                ElementMatrix* CreateKMatrix(Element* element);
     24                ElementMatrix* CreateKMatrixHO(Element* element);
     25                ElementMatrix* CreateKMatrixHOViscous(Element* element);
     26                ElementMatrix* CreateKMatrixHOFriction(Element* element);
    2427                ElementMatrix* CreateKMatrixSSA(Element* element);
    2528                ElementMatrix* CreateKMatrixSSAViscous(Element* element);
    … …  
    3639                ElementVector* CreatePVectorSSADrivingStress(Element* element);
    3740                ElementVector* CreatePVectorSSAFront(Element* element);
     41                void GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     42                void GetBHOprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    3843                void GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    3944                void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);

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

    r16827 r16829  
    144144                virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0;
    145145                virtual void   ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
     146                virtual void   ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
    146147                virtual void   ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
    147148                virtual int    VelocityInterpolation()=0;

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

    r16827 r16829  
    15431543/*}}}*/
    15441544/*FUNCTION Penta::GetStrainRate3dHO{{{*/
    1545 void Penta::GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){
     1545void Penta::GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input){
    15461546        /*Compute the 3d Blatter/HOStrain Rate (5 components):
    15471547         *
    … …  
    15721572/*}}}*/
    15731573/*FUNCTION Penta::GetStrainRate3d{{{*/
    1574 void Penta::GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
     1574void Penta::GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
    15751575        /*Compute the 3d Strain Rate (6 components):
    15761576         *
    … …  
    34823482        /*Assign output pointer*/
    34833483        *pphi = phi;
     3484}
     3485/*}}}*/
     3486/*FUNCTION Penta::ViscosityHO{{{*/
     3487void Penta::ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){
     3488
     3489        /*Intermediaries*/
     3490        IssmDouble viscosity;
     3491        IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
     3492
     3493        this->GetStrainRate3dHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     3494        material->GetViscosity3d(&viscosity, &epsilon[0]);
     3495
     3496        /*Assign output pointer*/
     3497        *pviscosity=viscosity;
    34843498}
    34853499/*}}}*/

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

    r16827 r16829  
    230230                void             GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
    231231                void           GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
    232                 void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
    233                 void           GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
     232                void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input);
     233                void           GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
    234234                Penta*         GetUpperElement(void);
    235235                void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum);
    … …  
    275275                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum){_error_("not implemented yet");};
    276276                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};
     277                void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    277278                void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    278279

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

    r16827 r16829  
    145145                ElementMatrix* NewElementMatrix(int approximation_enum){_error_("not implemented yet");};
    146146                void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
     147                void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    147148                void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    148149                #ifdef _HAVE_THERMAL_

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

    r16827 r16829  
    305305                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};
    306306                void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
     307                void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented yet");};
    307308                void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    308309
Note: See TracChangeset for help on using the changeset viewer.