Changeset 16831


Ignore:
Timestamp:
11/19/13 13:14:35 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added FS KMatrix Shelf

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

Legend:

Unmodified
Added
Removed

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

    r16805 r16831  
    5656/*Finite Element Analysis*/
    5757ElementMatrix* L2ProjectionBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/
    58         _error_("not implemented yet");
     58
     59        /*Intermediaries*/
     60        int      meshtype;
     61        Element* basalelement;
     62
     63        /*Get basal element*/
     64        element->FindParam(&meshtype,MeshTypeEnum);
     65        switch(meshtype){
     66                case Mesh2DhorizontalEnum:
     67                        basalelement = element;
     68                        break;
     69                case Mesh3DEnum:
     70                        if(!element->IsOnBed()) return NULL;
     71                        basalelement = element->SpawnBasalElement();
     72                        break;
     73                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     74        }
     75
     76        /*Intermediaries */
     77        IssmDouble  D,Jdet;
     78        IssmDouble *xyz_list  = NULL;
     79
     80        /*Fetch number of nodes and dof for this finite element*/
     81        int numnodes = basalelement->GetNumberOfNodes();
     82
     83        /*Initialize Element vector*/
     84        ElementMatrix* Ke    = basalelement->NewElementMatrix(NoneApproximationEnum);
     85        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     86
     87        /*Retrieve all inputs and parameters*/
     88        basalelement->GetVerticesCoordinates(&xyz_list);
     89
     90        /* Start  looping on the number of gaussian points: */
     91        Gauss* gauss=basalelement->NewGauss(2);
     92        for(int ig=gauss->begin();ig<gauss->end();ig++){
     93                gauss->GaussPoint(ig);
     94
     95                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     96                basalelement->NodalFunctions(basis,gauss);
     97                D=gauss->weight*Jdet;
     98
     99                TripleMultiply(basis,1,numnodes,1,
     100                                        &D,1,1,0,
     101                                        basis,1,numnodes,0,
     102                                        &Ke->values[0],1);
     103        }
     104
     105        /*Clean up and return*/
     106        xDelete<IssmDouble>(xyz_list);
     107        xDelete<IssmDouble>(basis);
     108        delete gauss;
     109        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     110        return Ke;
    59111}/*}}}*/
    60112ElementVector* L2ProjectionBaseAnalysis::CreatePVector(Element* element){/*{{{*/

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

    r16829 r16831  
    820820                case HOApproximationEnum:
    821821                        return CreateKMatrixHO(element);
     822                case FSApproximationEnum:
     823                        return CreateKMatrixFS(element);
    822824                case NoneApproximationEnum:
    823825                        return NULL;
    … …  
    16581660
    16591661/*FS*/
     1662ElementMatrix* StressbalanceAnalysis::CreateKMatrixFS(Element* element){/*{{{*/
     1663
     1664        /*compute all stiffness matrices for this element*/
     1665        ElementMatrix* Ke1=CreateKMatrixFSViscous(element);
     1666        ElementMatrix* Ke2=CreateKMatrixFSFriction(element);
     1667        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     1668
     1669        /*clean-up and return*/
     1670        delete Ke1;
     1671        delete Ke2;
     1672        return Ke;
     1673}/*}}}*/
     1674ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscous(Element* element){/*{{{*/
     1675
     1676        /*Intermediaries*/
     1677        int         i,meshtype,dim,epssize;
     1678        IssmDouble  viscosity,FSreconditioning,Jdet;
     1679        IssmDouble *xyz_list = NULL;
     1680
     1681        /*Get problem dimension*/
     1682        element->FindParam(&meshtype,MeshTypeEnum);
     1683        switch(meshtype){
     1684                case Mesh2DverticalEnum: dim = 2; epssize = 3; break;
     1685                case Mesh3DEnum:         dim = 3; epssize = 6; break;
     1686                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1687        }
     1688
     1689        /*Fetch number of nodes and dof for this finite element*/
     1690        int vnumnodes = element->GetNumberOfNodesVelocity();
     1691        int pnumnodes = element->GetNumberOfNodesPressure();
     1692        int numdof    = vnumnodes*dim + pnumnodes;
     1693        int bsize     = epssize + 2;
     1694
     1695        /*Prepare coordinate system list*/
     1696        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     1697        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     1698        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     1699        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     1700
     1701        /*Initialize Element matrix and vectors*/
     1702        ElementMatrix* Ke     = element->NewElementMatrix(FSvelocityEnum);
     1703        IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
     1704        IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
     1705        IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
     1706
     1707        /*Retrieve all inputs and parameters*/
     1708        element->GetVerticesCoordinates(&xyz_list);
     1709        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     1710        Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
     1711        Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
     1712        Input* vz_input;
     1713        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
     1714
     1715        /* Start  looping on the number of gaussian points: */
     1716        Gauss* gauss = element->NewGauss(5);
     1717        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1718                gauss->GaussPoint(ig);
     1719
     1720                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1721                this->GetBFS(B,element,dim,xyz_list,gauss);
     1722                this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
     1723
     1724                element->ViscosityFS(&viscosity,xyz_list,gauss,vx_input,vy_input,vz_input);
     1725                for(i=0;i<epssize;i++)     D[i*bsize+i] = + 2.*viscosity*gauss->weight*Jdet;
     1726                for(i=epssize;i<bsize;i++) D[i*bsize+i] = - FSreconditioning*gauss->weight*Jdet;
     1727
     1728                TripleMultiply(B,bsize,numdof,1,
     1729                                        D,bsize,bsize,0,
     1730                                        Bprime,bsize,numdof,0,
     1731                                        &Ke->values[0],1);
     1732        }
     1733
     1734        /*Transform Coordinate System*/
     1735        element->TransformStiffnessMatrixCoord(Ke,cs_list);
     1736
     1737        /*Clean up and return*/
     1738        delete gauss;
     1739        xDelete<IssmDouble>(xyz_list);
     1740        xDelete<IssmDouble>(D);
     1741        xDelete<IssmDouble>(Bprime);
     1742        xDelete<IssmDouble>(B);
     1743        xDelete<int>(cs_list);
     1744        return Ke;
     1745}/*}}}*/
     1746ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
     1747        return NULL;
     1748}/*}}}*/
    16601749ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
    16611750
    … …  
    18161905
    18171906        return NULL;
     1907}/*}}}*/
     1908void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1909        /*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
     1910         * For node i, Bvi can be expressed in the actual coordinate system
     1911         * by:     Bvi=[ dphi/dx          0        ]
     1912         *                                       [   0           dphi/dy     ]
     1913         *                                       [ 1/2*dphi/dy    1/2*dphi/dx]
     1914         *                                       [   0             0         ]
     1915         *                                       [ dphi/dx         dphi/dy   ]
     1916         *
     1917         *         Bpi=[  0    ]
     1918         *                                      [  0    ]
     1919         *                                      [  0    ]
     1920         *                                      [ phi_p ]
     1921         *                                      [  0    ]
     1922         *
     1923         *      In 3d:
     1924         *         Bvi=[ dh/dx          0             0      ]
     1925         *                                      [   0           dh/dy           0      ]
     1926         *                                      [   0             0           dh/dz    ]
     1927         *                                      [ 1/2*dh/dy    1/2*dh/dx        0      ]
     1928         *                                      [ 1/2*dh/dz       0         1/2*dh/dx  ]
     1929         *                                      [   0          1/2*dh/dz    1/2*dh/dy  ]
     1930         *                                      [   0             0             0      ]
     1931         *                                      [ dh/dx         dh/dy         dh/dz    ]
     1932         *
     1933         *         Bpi=[ 0 ]
     1934         *                                      [ 0 ]
     1935         *                                      [ 0 ]
     1936         *                                      [ 0 ]
     1937         *                                      [ 0 ]
     1938         *                                      [ 0 ]
     1939         *                                      [ h ]
     1940         *                                      [ 0 ]
     1941         *      where phi is the finiteelement function for node i.
     1942         *      Same thing for Bb except the last column that does not exist.
     1943         */
     1944
     1945        /*Fetch number of nodes for this finite element*/
     1946        int pnumnodes = element->NumberofNodesPressure();
     1947        int vnumnodes = element->NumberofNodesVelocity();
     1948
     1949        /*Get nodal functions derivatives*/
     1950        IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
     1951        IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
     1952        element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     1953        element->NodalFunctionsPressure(pbasis,gauss);
     1954
     1955        /*Build B: */
     1956        if(dim==2){
     1957                for(int i=0;i<vnumnodes;i++){
     1958                        B[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
     1959                        B[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
     1960                        B[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
     1961                        B[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
     1962                        B[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
     1963                        B[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
     1964                        B[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = 0.;
     1965                        B[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = 0.;
     1966                        B[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = vdbasis[0*vnumnodes+i];
     1967                        B[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = vdbasis[1*vnumnodes+i];
     1968                }
     1969                for(int i=0;i<pnumnodes;i++){
     1970                        B[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
     1971                        B[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
     1972                        B[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
     1973                        B[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = pbasis[i];
     1974                        B[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = 0.;
     1975                }
     1976        }
     1977        else{
     1978                for(int i=0;i<vnumnodes;i++){
     1979                        B[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
     1980                        B[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
     1981                        B[(dim*vnumnodes+pnumnodes)*0+dim*i+2] = 0.;
     1982                        B[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
     1983                        B[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
     1984                        B[(dim*vnumnodes+pnumnodes)*1+dim*i+2] = 0.;
     1985                        B[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = 0.;
     1986                        B[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = 0.;
     1987                        B[(dim*vnumnodes+pnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
     1988                        B[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
     1989                        B[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
     1990                        B[(dim*vnumnodes+pnumnodes)*3+dim*i+2] = 0.;
     1991                        B[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = .5*vdbasis[2*vnumnodes+i];
     1992                        B[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = 0.;
     1993                        B[(dim*vnumnodes+pnumnodes)*4+dim*i+2] = .5*vdbasis[0*vnumnodes+i];
     1994                        B[(dim*vnumnodes+pnumnodes)*5+dim*i+0] = 0.;
     1995                        B[(dim*vnumnodes+pnumnodes)*5+dim*i+1] = .5*vdbasis[2*vnumnodes+i];
     1996                        B[(dim*vnumnodes+pnumnodes)*5+dim*i+2] = .5*vdbasis[1*vnumnodes+i];
     1997                        B[(dim*vnumnodes+pnumnodes)*6+dim*i+0] = 0.;
     1998                        B[(dim*vnumnodes+pnumnodes)*6+dim*i+1] = 0.;
     1999                        B[(dim*vnumnodes+pnumnodes)*6+dim*i+2] = 0.;
     2000                        B[(dim*vnumnodes+pnumnodes)*7+dim*i+0] = vdbasis[0*vnumnodes+i];
     2001                        B[(dim*vnumnodes+pnumnodes)*7+dim*i+1] = vdbasis[1*vnumnodes+i];
     2002                        B[(dim*vnumnodes+pnumnodes)*7+dim*i+2] = vdbasis[2*vnumnodes+i];
     2003                }
     2004                for(int i=0;i<pnumnodes;i++){
     2005                        B[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
     2006                        B[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
     2007                        B[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
     2008                        B[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = 0.;
     2009                        B[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = 0.;
     2010                        B[(dim*vnumnodes+pnumnodes)*5+(dim*vnumnodes)+i] = 0.;
     2011                        B[(dim*vnumnodes+pnumnodes)*6+(dim*vnumnodes)+i] = pbasis[i];
     2012                        B[(dim*vnumnodes+pnumnodes)*7+(dim*vnumnodes)+i] = 0.;
     2013                }
     2014        }
     2015
     2016        /*Clean up*/
     2017        xDelete<IssmDouble>(vdbasis);
     2018        xDelete<IssmDouble>(pbasis);
     2019}/*}}}*/
     2020void StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2021        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
     2022         *      For node i, Bi' can be expressed in the actual coordinate system
     2023         *      by:
     2024         *                      Bvi' = [  dphi/dx     0     ]
     2025         *                                       [     0      dphi/dy ]
     2026         *                                       [  dphi/dy   dphi/dx ]
     2027         *                                       [  dphi/dx   dphi/dy ]
     2028         *                                       [     0      0       ]
     2029         *
     2030         * by:    Bpi=[  0  ]
     2031         *                                      [  0  ]
     2032         *                                      [  0  ]
     2033         *                                      [  0  ]
     2034         *                                      [ phi ]
     2035         *
     2036         *      In 3d
     2037         *         Bvi=[ dh/dx          0             0      ]
     2038         *                                      [   0           dh/dy           0      ]
     2039         *                                      [   0             0           dh/dz    ]
     2040         *                                      [ 1/2*dh/dy    1/2*dh/dx        0      ]
     2041         *                                      [ 1/2*dh/dz       0         1/2*dh/dx  ]
     2042         *                                      [   0          1/2*dh/dz    1/2*dh/dy  ]
     2043         *                                      [   0             0             0      ]
     2044         *                                      [ dh/dx         dh/dy         dh/dz    ]
     2045         *
     2046         *         Bpi=[ 0 ]
     2047         *                                      [ 0 ]
     2048         *                                      [ 0 ]
     2049         *                                      [ 0 ]
     2050         *                                      [ 0 ]
     2051         *                                      [ 0 ]
     2052         *                                      [ h ]
     2053         *                                      [ 0 ]
     2054         *      where phi is the finiteelement function for node i.
     2055         *      In 3d:
     2056         */
     2057
     2058        /*Fetch number of nodes for this finite element*/
     2059        int pnumnodes = element->NumberofNodesPressure();
     2060        int vnumnodes = element->NumberofNodesVelocity();
     2061
     2062        /*Get nodal functions derivatives*/
     2063        IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
     2064        IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
     2065        element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     2066        element->NodalFunctionsPressure(pbasis,gauss);
     2067
     2068        /*Build B_prime: */
     2069        if(dim==2){
     2070                for(int i=0;i<vnumnodes;i++){
     2071                        Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
     2072                        Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
     2073                        Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
     2074                        Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
     2075                        Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = vdbasis[1*vnumnodes+i];
     2076                        Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = vdbasis[0*vnumnodes+i];
     2077                        Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = vdbasis[0*vnumnodes+i];
     2078                        Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = vdbasis[1*vnumnodes+i];
     2079                        Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = 0.;
     2080                        Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = 0.;
     2081                }
     2082                for(int i=0;i<pnumnodes;i++){
     2083                        Bprime[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
     2084                        Bprime[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
     2085                        Bprime[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
     2086                        Bprime[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = 0.;
     2087                        Bprime[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = pbasis[i];
     2088                }
     2089        }
     2090        else{
     2091                for(int i=0;i<vnumnodes;i++){
     2092                        Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
     2093                        Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
     2094                        Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+2] = 0.;
     2095                        Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
     2096                        Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
     2097                        Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+2] = 0.;
     2098                        Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = 0.;
     2099                        Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = 0.;
     2100                        Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
     2101                        Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = vdbasis[1*vnumnodes+i];
     2102                        Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = vdbasis[0*vnumnodes+i];
     2103                        Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+2] = 0.;
     2104                        Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = vdbasis[2*vnumnodes+i];
     2105                        Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = 0.;
     2106                        Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+2] = vdbasis[0*vnumnodes+i];
     2107                        Bprime[(dim*vnumnodes+pnumnodes)*5+dim*i+0] = 0.;
     2108                        Bprime[(dim*vnumnodes+pnumnodes)*5+dim*i+1] = vdbasis[2*vnumnodes+i];
     2109                        Bprime[(dim*vnumnodes+pnumnodes)*5+dim*i+2] = vdbasis[1*vnumnodes+i];
     2110                        Bprime[(dim*vnumnodes+pnumnodes)*6+dim*i+0] = vdbasis[0*vnumnodes+i];
     2111                        Bprime[(dim*vnumnodes+pnumnodes)*6+dim*i+1] = vdbasis[1*vnumnodes+i];
     2112                        Bprime[(dim*vnumnodes+pnumnodes)*6+dim*i+2] = vdbasis[2*vnumnodes+i];
     2113                        Bprime[(dim*vnumnodes+pnumnodes)*7+dim*i+0] = 0.;
     2114                        Bprime[(dim*vnumnodes+pnumnodes)*7+dim*i+1] = 0.;
     2115                        Bprime[(dim*vnumnodes+pnumnodes)*7+dim*i+2] = 0.;
     2116                }
     2117                for(int i=0;i<pnumnodes;i++){
     2118                        Bprime[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
     2119                        Bprime[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
     2120                        Bprime[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
     2121                        Bprime[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = 0.;
     2122                        Bprime[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = 0.;
     2123                        Bprime[(dim*vnumnodes+pnumnodes)*5+(dim*vnumnodes)+i] = 0.;
     2124                        Bprime[(dim*vnumnodes+pnumnodes)*6+(dim*vnumnodes)+i] = 0.;
     2125                        Bprime[(dim*vnumnodes+pnumnodes)*7+(dim*vnumnodes)+i] = pbasis[i];
     2126                }
     2127        }
     2128
     2129        /*Clean up*/
     2130        xDelete<IssmDouble>(vdbasis);
     2131        xDelete<IssmDouble>(pbasis);
    18182132}/*}}}*/
    18192133void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r16829 r16831  
    2222                /*Finite element Analysis*/
    2323                ElementMatrix* CreateKMatrix(Element* element);
     24                ElementVector* CreatePVector(Element* element);
     25                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
     26                void GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
     27                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
     28
     29                /*SSA*/
     30                ElementMatrix* CreateKMatrixSSA(Element* element);
     31                ElementMatrix* CreateKMatrixSSAViscous(Element* element);
     32                ElementMatrix* CreateKMatrixSSAFriction(Element* element);
     33                ElementVector* CreatePVectorSSA(Element* element);
     34                ElementVector* CreatePVectorSSADrivingStress(Element* element);
     35                ElementVector* CreatePVectorSSAFront(Element* element);
     36                void GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     37                void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     38                void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
     39                /*L1L2*/
     40                void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
     41                /*HO*/
    2442                ElementMatrix* CreateKMatrixHO(Element* element);
    2543                ElementMatrix* CreateKMatrixHOViscous(Element* element);
    2644                ElementMatrix* CreateKMatrixHOFriction(Element* element);
    27                 ElementMatrix* CreateKMatrixSSA(Element* element);
    28                 ElementMatrix* CreateKMatrixSSAViscous(Element* element);
    29                 ElementMatrix* CreateKMatrixSSAFriction(Element* element);
    30                 ElementVector* CreatePVector(Element* element);
     45                ElementVector* CreatePVectorHO(Element* element);
     46                ElementVector* CreatePVectorHODrivingStress(Element* element);
     47                ElementVector* CreatePVectorHOFront(Element* element);
     48                void GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     49                void GetBHOprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     50                void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
     51                /*FS*/
     52                ElementMatrix* CreateKMatrixFS(Element* element);
     53                ElementMatrix* CreateKMatrixFSViscous(Element* element);
     54                ElementMatrix* CreateKMatrixFSFriction(Element* element);
    3155                ElementVector* CreatePVectorFS(Element* element);
    3256                ElementVector* CreatePVectorFSViscous(Element* element);
    3357                ElementVector* CreatePVectorFSShelf(Element* element);
    3458                ElementVector* CreatePVectorFSFront(Element* element);
    35                 ElementVector* CreatePVectorHO(Element* element);
    36                 ElementVector* CreatePVectorHODrivingStress(Element* element);
    37                 ElementVector* CreatePVectorHOFront(Element* element);
    38                 ElementVector* CreatePVectorSSA(Element* element);
    39                 ElementVector* CreatePVectorSSADrivingStress(Element* element);
    40                 ElementVector* CreatePVectorSSAFront(Element* element);
    41                 void GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    42                 void GetBHOprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    43                 void GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    44                 void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    45                 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
     59                void GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     60                void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    4661                void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
    47                 void GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
    48                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    4962                void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
    50                 void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
     63                /*Coupling*/
    5164                void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
    52                 void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
    53                 void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
    5465                void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);
    5566                void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);

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

    r16829 r16831  
    5757                virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0;
    5858                virtual void   NodalFunctions(IssmDouble* basis,Gauss* gauss)=0;
     59                virtual void   NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss)=0;
     60                virtual void   NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss)=0;
    5961                virtual void   NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
    60                 virtual void   NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss)=0;
    61                 virtual void   NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss)=0;
     62                virtual void   NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
    6263                virtual void   NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0;
    6364                virtual void   NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0;
    … …  
    144145                virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0;
    145146                virtual void   ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
     147                virtual void   ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
    146148                virtual void   ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
    147149                virtual void   ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;

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

    r16829 r16831  
    25682568}
    25692569/*}}}*/
     2570/*FUNCTION Penta::NodalFunctionsDerivativesVelocity{{{*/
     2571void Penta::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){
     2572
     2573        _assert_(gauss->Enum()==GaussPentaEnum);
     2574        this->GetNodalFunctionsDerivativesVelocity(dbasis,xyz_list,(GaussPenta*)gauss);
     2575
     2576}
     2577/*}}}*/
    25702578/*FUNCTION Penta::NodalFunctionsVelocity{{{*/
    25712579void Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){
    … …  
    32153223
    32163224        ::TransformStiffnessMatrixCoord(Ke,this->nodes,this->NumberofNodes(),transformenum);
     3225
     3226}
     3227/*}}}*/
     3228/*FUNCTION Penta::TransformStiffnessMatrixCoord(ElementMatrix* pe,int* transformenum_list){{{*/
     3229void Penta::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){
     3230
     3231        ::TransformStiffnessMatrixCoord(Ke,this->nodes,this->NumberofNodes(),transformenum_list);
    32173232
    32183233}
    … …  
    34823497        /*Assign output pointer*/
    34833498        *pphi = phi;
     3499}
     3500/*}}}*/
     3501/*FUNCTION Penta::ViscosityFS{{{*/
     3502void Penta::ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){
     3503
     3504        /*Intermediaries*/
     3505        IssmDouble viscosity;
     3506        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     3507
     3508        this->GetStrainRate3d(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     3509        material->GetViscosity3dFS(&viscosity, &epsilon[0]);
     3510
     3511        /*Assign output pointer*/
     3512        *pviscosity=viscosity;
    34843513}
    34853514/*}}}*/

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

    r16829 r16831  
    256256                void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
    257257                void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
     258                void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    258259                void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
    259260                void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
    … …  
    272273                void           TransformSolutionCoord(IssmDouble* values,int numnodes,int* transformenum_list);/*Tiling only*/
    273274                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum);
    274                 void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){_error_("not implemented yet");};
     275                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list);
    275276                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum){_error_("not implemented yet");};
    276277                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};
     278                void           ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    277279                void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    278280                void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};

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

    r16829 r16831  
    120120                void        NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
    121121                void        NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
     122                void        NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    122123                bool        NoIceInElement(){_error_("not implemented yet");};
    123124                void        NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    … …  
    145146                ElementMatrix* NewElementMatrix(int approximation_enum){_error_("not implemented yet");};
    146147                void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
     148                void           ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");};
    147149                void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    148150                void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};

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

    r16827 r16831  
    22432243}
    22442244/*}}}*/
     2245/*FUNCTION Tria::NodalFunctionsDerivativesVelocity{{{*/
     2246void Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){
     2247
     2248        _assert_(gauss->Enum()==GaussTriaEnum);
     2249        this->GetNodalFunctionsDerivativesVelocity(dbasis,xyz_list,(GaussTria*)gauss);
     2250
     2251}
     2252/*}}}*/
    22452253/*FUNCTION Tria::NodalFunctionsVelocity{{{*/
    22462254void Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){
    … …  
    27632771
    27642772        ::TransformStiffnessMatrixCoord(Ke,this->nodes,this->NumberofNodes(),transformenum);
     2773
     2774}
     2775/*}}}*/
     2776/*FUNCTION Tria::TransformStiffnessMatrixCoord(ElementMatrix* pe,int* transformenum_list){{{*/
     2777void Tria::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){
     2778
     2779        ::TransformStiffnessMatrixCoord(Ke,this->nodes,this->NumberofNodes(),transformenum_list);
    27652780
    27662781}
    … …  
    29082923        delete gauss;
    29092924
     2925}
     2926/*}}}*/
     2927/*FUNCTION Tria::ViscosityFS{{{*/
     2928void Tria::ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){
     2929
     2930        /*Intermediaries*/
     2931        IssmDouble viscosity;
     2932        IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];    */
     2933
     2934        this->GetStrainRate2d(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     2935        material->GetViscosity2dvertical(&viscosity, &epsilon[0]);
     2936
     2937        /*Assign output pointer*/
     2938        *pviscosity=viscosity;
    29102939}
    29112940/*}}}*/

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

    r16829 r16831  
    286286                void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
    287287                void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
     288                void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
    288289                void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
    289290                void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
    … …  
    301302                void           TransformSolutionCoord(IssmDouble* values,int numnodes,int* transformenum_list){_error_("not implemented yet");};
    302303                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum);
    303                 void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){_error_("not implemented yet");};
     304                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list);
    304305                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum){_error_("not implemented yet");};
    305306                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};
    306307                void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
     308                void           ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    307309                void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented yet");};
    308310                void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
Note: See TracChangeset for help on using the changeset viewer.