Changeset 16892


Ignore:
Timestamp:
11/22/13 12:01:07 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: KMatrix L1L2 (not working)

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

Legend:

Unmodified
Added
Removed

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

    r16882 r16892  
    820820                case SSAApproximationEnum:
    821821                        return CreateKMatrixSSA(element);
     822                case L1L2ApproximationEnum:
     823                        return CreateKMatrixL1L2(element);
    822824                case HOApproximationEnum:
    823825                        return CreateKMatrixHO(element);
    … …  
    14411443
    14421444/*L1L2*/
     1445ElementMatrix* StressbalanceAnalysis::CreateKMatrixL1L2(Element* element){/*{{{*/
     1446
     1447        /*compute all stiffness matrices for this element*/
     1448        ElementMatrix* Ke1=CreateKMatrixL1L2Viscous(element);
     1449        ElementMatrix* Ke2=CreateKMatrixL1L2Friction(element);
     1450        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     1451
     1452        /*clean-up and return*/
     1453        delete Ke1;
     1454        delete Ke2;
     1455        return Ke;
     1456}/*}}}*/
     1457ElementMatrix* StressbalanceAnalysis::CreateKMatrixL1L2Friction(Element* element){/*{{{*/
     1458
     1459        if(!element->IsOnBed() || !element->IsFloating()) return NULL;
     1460        Element* basalelement = element->SpawnBasalElement();
     1461        ElementMatrix* Ke = CreateKMatrixSSAFriction(basalelement);
     1462
     1463        /*clean-up and return*/
     1464        basalelement->DeleteMaterials(); delete basalelement;
     1465        return Ke;
     1466}/*}}}*/
     1467ElementMatrix* StressbalanceAnalysis::CreateKMatrixL1L2Viscous(Element* element){/*{{{*/
     1468
     1469        /*Intermediaries*/
     1470        IssmDouble  viscosity,Jdet;
     1471        IssmDouble *xyz_list = NULL;
     1472
     1473        /*Get element on base*/
     1474        Element* basalelement = element->GetBasalElement()->SpawnBasalElement();
     1475
     1476        /*Fetch number of nodes and dof for this finite element*/
     1477        int numnodes = basalelement->GetNumberOfNodes();
     1478        int numdof   = numnodes*2;
     1479
     1480        /*Initialize Element matrix and vectors*/
     1481        ElementMatrix* Ke     = basalelement->NewElementMatrix(L1L2ApproximationEnum);
     1482        IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
     1483        IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
     1484        IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
     1485
     1486        /*Retrieve all inputs and parameters*/
     1487        element->GetVerticesCoordinates(&xyz_list);
     1488        Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
     1489        Input* vx_input      = element->GetInput(VxEnum);      _assert_(vx_input);
     1490        Input* vy_input      = element->GetInput(VyEnum);      _assert_(vy_input);
     1491
     1492        /* Start  looping on the number of gaussian points: */
     1493        Gauss* gauss      = element->NewGauss(5);
     1494        Gauss* gauss_base = basalelement->NewGauss();
     1495        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1496
     1497                gauss->GaussPoint(ig);
     1498                gauss->SynchronizeGaussBase(gauss_base);
     1499
     1500                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1501                this->GetBSSA(B,basalelement,xyz_list,gauss_base);
     1502                this->GetBSSAprime(Bprime,basalelement,xyz_list,gauss_base);
     1503
     1504                element->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
     1505
     1506                for(int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->weight*Jdet;
     1507
     1508                TripleMultiply(B,3,numdof,1,
     1509                                        D,3,3,0,
     1510                                        Bprime,3,numdof,0,
     1511                                        &Ke->values[0],1);
     1512        }
     1513
     1514        /*Transform Coordinate System*/
     1515        basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum);
     1516
     1517        /*Clean up and return*/
     1518        delete gauss;
     1519        delete gauss_base;
     1520        basalelement->DeleteMaterials(); delete basalelement;
     1521        xDelete<IssmDouble>(xyz_list);
     1522        xDelete<IssmDouble>(D);
     1523        xDelete<IssmDouble>(Bprime);
     1524        xDelete<IssmDouble>(B);
     1525        return Ke;
     1526}/*}}}*/
    14431527ElementVector* StressbalanceAnalysis::CreatePVectorL1L2(Element* element){/*{{{*/
    14441528

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

    r16882 r16892  
    3939                void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
    4040                /*L1L2*/
     41                ElementMatrix* CreateKMatrixL1L2(Element* element);
     42                ElementMatrix* CreateKMatrixL1L2Viscous(Element* element);
     43                ElementMatrix* CreateKMatrixL1L2Friction(Element* element);
    4144                ElementVector* CreatePVectorL1L2(Element* element);
    4245                ElementVector* CreatePVectorL1L2DrivingStress(Element* element);

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

    r16886 r16892  
    170170                virtual void   ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
    171171                virtual void   ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
     172                virtual void   ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0;
    172173                virtual void   ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
    173174                virtual void   ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;

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

    r16886 r16892  
    35403540        /*Assign output pointer*/
    35413541        *pviscosity=viscosity;
     3542}
     3543/*}}}*/
     3544/*FUNCTION Penta::ViscosityL1L2{{{*/
     3545void Penta::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){
     3546        /*Compute the L1L2 viscosity
     3547         *
     3548         *      1
     3549         * mu = - A^-1 (sigma'_e)^(1-n)
     3550         *      2
     3551         *
     3552         * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18)
     3553         *
     3554         * L1L2 assumptions:
     3555         *
     3556         * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//|
     3557         * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2
     3558         *
     3559         * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|)
     3560         *
     3561         * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0     */
     3562
     3563        int        i;
     3564        IssmDouble z,s,viscosity,p,q,delta;
     3565        IssmDouble tau_perp,tau_par,eps_b,A;
     3566        IssmDouble epsilonvx[5]; /*exx eyy exy exz eyz*/
     3567        IssmDouble epsilonvy[5]; /*exx eyy exy exz eyz*/
     3568        IssmDouble epsilon[5];   /*exx eyy exy exz eyz*/
     3569        IssmDouble z_list[NUMVERTICES];
     3570        IssmDouble slope[3];
     3571
     3572        /*Check that both inputs have been found*/
     3573        if (!vx_input || !vy_input || !surface_input) _error_("Input missing");
     3574
     3575        /*Get tau_perp*/
     3576        for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[3*i+2];
     3577        surface_input->GetInputValue(&s,gauss);
     3578        surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
     3579        PentaRef::GetInputValue(&z,&z_list[0],gauss);
     3580        tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
     3581
     3582        /* Get eps_b*/
     3583        vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
     3584        vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
     3585        for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
     3586        eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
     3587        if(eps_b==0.){
     3588                *pviscosity = 2.5e+17;
     3589                return;
     3590        }
     3591
     3592        /*Get A*/
     3593        _assert_(material->GetN()==3.0);
     3594        A=material->GetA();
     3595
     3596        /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/
     3597        p     = tau_perp *tau_perp;
     3598        q     = - eps_b/A;
     3599        delta = q *q + p*p*p*4./27.;
     3600        _assert_(delta>0);
     3601        tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.);
     3602
     3603        /*Viscosity*/
     3604        viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp));
     3605        _assert_(!xIsNan(viscosity));
     3606        _assert_(viscosity > 0.);
     3607
     3608        /*Assign output pointer*/
     3609        *pviscosity = viscosity;
    35423610}
    35433611/*}}}*/
    … …  
    68496917
    68506918                gauss->GaussPoint(ig);
    6851                 gauss->SynchronizeGaussTria(gauss_tria);
     6919                gauss->SynchronizeGaussBase(gauss_tria);
    68526920
    68536921                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    … …  
    70547122
    70557123                gauss->GaussPoint(ig);
    7056                 gauss->SynchronizeGaussTria(gauss_tria);
     7124                gauss->SynchronizeGaussBase(gauss_tria);
    70577125
    70587126                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    … …  
    74707538
    74717539                gauss->GaussPoint(ig);
    7472                 gauss->SynchronizeGaussTria(gauss_tria);
     7540                gauss->SynchronizeGaussBase(gauss_tria);
    74737541
    74747542                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    … …  
    76187686
    76197687                gauss->GaussPoint(ig);
    7620                 gauss->SynchronizeGaussTria(gauss_tria);
     7688                gauss->SynchronizeGaussBase(gauss_tria);
    76217689
    76227690                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    … …  
    76257693
    76267694                /*Get viscosity for L1L2 model*/
    7627                 GetL1L2Viscosity(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input,surf_input);
     7695                ViscosityL1L2(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input,surf_input);
    76287696
    76297697                for(i=0;i<3;i++) D[i][i]=2*viscosity*gauss->weight*Jdet;
    … …  
    95769644        delete gauss;
    95779645        return Ke;
    9578 }
    9579 /*}}}*/
    9580 /*FUNCTION Penta::GetL1L2Viscosity{{{*/
    9581 void Penta::GetL1L2Viscosity(IssmDouble* pviscosity,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input,Input* surface_input){
    9582         /*Compute the L1L2 viscosity
    9583          *
    9584          *      1
    9585          * mu = - A^-1 (sigma'_e)^(1-n)
    9586          *      2
    9587          *
    9588          * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18)
    9589          *
    9590          * L1L2 assumptions:
    9591          *
    9592          * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//|
    9593          * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2
    9594          *
    9595          * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|)
    9596          *
    9597          * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0     */
    9598 
    9599         int        i;
    9600         IssmDouble z,s,viscosity,p,q,delta;
    9601         IssmDouble tau_perp,tau_par,eps_b,A;
    9602         IssmDouble epsilonvx[5]; /*exx eyy exy exz eyz*/
    9603         IssmDouble epsilonvy[5]; /*exx eyy exy exz eyz*/
    9604         IssmDouble epsilon[5];   /*exx eyy exy exz eyz*/
    9605         IssmDouble z_list[NUMVERTICES];
    9606         IssmDouble slope[3];
    9607 
    9608         /*Check that both inputs have been found*/
    9609         if (!vx_input || !vy_input || !surface_input) _error_("Input missing");
    9610 
    9611         /*Get tau_perp*/
    9612         for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[3*i+2];
    9613         surface_input->GetInputValue(&s,gauss);
    9614         surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
    9615         PentaRef::GetInputValue(&z,&z_list[0],gauss);
    9616         tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
    9617 
    9618         /* Get eps_b*/
    9619         vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
    9620         vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
    9621         for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
    9622         eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
    9623         if(eps_b==0.){
    9624                 *pviscosity = 2.5e+17;
    9625                 return;
    9626         }
    9627 
    9628         /*Get A*/
    9629         _assert_(material->GetN()==3.0);
    9630         A=material->GetA();
    9631 
    9632         /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/
    9633         p     = tau_perp *tau_perp;
    9634         q     = - eps_b/A;
    9635         delta = q *q + p*p*p*4./27.;
    9636         _assert_(delta>0);
    9637         tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.);
    9638 
    9639         /*Viscosity*/
    9640         viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp));
    9641         _assert_(!xIsNan(viscosity));
    9642         _assert_(viscosity > 0.);
    9643 
    9644         /*Assign output pointer*/
    9645         *pviscosity = viscosity;
    9646         return;
    96479646}
    96489647/*}}}*/

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

    r16886 r16892  
    297297                void           ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    298298                void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     299                void           ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
    299300                void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    300301                void           ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
    … …  
    366367                ElementVector* CreatePVectorStressbalanceVertVolume(void);
    367368                ElementVector* CreatePVectorStressbalanceVertBase(void);
    368                 void           GetL1L2Viscosity(IssmDouble*, IssmDouble*, GaussPenta*, Input*, Input*, Input*);
    369369                #endif
    370370

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

    r16886 r16892  
    172172                void        ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");};
    173173                void        ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
     174                void        ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not implemented yet");};
    174175                void        ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    175176                void        ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");};

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

    r16890 r16892  
    329329                void           ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    330330                void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented yet");};
     331                void           ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not implemented yet");};
    331332                void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    332333                void           ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");};

  • issm/trunk-jpl/src/c/classes/gauss/Gauss.h

    r16799 r16892  
    1919                virtual void GaussVertex(int iv)=0;
    2020                virtual void GaussNode(int finitelement,int iv)=0;
     21                virtual void SynchronizeGaussBase(Gauss* gauss)=0;
    2122
    2223};

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

    r16675 r16892  
    674674}
    675675/*}}}*/
    676 /*FUNCTION GaussPenta::SynchronizeGaussTria{{{*/
    677 void GaussPenta::SynchronizeGaussTria(GaussTria* gauss_tria){
     676/*FUNCTION GaussPenta::SynchronizeGaussBase{{{*/
     677void GaussPenta::SynchronizeGaussBase(Gauss* gauss){
     678
     679        _assert_(gauss->Enum()==GaussTriaEnum);
     680        GaussTria* gauss_tria = dynamic_cast<GaussTria*>(gauss);
    678681
    679682        gauss_tria->coord1=this->coord1;

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h

    r16801 r16892  
    4848                void GaussNode(int finitelement,int iv);
    4949                void GaussFaceTria(int index1, int index2, int index3, int order);
    50                 void SynchronizeGaussTria(GaussTria* gauss_tria);
     50                void SynchronizeGaussBase(Gauss* gauss);
    5151};
    5252#endif

  • issm/trunk-jpl/src/c/classes/gauss/GaussSeg.cpp

    r16675 r16892  
    123123}
    124124/*}}}*/
     125/*FUNCTION GaussSeg::SynchronizeGaussBase{{{*/
     126void GaussSeg::SynchronizeGaussBase(Gauss* gauss){
     127
     128        _error_("not supported");
     129}
     130/*}}}*/

  • issm/trunk-jpl/src/c/classes/gauss/GaussSeg.h

    r16801 r16892  
    3434                void GaussVertex(int iv);
    3535                void GaussNode(int finitelement,int iv);
     36                void SynchronizeGaussBase(Gauss* gauss);
    3637};
    3738#endif

  • issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp

    r16675 r16892  
    503503}
    504504/*}}}*/
     505/*FUNCTION GaussTria::SynchronizeGaussBase{{{*/
     506void GaussTria::SynchronizeGaussBase(Gauss* gauss){
     507
     508        _error_("not supported");
     509}
     510/*}}}*/

  • issm/trunk-jpl/src/c/classes/gauss/GaussTria.h

    r16801 r16892  
    4444                void GaussNode(int finitelement,int iv);
    4545                void GaussEdgeCenter(int index1,int index2);
     46                void SynchronizeGaussBase(Gauss* gauss);
    4647};
    4748#endif  /* _GAUSSTRIA_H_ */
Note: See TracChangeset for help on using the changeset viewer.