Changeset 13051


Ignore:
Timestamp:
08/15/12 17:15:45 (13 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added support l1l2 equations (still under development)

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r13046 r13051  
    203203                                        ./shared/Numerics/isnan.h\
    204204                                        ./shared/Numerics/isnan.cpp\
     205                                        ./shared/Numerics/cubic.cpp\
    205206                                        ./shared/Numerics/extrema.cpp\
    206207                                        ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\

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

    r13042 r13051  
    16771677        //Need to know the type of approximation for this element
    16781678        if(iomodel->Data(FlowequationElementEquationEnum)){
    1679                 if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealApproximationEnum){
     1679                if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealApproximationEnum){
    16801680                        this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealApproximationEnum));
    16811681                }
    1682                 else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==PattynApproximationEnum){
     1682                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==PattynApproximationEnum){
    16831683                        this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum));
    16841684                }
    1685                 else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealPattynApproximationEnum){
     1685                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealPattynApproximationEnum){
    16861686                        this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum));
    16871687                }
    1688                 else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==HutterApproximationEnum){
     1688                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HutterApproximationEnum){
    16891689                        this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum));
    16901690                }
    1691                 else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==StokesApproximationEnum){
     1691                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==L1L2ApproximationEnum){
     1692                        this->inputs->AddInput(new IntInput(ApproximationEnum,L1L2ApproximationEnum));
     1693                }
     1694                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==StokesApproximationEnum){
    16921695                        this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum));
    16931696                }
    1694                 else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealStokesApproximationEnum){
     1697                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealStokesApproximationEnum){
    16951698                        this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealStokesApproximationEnum));
    16961699                }
    1697                 else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==PattynStokesApproximationEnum){
     1700                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==PattynStokesApproximationEnum){
    16981701                        this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum));
    16991702                }
    1700                 else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==NoneApproximationEnum){
     1703                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==NoneApproximationEnum){
    17011704                        this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum));
    17021705                }
    17031706                else{
    1704                         _error_("Approximation type " << EnumToStringx((int)*(iomodel->Data(FlowequationElementEquationEnum)+index)) << " not supported yet");
     1707                        _error_("Approximation type " << EnumToStringx((int)iomodel->Data(FlowequationElementEquationEnum)[index]) << " not supported yet");
    17051708                }
    17061709        }
    … …  
    28012804
    28022805        /*Intermediaries*/
    2803         IssmInt i,j;
    2804         int     penta_type;
    2805         int     penta_node_ids[6];
    2806         int     penta_vertex_ids[6];
    2807         IssmDouble  nodeinputs[6];
    2808         IssmDouble  yts;
    2809         int     stabilization;
    2810         bool    dakota_analysis;
    2811         bool    isstokes;
    2812         IssmDouble  beta,heatcapacity,referencetemperature,meltingpoint,latentheat;
     2806        IssmInt    i,j;
     2807        int        penta_type;
     2808        int        penta_node_ids[6];
     2809        int        penta_vertex_ids[6];
     2810        IssmDouble nodeinputs[6];
     2811        IssmDouble yts;
     2812        int        stabilization;
     2813        bool       dakota_analysis;
     2814        bool       isstokes;
     2815        IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat;
    28132816
    28142817        /*Fetch parameters: */
    … …  
    42324235        const int    numdof=NDOF1*NUMVERTICES;
    42334236
    4234         int          i;
    4235         int*         doflist=NULL;
    4236         IssmDouble       values[numdof];
    4237         IssmDouble       temp;
    4238         GaussPenta   *gauss=NULL;
     4237        int         i;
     4238        int        *doflist = NULL;
     4239        IssmDouble  values[numdof];
     4240        IssmDouble  temp;
     4241        GaussPenta *gauss = NULL;
    42394242
    42404243        /*Get dof list: */
    … …  
    62006203                case MacAyealApproximationEnum:
    62016204                        return CreateKMatrixDiagnosticMacAyeal2d();
     6205                case L1L2ApproximationEnum:
     6206                        return CreateKMatrixDiagnosticL1L2();
    62026207                case PattynApproximationEnum:
    62036208                        return CreateKMatrixDiagnosticPattyn();
    … …  
    62266231
    62276232        /*Intermediaries*/
    6228         int       connectivity[2];
    6229         int       i,i0,i1,j0,j1;
    6230         IssmDouble    one0,one1;
     6233        int         connectivity[2];
     6234        int         i,i0,i1,j0,j1;
     6235        IssmDouble  one0,one1;
    62316236
    62326237        /*Initialize Element matrix*/
    … …  
    64536458        delete Ke2;
    64546459        delete Ke3;
     6460        return Ke;
     6461}
     6462/*}}}*/
     6463/*FUNCTION Penta::CreateKMatrixDiagnosticL1L2{{{*/
     6464ElementMatrix* Penta::CreateKMatrixDiagnosticL1L2(void){
     6465
     6466        /*compute all stiffness matrices for this element*/
     6467        ElementMatrix* Ke1=CreateKMatrixDiagnosticL1L2Viscous();
     6468        ElementMatrix* Ke2=CreateKMatrixDiagnosticL1L2Friction();
     6469        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     6470
     6471        /*clean-up and return*/
     6472        delete Ke1;
     6473        delete Ke2;
     6474        return Ke;
     6475}
     6476/*}}}*/
     6477/*FUNCTION Penta::CreateKMatrixDiagnosticL1L2Viscous{{{*/
     6478ElementMatrix* Penta::CreateKMatrixDiagnosticL1L2Viscous(void){
     6479
     6480        /*Constants*/
     6481        const int    numdof2d=2*NUMVERTICES2D;
     6482
     6483        /*Intermediaries */
     6484        int         i,j;
     6485        IssmDouble  Jdet,viscosity;
     6486        IssmDouble  epsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
     6487        IssmDouble  xyz_list[NUMVERTICES][3];
     6488        IssmDouble  B[3][numdof2d];
     6489        IssmDouble  Bprime[3][numdof2d];
     6490        IssmDouble  Ke_gg_gaussian[numdof2d][numdof2d];
     6491        IssmDouble  D[3][3]= {0.0};                 // material matrix, simple scalar matrix.
     6492        Tria       *tria       = NULL;
     6493        Penta      *pentabase  = NULL;
     6494        GaussPenta *gauss      = NULL;
     6495        GaussTria  *gauss_tria = NULL;
     6496
     6497        /*Find penta on bed as this is a macayeal elements: */
     6498        pentabase=GetBasalElement();
     6499        tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     6500
     6501        /*Initialize Element matrix*/
     6502        ElementMatrix* Ke=new ElementMatrix(tria->nodes,NUMVERTICES2D,this->parameters,L1L2ApproximationEnum);
     6503
     6504        /*Retrieve all inputs and parameters*/
     6505        GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES);
     6506        Input* vx_input=inputs->GetInput(VxEnum);        _assert_(vx_input);
     6507        Input* vy_input=inputs->GetInput(VyEnum);        _assert_(vy_input);
     6508        Input* surf_input=inputs->GetInput(SurfaceEnum); _assert_(surf_input);
     6509
     6510        /* Start  looping on the number of gaussian points: */
     6511        gauss=new GaussPenta(5,5);
     6512        gauss_tria=new GaussTria();
     6513        for(int ig=gauss->begin();ig<gauss->end();ig++){
     6514
     6515                gauss->GaussPoint(ig);
     6516                gauss->SynchronizeGaussTria(gauss_tria);
     6517
     6518                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     6519                tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria);
     6520                tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
     6521
     6522                /*Get viscosity for L1L2 model*/
     6523                GetL1L2Viscosity(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input,surf_input);
     6524
     6525                for(i=0;i<3;i++) D[i][i]=2*viscosity*gauss->weight*Jdet;
     6526
     6527                TripleMultiply( &B[0][0],3,numdof2d,1,
     6528                                        &D[0][0],3,3,0,
     6529                                        &Bprime[0][0],3,numdof2d,0,
     6530                                        &Ke_gg_gaussian[0][0],0);
     6531                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg_gaussian[i][j];
     6532        }
     6533
     6534        /*Transform Coordinate System*/
     6535        TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,XYEnum);
     6536
     6537        /*Clean up and return*/
     6538        delete tria->matice;
     6539        delete tria;
     6540        delete gauss_tria;
     6541        delete gauss;
     6542        return Ke;
     6543}
     6544/*}}}*/
     6545/*FUNCTION Penta::CreateKMatrixDiagnosticL1L2Friction{{{*/
     6546ElementMatrix* Penta::CreateKMatrixDiagnosticL1L2Friction(void){
     6547
     6548        /*Initialize Element matrix and return if necessary*/
     6549        if(IsFloating() || !IsOnBed()) return NULL;
     6550
     6551        /*Build a tria element using the 3 nodes of the base of the penta. Then use
     6552         * the tria functionality to build a friction stiffness matrix on these 3
     6553         * nodes: */
     6554        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     6555        ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyealFriction();
     6556        delete tria->matice; delete tria;
     6557
     6558        /*clean-up and return*/
    64556559        return Ke;
    64566560}
    … …  
    71937297                case PattynApproximationEnum:
    71947298                        return CreatePVectorDiagnosticPattyn();
     7299                case L1L2ApproximationEnum:
     7300                        return CreatePVectorDiagnosticL1L2();
    71957301                case HutterApproximationEnum:
    71967302                        return NULL;
    … …  
    73467452/*FUNCTION Penta::CreatePVectorDiagnosticMacAyeal{{{*/
    73477453ElementVector* Penta::CreatePVectorDiagnosticMacAyeal(void){
     7454
     7455        if (!IsOnBed()) return NULL;
     7456
     7457        /*Call Tria function*/
     7458        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     7459        ElementVector* pe=tria->CreatePVectorDiagnosticMacAyeal();
     7460        delete tria->matice; delete tria;
     7461
     7462        /*Clean up and return*/
     7463        return pe;
     7464}
     7465/*}}}*/
     7466/*FUNCTION Penta::CreatePVectorDiagnosticL1L2{{{*/
     7467ElementVector* Penta::CreatePVectorDiagnosticL1L2(void){
    73487468
    73497469        if (!IsOnBed()) return NULL;
    … …  
    78928012        const int    numdof=NDOF2*NUMVERTICES;
    78938013
    7894         int          i;
    7895         int          approximation;
    7896         int*         doflist=NULL;
    7897         IssmDouble       vx,vy;
    7898         IssmDouble       values[numdof];
    7899         GaussPenta*  gauss;
     8014        int         i;
     8015        int         approximation;
     8016        int        *doflist        = NULL;
     8017        IssmDouble  vx,vy;
     8018        IssmDouble  values[numdof];
     8019        GaussPenta *gauss;
    79008020
    79018021        /*Get approximation enum and dof list: */
    … …  
    80438163}
    80448164/*}}}*/
     8165/*FUNCTION Penta::GetL1L2Viscosity{{{*/
     8166void Penta::GetL1L2Viscosity(IssmDouble* pviscosity,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input,Input* surface_input){
     8167        /*Compute the L1L2 viscosity
     8168         *
     8169         *      1
     8170         * mu = - A^-1 (sigma'_e)^(1-n)
     8171         *      2
     8172         *
     8173         * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18)
     8174         *
     8175         * L1L2 assumptions:
     8176         *
     8177         * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//|
     8178         * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2
     8179         *
     8180         * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|)
     8181         *
     8182         * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0     */
     8183
     8184        _error_("Not supported yet");
     8185        int        i;
     8186        IssmDouble a,c,d,z,s,viscosity;
     8187        IssmDouble tau_perp,tau_par,eps_b,A;
     8188        IssmDouble epsilonvx[5]; /*exx eyy exy exz eyz*/
     8189        IssmDouble epsilonvy[5]; /*exx eyy exy exz eyz*/
     8190        IssmDouble epsilon[5];   /*exx eyy exy exz eyz*/
     8191        IssmDouble z_list[NUMVERTICES];
     8192        IssmDouble slope[3];
     8193
     8194        /*Check that both inputs have been found*/
     8195        if (!vx_input || !vy_input || !surface_input) _error_("Input missing");
     8196
     8197        /*Get tau_perp*/
     8198        for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[3*i+2];
     8199        surface_input->GetInputValue(&s,gauss);
     8200        surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
     8201        PentaRef::GetInputValue(&z,&z_list[0],gauss);
     8202        tau_perp = matpar->GetRhoIce() * matpar->GetG() * (s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
     8203
     8204        /* Get eps_b*/
     8205        vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
     8206        vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
     8207        for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
     8208        eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
     8209        if(eps_b==0.){
     8210                *pviscosity=2.5*pow(10.,17.);
     8211                return;
     8212        }
     8213
     8214        /*Get A*/
     8215        _assert_(matice->GetN()==3.0);
     8216        A=matice->GetA();
     8217
     8218        /*Solve for tau_par (http://en.wikipedia.org/wiki/Cubic_function)*/
     8219        a=A;
     8220        c=A*tau_perp*tau_perp;
     8221        d=-eps_b;
     8222        tau_par = -1./(3.*a) * pow(0.5*(27.*a*a*d + sqrt(27.*a*a*d*27.*a*a*d -4.*pow(-3.*a*c,3.))),1./3.);
     8223
     8224        //printf("=======================================\n");
     8225        //printf("tau_par = %g (%g=0?)\n",tau_par,a*tau_par*tau_par*tau_par+c*tau_par+d);
     8226        //printf("a=%g c=%g d=%g\n",a,c,d);
     8227        //IssmDouble coeff[4];
     8228        //coeff[0]=d;
     8229        //coeff[1]=c;
     8230        //coeff[2]=0.;
     8231        //coeff[3]=a;
     8232        //int numroots;
     8233        //IssmDouble roots[3];
     8234        //cubic(coeff,roots,&numroots);
     8235        //tau_par=roots[0];
     8236        //for(i=0;i<numroots;i++) printf(" %g ",roots[i]);
     8237        //printf(" (%g =0?)\n",a*roots[0]*roots[0]*roots[0]+c*roots[0]+d);
     8238       
     8239
     8240        /*Viscosity*/
     8241        viscosity = 1./(2.*A)*pow(tau_par*tau_par + tau_perp*tau_perp ,-2.);
     8242        _assert_(!isnan(viscosity));
     8243
     8244        /*Assign output pointer*/
     8245        *pviscosity = viscosity;
     8246        return;
     8247}
     8248/*}}}*/
    80458249/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{*/
    80468250void  Penta::InputUpdateFromSolutionDiagnosticHoriz(IssmDouble* solution){
    … …  
    80618265                        return;
    80628266                }
     8267        }
     8268        if (approximation==L1L2ApproximationEnum){
     8269                if (!IsOnBed()) return;
     8270                InputUpdateFromSolutionDiagnosticL1L2(solution);
     8271                return;
    80638272        }
    80648273        else if (approximation==PattynApproximationEnum){
    … …  
    83428551        xDelete<int>(doflistm);
    83438552        xDelete<int>(doflists);
     8553}
     8554/*}}}*/
     8555/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticL1L2 {{{*/
     8556void  Penta::InputUpdateFromSolutionDiagnosticL1L2(IssmDouble* solution){
     8557
     8558        const int    numdof=NDOF2*NUMVERTICES;
     8559
     8560        int     i;
     8561        IssmDouble  rho_ice,g;
     8562        IssmDouble  values[numdof];
     8563        IssmDouble  vx[NUMVERTICES];
     8564        IssmDouble  vy[NUMVERTICES];
     8565        IssmDouble  vz[NUMVERTICES];
     8566        IssmDouble  vel[NUMVERTICES];
     8567        IssmDouble  pressure[NUMVERTICES];
     8568        IssmDouble  surface[NUMVERTICES];
     8569        IssmDouble  xyz_list[NUMVERTICES][3];
     8570        int    *doflist = NULL;
     8571        Penta  *penta   = NULL;
     8572
     8573        /*Get dof list: */
     8574        GetDofList(&doflist,L1L2ApproximationEnum,GsetEnum);
     8575
     8576        /*Use the dof list to index into the solution vector: */
     8577        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
     8578
     8579        /*Transform solution in Cartesian Space*/
     8580        TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
     8581
     8582        /*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
     8583        for(i=0;i<3;i++){
     8584                vx[i]  =values[i*NDOF2+0];
     8585                vy[i]  =values[i*NDOF2+1];
     8586                vx[i+3]=vx[i];
     8587                vy[i+3]=vy[i];
     8588
     8589                /*Check solution*/
     8590                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
     8591                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     8592        }
     8593
     8594        /*Get parameters fro pressure computation*/
     8595        rho_ice=matpar->GetRhoIce();
     8596        g=matpar->GetG();
     8597
     8598        /*Start looping over all elements above current element and update all inputs*/
     8599        penta=this;
     8600        for(;;){
     8601
     8602                /*Get node data: */
     8603                GetVerticesCoordinates(&xyz_list[0][0],penta->nodes,NUMVERTICES);
     8604
     8605                /*Now Compute vel*/
     8606                GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
     8607                for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     8608
     8609                /*Now compute pressure*/
     8610                GetInputListOnVertices(&surface[0],SurfaceEnum);
     8611                for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
     8612
     8613                /*Now, we have to move the previous Vx and Vy inputs  to old
     8614                 * status, otherwise, we'll wipe them off: */
     8615                penta->inputs->ChangeEnum(VxEnum,VxPicardEnum);
     8616                penta->inputs->ChangeEnum(VyEnum,VyPicardEnum);
     8617                penta->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
     8618
     8619                /*Add vx and vy as inputs to the tria element: */
     8620                penta->inputs->AddInput(new PentaP1Input(VxEnum,vx));
     8621                penta->inputs->AddInput(new PentaP1Input(VyEnum,vy));
     8622                penta->inputs->AddInput(new PentaP1Input(VelEnum,vel));
     8623                penta->inputs->AddInput(new PentaP1Input(PressureEnum,pressure));
     8624
     8625                /*Stop if we have reached the surface*/
     8626                if (penta->IsOnSurface()) break;
     8627
     8628                /* get upper Penta*/
     8629                penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
     8630        }
     8631
     8632        /*Free ressources:*/
     8633        xDelete<int>(doflist);
    83448634}
    83458635/*}}}*/

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

    r13036 r13051  
    227227                ElementMatrix* CreateKMatrixDiagnosticMacAyealPattyn(void);
    228228                ElementMatrix* CreateKMatrixDiagnosticMacAyealStokes(void);
     229                ElementMatrix* CreateKMatrixDiagnosticL1L2(void);
     230                ElementMatrix* CreateKMatrixDiagnosticL1L2Viscous(void);
     231                ElementMatrix* CreateKMatrixDiagnosticL1L2Friction(void);
    229232                ElementMatrix* CreateKMatrixDiagnosticPattyn(void);
    230233                ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void);
    … …  
    245248                void           InputUpdateFromSolutionDiagnosticMacAyealPattyn( IssmDouble* solutiong);
    246249                void           InputUpdateFromSolutionDiagnosticMacAyealStokes( IssmDouble* solutiong);
     250                void           InputUpdateFromSolutionDiagnosticL1L2( IssmDouble* solutiong);
    247251                void           InputUpdateFromSolutionDiagnosticPattyn( IssmDouble* solutiong);
    248252                void           InputUpdateFromSolutionDiagnosticPattynStokes( IssmDouble* solutiong);
    … …  
    265269                ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void);
    266270                ElementVector* CreatePVectorDiagnosticMacAyealStokes(void);
     271                ElementVector* CreatePVectorDiagnosticL1L2(void);
    267272                ElementVector* CreatePVectorDiagnosticPattyn(void);
    268273                ElementVector* CreatePVectorDiagnosticPattynStokes(void);
    … …  
    273278                ElementVector* CreatePVectorDiagnosticVertVolume(void);
    274279                ElementVector* CreatePVectorDiagnosticVertBase(void);
     280                void GetL1L2Viscosity(IssmDouble*, IssmDouble*, GaussPenta*, Input*, Input*, Input*);
    275281                #endif
    276282

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

    r13036 r13051  
    232232        /*Intermediaries */
    233233        int        i,j,ig;
    234         IssmDouble     heatcapacity,latentheat;
    235         IssmDouble     Jdet,D_scalar;
    236         IssmDouble     xyz_list[NUMVERTICES][3];
    237         IssmDouble     L[3];
     234        IssmDouble heatcapacity,latentheat;
     235        IssmDouble Jdet,D_scalar;
     236        IssmDouble xyz_list[NUMVERTICES][3];
     237        IssmDouble L[3];
    238238        GaussTria *gauss=NULL;
    239239
    … …  
    14441444
    14451445        /*Intermediaries*/
    1446         int    i,j;
    1447         int    tria_vertex_ids[3];
     1446        int        i,j;
     1447        int        tria_vertex_ids[3];
    14481448        IssmDouble nodeinputs[3];
    14491449        IssmDouble cmmininputs[3];
    14501450        IssmDouble cmmaxinputs[3];
    1451         bool   control_analysis=false;
    1452         int    num_control_type;
     1451        bool       control_analysis   = false;
     1452        int        num_control_type;
    14531453        IssmDouble yts;
    1454         int    num_cm_responses;
     1454        int        num_cm_responses;
    14551455   
    14561456        /*Get parameters: */

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

    r13036 r13051  
    128128}
    129129/*}}}*/
     130/*FUNCTION Matice::GetA {{{*/
     131IssmDouble Matice::GetA(){
     132        /*
     133         * A = 1/B^n
     134         */
     135
     136        IssmDouble B,n;
     137
     138        inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
     139        n=this->GetN();
     140
     141        return pow(B,-n);
     142}
     143/*}}}*/
    130144/*FUNCTION Matice::GetB {{{*/
    131145IssmDouble Matice::GetB(){
    … …  
    236250                        if(A==0){
    237251                                /*Maxiviscositym viscosity for 0 shear areas: */
    238                                 viscosity=2.5*pow((IssmDouble)10,(IssmDouble)17);
     252                                viscosity=2.5*pow(10.,17.);
    239253                        }
    240254                        else{

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

    r12472 r13051  
    5858                /*}}}*/
    5959                /*Matice Numerics: {{{*/
    60                 void   SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin);
    61                 void   GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
    62                 void   GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
    63                 void   GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon);
    64                 void   GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
    65                 void   GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    66                 void   GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
     60                void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
     61                void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
     62                void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
     63                void       GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon);
     64                void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
     65                void       GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
     66                void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
     67                IssmDouble GetA();
    6768                IssmDouble GetB();
    6869                IssmDouble GetBbar();
    6970                IssmDouble GetN();
    70                 bool   IsInput(int name);
     71                bool       IsInput(int name);
    7172                /*}}}*/
    7273};

  • issm/trunk-jpl/src/c/classes/objects/Node.cpp

    r13036 r13051  
    9292                                for(k=1;k<=gsize;k++) this->FreezeDof(k);
    9393                        }
     94                        if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==L1L2ApproximationEnum && !iomodel->Data(MeshVertexonbedEnum)[io_index]){
     95                                for(k=1;k<=gsize;k++) this->FreezeDof(k);
     96                        }
    9497                        if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealPattynApproximationEnum && iomodel->Data(FlowequationBordermacayealEnum)[io_index]){
    9598                                if(!iomodel->Data(MeshVertexonbedEnum)[io_index]){
    … …  
    114117        /*Diagnostic Hutter*/
    115118        if (analysis_type==DiagnosticHutterAnalysisEnum){
    116                 if (!iomodel->Data(FlowequationVertexEquationEnum)) _error_("iomodel->vertices_type is NULL");
     119                _assert_(iomodel->Data(FlowequationVertexEquationEnum));
    117120                /*Constrain all nodes that are not Hutter*/
    118121                if (reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[io_index])!=HutterApproximationEnum){

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r12832 r13051  
    3636        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIshutterEnum));
    3737        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsmacayealpattynEnum));
     38        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsl1l2Enum));
    3839        parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsstokesEnum));
    3940        parameters->AddObject(iomodel->CopyConstantObject(SettingsOutputFrequencyEnum));

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp

    r13036 r13051  
    2020        double  rho_ice;
    2121        double  stokesreconditioning;
    22         bool    isstokes,ismacayealpattyn;
     22        bool    isstokes,isl1l2,ismacayealpattyn;
    2323   bool    spcpresent=false;
    2424        int Mx,Nx;
    … …  
    5656        iomodel->Constant(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
    5757        iomodel->Constant(&isstokes,FlowequationIsstokesEnum);
     58        iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum);
    5859        iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
    5960
    … …  
    6566       
    6667        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
    67         if (!ismacayealpattyn & !isstokes){
     68        if(!ismacayealpattyn & !isstokes & !isl1l2){
    6869                *pconstraints=constraints;
    6970                return;

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

    r12832 r13051  
    9494                        count++;
    9595                }
     96                else if ((int)*(elements_type+element)==(L1L2ApproximationEnum)){
     97                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,L1L2IceFrontEnum,DiagnosticHorizAnalysisEnum));
     98                        count++;
     99                }
    96100                else if ((int)*(elements_type+element)==(StokesApproximationEnum)){
    97101                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,StokesIceFrontEnum,DiagnosticHorizAnalysisEnum));

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp

    r12832 r13051  
    1919        bool   continuous_galerkin=true;
    2020        int    numberofvertices;
    21         bool   isstokes,ismacayealpattyn;
     21        bool   isstokes,isl1l2,ismacayealpattyn;
    2222
    2323        /*DataSets: */
    … …  
    2727        iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
    2828        iomodel->Constant(&isstokes,FlowequationIsstokesEnum);
     29        iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum);
    2930        iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
    3031
    … …  
    3637       
    3738        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
    38         if(!ismacayealpattyn & !isstokes){
     39        if(!ismacayealpattyn & !isstokes & !isl1l2){
    3940                *pnodes=nodes;
    4041                return;

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp

    r12832 r13051  
    1919        int    numberofelements;
    2020        bool   ismacayealpattyn;
     21        bool   isl1l2;
    2122        bool   isstokes;
    2223        bool   control_analysis;
    … …  
    2526        /*Fetch constants needed: */
    2627        iomodel->Constant(&isstokes,FlowequationIsstokesEnum);
     28        iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum);
    2729        iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
    2830        iomodel->Constant(&dim,MeshDimensionEnum);
    … …  
    3234
    3335        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
    34         if(!ismacayealpattyn & !isstokes) return;
     36        if(!ismacayealpattyn & !isstokes &!isl1l2) return;
    3537
    3638        /*Fetch data needed: */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp

    r13036 r13051  
    1919        if (analysis_type==DiagnosticHorizAnalysisEnum){
    2020                if (vertices_type[0]==MacAyealApproximationEnum){
     21                        numdofs=2;
     22                }
     23                else if (vertices_type[0]==L1L2ApproximationEnum){
    2124                        numdofs=2;
    2225                }

  • issm/trunk-jpl/src/c/solutions/convergence.cpp

    r13036 r13051  
    8181        //print
    8282        if(res<eps_res){
    83                 if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   mechanical equilibrium convergence criterion" << res*100 << " < " << eps_res*100 << " %");
     83                if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<"   mechanical equilibrium convergence criterion"<<res*100<< " < "<<eps_res*100<<" %");
    8484                converged=true;
    8585        }
    8686        else{
    87                 if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   mechanical equilibrium convergence criterion" << res*100 << " > " << eps_res*100 << " %");
     87                if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<"   mechanical equilibrium convergence criterion"<<res*100<<" > "<<eps_res*100<<" %");
    8888                converged=false;
    8989        }

  • issm/trunk-jpl/src/c/solutions/diagnostic_core.cpp

    r12832 r13051  
    2020        bool  ishutter          = false;
    2121        bool  ismacayealpattyn  = false;
     22        bool  isl1l2            = false;
    2223        bool  isstokes          = false;
    2324        bool  isnewton          = false;
    … …  
    3334        femmodel->parameters->FindParam(&ishutter,FlowequationIshutterEnum);
    3435        femmodel->parameters->FindParam(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
     36        femmodel->parameters->FindParam(&isl1l2,FlowequationIsl1l2Enum);
    3537        femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum);
    3638        femmodel->parameters->FindParam(&isnewton,DiagnosticIsnewtonEnum);
    … …  
    7072        }
    7173
    72         if (ismacayealpattyn ^ isstokes){ // ^ = xor
     74        if ((ismacayealpattyn || isl1l2) ^ isstokes){ // ^ = xor
    7375               
    7476                if(VerboseSolution()) _pprintLine_("   computing velocities");
Note: See TracChangeset for help on using the changeset viewer.