Changeset 17297


Ignore:
Timestamp:
02/18/14 16:13:16 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: preparing SSA 1d

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

Legend:

Unmodified
Added
Removed

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

    r17291 r17297  
    66#include "../solutionsequences/solutionsequences.h"
    77
    8 //#define FSANALYTICAL 5
     8#define FSANALYTICAL 5
    99
    1010/*Model processing*/
    … …  
    11761176
    11771177        /*Intermediaries*/
     1178        int         dim,meshtype;
    11781179        bool        mainlyfloating;
    11791180        int         migration_style,point1;
    … …  
    11831184        Gauss*      gauss     = NULL;
    11841185
     1186        /*Get problem dimension*/
     1187        element->FindParam(&meshtype,MeshTypeEnum);
     1188        switch(meshtype){
     1189                case Mesh2DverticalEnum:   dim = 1;break;
     1190                case Mesh2DhorizontalEnum: dim = 2;break;
     1191                case Mesh3DEnum:           dim = 2;break;
     1192                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1193        }
     1194
    11851195        /*Fetch number of nodes and dof for this finite element*/
    11861196        int numnodes = element->GetNumberOfNodes();
    … …  
    11881198
    11891199        /*Initialize Element matrix and vectors*/
    1190         ElementMatrix* Ke      = element->NewElementMatrix(SSAApproximationEnum);
    1191         IssmDouble*    B       = xNew<IssmDouble>(2*numdof);
    1192         IssmDouble     D[2][2] = {0.};
     1200        ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum);
     1201        IssmDouble*    B  = xNew<IssmDouble>(dim*numdof);
     1202        IssmDouble*    D  = xNewZeroInit<IssmDouble>(dim*dim);
    11931203
    11941204        /*Retrieve all inputs and parameters*/
    … …  
    11971207        Input* surface_input    = element->GetInput(SurfaceEnum); _assert_(surface_input);
    11981208        Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
    1199         Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
    12001209        Input* gllevelset_input = NULL;
     1210        Input* vy_input         = NULL;
     1211        if(dim==2){vy_input     = element->GetInput(VyEnum);      _assert_(vy_input);}
    12011212
    12021213        /*build friction object, used later on: */
    … …  
    12251236                }
    12261237
    1227                 this->GetBSSAFriction(B, element, xyz_list, gauss);
     1238                this->GetBSSAFriction(B,element,dim,xyz_list,gauss);
    12281239                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
    1229                 for(int i=0;i<2;i++) D[i][i]=alpha2*gauss->weight*Jdet;
    1230 
    1231                 TripleMultiply(B,2,numdof,1,
    1232                                         &D[0][0],2,2,0,
    1233                                         B,2,numdof,0,
     1240                for(int i=0;i<dim;i++) D[i*dim+i]=alpha2*gauss->weight*Jdet;
     1241
     1242                TripleMultiply(B,dim,numdof,1,
     1243                                        D,dim,dim,0,
     1244                                        B,dim,numdof,0,
    12341245                                        &Ke->values[0],1);
    12351246        }
    … …  
    12431254        xDelete<IssmDouble>(xyz_list);
    12441255        xDelete<IssmDouble>(B);
     1256        xDelete<IssmDouble>(D);
    12451257        return Ke;
    12461258}/*}}}*/
    … …  
    12511263
    12521264        /*Intermediaries*/
     1265        int         dim,meshtype,bsize;
    12531266        IssmDouble  viscosity,newviscosity,oldviscosity;
    12541267        IssmDouble  viscosity_overshoot,thickness,Jdet;
    … …  
    12561269        IssmDouble *xyz_list = NULL;
    12571270
     1271        /*Get problem dimension*/
     1272        element->FindParam(&meshtype,MeshTypeEnum);
     1273        switch(meshtype){
     1274                case Mesh2DverticalEnum:   dim = 1; bsize = 1; break;
     1275                case Mesh2DhorizontalEnum: dim = 2; bsize = 1; break;
     1276                case Mesh3DEnum:           dim = 2; bsize = 3; break;
     1277                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1278        }
     1279
    12581280        /*Fetch number of nodes and dof for this finite element*/
    12591281        int numnodes = element->GetNumberOfNodes();
    1260         int numdof   = numnodes*2;
     1282        int numdof   = numnodes*dim;
    12611283
    12621284        /*Initialize Element matrix and vectors*/
    12631285        ElementMatrix* Ke     = element->NewElementMatrix(SSAApproximationEnum);
    1264         IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
    1265         IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
    1266         IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
     1286        IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
     1287        IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
     1288        IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
    12671289
    12681290        /*Retrieve all inputs and parameters*/
    … …  
    12701292        Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
    12711293        Input* vx_input=element->GetInput(VxEnum);               _assert_(vx_input);
    1272         Input* vy_input=element->GetInput(VyEnum);               _assert_(vy_input);
    12731294        Input* vxold_input=element->GetInput(VxPicardEnum);      _assert_(vxold_input);
    1274         Input* vyold_input=element->GetInput(VyPicardEnum);      _assert_(vyold_input);
     1295        Input* vy_input    = NULL;
     1296        Input* vyold_input = NULL;
     1297        if(dim==3){
     1298                vy_input=element->GetInput(VyEnum);          _assert_(vy_input);
     1299                vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
     1300        }
    12751301        element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
    12761302
    … …  
    12811307
    12821308                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    1283                 this->GetBSSA(B,element,xyz_list,gauss);
    1284                 this->GetBSSAprime(Bprime,element,xyz_list,gauss);
    1285 
    1286                 element->ViscositySSA(&viscosity,xyz_list,gauss,vx_input,vy_input);
    1287                 element->ViscositySSA(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
     1309                this->GetBSSA(B,element,dim,xyz_list,gauss);
     1310                this->GetBSSAprime(Bprime,element,dim,xyz_list,gauss);
     1311
     1312                element->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
     1313                element->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
    12881314                thickness_input->GetInputValue(&thickness, gauss);
    12891315
    … …  
    12921318                for(int i=0;i<3;i++) D[i*3+i]=D_scalar;
    12931319
    1294                 TripleMultiply(B,3,numdof,1,
    1295                                         D,3,3,0,
    1296                                         Bprime,3,numdof,0,
     1320                TripleMultiply(B,bsize,numdof,1,
     1321                                        D,bsize,bsize,0,
     1322                                        Bprime,bsize,numdof,0,
    12971323                                        &Ke->values[0],1);
    12981324        }
    12991325
    13001326        /*Transform Coordinate System*/
    1301         element->TransformStiffnessMatrixCoord(Ke,XYEnum);
     1327        if(dim==2) element->TransformStiffnessMatrixCoord(Ke,XYEnum);
    13021328
    13031329        /*Clean up and return*/
    … …  
    13481374
    13491375        /*Intermediaries */
     1376        int         dim,meshtype;
    13501377        IssmDouble  thickness,Jdet,slope[2];
    13511378        IssmDouble* xyz_list = NULL;
     1379
     1380        /*Get problem dimension*/
     1381        element->FindParam(&meshtype,MeshTypeEnum);
     1382        switch(meshtype){
     1383                case Mesh2DverticalEnum:   dim = 1;break;
     1384                case Mesh2DhorizontalEnum: dim = 2;break;
     1385                case Mesh3DEnum:           dim = 2;break;
     1386                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1387        }
    13521388
    13531389        /*Fetch number of nodes and dof for this finite element*/
    … …  
    13771413
    13781414                for(int i=0;i<numnodes;i++){
    1379                         pe->values[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i];
    1380                         pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i];
     1415                        pe->values[i*dim+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i];
     1416                        if(dim==2) pe->values[i*dim+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i];
    13811417                }
    13821418        }
    … …  
    14001436
    14011437        /*Intermediaries*/
     1438        int         dim,meshtype;
    14021439        IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure;
    14031440        IssmDouble  surface_under_water,base_under_water,pressure;
    … …  
    14051442        IssmDouble *xyz_list_front = NULL;
    14061443        IssmDouble  normal[2];
     1444
     1445        /*Get problem dimension*/
     1446        element->FindParam(&meshtype,MeshTypeEnum);
     1447        switch(meshtype){
     1448                case Mesh2DverticalEnum:   dim = 1;break;
     1449                case Mesh2DhorizontalEnum: dim = 2;break;
     1450                case Mesh3DEnum:           dim = 2;break;
     1451                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1452        }
    14071453
    14081454        /*Fetch number of nodes for this finite element*/
    … …  
    14201466        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    14211467        element->GetVerticesCoordinates(&xyz_list);
    1422 //      element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
     1468        //element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
    14231469        element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
    14241470        element->NormalSection(&normal[0],xyz_list_front);
    … …  
    14401486                pressure = ice_pressure + water_pressure;
    14411487
    1442                 for (int i=0;i<numnodes;i++){
    1443                         pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
    1444                         pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
     1488                for(int i=0;i<numnodes;i++){
     1489                        pe->values[dim*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
     1490                        if(dim==2) pe->values[dim*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
    14451491                }
    14461492        }
    … …  
    14561502        return pe;
    14571503}/*}}}*/
    1458 void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1504void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    14591505        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    14601506         * For node i, Bi can be expressed in the actual coordinate system
    14611507         * by:
    1462          *       Bi=[ dN/dx           0    ]
    1463          *          [   0           dN/dy  ]
    1464          *          [ 1/2*dN/dy  1/2*dN/dx ]
     1508         *                   2D                      1D
     1509         *       Bi=[ dN/dx           0    ]   Bi=[ dN/dx ]
     1510         *          [   0           dN/dy  ]     
     1511         *          [ 1/2*dN/dy  1/2*dN/dx ]     
    14651512         * where N is the finiteelement function for node i.
    14661513         *
    … …  
    14721519
    14731520        /*Get nodal functions derivatives*/
    1474         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     1521        IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
    14751522        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    14761523
    14771524        /*Build B: */
    1478         for(int i=0;i<numnodes;i++){
    1479                 B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
    1480                 B[2*numnodes*0+2*i+1] = 0.;
    1481                 B[2*numnodes*1+2*i+0] = 0.;
    1482                 B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
    1483                 B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
    1484                 B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
     1525        if(dim==2){
     1526                for(int i=0;i<numnodes;i++){
     1527                        B[i] = dbasis[i];
     1528                }
     1529        }
     1530        else{
     1531                for(int i=0;i<numnodes;i++){
     1532                        B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
     1533                        B[2*numnodes*0+2*i+1] = 0.;
     1534                        B[2*numnodes*1+2*i+0] = 0.;
     1535                        B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
     1536                        B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
     1537                        B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
     1538                }
    14851539        }
    14861540
    … …  
    14881542        xDelete<IssmDouble>(dbasis);
    14891543}/*}}}*/
    1490 void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1544void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    14911545        /*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2.
    14921546         * For node i, Bi can be expressed in the actual coordinate system
    14931547         * by:
    1494          *                 Bi=[ N   0 ]
     1548         *                       2D             1D
     1549         *                 Bi=[ N   0 ]    Bi = N
    14951550         *                    [ 0   N ]
    14961551         * where N is the finiteelement function for node i.
    … …  
    15071562
    15081563        /*Build L: */
    1509         for(int i=0;i<numnodes;i++){
    1510                 B[2*numnodes*0+2*i+0] = basis[i];
    1511                 B[2*numnodes*0+2*i+1] = 0.;
    1512                 B[2*numnodes*1+2*i+0] = 0.;
    1513                 B[2*numnodes*1+2*i+1] = basis[i];
     1564        if(dim==2){
     1565                for(int i=0;i<numnodes;i++){
     1566                        B[2*numnodes*0+2*i+0] = basis[i];
     1567                        B[2*numnodes*0+2*i+1] = 0.;
     1568                        B[2*numnodes*1+2*i+0] = 0.;
     1569                        B[2*numnodes*1+2*i+1] = basis[i];
     1570                }
     1571        }
     1572        else{
     1573                for(int i=0;i<numnodes;i++){
     1574                        B[i] = basis[i];
     1575                }
    15141576        }
    15151577
    … …  
    15171579        xDelete<IssmDouble>(basis);
    15181580}/*}}}*/
    1519 void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     1581void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    15201582        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    15211583         * For node i, Bi' can be expressed in the actual coordinate system
    15221584         * by:
    1523          *       Bi_prime=[ 2*dN/dx    dN/dy ]
     1585         *                         2D                        1D
     1586         *       Bi_prime=[ 2*dN/dx    dN/dy ]     Bi_prime=[ 2*dN/dx ]
    15241587         *                [   dN/dx  2*dN/dy ]
    15251588         *                [   dN/dy    dN/dx ]
    … …  
    15331596
    15341597        /*Get nodal functions derivatives*/
    1535         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     1598        IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
    15361599        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    15371600
    15381601        /*Build B': */
    1539         for(int i=0;i<numnodes;i++){
    1540                 Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
    1541                 Bprime[2*numnodes*0+2*i+1] =    dbasis[1*numnodes+i];
    1542                 Bprime[2*numnodes*1+2*i+0] =    dbasis[0*numnodes+i];
    1543                 Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
    1544                 Bprime[2*numnodes*2+2*i+0] =    dbasis[1*numnodes+i];
    1545                 Bprime[2*numnodes*2+2*i+1] =    dbasis[0*numnodes+i];
     1602        if(dim==2){
     1603                for(int i=0;i<numnodes;i++){
     1604                        Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
     1605                        Bprime[2*numnodes*0+2*i+1] =    dbasis[1*numnodes+i];
     1606                        Bprime[2*numnodes*1+2*i+0] =    dbasis[0*numnodes+i];
     1607                        Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
     1608                        Bprime[2*numnodes*2+2*i+0] =    dbasis[1*numnodes+i];
     1609                        Bprime[2*numnodes*2+2*i+1] =    dbasis[0*numnodes+i];
     1610                }
     1611        }
     1612        else{
     1613                for(int i=0;i<numnodes;i++){
     1614                        Bprime[i] = 2.*dbasis[i];
     1615                }
    15461616        }
    15471617
    … …  
    15511621void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
    15521622
    1553         int         i,meshtype;
     1623        int         i,dim,meshtype;
    15541624        IssmDouble  rho_ice,g;
    15551625        int*        doflist=NULL;
    … …  
    15701640                        element->GetInputListOnVertices(thickness,ThicknessEnum);
    15711641                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
     1642                        dim=2;
    15721643                        break;
    15731644                case Mesh3DEnum:
    … …  
    15751646                        element->GetInputListOnVertices(surface,SurfaceEnum);
    15761647                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
     1648                        dim=2;
     1649                        break;
     1650                case Mesh2DverticalEnum:
     1651                        element->GetVerticesCoordinates(&xyz_list);
     1652                        element->GetInputListOnVertices(surface,SurfaceEnum);
     1653                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
     1654                        dim=1;
    15771655                        break;
    15781656                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    … …  
    15881666                        basalelement = element;
    15891667                        break;
    1590                 case Mesh3DEnum:
     1668                case Mesh3DEnum: case Mesh2DverticalEnum:
    15911669                        if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
    15921670                        basalelement=element->SpawnBasalElement();
    … …  
    15971675        /*Fetch number of nodes and dof for this finite element*/
    15981676        int numnodes = basalelement->GetNumberOfNodes();
    1599         int numdof   = numnodes*2;
     1677        int numdof   = numnodes*dim;
    16001678
    16011679        /*Fetch dof list and allocate solution vectors*/
    … …  
    16111689
    16121690        /*Transform solution in Cartesian Space*/
    1613         basalelement->TransformSolutionCoord(&values[0],XYEnum);
    1614         basalelement->FindParam(&meshtype,MeshTypeEnum);
     1691        if(dim==2) basalelement->TransformSolutionCoord(&values[0],XYEnum);
    16151692
    16161693        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    16171694        for(i=0;i<numnodes;i++){
    16181695                vx[i]=values[i*2+0];
    1619                 vy[i]=values[i*2+1];
    1620 
    1621                 /*Check solution*/
    16221696                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    1623                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     1697
     1698                if(dim==2){
     1699                        vy[i]=values[i*2+1];
     1700                        if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     1701                }
    16241702        }
    16251703
    16261704        /*Get Vz and compute vel*/
    1627         basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
    1628         for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     1705        if(dim==2){
     1706                basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
     1707                for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     1708        }
     1709        else{
     1710                element->GetInputListOnNodes(&vy[0],VyEnum,0.);
     1711                for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
     1712        }
    16291713
    16301714        /*Now, we have to move the previous Vx and Vy inputs  to old
    16311715         * status, otherwise, we'll wipe them off: */
    16321716        element->InputChangeName(VxEnum,VxPicardEnum);
    1633         element->InputChangeName(VyEnum,VyPicardEnum);
     1717        if(dim==2)element->InputChangeName(VyEnum,VyPicardEnum);
    16341718
    16351719        /*Add vx and vy as inputs to the tria element: */
    16361720        element->AddBasalInput(VxEnum,vx,P1Enum);
    1637         element->AddBasalInput(VyEnum,vy,P1Enum);
     1721        if(dim==2)element->AddBasalInput(VyEnum,vy,P1Enum);
    16381722        element->AddBasalInput(VelEnum,vel,P1Enum);
    16391723
    … …  
    17061790
    17071791                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    1708                 this->GetBSSA(B,basalelement,xyz_list,gauss_base);
    1709                 this->GetBSSAprime(Bprime,basalelement,xyz_list,gauss_base);
     1792                this->GetBSSA(B,basalelement,2,xyz_list,gauss_base);
     1793                this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base);
    17101794
    17111795                element->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
    … …  
    24282512                        Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
    24292513                }
    2430                 }
     2514        }
    24312515
    24322516        /*Clean-up*/
    … …  
    36963780
    36973781                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
    3698                 this->GetBSSA(&B[0][0],basaltria, xyz_list, gauss_tria);
    3699                 this->GetBSSAprime(&Bprime[0][0], basaltria,xyz_list, gauss_tria);
     3782                this->GetBSSA(&B[0][0],basaltria,2,xyz_list, gauss_tria);
     3783                this->GetBSSAprime(&Bprime[0][0],basaltria,2,xyz_list, gauss_tria);
    37003784
    37013785                if(approximation==SSAHOApproximationEnum){
    … …  
    39984082                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
    39994083                this->GetBSSAHO(B, element,xyz_list, gauss);
    4000                 this->GetBSSAprime(Bprime, basaltria,xyz_list, gauss_tria);
     4084                this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria);
    40014085                element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input);
    40024086                element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input);

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

    r17242 r17297  
    3939                ElementVector* CreatePVectorSSADrivingStress(Element* element);
    4040                ElementVector* CreatePVectorSSAFront(Element* element);
    41                 void GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    42                 void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    43                 void GetBSSAFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     41                void GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     42                void GetBSSAprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     43                void GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    4444                void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
    4545                /*L1L2*/
Note: See TracChangeset for help on using the changeset viewer.