Changeset 16886


Ignore:
Timestamp:
11/22/13 08:26:44 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: KMatrix Adjoint FS and HO

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

Legend:

Unmodified
Added
Removed

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

    r16876 r16886  
    127127}/*}}}*/
    128128ElementMatrix* AdjointHorizAnalysis::CreateKMatrixHO(Element* element){/*{{{*/
    129         _error_("S");
     129
     130        /*Intermediaries */
     131        bool        incomplete_adjoint;
     132        IssmDouble  Jdet,mu_prime;
     133        IssmDouble  eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij;
     134        IssmDouble  eps1[3],eps2[3],epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
     135        IssmDouble *xyz_list = NULL;
     136
     137        /*Fetch number of nodes and dof for this finite element*/
     138        int numnodes = element->GetNumberOfNodes();
     139
     140        /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
     141        element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
     142        StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
     143        ElementMatrix* Ke=analysis->CreateKMatrix(element);
     144        delete analysis;
     145        if(incomplete_adjoint) return Ke;
     146
     147        /*Retrieve all inputs and parameters*/
     148        element->GetVerticesCoordinates(&xyz_list);
     149        Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
     150        Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
     151
     152        /*Allocate dbasis*/
     153        IssmDouble* dbasis = xNew<IssmDouble>(3*numnodes);
     154
     155        /* Start  looping on the number of gaussian points: */
     156        Gauss* gauss=element->NewGauss(5);
     157        for(int ig=gauss->begin();ig<gauss->end();ig++){
     158                gauss->GaussPoint(ig);
     159
     160                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     161                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     162
     163                element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     164                element->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]);
     165                eps1[0]=2.*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
     166                eps1[1]=epsilon[2];                 eps2[1]=epsilon[0]+2.*epsilon[1];
     167                eps1[2]=epsilon[3];                 eps2[2]=epsilon[4];
     168
     169                for(int i=0;i<numnodes;i++){
     170                        for(int j=0;j<numnodes;j++){
     171                                eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i]+eps1[2]*dbasis[2*numnodes+i];
     172                                eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j]+eps1[2]*dbasis[2*numnodes+j];
     173                                eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i]+eps2[2]*dbasis[2*numnodes+i];
     174                                eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j]+eps2[2]*dbasis[2*numnodes+j];
     175
     176                                Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
     177                                Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
     178                                Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
     179                                Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
     180                        }
     181                }
     182        }
     183
     184        /*Transform Coordinate System*/
     185        element->TransformStiffnessMatrixCoord(Ke,XYEnum);
     186
     187        /*Clean up and return*/
     188        delete gauss;
     189        xDelete<IssmDouble>(dbasis);
     190        xDelete<IssmDouble>(xyz_list);
     191        return Ke;
    130192}/*}}}*/
    131193ElementMatrix* AdjointHorizAnalysis::CreateKMatrixFS(Element* element){/*{{{*/
    132         _error_("S");
     194
     195        /*Intermediaries */
     196        bool        incomplete_adjoint;
     197        IssmDouble  Jdet,mu_prime;
     198        IssmDouble  eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij,eps3dotdphii,eps3dotdphij;
     199        IssmDouble  eps1[3],eps2[3],eps3[3],epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
     200        IssmDouble *xyz_list = NULL;
     201
     202        /*Fetch number of nodes and dof for this finite element*/
     203        int vnumnodes = element->GetNumberOfNodesVelocity();
     204        int pnumnodes = element->GetNumberOfNodesPressure();
     205        int numdof    = vnumnodes*3 + pnumnodes;
     206
     207        /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
     208        element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
     209        StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
     210        ElementMatrix* Ke=analysis->CreateKMatrix(element);
     211        delete analysis;
     212        if(incomplete_adjoint) return Ke;
     213
     214        /*Retrieve all inputs and parameters*/
     215        element->GetVerticesCoordinates(&xyz_list);
     216        Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input);
     217        Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input);
     218        Input* vz_input = element->GetInput(VzEnum);_assert_(vz_input);
     219
     220        /*Allocate dbasis*/
     221        IssmDouble* dbasis = xNew<IssmDouble>(3*vnumnodes);
     222
     223        /* Start  looping on the number of gaussian points: */
     224        Gauss* gauss=element->NewGauss(5);
     225        for(int ig=gauss->begin();ig<gauss->end();ig++){
     226                gauss->GaussPoint(ig);
     227
     228                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     229                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     230
     231                element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     232                element->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]);
     233                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
     234                eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
     235                eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
     236
     237                for(int i=0;i<vnumnodes;i++){
     238                        for(int j=0;j<vnumnodes;j++){
     239                                eps1dotdphii=eps1[0]*dbasis[0*vnumnodes+i]+eps1[1]*dbasis[1*vnumnodes+i]+eps1[2]*dbasis[2*vnumnodes+i];
     240                                eps1dotdphij=eps1[0]*dbasis[0*vnumnodes+j]+eps1[1]*dbasis[1*vnumnodes+j]+eps1[2]*dbasis[2*vnumnodes+j];
     241                                eps2dotdphii=eps2[0]*dbasis[0*vnumnodes+i]+eps2[1]*dbasis[1*vnumnodes+i]+eps2[2]*dbasis[2*vnumnodes+i];
     242                                eps2dotdphij=eps2[0]*dbasis[0*vnumnodes+j]+eps2[1]*dbasis[1*vnumnodes+j]+eps2[2]*dbasis[2*vnumnodes+j];
     243                                eps3dotdphii=eps3[0]*dbasis[0*vnumnodes+i]+eps3[1]*dbasis[1*vnumnodes+i]+eps3[2]*dbasis[2*vnumnodes+i];
     244                                eps3dotdphij=eps3[0]*dbasis[0*vnumnodes+j]+eps3[1]*dbasis[1*vnumnodes+j]+eps3[2]*dbasis[2*vnumnodes+j];
     245
     246                                Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
     247                                Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
     248                                Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
     249
     250                                Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
     251                                Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
     252                                Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
     253
     254                                Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
     255                                Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
     256                                Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
     257                        }
     258                }
     259        }
     260
     261        /*Transform Coordinate System*/
     262        element->TransformStiffnessMatrixCoord(Ke,XYZEnum);
     263
     264        /*Clean up and return*/
     265        delete gauss;
     266        xDelete<IssmDouble>(dbasis);
     267        xDelete<IssmDouble>(xyz_list);
     268        return Ke;
    133269}/*}}}*/
    134270ElementVector* AdjointHorizAnalysis::CreatePVector(Element* element){/*{{{*/
    … …  
    156292        /*Intermediaries */
    157293        int        num_responses,i,meshtype,dim;
    158         IssmDouble thickness,Jdet,obs_velocity_mag,velocity_mag;
     294        IssmDouble Jdet,obs_velocity_mag,velocity_mag;
    159295        IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight;
    160296        IssmDouble scalex,scaley,scale,S;
    … …  
    355491        /*Intermediaries */
    356492        int        num_responses,i;
    357         IssmDouble thickness,Jdet,obs_velocity_mag,velocity_mag;
     493        IssmDouble Jdet,obs_velocity_mag,velocity_mag;
    358494        IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight;
    359495        IssmDouble scalex,scaley,scale,S;
    … …  
    553689        /*Intermediaries */
    554690        int         num_responses,i;
    555         IssmDouble  thickness,Jdet,obs_velocity_mag,velocity_mag;
     691        IssmDouble  Jdet,obs_velocity_mag,velocity_mag;
    556692        IssmDouble  vx,vy,vxobs,vyobs,dux,duy,weight;
    557693        IssmDouble scalex,scaley,scale,S;

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

    r16876 r16886  
    172172                virtual void   ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
    173173                virtual void   ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;
     174                virtual void   ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;
     175                virtual void   ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;
    174176                virtual void   StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
     177                virtual void   StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
     178                virtual void   StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
    175179                virtual int    VelocityInterpolation()=0;
    176180                virtual int    PressureInterpolation()=0;

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

    r16876 r16886  
    275275
    276276                /*Compute strain rate viscosity and pressure: */
    277                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     277                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    278278                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    279279                pressure_input->GetInputValue(&pressure,gauss);
    … …  
    341341
    342342                /*Compute strain rate viscosity and pressure: */
    343                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     343                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    344344                material->GetViscosity3d(&viscosity,&epsilon[0]);
    345345                pressure_input->GetInputValue(&pressure,gauss);
    … …  
    15781578
    15791579        return tau_parameter;
    1580 }
    1581 /*}}}*/
    1582 /*FUNCTION Penta::GetStrainRate3dHO{{{*/
    1583 void Penta::GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input){
    1584         /*Compute the 3d Blatter/HOStrain Rate (5 components):
    1585          *
    1586          * epsilon=[exx eyy exy exz eyz]
    1587          *
    1588          * with exz=1/2 du/dz
    1589          *      eyz=1/2 dv/dz
    1590          *
    1591          * the contribution of vz is neglected
    1592          */
    1593 
    1594         int i;
    1595         IssmDouble epsilonvx[5];
    1596         IssmDouble epsilonvy[5];
    1597 
    1598         /*Check that both inputs have been found*/
    1599         if (!vx_input || !vy_input){
    1600                 _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n");
    1601         }
    1602 
    1603         /*Get strain rate assuming that epsilon has been allocated*/
    1604         vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
    1605         vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
    1606 
    1607         /*Sum all contributions*/
    1608         for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
    1609 }
    1610 /*}}}*/
    1611 /*FUNCTION Penta::GetStrainRate3d{{{*/
    1612 void Penta::GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
    1613         /*Compute the 3d Strain Rate (6 components):
    1614          *
    1615          * epsilon=[exx eyy ezz exy exz eyz]
    1616          */
    1617 
    1618         IssmDouble epsilonvx[6];
    1619         IssmDouble epsilonvy[6];
    1620         IssmDouble epsilonvz[6];
    1621 
    1622         /*Check that both inputs have been found*/
    1623         if (!vx_input || !vy_input || !vz_input){
    1624                 _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vz: " << vz_input << "\n");
    1625         }
    1626 
    1627         /*Get strain rate assuming that epsilon has been allocated*/
    1628         vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss);
    1629         vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss);
    1630         vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss);
    1631 
    1632         /*Sum all contributions*/
    1633         for(int i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
    16341580}
    16351581/*}}}*/
    … …  
    35513497                thickness_input->GetInputValue(&thickness,gauss);
    35523498
    3553                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     3499                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    35543500                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    35553501                GetPhi(&phi, &epsilon[0], viscosity);
    … …  
    35743520
    35753521        _assert_(gauss->Enum()==GaussPentaEnum);
    3576         this->GetStrainRate3d(&epsilon[0],xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
     3522        this->StrainRateFS(&epsilon[0],xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
    35773523        material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    35783524        GetPhi(&phi,&epsilon[0],viscosity);
    … …  
    35893535        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    35903536
    3591         this->GetStrainRate3d(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     3537        this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    35923538        material->GetViscosity3dFS(&viscosity, &epsilon[0]);
    35933539
    … …  
    36033549        IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
    36043550
    3605         this->GetStrainRate3dHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     3551        this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    36063552        material->GetViscosity3d(&viscosity, &epsilon[0]);
    36073553
    36083554        /*Assign output pointer*/
    36093555        *pviscosity=viscosity;
     3556}
     3557/*}}}*/
     3558/*FUNCTION Penta::ViscosityHODerivativeEpsSquare{{{*/
     3559void Penta::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){
     3560        this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
     3561}
     3562/*}}}*/
     3563/*FUNCTION Penta::StrainRateHO{{{*/
     3564void Penta::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){
     3565        /*Compute the 3d Blatter/HOStrain Rate (5 components):
     3566         *
     3567         * epsilon=[exx eyy exy exz eyz]
     3568         *
     3569         * with exz=1/2 du/dz
     3570         *      eyz=1/2 dv/dz
     3571         *
     3572         * the contribution of vz is neglected
     3573         */
     3574
     3575        int i;
     3576        IssmDouble epsilonvx[5];
     3577        IssmDouble epsilonvy[5];
     3578
     3579        /*Check that both inputs have been found*/
     3580        if (!vx_input || !vy_input){
     3581                _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n");
     3582        }
     3583
     3584        /*Get strain rate assuming that epsilon has been allocated*/
     3585        vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
     3586        vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
     3587
     3588        /*Sum all contributions*/
     3589        for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
     3590}
     3591/*}}}*/
     3592/*FUNCTION Penta::ViscosityFSDerivativeEpsSquare{{{*/
     3593void Penta::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){
     3594        this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
     3595}
     3596/*}}}*/
     3597/*FUNCTION Penta::StrainRateFS{{{*/
     3598void Penta::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){
     3599        /*Compute the 3d Strain Rate (6 components):
     3600         *
     3601         * epsilon=[exx eyy ezz exy exz eyz]
     3602         */
     3603
     3604        IssmDouble epsilonvx[6];
     3605        IssmDouble epsilonvy[6];
     3606        IssmDouble epsilonvz[6];
     3607
     3608        /*Check that both inputs have been found*/
     3609        if (!vx_input || !vy_input || !vz_input){
     3610                _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vz: " << vz_input << "\n");
     3611        }
     3612
     3613        /*Get strain rate assuming that epsilon has been allocated*/
     3614        vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss);
     3615        vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss);
     3616        vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss);
     3617
     3618        /*Sum all contributions*/
     3619        for(int i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
    36103620}
    36113621/*}}}*/
    … …  
    44504460                GetNodalFunctionsP1(&L[0], gauss);
    44514461
    4452                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     4462                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    44534463                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    44544464                GetPhi(&phi, &epsilon[0], viscosity);
    … …  
    47204730                GetNodalFunctionsP1(&L[0], gauss);
    47214731
    4722                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     4732                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    47234733                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    47244734                GetPhi(&phi, &epsilon[0], viscosity);
    … …  
    53615371                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    53625372
    5363                 this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     5373                this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    53645374                material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    53655375                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    … …  
    54325442                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    54335443
    5434                 this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     5444                this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    54355445                material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    54365446                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
    … …  
    68456855                tria->GetBprimeSSA(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    68466856
    6847                 this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    6848                 this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
     6857                this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     6858                this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    68496859                material->GetViscosity3d(&viscosity, &epsilon[0]);
    68506860                material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    … …  
    70527062                GetBprimeSSAFS(&Bprime2[0][0], &xyz_list[0][0], gauss);
    70537063
    7054                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     7064                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    70557065                material->GetViscosity3dFS(&viscosity, &epsilon[0]);
    70567066
    … …  
    71697179                GetLprimeFSSSA(&LprimeFSSSA[0][0], &xyz_list[0][0], gauss);
    71707180
    7171                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     7181                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    71727182                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    71737183
    … …  
    74677477
    74687478                if(approximation==SSAHOApproximationEnum){
    7469                         this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    7470                         this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
     7479                        this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     7480                        this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    74717481                        material->GetViscosity3d(&viscosity, &epsilon[0]);
    74727482                        material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    … …  
    74757485                }
    74767486                else if (approximation==SSAFSApproximationEnum){
    7477                         this->GetStrainRate3d(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     7487                        this->StrainRateFS(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    74787488                        material->GetViscosity3dFS(&newviscosity,&epsilons[0]);
    74797489                }
    … …  
    77107720                GetBprimeHO(&Bprime[0], &xyz_list[0][0], gauss);
    77117721
    7712                 this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    7713                 this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
     7722                this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     7723                this->StrainRateHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    77147724                material->GetViscosity3d(&viscosity, &epsilon[0]);
    77157725                material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    … …  
    79217931                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    79227932
    7923                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     7933                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    79247934                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    79257935
    … …  
    80198029                GetBprimeFS(Bprime,&xyz_list[0][0],gauss);
    80208030
    8021                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     8031                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    80228032                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    80238033
    … …  
    82778287                vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    82788288
    8279                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     8289                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    82808290                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    82818291
    … …  
    83558365
    83568366                NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    8357                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     8367                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    83588368                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    83598369                friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
    … …  
    84398449                vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    84408450
    8441                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     8451                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    84428452                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    84438453
    … …  
    85208530
    85218531                NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    8522                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     8532                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    85238533                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    85248534                friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
    … …  
    90429052
    90439053                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    9044                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     9054                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    90459055                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    90469056
    … …  
    94469456                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    94479457
    9448                 this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     9458                this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    94499459                material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    94509460                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    … …  
    95239533                GetNodalFunctionsDerivativesVelocity(dbasis,&xyz_list[0][0],gauss);
    95249534
    9525                 //this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     9535                //this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    95269536                //material->GetViscosityDerivativeEpsSquareFS(&mu_prime,&epsilon[0]);
    95279537                //eps1[0]=epsilon[0];   eps2[0]=epsilon[3];   eps3[0]=epsilon[4];
    95289538                //eps1[1]=epsilon[3];   eps2[1]=epsilon[1];   eps3[1]=epsilon[5];
    95299539                //eps1[2]=epsilon[4];   eps2[2]=epsilon[5];   eps3[2]=epsilon[2];
    9530                 this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     9540                this->StrainRateHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    95319541                material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    95329542                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];

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

    r16876 r16886  
    298298                void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    299299                void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
     300                void           ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
     301                void           ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
    300302                void           ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");};
     303                void           StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
     304                void           StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    301305                void           StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    302306

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

    r16876 r16886  
    173173                void        ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    174174                void        ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
     175                void        ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");};
     176                void        ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");};
    175177                void        ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");};
     178                void        StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");};
     179                void        StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    176180                void        StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    177181                bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};

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

    r16876 r16886  
    330330                void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented yet");};
    331331                void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     332                void           ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");};
     333                void           ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");};
    332334                void           ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
     335                void           StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");};
     336                void           StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    333337                void           StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    334338
Note: See TracChangeset for help on using the changeset viewer.