source: issm/oecreview/Archive/16554-17801/ISSM-17547-17548.diff

Last change on this file was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 7.9 KB

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

     
    42574257        IssmDouble  dvx[3],dvy[3],dvz[3],B,n;
    42584258        IssmDouble *xyz_list = NULL;
    42594259        IssmDouble  Jdet,r;
     4260        Gauss*      gauss = NULL;
    42604261
    42614262        parameters->FindParam(&meshtype,MeshTypeEnum);
    42624263        switch(meshtype){
    … …  
    43084309                        sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
    43094310                }
    43104311
    4311                 Gauss* gauss=element->NewGauss(5);
     4312                gauss=element->NewGauss(5);
    43124313                for(int ig=gauss->begin();ig<gauss->end();ig++){
    43134314                        gauss->GaussPoint(ig);
    43144315                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    … …  
    43734374                                                Ke,1);
    43744375
    43754376                        /*Create Right hand sides*/
    4376                         for(int ii=0;ii<tnumnodes;ii++) pe_xx[ii] += sigmapxx*tbasis[ii]*gauss->weight*Jdet;
    4377                         for(int ii=0;ii<tnumnodes;ii++) pe_yy[ii] += sigmapyy*tbasis[ii]*gauss->weight*Jdet;
    4378                         for(int ii=0;ii<tnumnodes;ii++) pe_xy[ii] += sigmapxy*tbasis[ii]*gauss->weight*Jdet;
     4377                        for(int ii=0;ii<tnumnodes;ii++) pe_xx[ii] += (r*epsxx+sigmapxx)*tbasis[ii]*gauss->weight*Jdet;
     4378                        for(int ii=0;ii<tnumnodes;ii++) pe_yy[ii] += (r*epsyy+sigmapyy)*tbasis[ii]*gauss->weight*Jdet;
     4379                        for(int ii=0;ii<tnumnodes;ii++) pe_xy[ii] += (r*epsxy+sigmapxy)*tbasis[ii]*gauss->weight*Jdet;
    43794380                        if(dim==3){
    4380                                 for(int ii=0;ii<tnumnodes;ii++) pe_zz[ii] += sigmapzz*tbasis[ii]*gauss->weight*Jdet;
    4381                                 for(int ii=0;ii<tnumnodes;ii++) pe_xz[ii] += sigmapxz*tbasis[ii]*gauss->weight*Jdet;
    4382                                 for(int ii=0;ii<tnumnodes;ii++) pe_yz[ii] += sigmapyz*tbasis[ii]*gauss->weight*Jdet;
     4381                                for(int ii=0;ii<tnumnodes;ii++) pe_zz[ii] += (r*epszz+sigmapzz)*tbasis[ii]*gauss->weight*Jdet;
     4382                                for(int ii=0;ii<tnumnodes;ii++) pe_xz[ii] += (r*epsxz+sigmapxz)*tbasis[ii]*gauss->weight*Jdet;
     4383                                for(int ii=0;ii<tnumnodes;ii++) pe_yz[ii] += (r*epsyz+sigmapyz)*tbasis[ii]*gauss->weight*Jdet;
    43834384                        }
    43844385                }
    43854386
    43864387                /*Solve the systems*/
    4387                 _error_("S");
     4388                IssmDouble* d_xx = xNew<IssmDouble>(tnumnodes);
     4389                IssmDouble* d_yy = xNew<IssmDouble>(tnumnodes);
     4390                IssmDouble* d_xy = xNew<IssmDouble>(tnumnodes);
     4391                IssmDouble* d_zz = NULL;
     4392                IssmDouble* d_xz = NULL;
     4393                IssmDouble* d_yz = NULL;
     4394                if(dim==2){
     4395                        _assert_(tnumnodes==3);
     4396                        Matrix3x3Solve(&d_xx[0],Ke,pe_xx);
     4397                        Matrix3x3Solve(&d_yy[0],Ke,pe_yy);
     4398                        Matrix3x3Solve(&d_xy[0],Ke,pe_xy);
     4399                        element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
     4400                        element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
     4401                        element->AddInput(StrainRatexyEnum,d_xy,P1DGEnum);
     4402                }
     4403                else{
     4404                        _assert_(tnumnodes==4);
     4405                        Matrix3x3Solve(&d_xx[0],Ke,pe_xx);
     4406                        Matrix3x3Solve(&d_yy[0],Ke,pe_yy);
     4407                        Matrix3x3Solve(&d_xy[0],Ke,pe_xy);
     4408                        Matrix3x3Solve(&d_zz[0],Ke,pe_zz);
     4409                        Matrix3x3Solve(&d_xz[0],Ke,pe_xz);
     4410                        Matrix3x3Solve(&d_yz[0],Ke,pe_yz);
     4411                        element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
     4412                        element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
     4413                        element->AddInput(StrainRatexyEnum,d_xy,P1DGEnum);
     4414                        element->AddInput(StrainRatezzEnum,d_zz,P1DGEnum);
     4415                        element->AddInput(StrainRatexzEnum,d_xz,P1DGEnum);
     4416                        element->AddInput(StrainRateyzEnum,d_yz,P1DGEnum);
     4417                }
    43884418
     4419                /*Update tau accordingly*/
     4420                IssmDouble* tau_xx = xNew<IssmDouble>(tnumnodes);
     4421                IssmDouble* tau_yy = xNew<IssmDouble>(tnumnodes);
     4422                IssmDouble* tau_xy = xNew<IssmDouble>(tnumnodes);
     4423                IssmDouble* tau_zz = NULL;
     4424                IssmDouble* tau_xz = NULL;
     4425                IssmDouble* tau_yz = NULL;
     4426                Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
     4427                Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
     4428                Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
     4429                Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
     4430                if(dim==3){
     4431                        epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
     4432                        epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
     4433                        epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
     4434                }
     4435                Gauss* gauss = element->NewGauss();
     4436                for(int ig=0;ig<tnumnodes;ig++){
     4437                        gauss->GaussNode(P1DGEnum,ig);
     4438
     4439                        /*Get D(u)*/
     4440                        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     4441                        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     4442                        if(dim==3){
     4443                                vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
     4444                        }
     4445                        epsxx = dvx[0];
     4446                        epsyy = dvy[1];
     4447                        epsxy = 0.5*(dvx[1] + dvy[0]);
     4448                        if(dim==3){
     4449                                epszz = dvz[2];
     4450                                epsxz = 0.5*(dvx[2] + dvz[0]);
     4451                                epsyz = 0.5*(dvy[2] + dvz[1]);
     4452                        }
     4453
     4454                        /*Get tau^(n-1) from inputs*/
     4455                        sigmapxx_input->GetInputValue(&sigmapxx,gauss);
     4456                        sigmapyy_input->GetInputValue(&sigmapyy,gauss);
     4457                        sigmapxy_input->GetInputValue(&sigmapxy,gauss);
     4458                        if(dim==3){
     4459                                sigmapzz_input->GetInputValue(&sigmapzz,gauss);
     4460                                sigmapxz_input->GetInputValue(&sigmapxz,gauss);
     4461                                sigmapyz_input->GetInputValue(&sigmapyz,gauss);
     4462                        }
     4463
     4464                        /*Get d and update tau accordingly*/
     4465                        tau_xx[ig] = sigmapxx + r*(epsxx - d_xx[ig]);
     4466                        tau_yy[ig] = sigmapyy + r*(epsyy - d_yy[ig]);
     4467                        tau_xy[ig] = sigmapxy + r*(epsxy - d_xy[ig]);
     4468                        if(dim==3){
     4469                                tau_zz[ig] = sigmapzz + r*(epszz - d_zz[ig]);
     4470                                tau_xz[ig] = sigmapxz + r*(epsxz - d_xz[ig]);
     4471                                tau_yz[ig] = sigmapyz + r*(epsyz - d_yz[ig]);
     4472                        }
     4473                }
     4474
     4475                /*Add inputs*/
     4476                element->AddInput(StrainRatexxEnum,tau_xx,P1DGEnum);
     4477                element->AddInput(StrainRateyyEnum,tau_yy,P1DGEnum);
     4478                element->AddInput(StrainRatexyEnum,tau_xy,P1DGEnum);
     4479                if(dim==3){
     4480                        element->AddInput(StrainRatezzEnum,tau_zz,P1DGEnum);
     4481                        element->AddInput(StrainRatexzEnum,tau_xz,P1DGEnum);
     4482                        element->AddInput(StrainRateyzEnum,tau_yz,P1DGEnum);
     4483                }
     4484
    43894485                /*Clean up and */
     4486                xDelete<IssmDouble>(xyz_list);
    43904487                xDelete<IssmDouble>(tbasis);
    43914488                xDelete<IssmDouble>(Ke);
    4392                 xDelete<IssmDouble>(pe_xx);
    4393                 xDelete<IssmDouble>(pe_yy);
    4394                 xDelete<IssmDouble>(pe_zz);
    4395                 xDelete<IssmDouble>(pe_xy);
    4396                 xDelete<IssmDouble>(pe_xz);
    4397                 xDelete<IssmDouble>(pe_yz);
     4489                xDelete<IssmDouble>(pe_xx); xDelete<IssmDouble>(d_xx);   xDelete<IssmDouble>(tau_xx);
     4490                xDelete<IssmDouble>(pe_yy); xDelete<IssmDouble>(d_yy);   xDelete<IssmDouble>(tau_yy);
     4491                xDelete<IssmDouble>(pe_zz); xDelete<IssmDouble>(d_zz);   xDelete<IssmDouble>(tau_zz);
     4492                xDelete<IssmDouble>(pe_xy); xDelete<IssmDouble>(d_xy);   xDelete<IssmDouble>(tau_xy);
     4493                xDelete<IssmDouble>(pe_xz); xDelete<IssmDouble>(d_xz);   xDelete<IssmDouble>(tau_xz);
     4494                xDelete<IssmDouble>(pe_yz); xDelete<IssmDouble>(d_yz);   xDelete<IssmDouble>(tau_yz);
    43984495        }
    43994496
     4497        delete gauss;
     4498
    44004499}/*}}}*/
    44014500
    44024501/*Coupling (Tiling)*/

  • ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

     
    366366        Ainv[8]=(A[0]*A[4]-A[1]*A[3])/det; /* = (a*e-b*d)/det */
    367367
    368368}/*}}}*/
     369/*FUNCTION Matrix3x3Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B) {{{*/
     370void Matrix3x3Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B){
     371
     372        IssmDouble Ainv[3][3];
     373
     374        Matrix3x3Invert(&Ainv[0][0],A);
     375        for(int i=0;i<3;i++) X[i]=Ainv[i][0]*B[0] + Ainv[i][1]*B[1] + Ainv[i][2]*B[2];
     376
     377}/*}}}*/
     378/*FUNCTION Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble *B) {{{*/
     379void Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble *B){
     380        _error_("STOP");
     381}/*}}}*/

  • ../trunk-jpl/src/c/shared/Matrix/matrix.h

     
    1414void Matrix2x2Determinant(IssmDouble* Adet,IssmDouble* A);
    1515void Matrix3x3Invert(IssmDouble* Ainv, IssmDouble* A);
    1616void Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A);
     17void Matrix3x3Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B);
     18void Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B);
    1719
    1820#endif //ifndef _MATRIXUTILS_H_
Note: See TracBrowser for help on using the repository browser.