Changeset 17242


Ignore:
Timestamp:
02/07/14 16:57:05 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: preparing HO 2dvertical

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

Legend:

Unmodified
Added
Removed

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

    r17241 r17242  
    1818                case SSAApproximationEnum:  numdofs =2; break;
    1919                case L1L2ApproximationEnum: numdofs =2; break;
    20                 case HOApproximationEnum:   numdofs =2; break;
     20                case HOApproximationEnum:   
     21                         switch(meshtype){
     22                                 case Mesh3DEnum:         numdofs=2; break;
     23                                 case Mesh2DverticalEnum: numdofs=1; break;
     24                                 default: _error_("mesh type not supported yet");
     25                         }
     26                         break;
    2127                case SIAApproximationEnum:  numdofs =2; break;
    2228                case FSvelocityEnum:
    … …  
    478484                else{
    479485                        IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1);
    480                         IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2);
     486                        if(iomodel->meshtype==Mesh3DEnum){
     487                                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceAnalysisEnum,finiteelement,2);
     488                        }
    481489                }
    482490
    … …  
    955963
    956964        IssmDouble   vx,vy;
    957         int          approximation;
     965        int          meshtype,dim,approximation;
    958966        int*         doflist = NULL;
     967
     968        /*Get some parameters*/
     969        element->FindParam(&meshtype,MeshTypeEnum);
     970        switch(meshtype){
     971                case Mesh2DverticalEnum: dim = 2; break;
     972                case Mesh3DEnum:         dim = 3; break;
     973                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     974        }
    959975
    960976        /*Fetch number of nodes and dof for this finite element*/
    961977        int numnodes = element->GetNumberOfNodes();
    962         int numdof   = numnodes*2;
     978        int numdof   = numnodes*(dim-1);
    963979        element->GetInputValue(&approximation,ApproximationEnum);
    964980
    … …  
    969985        /*Get inputs*/
    970986        Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
    971         Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
     987        Input* vy_input=NULL;
     988        if(dim==3){vy_input=element->GetInput(VyEnum); _assert_(vy_input);}
    972989
    973990        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    978995                /*Recover vx and vy*/
    979996                vx_input->GetInputValue(&vx,gauss);
    980                 vy_input->GetInputValue(&vy,gauss);
    981                 values[i*2+0]=vx;
    982                 values[i*2+1]=vy;
     997                values[i*(dim-1)+0]=vx;
     998                if(dim==3){
     999                        vy_input->GetInputValue(&vy,gauss);
     1000                        values[i*(dim-1)+1]=vy;
     1001                }
    9831002        }
    9841003
    … …  
    20212040
    20222041        /*Intermediaries*/
     2042        int         dim,meshtype,bsize;
    20232043        IssmDouble  viscosity,newviscosity,oldviscosity;
    20242044        IssmDouble  viscosity_overshoot,thickness,Jdet;
    … …  
    20262046        IssmDouble *xyz_list = NULL;
    20272047
     2048        /*Get problem dimension*/
     2049        element->FindParam(&meshtype,MeshTypeEnum);
     2050        switch(meshtype){
     2051                case Mesh2DverticalEnum: dim = 2; bsize = 2; break;
     2052                case Mesh3DEnum:         dim = 3; bsize = 5; break;
     2053                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     2054        }
     2055
    20282056        /*Fetch number of nodes and dof for this finite element*/
    20292057        int numnodes = element->GetNumberOfNodes();
    … …  
    20322060        /*Initialize Element matrix and vectors*/
    20332061        ElementMatrix* Ke     = element->NewElementMatrix(HOApproximationEnum);
    2034         IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
    2035         IssmDouble*    Bprime = xNew<IssmDouble>(5*numdof);
    2036         IssmDouble*    D      = xNewZeroInit<IssmDouble>(5*5);
     2062        IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
     2063        IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
     2064        IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
    20372065
    20382066        /*Retrieve all inputs and parameters*/
    20392067        element->GetVerticesCoordinates(&xyz_list);
    2040         Input* vx_input=element->GetInput(VxEnum);               _assert_(vx_input);
    2041         Input* vy_input=element->GetInput(VyEnum);               _assert_(vy_input);
    2042         Input* vxold_input=element->GetInput(VxPicardEnum);      _assert_(vxold_input);
    2043         Input* vyold_input=element->GetInput(VyPicardEnum);      _assert_(vyold_input);
     2068        Input* vx_input    = element->GetInput(VxEnum);       _assert_(vx_input);
     2069        Input* vxold_input = element->GetInput(VxPicardEnum); _assert_(vxold_input);
     2070        Input* vy_input    = NULL;
     2071        Input* vyold_input = NULL;
     2072        if(dim==3){
     2073                vy_input=element->GetInput(VyEnum);          _assert_(vy_input);
     2074                vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
     2075        }
    20442076        element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
    20452077
    … …  
    20502082
    20512083                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    2052                 this->GetBHO(B,element,xyz_list,gauss);
    2053                 this->GetBHOprime(Bprime,element,xyz_list,gauss);
    2054 
    2055                 element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input);
    2056                 element->ViscosityHO(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
     2084                this->GetBHO(B,element,dim,xyz_list,gauss);
     2085                this->GetBHOprime(Bprime,element,dim,xyz_list,gauss);
     2086
     2087                element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
     2088                element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
    20572089
    20582090                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    20592091                D_scalar=2.*newviscosity*gauss->weight*Jdet;
    2060                 for(int i=0;i<5;i++) D[i*5+i]=D_scalar;
    2061 
    2062                 TripleMultiply(B,5,numdof,1,
    2063                                         D,5,5,0,
    2064                                         Bprime,5,numdof,0,
     2092                for(int i=0;i<bsize;i++) D[i*bsize+i]=D_scalar;
     2093
     2094                TripleMultiply(B,bsize,numdof,1,
     2095                                        D,bsize,bsize,0,
     2096                                        Bprime,bsize,numdof,0,
    20652097                                        &Ke->values[0],1);
    20662098        }
    … …  
    20822114
    20832115        /*Intermediaries*/
     2116        int         dim,meshtype;
    20842117        bool        mainlyfloating;
    20852118        int         migration_style,point1;
    … …  
    20892122        Gauss*      gauss         = NULL;
    20902123
     2124        /*Get problem dimension*/
     2125        element->FindParam(&meshtype,MeshTypeEnum);
     2126        switch(meshtype){
     2127                case Mesh2DverticalEnum: dim = 2; break;
     2128                case Mesh3DEnum:         dim = 3; break;
     2129                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     2130        }
     2131
    20912132        /*Fetch number of nodes and dof for this finite element*/
    20922133        int numnodes = element->GetNumberOfNodes();
    2093         int numdof   = numnodes*2;
     2134        int numdof   = numnodes*(dim-1);
    20942135
    20952136        /*Initialize Element matrix and vectors*/
    2096         ElementMatrix* Ke      = element->NewElementMatrix(HOApproximationEnum);
    2097         IssmDouble*    B       = xNew<IssmDouble>(2*numdof);
    2098         IssmDouble     D[2][2] = {0.};
     2137        ElementMatrix* Ke = element->NewElementMatrix(HOApproximationEnum);
     2138        IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
     2139        IssmDouble*    D  = xNew<IssmDouble>((dim-1)*(dim-1));
    20992140
    21002141        /*Retrieve all inputs and parameters*/
    21012142        element->GetVerticesCoordinatesBase(&xyz_list_base);
    21022143        element->FindParam(&migration_style,GroundinglineMigrationEnum);
    2103         Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
    2104         Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
    2105         Input* vz_input         = element->GetInput(VzEnum);      _assert_(vz_input);
     2144        Input* vx_input = element->GetInput(VxEnum);   _assert_(vx_input);
     2145        Input* vy_input = element->GetInput(VyEnum);   _assert_(vy_input);
     2146        Input* vz_input = NULL;
    21062147        Input* gllevelset_input = NULL;
     2148        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
    21072149
    21082150        /*build friction object, used later on: */
    … …  
    21322174                }
    21332175
    2134                 this->GetBHOFriction(B,element,xyz_list_base,gauss);
     2176                this->GetBHOFriction(B,element,dim,xyz_list_base,gauss);
    21352177                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    2136                 for(int i=0;i<2;i++) D[i][i]=alpha2*gauss->weight*Jdet;
    2137 
    2138                 TripleMultiply(B,2,numdof,1,
    2139                                         &D[0][0],2,2,0,
    2140                                         B,2,numdof,0,
     2178                for(int i=0;i<dim-1;i++) D[i*(dim-1)+i]=alpha2*gauss->weight*Jdet;
     2179
     2180                TripleMultiply(B,dim-1,numdof,1,
     2181                                        D,dim-1,dim-1,0,
     2182                                        B,dim-1,numdof,0,
    21412183                                        &Ke->values[0],1);
    21422184        }
    … …  
    21502192        xDelete<IssmDouble>(xyz_list_base);
    21512193        xDelete<IssmDouble>(B);
     2194        xDelete<IssmDouble>(D);
    21522195        return Ke;
    21532196}/*}}}*/
    … …  
    21672210
    21682211        /*Intermediaries */
     2212        int         dim,meshtype;
    21692213        IssmDouble  Jdet,slope[3];
    21702214        IssmDouble* xyz_list = NULL;
     2215
     2216        /*Get problem dimension*/
     2217        element->FindParam(&meshtype,MeshTypeEnum);
     2218        switch(meshtype){
     2219                case Mesh2DverticalEnum: dim = 2; break;
     2220                case Mesh3DEnum:         dim = 3; break;
     2221                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     2222        }
    21712223
    21722224        /*Fetch number of nodes and dof for this finite element*/
    … …  
    21922244
    21932245                for(int i=0;i<numnodes;i++){
    2194                         pe->values[i*2+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i];
    2195                         pe->values[i*2+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i];
     2246                        pe->values[i*(dim-1)+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i];
     2247                        if(dim==3) pe->values[i*(dim-1)+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i];
    21962248                }
    21972249        }
    … …  
    22122264
    22132265        /*Intermediaries*/
     2266        int         dim,meshtype;
    22142267        IssmDouble  Jdet,surface,z,water_pressure,ice_pressure;
    22152268        IssmDouble  surface_under_water,base_under_water,pressure;
    … …  
    22182271        IssmDouble  normal[3];
    22192272        Gauss*      gauss = NULL;
     2273
     2274        /*Get problem dimension*/
     2275        element->FindParam(&meshtype,MeshTypeEnum);
     2276        switch(meshtype){
     2277                case Mesh2DverticalEnum: dim = 2; break;
     2278                case Mesh3DEnum:         dim = 3; break;
     2279                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     2280        }
    22202281
    22212282        /*Fetch number of nodes and dof for this finite element*/
    … …  
    22372298
    22382299        /*Initialize gauss points*/
    2239         IssmDouble zmax=xyz_list[0*3+2]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+2]>zmax) zmax=xyz_list[i*3+2];
    2240         IssmDouble zmin=xyz_list[0*3+2]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+2]<zmin) zmin=xyz_list[i*3+2];
    2241         if(zmax>0 && zmin<0) gauss=element->NewGauss(xyz_list,xyz_list_front,3,10);//refined in vertical because of the sea level discontinuity
    2242         else                 gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
     2300        IssmDouble zmax=xyz_list[0*3+(dim-1)]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+(dim-1)]>zmax) zmax=xyz_list[i*3+(dim-1)];
     2301        IssmDouble zmin=xyz_list[0*3+(dim-1)]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+(dim-1)]<zmin) zmin=xyz_list[i*3+(dim-1)];
     2302        if(zmax>0. && zmin<0.) gauss=element->NewGauss(xyz_list,xyz_list_front,3,10);//refined in vertical because of the sea level discontinuity
     2303        else                   gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
    22432304
    22442305        /* Start  looping on the number of gaussian points: */
    … …  
    22562317
    22572318                for (int i=0;i<numnodes;i++){
    2258                         pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
    2259                         pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
     2319                        pe->values[(dim-1)*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
     2320                        if(dim==3) pe->values[(dim-1)*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
    22602321                }
    22612322        }
    … …  
    22712332        return pe;
    22722333}/*}}}*/
    2273 void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2334void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    22742335        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    22752336         * For node i, Bi can be expressed in the actual coordinate system
    22762337         * by:
    2277          *       Bi=[ dh/dx          0      ]
    2278          *          [   0           dh/dy   ]
    2279          *          [ 1/2*dh/dy  1/2*dh/dx  ]
    2280          *          [ 1/2*dh/dz      0      ]
    2281          *          [  0         1/2*dh/dz  ]
     2338         *                   3D                        2D
     2339         *
     2340         *       Bi=[ dh/dx          0      ]  Bi=[ dh/dx]
     2341         *          [   0           dh/dy   ]     [ dh/dy]
     2342         *          [ 1/2*dh/dy  1/2*dh/dx  ]     
     2343         *          [ 1/2*dh/dz      0      ]   
     2344         *          [  0         1/2*dh/dz  ]   
    22822345         * where h is the interpolation function for node i.
    22832346         *
    … …  
    22892352
    22902353        /*Get nodal functions derivatives*/
    2291         IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
     2354        IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
    22922355        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    22932356
    22942357        /*Build B: */
    2295         for(int i=0;i<numnodes;i++){
    2296                 B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
    2297                 B[2*numnodes*0+2*i+1] = 0.;
    2298                 B[2*numnodes*1+2*i+0] = 0.;
    2299                 B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
    2300                 B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
    2301                 B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
    2302                 B[2*numnodes*3+2*i+0] = .5*dbasis[2*numnodes+i];
    2303                 B[2*numnodes*3+2*i+1] = 0.;
    2304                 B[2*numnodes*4+2*i+0] = 0.;
    2305                 B[2*numnodes*4+2*i+1] = .5*dbasis[2*numnodes+i];
    2306         }
     2358        if(dim==2){
     2359                for(int i=0;i<numnodes;i++){
     2360                        B[numnodes*0+i] = dbasis[0*numnodes+i];
     2361                        B[numnodes*1+i] = .5*dbasis[1*numnodes+i];
     2362                }
     2363        }
     2364        else{
     2365                for(int i=0;i<numnodes;i++){
     2366                        B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
     2367                        B[2*numnodes*0+2*i+1] = 0.;
     2368                        B[2*numnodes*1+2*i+0] = 0.;
     2369                        B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
     2370                        B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
     2371                        B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
     2372                        B[2*numnodes*3+2*i+0] = .5*dbasis[2*numnodes+i];
     2373                        B[2*numnodes*3+2*i+1] = 0.;
     2374                        B[2*numnodes*4+2*i+0] = 0.;
     2375                        B[2*numnodes*4+2*i+1] = .5*dbasis[2*numnodes+i];
     2376                }
     2377        }
     2378
    23072379
    23082380        /*Clean-up*/
    23092381        xDelete<IssmDouble>(dbasis);
    23102382}/*}}}*/
    2311 void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2383void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    23122384        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    23132385         * For node i, Bi' can be expressed in the actual coordinate system
    23142386         * by:
    2315          *       Bi_prime=[ 2*dN/dx    dN/dy ]
    2316          *                [   dN/dx  2*dN/dy ]
    2317          *                [   dN/dy    dN/dx ]
     2387         *                          3D                      2D
     2388         *       Bi_prime=[ 2*dN/dx    dN/dy ] Bi_prime=[ 2*dN/dx ]
     2389         *                [   dN/dx  2*dN/dy ]          [   dN/dy ]
     2390         *                [   dN/dy    dN/dx ] 
    23182391         * where hNis the finiteelement function for node i.
    23192392         *
    … …  
    23252398
    23262399        /*Get nodal functions derivatives*/
    2327         IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
     2400        IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
    23282401        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    23292402
    23302403        /*Build B': */
    2331         for(int i=0;i<numnodes;i++){
    2332                 Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
    2333                 Bprime[2*numnodes*0+2*i+1] = dbasis[1*numnodes+i];
    2334                 Bprime[2*numnodes*1+2*i+0] = dbasis[0*numnodes+i];
    2335                 Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
    2336                 Bprime[2*numnodes*2+2*i+0] = dbasis[1*numnodes+i];
    2337                 Bprime[2*numnodes*2+2*i+1] = dbasis[0*numnodes+i];
    2338                 Bprime[2*numnodes*3+2*i+0] = dbasis[2*numnodes+i];
    2339                 Bprime[2*numnodes*3+2*i+1] = 0.;
    2340                 Bprime[2*numnodes*4+2*i+0] = 0.;
    2341                 Bprime[2*numnodes*4+2*i+1] = dbasis[2*numnodes+i];
    2342         }
     2404        if(dim==3){
     2405                for(int i=0;i<numnodes;i++){
     2406                        Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
     2407                        Bprime[2*numnodes*0+2*i+1] = dbasis[1*numnodes+i];
     2408                        Bprime[2*numnodes*1+2*i+0] = dbasis[0*numnodes+i];
     2409                        Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
     2410                        Bprime[2*numnodes*2+2*i+0] = dbasis[1*numnodes+i];
     2411                        Bprime[2*numnodes*2+2*i+1] = dbasis[0*numnodes+i];
     2412                        Bprime[2*numnodes*3+2*i+0] = dbasis[2*numnodes+i];
     2413                        Bprime[2*numnodes*3+2*i+1] = 0.;
     2414                        Bprime[2*numnodes*4+2*i+0] = 0.;
     2415                        Bprime[2*numnodes*4+2*i+1] = dbasis[2*numnodes+i];
     2416                }
     2417        }
     2418        else{
     2419                for(int i=0;i<numnodes;i++){
     2420                        Bprime[numnodes*0+i] = 2.*dbasis[0*numnodes+i];
     2421                        Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
     2422                }
     2423                }
    23432424
    23442425        /*Clean-up*/
    23452426        xDelete<IssmDouble>(dbasis);
    23462427}/*}}}*/
    2347 void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     2428void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    23482429        /*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2.
    23492430         * For node i, Bi can be expressed in the actual coordinate system
    23502431         * by:
    2351          *                 Bi=[ N   0 ]
     2432         *                       3D           2D
     2433         *                 Bi=[ N   0 ]    Bi=N
    23522434         *                    [ 0   N ]
    23532435         * where N is the finiteelement function for node i.
    … …  
    23642446
    23652447        /*Build L: */
    2366         for(int i=0;i<numnodes;i++){
    2367                 B[2*numnodes*0+2*i+0] = basis[i];
    2368                 B[2*numnodes*0+2*i+1] = 0.;
    2369                 B[2*numnodes*1+2*i+0] = 0.;
    2370                 B[2*numnodes*1+2*i+1] = basis[i];
     2448        if(dim==3){
     2449                for(int i=0;i<numnodes;i++){
     2450                        B[2*numnodes*0+2*i+0] = basis[i];
     2451                        B[2*numnodes*0+2*i+1] = 0.;
     2452                        B[2*numnodes*1+2*i+0] = 0.;
     2453                        B[2*numnodes*1+2*i+1] = basis[i];
     2454                }
     2455        }
     2456        else{
     2457                for(int i=0;i<numnodes;i++){
     2458                        B[i] = basis[i];
     2459                }
    23712460        }
    23722461
    … …  
    23762465void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
    23772466
    2378         int         i;
     2467        int         i,dim,meshtype;
    23792468        int*        doflist=NULL;
    23802469        IssmDouble* xyz_list=NULL;
     2470
     2471        /*Get mesh dimension*/
     2472        element->FindParam(&meshtype,MeshTypeEnum);
     2473        switch(meshtype){
     2474                case Mesh2DverticalEnum: dim = 2; break;
     2475                case Mesh3DEnum:         dim = 3; break;
     2476                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     2477        }
    23812478
    23822479        /*Deal with pressure first*/
    … …  
    23952492        /*Fetch number of nodes and dof for this finite element*/
    23962493        int numnodes = element->GetNumberOfNodes();
    2397         int numdof   = numnodes*2;
     2494        int numdof   = numnodes*(dim-1);
    23982495
    23992496        /*Fetch dof list and allocate solution vectors*/
    24002497        element->GetDofList(&doflist,HOApproximationEnum,GsetEnum);
    2401         IssmDouble* values    = xNew<IssmDouble>(numdof);
    2402         IssmDouble* vx        = xNew<IssmDouble>(numnodes);
    2403         IssmDouble* vy        = xNew<IssmDouble>(numnodes);
    2404         IssmDouble* vz        = xNew<IssmDouble>(numnodes);
    2405         IssmDouble* vel       = xNew<IssmDouble>(numnodes);
     2498        IssmDouble* values = xNew<IssmDouble>(numdof);
     2499        IssmDouble* vx     = xNew<IssmDouble>(numnodes);
     2500        IssmDouble* vy     = xNew<IssmDouble>(numnodes);
     2501        IssmDouble* vz     = xNew<IssmDouble>(numnodes);
     2502        IssmDouble* vel    = xNew<IssmDouble>(numnodes);
    24062503
    24072504        /*Use the dof list to index into the solution vector: */
    … …  
    24132510        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    24142511        for(i=0;i<numnodes;i++){
    2415                 vx[i]=values[i*2+0];
    2416                 vy[i]=values[i*2+1];
    2417 
    2418                 /*Check solution*/
     2512                vx[i]=values[i*(dim-1)+0];
    24192513                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    2420                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     2514                if(dim==3){
     2515                        vy[i]=values[i*(dim-1)+1];
     2516                        if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     2517                }
    24212518        }
    24222519
    24232520        /*Get Vz and compute vel*/
    2424         element->GetInputListOnNodes(&vz[0],VzEnum,0.);
    2425         for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     2521        if(dim==3){
     2522                element->GetInputListOnNodes(&vz[0],VzEnum,0.);
     2523                for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     2524        }
     2525        else{
     2526                element->GetInputListOnNodes(&vy[0],VyEnum,0.);
     2527                for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
     2528        }
    24262529
    24272530        /*Now, we have to move the previous Vx and Vy inputs  to old
    24282531         * status, otherwise, we'll wipe them off: */
    24292532        element->InputChangeName(VxEnum,VxPicardEnum);
    2430         element->InputChangeName(VyEnum,VyPicardEnum);
    2431         element->InputChangeName(PressureEnum,PressurePicardEnum);
     2533        if(dim==3)element->InputChangeName(VyEnum,VyPicardEnum);
    24322534
    24332535        /*Add vx and vy as inputs to the element: */
    … …  
    24362538        //element->AddInput(VelEnum,vel,element->GetElementType());
    24372539        element->AddInput(VxEnum,vx,P1Enum);
    2438         element->AddInput(VyEnum,vy,P1Enum);
     2540        if(dim==3)element->AddInput(VyEnum,vy,P1Enum);
    24392541        element->AddInput(VelEnum,vel,P1Enum);
    24402542
    … …  
    35793681
    35803682                if(approximation==SSAHOApproximationEnum){
    3581                         element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input);
    3582                         element->ViscosityHO(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
     3683                        element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
     3684                        element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
    35833685                        newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    35843686                }
    … …  
    37933895                friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    37943896                element->JacobianDeterminantBase(&Jdet2d, xyz_list_tria,gauss);
    3795                 this->GetBHOFriction(L,element,xyz_list_tria,gauss);
     3897                this->GetBHOFriction(L,element,3,xyz_list_tria,gauss);
    37963898
    37973899                DL_scalar=alpha2*gauss->weight*Jdet2d;
    … …  
    38783980                this->GetBSSAHO(B, element,xyz_list, gauss);
    38793981                this->GetBSSAprime(Bprime, basaltria,xyz_list, gauss_tria);
    3880                 element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input);
    3881                 element->ViscosityHO(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
     3982                element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input);
     3983                element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input);
    38823984
    38833985                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);

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

    r17212 r17242  
    5959                ElementVector* CreatePVectorHODrivingStress(Element* element);
    6060                ElementVector* CreatePVectorHOFront(Element* element);
    61                 void GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    62                 void GetBHOprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    63                 void GetBHOFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     61                void GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     62                void GetBHOprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     63                void GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    6464                void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
    6565                /*FS*/

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

    r17100 r17242  
    378378        *pviscosity = viscosity;
    379379}/*}}}*/
    380 void Element::ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     380void Element::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
    381381
    382382        /*Intermediaries*/
    383383        IssmDouble viscosity;
    384         IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
    385 
    386         this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    387         material->GetViscosity3d(&viscosity, &epsilon[0]);
     384        IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
     385        IssmDouble epsilon2d[2]; /* epsilon=[exx,exy];           */
     386
     387        if(dim==3){
     388                this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
     389                material->GetViscosity3d(&viscosity, &epsilon3d[0]);
     390        }
     391        else{
     392                this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
     393                material->GetViscosity2dverticalHO(&viscosity, &epsilon2d[0]);
     394        }
    388395
    389396        /*Assign output pointer*/
    … …  
    469476
    470477}/*}}}*/
     478void Element::StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     479        /*Compute the 2d Blatter/HOStrain Rate (2 components):
     480         *
     481         * epsilon=[exx exz]
     482         *
     483         * with exz=1/2 du/dz
     484         *
     485         * the contribution of vz is neglected
     486         */
     487
     488        /*Intermediaries*/
     489        IssmDouble dvx[3];
     490
     491        /*Check that both inputs have been found*/
     492        if (!vx_input){
     493                _error_("Input missing. Here are the input pointers we have for vx: " << vx_input <<"\n");
     494        }
     495
     496        /*Get strain rate assuming that epsilon has been allocated*/
     497        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     498        epsilon[0] = dvx[0];
     499        epsilon[1] = 0.5*dvx[1];
     500
     501}/*}}}*/
    471502void Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
    472503

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

    r17236 r17242  
    7070                void       StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    7171                void       StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     72                void       StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    7273                void       StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    7374                IssmDouble TMeltingPoint(IssmDouble pressure);
    7475                void       ViscousHeatingCreateInput(void);
    7576                void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    76                 void       ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     77                void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    7778                void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
    7879                void       ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);

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

    r17226 r17242  
    451451
    452452        switch(finiteelement){
     453                case P0Enum:
     454                        /*Nodal function 1*/
     455                        dbasis[NUMNODESP0*0+0] = 0.;
     456                        dbasis[NUMNODESP0*1+0] = 0.;
     457                        return;
    453458                case P1Enum: case P1DGEnum:
    454459                        /*Nodal function 1*/

  • issm/trunk-jpl/src/c/classes/Materials/Material.h

    r16638 r17242  
    2828                virtual void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
    2929                virtual void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
     30                virtual void       GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
    3031                virtual void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon)=0;
    3132                virtual void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon)=0;

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r17236 r17242  
    351351                        /*Build viscosity: viscosity=B/(2*A^e) */
    352352                        A=exx*exx+eyy*eyy+2.*exy*exy;
     353                        if(A==0.){
     354                                /*Maxiviscositym viscosity for 0 shear areas: */
     355                                viscosity=2.5*2.e+17;
     356                        }
     357                        else{
     358                                e=(n-1.)/(2.*n);
     359                                viscosity=(1.-D)*B/(2.*pow(A,e));
     360                        }
     361                }
     362        }
     363
     364        /*Checks in debugging mode*/
     365        if(viscosity<=0) _error_("Negative viscosity");
     366        _assert_(B>0);
     367        _assert_(n>0);
     368        _assert_(D>=0 && D<1);
     369
     370        /*Return: */
     371        *pviscosity=viscosity;
     372}
     373/*}}}*/
     374/*FUNCTION Matice::GetViscosity2dverticalHO {{{*/
     375void  Matice::GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* epsilon){
     376        /*From a string tensor and a material object, return viscosity, using Glen's flow law.
     377                                                                           (1-D) B
     378          viscosity= --------------------------------------
     379                                                  2[ 2exx^2 + 2exy^2]^[(n-1)/2n]
     380
     381          where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity
     382          vector, and n the flow law exponent.
     383
     384          If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we
     385          return 10^14, initial viscosity.
     386          */
     387
     388        /*output: */
     389        IssmDouble viscosity;
     390
     391        /*input strain rate: */
     392        IssmDouble exx,exy;
     393
     394        /*Intermediary: */
     395        IssmDouble A,e;
     396        IssmDouble B,D,n;
     397
     398        /*Get B and n*/
     399        B=GetB();
     400        n=GetN();
     401        D=GetD();
     402
     403        if (n==1){
     404                /*Viscous behaviour! viscosity=B: */
     405                viscosity=(1-D)*B/2;
     406        }
     407        else{
     408                if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
     409                        viscosity=0.5*pow(10.,14);
     410                }
     411                else{
     412                        /*Retrive strain rate components: */
     413                        exx=epsilon[0];
     414                        exy=epsilon[1];
     415
     416                        /*Build viscosity: viscosity=B/(2*A^e) */
     417                        A=2*exx*exx+2.*exy*exy;
    353418                        if(A==0.){
    354419                                /*Maxiviscositym viscosity for 0 shear areas: */

  • issm/trunk-jpl/src/c/classes/Materials/Matice.h

    r16783 r17242  
    5454                void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
    5555                void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon);
     56                void       GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* pepsilon);
    5657                void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
    5758                void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon);

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r17021 r17242  
    8484                void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
    8585                void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
     86                void       GetViscosity2dverticalHO(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
    8687                void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon){_error_("not supported");};
    8788                void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon){_error_("not supported");};
Note: See TracChangeset for help on using the changeset viewer.