Changeset 16993


Ignore:
Timestamp:
12/03/13 10:36:26 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: moving all Stiffness matrices and Jacobian matrices to analyses

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

Legend:

Unmodified
Added
Removed

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

    r16992 r16993  
    813813/*Finite Element Analysis*/
    814814ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
    815 _error_("Not implemented");
     815
     816        int approximation;
     817        element->GetInputValue(&approximation,ApproximationEnum);
     818        switch(approximation){
     819                case SSAApproximationEnum:
     820                        return CreateJacobianMatrixSSA(element);
     821                case HOApproximationEnum:
     822                        return CreateJacobianMatrixHO(element);
     823                case FSApproximationEnum:
     824                        return CreateJacobianMatrixFS(element);
     825                case NoneApproximationEnum:
     826                        return NULL;
     827                default:
     828                        _error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
     829        }
    816830}/*}}}*/
    817831ElementMatrix* StressbalanceAnalysis::CreateKMatrix(Element* element){/*{{{*/
    … …  
    961975
    962976/*SSA*/
     977ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrixSSA(Element* element){/*{{{*/
     978
     979        /*Intermediaries*/
     980        int      meshtype;
     981        Element* basalelement;
     982
     983        /*Get basal element*/
     984        element->FindParam(&meshtype,MeshTypeEnum);
     985        switch(meshtype){
     986                case Mesh2DhorizontalEnum:
     987                        basalelement = element;
     988                        break;
     989                case Mesh3DEnum:
     990                        if(!element->IsOnBed()) return NULL;
     991                        basalelement = element->SpawnBasalElement();
     992                        break;
     993                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     994        }
     995
     996        /*Intermediaries */
     997        IssmDouble Jdet,thickness;
     998        IssmDouble eps1dotdphii,eps1dotdphij;
     999        IssmDouble eps2dotdphii,eps2dotdphij;
     1000        IssmDouble mu_prime;
     1001        IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/
     1002        IssmDouble eps1[2],eps2[2];
     1003        IssmDouble *xyz_list = NULL;
     1004
     1005        /*Fetch number of nodes and dof for this finite element*/
     1006        int numnodes = basalelement->GetNumberOfNodes();
     1007
     1008        /*Initialize Element matrix, vectors and Gaussian points*/
     1009        ElementMatrix* Ke=this->CreateKMatrixSSA(element); //Initialize Jacobian with regular SSA (first part of the Gateau derivative)
     1010        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
     1011
     1012        /*Retrieve all inputs and parameters*/
     1013        element->GetVerticesCoordinates(&xyz_list);
     1014        Input* thickness_input = basalelement->GetInput(ThicknessEnum);_assert_(thickness_input);
     1015        Input* vx_input        = basalelement->GetInput(VxEnum);       _assert_(vx_input);
     1016        Input* vy_input        = basalelement->GetInput(VyEnum);       _assert_(vy_input);
     1017
     1018        /* Start  looping on the number of gaussian points: */
     1019        Gauss* gauss = basalelement->NewGauss(2);
     1020        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1021                gauss->GaussPoint(ig);
     1022
     1023                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1024                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1025
     1026                thickness_input->GetInputValue(&thickness, gauss);
     1027                basalelement->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     1028                basalelement->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]);
     1029                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
     1030                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
     1031
     1032                for(int i=0;i<numnodes;i++){
     1033                        for(int j=0;j<numnodes;j++){
     1034                                eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
     1035                                eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
     1036                                eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
     1037                                eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];
     1038
     1039                                Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
     1040                                Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
     1041                                Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
     1042                                Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
     1043                        }
     1044                }
     1045        }
     1046
     1047        /*Transform Coordinate System*/
     1048        basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum);
     1049
     1050        /*Clean up and return*/
     1051        delete gauss;
     1052        xDelete<IssmDouble>(xyz_list);
     1053        xDelete<IssmDouble>(dbasis);
     1054
     1055        /*clean-up and return*/
     1056        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     1057        return Ke;
     1058
     1059}/*}}}*/
    9631060ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA(Element* element){/*{{{*/
    9641061
    … …  
    17781875
    17791876/*HO*/
     1877ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrixHO(Element* element){/*{{{*/
     1878
     1879        /*Intermediaries */
     1880        IssmDouble Jdet;
     1881        IssmDouble eps1dotdphii,eps1dotdphij;
     1882        IssmDouble eps2dotdphii,eps2dotdphij;
     1883        IssmDouble mu_prime;
     1884        IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     1885        IssmDouble eps1[3],eps2[3];
     1886        IssmDouble *xyz_list = NULL;
     1887
     1888        /*Fetch number of nodes and dof for this finite element*/
     1889        int numnodes = element->GetNumberOfNodes();
     1890
     1891        /*Initialize Element matrix, vectors and Gaussian points*/
     1892        ElementMatrix* Ke=this->CreateKMatrixHO(element); //Initialize Jacobian with regular HO (first part of the Gateau derivative)
     1893        IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
     1894
     1895        /*Retrieve all inputs and parameters*/
     1896        element->GetVerticesCoordinates(&xyz_list);
     1897        Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
     1898        Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
     1899
     1900        /* Start  looping on the number of gaussian points: */
     1901        Gauss* gauss = element->NewGauss(5);
     1902        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1903                gauss->GaussPoint(ig);
     1904
     1905                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1906                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     1907
     1908                element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     1909                element->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]);
     1910                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
     1911                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
     1912                eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
     1913
     1914                for(int i=0;i<numnodes;i++){
     1915                        for(int j=0;j<numnodes;j++){
     1916                                eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i]+eps1[2]*dbasis[2*numnodes+i];
     1917                                eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j]+eps1[2]*dbasis[2*numnodes+j];
     1918                                eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i]+eps2[2]*dbasis[2*numnodes+i];
     1919                                eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j]+eps2[2]*dbasis[2*numnodes+j];
     1920
     1921                                Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*eps1dotdphij*eps1dotdphii;
     1922                                Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*eps2dotdphij*eps1dotdphii;
     1923                                Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2.*mu_prime*eps1dotdphij*eps2dotdphii;
     1924                                Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2.*mu_prime*eps2dotdphij*eps2dotdphii;
     1925                        }
     1926                }
     1927        }
     1928
     1929        /*Transform Coordinate System*/
     1930        element->TransformStiffnessMatrixCoord(Ke,XYEnum);
     1931
     1932        /*Clean up and return*/
     1933        delete gauss;
     1934        xDelete<IssmDouble>(xyz_list);
     1935        xDelete<IssmDouble>(dbasis);
     1936        return Ke;
     1937}/*}}}*/
    17801938ElementMatrix* StressbalanceAnalysis::CreateKMatrixHO(Element* element){/*{{{*/
    17811939
    … …  
    22192377
    22202378/*FS*/
     2379ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrixFS(Element* element){/*{{{*/
     2380
     2381        /*Intermediaries */
     2382        int        i,j;
     2383        IssmDouble Jdet;
     2384        IssmDouble eps1dotdphii,eps1dotdphij;
     2385        IssmDouble eps2dotdphii,eps2dotdphij;
     2386        IssmDouble eps3dotdphii,eps3dotdphij;
     2387        IssmDouble mu_prime;
     2388        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     2389        IssmDouble eps1[3],eps2[3],eps3[3];
     2390        IssmDouble *xyz_list = NULL;
     2391
     2392        /*Fetch number of nodes and dof for this finite element*/
     2393        int vnumnodes = element->GetNumberOfNodesVelocity();
     2394        int pnumnodes = element->GetNumberOfNodesPressure();
     2395        int numdof    = vnumnodes*NDOF3 + pnumnodes*NDOF1;
     2396
     2397        /*Prepare coordinate system list*/
     2398        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     2399        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     2400        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     2401
     2402        /*Initialize Element matrix, vectors and Gaussian points*/
     2403        ElementMatrix* Ke=this->CreateKMatrixFS(element); //Initialize Jacobian with regular FS (first part of the Gateau derivative)
     2404        IssmDouble*    dbasis = xNew<IssmDouble>(3*vnumnodes);
     2405
     2406        /*Retrieve all inputs and parameters*/
     2407        element->GetVerticesCoordinates(&xyz_list);
     2408        Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
     2409        Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
     2410        Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input);
     2411
     2412        /* Start  looping on the number of gaussian points: */
     2413        Gauss* gauss = element->NewGauss(5);
     2414        for(int ig=gauss->begin();ig<gauss->end();ig++){
     2415                gauss->GaussPoint(ig);
     2416
     2417                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     2418                element->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss);
     2419
     2420                //element->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     2421                //eps1[0]=epsilon[0];   eps2[0]=epsilon[3];   eps3[0]=epsilon[4];
     2422                //eps1[1]=epsilon[3];   eps2[1]=epsilon[1];   eps3[1]=epsilon[5];
     2423                //eps1[2]=epsilon[4];   eps2[2]=epsilon[5];   eps3[2]=epsilon[2];
     2424                element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     2425                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
     2426                eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
     2427                eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
     2428                element->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]);
     2429
     2430                for(i=0;i<vnumnodes;i++){
     2431                        for(j=0;j<vnumnodes;j++){
     2432                                eps1dotdphii=eps1[0]*dbasis[0*vnumnodes+i]+eps1[1]*dbasis[1*vnumnodes+i]+eps1[2]*dbasis[2*vnumnodes+i];
     2433                                eps1dotdphij=eps1[0]*dbasis[0*vnumnodes+j]+eps1[1]*dbasis[1*vnumnodes+j]+eps1[2]*dbasis[2*vnumnodes+j];
     2434                                eps2dotdphii=eps2[0]*dbasis[0*vnumnodes+i]+eps2[1]*dbasis[1*vnumnodes+i]+eps2[2]*dbasis[2*vnumnodes+i];
     2435                                eps2dotdphij=eps2[0]*dbasis[0*vnumnodes+j]+eps2[1]*dbasis[1*vnumnodes+j]+eps2[2]*dbasis[2*vnumnodes+j];
     2436                                eps3dotdphii=eps3[0]*dbasis[0*vnumnodes+i]+eps3[1]*dbasis[1*vnumnodes+i]+eps3[2]*dbasis[2*vnumnodes+i];
     2437                                eps3dotdphij=eps3[0]*dbasis[0*vnumnodes+j]+eps3[1]*dbasis[1*vnumnodes+j]+eps3[2]*dbasis[2*vnumnodes+j];
     2438
     2439                                Ke->values[numdof*(3*i+0)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
     2440                                Ke->values[numdof*(3*i+0)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
     2441                                Ke->values[numdof*(3*i+0)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
     2442
     2443                                Ke->values[numdof*(3*i+1)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
     2444                                Ke->values[numdof*(3*i+1)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
     2445                                Ke->values[numdof*(3*i+1)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
     2446
     2447                                Ke->values[numdof*(3*i+2)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
     2448                                Ke->values[numdof*(3*i+2)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
     2449                                Ke->values[numdof*(3*i+2)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
     2450                        }
     2451                }
     2452        }
     2453
     2454        /*Transform Coordinate System*/
     2455        element->TransformStiffnessMatrixCoord(Ke,cs_list);
     2456
     2457        /*Clean up and return*/
     2458        delete gauss;
     2459        xDelete<IssmDouble>(xyz_list);
     2460        xDelete<IssmDouble>(dbasis);
     2461        xDelete<int>(cs_list);
     2462        return Ke;
     2463}/*}}}*/
    22212464ElementMatrix* StressbalanceAnalysis::CreateKMatrixFS(Element* element){/*{{{*/
    22222465

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

    r16992 r16993  
    2929
    3030                /*SSA*/
     31                ElementMatrix* CreateJacobianMatrixSSA(Element* element);
    3132                ElementMatrix* CreateKMatrixSSA(Element* element);
    3233                ElementMatrix* CreateKMatrixSSAViscous(Element* element);
    … …  
    4849                void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
    4950                /*HO*/
     51                ElementMatrix* CreateJacobianMatrixHO(Element* element);
    5052                ElementMatrix* CreateKMatrixHO(Element* element);
    5153                ElementMatrix* CreateKMatrixHOViscous(Element* element);
    … …  
    5961                void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
    6062                /*FS*/
     63                ElementMatrix* CreateJacobianMatrixFS(Element* element);
    6164                ElementMatrix* CreateKMatrixFS(Element* element);
    6265                ElementMatrix* CreateKMatrixFSViscous(Element* element);

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

    r16982 r16993  
    104104                virtual void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0;
    105105                virtual void       SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
    106                 virtual ElementMatrix* CreateKMatrix(void)=0;
    107106                virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
    108                 virtual void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
    109107                virtual void   ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
    110108                virtual void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0;

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

    r16982 r16993  
    379379}
    380380/*}}}*/
    381 /*FUNCTION Penta::CreateKMatrix(void){{{*/
    382 ElementMatrix* Penta::CreateKMatrix(void){
    383 
    384         /*retrieve parameters: */
    385         int analysis_type;
    386         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    387 
    388         /*Checks in debugging {{{*/
    389         _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
    390         /*}}}*/
    391 
    392         /*Skip if water element*/
    393         if(NoIceInElement()) return NULL;
    394 
    395         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    396         switch(analysis_type){
    397                 #ifdef _HAVE_STRESSBALANCE_
    398                 case StressbalanceAnalysisEnum:
    399                         return CreateKMatrixStressbalanceHoriz();
    400                         break;
    401                 case AdjointHorizAnalysisEnum:
    402                         return CreateKMatrixAdjointHoriz();
    403                         break;
    404                 case StressbalanceSIAAnalysisEnum:
    405                         return CreateKMatrixStressbalanceSIA();
    406                         break;
    407                 case StressbalanceVerticalAnalysisEnum:
    408                         return CreateKMatrixStressbalanceVert();
    409                         break;
    410                 #endif
    411                 case L2ProjectionBaseAnalysisEnum:
    412                         return CreateBasalMassMatrix();
    413                         break;
    414                 case MasstransportAnalysisEnum:
    415                         return CreateKMatrixMasstransport();
    416                         break;
    417                 case FreeSurfaceTopAnalysisEnum:
    418                         return CreateKMatrixFreeSurfaceTop();
    419                         break;
    420                 case FreeSurfaceBaseAnalysisEnum:
    421                         return CreateKMatrixFreeSurfaceBase();
    422                         break;
    423                 #ifdef _HAVE_BALANCED_
    424                 case BalancethicknessAnalysisEnum:
    425                         return CreateKMatrixBalancethickness();
    426                         break;
    427                 #endif
    428                 #ifdef _HAVE_THERMAL_
    429                 case ThermalAnalysisEnum:
    430                         return CreateKMatrixThermal();
    431                         break;
    432                 case EnthalpyAnalysisEnum:
    433                         return CreateKMatrixEnthalpy();
    434                         break;
    435                 case MeltingAnalysisEnum:
    436                         return CreateKMatrixMelting();
    437                         break;
    438                 #endif
    439                 #ifdef _HAVE_HYDROLOGY_
    440           case HydrologyDCInefficientAnalysisEnum:
    441                         return CreateKMatrixHydrologyDCInefficient();
    442                         break;
    443           case HydrologyDCEfficientAnalysisEnum:
    444                         return CreateKMatrixHydrologyDCEfficient();
    445                         break;
    446           case L2ProjectionEPLAnalysisEnum:
    447                         return CreateEPLDomainMassMatrix();
    448                         break;
    449                 #endif
    450                 default:
    451                         _error_("analysis " << EnumToStringx(analysis_type) << " not supported yet");
    452         }
    453 
    454         /*Make compiler happy*/
    455         return NULL;
    456 }
    457 /*}}}*/
    458 /*FUNCTION Penta::CreateKMatrixMasstransport {{{*/
    459 ElementMatrix* Penta::CreateKMatrixMasstransport(void){
    460 
    461         if (!IsOnBed()) return NULL;
    462 
    463         /*Depth Averaging Vx and Vy*/
    464         this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
    465         this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
    466 
    467         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    468         ElementMatrix* Ke=tria->CreateKMatrixMasstransport();
    469         delete tria->material; delete tria;
    470 
    471         /*Delete Vx and Vy averaged*/
    472         this->inputs->DeleteInput(VxAverageEnum);
    473         this->inputs->DeleteInput(VyAverageEnum);
    474 
    475         /*clean up and return*/
    476         return Ke;
    477 }
    478 /*}}}*/
    479 /*FUNCTION Penta::CreateKMatrixFreeSurfaceTop {{{*/
    480 ElementMatrix* Penta::CreateKMatrixFreeSurfaceTop(void){
    481 
    482         if(!IsOnSurface()) return NULL;
    483 
    484         Tria* tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
    485         ElementMatrix* Ke=tria->CreateKMatrixFreeSurfaceTop();
    486         delete tria->material; delete tria;
    487 
    488         /*clean up and return*/
    489         return Ke;
    490 }
    491 /*}}}*/
    492 /*FUNCTION Penta::CreateKMatrixFreeSurfaceBase {{{*/
    493 ElementMatrix* Penta::CreateKMatrixFreeSurfaceBase(void){
    494 
    495         if(!IsOnBed()) return NULL;
    496 
    497         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    498         ElementMatrix* Ke=tria->CreateKMatrixFreeSurfaceBase();
    499         delete tria->material; delete tria;
    500 
    501         /*clean up and return*/
    502         return Ke;
    503 }
    504 /*}}}*/
    505 /*FUNCTION Penta::CreateBasalMassMatrix{{{*/
    506 ElementMatrix* Penta::CreateBasalMassMatrix(void){
    507 
    508         if (!IsOnBed()) return NULL;
    509 
    510         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    511         ElementMatrix* Ke=tria->CreateMassMatrix();
    512         delete tria->material; delete tria;
    513 
    514         /*clean up and return*/
    515         return Ke;
    516 }
    517 /*}}}*/
    518381/*FUNCTION Penta::CreateDVector {{{*/
    519382void  Penta::CreateDVector(Vector<IssmDouble>* df){
    … …  
    539402                De->InsertIntoGlobal(df);
    540403                delete De;
    541         }
    542 }
    543 /*}}}*/
    544 /*FUNCTION Penta::CreateJacobianMatrix{{{*/
    545 void  Penta::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){
    546 
    547         /*retrieve parameters: */
    548         ElementMatrix* Ke=NULL;
    549         int analysis_type;
    550         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    551 
    552         /*Checks in debugging {{{*/
    553         _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
    554         /*}}}*/
    555 
    556         /*Skip if water element*/
    557         if(NoIceInElement()) return;
    558 
    559         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    560         switch(analysis_type){
    561 #ifdef _HAVE_STRESSBALANCE_
    562                 case StressbalanceAnalysisEnum:
    563                         Ke=CreateJacobianStressbalanceHoriz();
    564                         break;
    565 #endif
    566                 default:
    567                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    568         }
    569 
    570         /*Add to global matrix*/
    571         if(Ke){
    572                 /*Condense if requested*/
    573                 if(this->element_type==MINIcondensedEnum){
    574                         int indices[3]={18,19,20};
    575                         Ke->StaticCondensation(3,&indices[0]);
    576                 }
    577                 else if(this->element_type==P1bubblecondensedEnum){
    578                         int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
    579                         int offset = 0;
    580                         for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
    581                         int* indices=xNew<int>(size);
    582                         for(int i=0;i<size;i++) indices[i] = offset+i;
    583                         Ke->StaticCondensation(size,indices);
    584                         xDelete<int>(indices);
    585                 }
    586 
    587                 Ke->AddToGlobal(Jff);
    588                 delete Ke;
    589404        }
    590405}
    … …  
    34373252
    34383253#ifdef _HAVE_THERMAL_
    3439 /*FUNCTION Penta::CreateKMatrixEnthalpy {{{*/
    3440 ElementMatrix* Penta::CreateKMatrixEnthalpy(void){
    3441 
    3442         /*compute all stiffness matrices for this element*/
    3443         ElementMatrix* Ke1=CreateKMatrixEnthalpyVolume();
    3444         ElementMatrix* Ke2=CreateKMatrixEnthalpyShelf();
    3445         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    3446 
    3447         /*clean-up and return*/
    3448         delete Ke1;
    3449         delete Ke2;
    3450         return Ke;
    3451 }
    3452 /*}}}*/
    3453 /*FUNCTION Penta::CreateKMatrixEnthalpyVolume {{{*/
    3454 ElementMatrix* Penta::CreateKMatrixEnthalpyVolume(void){
    3455 
    3456         /*Constants*/
    3457         const int    numdof=NDOF1*NUMVERTICES;
    3458 
    3459         /*Intermediaries */
    3460         int        stabilization;
    3461         int        i,j;
    3462         IssmDouble Jdet,u,v,w,um,vm,wm;
    3463         IssmDouble h,hx,hy,hz,vx,vy,vz,vel;
    3464         IssmDouble gravity,rho_ice,rho_water;
    3465         IssmDouble epsvel=2.220446049250313e-16;
    3466         IssmDouble heatcapacity,thermalconductivity,dt;
    3467         IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES];
    3468         IssmDouble latentheat,kappa;
    3469         IssmDouble tau_parameter,diameter;
    3470         IssmDouble xyz_list[NUMVERTICES][3];
    3471         IssmDouble B_conduct[3][numdof];
    3472         IssmDouble B_advec[3][numdof];
    3473         IssmDouble Bprime_advec[3][numdof];
    3474         IssmDouble L[numdof];
    3475         IssmDouble dbasis[3][6];
    3476         IssmDouble D_scalar_conduct,D_scalar_advec;
    3477         IssmDouble D_scalar_trans,D_scalar_stab;
    3478         IssmDouble D[3][3];
    3479         IssmDouble K[3][3]={0.0};
    3480         GaussPenta *gauss=NULL;
    3481 
    3482         /*Initialize Element matrix*/
    3483         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    3484 
    3485         /*Retrieve all inputs and parameters*/
    3486         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3487         rho_water=matpar->GetRhoWater();
    3488         rho_ice=matpar->GetRhoIce();
    3489         gravity=matpar->GetG();
    3490         heatcapacity=matpar->GetHeatCapacity();
    3491         latentheat=matpar->GetLatentHeat();
    3492         thermalconductivity=matpar->GetThermalConductivity();
    3493         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    3494         this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    3495         Input* vx_input=inputs->GetInput(VxEnum);                  _assert_(vx_input);
    3496         Input* vy_input=inputs->GetInput(VyEnum);                  _assert_(vy_input);
    3497         Input* vz_input=inputs->GetInput(VzEnum);                  _assert_(vz_input);
    3498         Input* vxm_input=inputs->GetInput(VxMeshEnum);             _assert_(vxm_input);
    3499         Input* vym_input=inputs->GetInput(VyMeshEnum);             _assert_(vym_input);
    3500         Input* vzm_input=inputs->GetInput(VzMeshEnum);             _assert_(vzm_input);
    3501         if (stabilization==2) diameter=MinEdgeLength(&xyz_list[0][0]);
    3502 
    3503         /* Start  looping on the number of gaussian points: */
    3504         gauss=new GaussPenta(2,2);
    3505         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3506 
    3507                 gauss->GaussPoint(ig);
    3508 
    3509                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3510 
    3511                 /*Conduction: */ 
    3512                 GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss);
    3513                 GetInputListOnVertices(&enthalpy[0],EnthalpyPicardEnum);
    3514                 GetInputListOnVertices(&pressure[0],PressureEnum);
    3515                 kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,enthalpy,pressure); _assert_(kappa>0.);
    3516 
    3517                 D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice;
    3518                 if(reCast<bool,IssmDouble>(dt)) D_scalar_conduct=D_scalar_conduct*dt;
    3519 
    3520                 D[0][0]=D_scalar_conduct; D[0][1]=0; D[0][2]=0;
    3521                 D[1][0]=0; D[1][1]=D_scalar_conduct; D[1][2]=0;
    3522                 D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_conduct;
    3523 
    3524                 TripleMultiply(&B_conduct[0][0],3,numdof,1,
    3525                                         &D[0][0],3,3,0,
    3526                                         &B_conduct[0][0],3,numdof,0,
    3527                                         &Ke->values[0],1);
    3528 
    3529                 /*Advection: */
    3530                 GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss);
    3531                 GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
    3532 
    3533                 vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
    3534                 vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
    3535                 vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
    3536 
    3537                 D_scalar_advec=gauss->weight*Jdet;
    3538                 if(reCast<bool,IssmDouble>(dt)) D_scalar_advec=D_scalar_advec*dt;
    3539 
    3540                 D[0][0]=D_scalar_advec*vx;D[0][1]=0;                D[0][2]=0;
    3541                 D[1][0]=0;                D[1][1]=D_scalar_advec*vy;D[1][2]=0;
    3542                 D[2][0]=0;                D[2][1]=0;                D[2][2]=D_scalar_advec*vz;
    3543 
    3544                 TripleMultiply(&B_advec[0][0],3,numdof,1,
    3545                                         &D[0][0],3,3,0,
    3546                                         &Bprime_advec[0][0],3,numdof,0,
    3547                                         &Ke->values[0],1);
    3548 
    3549                 /*Transient: */
    3550                 if(reCast<bool,IssmDouble>(dt)){
    3551                         GetNodalFunctionsP1(&L[0], gauss);
    3552                         D_scalar_trans=gauss->weight*Jdet;
    3553 
    3554                         TripleMultiply(&L[0],numdof,1,0,
    3555                                                 &D_scalar_trans,1,1,0,
    3556                                                 &L[0],1,numdof,0,
    3557                                                 &Ke->values[0],1);
    3558                 }
    3559 
    3560                 /*Artificial diffusivity*/
    3561                 if(stabilization==1){
    3562                         /*Build K: */
    3563                         ElementSizes(&hx,&hy,&hz);
    3564                         vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
    3565                         h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
    3566 
    3567                         K[0][0]=h/(2*vel)*vx*vx;  K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz;
    3568                         K[1][0]=h/(2*vel)*vy*vx;  K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz;
    3569                         K[2][0]=h/(2*vel)*vz*vx;  K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz;
    3570 
    3571                         D_scalar_stab=gauss->weight*Jdet;
    3572                         if(reCast<bool,IssmDouble>(dt)) D_scalar_stab=D_scalar_stab*dt;
    3573                         for(i=0;i<3;i++) for(j=0;j<3;j++) K[i][j] = D_scalar_stab*K[i][j];
    3574 
    3575                         GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
    3576 
    3577                         TripleMultiply(&Bprime_advec[0][0],3,numdof,1,
    3578                                                 &K[0][0],3,3,0,
    3579                                                 &Bprime_advec[0][0],3,numdof,0,
    3580                                                 &Ke->values[0],1);
    3581                 }
    3582                 else if(stabilization==2){
    3583                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    3584                         tau_parameter=StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice);
    3585 
    3586                         for(i=0;i<numdof;i++){
    3587                                 for(j=0;j<numdof;j++){
    3588                                         Ke->values[i*numdof+j]+=tau_parameter*D_scalar_advec*
    3589                                           ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]);
    3590                                 }
    3591                         }
    3592                         if(reCast<bool,IssmDouble>(dt)){
    3593                                 for(i=0;i<numdof;i++){
    3594                                         for(j=0;j<numdof;j++){
    3595                                                 Ke->values[i*numdof+j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);
    3596                                         }
    3597                                 }
    3598                         }
    3599                 }
    3600         }
    3601 
    3602         /*Clean up and return*/
    3603         delete gauss;
    3604         return Ke;
    3605 }
    3606 /*}}}*/
    3607 /*FUNCTION Penta::CreateKMatrixEnthalpyShelf {{{*/
    3608 ElementMatrix* Penta::CreateKMatrixEnthalpyShelf(void){
    3609 
    3610         /*Constants*/
    3611         const int    numdof=NDOF1*NUMVERTICES;
    3612 
    3613         /*Intermediaries */
    3614         int        i,j;
    3615         IssmDouble mixed_layer_capacity,thermal_exchange_velocity;
    3616         IssmDouble rho_ice,rho_water,heatcapacity;
    3617         IssmDouble Jdet2d,dt;
    3618         IssmDouble xyz_list[NUMVERTICES][3];
    3619         IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    3620         IssmDouble basis[NUMVERTICES];
    3621         IssmDouble D_scalar;
    3622         GaussPenta *gauss=NULL;
    3623 
    3624         /*Initialize Element matrix and return if necessary*/
    3625         if (!IsOnBed() || !IsFloating()) return NULL;
    3626         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    3627 
    3628         /*Retrieve all inputs and parameters*/
    3629         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    3630         mixed_layer_capacity=matpar->GetMixedLayerCapacity();
    3631         thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
    3632         rho_water=matpar->GetRhoWater();
    3633         rho_ice=matpar->GetRhoIce();
    3634         heatcapacity=matpar->GetHeatCapacity();
    3635         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3636         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    3637 
    3638         /* Start looping on the number of gauss (nodes on the bedrock) */
    3639         gauss=new GaussPenta(0,1,2,2);
    3640         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3641 
    3642                 gauss->GaussPoint(ig);
    3643 
    3644                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3645                 GetNodalFunctionsP1(&basis[0], gauss);
    3646 
    3647                 D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(rho_ice*heatcapacity);
    3648                 if(reCast<bool,IssmDouble>(dt)) D_scalar=dt*D_scalar;
    3649 
    3650                 TripleMultiply(&basis[0],numdof,1,0,
    3651                                         &D_scalar,1,1,0,
    3652                                         &basis[0],1,numdof,0,
    3653                                         &Ke->values[0],1);
    3654         }
    3655 
    3656         /*Clean up and return*/
    3657         delete gauss;
    3658         return Ke;
    3659 }
    3660 /*}}}*/
    3661 /*FUNCTION Penta::CreateKMatrixMelting {{{*/
    3662 ElementMatrix* Penta::CreateKMatrixMelting(void){
    3663 
    3664         if (!IsOnBed()) return NULL;
    3665 
    3666         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    3667         ElementMatrix* Ke=tria->CreateKMatrixMelting();
    3668 
    3669         delete tria->material; delete tria;
    3670         return Ke;
    3671 }
    3672 /*}}}*/
    3673 /*FUNCTION Penta::CreateKMatrixThermal {{{*/
    3674 ElementMatrix* Penta::CreateKMatrixThermal(void){
    3675 
    3676         /*compute all stiffness matrices for this element*/
    3677         ElementMatrix* Ke1=CreateKMatrixThermalVolume();
    3678         ElementMatrix* Ke2=CreateKMatrixThermalShelf();
    3679         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    3680 
    3681         /*clean-up and return*/
    3682         delete Ke1;
    3683         delete Ke2;
    3684         return Ke;
    3685 }
    3686 /*}}}*/
    3687 /*FUNCTION Penta::CreateKMatrixThermalVolume {{{*/
    3688 ElementMatrix* Penta::CreateKMatrixThermalVolume(void){
    3689 
    3690         /*Constants*/
    3691         const int    numdof=NDOF1*NUMVERTICES;
    3692 
    3693         /*Intermediaries */
    3694         int        stabilization;
    3695         int        i,j;
    3696         IssmDouble Jdet,u,v,w,um,vm,wm,vel;
    3697         IssmDouble h,hx,hy,hz,vx,vy,vz;
    3698         IssmDouble gravity,rho_ice,rho_water,kappa;
    3699         IssmDouble heatcapacity,thermalconductivity,dt;
    3700         IssmDouble tau_parameter,diameter;
    3701         IssmDouble xyz_list[NUMVERTICES][3];
    3702         IssmDouble B_conduct[3][numdof];
    3703         IssmDouble B_advec[3][numdof];
    3704         IssmDouble Bprime_advec[3][numdof];
    3705         IssmDouble L[numdof];
    3706         IssmDouble dbasis[3][6];
    3707         IssmDouble D_scalar_conduct,D_scalar_advec;
    3708         IssmDouble D_scalar_trans,D_scalar_stab;
    3709         IssmDouble D[3][3];
    3710         IssmDouble K[3][3]={0.0};
    3711         GaussPenta *gauss=NULL;
    3712 
    3713         /*Initialize Element matrix*/
    3714         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    3715 
    3716         /*Retrieve all inputs and parameters*/
    3717         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3718         rho_water=matpar->GetRhoWater();
    3719         rho_ice=matpar->GetRhoIce();
    3720         gravity=matpar->GetG();
    3721         heatcapacity=matpar->GetHeatCapacity();
    3722         thermalconductivity=matpar->GetThermalConductivity();
    3723         kappa=thermalconductivity/(rho_ice*heatcapacity);
    3724         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    3725         this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    3726         Input* vx_input=inputs->GetInput(VxEnum);      _assert_(vx_input);
    3727         Input* vy_input=inputs->GetInput(VyEnum);      _assert_(vy_input);
    3728         Input* vz_input=inputs->GetInput(VzEnum);      _assert_(vz_input);
    3729         Input* vxm_input=inputs->GetInput(VxMeshEnum); _assert_(vxm_input);
    3730         Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input);
    3731         Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input);
    3732         if (stabilization==2) diameter=MinEdgeLength(&xyz_list[0][0]);
    3733 
    3734         /* Start  looping on the number of gaussian points: */
    3735         gauss=new GaussPenta(2,2);
    3736         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3737 
    3738                 gauss->GaussPoint(ig);
    3739 
    3740                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3741 
    3742                 /*Conduction: */
    3743 
    3744                 GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss);
    3745 
    3746                 D_scalar_conduct=gauss->weight*Jdet*kappa;
    3747                 if(reCast<bool,IssmDouble>(dt)) D_scalar_conduct=D_scalar_conduct*dt;
    3748 
    3749                 D[0][0]=D_scalar_conduct; D[0][1]=0; D[0][2]=0;
    3750                 D[1][0]=0; D[1][1]=D_scalar_conduct; D[1][2]=0;
    3751                 D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_conduct;
    3752 
    3753                 TripleMultiply(&B_conduct[0][0],3,numdof,1,
    3754                                         &D[0][0],3,3,0,
    3755                                         &B_conduct[0][0],3,numdof,0,
    3756                                         &Ke->values[0],1);
    3757 
    3758                 /*Advection: */
    3759                 GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss);
    3760                 GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
    3761 
    3762                 vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
    3763                 vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
    3764                 vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
    3765 
    3766                 D_scalar_advec=gauss->weight*Jdet;
    3767                 if(reCast<bool,IssmDouble>(dt)) D_scalar_advec=D_scalar_advec*dt;
    3768 
    3769                 D[0][0]=D_scalar_advec*vx;    D[0][1]=0;                    D[0][2]=0;
    3770                 D[1][0]=0;                    D[1][1]=D_scalar_advec*vy;    D[1][2]=0;
    3771                 D[2][0]=0;                    D[2][1]=0;                    D[2][2]=D_scalar_advec*vz;
    3772 
    3773                 TripleMultiply(&B_advec[0][0],3,numdof,1,
    3774                                         &D[0][0],3,3,0,
    3775                                         &Bprime_advec[0][0],3,numdof,0,
    3776                                         &Ke->values[0],1);
    3777 
    3778                 /*Transient: */
    3779                 if(reCast<bool,IssmDouble>(dt)){
    3780                         GetNodalFunctionsP1(&L[0], gauss);
    3781                         D_scalar_trans=gauss->weight*Jdet;
    3782 
    3783                         TripleMultiply(&L[0],numdof,1,0,
    3784                                                 &D_scalar_trans,1,1,0,
    3785                                                 &L[0],1,numdof,0,
    3786                                                 &Ke->values[0],1);
    3787                 }
    3788 
    3789                 /*Artifficial diffusivity*/
    3790                 if(stabilization==1){
    3791                         /*Build K: */
    3792                         ElementSizes(&hx,&hy,&hz);
    3793                         vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
    3794                         h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
    3795 
    3796                         K[0][0]=h/(2*vel)*fabs(vx*vx);  K[0][1]=h/(2*vel)*fabs(vx*vy); K[0][2]=h/(2*vel)*fabs(vx*vz);
    3797                         K[1][0]=h/(2*vel)*fabs(vy*vx);  K[1][1]=h/(2*vel)*fabs(vy*vy); K[1][2]=h/(2*vel)*fabs(vy*vz);
    3798                         K[2][0]=h/(2*vel)*fabs(vz*vx);  K[2][1]=h/(2*vel)*fabs(vz*vy); K[2][2]=h/(2*vel)*fabs(vz*vz);
    3799 
    3800                         D_scalar_stab=gauss->weight*Jdet;
    3801                         if(reCast<bool,IssmDouble>(dt)) D_scalar_stab=D_scalar_stab*dt;
    3802                         for(i=0;i<3;i++) for(j=0;j<3;j++) K[i][j] = D_scalar_stab*K[i][j];
    3803 
    3804                         GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
    3805 
    3806                         TripleMultiply(&Bprime_advec[0][0],3,numdof,1,
    3807                                                 &K[0][0],3,3,0,
    3808                                                 &Bprime_advec[0][0],3,numdof,0,
    3809                                                 &Ke->values[0],1);
    3810                 }
    3811                 else if(stabilization==2){
    3812                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    3813                         tau_parameter=StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
    3814 
    3815                         for(i=0;i<numdof;i++){
    3816                                 for(j=0;j<numdof;j++){
    3817                                         Ke->values[i*numdof+j]+=tau_parameter*D_scalar_advec*
    3818                                           ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]);
    3819                                 }
    3820                         }
    3821                         if(reCast<bool,IssmDouble>(dt)){
    3822                                 for(i=0;i<numdof;i++){
    3823                                         for(j=0;j<numdof;j++){
    3824                                                 Ke->values[i*numdof+j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);
    3825                                         }
    3826                                 }
    3827                         }
    3828                 }
    3829         }
    3830 
    3831         /*Clean up and return*/
    3832         delete gauss;
    3833         return Ke;
    3834 }
    3835 /*}}}*/
    3836 /*FUNCTION Penta::CreateKMatrixThermalShelf {{{*/
    3837 ElementMatrix* Penta::CreateKMatrixThermalShelf(void){
    3838 
    3839         /*Constants*/
    3840         const int    numdof=NDOF1*NUMVERTICES;
    3841 
    3842         /*Intermediaries */
    3843         int       i,j;
    3844         IssmDouble mixed_layer_capacity,thermal_exchange_velocity;
    3845         IssmDouble rho_ice,rho_water,heatcapacity;
    3846         IssmDouble Jdet2d,dt;
    3847         IssmDouble xyz_list[NUMVERTICES][3];
    3848         IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    3849         IssmDouble basis[NUMVERTICES];
    3850         IssmDouble D_scalar;
    3851         GaussPenta *gauss=NULL;
    3852 
    3853         /*Initialize Element matrix and return if necessary*/
    3854         if (!IsOnBed() || !IsFloating()) return NULL;
    3855         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    3856 
    3857         /*Retrieve all inputs and parameters*/
    3858         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    3859         mixed_layer_capacity=matpar->GetMixedLayerCapacity();
    3860         thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
    3861         rho_water=matpar->GetRhoWater();
    3862         rho_ice=matpar->GetRhoIce();
    3863         heatcapacity=matpar->GetHeatCapacity();
    3864         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3865         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    3866 
    3867         /* Start looping on the number of gauss (nodes on the bedrock) */
    3868         gauss=new GaussPenta(0,1,2,2);
    3869         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3870 
    3871                 gauss->GaussPoint(ig);
    3872 
    3873                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    3874                 GetNodalFunctionsP1(&basis[0], gauss);
    3875 
    3876                 D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);
    3877                 if(reCast<bool,IssmDouble>(dt)) D_scalar=dt*D_scalar;
    3878 
    3879                 TripleMultiply(&basis[0],numdof,1,0,
    3880                                         &D_scalar,1,1,0,
    3881                                         &basis[0],1,numdof,0,
    3882                                         &Ke->values[0],1);
    3883         }
    3884 
    3885         /*Clean up and return*/
    3886         delete gauss;
    3887         return Ke;
    3888 }
    3889 /*}}}*/
    38903254/*FUNCTION Penta::UpdateBasalConstraintsEnthalpy{{{*/
    38913255void  Penta::UpdateBasalConstraintsEnthalpy(void){
    … …  
    42823646
    42833647}/*}}}*/
    4284 /*FUNCTION Penta::CreateKMatrixAdjointHoriz{{{*/
    4285 ElementMatrix* Penta::CreateKMatrixAdjointHoriz(void){
    4286 
    4287         int approximation;
    4288         inputs->GetInputValue(&approximation,ApproximationEnum);
    4289 
    4290         switch(approximation){
    4291                 case SSAApproximationEnum:
    4292                         return CreateKMatrixAdjointSSA2d();
    4293                 case HOApproximationEnum:
    4294                         return CreateKMatrixAdjointHO();
    4295                 case FSApproximationEnum:
    4296                         return CreateKMatrixAdjointFS();
    4297                 case NoneApproximationEnum:
    4298                         return NULL;
    4299                 default:
    4300                         _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
    4301         }
    4302 }
    4303 /*}}}*/
    4304 /*FUNCTION Penta::CreateKMatrixAdjointSSA2d{{{*/
    4305 ElementMatrix* Penta::CreateKMatrixAdjointSSA2d(void){
    4306 
    4307         /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the
    4308           bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build
    4309           the stiffness matrix. */
    4310         if (!IsOnBed()) return NULL;
    4311 
    4312         /*Depth average some fields*/
    4313         switch(this->material->ObjectEnum()){
    4314                 case MaticeEnum:
    4315                         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    4316                         this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
    4317                         break;
    4318                 default:
    4319                         _error_("material "<<EnumToStringx(this->material->ObjectEnum())<<" not supported");
    4320         }
    4321 
    4322         /*Call Tria function*/
    4323         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    4324         ElementMatrix* Ke=tria->CreateKMatrixAdjointSSA();
    4325         delete tria->material; delete tria;
    4326 
    4327         /*Delete averaged fields*/
    4328         this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    4329         this->material->inputs->DeleteInput(DamageDbarEnum);
    4330 
    4331         /*clean up and return*/
    4332         return Ke;
    4333 }
    4334 /*}}}*/
    4335 /*FUNCTION Penta::CreateKMatrixAdjointHO{{{*/
    4336 ElementMatrix* Penta::CreateKMatrixAdjointHO(void){
    4337 
    4338         /*Intermediaries */
    4339         int        i,j;
    4340         bool       incomplete_adjoint;
    4341         IssmDouble xyz_list[NUMVERTICES][3];
    4342         IssmDouble Jdet;
    4343         IssmDouble eps1dotdphii,eps1dotdphij;
    4344         IssmDouble eps2dotdphii,eps2dotdphij;
    4345         IssmDouble mu_prime;
    4346         IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    4347         IssmDouble eps1[3],eps2[3];
    4348         IssmDouble dphi[3][NUMVERTICES];
    4349         GaussPenta *gauss=NULL;
    4350 
    4351         /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
    4352         parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
    4353         ElementMatrix* Ke=CreateKMatrixStressbalanceHO();
    4354         if(incomplete_adjoint) return Ke;
    4355 
    4356         /*Retrieve all inputs and parameters*/
    4357         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4358         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    4359         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    4360 
    4361         /* Start  looping on the number of gaussian points: */
    4362         gauss=new GaussPenta(5,5);
    4363         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4364 
    4365                 gauss->GaussPoint(ig);
    4366 
    4367                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    4368                 GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    4369 
    4370                 this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    4371                 material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    4372                 eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    4373                 eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
    4374                 eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
    4375 
    4376                 for(i=0;i<6;i++){
    4377                         for(j=0;j<6;j++){
    4378                                 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
    4379                                 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
    4380                                 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
    4381                                 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
    4382 
    4383                                 Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
    4384                                 Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
    4385                                 Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
    4386                                 Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
    4387                         }
    4388                 }
    4389         }
    4390 
    4391         /*Transform Coordinate System*/
    4392         ::TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
    4393 
    4394         /*Clean up and return*/
    4395         delete gauss;
    4396         return Ke;
    4397 }
    4398 /*}}}*/
    4399 /*FUNCTION Penta::CreateKMatrixAdjointFS{{{*/
    4400 ElementMatrix* Penta::CreateKMatrixAdjointFS(void){
    4401 
    4402         /*Intermediaries */
    4403         int        i,j;
    4404         bool       incomplete_adjoint;
    4405         IssmDouble xyz_list[NUMVERTICES][3];
    4406         IssmDouble Jdet;
    4407         IssmDouble eps1dotdphii,eps1dotdphij;
    4408         IssmDouble eps2dotdphii,eps2dotdphij;
    4409         IssmDouble eps3dotdphii,eps3dotdphij;
    4410         IssmDouble mu_prime;
    4411         IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    4412         IssmDouble eps1[3],eps2[3],eps3[3];
    4413         IssmDouble dphi[3][NUMVERTICES];
    4414         GaussPenta *gauss=NULL;
    4415 
    4416         /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
    4417         parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
    4418         ElementMatrix* Ke=CreateKMatrixStressbalanceFS();
    4419         if(incomplete_adjoint) return Ke;
    4420 
    4421         /*Fetch number of nodes and dof for this finite element*/
    4422         int vnumnodes = this->NumberofNodesVelocity();
    4423         int pnumnodes = this->NumberofNodesPressure();
    4424         int numdof    = vnumnodes*NDOF3 + pnumnodes*NDOF1;
    4425 
    4426         /*Retrieve all inputs and parameters*/
    4427         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4428         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    4429         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    4430         Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
    4431 
    4432         /* Start  looping on the number of gaussian points: */
    4433         gauss=new GaussPenta(5,5);
    4434         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4435 
    4436                 gauss->GaussPoint(ig);
    4437 
    4438                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    4439                 GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    4440 
    4441                 this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    4442                 material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    4443                 eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
    4444                 eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
    4445                 eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
    4446 
    4447                 for(i=0;i<6;i++){
    4448                         for(j=0;j<6;j++){
    4449                                 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
    4450                                 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
    4451                                 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
    4452                                 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
    4453                                 eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i];
    4454                                 eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j];
    4455 
    4456                                 Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
    4457                                 Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
    4458                                 Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
    4459 
    4460                                 Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
    4461                                 Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
    4462                                 Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
    4463 
    4464                                 Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
    4465                                 Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
    4466                                 Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
    4467                         }
    4468                 }
    4469         }
    4470 
    4471         /*Transform Coordinate System*/
    4472         ::TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZEnum);
    4473 
    4474         /*Clean up and return*/
    4475         delete gauss;
    4476         return Ke;
    4477 }
    4478 /*}}}*/
    44793648/*FUNCTION Penta::GradientIndexing{{{*/
    44803649void Penta::GradientIndexing(int* indexing,int control_index){
    … …  
    53514520}
    53524521/*}}}*/
    5353 /*FUNCTION Penta::CreateKMatrixCouplingSSAHO{{{*/
    5354 ElementMatrix* Penta::CreateKMatrixCouplingSSAHO(void){
    5355 
    5356         /*compute all stiffness matrices for this element*/
    5357         ElementMatrix* Ke1=CreateKMatrixCouplingSSAHOViscous();
    5358         ElementMatrix* Ke2=CreateKMatrixCouplingSSAHOFriction();
    5359         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    5360 
    5361         /*clean-up and return*/
    5362         delete Ke1;
    5363         delete Ke2;
    5364         return Ke;
    5365 }
    5366 /*}}}*/
    5367 /*FUNCTION Penta::CreateKMatrixCouplingSSAHOViscous{{{*/
    5368 ElementMatrix* Penta::CreateKMatrixCouplingSSAHOViscous(void){
    5369 
    5370         /*Constants*/
    5371         const int numnodes    = 2 *NUMVERTICES;
    5372         const int numdofm     = NDOF2 *NUMVERTICES2D;
    5373         const int numdofp     = NDOF2 *NUMVERTICES;
    5374         const int numdoftotal = 2 *NDOF2*NUMVERTICES;
    5375 
    5376         /*Intermediaries */
    5377         int         i,j;
    5378         IssmDouble  Jdet;
    5379         IssmDouble  viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
    5380         IssmDouble  epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    5381         IssmDouble  xyz_list[NUMVERTICES][3];
    5382         IssmDouble  B[3][numdofp];
    5383         IssmDouble  Bprime[3][numdofm];
    5384         IssmDouble  D[3][3]={0.0};            // material matrix, simple scalar matrix.
    5385         IssmDouble  D_scalar;
    5386         IssmDouble  Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix
    5387         IssmDouble  Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point.
    5388         GaussPenta *gauss=NULL;
    5389         GaussTria  *gauss_tria=NULL;
    5390         Node       *node_list[numnodes];
    5391         int         cs_list[numnodes];
    5392 
    5393         /*Find penta on bed as HO must be coupled to the dofs on the bed: */
    5394         Penta* pentabase=(Penta*)GetBasalElement();
    5395         Tria*  tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
    5396 
    5397         /*prepare node list*/
    5398         for(i=0;i<NUMVERTICES;i++){
    5399                 node_list[i+0*NUMVERTICES] = pentabase->nodes[i];
    5400                 node_list[i+1*NUMVERTICES] = this->nodes[i];
    5401                 cs_list[i+0*NUMVERTICES] = XYEnum;
    5402                 cs_list[i+1*NUMVERTICES] = XYEnum;
    5403         }
    5404 
    5405         /*Initialize Element matrix*/
    5406         ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,SSAApproximationEnum);
    5407         ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,HOApproximationEnum);
    5408         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    5409         delete Ke1; delete Ke2;
    5410 
    5411         /* Get node coordinates and dof list: */
    5412         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5413         this->parameters->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
    5414         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    5415         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    5416         Input* vxold_input=inputs->GetInput(VxPicardEnum); _assert_(vxold_input);
    5417         Input* vyold_input=inputs->GetInput(VyPicardEnum); _assert_(vyold_input);
    5418 
    5419         /* Start  looping on the number of gaussian points: */
    5420         gauss=new GaussPenta(5,5);
    5421         gauss_tria=new GaussTria();
    5422         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5423 
    5424                 gauss->GaussPoint(ig);
    5425                 gauss->SynchronizeGaussBase(gauss_tria);
    5426 
    5427                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    5428                 GetBSSAHO(&B[0][0], &xyz_list[0][0], gauss);
    5429                 tria->GetBprimeSSA(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    5430 
    5431                 this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    5432                 this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    5433                 material->GetViscosity3d(&viscosity, &epsilon[0]);
    5434                 material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    5435 
    5436                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    5437                 D_scalar=2*newviscosity*gauss->weight*Jdet;
    5438                 for (i=0;i<3;i++) D[i][i]=D_scalar;
    5439 
    5440                 TripleMultiply( &B[0][0],3,numdofp,1,
    5441                                         &D[0][0],3,3,0,
    5442                                         &Bprime[0][0],3,numdofm,0,
    5443                                         &Ke_gg_gaussian[0][0],0);
    5444 
    5445                 for( i=0; i<numdofp; i++) for(j=0;j<numdofm; j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    5446         }
    5447         for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j];
    5448         for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j][i];
    5449 
    5450         /*Transform Coordinate System*/
    5451         ::TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
    5452 
    5453         /*Clean-up and return*/
    5454         delete tria->material; delete tria;
    5455         delete gauss;
    5456         delete gauss_tria;
    5457         return Ke;
    5458 }
    5459 /*}}}*/
    5460 /*FUNCTION Penta::CreateKMatrixCouplingSSAHOFriction{{{*/
    5461 ElementMatrix* Penta::CreateKMatrixCouplingSSAHOFriction(void){
    5462 
    5463         /*Constants*/
    5464         const int numnodes    = 2 *NUMVERTICES;
    5465         const int numdof      = NDOF2 *NUMVERTICES;
    5466         const int numdoftotal = NDOF4 *NUMVERTICES;
    5467 
    5468         /*Intermediaries */
    5469         int       i,j,analysis_type;
    5470         IssmDouble Jdet2d,alpha2;
    5471         IssmDouble xyz_list[NUMVERTICES][3];
    5472         IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
    5473         IssmDouble L[2][numdof];
    5474         IssmDouble DL[2][2]                  ={{ 0,0 },{0,0}}; //for basal drag
    5475         IssmDouble DL_scalar;
    5476         IssmDouble Ke_gg[numdof][numdof]     ={0.0};
    5477         IssmDouble Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    5478         Friction  *friction = NULL;
    5479         GaussPenta *gauss=NULL;
    5480         Node       *node_list[numnodes];
    5481         int         cs_list[numnodes];
    5482 
    5483         /*Initialize Element matrix and return if necessary*/
    5484         if(IsFloating() || !IsOnBed()) return NULL;
    5485         ElementMatrix* Ke1=new ElementMatrix(nodes,NUMVERTICES,this->parameters,SSAApproximationEnum);
    5486         ElementMatrix* Ke2=new ElementMatrix(nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
    5487         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    5488         delete Ke1; delete Ke2;
    5489 
    5490         /*Prepare node list*/
    5491         for(i=0;i<NUMVERTICES;i++){
    5492                 node_list[i+0*NUMVERTICES] = this->nodes[i];
    5493                 node_list[i+1*NUMVERTICES] = this->nodes[i];
    5494                 cs_list[i+0*NUMVERTICES] = XYEnum;
    5495                 cs_list[i+1*NUMVERTICES] = XYEnum;
    5496         }
    5497 
    5498         /*retrieve inputs :*/
    5499         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5500         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    5501         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    5502         Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
    5503         Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
    5504         Input* vz_input=inputs->GetInput(VzEnum);           _assert_(vz_input);
    5505 
    5506         /*build friction object, used later on: */
    5507         friction=new Friction(this,2);
    5508 
    5509         /* Start  looping on the number of gaussian points: */
    5510         gauss=new GaussPenta(0,1,2,2);
    5511         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5512 
    5513                 gauss->GaussPoint(ig);
    5514 
    5515                 /*Friction: */
    5516                 friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    5517 
    5518                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    5519                 GetBHOFriction(&L[0][0],gauss);
    5520 
    5521                 DL_scalar=alpha2*gauss->weight*Jdet2d;
    5522                 for (i=0;i<2;i++) DL[i][i]=DL_scalar;
    5523 
    5524                 /*  Do the triple producte tL*D*L: */
    5525                 TripleMultiply( &L[0][0],2,numdof,1,
    5526                                         &DL[0][0],2,2,0,
    5527                                         &L[0][0],2,numdof,0,
    5528                                         &Ke_gg_gaussian[0][0],0);
    5529 
    5530                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    5531         }
    5532 
    5533         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i][j];
    5534         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i][j];
    5535 
    5536         /*Transform Coordinate System*/
    5537         ::TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
    5538 
    5539         /*Clean up and return*/
    5540         delete gauss;
    5541         delete friction;
    5542         return Ke;
    5543 }
    5544 /*}}}*/
    5545 /*FUNCTION Penta::CreateKMatrixCouplingSSAFS{{{*/
    5546 ElementMatrix* Penta::CreateKMatrixCouplingSSAFS(void){
    5547 
    5548         /*compute all stiffness matrices for this element*/
    5549         ElementMatrix* Ke1=CreateKMatrixCouplingSSAFSViscous();
    5550         ElementMatrix* Ke2=CreateKMatrixCouplingSSAFSFriction();
    5551         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    5552 
    5553         /*clean-up and return*/
    5554         delete Ke1;
    5555         delete Ke2;
    5556         return Ke;
    5557 }
    5558 /*}}}*/
    5559 /*FUNCTION Penta::CreateKMatrixCouplingSSAFSViscous{{{*/
    5560 ElementMatrix* Penta::CreateKMatrixCouplingSSAFSViscous(void){
    5561 
    5562         /*Constants*/
    5563         const int numdofm     = NDOF2 *NUMVERTICES2D;
    5564         const int numdofs     = NDOF4 *NUMVERTICES + NDOF3;
    5565         const int numdoftotal = 2 *numdofm+numdofs;
    5566 
    5567         /*Intermediaries */
    5568         int        i,j;
    5569         IssmDouble Jdet;
    5570         IssmDouble viscosity,FSreconditioning; //viscosity
    5571         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    5572         IssmDouble xyz_list[NUMVERTICES][3];
    5573         IssmDouble B[4][numdofs];
    5574         IssmDouble Bprime[4][numdofm];
    5575         IssmDouble B2[3][numdofm];
    5576         IssmDouble Bprime2[3][numdofs];
    5577         IssmDouble D[4][4]={0.0};            // material matrix, simple scalar matrix.
    5578         IssmDouble D2[3][3]={0.0};            // material matrix, simple scalar matrix.
    5579         IssmDouble D_scalar;
    5580         IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix
    5581         IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix
    5582         IssmDouble Ke_gg_gaussian[numdofs][numdofm]; //stiffness matrix evaluated at the gaussian point.
    5583         IssmDouble Ke_gg_gaussian2[numdofm][numdofs]; //stiffness matrix evaluated at the gaussian point.
    5584         GaussPenta *gauss=NULL;
    5585         GaussTria  *gauss_tria=NULL;
    5586         Node       *node_list[20];
    5587 
    5588         /*Find penta on bed as FS must be coupled to the dofs on the bed: */
    5589         Penta* pentabase=(Penta*)GetBasalElement();
    5590         Tria* tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
    5591 
    5592         int vnumnodes = this->NumberofNodesVelocity();
    5593         int pnumnodes = this->NumberofNodesPressure();
    5594         int numnodes  = 2*vnumnodes-1+pnumnodes;
    5595 
    5596         /*Prepare node list*/
    5597         int* cs_list = xNew<int>(2*vnumnodes+pnumnodes);
    5598         for(i=0;i<vnumnodes-1;i++){
    5599                 node_list[i] = pentabase->nodes[i];
    5600                 cs_list[i] = XYEnum;
    5601         }
    5602         for(i=0;i<vnumnodes;i++){
    5603                 node_list[i+vnumnodes-1] = this->nodes[i];
    5604                 cs_list[i+vnumnodes-1] = XYZEnum;
    5605         }
    5606         for(i=0;i<pnumnodes;i++){
    5607                 node_list[2*vnumnodes-1+i] = this->nodes[vnumnodes+i];
    5608                 cs_list[2*vnumnodes-1+i] = PressureEnum;
    5609         }
    5610 
    5611         /*Initialize Element matrix and return if necessary*/
    5612         ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,    this->parameters,SSAApproximationEnum);
    5613         ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,2*NUMVERTICES+1,this->parameters,FSvelocityEnum);
    5614         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    5615         delete Ke1; delete Ke2;
    5616 
    5617         /* Get node coordinates and dof list: */
    5618         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5619         parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    5620         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    5621         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    5622         Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
    5623 
    5624         /* Start  looping on the number of gaussian points: */
    5625         gauss=new GaussPenta(5,5);
    5626         gauss_tria=new GaussTria();
    5627         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5628 
    5629                 gauss->GaussPoint(ig);
    5630                 gauss->SynchronizeGaussBase(gauss_tria);
    5631 
    5632                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    5633                 GetBSSAFS(&B[0][0], &xyz_list[0][0], gauss);
    5634                 tria->GetBprimeSSAFS(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    5635                 tria->GetBSSAFS(&B2[0][0], &xyz_list[0][0], gauss_tria);
    5636                 GetBprimeSSAFS(&Bprime2[0][0], &xyz_list[0][0], gauss);
    5637 
    5638                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    5639                 material->GetViscosity3dFS(&viscosity, &epsilon[0]);
    5640 
    5641                 D_scalar=2*viscosity*gauss->weight*Jdet;
    5642                 for (i=0;i<3;i++) D[i][i]=D_scalar;
    5643                 D[3][3]=-gauss->weight*Jdet*FSreconditioning;
    5644                 for (i=0;i<3;i++) D2[i][i]=D_scalar;
    5645 
    5646                 TripleMultiply( &B[0][0],4,numdofs,1,
    5647                                         &D[0][0],4,4,0,
    5648                                         &Bprime[0][0],4,numdofm,0,
    5649                                         &Ke_gg_gaussian[0][0],0);
    5650 
    5651                 TripleMultiply( &B2[0][0],3,numdofm,1,
    5652                                         &D2[0][0],3,3,0,
    5653                                         &Bprime2[0][0],3,numdofs,0,
    5654                                         &Ke_gg_gaussian2[0][0],0);
    5655 
    5656                 for( i=0; i<numdofs; i++) for(j=0;j<numdofm; j++)      Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    5657                 for( i=0; i<numdofm; i++)      for(j=0;j<numdofs; j++) Ke_gg2[i][j]+=Ke_gg_gaussian2[i][j];
    5658         }
    5659         for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++)      Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j];
    5660         for(i=0;i<numdofm;i++)      for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j];
    5661 
    5662         /*Transform Coordinate System*/
    5663         ::TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
    5664 
    5665         /*Clean-up and return*/
    5666         xDelete<int>(cs_list);
    5667         delete tria->material; delete tria;
    5668         delete gauss;
    5669         delete gauss_tria;
    5670         return Ke;
    5671 }
    5672 /*}}}*/
    5673 /*FUNCTION Penta::CreateKMatrixCouplingSSAFSFriction {{{*/
    5674 ElementMatrix* Penta::CreateKMatrixCouplingSSAFSFriction(void){
    5675         /*Constants*/
    5676         const int numdofs   = (NUMVERTICES+1)*NDOF3 + NUMVERTICES*NDOF1;
    5677         const int numdofm   = NUMVERTICES *NDOF2;
    5678         const int numdof2d  = NUMVERTICES2D *NDOF3;
    5679         const int numdof2dm = NUMVERTICES2D *NDOF2;
    5680         const int numdoftot = NUMVERTICES*2 + (NUMVERTICES+1)*3 +NUMVERTICES; // HO + FS vel + FS Pressure
    5681 
    5682         /*Intermediaries */
    5683         int        i,j;
    5684         int        analysis_type,approximation;
    5685         IssmDouble FSreconditioning;
    5686         IssmDouble viscosity,alpha2_gauss,Jdet2d;
    5687         IssmDouble bed_normal[3];
    5688         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    5689         IssmDouble xyz_list[NUMVERTICES][3];
    5690         IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    5691         IssmDouble LSSAFS[8][numdof2dm];
    5692         IssmDouble LprimeSSAFS[8][numdofs];
    5693         IssmDouble DLSSAFS[8][8]={0.0};
    5694         IssmDouble LFSSSA[4][numdof2d];
    5695         IssmDouble LprimeFSSSA[4][numdof2dm];
    5696         IssmDouble DLFSSSA[4][4]={0.0};
    5697         IssmDouble Ke_drag_gaussian[numdof2dm][numdofs];
    5698         IssmDouble Ke_drag_gaussian2[numdof2d][numdof2dm];
    5699         Friction*  friction=NULL;
    5700         GaussPenta *gauss=NULL;
    5701         Node       *node_list[20];
    5702 
    5703         /*If on water or not FS, skip stiffness: */
    5704         inputs->GetInputValue(&approximation,ApproximationEnum);
    5705         if(IsFloating() || !IsOnBed()) return NULL;
    5706 
    5707         int vnumnodes = this->NumberofNodesVelocity();
    5708         int pnumnodes = this->NumberofNodesPressure();
    5709         int numnodes  = 2*vnumnodes-1+pnumnodes;
    5710 
    5711         /*Prepare node list*/
    5712         int* cs_list = xNew<int>(2*vnumnodes+pnumnodes);
    5713         for(i=0;i<vnumnodes-1;i++){
    5714                 node_list[i] = this->nodes[i];
    5715                 cs_list[i] = XYEnum;
    5716         }
    5717         for(i=0;i<vnumnodes;i++){
    5718                 node_list[i+vnumnodes-1] = this->nodes[i];
    5719                 cs_list[i+vnumnodes-1] = XYZEnum;
    5720         }
    5721         for(i=0;i<pnumnodes;i++){
    5722                 node_list[2*vnumnodes-1+i] = this->nodes[vnumnodes+i];
    5723                 cs_list[2*vnumnodes-1+i] = PressureEnum;
    5724         }
    5725 
    5726         ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,        this->parameters,SSAApproximationEnum);
    5727         ElementMatrix* Ke2=new ElementMatrix(this->nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    5728         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    5729         delete Ke1; delete Ke2;
    5730 
    5731         /*Retrieve all inputs and parameters*/
    5732         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5733         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    5734         parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    5735         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    5736         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    5737         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    5738         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    5739 
    5740         /*build friction object, used later on: */
    5741         friction=new Friction(this,3);
    5742 
    5743         /* Start  looping on the number of gaussian points: */
    5744         gauss=new GaussPenta(0,1,2,2);
    5745         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5746 
    5747                 gauss->GaussPoint(ig);
    5748 
    5749                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    5750                 GetLSSAFS(&LSSAFS[0][0], gauss);
    5751                 GetLprimeSSAFS(&LprimeSSAFS[0][0], &xyz_list[0][0], gauss);
    5752                 GetLFSSSA(&LFSSSA[0][0], gauss);
    5753                 GetLprimeFSSSA(&LprimeFSSSA[0][0], &xyz_list[0][0], gauss);
    5754 
    5755                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    5756                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    5757 
    5758                 NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    5759                 friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
    5760 
    5761                 DLSSAFS[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
    5762                 DLSSAFS[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
    5763                 DLSSAFS[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
    5764                 DLSSAFS[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
    5765                 DLSSAFS[4][4]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
    5766                 DLSSAFS[5][5]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
    5767                 DLSSAFS[6][6]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[0];
    5768                 DLSSAFS[7][7]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[1];
    5769 
    5770                 DLFSSSA[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
    5771                 DLFSSSA[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
    5772                 DLFSSSA[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
    5773                 DLFSSSA[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
    5774 
    5775                 TripleMultiply( &LSSAFS[0][0],8,numdof2dm,1,
    5776                                         &DLSSAFS[0][0],8,8,0,
    5777                                         &LprimeSSAFS[0][0],8,numdofs,0,
    5778                                         &Ke_drag_gaussian[0][0],0);
    5779 
    5780                 TripleMultiply( &LFSSSA[0][0],4,numdof2d,1,
    5781                                         &DLFSSSA[0][0],4,4,0,
    5782                                         &LprimeFSSSA[0][0],4,numdof2dm,0,
    5783                                         &Ke_drag_gaussian2[0][0],0);
    5784                 for(i=0;i<numdof2dm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftot+j+numdofm]+=Ke_drag_gaussian[i][j];
    5785                 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2dm;j++) Ke->values[(i+numdofm)*numdoftot+j]+=Ke_drag_gaussian2[i][j];
    5786         }
    5787 
    5788         /*Transform Coordinate System*/
    5789         ::TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
    5790 
    5791         /*Clean up and return*/
    5792         xDelete<int>(cs_list);
    5793         delete gauss;
    5794         delete friction;
    5795         return Ke;
    5796 }
    5797 /*}}}*/
    5798 /*FUNCTION Penta::CreateKMatrixCouplingHOFS{{{*/
    5799 ElementMatrix* Penta::CreateKMatrixCouplingHOFS(void){
    5800 
    5801         /*Constants*/
    5802         const int numnodes       = 3 *NUMVERTICES+1;
    5803         const int numdofp        = NDOF2 *NUMVERTICES;
    5804         const int numdofs        = NDOF4 * 6 + NDOF3;
    5805         const int numdoftotal    = (NDOF2+NDOF4) *NUMVERTICES + NDOF3;
    5806 
    5807         /*Intermediaries*/
    5808         int   i,j,init;
    5809         Node  *node_list[NUMVERTICES*3+1];
    5810         int   cs_list[NUMVERTICES*3+1];
    5811         int   cs_list2[NUMVERTICES*2+1];
    5812 
    5813         /*Some parameters needed*/
    5814         init=this->element_type;
    5815 
    5816         /*prepare node list*/
    5817         for(i=0;i<NUMVERTICES+1;i++){
    5818                 node_list[i+NUMVERTICES] = this->nodes[i];
    5819                 cs_list[i+NUMVERTICES]   = XYZEnum;
    5820                 cs_list2[i]              = XYZEnum;
    5821         }
    5822         for(i=0;i<NUMVERTICES;i++){
    5823                 node_list[i]                 = this->nodes[i];
    5824                 node_list[i+2*NUMVERTICES+1] = this->nodes[i+NUMVERTICES+1];
    5825                 cs_list[i]                   = XYEnum;
    5826                 cs_list[i+2*NUMVERTICES+1]   = PressureEnum;
    5827                 cs_list2[i+NUMVERTICES+1]    = PressureEnum;
    5828         }
    5829 
    5830         /*compute all stiffness matrices for this element*/
    5831         ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
    5832         ElementMatrix* Ke2=new ElementMatrix(this->nodes,2*NUMVERTICES+1,this->parameters,FSvelocityEnum);
    5833         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    5834         delete Ke1; delete Ke2;
    5835 
    5836         /*Compute HO Matrix with P1 element type\n");*/
    5837         this->element_type=P1Enum;
    5838         Ke1=CreateKMatrixStressbalanceHO(); TransformInvStiffnessMatrixCoord(Ke1,XYEnum);
    5839         this->element_type=init;
    5840         /*Compute FS Matrix and condense it \n");*/
    5841         Ke2=CreateKMatrixStressbalanceFS(); TransformInvStiffnessMatrixCoord(Ke2,node_list,2*NUMVERTICES+1,cs_list2);
    5842         int indices[3]={18,19,20};
    5843         Ke2->StaticCondensation(3,&indices[0]);
    5844 
    5845         for(i=0;i<numdofs;i++) for(j=0;j<NUMVERTICES;j++){
    5846                 Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+0]+=Ke2->values[i*numdofs+NDOF3*j+0];
    5847                 Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+1]+=Ke2->values[i*numdofs+NDOF3*j+1];
    5848         }
    5849         for(i=0;i<numdofp;i++) for(j=0;j<NUMVERTICES;j++){
    5850                 Ke->values[i*numdoftotal+numdofp+NDOF3*j+0]+=Ke1->values[i*numdofp+NDOF2*j+0];
    5851                 Ke->values[i*numdoftotal+numdofp+NDOF3*j+1]+=Ke1->values[i*numdofp+NDOF2*j+1];
    5852         }
    5853 
    5854         /*Transform Coordinate System*/ //Do not transform, already done in the matrices
    5855         ::TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
    5856 
    5857         /*clean-up and return*/
    5858         delete Ke1;
    5859         delete Ke2;
    5860         return Ke;
    5861 }
    5862 //*}}}*/
    5863 /*FUNCTION Penta::CreateKMatrixStressbalanceHoriz {{{*/
    5864 ElementMatrix* Penta::CreateKMatrixStressbalanceHoriz(void){
    5865 
    5866         int approximation;
    5867         inputs->GetInputValue(&approximation,ApproximationEnum);
    5868         switch(approximation){
    5869                 case SSAApproximationEnum:
    5870                         return CreateKMatrixStressbalanceSSA2d();
    5871                 case L1L2ApproximationEnum:
    5872                         return CreateKMatrixStressbalanceL1L2();
    5873                 case HOApproximationEnum:
    5874                         return CreateKMatrixStressbalanceHO();
    5875                 case FSApproximationEnum:
    5876                         return CreateKMatrixStressbalanceFS();
    5877                 case SIAApproximationEnum:
    5878                         return NULL;
    5879                 case NoneApproximationEnum:
    5880                         return NULL;
    5881                 case SSAHOApproximationEnum:
    5882                         return CreateKMatrixStressbalanceSSAHO();
    5883                 case SSAFSApproximationEnum:
    5884                         return CreateKMatrixStressbalanceSSAFS();
    5885                 case HOFSApproximationEnum:
    5886                         return CreateKMatrixStressbalanceHOFS();
    5887                 default:
    5888                         _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
    5889         }
    5890 }
    5891 /*}}}*/
    5892 /*FUNCTION Penta::CreateKMatrixStressbalanceSIA{{{*/
    5893 ElementMatrix* Penta::CreateKMatrixStressbalanceSIA(void){
    5894 
    5895         /*Intermediaries*/
    5896         IssmDouble connectivity[2];
    5897         IssmDouble one0,one1;
    5898         int        i,i0,i1,j0,j1;
    5899 
    5900         /*Fetch number of nodes and dof for this finite element*/
    5901         int numnodes = this->NumberofNodes(); _assert_(numnodes==6);
    5902         int numdof   = numnodes*NDOF2;
    5903 
    5904         /*Initialize Element matrix*/
    5905         ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    5906 
    5907         /*3 vertical edges*/
    5908         for(i=0;i<3;i++){
    5909 
    5910                 /*2 dofs of first node*/
    5911                 i0=2*i;     i1=2*i+1;
    5912                 /*2 dofs of second node*/
    5913                 j0=2*(i+3); j1=2*(i+3)+1;
    5914 
    5915                 /*Find connectivity for the two nodes*/
    5916                 connectivity[0]=(IssmDouble)vertices[i]->Connectivity();
    5917                 connectivity[1]=(IssmDouble)vertices[i+3]->Connectivity();
    5918                 one0=1./connectivity[0];
    5919                 one1=1./connectivity[1];
    5920 
    5921                 /*Create matrix for these two nodes*/
    5922                 if (IsOnBed() && IsOnSurface()){
    5923                         Ke->values[i0*numdof+i0] = +one0;
    5924                         Ke->values[i1*numdof+i1] = +one0;
    5925                         Ke->values[j0*numdof+i0] = -one1;
    5926                         Ke->values[j0*numdof+j0] = +one1;
    5927                         Ke->values[j1*numdof+i1] = -one1;
    5928                         Ke->values[j1*numdof+j1] = +one1;
    5929                 }
    5930                 else if (IsOnBed()){
    5931                         Ke->values[i0*numdof+i0] = one0;
    5932                         Ke->values[i1*numdof+i1] = one0;
    5933                         Ke->values[j0*numdof+i0] = -2.*one1;
    5934                         Ke->values[j0*numdof+j0] = +2.*one1;
    5935                         Ke->values[j1*numdof+i1] = -2.*one1;
    5936                         Ke->values[j1*numdof+j1] = +2.*one1;
    5937                 }
    5938                 else if (IsOnSurface()){
    5939                         Ke->values[j0*numdof+i0] = -one1;
    5940                         Ke->values[j0*numdof+j0] = +one1;
    5941                         Ke->values[j1*numdof+i1] = -one1;
    5942                         Ke->values[j1*numdof+j1] = +one1;
    5943                 }
    5944                 else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity
    5945                         Ke->values[j0*numdof+i0] = -2.*one1;
    5946                         Ke->values[j0*numdof+j0] = +2.*one1;
    5947                         Ke->values[j1*numdof+i1] = -2.*one1;
    5948                         Ke->values[j1*numdof+j1] = +2.*one1;
    5949                 }
    5950         }
    5951 
    5952         /*Clean up and return*/
    5953         return Ke;
    5954 }/*}}}*/
    5955 /*FUNCTION Penta::CreateKMatrixStressbalanceSSA2d{{{*/
    5956 ElementMatrix* Penta::CreateKMatrixStressbalanceSSA2d(void){
    5957 
    5958         /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the
    5959           bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build
    5960           the stiffness matrix. */
    5961         if (!IsOnBed()) return NULL;
    5962 
    5963         /*Depth average some fields*/
    5964         switch(this->material->ObjectEnum()){
    5965                 case MaticeEnum:
    5966                         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    5967                         this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
    5968                         break;
    5969                 default:
    5970                         _error_("material "<<EnumToStringx(this->material->ObjectEnum())<<" not supported");
    5971         }
    5972 
    5973         /*Call Tria function*/
    5974         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    5975         ElementMatrix* Ke=tria->CreateKMatrixStressbalanceSSA();
    5976         delete tria->material; delete tria;
    5977 
    5978         /*Delete averaged fields*/
    5979         this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    5980         this->material->inputs->DeleteInput(DamageDbarEnum);
    5981 
    5982         /*clean up and return*/
    5983         return Ke;
    5984 }
    5985 /*}}}*/
    5986 /*FUNCTION Penta::CreateKMatrixStressbalanceSSA3d{{{*/
    5987 ElementMatrix* Penta::CreateKMatrixStressbalanceSSA3d(void){
    5988 
    5989         /*compute all stiffness matrices for this element*/
    5990         ElementMatrix* Ke1=CreateKMatrixStressbalanceSSA3dViscous();
    5991         ElementMatrix* Ke2=CreateKMatrixStressbalanceSSA3dFriction();
    5992         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    5993 
    5994         /*clean-up and return*/
    5995         delete Ke1;
    5996         delete Ke2;
    5997         return Ke;
    5998 }
    5999 /*}}}*/
    6000 /*FUNCTION Penta::CreateKMatrixStressbalanceSSA3dViscous{{{*/
    6001 ElementMatrix* Penta::CreateKMatrixStressbalanceSSA3dViscous(void){
    6002 
    6003         /*Constants*/
    6004         const int    numdof2d=2*NUMVERTICES2D;
    6005 
    6006         /*Intermediaries */
    6007         int         i,j,approximation;
    6008         IssmDouble  Jdet;
    6009         IssmDouble  viscosity , oldviscosity, newviscosity, viscosity_overshoot;
    6010         IssmDouble  epsilon[5],oldepsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
    6011         IssmDouble  epsilons[6];                    //6 for FS
    6012         IssmDouble  xyz_list[NUMVERTICES][3];
    6013         IssmDouble  B[3][numdof2d];
    6014         IssmDouble  Bprime[3][numdof2d];
    6015         IssmDouble  D[3][3]= {0.0};                 // material matrix, simple scalar matrix.
    6016         IssmDouble  D_scalar;
    6017         IssmDouble  Ke_gg_gaussian[numdof2d][numdof2d];
    6018         Tria       *tria       = NULL;
    6019         Penta      *pentabase  = NULL;
    6020         GaussPenta *gauss      = NULL;
    6021         GaussTria  *gauss_tria = NULL;
    6022 
    6023         /*Find penta on bed as this is a SSA elements: */
    6024         pentabase=(Penta*)GetBasalElement();
    6025         tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
    6026 
    6027         /*Initialize Element matrix*/
    6028         ElementMatrix* Ke=new ElementMatrix(tria->nodes,NUMVERTICES2D,this->parameters,SSAApproximationEnum);
    6029         inputs->GetInputValue(&approximation,ApproximationEnum);
    6030 
    6031         /*Retrieve all inputs and parameters*/
    6032         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6033         this->parameters->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
    6034         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    6035         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    6036         Input* vxold_input=inputs->GetInput(VxPicardEnum); _assert_(vxold_input);
    6037         Input* vyold_input=inputs->GetInput(VyPicardEnum); _assert_(vyold_input);
    6038         Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
    6039 
    6040         /* Start  looping on the number of gaussian points: */
    6041         gauss=new GaussPenta(5,5);
    6042         gauss_tria=new GaussTria();
    6043         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6044 
    6045                 gauss->GaussPoint(ig);
    6046                 gauss->SynchronizeGaussBase(gauss_tria);
    6047 
    6048                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6049                 tria->GetBSSA(&B[0][0], &xyz_list[0][0], gauss_tria);
    6050                 tria->GetBprimeSSA(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    6051 
    6052                 if(approximation==SSAHOApproximationEnum){
    6053                         this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    6054                         this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    6055                         material->GetViscosity3d(&viscosity, &epsilon[0]);
    6056                         material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    6057 
    6058                         newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    6059                 }
    6060                 else if (approximation==SSAFSApproximationEnum){
    6061                         this->StrainRateFS(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    6062                         material->GetViscosity3dFS(&newviscosity,&epsilons[0]);
    6063                 }
    6064                 else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet");
    6065 
    6066                 D_scalar=2*newviscosity*gauss->weight*Jdet;
    6067                 for (i=0;i<3;i++) D[i][i]=D_scalar;
    6068 
    6069                 TripleMultiply( &B[0][0],3,numdof2d,1,
    6070                                         &D[0][0],3,3,0,
    6071                                         &Bprime[0][0],3,numdof2d,0,
    6072                                         &Ke_gg_gaussian[0][0],0);
    6073 
    6074                 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg_gaussian[i][j];
    6075         }
    6076 
    6077         /*Transform Coordinate System*/
    6078         ::TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,XYEnum);
    6079 
    6080         /*Clean up and return*/
    6081         delete tria->material;
    6082         delete tria;
    6083         delete gauss_tria;
    6084         delete gauss;
    6085         return Ke;
    6086 }
    6087 /*}}}*/
    6088 /*FUNCTION Penta::CreateKMatrixStressbalanceSSA3dFriction{{{*/
    6089 ElementMatrix* Penta::CreateKMatrixStressbalanceSSA3dFriction(void){
    6090 
    6091         /*Initialize Element matrix and return if necessary*/
    6092         if(IsFloating() || !IsOnBed()) return NULL;
    6093 
    6094         /*Build a tria element using the 3 nodes of the base of the penta. Then use
    6095          * the tria functionality to build a friction stiffness matrix on these 3
    6096          * nodes: */
    6097         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    6098         ElementMatrix* Ke=tria->CreateKMatrixStressbalanceSSAFriction();
    6099         delete tria->material; delete tria;
    6100 
    6101         /*clean-up and return*/
    6102         return Ke;
    6103 }
    6104 /*}}}*/
    6105 /*FUNCTION Penta::CreateKMatrixStressbalanceSSAHO{{{*/
    6106 ElementMatrix* Penta::CreateKMatrixStressbalanceSSAHO(void){
    6107 
    6108         /*compute all stiffness matrices for this element*/
    6109         ElementMatrix* Ke1=CreateKMatrixStressbalanceSSA3d();
    6110         ElementMatrix* Ke2=CreateKMatrixStressbalanceHO();
    6111         ElementMatrix* Ke3=CreateKMatrixCouplingSSAHO();
    6112         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    6113 
    6114         /*clean-up and return*/
    6115         delete Ke1;
    6116         delete Ke2;
    6117         delete Ke3;
    6118         return Ke;
    6119 }
    6120 /*}}}*/
    6121 /*FUNCTION Penta::CreateKMatrixStressbalanceSSAFS{{{*/
    6122 ElementMatrix* Penta::CreateKMatrixStressbalanceSSAFS(void){
    6123 
    6124         /*compute all stiffness matrices for this element*/
    6125         ElementMatrix* Ke1=CreateKMatrixStressbalanceFS();
    6126         int indices[3]={18,19,20};
    6127         Ke1->StaticCondensation(3,&indices[0]);
    6128         int init = this->element_type;
    6129         this->element_type=P1Enum; //P1 needed for HO
    6130         ElementMatrix* Ke2=CreateKMatrixStressbalanceSSA3d();
    6131         this->element_type=init;
    6132         ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS();
    6133         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    6134 
    6135         /*clean-up and return*/
    6136         delete Ke1;
    6137         delete Ke2;
    6138         delete Ke3;
    6139         return Ke;
    6140 }
    6141 /*}}}*/
    6142 /*FUNCTION Penta::CreateKMatrixStressbalanceL1L2{{{*/
    6143 ElementMatrix* Penta::CreateKMatrixStressbalanceL1L2(void){
    6144 
    6145         /*compute all stiffness matrices for this element*/
    6146         ElementMatrix* Ke1=CreateKMatrixStressbalanceL1L2Viscous();
    6147         ElementMatrix* Ke2=CreateKMatrixStressbalanceL1L2Friction();
    6148         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    6149 
    6150         /*clean-up and return*/
    6151         delete Ke1;
    6152         delete Ke2;
    6153         return Ke;
    6154 }
    6155 /*}}}*/
    6156 /*FUNCTION Penta::CreateKMatrixStressbalanceL1L2Viscous{{{*/
    6157 ElementMatrix* Penta::CreateKMatrixStressbalanceL1L2Viscous(void){
    6158 
    6159         /*Constants*/
    6160         const int    numdof2d=2*NUMVERTICES2D;
    6161 
    6162         /*Intermediaries */
    6163         int         i,j;
    6164         IssmDouble  Jdet,viscosity;
    6165         IssmDouble  xyz_list[NUMVERTICES][3];
    6166         IssmDouble  B[3][numdof2d];
    6167         IssmDouble  Bprime[3][numdof2d];
    6168         IssmDouble  Ke_gg_gaussian[numdof2d][numdof2d];
    6169         IssmDouble  D[3][3]= {0.0};                 // material matrix, simple scalar matrix.
    6170         Tria       *tria       = NULL;
    6171         Penta      *pentabase  = NULL;
    6172         GaussPenta *gauss      = NULL;
    6173         GaussTria  *gauss_tria = NULL;
    6174 
    6175         /*Find penta on bed as this is a SSA elements: */
    6176         pentabase=(Penta*)GetBasalElement();
    6177         tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
    6178 
    6179         /*Initialize Element matrix*/
    6180         ElementMatrix* Ke=new ElementMatrix(tria->nodes,NUMVERTICES2D,this->parameters,L1L2ApproximationEnum);
    6181 
    6182         /*Retrieve all inputs and parameters*/
    6183         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6184         Input* vx_input=inputs->GetInput(VxEnum);        _assert_(vx_input);
    6185         Input* vy_input=inputs->GetInput(VyEnum);        _assert_(vy_input);
    6186         Input* surf_input=inputs->GetInput(SurfaceEnum); _assert_(surf_input);
    6187 
    6188         /* Start  looping on the number of gaussian points: */
    6189         gauss=new GaussPenta(5,5);
    6190         gauss_tria=new GaussTria();
    6191         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6192 
    6193                 gauss->GaussPoint(ig);
    6194                 gauss->SynchronizeGaussBase(gauss_tria);
    6195 
    6196                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6197                 tria->GetBSSA(&B[0][0], &xyz_list[0][0], gauss_tria);
    6198                 tria->GetBprimeSSA(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    6199 
    6200                 /*Get viscosity for L1L2 model*/
    6201                 ViscosityL1L2(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input,surf_input);
    6202 
    6203                 for(i=0;i<3;i++) D[i][i]=2*viscosity*gauss->weight*Jdet;
    6204 
    6205                 TripleMultiply( &B[0][0],3,numdof2d,1,
    6206                                         &D[0][0],3,3,0,
    6207                                         &Bprime[0][0],3,numdof2d,0,
    6208                                         &Ke_gg_gaussian[0][0],0);
    6209                 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg_gaussian[i][j];
    6210         }
    6211 
    6212         /*Transform Coordinate System*/
    6213         ::TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,XYEnum);
    6214 
    6215         /*Clean up and return*/
    6216         delete tria->material;
    6217         delete tria;
    6218         delete gauss_tria;
    6219         delete gauss;
    6220         return Ke;
    6221 }
    6222 /*}}}*/
    6223 /*FUNCTION Penta::CreateKMatrixStressbalanceL1L2Friction{{{*/
    6224 ElementMatrix* Penta::CreateKMatrixStressbalanceL1L2Friction(void){
    6225 
    6226         /*Initialize Element matrix and return if necessary*/
    6227         if(IsFloating() || !IsOnBed()) return NULL;
    6228 
    6229         /*Build a tria element using the 3 nodes of the base of the penta. Then use
    6230          * the tria functionality to build a friction stiffness matrix on these 3
    6231          * nodes: */
    6232         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    6233         ElementMatrix* Ke=tria->CreateKMatrixStressbalanceSSAFriction();
    6234         delete tria->material; delete tria;
    6235 
    6236         /*clean-up and return*/
    6237         return Ke;
    6238 }
    6239 /*}}}*/
    6240 /*FUNCTION Penta::CreateKMatrixStressbalanceHO{{{*/
    6241 ElementMatrix* Penta::CreateKMatrixStressbalanceHO(void){
    6242 
    6243         /*compute all stiffness matrices for this element*/
    6244         ElementMatrix* Ke1=CreateKMatrixStressbalanceHOViscous();
    6245         ElementMatrix* Ke2=CreateKMatrixStressbalanceHOFriction();
    6246         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    6247 
    6248         /*clean-up and return*/
    6249         delete Ke1;
    6250         delete Ke2;
    6251         return Ke;
    6252 
    6253 }
    6254 /*}}}*/
    6255 /*FUNCTION Penta::CreateKMatrixStressbalanceHOViscous{{{*/
    6256 ElementMatrix* Penta::CreateKMatrixStressbalanceHOViscous(void){
    6257 
    6258         /*Intermediaries */
    6259         int         approximation;
    6260         IssmDouble  xyz_list[NUMVERTICES][3];
    6261         IssmDouble  Jdet;
    6262         IssmDouble  viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
    6263         IssmDouble  epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    6264         IssmDouble  D_scalar;
    6265         GaussPenta *gauss=NULL;
    6266 
    6267         /*Fetch number of nodes and dof for this finite element*/
    6268         int numnodes = this->NumberofNodes();
    6269         int numdof   = numnodes*NDOF2;
    6270 
    6271         /*Initialize Element matrix and vectors*/
    6272         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,HOApproximationEnum);
    6273         IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
    6274         IssmDouble*    Bprime = xNew<IssmDouble>(5*numdof);
    6275         IssmDouble*    D      = xNewZeroInit<IssmDouble>(5*5);
    6276 
    6277         /*Retrieve all inputs and parameters*/
    6278         inputs->GetInputValue(&approximation,ApproximationEnum);
    6279         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6280         this->parameters->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
    6281         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    6282         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    6283         Input* vxold_input=inputs->GetInput(VxPicardEnum); _assert_(vxold_input);
    6284         Input* vyold_input=inputs->GetInput(VyPicardEnum); _assert_(vyold_input);
    6285 
    6286         /* Start  looping on the number of gaussian points: */
    6287         gauss=new GaussPenta(5,5);
    6288         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6289 
    6290                 gauss->GaussPoint(ig);
    6291 
    6292                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6293                 GetBHO(&B[0], &xyz_list[0][0], gauss);
    6294                 GetBprimeHO(&Bprime[0], &xyz_list[0][0], gauss);
    6295 
    6296                 this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    6297                 this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    6298                 material->GetViscosity3d(&viscosity, &epsilon[0]);
    6299                 material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    6300                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    6301 
    6302                 D_scalar=2*newviscosity*gauss->weight*Jdet;
    6303                 for(int i=0;i<5;i++) D[i*5+i]=D_scalar;
    6304 
    6305                 TripleMultiply(B,5,numdof,1,
    6306                                         D,5,5,0,
    6307                                         Bprime,5,numdof,0,
    6308                                         &Ke->values[0],1);
    6309         }
    6310 
    6311         /*Transform Coordinate System*/
    6312         ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);
    6313 
    6314         /*Clean up and return*/
    6315         delete gauss;
    6316         xDelete<IssmDouble>(D);
    6317         xDelete<IssmDouble>(Bprime);
    6318         xDelete<IssmDouble>(B);
    6319         return Ke;
    6320 }
    6321 /*}}}*/
    6322 /*FUNCTION Penta::CreateKMatrixStressbalanceHOFriction{{{*/
    6323 ElementMatrix* Penta::CreateKMatrixStressbalanceHOFriction(void){
    6324 
    6325         /*Intermediaries */
    6326         bool       mainlyfloating;
    6327         int         i,j;
    6328         int         analysis_type,migration_style;
    6329         int         point1;
    6330         IssmDouble  xyz_list[NUMVERTICES][3];
    6331         IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
    6332         IssmDouble  alpha2,Jdet;
    6333         IssmDouble fraction1,fraction2;
    6334         IssmDouble gllevelset;
    6335         IssmDouble  phi=1.0;
    6336         IssmDouble  DL_scalar;
    6337         Friction   *friction = NULL;
    6338         GaussPenta *gauss    = NULL;
    6339 
    6340         /*Initialize Element matrix and return if necessary*/
    6341         if(IsFloating() || !IsOnBed()) return NULL;
    6342 
    6343         /*Fetch number of nodes and dof for this finite element*/
    6344         int numnodes = this->NumberofNodes();
    6345         int numdof   = numnodes*NDOF2;
    6346 
    6347         /*Initialize Element matrix and vectors*/
    6348         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,HOApproximationEnum);
    6349         IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
    6350         IssmDouble*    D      = xNewZeroInit<IssmDouble>(2*2);
    6351 
    6352         /*Retrieve all inputs and parameters*/
    6353         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6354         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    6355         parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
    6356         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    6357         Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    6358         Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
    6359         Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
    6360         Input* vz_input=inputs->GetInput(VzEnum);           _assert_(vz_input);
    6361         Input* gllevelset_input=NULL;
    6362 
    6363         /*build friction object, used later on: */
    6364         friction=new Friction(this,2);
    6365 
    6366         /*Recover portion of element that is grounded*/
    6367         if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list_tria[0][0]);
    6368         if(migration_style==SubelementMigration2Enum){
    6369                 gllevelset_input=inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
    6370                 this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    6371                 //gauss=new GaussPenta(point1,fraction1,fraction2,mainlyfloating,2);
    6372                 gauss=new GaussPenta(0,1,2,2);
    6373         }
    6374         else{
    6375                 gauss=new GaussPenta(0,1,2,2);
    6376         }
    6377 
    6378         /* Start  looping on the number of gaussian points: */
    6379         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6380 
    6381                 gauss->GaussPoint(ig);
    6382 
    6383                 GetTriaJacobianDeterminant(&Jdet,&xyz_list_tria[0][0],gauss);
    6384                 GetBHOFriction(&B[0],gauss);
    6385 
    6386                 friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    6387                 if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
    6388                 if(migration_style==SubelementMigration2Enum){
    6389                         gllevelset_input->GetInputValue(&gllevelset, gauss);
    6390                         if(gllevelset<0) alpha2=0;
    6391                 }
    6392 
    6393                 DL_scalar=alpha2*gauss->weight*Jdet;
    6394                 for (i=0;i<2;i++) D[i*2+i]=DL_scalar;
    6395 
    6396                 TripleMultiply(B,2,numdof,1,
    6397                                         D,2,2,0,
    6398                                         B,2,numdof,0,
    6399                                         &Ke->values[0],1);
    6400         }
    6401 
    6402         /*Transform Coordinate System*/
    6403         ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);
    6404 
    6405         /*Clean up and return*/
    6406         delete gauss;
    6407         xDelete<IssmDouble>(D);
    6408         xDelete<IssmDouble>(B);
    6409         delete friction;
    6410         return Ke;
    6411 }
    6412 /*}}}*/
    6413 /*FUNCTION Penta::CreateKMatrixStressbalanceHOFS{{{*/
    6414 ElementMatrix* Penta::CreateKMatrixStressbalanceHOFS(void){
    6415 
    6416         /*compute all stiffness matrices for this element*/
    6417         ElementMatrix* Ke1=CreateKMatrixStressbalanceFS();
    6418         int indices[3]={18,19,20};
    6419         Ke1->StaticCondensation(3,&indices[0]);
    6420         int init = this->element_type;
    6421         this->element_type=P1Enum; //P1 needed for HO
    6422         ElementMatrix* Ke2=CreateKMatrixStressbalanceHO();
    6423         this->element_type=init;
    6424         ElementMatrix* Ke3=CreateKMatrixCouplingHOFS();
    6425         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    6426 
    6427         /*clean-up and return*/
    6428         delete Ke1;
    6429         delete Ke2;
    6430         delete Ke3;
    6431         return Ke;
    6432 }
    6433 /*}}}*/
    6434 /*FUNCTION Penta::CreateKMatrixStressbalanceFS{{{*/
    6435 ElementMatrix* Penta::CreateKMatrixStressbalanceFS(void){
    6436 
    6437         ElementMatrix* Ke1 = NULL;
    6438         ElementMatrix* Ke2 = NULL;
    6439         ElementMatrix* Ke  = NULL;
    6440 
    6441         /*compute all stiffness matrices for this element*/
    6442         Ke1=CreateKMatrixStressbalanceFSViscous();
    6443         Ke2=CreateKMatrixStressbalanceFSFriction();
    6444         Ke =new ElementMatrix(Ke1,Ke2);
    6445 
    6446         /*clean-up and return*/
    6447         delete Ke1;
    6448         delete Ke2;
    6449         return Ke;
    6450 
    6451 }
    6452 /*}}}*/
    64534522/*FUNCTION Penta::KMatrixGLSstabilization{{{*/
    64544523void Penta::KMatrixGLSstabilization(ElementMatrix* Ke){
    … …  
    65564625}
    65574626/*}}}*/
    6558 /*FUNCTION Penta::CreateKMatrixStressbalanceFSViscous {{{*/
    6559 ElementMatrix* Penta::CreateKMatrixStressbalanceFSViscous(void){
    6560 
    6561         /*Intermediaries */
    6562         int        i,approximation;
    6563         IssmDouble Jdet,viscosity,FSreconditioning,D_scalar;
    6564         IssmDouble xyz_list[NUMVERTICES][3];
    6565         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    6566         GaussPenta *gauss=NULL;
    6567 
    6568         /*If on water or not FS, skip stiffness: */
    6569         inputs->GetInputValue(&approximation,ApproximationEnum);
    6570         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    6571 
    6572         /*Fetch number of nodes and dof for this finite element*/
    6573         int vnumnodes = this->NumberofNodesVelocity();
    6574         int pnumnodes = this->NumberofNodesPressure();
    6575         int numdof    = vnumnodes*NDOF3 + pnumnodes*NDOF1;
    6576 
    6577         /*Prepare coordinate system list*/
    6578         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    6579         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    6580         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    6581 
    6582         /*Initialize Element matrix and vectors*/
    6583         ElementMatrix* Ke     = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    6584         IssmDouble*    B      = xNew<IssmDouble>(8*numdof);
    6585         IssmDouble*    Bprime = xNew<IssmDouble>(8*numdof);
    6586         IssmDouble*    D      = xNewZeroInit<IssmDouble>(8*8);
    6587 
    6588         /*Retrieve all inputs and parameters*/
    6589         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6590         parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    6591         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    6592         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    6593         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    6594 
    6595         /* Start  looping on the number of gaussian points: */
    6596         gauss=new GaussPenta(5,5);
    6597         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6598 
    6599                 gauss->GaussPoint(ig);
    6600 
    6601                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6602                 GetBFS(B,&xyz_list[0][0],gauss);
    6603                 GetBprimeFS(Bprime,&xyz_list[0][0],gauss);
    6604 
    6605                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    6606                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    6607 
    6608                 D_scalar=gauss->weight*Jdet;
    6609                 for(i=0;i<6;i++) D[i*8+i] = +D_scalar*2.*viscosity;
    6610                 for(i=6;i<8;i++) D[i*8+i] = -D_scalar*FSreconditioning;
    6611 
    6612                 TripleMultiply(B,8,numdof,1,
    6613                                         D,8,8,0,
    6614                                         Bprime,8,numdof,0,
    6615                                         Ke->values,1);
    6616         }
    6617 
    6618         /*Transform Coordinate System*/
    6619         ::TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);
    6620 
    6621         /*Clean up and return*/
    6622         delete gauss;
    6623         xDelete<int>(cs_list);
    6624         xDelete<IssmDouble>(B);
    6625         xDelete<IssmDouble>(Bprime);
    6626         xDelete<IssmDouble>(D);
    6627         return Ke;
    6628 }
    6629 /*}}}*/
    6630 /*FUNCTION Penta::CreateKMatrixStressbalanceFSFriction{{{*/
    6631 ElementMatrix* Penta::CreateKMatrixStressbalanceFSFriction(void){
    6632 
    6633         /*Intermediaries */
    6634         int         i,j;
    6635         int         analysis_type,approximation;
    6636         IssmDouble  alpha2,Jdet2d;
    6637         IssmDouble  xyz_list[NUMVERTICES][3];
    6638         IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
    6639         Friction   *friction = NULL;
    6640         GaussPenta *gauss    = NULL;
    6641 
    6642         /*Initialize Element matrix and return if necessary*/
    6643         if(IsFloating() || !IsOnBed()) return NULL;
    6644 
    6645         /*If on water or not FS, skip stiffness: */
    6646         inputs->GetInputValue(&approximation,ApproximationEnum);
    6647         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    6648 
    6649         /*Fetch number of nodes and dof for this finite element*/
    6650         int vnumnodes = this->NumberofNodesVelocity();
    6651         int pnumnodes = this->NumberofNodesPressure();
    6652         int  numdof   = vnumnodes*NDOF3 + pnumnodes*NDOF1;
    6653 
    6654         /*Initialize Element matrix and vectors*/
    6655         ElementMatrix* Ke        = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    6656         IssmDouble*    BFriction = xNew<IssmDouble>(2*numdof);
    6657         IssmDouble*    D         = xNewZeroInit<IssmDouble>(2*2);
    6658 
    6659         /*Retrieve all inputs and parameters*/
    6660         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6661         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    6662         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    6663         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    6664         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    6665         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    6666 
    6667         /*build friction object, used later on: */
    6668         friction=new Friction(this,3);
    6669 
    6670         /* Start  looping on the number of gaussian points: */
    6671         gauss=new GaussPenta(0,1,2,3);
    6672         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6673 
    6674                 gauss->GaussPoint(ig);
    6675 
    6676                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    6677                 GetLFS(BFriction,gauss);
    6678 
    6679                 friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    6680 
    6681                 D[0*2+0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
    6682                 D[1*2+1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy
    6683 
    6684                 TripleMultiply(BFriction,2,numdof,1,
    6685                                         D,2,2,0,
    6686                                         BFriction,2,numdof,0,
    6687                                         Ke->values,1);
    6688         }
    6689 
    6690         /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
    6691         //::TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZEnum);
    6692 
    6693         /*Clean up and return*/
    6694         delete gauss;
    6695         delete friction;
    6696         xDelete<IssmDouble>(BFriction);
    6697         xDelete<IssmDouble>(D);
    6698         return Ke;
    6699 }
    6700 /*}}}*/
    6701 /*FUNCTION Penta::CreateKMatrixStressbalanceVert {{{*/
    6702 ElementMatrix* Penta::CreateKMatrixStressbalanceVert(void){
    6703 
    6704         /*compute all stiffness matrices for this element*/
    6705         ElementMatrix* Ke1=CreateKMatrixStressbalanceVertVolume();
    6706         ElementMatrix* Ke2=CreateKMatrixStressbalanceVertSurface();
    6707         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    6708 
    6709         /*clean-up and return*/
    6710         delete Ke1;
    6711         delete Ke2;
    6712         return Ke;
    6713 
    6714 }
    6715 /*}}}*/
    6716 /*FUNCTION Penta::CreateKMatrixStressbalanceVertVolume {{{*/
    6717 ElementMatrix* Penta::CreateKMatrixStressbalanceVertVolume(void){
    6718 
    6719         /*Intermediaries */
    6720         IssmDouble  Jdet;
    6721         IssmDouble  xyz_list[NUMVERTICES][3];
    6722         IssmDouble  B[NDOF1][NUMVERTICES];
    6723         IssmDouble  Bprime[NDOF1][NUMVERTICES];
    6724         IssmDouble  DL_scalar;
    6725         GaussPenta  *gauss=NULL;
    6726 
    6727         /*Initialize Element matrix*/
    6728         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    6729 
    6730         /*Retrieve all inputs and parameters*/
    6731         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6732 
    6733         /* Start  looping on the number of gaussian points: */
    6734         gauss=new GaussPenta(2,2);
    6735         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6736 
    6737                 gauss->GaussPoint(ig);
    6738 
    6739                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6740                 GetBVert(&B[0][0], &xyz_list[0][0], gauss);
    6741                 GetBprimeVert(&Bprime[0][0], &xyz_list[0][0], gauss);
    6742 
    6743                 DL_scalar=gauss->weight*Jdet;
    6744 
    6745                 TripleMultiply( &B[0][0],1,NUMVERTICES,1,
    6746                                         &DL_scalar,1,1,0,
    6747                                         &Bprime[0][0],1,NUMVERTICES,0,
    6748                                         &Ke->values[0],1);
    6749         }
    6750 
    6751         /*Clean up and return*/
    6752         delete gauss;
    6753         return Ke;
    6754 }
    6755 /*}}}*/
    6756 /*FUNCTION Penta::CreateKMatrixStressbalanceVertSurface {{{*/
    6757 ElementMatrix* Penta::CreateKMatrixStressbalanceVertSurface(void){
    6758 
    6759         if (!IsOnSurface()) return NULL;
    6760 
    6761         /*Constants*/
    6762         const int numdof=NDOF1*NUMVERTICES;
    6763 
    6764         /*Intermediaries */
    6765         int       i,j;
    6766         IssmDouble xyz_list[NUMVERTICES][3];
    6767         IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    6768         IssmDouble surface_normal[3];
    6769         IssmDouble Jdet2d,DL_scalar;
    6770         IssmDouble basis[NUMVERTICES];
    6771         GaussPenta *gauss=NULL;
    6772 
    6773         /*Initialize Element matrix*/
    6774         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    6775 
    6776         /*Retrieve all inputs and parameters*/
    6777         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6778         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i+3][j];
    6779         NormalTop(&surface_normal[0],&xyz_list_tria[0][0]);
    6780 
    6781         /* Start  looping on the number of gaussian points: */
    6782         gauss=new GaussPenta(3,4,5,2);
    6783         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6784 
    6785                 gauss->GaussPoint(ig);
    6786 
    6787                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    6788                 GetNodalFunctionsP1(&basis[0], gauss);
    6789 
    6790                 DL_scalar= - gauss->weight*Jdet2d*surface_normal[2];
    6791 
    6792                 TripleMultiply( basis,1,numdof,1,
    6793                                         &DL_scalar,1,1,0,
    6794                                         basis,1,numdof,0,
    6795                                         &Ke->values[0],1);
    6796         }
    6797 
    6798         /*Clean up and return*/
    6799         delete gauss;
    6800         return Ke;
    6801 }
    6802 /*}}}*/
    6803 /*FUNCTION Penta::CreateJacobianStressbalanceHoriz{{{*/
    6804 ElementMatrix* Penta::CreateJacobianStressbalanceHoriz(void){
    6805 
    6806         int approximation;
    6807         inputs->GetInputValue(&approximation,ApproximationEnum);
    6808 
    6809         switch(approximation){
    6810                 case SSAApproximationEnum:
    6811                         return CreateJacobianStressbalanceSSA2d();
    6812                 case HOApproximationEnum:
    6813                         return CreateJacobianStressbalanceHO();
    6814                 case FSApproximationEnum:
    6815                         return CreateJacobianStressbalanceFS();
    6816                 case NoneApproximationEnum:
    6817                         return NULL;
    6818                 default:
    6819                         _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
    6820         }
    6821 }
    6822 /*}}}*/
    6823 /*FUNCTION Penta::CreateJacobianStressbalanceSSA2d{{{*/
    6824 ElementMatrix* Penta::CreateJacobianStressbalanceSSA2d(void){
    6825 
    6826         /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the
    6827           bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build
    6828           the stiffness matrix. */
    6829         if (!IsOnBed()) return NULL;
    6830 
    6831         /*Depth average some fields*/
    6832         switch(this->material->ObjectEnum()){
    6833                 case MaticeEnum:
    6834                         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    6835                         this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
    6836                         break;
    6837                 default:
    6838                         _error_("material "<<EnumToStringx(this->material->ObjectEnum())<<" not supported");
    6839         }
    6840 
    6841         /*Call Tria function*/
    6842         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    6843         ElementMatrix* Ke=tria->CreateJacobianStressbalanceSSA();
    6844         delete tria->material; delete tria;
    6845 
    6846         /*Delete averaged inputs*/
    6847         this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    6848         this->material->inputs->DeleteInput(DamageDbarEnum);
    6849 
    6850         /*clean up and return*/
    6851         return Ke;
    6852 }
    6853 /*}}}*/
    6854 /*FUNCTION Penta::CreateJacobianStressbalanceHO{{{*/
    6855 ElementMatrix* Penta::CreateJacobianStressbalanceHO(void){
    6856 
    6857         /*Intermediaries */
    6858         int        i,j;
    6859         IssmDouble xyz_list[NUMVERTICES][3];
    6860         IssmDouble Jdet;
    6861         IssmDouble eps1dotdphii,eps1dotdphij;
    6862         IssmDouble eps2dotdphii,eps2dotdphij;
    6863         IssmDouble mu_prime;
    6864         IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    6865         IssmDouble eps1[3],eps2[3];
    6866         IssmDouble dphi[3][NUMVERTICES];
    6867         GaussPenta *gauss=NULL;
    6868 
    6869         /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
    6870         ElementMatrix* Ke=CreateKMatrixStressbalanceHO();
    6871 
    6872         /*Retrieve all inputs and parameters*/
    6873         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6874         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    6875         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    6876 
    6877         /* Start  looping on the number of gaussian points: */
    6878         gauss=new GaussPenta(5,5);
    6879         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6880 
    6881                 gauss->GaussPoint(ig);
    6882 
    6883                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6884                 GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    6885 
    6886                 this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    6887                 material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    6888                 eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    6889                 eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
    6890                 eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
    6891 
    6892                 for(i=0;i<6;i++){
    6893                         for(j=0;j<6;j++){
    6894                                 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
    6895                                 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
    6896                                 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
    6897                                 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
    6898 
    6899                                 Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
    6900                                 Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
    6901                                 Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
    6902                                 Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
    6903                         }
    6904                 }
    6905         }
    6906 
    6907         /*Transform Coordinate System*/
    6908         ::TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
    6909 
    6910         /*Clean up and return*/
    6911         delete gauss;
    6912         return Ke;
    6913 }
    6914 /*}}}*/
    6915 /*FUNCTION Penta::CreateJacobianStressbalanceFS{{{*/
    6916 ElementMatrix* Penta::CreateJacobianStressbalanceFS(void){
    6917 
    6918         /*Intermediaries */
    6919         int        i,j;
    6920         IssmDouble xyz_list[NUMVERTICES][3];
    6921         IssmDouble Jdet;
    6922         IssmDouble eps1dotdphii,eps1dotdphij;
    6923         IssmDouble eps2dotdphii,eps2dotdphij;
    6924         IssmDouble eps3dotdphii,eps3dotdphij;
    6925         IssmDouble mu_prime;
    6926         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    6927         IssmDouble eps1[3],eps2[3],eps3[3];
    6928         GaussPenta *gauss=NULL;
    6929 
    6930         /*If on water or not FS, skip stiffness: */
    6931         //inputs->GetInputValue(&approximation,ApproximationEnum);
    6932         //if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    6933 
    6934         /*Fetch number of nodes and dof for this finite element*/
    6935         int vnumnodes = this->NumberofNodesVelocity();
    6936         int pnumnodes = this->NumberofNodesPressure();
    6937         int numdof    = vnumnodes*NDOF3 + pnumnodes*NDOF1;
    6938 
    6939         /*Prepare coordinate system list*/
    6940         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    6941         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    6942         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    6943 
    6944         /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
    6945         ElementMatrix* Ke=CreateKMatrixStressbalanceFS();
    6946         IssmDouble*    dbasis = xNew<IssmDouble>(3*vnumnodes);
    6947 
    6948         /*Retrieve all inputs and parameters*/
    6949         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6950         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    6951         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    6952         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    6953 
    6954         /* Start  looping on the number of gaussian points: */
    6955         gauss=new GaussPenta(5,5);
    6956         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6957 
    6958                 gauss->GaussPoint(ig);
    6959 
    6960                 GetJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    6961                 GetNodalFunctionsDerivativesVelocity(dbasis,&xyz_list[0][0],gauss);
    6962 
    6963                 //this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    6964                 //material->GetViscosityDerivativeEpsSquareFS(&mu_prime,&epsilon[0]);
    6965                 //eps1[0]=epsilon[0];   eps2[0]=epsilon[3];   eps3[0]=epsilon[4];
    6966                 //eps1[1]=epsilon[3];   eps2[1]=epsilon[1];   eps3[1]=epsilon[5];
    6967                 //eps1[2]=epsilon[4];   eps2[2]=epsilon[5];   eps3[2]=epsilon[2];
    6968                 this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    6969                 material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    6970                 eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
    6971                 eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
    6972                 eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
    6973 
    6974                 for(i=0;i<vnumnodes-1;i++){
    6975                         for(j=0;j<vnumnodes-1;j++){
    6976                                 eps1dotdphii=eps1[0]*dbasis[0*vnumnodes+i]+eps1[1]*dbasis[1*vnumnodes+i]+eps1[2]*dbasis[2*vnumnodes+i];
    6977                                 eps1dotdphij=eps1[0]*dbasis[0*vnumnodes+j]+eps1[1]*dbasis[1*vnumnodes+j]+eps1[2]*dbasis[2*vnumnodes+j];
    6978                                 eps2dotdphii=eps2[0]*dbasis[0*vnumnodes+i]+eps2[1]*dbasis[1*vnumnodes+i]+eps2[2]*dbasis[2*vnumnodes+i];
    6979                                 eps2dotdphij=eps2[0]*dbasis[0*vnumnodes+j]+eps2[1]*dbasis[1*vnumnodes+j]+eps2[2]*dbasis[2*vnumnodes+j];
    6980                                 eps3dotdphii=eps3[0]*dbasis[0*vnumnodes+i]+eps3[1]*dbasis[1*vnumnodes+i]+eps3[2]*dbasis[2*vnumnodes+i];
    6981                                 eps3dotdphij=eps3[0]*dbasis[0*vnumnodes+j]+eps3[1]*dbasis[1*vnumnodes+j]+eps3[2]*dbasis[2*vnumnodes+j];
    6982 
    6983                                 Ke->values[numdof*(3*i+0)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
    6984                                 Ke->values[numdof*(3*i+0)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
    6985                                 Ke->values[numdof*(3*i+0)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
    6986 
    6987                                 Ke->values[numdof*(3*i+1)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
    6988                                 Ke->values[numdof*(3*i+1)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
    6989                                 Ke->values[numdof*(3*i+1)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
    6990 
    6991                                 Ke->values[numdof*(3*i+2)+3*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
    6992                                 Ke->values[numdof*(3*i+2)+3*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
    6993                                 Ke->values[numdof*(3*i+2)+3*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
    6994                         }
    6995                 }
    6996         }
    6997 
    6998         /*Transform Coordinate System*/
    6999         ::TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);
    7000 
    7001         /*Clean up and return*/
    7002         xDelete<int>(cs_list);
    7003         xDelete<IssmDouble>(dbasis);
    7004         delete gauss;
    7005         return Ke;
    7006 }
    7007 /*}}}*/
    70084627#endif
    70094628
    7010 #ifdef _HAVE_BALANCED_
    7011 /*FUNCTION Penta::CreateKMatrixBalancethickness {{{*/
    7012 ElementMatrix* Penta::CreateKMatrixBalancethickness(void){
    7013 
    7014         /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the
    7015           bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build
    7016           the stiffness matrix. */
    7017         if (!IsOnBed()) return NULL;
    7018 
    7019         /*Depth Averaging Vx and Vy*/
    7020         this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
    7021         this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
    7022 
    7023         /*Spawn Tria element from the base of the Penta: */
    7024         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    7025         ElementMatrix* Ke=tria->CreateKMatrixBalancethickness();
    7026         delete tria->material; delete tria;
    7027 
    7028         /*Delete Vx and Vy averaged*/
    7029         this->inputs->DeleteInput(VxAverageEnum);
    7030         this->inputs->DeleteInput(VyAverageEnum);
    7031 
    7032         /*clean up and return*/
    7033         return Ke;
    7034 }
    7035 /*}}}*/
    7036 #endif
    7037 
    70384629#ifdef _HAVE_HYDROLOGY_
    7039 /*FUNCTION Penta::CreateKMatrixHydrologyDCInefficient {{{*/
    7040 ElementMatrix* Penta::CreateKMatrixHydrologyDCInefficient(void){
    7041 
    7042         if (!IsOnBed()) return NULL;
    7043 
    7044         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    7045         ElementMatrix* Ke=tria->CreateKMatrixHydrologyDCInefficient();
    7046         delete tria->material; delete tria;
    7047 
    7048         /*clean up and return*/
    7049         return Ke;
    7050 }
    7051 /*}}}*/
    7052 /*FUNCTION Penta::CreateKMatrixHydrologyDCEfficient {{{*/
    7053 ElementMatrix* Penta::CreateKMatrixHydrologyDCEfficient(void){
    7054 
    7055         if (!IsOnBed()) return NULL;
    7056 
    7057         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    7058         ElementMatrix* Ke=tria->CreateKMatrixHydrologyDCEfficient();
    7059         delete tria->material; delete tria;
    7060 
    7061         /*clean up and return*/
    7062         return Ke;
    7063 }
    7064 /*}}}*/
    70654630/*FUNCTION Penta::CreateEPLDomainMassMatrix {{{*/
    70664631ElementMatrix* Penta::CreateEPLDomainMassMatrix(void){

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

    r16982 r16993  
    7272                void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
    7373                void   CreateDVector(Vector<IssmDouble>* df);
    74                 void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
    7574                void   Delta18oParameterization(void);
    7675                void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
    … …  
    195194                void             NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list);
    196195                ElementMatrix* CreateBasalMassMatrix(void);
    197                 ElementMatrix* CreateKMatrix(void);
    198                 ElementMatrix* CreateKMatrixMasstransport(void);
    199                 ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
    200                 ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
    201196                IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
    202197                IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
    … …  
    260255
    261256                #ifdef _HAVE_STRESSBALANCE_
    262                 ElementMatrix* CreateKMatrixCouplingSSAHO(void);
    263                 ElementMatrix* CreateKMatrixCouplingSSAHOViscous(void);
    264                 ElementMatrix* CreateKMatrixCouplingSSAHOFriction(void);
    265                 ElementMatrix* CreateKMatrixCouplingSSAFS(void);
    266                 ElementMatrix* CreateKMatrixCouplingSSAFSViscous(void);
    267                 ElementMatrix* CreateKMatrixCouplingSSAFSFriction(void);
    268                 ElementMatrix* CreateKMatrixCouplingHOFS(void);
    269                 ElementMatrix* CreateKMatrixStressbalanceHoriz(void);
    270                 ElementMatrix* CreateKMatrixAdjointHoriz(void);
    271257                ElementVector* CreateDVectorStressbalanceHoriz(void);
    272258                ElementVector* CreateDVectorStressbalanceFS(void);
    273                 ElementMatrix* CreateKMatrixStressbalanceSIA(void);
    274                 ElementMatrix* CreateKMatrixStressbalanceSSA2d(void);
    275                 ElementMatrix* CreateKMatrixStressbalanceSSA3d(void);
    276                 ElementMatrix* CreateKMatrixStressbalanceSSA3dViscous(void);
    277                 ElementMatrix* CreateKMatrixStressbalanceSSA3dFriction(void);
    278                 ElementMatrix* CreateKMatrixStressbalanceSSAHO(void);
    279                 ElementMatrix* CreateKMatrixStressbalanceSSAFS(void);
    280                 ElementMatrix* CreateKMatrixStressbalanceL1L2(void);
    281                 ElementMatrix* CreateKMatrixStressbalanceL1L2Viscous(void);
    282                 ElementMatrix* CreateKMatrixStressbalanceL1L2Friction(void);
    283                 ElementMatrix* CreateKMatrixStressbalanceHO(void);
    284                 ElementMatrix* CreateKMatrixStressbalanceHOViscous(void);
    285                 ElementMatrix* CreateKMatrixStressbalanceHOFriction(void);
    286                 ElementMatrix* CreateKMatrixStressbalanceHOFS(void);
    287                 ElementMatrix* CreateKMatrixStressbalanceFS(void);
    288                 ElementMatrix* CreateKMatrixStressbalanceFSViscous(void);
    289259                void           KMatrixGLSstabilization(ElementMatrix* Ke);
    290                 ElementMatrix* CreateKMatrixStressbalanceFSFriction(void);
    291                 ElementMatrix* CreateKMatrixStressbalanceVert(void);
    292                 ElementMatrix* CreateKMatrixStressbalanceVertVolume(void);
    293                 ElementMatrix* CreateKMatrixStressbalanceVertSurface(void);
    294                 ElementMatrix* CreateJacobianStressbalanceHoriz(void);
    295                 ElementMatrix* CreateJacobianStressbalanceSSA2d(void);
    296                 ElementMatrix* CreateJacobianStressbalanceHO(void);
    297                 ElementMatrix* CreateJacobianStressbalanceFS(void);
    298                 #endif
    299 
    300                 #ifdef _HAVE_CONTROL_
    301                 ElementMatrix* CreateKMatrixAdjointSSA2d(void);
    302                 ElementMatrix* CreateKMatrixAdjointHO(void);
    303                 ElementMatrix* CreateKMatrixAdjointFS(void);
    304260                #endif
    305261
    306262                #ifdef _HAVE_HYDROLOGY_
    307                 ElementMatrix* CreateKMatrixHydrologyDCInefficient(void);
    308                 ElementMatrix* CreateKMatrixHydrologyDCEfficient(void);
    309263                ElementMatrix* CreateEPLDomainMassMatrix(void);
    310264                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    … …  
    319273                void           UpdateConstraintsExtrudeFromTop(void){_error_("not implemented yet");};
    320274                #ifdef _HAVE_THERMAL_
    321                 ElementMatrix* CreateKMatrixEnthalpy(void);
    322                 ElementMatrix* CreateKMatrixEnthalpyVolume(void);
    323                 ElementMatrix* CreateKMatrixEnthalpyShelf(void);
    324                 ElementMatrix* CreateKMatrixThermal(void);
    325                 ElementMatrix* CreateKMatrixMelting(void);
    326                 ElementMatrix* CreateKMatrixThermalVolume(void);
    327                 ElementMatrix* CreateKMatrixThermalShelf(void);
    328275                void           UpdateBasalConstraintsEnthalpy(void);
    329276                void           ComputeBasalMeltingrate(void);
    330277                void           DrainWaterfraction(IssmDouble* drainrate_element);
    331278                #endif
    332 
    333                 #ifdef _HAVE_BALANCED_
    334                 ElementMatrix* CreateKMatrixBalancethickness(void);
    335                 #endif
    336279                /*}}}*/
    337280};

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

    r16982 r16993  
    138138}
    139139/*}}}*/
    140 /*FUNCTION Seg::CreateKMatrixFreeSurfaceTop {{{*/
    141 ElementMatrix* Seg::CreateKMatrixFreeSurfaceTop(void){
    142 
    143         /*Intermediaries */
    144         int        stabilization;
    145         IssmDouble Jdet,D_scalar,dt,h;
    146         IssmDouble vx,vel;
    147         IssmDouble xyz_list[NUMVERTICES][3];
    148 
    149         /*Fetch number of nodes for this finite element*/
    150         int numnodes = this->NumberofNodes();
    151 
    152         /*Initialize Element matrix and vectors*/
    153         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    154         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    155         IssmDouble*    B      = xNew<IssmDouble>(1*numnodes);
    156         IssmDouble*    Bprime = xNew<IssmDouble>(1*numnodes);
    157 
    158         /*Retrieve all inputs and parameters*/
    159         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    160         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    161         this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
    162         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    163         h=this->GetSize();
    164 
    165         /* Start  looping on the number of gaussian points: */
    166         GaussSeg *gauss=new GaussSeg(2);
    167         for(int ig=gauss->begin();ig<gauss->end();ig++){
    168 
    169                 gauss->GaussPoint(ig);
    170 
    171                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    172                 GetNodalFunctions(basis,gauss);
    173 
    174                 vx_input->GetInputValue(&vx,gauss);
    175 
    176                 D_scalar=gauss->weight*Jdet;
    177 
    178                 TripleMultiply(basis,1,numnodes,1,
    179                                         &D_scalar,1,1,0,
    180                                         basis,1,numnodes,0,
    181                                         &Ke->values[0],1);
    182 
    183                 GetNodalFunctions(B,gauss);
    184                 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
    185 
    186                 D_scalar=dt*gauss->weight*Jdet*vx;
    187                 TripleMultiply(B,1,numnodes,1,
    188                                         &D_scalar,1,1,0,
    189                                         Bprime,1,numnodes,0,
    190                                         &Ke->values[0],1);
    191 
    192                 if(stabilization==2){
    193                         /*Streamline upwinding*/
    194                         vel=fabs(vx)+1.e-8;
    195                         D_scalar=dt*gauss->weight*Jdet*h/(2.*vel)*vx;
    196                 }
    197                 else if(stabilization==1){
    198                         /*SSA*/
    199                         vx_input->GetInputAverage(&vx);
    200                         D_scalar=dt*gauss->weight*Jdet*h/2.*fabs(vx);
    201                 }
    202                 if(stabilization==1 || stabilization==2){
    203                         TripleMultiply(Bprime,1,numnodes,1,
    204                                                 &D_scalar,1,1,0,
    205                                                 Bprime,1,numnodes,0,
    206                                                 &Ke->values[0],1);
    207                 }
    208         }
    209 
    210         /*Clean up and return*/
    211         xDelete<IssmDouble>(basis);
    212         xDelete<IssmDouble>(B);
    213         xDelete<IssmDouble>(Bprime);
    214         delete gauss;
    215         return Ke;
    216 }
    217 /*}}}*/
    218 /*FUNCTION Seg::CreateKMatrixFreeSurfaceBase {{{*/
    219 ElementMatrix* Seg::CreateKMatrixFreeSurfaceBase(void){
    220 
    221         /*Intermediaries */
    222         int        stabilization;
    223         IssmDouble Jdet,D_scalar,dt,h;
    224         IssmDouble vx,vel;
    225         IssmDouble xyz_list[NUMVERTICES][3];
    226 
    227         /*Fetch number of nodes for this finite element*/
    228         int numnodes = this->NumberofNodes();
    229 
    230         /*Initialize Element matrix and vectors*/
    231         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    232         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    233         IssmDouble*    B      = xNew<IssmDouble>(1*numnodes);
    234         IssmDouble*    Bprime = xNew<IssmDouble>(1*numnodes);
    235 
    236         /*Retrieve all inputs and parameters*/
    237         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    238         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    239         this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
    240         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    241         h=this->GetSize();
    242 
    243         /* Start  looping on the number of gaussian points: */
    244         GaussSeg *gauss=new GaussSeg(2);
    245         for(int ig=gauss->begin();ig<gauss->end();ig++){
    246 
    247                 gauss->GaussPoint(ig);
    248 
    249                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    250                 GetNodalFunctions(basis,gauss);
    251 
    252                 vx_input->GetInputValue(&vx,gauss);
    253 
    254                 D_scalar=gauss->weight*Jdet;
    255 
    256                 TripleMultiply(basis,1,numnodes,1,
    257                                         &D_scalar,1,1,0,
    258                                         basis,1,numnodes,0,
    259                                         &Ke->values[0],1);
    260 
    261                 GetNodalFunctions(B,gauss);
    262                 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
    263 
    264                 D_scalar=dt*gauss->weight*Jdet*vx;
    265                 TripleMultiply(B,1,numnodes,1,
    266                                         &D_scalar,1,1,0,
    267                                         Bprime,1,numnodes,0,
    268                                         &Ke->values[0],1);
    269 
    270                 if(stabilization==2){
    271                         /*Streamline upwinding*/
    272                         vel=fabs(vx)+1.e-8;
    273                         D_scalar=dt*gauss->weight*Jdet*h/(2.*vel)*vx;
    274                 }
    275                 else if(stabilization==1){
    276                         /*SSA*/
    277                         vx_input->GetInputAverage(&vx);
    278                         D_scalar=dt*gauss->weight*Jdet*h/2.*fabs(vx);
    279                 }
    280                 if(stabilization==1 || stabilization==2){
    281                         TripleMultiply(Bprime,1,numnodes,1,
    282                                                 &D_scalar,1,1,0,
    283                                                 Bprime,1,numnodes,0,
    284                                                 &Ke->values[0],1);
    285                 }
    286         }
    287 
    288         /*Clean up and return*/
    289         xDelete<IssmDouble>(basis);
    290         xDelete<IssmDouble>(B);
    291         xDelete<IssmDouble>(Bprime);
    292         delete gauss;
    293         return Ke;
    294 }
    295 /*}}}*/
    296 /*FUNCTION Seg::CreateMassMatrix {{{*/
    297 ElementMatrix* Seg::CreateMassMatrix(void){
    298 
    299         /* Intermediaries */
    300         IssmDouble  D,Jdet;
    301         IssmDouble  xyz_list[NUMVERTICES][3];
    302 
    303         /*Fetch number of nodes and dof for this finite element*/
    304         int numnodes = this->NumberofNodes();
    305 
    306         /*Initialize Element matrix and vectors*/
    307         ElementMatrix* Ke    = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    308         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    309 
    310         /*Retrieve all inputs and parameters*/
    311         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    312 
    313         /* Start looping on the number of gaussian points: */
    314         GaussSeg* gauss=new GaussSeg(2);
    315         for(int ig=gauss->begin();ig<gauss->end();ig++){
    316 
    317                 gauss->GaussPoint(ig);
    318 
    319                 GetNodalFunctions(basis,gauss);
    320                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    321                 D=gauss->weight*Jdet;
    322 
    323                 TripleMultiply(basis,1,numnodes,1,
    324                                         &D,1,1,0,
    325                                         basis,1,numnodes,0,
    326                                         &Ke->values[0],1);
    327         }
    328 
    329         /*Clean up and return*/
    330         delete gauss;
    331         xDelete<IssmDouble>(basis);
    332         return Ke;
    333 }
    334 /*}}}*/
    335140/*FUNCTION Seg::GetNumberOfNodes;{{{*/
    336141int Seg::GetNumberOfNodes(void){

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

    r16982 r16993  
    6969                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
    7070                void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");};
    71                 ElementMatrix* CreateKMatrix(void){_error_("not implemented yet");};
    7271                void        CreateDVector(Vector<IssmDouble>* df){_error_("not implemented yet");};
    73                 void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("not implemented yet");};
    7472                void        Delta18oParameterization(void){_error_("not implemented yet");};
    7573                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
    … …  
    257255                /*}}}*/
    258256                /*Seg specific routines:*/
    259                 ElementMatrix* CreateMassMatrix(void);
    260                 ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
    261                 ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
    262257                IssmDouble     GetSize(void);
    263258};

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

    r16982 r16993  
    191191}
    192192/*}}}*/
    193 /*FUNCTION Tria::CreateKMatrix(void) {{{*/
    194 ElementMatrix* Tria::CreateKMatrix(void){
    195 
    196         /*retreive parameters: */
    197         int analysis_type;
    198         int meshtype;
    199         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    200 
    201         /*Checks in debugging mode{{{*/
    202         _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs);
    203         /*}}}*/
    204 
    205         /*Skip if water element*/
    206         if(NoIceInElement()) return NULL;
    207 
    208         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    209         switch(analysis_type){
    210                 #ifdef _HAVE_STRESSBALANCE_
    211                 case StressbalanceAnalysisEnum:
    212                         int approximation;
    213                         inputs->GetInputValue(&approximation,ApproximationEnum);
    214                         switch(approximation){
    215                                 case SSAApproximationEnum:
    216                                         return CreateKMatrixStressbalanceSSA();
    217                                 case FSApproximationEnum:
    218                                         return CreateKMatrixStressbalanceFS();
    219                                 case SIAApproximationEnum:
    220                                         return NULL;
    221                                 case NoneApproximationEnum:
    222                                         return NULL;
    223                                 default:
    224                                         _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
    225                         }
    226                         break;
    227                 case StressbalanceSIAAnalysisEnum:
    228                         return CreateKMatrixStressbalanceSIA();
    229                         break;
    230                 #endif
    231                 #ifdef _HAVE_DAMAGE_
    232                 case DamageEvolutionAnalysisEnum:
    233                         return CreateKMatrixDamageEvolution();
    234                         break;
    235                 #endif
    236                 case L2ProjectionBaseAnalysisEnum:
    237                         parameters->FindParam(&meshtype,MeshTypeEnum);
    238                         if(meshtype==Mesh2DverticalEnum){
    239                                 return CreateBasalMassMatrix();
    240                         }
    241                         else{
    242                                 return CreateMassMatrix();
    243                         }
    244                         break;
    245                 #ifdef _HAVE_MASSTRANSPORT_
    246                 case MasstransportAnalysisEnum:
    247                         return CreateKMatrixMasstransport();
    248                         break;
    249                 case FreeSurfaceTopAnalysisEnum:
    250                         parameters->FindParam(&meshtype,MeshTypeEnum);
    251                         switch(meshtype){
    252                                 case Mesh2DverticalEnum:
    253                                         return CreateKMatrixFreeSurfaceTop1D();
    254                                 case Mesh3DEnum:
    255                                         return CreateKMatrixFreeSurfaceTop();
    256                                 default:
    257                                         _error_("Mesh not supported yet");
    258                         }
    259                         break;
    260                 case FreeSurfaceBaseAnalysisEnum:
    261                         parameters->FindParam(&meshtype,MeshTypeEnum);
    262                         switch(meshtype){
    263                                 case Mesh2DverticalEnum:
    264                                         return CreateKMatrixFreeSurfaceBase1D();
    265                                 case Mesh3DEnum:
    266                                         return CreateKMatrixFreeSurfaceBase();
    267                                 default:
    268                                         _error_("Mesh not supported yet");
    269                         }
    270                         break;
    271                 case ExtrudeFromBaseAnalysisEnum: case ExtrudeFromTopAnalysisEnum:
    272                         return CreateKMatrixExtrusion();
    273                 #endif
    274                 #ifdef _HAVE_HYDROLOGY_
    275                 case HydrologyShreveAnalysisEnum:
    276                         return CreateKMatrixHydrologyShreve();
    277                         break;
    278                 case HydrologyDCInefficientAnalysisEnum:
    279                         return CreateKMatrixHydrologyDCInefficient();
    280                         break;
    281                 case HydrologyDCEfficientAnalysisEnum:
    282                         return CreateKMatrixHydrologyDCEfficient();
    283                         break;
    284           case L2ProjectionEPLAnalysisEnum:
    285                         return CreateEPLDomainMassMatrix();
    286                         break;
    287                 #endif
    288                 #ifdef _HAVE_BALANCED_
    289                 case BalancethicknessAnalysisEnum:
    290                         return CreateKMatrixBalancethickness();
    291                         break;
    292                 case BalancevelocityAnalysisEnum:
    293                         return CreateKMatrixBalancevelocity();
    294                         break;
    295                 case SmoothedSurfaceSlopeXAnalysisEnum: case SmoothedSurfaceSlopeYAnalysisEnum:
    296                         return CreateKMatrixSmoothedSlope();
    297                         break;
    298                 #endif
    299                 #ifdef _HAVE_CONTROL_
    300                 case AdjointBalancethicknessAnalysisEnum:
    301                         return CreateKMatrixAdjointBalancethickness();
    302                         break;
    303                 case AdjointHorizAnalysisEnum:
    304                         return CreateKMatrixAdjointSSA();
    305                         break;
    306                 #endif
    307                 default:
    308                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    309         }
    310 
    311         /*Make compiler happy*/
    312         return NULL;
    313 }
    314 /*}}}*/
    315 /*FUNCTION Tria::CreateMassMatrix {{{*/
    316 ElementMatrix* Tria::CreateMassMatrix(void){
    317 
    318         /* Intermediaries */
    319         IssmDouble  D,Jdet;
    320         IssmDouble  xyz_list[NUMVERTICES][3];
    321 
    322         /*Fetch number of nodes and dof for this finite element*/
    323         int numnodes = this->NumberofNodes();
    324 
    325         /*Initialize Element matrix and vectors*/
    326         ElementMatrix* Ke    = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    327         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    328 
    329         /*Retrieve all inputs and parameters*/
    330         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    331 
    332         /* Start looping on the number of gaussian points: */
    333         GaussTria* gauss=new GaussTria(2);
    334         for(int ig=gauss->begin();ig<gauss->end();ig++){
    335 
    336                 gauss->GaussPoint(ig);
    337 
    338                 GetNodalFunctions(basis,gauss);
    339                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    340                 D=gauss->weight*Jdet;
    341 
    342                 TripleMultiply(basis,1,numnodes,1,
    343                                         &D,1,1,0,
    344                                         basis,1,numnodes,0,
    345                                         &Ke->values[0],1);
    346         }
    347 
    348         /*Clean up and return*/
    349         delete gauss;
    350         xDelete<IssmDouble>(basis);
    351         return Ke;
    352 }
    353 /*}}}*/
    354 /*FUNCTION Tria::CreateBasalMassMatrix{{{*/
    355 ElementMatrix* Tria::CreateBasalMassMatrix(void){
    356 
    357         if(!HasEdgeOnBed()) return NULL;
    358 
    359         int index1,index2;
    360         this->EdgeOnBedIndices(&index1,&index2);
    361 
    362         Seg* seg=(Seg*)SpawnSeg(index1,index2);
    363         ElementMatrix* Ke=seg->CreateMassMatrix();
    364         delete seg->material; delete seg;
    365 
    366         /*clean up and return*/
    367         return Ke;
    368 }
    369 /*}}}*/
    370193/*FUNCTION Tria::CreateDVector {{{*/
    371194void  Tria::CreateDVector(Vector<IssmDouble>* df){
    372195
    373196        /*Nothing done yet*/
    374 }
    375 /*}}}*/
    376 /*FUNCTION Tria::CreateJacobianMatrix{{{*/
    377 void  Tria::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){
    378 
    379         /*retrieve parameters: */
    380         ElementMatrix* Ke=NULL;
    381         int analysis_type;
    382         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    383 
    384         /*Checks in debugging {{{*/
    385         _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs);
    386         /*}}}*/
    387 
    388         /*Skip if water element*/
    389         if(NoIceInElement()) return;
    390 
    391         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    392         switch(analysis_type){
    393 #ifdef _HAVE_STRESSBALANCE_
    394                 case StressbalanceAnalysisEnum:
    395                         Ke=CreateJacobianStressbalanceSSA();
    396                         break;
    397 #endif
    398                 default:
    399                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    400         }
    401 
    402         /*Add to global matrix*/
    403         if(Ke){
    404                 Ke->AddToGlobal(Jff);
    405                 delete Ke;
    406         }
    407197}
    408198/*}}}*/
    … …  
    32733063
    32743064#ifdef _HAVE_STRESSBALANCE_
    3275 /*FUNCTION Tria::CreateKMatrixStressbalanceFS{{{*/
    3276 ElementMatrix* Tria::CreateKMatrixStressbalanceFS(void){
    3277 
    3278         ElementMatrix* Ke1 = NULL;
    3279         ElementMatrix* Ke2 = NULL;
    3280         ElementMatrix* Ke  = NULL;
    3281 
    3282         /*compute all stiffness matrices for this element*/
    3283         Ke1=CreateKMatrixStressbalanceFSViscous();
    3284         Ke2=CreateKMatrixStressbalanceFSFriction();
    3285         Ke =new ElementMatrix(Ke1,Ke2);
    3286 
    3287         /*clean-up and return*/
    3288         delete Ke1;
    3289         delete Ke2;
    3290         return Ke;
    3291 
    3292 }
    3293 /*}}}*/
    3294 /*FUNCTION Tria::CreateKMatrixStressbalanceFSViscous {{{*/
    3295 ElementMatrix* Tria::CreateKMatrixStressbalanceFSViscous(void){
    3296 
    3297         /*Intermediaries */
    3298         int        i,approximation;
    3299         IssmDouble Jdet,viscosity,FSreconditioning,D_scalar;
    3300         IssmDouble xyz_list[NUMVERTICES][3];
    3301         IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    3302         GaussTria *gauss=NULL;
    3303 
    3304         /*Fetch number of nodes and dof for this finite element*/
    3305         int vnumnodes = this->NumberofNodesVelocity();
    3306         int pnumnodes = this->NumberofNodesPressure();
    3307         int numdof    = vnumnodes*NDOF2 + pnumnodes*NDOF1;
    3308 
    3309         /*Prepare coordinate system list*/
    3310         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    3311         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
    3312         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    3313 
    3314         /*Initialize Element matrix and vectors*/
    3315         ElementMatrix* Ke     = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    3316         IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
    3317         IssmDouble*    Bprime = xNew<IssmDouble>(5*numdof);
    3318         IssmDouble*    D      = xNewZeroInit<IssmDouble>(5*5);
    3319 
    3320         /*Retrieve all inputs and parameters*/
    3321         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3322         parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    3323         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    3324         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    3325 
    3326         /* Start  looping on the number of gaussian points: */
    3327         gauss=new GaussTria(5);
    3328         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3329 
    3330                 gauss->GaussPoint(ig);
    3331 
    3332                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3333                 GetBFS(B,&xyz_list[0][0],gauss);
    3334                 GetBprimeFS(Bprime,&xyz_list[0][0],gauss);
    3335 
    3336                 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    3337                 material->GetViscosity2dvertical(&viscosity,&epsilon[0]);
    3338 
    3339                 D_scalar=gauss->weight*Jdet;
    3340                 for(i=0;i<3;i++) D[i*5+i] = +D_scalar*2.*viscosity;
    3341                 for(i=3;i<5;i++) D[i*5+i] = -D_scalar*FSreconditioning;
    3342 
    3343                 TripleMultiply(B,5,numdof,1,
    3344                                         D,5,5,0,
    3345                                         Bprime,5,numdof,0,
    3346                                         Ke->values,1);
    3347         }
    3348 
    3349         /*Transform Coordinate System*/
    3350         ::TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);
    3351 
    3352         /*Clean up and return*/
    3353         delete gauss;
    3354         xDelete<int>(cs_list);
    3355         xDelete<IssmDouble>(B);
    3356         xDelete<IssmDouble>(Bprime);
    3357         xDelete<IssmDouble>(D);
    3358         return Ke;
    3359 }
    3360 /*}}}*/
    3361 /*FUNCTION Tria::CreateKMatrixStressbalanceFSFriction{{{*/
    3362 ElementMatrix* Tria::CreateKMatrixStressbalanceFSFriction(void){
    3363 
    3364         /*Intermediaries */
    3365         int         i,j;
    3366         int         analysis_type;
    3367         int         indices[2];
    3368         IssmDouble  alpha2,Jdet;
    3369         IssmDouble  xyz_list[NUMVERTICES][3];
    3370         IssmDouble  xyz_list_seg[NUMVERTICES1D][3];
    3371         Friction   *friction = NULL;
    3372 
    3373         /*Initialize Element matrix and return if necessary*/
    3374         if(IsFloating() || !HasEdgeOnBed()) return NULL;
    3375 
    3376         /*Fetch number of nodes and dof for this finite element*/
    3377         int vnumnodes = this->NumberofNodesVelocity();
    3378         int pnumnodes = this->NumberofNodesPressure();
    3379         int  numdof   = vnumnodes*NDOF2 + pnumnodes*NDOF1;
    3380 
    3381         /*Prepare coordinate system list*/
    3382         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    3383         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
    3384         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    3385 
    3386         /*Initialize Element matrix and vectors*/
    3387         ElementMatrix* Ke        = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    3388         IssmDouble*    BFriction = xNew<IssmDouble>(1*numdof);
    3389         IssmDouble     D;
    3390 
    3391         /*Retrieve all inputs and parameters*/
    3392         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3393         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3394         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    3395         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    3396 
    3397         /*Get vertex indices that lie on bed*/
    3398         this->EdgeOnBedIndices(&indices[0],&indices[1]);
    3399         for(i=0;i<NUMVERTICES1D;i++) for(j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
    3400 
    3401         /*build friction object, used later on: */
    3402         friction=new Friction(this,1);
    3403 
    3404         /* Start looping on the number of gauss 1d (nodes on the bedrock) */
    3405         GaussTria* gauss=new GaussTria(indices[0],indices[1],3);
    3406         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3407 
    3408                 gauss->GaussPoint(ig);
    3409 
    3410                 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);
    3411                 GetLFS(BFriction,gauss);
    3412 
    3413                 friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,NULL);
    3414                 D = +alpha2*gauss->weight*Jdet; //taub_x = -alpha2 vx
    3415 
    3416                 TripleMultiply(BFriction,1,numdof,1,
    3417                                         &D,1,1,0,
    3418                                         BFriction,1,numdof,0,
    3419                                         Ke->values,1);
    3420         }
    3421 
    3422         /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
    3423         //::TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);
    3424 
    3425         /*Clean up and return*/
    3426         delete gauss;
    3427         delete friction;
    3428         xDelete<IssmDouble>(BFriction);
    3429         xDelete<int>(cs_list);
    3430         return Ke;
    3431 }
    3432 /*}}}*/
    3433 /*FUNCTION Tria::CreateKMatrixStressbalanceSSA {{{*/
    3434 ElementMatrix* Tria::CreateKMatrixStressbalanceSSA(void){
    3435 
    3436         /*compute all stiffness matrices for this element*/
    3437         ElementMatrix* Ke1=CreateKMatrixStressbalanceSSAViscous();
    3438         ElementMatrix* Ke2=CreateKMatrixStressbalanceSSAFriction();
    3439         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    3440 
    3441         /*clean-up and return*/
    3442         delete Ke1;
    3443         delete Ke2;
    3444         return Ke;
    3445 }
    3446 /*}}}*/
    3447 /*FUNCTION Tria::CreateKMatrixStressbalanceSSAViscous{{{*/
    3448 ElementMatrix* Tria::CreateKMatrixStressbalanceSSAViscous(void){
    3449 
    3450         /*Intermediaries*/
    3451         IssmDouble xyz_list[NUMVERTICES][3];
    3452         IssmDouble viscosity,newviscosity,oldviscosity;
    3453         IssmDouble viscosity_overshoot,thickness,Jdet;
    3454         IssmDouble epsilon[3],oldepsilon[3];    /* epsilon=[exx,eyy,exy];    */
    3455         IssmDouble D_scalar;
    3456 
    3457         /*Fetch number of nodes and dof for this finite element*/
    3458         int numnodes = this->NumberofNodes();
    3459         int numdof   = numnodes*NDOF2;
    3460 
    3461         /*Initialize Element matrix and vectors*/
    3462         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,SSAApproximationEnum);
    3463         IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
    3464         IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
    3465         IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
    3466 
    3467         /*Retrieve all inputs and parameters*/
    3468         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3469         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    3470         Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    3471         Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    3472         Input* vxold_input=inputs->GetInput(VxPicardEnum);      _assert_(vxold_input);
    3473         Input* vyold_input=inputs->GetInput(VyPicardEnum);      _assert_(vyold_input);
    3474         this->parameters->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
    3475 
    3476         /* Start  looping on the number of gaussian points: */
    3477         GaussTria* gauss  = new GaussTria(2);
    3478         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3479 
    3480                 gauss->GaussPoint(ig);
    3481 
    3482                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3483                 GetBSSA(&B[0], &xyz_list[0][0], gauss);
    3484                 GetBprimeSSA(&Bprime[0], &xyz_list[0][0], gauss);
    3485 
    3486                 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    3487                 this->StrainRateSSA(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    3488                 material->GetViscosity2d(&viscosity, &epsilon[0]);
    3489                 material->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
    3490                 thickness_input->GetInputValue(&thickness, gauss);
    3491 
    3492                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    3493                 D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet;
    3494                 for(int i=0;i<3;i++) D[i*3+i]=D_scalar;
    3495 
    3496                 TripleMultiply(B,3,numdof,1,
    3497                                         D,3,3,0,
    3498                                         Bprime,3,numdof,0,
    3499                                         &Ke->values[0],1);
    3500         }
    3501 
    3502         /*Transform Coordinate System*/
    3503         ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);
    3504 
    3505         /*Clean up and return*/
    3506         delete gauss;
    3507         xDelete<IssmDouble>(D);
    3508         xDelete<IssmDouble>(Bprime);
    3509         xDelete<IssmDouble>(B);
    3510         return Ke;
    3511 }
    3512 /*}}}*/
    3513 /*FUNCTION Tria::CreateKMatrixStressbalanceSSAFriction {{{*/
    3514 ElementMatrix* Tria::CreateKMatrixStressbalanceSSAFriction(void){
    3515 
    3516         /*Intermediaries*/
    3517         bool       mainlyfloating;
    3518         int        analysis_type,migration_style;
    3519         int        point1;
    3520         IssmDouble alpha2,Jdet;
    3521         IssmDouble fraction1,fraction2;
    3522         IssmDouble phi=1.0;
    3523         IssmDouble D_scalar;
    3524         IssmDouble gllevelset;
    3525         IssmDouble xyz_list[NUMVERTICES][3];
    3526         Friction  *friction = NULL;
    3527         GaussTria *gauss    = NULL;
    3528 
    3529         /*Return if element is inactive*/
    3530         if(IsFloating()) return NULL;
    3531 
    3532         /*Fetch number of nodes and dof for this finite element*/
    3533         int numnodes = this->NumberofNodes();
    3534         int numdof   = numnodes*NDOF2;
    3535 
    3536         /*Initialize Element matrix and vectors*/
    3537         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,SSAApproximationEnum);
    3538         IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
    3539         IssmDouble*    D      = xNewZeroInit<IssmDouble>(2*2);
    3540 
    3541         /*Retrieve all inputs and parameters*/
    3542         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3543         parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
    3544         Input* surface_input=inputs->GetInput(SurfaceEnum);       _assert_(surface_input);
    3545         Input* vx_input=inputs->GetInput(VxEnum);                 _assert_(vx_input);
    3546         Input* vy_input=inputs->GetInput(VyEnum);                 _assert_(vy_input);
    3547         Input* gllevelset_input=NULL;
    3548         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3549 
    3550         /*build friction object, used later on: */
    3551         friction=new Friction(this,2);
    3552 
    3553         /*Recover portion of element that is grounded*/
    3554         if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list[0][0]);
    3555         if(migration_style==SubelementMigration2Enum){
    3556                 gllevelset_input=inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
    3557                 this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    3558                 gauss=new GaussTria(point1,fraction1,fraction2,mainlyfloating,2);
    3559         }
    3560         else{
    3561                 gauss=new GaussTria(2);
    3562         }
    3563 
    3564         /* Start  looping on the number of gaussian points: */
    3565         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3566 
    3567                 gauss->GaussPoint(ig);
    3568 
    3569                 friction->GetAlpha2(&alpha2, gauss,vx_input,vy_input,NULL);
    3570                 if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
    3571                 if(migration_style==SubelementMigration2Enum){
    3572                         gllevelset_input->GetInputValue(&gllevelset, gauss);
    3573                         if(gllevelset<0) alpha2=0;
    3574                 }
    3575 
    3576                 GetBSSAFriction(&B[0], &xyz_list[0][0], gauss);
    3577                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3578                 D_scalar=alpha2*gauss->weight*Jdet;
    3579                 for(int i=0;i<2;i++) D[i*2+i]=D_scalar;
    3580 
    3581                 TripleMultiply(B,2,numdof,1,
    3582                                         D,2,2,0,
    3583                                         B,2,numdof,0,
    3584                                         &Ke->values[0],1);
    3585         }
    3586 
    3587         /*Transform Coordinate System*/
    3588         ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);
    3589 
    3590         /*Clean up and return*/
    3591         delete gauss;
    3592         delete friction;
    3593         xDelete<IssmDouble>(D);
    3594         xDelete<IssmDouble>(B);
    3595         return Ke;
    3596 }
    3597 /*}}}*/
    3598 /*FUNCTION Tria::CreateKMatrixStressbalanceSIA{{{*/
    3599 ElementMatrix* Tria::CreateKMatrixStressbalanceSIA(void){
    3600 
    3601         /*Intermediaries*/
    3602         IssmDouble connectivity;
    3603 
    3604         /*Fetch number of nodes and dof for this finite element*/
    3605         int numnodes = this->NumberofNodes(); _assert_(numnodes==3);
    3606         int numdof   = numnodes*NDOF2;
    3607 
    3608         /*Initialize Element matrix*/
    3609         ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    3610 
    3611         /*Create Element matrix*/
    3612         for(int i=0;i<numnodes;i++){
    3613                 connectivity=(IssmDouble)vertices[i]->Connectivity();
    3614                 Ke->values[(2*i+0)*numdof+(2*i+0)]=1./connectivity;
    3615                 Ke->values[(2*i+1)*numdof+(2*i+1)]=1./connectivity;
    3616         }
    3617 
    3618         /*Clean up and return*/
    3619         return Ke;
    3620 }
    3621 /*}}}*/
    3622 /*FUNCTION Tria::CreateJacobianStressbalanceSSA{{{*/
    3623 ElementMatrix* Tria::CreateJacobianStressbalanceSSA(void){
    3624 
    3625         /*Intermediaries */
    3626         int        i,j;
    3627         IssmDouble xyz_list[NUMVERTICES][3];
    3628         IssmDouble Jdet,thickness;
    3629         IssmDouble eps1dotdphii,eps1dotdphij;
    3630         IssmDouble eps2dotdphii,eps2dotdphij;
    3631         IssmDouble mu_prime;
    3632         IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/
    3633         IssmDouble eps1[2],eps2[2];
    3634         GaussTria* gauss = NULL;
    3635 
    3636         /*Fetch number of nodes and dof for this finite element*/
    3637         int numnodes = this->NumberofNodes();
    3638 
    3639         /*Initialize Element matrix, vectors and Gaussian points*/
    3640         ElementMatrix* Ke=CreateKMatrixStressbalanceSSA(); //Initialize Jacobian with regular SSA (first part of the Gateau derivative)
    3641         IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    3642 
    3643         /*Retrieve all inputs and parameters*/
    3644         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3645         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    3646         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    3647         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    3648 
    3649         /* Start  looping on the number of gaussian points: */
    3650         gauss=new GaussTria(2);
    3651         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3652 
    3653                 gauss->GaussPoint(ig);
    3654 
    3655                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3656                 GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
    3657 
    3658                 thickness_input->GetInputValue(&thickness, gauss);
    3659                 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    3660                 material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    3661                 eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    3662                 eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
    3663 
    3664                 for(i=0;i<numnodes;i++){
    3665                         for(j=0;j<numnodes;j++){
    3666                                 eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
    3667                                 eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
    3668                                 eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
    3669                                 eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];
    3670 
    3671                                 Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
    3672                                 Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
    3673                                 Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
    3674                                 Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
    3675                         }
    3676                 }
    3677         }
    3678 
    3679         /*Transform Coordinate System*/
    3680         ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);
    3681 
    3682         /*Clean up and return*/
    3683         xDelete<IssmDouble>(dbasis);
    3684         delete gauss;
    3685         return Ke;
    3686 }
    3687 /*}}}*/
    36883065/*FUNCTION Tria::GetYcoord {{{*/
    36893066IssmDouble Tria::GetYcoord(GaussTria* gauss){
    … …  
    50374414}
    50384415/*}}}*/
    5039 /*FUNCTION Tria::CreateKMatrixAdjointBalancethickness {{{*/
    5040 ElementMatrix* Tria::CreateKMatrixAdjointBalancethickness(void){
    5041 
    5042         ElementMatrix* Ke=NULL;
    5043 
    5044         /*Get Element Matrix of the forward model*/
    5045         switch(GetElementType()){
    5046                 case P1Enum:
    5047                         Ke=CreateKMatrixBalancethickness_CG();
    5048                         break;
    5049                 case P1DGEnum:
    5050                         Ke=CreateKMatrixBalancethickness_DG();
    5051                         break;
    5052                 default:
    5053                         _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
    5054         }
    5055 
    5056         /*Transpose and return Ke*/
    5057         Ke->Transpose();
    5058         return Ke;
    5059 }
    5060 /*}}}*/
    5061 /*FUNCTION Tria::CreateKMatrixAdjointSSA{{{*/
    5062 ElementMatrix* Tria::CreateKMatrixAdjointSSA(void){
    5063 
    5064         /*Intermediaries */
    5065         int         i,j;
    5066         bool        incomplete_adjoint;
    5067         IssmDouble  xyz_list[NUMVERTICES][3];
    5068         IssmDouble  Jdet,thickness;
    5069         IssmDouble  eps1dotdphii,eps1dotdphij;
    5070         IssmDouble  eps2dotdphii,eps2dotdphij;
    5071         IssmDouble  mu_prime;
    5072         IssmDouble  epsilon[3];/* epsilon=[exx,eyy,exy];*/
    5073         IssmDouble  eps1[2],eps2[2];
    5074         GaussTria  *gauss=NULL;
    5075 
    5076         /*Fetch number of nodes and dof for this finite element*/
    5077         int numnodes = this->NumberofNodes();
    5078 
    5079         /*Initialize Jacobian with regular SSA (first part of the Gateau derivative)*/
    5080         parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
    5081         ElementMatrix* Ke=CreateKMatrixStressbalanceSSA();
    5082         if(incomplete_adjoint) return Ke;
    5083 
    5084         /*Retrieve all inputs and parameters*/
    5085         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5086         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    5087         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    5088         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    5089 
    5090         /*Allocate dbasis*/
    5091         IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);
    5092 
    5093         /* Start  looping on the number of gaussian points: */
    5094         gauss=new GaussTria(2);
    5095         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5096 
    5097                 gauss->GaussPoint(ig);
    5098 
    5099                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    5100                 GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
    5101 
    5102                 thickness_input->GetInputValue(&thickness, gauss);
    5103                 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    5104                 material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    5105                 eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    5106                 eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
    5107 
    5108                 for(i=0;i<numnodes;i++){
    5109                         for(j=0;j<numnodes;j++){
    5110                                 eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];
    5111                                 eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];
    5112                                 eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];
    5113                                 eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];
    5114 
    5115                                 Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
    5116                                 Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
    5117                                 Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
    5118                                 Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
    5119                         }
    5120                 }
    5121         }
    5122 
    5123         /*Transform Coordinate System*/
    5124         ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);
    5125 
    5126         /*Clean up and return*/
    5127         delete gauss;
    5128         xDelete<IssmDouble>(dbasis);
    5129         //Ke->Transpose();
    5130         return Ke;
    5131 }
    5132 /*}}}*/
    51334416/*FUNCTION Tria::GetVectorFromControlInputs{{{*/
    51344417void  Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){
    … …  
    51914474
    51924475        ((ControlInput*)input)->SetInput(new_input);
    5193 }
    5194 /*}}}*/
    5195 #endif
    5196 
    5197 #ifdef _HAVE_THERMAL_
    5198 /*FUNCTION Tria::CreateKMatrixMelting {{{*/
    5199 ElementMatrix* Tria::CreateKMatrixMelting(void){
    5200 
    5201         /*Intermediaries */
    5202         IssmDouble heatcapacity,latentheat;
    5203         IssmDouble Jdet,D_scalar;
    5204         IssmDouble xyz_list[NUMVERTICES][3];
    5205 
    5206         /*Fetch number of nodes and dof for this finite element*/
    5207         int numnodes = this->NumberofNodes();
    5208 
    5209         /*Initialize Element vector and vectors*/
    5210         ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    5211         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    5212 
    5213         /*Retrieve all inputs and parameters*/
    5214         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5215         latentheat   = matpar->GetLatentHeat();
    5216         heatcapacity = matpar->GetHeatCapacity();
    5217 
    5218         /* Start looping on the number of gauss  (nodes on the bedrock) */
    5219         GaussTria* gauss=new GaussTria(2);
    5220         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5221 
    5222                 gauss->GaussPoint(ig);
    5223 
    5224                 GetNodalFunctions(&basis[0],gauss);
    5225                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0], gauss);
    5226 
    5227                 D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
    5228 
    5229                 TripleMultiply(basis,numnodes,1,0,
    5230                                         &D_scalar,1,1,0,
    5231                                         basis,1,numnodes,0,
    5232                                         &Ke->values[0],1);
    5233         }
    5234 
    5235         /*Clean up and return*/
    5236         xDelete<IssmDouble>(basis);
    5237         delete gauss;
    5238         return Ke;
    52394476}
    52404477/*}}}*/
    … …  
    52944531        this->inputs->AddInput(new TriaInput(HydrologyWaterVxEnum,vx,P1Enum));
    52954532        this->inputs->AddInput(new TriaInput(HydrologyWaterVyEnum,vy,P1Enum));
    5296 }
    5297 /*}}}*/
    5298 /*FUNCTION Tria::CreateKMatrixHydrologyShreve{{{*/
    5299 ElementMatrix* Tria::CreateKMatrixHydrologyShreve(void){
    5300 
    5301         /*Intermediaries */
    5302         IssmDouble diffusivity;
    5303         IssmDouble Jdet,D_scalar,dt,h;
    5304         IssmDouble vx,vy,vel,dvxdx,dvydy;
    5305         IssmDouble dvx[2],dvy[2];
    5306         IssmDouble xyz_list[NUMVERTICES][3];
    5307 
    5308         /*Skip if water or ice shelf element*/
    5309         if(NoIceInElement() || IsFloating()) return NULL;
    5310 
    5311         /*Fetch number of nodes and dof for this finite element*/
    5312         int numnodes = this->NumberofNodes();
    5313 
    5314         /*Initialize Element matrix and vectors*/
    5315         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    5316         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    5317         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    5318         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    5319         IssmDouble     D[2][2];
    5320 
    5321         /*Create water velocity vx and vy from current inputs*/
    5322         CreateHydrologyWaterVelocityInput();
    5323 
    5324         /*Retrieve all inputs and parameters*/
    5325         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5326         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    5327         this->parameters->FindParam(&diffusivity,HydrologyshreveStabilizationEnum);
    5328         Input* vx_input=inputs->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
    5329         Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
    5330         h=sqrt(2*this->GetArea());
    5331 
    5332         /* Start  looping on the number of gaussian points: */
    5333         GaussTria* gauss=new GaussTria(2);
    5334         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5335 
    5336                 gauss->GaussPoint(ig);
    5337 
    5338                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    5339                 GetNodalFunctions(basis,gauss);
    5340 
    5341                 vx_input->GetInputValue(&vx,gauss);
    5342                 vy_input->GetInputValue(&vy,gauss);
    5343                 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    5344                 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    5345 
    5346                 D_scalar=gauss->weight*Jdet;
    5347 
    5348                 TripleMultiply(basis,1,numnodes,1,
    5349                                         &D_scalar,1,1,0,
    5350                                         basis,1,numnodes,0,
    5351                                         Ke->values,1);
    5352 
    5353                 GetBMasstransport(B,&xyz_list[0][0], gauss);
    5354                 GetBprimeMasstransport(Bprime,&xyz_list[0][0], gauss);
    5355 
    5356                 dvxdx=dvx[0];
    5357                 dvydy=dvy[1];
    5358                 D_scalar=dt*gauss->weight*Jdet;
    5359 
    5360                 D[0][0]=D_scalar*dvxdx;
    5361                 D[0][1]=0.;
    5362                 D[1][0]=0.;
    5363                 D[1][1]=D_scalar*dvydy;
    5364                 TripleMultiply(B,2,numnodes,1,
    5365                                         &D[0][0],2,2,0,
    5366                                         B,2,numnodes,0,
    5367                                         &Ke->values[0],1);
    5368 
    5369                 D[0][0]=D_scalar*vx;
    5370                 D[1][1]=D_scalar*vy;
    5371                 TripleMultiply(B,2,numnodes,1,
    5372                                         &D[0][0],2,2,0,
    5373                                         Bprime,2,numnodes,0,
    5374                                         &Ke->values[0],1);
    5375 
    5376                 /*Artificial diffusivity*/
    5377                 vel=sqrt(vx*vx+vy*vy);
    5378                 D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx;
    5379                 D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx;
    5380                 D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;
    5381                 D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;
    5382                 TripleMultiply(Bprime,2,numnodes,1,
    5383                                         &D[0][0],2,2,0,
    5384                                         Bprime,2,numnodes,0,
    5385                                         &Ke->values[0],1);
    5386         }
    5387 
    5388         /*Clean up and return*/
    5389         xDelete<IssmDouble>(basis);
    5390         xDelete<IssmDouble>(B);
    5391         xDelete<IssmDouble>(Bprime);
    5392         delete gauss;
    5393         return Ke;
    5394 }
    5395 /*}}}*/
    5396 /*FUNCTION Tria::CreateKMatrixHydrologyDCInefficient{{{*/
    5397 ElementMatrix* Tria::CreateKMatrixHydrologyDCInefficient(void){
    5398 
    5399         /* Intermediaries */
    5400         IssmDouble  D_scalar,Jdet;
    5401         IssmDouble      sediment_transmitivity,dt;
    5402         IssmDouble  sediment_storing;
    5403         IssmDouble  xyz_list[NUMVERTICES][3];
    5404 
    5405         /*Fetch number of nodes and dof for this finite element*/
    5406         int numnodes = this->NumberofNodes();
    5407 
    5408         /*Initialize Element matrix and vectors*/
    5409         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
    5410         IssmDouble*    B      = xNew<IssmDouble>(5*numnodes);
    5411         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    5412         IssmDouble     D[2][2];
    5413 
    5414         /*Retrieve all inputs and parameters*/
    5415         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5416         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    5417         sediment_storing       = matpar->GetSedimentStoring();
    5418         sediment_transmitivity = matpar->GetSedimentTransmitivity();
    5419 
    5420         /* Start looping on the number of gaussian points: */
    5421         GaussTria* gauss=new GaussTria(2);
    5422         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5423 
    5424                 gauss->GaussPoint(ig);
    5425 
    5426                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    5427 
    5428                 /*Diffusivity*/
    5429                 D_scalar=sediment_transmitivity*gauss->weight*Jdet;
    5430                 if(reCast<bool,IssmDouble>(dt)) D_scalar=D_scalar*dt;
    5431                 D[0][0]=D_scalar; D[0][1]=0.;
    5432                 D[1][0]=0.;       D[1][1]=D_scalar;
    5433                 GetBHydro(B,&xyz_list[0][0],gauss);
    5434                 TripleMultiply(B,2,numnodes,1,
    5435                                         &D[0][0],2,2,0,
    5436                                         B,2,numnodes,0,
    5437                                         &Ke->values[0],1);
    5438 
    5439                 /*Transient*/
    5440                 if(reCast<bool,IssmDouble>(dt)){
    5441                         GetNodalFunctions(&basis[0],gauss);
    5442                         D_scalar=sediment_storing*gauss->weight*Jdet;
    5443 
    5444                         TripleMultiply(basis,numnodes,1,0,
    5445                                                 &D_scalar,1,1,0,
    5446                                                 basis,1,numnodes,0,
    5447                                                 &Ke->values[0],1);
    5448                 }
    5449         }
    5450         /*Clean up and return*/
    5451         xDelete<IssmDouble>(basis);
    5452         xDelete<IssmDouble>(B);
    5453         delete gauss;
    5454         return Ke;
    5455 }
    5456 /*}}}*/
    5457 /*FUNCTION Tria::CreateKMatrixHydrologyDCEfficient{{{*/
    5458 ElementMatrix* Tria::CreateKMatrixHydrologyDCEfficient(void){
    5459 
    5460         /* Intermediaries */
    5461         IssmDouble  D_scalar,Jdet,dt;
    5462         IssmDouble  epl_thickness;
    5463         IssmDouble      epl_conductivity;
    5464         IssmDouble  epl_specificstoring;
    5465         IssmDouble  xyz_list[NUMVERTICES][3];
    5466 
    5467         /*Check that all nodes are active, else return empty matrix*/
    5468         if(!this->AllActive()){
    5469                 return NULL;
    5470         }
    5471 
    5472         /*Fetch number of nodes and dof for this finite element*/
    5473         int numnodes = this->NumberofNodes();
    5474 
    5475         /*Initialize Element matrix and vectors*/
    5476         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
    5477         IssmDouble*    B      = xNew<IssmDouble>(5*numnodes);
    5478         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    5479         IssmDouble     D[2][2];
    5480 
    5481         /*Retrieve all inputs and parameters*/
    5482         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5483         Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input);
    5484         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    5485         epl_specificstoring = matpar->GetEplSpecificStoring();
    5486         epl_conductivity    = matpar->GetEplConductivity();
    5487 
    5488 
    5489         /* Start looping on the number of gaussian points: */
    5490         GaussTria* gauss=new GaussTria(2);
    5491         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5492                
    5493                
    5494                 gauss->GaussPoint(ig);
    5495                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    5496                 thickness_input->GetInputValue(&epl_thickness,gauss);
    5497 
    5498                 /*Diffusivity*/
    5499                 D_scalar=epl_conductivity*epl_thickness*gauss->weight*Jdet;
    5500                 if(reCast<bool,IssmDouble>(dt)) D_scalar=D_scalar*dt;
    5501                 D[0][0]=D_scalar; D[0][1]=0.;
    5502                 D[1][0]=0.;       D[1][1]=D_scalar;
    5503                 GetBHydro(B,&xyz_list[0][0],gauss);
    5504                 TripleMultiply(B,2,numnodes,1,
    5505                                         &D[0][0],2,2,0,
    5506                                         B,2,numnodes,0,
    5507                                         &Ke->values[0],1);
    5508 
    5509                 /*Transient*/
    5510                 if(reCast<bool,IssmDouble>(dt)){
    5511                         GetNodalFunctions(basis,gauss);
    5512                         D_scalar=epl_specificstoring*epl_thickness*gauss->weight*Jdet;
    5513 
    5514                         TripleMultiply(basis,numnodes,1,0,
    5515                                                 &D_scalar,1,1,0,
    5516                                                 basis,1,numnodes,0,
    5517                                                 &Ke->values[0],1);
    5518                 }
    5519         }
    5520 
    5521         /*Clean up and return*/
    5522         xDelete<IssmDouble>(basis);
    5523         xDelete<IssmDouble>(B);
    5524         delete gauss;
    5525         return Ke;
    55264533}
    55274534/*}}}*/
    … …  
    58984905#endif
    58994906
    5900 #ifdef _HAVE_MASSTRANSPORT_
    5901 /*FUNCTION Tria::CreateKMatrixExtrusion {{{*/
    5902 ElementMatrix* Tria::CreateKMatrixExtrusion(void){
    5903 
    5904         /*compute all stiffness matrices for this element*/
    5905         ElementMatrix* Ke1=CreateKMatrixExtrusionVolume();
    5906         ElementMatrix* Ke2=CreateKMatrixExtrusionSurface();
    5907         ElementMatrix* Ke3=CreateKMatrixExtrusionBed();
    5908         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    5909 
    5910         /*clean-up and return*/
    5911         delete Ke1;
    5912         delete Ke2;
    5913         delete Ke3;
    5914         return Ke;
    5915 
    5916 }
    5917 /*}}}*/
    5918 /*FUNCTION Tria::CreateKMatrixExtrusionVolume {{{*/
    5919 ElementMatrix* Tria::CreateKMatrixExtrusionVolume(void){
    5920 
    5921         /*Intermediaries */
    5922         IssmDouble  Jdet;
    5923         IssmDouble  xyz_list[NUMVERTICES][3];
    5924         IssmDouble  B[NDOF1][NUMVERTICES];
    5925         IssmDouble  Bprime[NDOF1][NUMVERTICES];
    5926         IssmDouble  DL_scalar;
    5927         GaussTria  *gauss=NULL;
    5928 
    5929         /*Initialize Element matrix*/
    5930         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    5931 
    5932         /*Retrieve all inputs and parameters*/
    5933         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5934 
    5935         /* Start  looping on the number of gaussian points: */
    5936         gauss=new GaussTria(2);
    5937         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5938 
    5939                 gauss->GaussPoint(ig);
    5940 
    5941                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    5942                 GetBExtrusion(&B[0][0], &xyz_list[0][0], gauss);
    5943                 GetNodalFunctions(&Bprime[0][0],gauss);
    5944 
    5945                 DL_scalar=gauss->weight*Jdet;
    5946 
    5947                 TripleMultiply(&B[0][0],1,NUMVERTICES,1,
    5948                                         &DL_scalar,1,1,0,
    5949                                         &Bprime[0][0],1,NUMVERTICES,0,
    5950                                         &Ke->values[0],1);
    5951         }
    5952 
    5953         /*Clean up and return*/
    5954         delete gauss;
    5955         return Ke;
    5956 }
    5957 /*}}}*/
    5958 /*FUNCTION Tria::CreateKMatrixExtrusionSurface {{{*/
    5959 ElementMatrix* Tria::CreateKMatrixExtrusionSurface(void){
    5960 
    5961         if (!HasEdgeOnSurface()) return NULL;
    5962 
    5963         /*Constants*/
    5964         const int numdof=NDOF1*NUMVERTICES;
    5965 
    5966         /*Intermediaries */
    5967         int indices[2];
    5968         IssmDouble xyz_list[NUMVERTICES][3];
    5969         IssmDouble xyz_list_seg[NUMVERTICES1D][3];
    5970         IssmDouble normal[3];
    5971         IssmDouble Jdet,DL_scalar;
    5972         IssmDouble basis[NUMVERTICES];
    5973         GaussTria *gauss=NULL;
    5974 
    5975         /*Initialize Element matrix*/
    5976         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    5977 
    5978         /*Retrieve all inputs and parameters*/
    5979         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5980 
    5981         /*Get vertex indices that lie on bed*/
    5982         this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
    5983         for(int i=0;i<NUMVERTICES1D;i++) for(int j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
    5984         NormalSection(&normal[0],&xyz_list_seg[0][0]);
    5985 
    5986         /* Start  looping on the number of gaussian points: */
    5987         gauss=new GaussTria(indices[0],indices[1],2);
    5988         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5989 
    5990                 gauss->GaussPoint(ig);
    5991 
    5992                 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);
    5993                 GetNodalFunctions(&basis[0], gauss);
    5994 
    5995                 DL_scalar= - gauss->weight*Jdet*normal[1];
    5996 
    5997                 TripleMultiply( basis,1,numdof,1,
    5998                                         &DL_scalar,1,1,0,
    5999                                         basis,1,numdof,0,
    6000                                         &Ke->values[0],1);
    6001         }
    6002 
    6003         /*Clean up and return*/
    6004         delete gauss;
    6005         return Ke;
    6006 }
    6007 /*}}}*/
    6008 /*FUNCTION Tria::CreateKMatrixExtrusionBed {{{*/
    6009 ElementMatrix* Tria::CreateKMatrixExtrusionBed(void){
    6010 
    6011         if (!HasEdgeOnBed()) return NULL;
    6012 
    6013         /*Constants*/
    6014         const int numdof=NDOF1*NUMVERTICES;
    6015 
    6016         /*Intermediaries */
    6017         int indices[2];
    6018         IssmDouble xyz_list[NUMVERTICES][3];
    6019         IssmDouble xyz_list_seg[NUMVERTICES1D][3];
    6020         IssmDouble normal[3];
    6021         IssmDouble Jdet,DL_scalar;
    6022         IssmDouble basis[NUMVERTICES];
    6023         GaussTria *gauss=NULL;
    6024 
    6025         /*Initialize Element matrix*/
    6026         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    6027 
    6028         /*Retrieve all inputs and parameters*/
    6029         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6030 
    6031         /*Get vertex indices that lie on bed*/
    6032         this->EdgeOnBedIndices(&indices[0],&indices[1]);
    6033         for(int i=0;i<NUMVERTICES1D;i++) for(int j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
    6034         NormalSection(&normal[0],&xyz_list_seg[0][0]);
    6035 
    6036         /* Start  looping on the number of gaussian points: */
    6037         gauss=new GaussTria(indices[0],indices[1],2);
    6038         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6039 
    6040                 gauss->GaussPoint(ig);
    6041 
    6042                 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);
    6043                 GetNodalFunctions(&basis[0], gauss);
    6044 
    6045                 DL_scalar= - gauss->weight*Jdet*normal[1];
    6046 
    6047                 TripleMultiply( basis,1,numdof,1,
    6048                                         &DL_scalar,1,1,0,
    6049                                         basis,1,numdof,0,
    6050                                         &Ke->values[0],1);
    6051         }
    6052 
    6053         /*Clean up and return*/
    6054         delete gauss;
    6055         return Ke;
    6056 }
    6057 /*}}}*/
    6058 /*FUNCTION Tria::CreateKMatrixMasstransport {{{*/
    6059 ElementMatrix* Tria::CreateKMatrixMasstransport(void){
    6060 
    6061         switch(GetElementType()){
    6062                 case P1Enum: case P2Enum:
    6063                         return CreateKMatrixMasstransport_CG();
    6064                 case P1DGEnum:
    6065                         return CreateKMatrixMasstransport_DG();
    6066                 default:
    6067                         _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
    6068         }
    6069 }
    6070 /*}}}*/
    6071 /*FUNCTION Tria::CreateKMatrixMasstransport_CG {{{*/
    6072 ElementMatrix* Tria::CreateKMatrixMasstransport_CG(void){
    6073 
    6074         /*Intermediaries */
    6075         int        stabilization;
    6076         int        meshtype;
    6077         IssmDouble Jdet,D_scalar,dt,h;
    6078         IssmDouble vel,vx,vy,dvxdx,dvydy;
    6079         IssmDouble dvx[2],dvy[2];
    6080         IssmDouble xyz_list[NUMVERTICES][3];
    6081 
    6082         /*Fetch number of nodes for this finite element*/
    6083         int numnodes = this->NumberofNodes();
    6084 
    6085         /*Initialize Element matrix and vectors*/
    6086         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    6087         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    6088         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    6089         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    6090         IssmDouble     D[2][2];
    6091 
    6092         /*Retrieve all inputs and parameters*/
    6093         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6094         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6095         this->parameters->FindParam(&meshtype,MeshTypeEnum);
    6096         this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
    6097         Input* vxaverage_input=NULL;
    6098         Input* vyaverage_input=NULL;
    6099         if(meshtype==Mesh2DhorizontalEnum){
    6100                 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    6101                 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    6102         }
    6103         else{
    6104                 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    6105                 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
    6106         }
    6107         h=sqrt(2*this->GetArea());
    6108 
    6109         /* Start  looping on the number of gaussian points: */
    6110         GaussTria *gauss=new GaussTria(2);
    6111         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6112 
    6113                 gauss->GaussPoint(ig);
    6114 
    6115                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6116                 GetNodalFunctions(basis,gauss);
    6117 
    6118                 vxaverage_input->GetInputValue(&vx,gauss);
    6119                 vyaverage_input->GetInputValue(&vy,gauss);
    6120                 vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    6121                 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    6122 
    6123                 D_scalar=gauss->weight*Jdet;
    6124 
    6125                 TripleMultiply(basis,1,numnodes,1,
    6126                                         &D_scalar,1,1,0,
    6127                                         basis,1,numnodes,0,
    6128                                         &Ke->values[0],1);
    6129 
    6130                 GetBMasstransport(B,&xyz_list[0][0],gauss);
    6131                 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
    6132 
    6133                 dvxdx=dvx[0];
    6134                 dvydy=dvy[1];
    6135                 D_scalar=dt*gauss->weight*Jdet;
    6136 
    6137                 D[0][0]=D_scalar*dvxdx;
    6138                 D[0][1]=0.;
    6139                 D[1][0]=0.;
    6140                 D[1][1]=D_scalar*dvydy;
    6141                 TripleMultiply(B,2,numnodes,1,
    6142                                         &D[0][0],2,2,0,
    6143                                         B,2,numnodes,0,
    6144                                         &Ke->values[0],1);
    6145 
    6146                 D[0][0]=D_scalar*vx;
    6147                 D[1][1]=D_scalar*vy;
    6148                 TripleMultiply(B,2,numnodes,1,
    6149                                         &D[0][0],2,2,0,
    6150                                         Bprime,2,numnodes,0,
    6151                                         &Ke->values[0],1);
    6152 
    6153                 if(stabilization==2){
    6154                         /*Streamline upwinding*/
    6155                         vel=sqrt(vx*vx+vy*vy)+1.e-8;
    6156                         D[0][0]=h/(2*vel)*vx*vx;
    6157                         D[1][0]=h/(2*vel)*vy*vx;
    6158                         D[0][1]=h/(2*vel)*vx*vy;
    6159                         D[1][1]=h/(2*vel)*vy*vy;
    6160                 }
    6161                 else if(stabilization==1){
    6162                         /*SSA*/
    6163                         vxaverage_input->GetInputAverage(&vx);
    6164                         vyaverage_input->GetInputAverage(&vy);
    6165                         D[0][0]=h/2.0*fabs(vx);
    6166                         D[0][1]=0.;
    6167                         D[1][0]=0.;
    6168                         D[1][1]=h/2.0*fabs(vy);
    6169                 }
    6170                 if(stabilization==1 || stabilization==2){
    6171                         D[0][0]=D_scalar*D[0][0];
    6172                         D[1][0]=D_scalar*D[1][0];
    6173                         D[0][1]=D_scalar*D[0][1];
    6174                         D[1][1]=D_scalar*D[1][1];
    6175                         TripleMultiply(Bprime,2,numnodes,1,
    6176                                                 &D[0][0],2,2,0,
    6177                                                 Bprime,2,numnodes,0,
    6178                                                 &Ke->values[0],1);
    6179                 }
    6180         }
    6181 
    6182         /*Clean up and return*/
    6183         xDelete<IssmDouble>(basis);
    6184         xDelete<IssmDouble>(B);
    6185         xDelete<IssmDouble>(Bprime);
    6186         delete gauss;
    6187         return Ke;
    6188 }
    6189 /*}}}*/
    6190 /*FUNCTION Tria::CreateKMatrixMasstransport_DG {{{*/
    6191 ElementMatrix* Tria::CreateKMatrixMasstransport_DG(void){
    6192 
    6193         /*Intermediaries */
    6194         int        meshtype;
    6195         IssmDouble xyz_list[NUMVERTICES][3];
    6196         IssmDouble Jdet,D_scalar,dt,vx,vy;
    6197 
    6198         /*Fetch number of nodes for this finite element*/
    6199         int numnodes = this->NumberofNodes();
    6200 
    6201         /*Initialize Element matrix and vectors*/
    6202         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    6203         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    6204         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    6205         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    6206         IssmDouble     D[2][2];
    6207 
    6208         /*Retrieve all inputs and parameters*/
    6209         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6210         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6211         this->parameters->FindParam(&meshtype,MeshTypeEnum);
    6212         Input* vxaverage_input=NULL;
    6213         Input* vyaverage_input=NULL;
    6214         if(meshtype==Mesh2DhorizontalEnum){
    6215                 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    6216                 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    6217         }
    6218         else{
    6219                 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    6220                 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
    6221         }
    6222 
    6223         /* Start  looping on the number of gaussian points: */
    6224         GaussTria* gauss=new GaussTria(2);
    6225         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6226 
    6227                 gauss->GaussPoint(ig);
    6228 
    6229                 vxaverage_input->GetInputValue(&vx,gauss);
    6230                 vyaverage_input->GetInputValue(&vy,gauss);
    6231 
    6232                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6233                 GetNodalFunctions(basis,gauss);
    6234 
    6235                 D_scalar=gauss->weight*Jdet;
    6236 
    6237                 TripleMultiply(basis,1,numnodes,1,
    6238                                         &D_scalar,1,1,0,
    6239                                         basis,1,numnodes,0,
    6240                                         &Ke->values[0],1);
    6241 
    6242                 /*WARNING: B and Bprime are inverted compared to usual masstransport!!!!*/
    6243                 GetBMasstransport(Bprime, &xyz_list[0][0], gauss);
    6244                 GetBprimeMasstransport(B, &xyz_list[0][0], gauss);
    6245 
    6246                 D_scalar=-dt*gauss->weight*Jdet;
    6247                 D[0][0]=D_scalar*vx;
    6248                 D[0][1]=0.;
    6249                 D[1][0]=0.;
    6250                 D[1][1]=D_scalar*vy;
    6251 
    6252                 TripleMultiply(B,2,numnodes,1,
    6253                                         &D[0][0],2,2,0,
    6254                                         Bprime,2,numnodes,0,
    6255                                         &Ke->values[0],1);
    6256         }
    6257 
    6258         /*Clean up and return*/
    6259         xDelete<IssmDouble>(basis);
    6260         xDelete<IssmDouble>(B);
    6261         xDelete<IssmDouble>(Bprime);
    6262         delete gauss;
    6263         return Ke;
    6264 }
    6265 /*}}}*/
    6266 /*FUNCTION Tria::CreateKMatrixFreeSurfaceTop {{{*/
    6267 ElementMatrix* Tria::CreateKMatrixFreeSurfaceTop(void){
    6268 
    6269         /*Intermediaries */
    6270         int        stabilization;
    6271         IssmDouble Jdet,D_scalar,dt,h;
    6272         IssmDouble vel,vx,vy;
    6273         IssmDouble xyz_list[NUMVERTICES][3];
    6274 
    6275         /*Fetch number of nodes for this finite element*/
    6276         int numnodes = this->NumberofNodes();
    6277 
    6278         /*Initialize Element matrix and vectors*/
    6279         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    6280         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    6281         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    6282         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    6283         IssmDouble     D[2][2];
    6284 
    6285         /*Retrieve all inputs and parameters*/
    6286         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6287         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6288         this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
    6289         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    6290         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    6291         h=sqrt(2*this->GetArea());
    6292 
    6293         /* Start  looping on the number of gaussian points: */
    6294         GaussTria *gauss=new GaussTria(2);
    6295         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6296 
    6297                 gauss->GaussPoint(ig);
    6298 
    6299                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6300                 GetNodalFunctions(basis,gauss);
    6301 
    6302                 vx_input->GetInputValue(&vx,gauss);
    6303                 vy_input->GetInputValue(&vy,gauss);
    6304 
    6305                 D_scalar=gauss->weight*Jdet;
    6306 
    6307                 TripleMultiply(basis,1,numnodes,1,
    6308                                         &D_scalar,1,1,0,
    6309                                         basis,1,numnodes,0,
    6310                                         &Ke->values[0],1);
    6311 
    6312                 GetBMasstransport(B,&xyz_list[0][0],gauss);
    6313                 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
    6314 
    6315                 D_scalar=dt*gauss->weight*Jdet;
    6316 
    6317                 D[0][0]=D_scalar*vx;
    6318                 D[0][1]=0.;
    6319                 D[1][0]=0.;
    6320                 D[1][1]=D_scalar*vy;
    6321                 TripleMultiply(B,2,numnodes,1,
    6322                                         &D[0][0],2,2,0,
    6323                                         Bprime,2,numnodes,0,
    6324                                         &Ke->values[0],1);
    6325 
    6326                 if(stabilization==2){
    6327                         /*Streamline upwinding*/
    6328                         vel=sqrt(vx*vx+vy*vy)+1.e-8;
    6329                         D[0][0]=h/(2*vel)*vx*vx;
    6330                         D[1][0]=h/(2*vel)*vy*vx;
    6331                         D[0][1]=h/(2*vel)*vx*vy;
    6332                         D[1][1]=h/(2*vel)*vy*vy;
    6333                 }
    6334                 else if(stabilization==1){
    6335                         /*SSA*/
    6336                         vx_input->GetInputAverage(&vx);
    6337                         vy_input->GetInputAverage(&vy);
    6338                         D[0][0]=h/2.0*fabs(vx);
    6339                         D[0][1]=0.;
    6340                         D[1][0]=0.;
    6341                         D[1][1]=h/2.0*fabs(vy);
    6342                 }
    6343                 if(stabilization==1 || stabilization==2){
    6344                         D[0][0]=D_scalar*D[0][0];
    6345                         D[1][0]=D_scalar*D[1][0];
    6346                         D[0][1]=D_scalar*D[0][1];
    6347                         D[1][1]=D_scalar*D[1][1];
    6348                         TripleMultiply(Bprime,2,numnodes,1,
    6349                                                 &D[0][0],2,2,0,
    6350                                                 Bprime,2,numnodes,0,
    6351                                                 &Ke->values[0],1);
    6352                 }
    6353         }
    6354 
    6355         /*Clean up and return*/
    6356         xDelete<IssmDouble>(basis);
    6357         xDelete<IssmDouble>(B);
    6358         xDelete<IssmDouble>(Bprime);
    6359         delete gauss;
    6360         return Ke;
    6361 }
    6362 /*}}}*/
    6363 /*FUNCTION Tria::CreateKMatrixFreeSurfaceTop1D {{{*/
    6364 ElementMatrix* Tria::CreateKMatrixFreeSurfaceTop1D(void){
    6365 
    6366         if(!HasEdgeOnSurface()) return NULL;
    6367 
    6368         int index1,index2;
    6369         this->EdgeOnSurfaceIndices(&index1,&index2);
    6370 
    6371         Seg* seg=(Seg*)SpawnSeg(index1,index2);
    6372         ElementMatrix* Ke=seg->CreateKMatrixFreeSurfaceTop();
    6373         delete seg->material; delete seg;
    6374 
    6375         /*clean up and return*/
    6376         return Ke;
    6377 }
    6378 /*}}}*/
    6379 /*FUNCTION Tria::CreateKMatrixFreeSurfaceBase {{{*/
    6380 ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase(void){
    6381 
    6382         /*Intermediaries */
    6383         int        stabilization;
    6384         IssmDouble Jdet,D_scalar,dt,h;
    6385         IssmDouble vel,vx,vy;
    6386         IssmDouble xyz_list[NUMVERTICES][3];
    6387 
    6388         /*Fetch number of nodes for this finite element*/
    6389         int numnodes = this->NumberofNodes();
    6390 
    6391         /*Initialize Element matrix and vectors*/
    6392         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    6393         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    6394         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    6395         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    6396         IssmDouble     D[2][2];
    6397 
    6398         /*Retrieve all inputs and parameters*/
    6399         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6400         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6401         this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
    6402         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    6403         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    6404         h=sqrt(2*this->GetArea());
    6405 
    6406         /* Start  looping on the number of gaussian points: */
    6407         GaussTria *gauss=new GaussTria(2);
    6408         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6409 
    6410                 gauss->GaussPoint(ig);
    6411 
    6412                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6413                 GetNodalFunctions(basis,gauss);
    6414 
    6415                 vx_input->GetInputValue(&vx,gauss);
    6416                 vy_input->GetInputValue(&vy,gauss);
    6417 
    6418                 D_scalar=gauss->weight*Jdet;
    6419 
    6420                 TripleMultiply(basis,1,numnodes,1,
    6421                                         &D_scalar,1,1,0,
    6422                                         basis,1,numnodes,0,
    6423                                         &Ke->values[0],1);
    6424 
    6425                 GetBMasstransport(B,&xyz_list[0][0],gauss);
    6426                 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
    6427 
    6428                 D_scalar=dt*gauss->weight*Jdet;
    6429 
    6430                 D[0][0]=D_scalar*vx;
    6431                 D[0][1]=0.;
    6432                 D[1][0]=0.;
    6433                 D[1][1]=D_scalar*vy;
    6434                 TripleMultiply(B,2,numnodes,1,
    6435                                         &D[0][0],2,2,0,
    6436                                         Bprime,2,numnodes,0,
    6437                                         &Ke->values[0],1);
    6438 
    6439                 if(stabilization==2){
    6440                         /*Streamline upwinding*/
    6441                         vel=sqrt(vx*vx+vy*vy)+1.e-8;
    6442                         D[0][0]=h/(2*vel)*vx*vx;
    6443                         D[1][0]=h/(2*vel)*vy*vx;
    6444                         D[0][1]=h/(2*vel)*vx*vy;
    6445                         D[1][1]=h/(2*vel)*vy*vy;
    6446                 }
    6447                 else if(stabilization==1){
    6448                         /*SSA*/
    6449                         vx_input->GetInputAverage(&vx);
    6450                         vy_input->GetInputAverage(&vy);
    6451                         D[0][0]=h/2.0*fabs(vx);
    6452                         D[0][1]=0.;
    6453                         D[1][0]=0.;
    6454                         D[1][1]=h/2.0*fabs(vy);
    6455                 }
    6456                 if(stabilization==1 || stabilization==2){
    6457                         D[0][0]=D_scalar*D[0][0];
    6458                         D[1][0]=D_scalar*D[1][0];
    6459                         D[0][1]=D_scalar*D[0][1];
    6460                         D[1][1]=D_scalar*D[1][1];
    6461                         TripleMultiply(Bprime,2,numnodes,1,
    6462                                                 &D[0][0],2,2,0,
    6463                                                 Bprime,2,numnodes,0,
    6464                                                 &Ke->values[0],1);
    6465                 }
    6466         }
    6467 
    6468         /*Clean up and return*/
    6469         xDelete<IssmDouble>(basis);
    6470         xDelete<IssmDouble>(B);
    6471         xDelete<IssmDouble>(Bprime);
    6472         delete gauss;
    6473         return Ke;
    6474 }
    6475 /*}}}*/
    6476 /*FUNCTION Tria::CreateKMatrixFreeSurfaceBase1D {{{*/
    6477 ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase1D(void){
    6478 
    6479         if(!HasEdgeOnBed()) return NULL;
    6480 
    6481         int index1,index2;
    6482         this->EdgeOnBedIndices(&index1,&index2);
    6483 
    6484         Seg* seg=(Seg*)SpawnSeg(index1,index2);
    6485         ElementMatrix* Ke=seg->CreateKMatrixFreeSurfaceBase();
    6486         delete seg->material; delete seg;
    6487 
    6488         /*clean up and return*/
    6489         return Ke;
    6490 }
    6491 /*}}}*/
    6492 #endif
    6493 
    64944907#ifdef _HAVE_DAMAGE_
    6495 /*FUNCTION Tria::CreateKMatrixDamageEvolution {{{*/
    6496 ElementMatrix* Tria::CreateKMatrixDamageEvolution(void){
    6497 
    6498         switch(GetElementType()){
    6499                 case P1Enum: case P2Enum:
    6500                         return CreateKMatrixDamageEvolution_CG();
    6501                 case P1DGEnum:
    6502                         _error_("DG not implemented yet!");break;
    6503                 default:
    6504                         _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
    6505         }
    6506 }
    6507 /*}}}*/
    6508 /*FUNCTION Tria::CreateKMatrixDamageEvolution_CG {{{*/
    6509 ElementMatrix* Tria::CreateKMatrixDamageEvolution_CG(void){
    6510 
    6511         /*Intermediaries */
    6512         int        stabilization;
    6513         int        meshtype;
    6514         IssmDouble Jdet,D_scalar,dt,h;
    6515         IssmDouble vel,vx,vy,dvxdx,dvydy;
    6516         IssmDouble dvx[2],dvy[2];
    6517         IssmDouble xyz_list[NUMVERTICES][3];
    6518 
    6519         /*Fetch number of nodes for this finite element*/
    6520         int numnodes = this->NumberofNodes();
    6521 
    6522         /*Initialize Element matrix and vectors*/
    6523         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    6524         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    6525         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    6526         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    6527         IssmDouble     D[2][2];
    6528 
    6529         /*Retrieve all inputs and parameters*/
    6530         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6531         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6532         this->parameters->FindParam(&meshtype,MeshTypeEnum);
    6533         this->parameters->FindParam(&stabilization,DamageStabilizationEnum);
    6534         Input* vxaverage_input=NULL;
    6535         Input* vyaverage_input=NULL;
    6536         if(meshtype==Mesh2DhorizontalEnum){
    6537                 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    6538                 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    6539         }
    6540         else{
    6541                 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    6542                 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
    6543         }
    6544         h=sqrt(2*this->GetArea());
    6545 
    6546         /* Start  looping on the number of gaussian points: */
    6547         GaussTria *gauss=new GaussTria(2);
    6548         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6549 
    6550                 gauss->GaussPoint(ig);
    6551 
    6552                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6553                 GetNodalFunctions(basis,gauss);
    6554 
    6555                 vxaverage_input->GetInputValue(&vx,gauss);
    6556                 vyaverage_input->GetInputValue(&vy,gauss);
    6557                 vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    6558                 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    6559 
    6560                 D_scalar=gauss->weight*Jdet;
    6561 
    6562                 TripleMultiply(basis,1,numnodes,1,
    6563                                         &D_scalar,1,1,0,
    6564                                         basis,1,numnodes,0,
    6565                                         &Ke->values[0],1);
    6566                 GetBMasstransport(B,&xyz_list[0][0],gauss);
    6567                 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
    6568 
    6569                 dvxdx=dvx[0];
    6570                 dvydy=dvy[1];
    6571                 D_scalar=dt*gauss->weight*Jdet;
    6572 
    6573                 D[0][0]=D_scalar*dvxdx;
    6574                 D[0][1]=0.;
    6575                 D[1][0]=0.;
    6576                 D[1][1]=D_scalar*dvydy;
    6577                 TripleMultiply(B,2,numnodes,1,
    6578                                         &D[0][0],2,2,0,
    6579                                         B,2,numnodes,0,
    6580                                         &Ke->values[0],1);
    6581 
    6582                 D[0][0]=D_scalar*vx;
    6583                 D[1][1]=D_scalar*vy;
    6584                 TripleMultiply(B,2,numnodes,1,
    6585                                         &D[0][0],2,2,0,
    6586                                         Bprime,2,numnodes,0,
    6587                                         &Ke->values[0],1);
    6588 
    6589                 if(stabilization==2){
    6590                         /*Streamline upwinding*/
    6591                         vel=sqrt(vx*vx+vy*vy)+1.e-8;
    6592                         D[0][0]=h/(2*vel)*vx*vx;
    6593                         D[1][0]=h/(2*vel)*vy*vx;
    6594                         D[0][1]=h/(2*vel)*vx*vy;
    6595                         D[1][1]=h/(2*vel)*vy*vy;
    6596                 }
    6597                 else if(stabilization==1){
    6598                         /*SSA*/
    6599                         vxaverage_input->GetInputAverage(&vx);
    6600                         vyaverage_input->GetInputAverage(&vy);
    6601                         D[0][0]=h/2.0*fabs(vx);
    6602                         D[0][1]=0.;
    6603                         D[1][0]=0.;
    6604                         D[1][1]=h/2.0*fabs(vy);
    6605                 }
    6606                 if(stabilization==1 || stabilization==2){
    6607                         D[0][0]=D_scalar*D[0][0];
    6608                         D[1][0]=D_scalar*D[1][0];
    6609                         D[0][1]=D_scalar*D[0][1];
    6610                         D[1][1]=D_scalar*D[1][1];
    6611                         TripleMultiply(Bprime,2,numnodes,1,
    6612                                                 &D[0][0],2,2,0,
    6613                                                 Bprime,2,numnodes,0,
    6614                                                 &Ke->values[0],1);
    6615                 }
    6616         }
    6617 
    6618         /*Clean up and return*/
    6619         xDelete<IssmDouble>(basis);
    6620         xDelete<IssmDouble>(B);
    6621         xDelete<IssmDouble>(Bprime);
    6622         delete gauss;
    6623         return Ke;
    6624 }
    6625 /*}}}*/
    66264908/*FUNCTION Tria::DamageEvolutionF{{{*/
    66274909void Tria::DamageEvolutionF(IssmDouble* f){
    … …  
    68195101        }
    68205102
    6821 }
    6822 /*}}}*/
    6823 #endif
    6824 
    6825 #ifdef _HAVE_BALANCED_
    6826 /*FUNCTION Tria::CreateKMatrixBalancethickness {{{*/
    6827 ElementMatrix* Tria::CreateKMatrixBalancethickness(void){
    6828 
    6829         switch(GetElementType()){
    6830                 case P1Enum:
    6831                         return CreateKMatrixBalancethickness_CG();
    6832                 case P1DGEnum:
    6833                         return CreateKMatrixBalancethickness_DG();
    6834                 default:
    6835                         _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
    6836         }
    6837 
    6838 }
    6839 /*}}}*/
    6840 /*FUNCTION Tria::CreateKMatrixBalancethickness_CG {{{*/
    6841 ElementMatrix* Tria::CreateKMatrixBalancethickness_CG(void){
    6842 
    6843         /*Intermediaries */
    6844         int        stabilization,meshtype;
    6845         IssmDouble Jdet,vx,vy,dvxdx,dvydy,vel,h;
    6846         IssmDouble D_scalar;
    6847         IssmDouble dvx[2],dvy[2];
    6848         IssmDouble xyz_list[NUMVERTICES][3];
    6849 
    6850         /*Fetch number of nodes and dof for this finite element*/
    6851         int numnodes = this->NumberofNodes();
    6852 
    6853         /*Initialize Element matrix and vectors*/
    6854         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    6855         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    6856         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    6857         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    6858         IssmDouble     D[2][2];
    6859 
    6860         /*Retrieve all Inputs and parameters: */
    6861         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6862         this->parameters->FindParam(&stabilization,BalancethicknessStabilizationEnum);
    6863         this->parameters->FindParam(&meshtype,MeshTypeEnum);
    6864         Input* vxaverage_input=NULL;
    6865         Input* vyaverage_input=NULL;
    6866         if(meshtype==Mesh2DhorizontalEnum){
    6867                 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
    6868                 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
    6869         }
    6870         else{
    6871                 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    6872                 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
    6873         }
    6874         h=sqrt(2.*this->GetArea());
    6875 
    6876         /*Start looping on the number of gaussian points:*/
    6877         GaussTria* gauss=new GaussTria(2);
    6878         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6879 
    6880                 gauss->GaussPoint(ig);
    6881 
    6882                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6883                 GetBMasstransport(B,&xyz_list[0][0],gauss);
    6884                 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
    6885 
    6886                 vxaverage_input->GetInputValue(&vx,gauss);
    6887                 vyaverage_input->GetInputValue(&vy,gauss);
    6888                 vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    6889                 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    6890 
    6891                 dvxdx=dvx[0];
    6892                 dvydy=dvy[1];
    6893                 D_scalar=gauss->weight*Jdet;
    6894 
    6895                 D[0][0]=D_scalar*dvxdx;
    6896                 D[0][1]=0.;
    6897                 D[1][0]=0.;
    6898                 D[1][1]=D_scalar*dvydy;
    6899                 TripleMultiply(B,2,numnodes,1,
    6900                                         &D[0][0],2,2,0,
    6901                                         B,2,numnodes,0,
    6902                                         &Ke->values[0],1);
    6903 
    6904                 D[0][0]=D_scalar*vx;
    6905                 D[1][1]=D_scalar*vy;
    6906                 TripleMultiply(B,2,numnodes,1,
    6907                                         &D[0][0],2,2,0,
    6908                                         Bprime,2,numnodes,0,
    6909                                         &Ke->values[0],1);
    6910 
    6911                 if(stabilization==1){
    6912                         /*Streamline upwinding*/
    6913                         vel=sqrt(vx*vx+vy*vy);
    6914                         D[0][0]=h/(2*vel)*vx*vx;
    6915                         D[1][0]=h/(2*vel)*vy*vx;
    6916                         D[0][1]=h/(2*vel)*vx*vy;
    6917                         D[1][1]=h/(2*vel)*vy*vy;
    6918                 }
    6919                 else if(stabilization==2){
    6920                         /*SSA*/
    6921                         vxaverage_input->GetInputAverage(&vx);
    6922                         vyaverage_input->GetInputAverage(&vy);
    6923                         D[0][0]=h/2.0*fabs(vx);
    6924                         D[0][1]=0.;
    6925                         D[1][0]=0.;
    6926                         D[1][1]=h/2.0*fabs(vy);
    6927                 }
    6928                 if(stabilization==1 || stabilization==2){
    6929                         D[0][0]=D_scalar*D[0][0];
    6930                         D[1][0]=D_scalar*D[1][0];
    6931                         D[0][1]=D_scalar*D[0][1];
    6932                         D[1][1]=D_scalar*D[1][1];
    6933                         TripleMultiply(Bprime,2,numnodes,1,
    6934                                                 &D[0][0],2,2,0,
    6935                                                 Bprime,2,numnodes,0,
    6936                                                 &Ke->values[0],1);
    6937                 }
    6938         }
    6939 
    6940         /*Clean up and return*/
    6941         xDelete<IssmDouble>(basis);
    6942         xDelete<IssmDouble>(B);
    6943         xDelete<IssmDouble>(Bprime);
    6944         delete gauss;
    6945         return Ke;
    6946 }
    6947 /*}}}*/
    6948 /*FUNCTION Tria::CreateKMatrixBalancethickness_DG {{{*/
    6949 ElementMatrix* Tria::CreateKMatrixBalancethickness_DG(void){
    6950 
    6951         /*Intermediaries*/
    6952         IssmDouble vx,vy,D_scalar,Jdet;
    6953         IssmDouble xyz_list[NUMVERTICES][3];
    6954 
    6955         /*Fetch number of nodes for this finite element*/
    6956         int numnodes = this->NumberofNodes();
    6957 
    6958         /*Initialize Element matrix and vectors*/
    6959         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    6960         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    6961         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    6962         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    6963         IssmDouble     D[2][2];
    6964 
    6965         /*Retrieve all inputs and parameters*/
    6966         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6967         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    6968         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    6969 
    6970         /*Start looping on the number of gaussian points:*/
    6971         GaussTria* gauss=new GaussTria(2);
    6972         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6973 
    6974                 gauss->GaussPoint(ig);
    6975 
    6976                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6977                 /*WARNING: B and Bprime are inverted compared to usual masstransport!!!!*/
    6978                 GetBMasstransport(Bprime,&xyz_list[0][0],gauss);
    6979                 GetBprimeMasstransport(B,&xyz_list[0][0],gauss);
    6980 
    6981                 vx_input->GetInputValue(&vx,gauss);
    6982                 vy_input->GetInputValue(&vy,gauss);
    6983 
    6984                 D_scalar=-gauss->weight*Jdet;
    6985                 D[0][0]=D_scalar*vx;
    6986                 D[0][1]=0.;
    6987                 D[1][0]=0.;
    6988                 D[1][1]=D_scalar*vy;
    6989 
    6990                 TripleMultiply(B,2,numnodes,1,
    6991                                         &D[0][0],2,2,0,
    6992                                         Bprime,2,numnodes,0,
    6993                                         &Ke->values[0],1);
    6994         }
    6995 
    6996         /*Clean up and return*/
    6997         xDelete<IssmDouble>(basis);
    6998         xDelete<IssmDouble>(B);
    6999         xDelete<IssmDouble>(Bprime);
    7000         delete gauss;
    7001         return Ke;
    7002 }
    7003 /*}}}*/
    7004 /*FUNCTION Tria::CreateKMatrixBalancevelocity{{{*/
    7005 ElementMatrix* Tria::CreateKMatrixBalancevelocity(void){
    7006 
    7007         /*Intermediaries */
    7008         IssmDouble xyz_list[NUMVERTICES][3];
    7009         IssmDouble dhdt_g,mb_g,ms_g,Jdet;
    7010         IssmDouble h,gamma,thickness;
    7011         IssmDouble hnx,hny,dhnx[2],dhny[2];
    7012 
    7013         /*Fetch number of nodes and dof for this finite element*/
    7014         int numnodes = this->NumberofNodes();
    7015 
    7016         /*Initialize Element matrix and vectors*/
    7017         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    7018         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    7019         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    7020         IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    7021         IssmDouble*    HNx    = xNew<IssmDouble>(numnodes);
    7022         IssmDouble*    HNy    = xNew<IssmDouble>(numnodes);
    7023         IssmDouble*    H      = xNew<IssmDouble>(numnodes);
    7024         IssmDouble*    Nx     = xNew<IssmDouble>(numnodes);
    7025         IssmDouble*    Ny     = xNew<IssmDouble>(numnodes);
    7026 
    7027         /*Retrieve all Inputs and parameters: */
    7028         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7029         Input* H_input =inputs->GetInput(ThicknessEnum); _assert_(H_input);
    7030         h=sqrt(2.*this->GetArea());
    7031 
    7032         /*Get vector N for all nodes and build HNx and HNy*/
    7033         GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
    7034         GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
    7035         GetInputListOnNodes(H,ThicknessEnum);
    7036         for(int i=0;i<numnodes;i++){
    7037                 IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
    7038                 HNx[i] = -H[i]*Nx[i]/norm;
    7039                 HNy[i] = -H[i]*Ny[i]/norm;
    7040         }
    7041 
    7042         /*Start looping on the number of gaussian points:*/
    7043         GaussTria* gauss=new GaussTria(2);
    7044         for(int ig=gauss->begin();ig<gauss->end();ig++){
    7045 
    7046                 gauss->GaussPoint(ig);
    7047 
    7048                 H_input->GetInputValue(&thickness,gauss);
    7049                 if(thickness<50.) thickness=50.;
    7050                 TriaRef::GetInputDerivativeValue(&dhnx[0],HNx,&xyz_list[0][0],gauss);
    7051                 TriaRef::GetInputDerivativeValue(&dhny[0],HNy,&xyz_list[0][0],gauss);
    7052                 TriaRef::GetInputValue(&hnx,HNx,gauss);
    7053                 TriaRef::GetInputValue(&hny,HNy,gauss);
    7054 
    7055                 gamma=h/(2.*thickness+1.e-10);
    7056 
    7057                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7058                 GetNodalFunctions(basis,gauss);
    7059                 GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
    7060 
    7061                 for(int i=0;i<numnodes;i++){
    7062                         for(int j=0;j<numnodes;j++){
    7063                                 Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
    7064                                                         (basis[i]+gamma*(basis[i]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+i]*hnx + dbasis[1*numnodes+i]*hny))*
    7065                                                         (basis[j]*(dhnx[0]+dhny[1])  + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)
    7066                                                         );
    7067                         }
    7068                 }
    7069         }
    7070 
    7071         /*Clean up and return*/
    7072         xDelete<IssmDouble>(basis);
    7073         xDelete<IssmDouble>(dbasis);
    7074         xDelete<IssmDouble>(H);
    7075         xDelete<IssmDouble>(Nx);
    7076         xDelete<IssmDouble>(Ny);
    7077         xDelete<IssmDouble>(HNx);
    7078         xDelete<IssmDouble>(HNy);
    7079         xDelete<IssmDouble>(B);
    7080         delete gauss;
    7081         return Ke;
    7082 }
    7083 /*}}}*/
    7084 /*FUNCTION Tria::CreateKMatrixSmoothedSlope {{{*/
    7085 ElementMatrix* Tria::CreateKMatrixSmoothedSlope(void){
    7086 
    7087         /* Intermediaries */
    7088         IssmDouble D_scalar,Jdet,thickness;
    7089         IssmDouble xyz_list[NUMVERTICES][3];
    7090         IssmDouble l=8.;
    7091 
    7092         /*Fetch number of nodes and dof for this finite element*/
    7093         int numnodes = this->NumberofNodes();
    7094 
    7095         /*Initialize Element matrix and vectors*/
    7096         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
    7097         IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    7098         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    7099 
    7100         /*Retrieve all inputs and parameters*/
    7101         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7102         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    7103 
    7104         /* Start looping on the number of gaussian points: */
    7105         GaussTria* gauss=new GaussTria(2);
    7106         for(int ig=gauss->begin();ig<gauss->end();ig++){
    7107 
    7108                 gauss->GaussPoint(ig);
    7109 
    7110                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7111                 thickness_input->GetInputValue(&thickness,gauss);
    7112                 if(thickness<50.) thickness=50.;
    7113 
    7114                 GetNodalFunctions(basis,gauss);
    7115                 GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
    7116 
    7117                 for(int i=0;i<numnodes;i++){
    7118                         for(int j=0;j<numnodes;j++){
    7119                                 Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
    7120                                                         basis[i]*basis[j]
    7121                                                         +(l*thickness)*(l*thickness)*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
    7122                                                         );
    7123                         }
    7124                 }
    7125         }
    7126 
    7127         /*Clean up and return*/
    7128         delete gauss;
    7129         xDelete<IssmDouble>(dbasis);
    7130         xDelete<IssmDouble>(basis);
    7131         return Ke;
    71325103}
    71335104/*}}}*/

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

    r16982 r16993  
    6868                void        SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
    6969                void        CreateDVector(Vector<IssmDouble>* df);
    70                 void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
    7170                void        Delta18oParameterization(void);
    7271                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
    … …  
    202201                void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
    203202                void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);
    204                 ElementMatrix* CreateKMatrix(void);
    205                 ElementMatrix* CreateKMatrixBalancethickness(void);
    206                 ElementMatrix* CreateKMatrixBalancethickness_DG(void);
    207                 ElementMatrix* CreateKMatrixBalancethickness_CG(void);
    208                 ElementMatrix* CreateKMatrixBalancevelocity(void);
    209                 ElementMatrix* CreateKMatrixSmoothedSlope(void);
    210                 ElementMatrix* CreateKMatrixMelting(void);
    211                 ElementMatrix* CreateKMatrixMasstransport(void);
    212                 ElementMatrix* CreateKMatrixMasstransport_CG(void);
    213                 ElementMatrix* CreateKMatrixMasstransport_DG(void);
    214                 ElementMatrix* CreateKMatrixExtrusion(void);
    215                 ElementMatrix* CreateKMatrixExtrusionVolume(void);
    216                 ElementMatrix* CreateKMatrixExtrusionSurface(void);
    217                 ElementMatrix* CreateKMatrixExtrusionBed(void);
    218                 ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
    219                 ElementMatrix* CreateKMatrixFreeSurfaceTop1D(void);
    220                 ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
    221                 ElementMatrix* CreateKMatrixFreeSurfaceBase1D(void);
    222                 ElementMatrix* CreateMassMatrix(void);
    223                 ElementMatrix* CreateBasalMassMatrix(void);
    224203                IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
    225204                IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
    … …  
    283262
    284263                #ifdef _HAVE_STRESSBALANCE_
    285                 ElementMatrix* CreateKMatrixStressbalanceSSA(void);
    286                 ElementMatrix* CreateKMatrixStressbalanceSSAViscous(void);
    287                 ElementMatrix* CreateKMatrixStressbalanceSSAFriction(void);
    288                 ElementMatrix* CreateKMatrixStressbalanceSIA(void);
    289                 ElementMatrix* CreateKMatrixStressbalanceFS(void);
    290                 ElementMatrix* CreateKMatrixStressbalanceFSViscous(void);
    291                 ElementMatrix* CreateKMatrixStressbalanceFSFriction(void);
    292264                void           PVectorGLSstabilization(ElementVector* pe);
    293                 ElementMatrix* CreateJacobianStressbalanceSSA(void);
    294265                IssmDouble     GetYcoord(GaussTria* gauss);
    295                 #endif
    296 
    297                 #ifdef _HAVE_CONTROL_
    298                 ElementMatrix* CreateKMatrixAdjointBalancethickness(void);
    299                 ElementMatrix* CreateKMatrixAdjointSSA(void);
    300266                #endif
    301267
    … …  
    309275
    310276                #ifdef _HAVE_HYDROLOGY_
    311                 ElementMatrix* CreateKMatrixHydrologyShreve(void);
    312                 ElementMatrix* CreateKMatrixHydrologyDCInefficient(void);
    313                 ElementMatrix* CreateKMatrixHydrologyDCEfficient(void);
    314277                ElementMatrix* CreateEPLDomainMassMatrix(void);
    315278                void           CreateHydrologyWaterVelocityInput(void);
    … …  
    324287
    325288                #ifdef _HAVE_DAMAGE_
    326                 ElementMatrix* CreateKMatrixDamageEvolution(void);
    327                 ElementMatrix* CreateKMatrixDamageEvolution_CG(void);
    328289                void           DamageEvolutionF(IssmDouble* flist);
    329290                #endif

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

    r16164 r16993  
    1010void CreateJacobianMatrixx(Matrix<IssmDouble>** pJff,FemModel* femmodel,IssmDouble kmax){
    1111
    12         int      i,connectivity;
    13         int      configuration_type;
     12        int      i;
     13        int      configuration_type,analysisenum;
    1414        Element *element = NULL;
    1515        Load    *load    = NULL;
    … …  
    2121        /*Recover some parameters*/
    2222        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    23         femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
     23        femmodel->parameters->FindParam(&analysisenum,AnalysisTypeEnum);
     24        Analysis* analysis = EnumToAnalysis(analysisenum);
    2425
    2526        /*Initialize Jacobian Matrix*/
    … …  
    2930        for(i=0;i<femmodel->elements->Size();i++){
    3031                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    31                 element->CreateJacobianMatrix(Jff);
     32                ElementMatrix* Je = analysis->CreateJacobianMatrix(element);
     33                if(Je) Je->AddToGlobal(Jff);
     34                delete Je;
    3235        }
    3336        for (i=0;i<femmodel->loads->Size();i++){
    … …  
    3942
    4043        /*Assign output pointer*/
     44        delete analysis;
    4145        *pJff=Jff;
    4246

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

    r16126 r16993  
    66
    77#include "../../classes/classes.h"
     8#include "../../analyses/analyses.h"
    89
    910/* local prototypes: */

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

    r16982 r16993  
    4141                for (i=0;i<femmodel->elements->Size();i++){
    4242                        element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    43                         ElementMatrix* Ke = element->CreateKMatrix();
    44                         //ElementVector* pe = element->CreatePVector();
     43                        ElementMatrix* Ke = analysis->CreateKMatrix(element);
    4544                        ElementVector* pe = analysis->CreatePVector(element);
    4645                        element->ReduceMatrices(Ke,pe);
    … …  
    7372        for (i=0;i<femmodel->elements->Size();i++){
    7473                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    75                 ElementMatrix* Ke = element->CreateKMatrix();
    76                 //ElementMatrix* Ke = analysis->CreateKMatrix(element);
    77                 //ElementVector* pe = element->CreatePVector();
     74                ElementMatrix* Ke = analysis->CreateKMatrix(element);
    7875                ElementVector* pe = analysis->CreatePVector(element);
    7976                element->ReduceMatrices(Ke,pe);
Note: See TracChangeset for help on using the changeset viewer.