Changeset 18385


Ignore:
Timestamp:
08/13/14 14:47:40 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added maximum stress eigen value

Location:
issm/trunk-jpl/src
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

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

    r18363 r18385  
    11581158                                input=this->inputs->GetInput(output_enum);
    11591159                                break;
     1160                        case StressMaxPrincipalEnum:
     1161                                this->StressMaxPrincipalCreateInput();
     1162                                input=this->inputs->GetInput(output_enum);
     1163                                break;
    11601164                        case StressTensorxxEnum:
    11611165                        case StressTensorxyEnum:
    … …  
    13381342        return this->matpar->TMeltingPoint(pressure);
    13391343}/*}}}*/
     1344void       Element::StressMaxPrincipalCreateInput(void){/*{{{*/
     1345
     1346        /*Intermediaries*/
     1347        IssmDouble *xyz_list = NULL;
     1348        IssmDouble  sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz;
     1349        IssmDouble  a,b,c,d,x[3],max;
     1350        int         dim,numroots;
     1351
     1352        /*First: get stress tensor*/
     1353        this->ComputeStressTensor();
     1354
     1355        /*Get domain dimension*/
     1356        this->FindParam(&dim,DomainDimensionEnum);
     1357
     1358        /*Fetch number vertices and allocate memory*/
     1359        int         numvertices  = this->GetNumberOfVertices();
     1360        IssmDouble* maxprincipal = xNew<IssmDouble>(numvertices);
     1361
     1362        /*Retrieve all inputs and parameters*/
     1363        this->GetVerticesCoordinatesBase(&xyz_list);
     1364        Input* sigma_xx_input  = this->GetInput(StressTensorxxEnum); _assert_(sigma_xx_input);
     1365        Input* sigma_yy_input  = this->GetInput(StressTensoryyEnum); _assert_(sigma_yy_input);
     1366        Input* sigma_xy_input  = this->GetInput(StressTensorxyEnum); _assert_(sigma_xy_input);
     1367        Input* sigma_xz_input  = NULL;
     1368        Input* sigma_yz_input  = NULL;
     1369        Input* sigma_zz_input  = NULL;
     1370        if(dim==3){
     1371                sigma_xz_input  = this->GetInput(StressTensorxzEnum); _assert_(sigma_xz_input);
     1372                sigma_yz_input  = this->GetInput(StressTensoryzEnum); _assert_(sigma_yz_input);
     1373                sigma_zz_input  = this->GetInput(StressTensorzzEnum); _assert_(sigma_zz_input);
     1374        }
     1375
     1376        /*loop over vertices: */
     1377        Gauss* gauss=this->NewGauss();
     1378        for (int iv=0;iv<numvertices;iv++){
     1379                gauss->GaussVertex(iv);
     1380
     1381                sigma_xx_input->GetInputValue(&sigma_xx,gauss);
     1382                sigma_yy_input->GetInputValue(&sigma_yy,gauss);
     1383                sigma_xy_input->GetInputValue(&sigma_xy,gauss);
     1384                if(dim==3){
     1385                        sigma_xz_input->GetInputValue(&sigma_xz,gauss);
     1386                        sigma_yz_input->GetInputValue(&sigma_yz,gauss);
     1387                        sigma_zz_input->GetInputValue(&sigma_zz,gauss);
     1388                }
     1389
     1390                if(dim==2){
     1391                        a = 0.;
     1392                        b = 1.;
     1393                        c = -sigma_yy -sigma_xx;
     1394                        d = sigma_xx*sigma_yy - sigma_xy*sigma_xy;
     1395                }
     1396                else{
     1397                        a = -1.;
     1398                        b = sigma_xx+sigma_yy+sigma_zz;
     1399                        c = -sigma_xx*sigma_yy -sigma_xx*sigma_zz -sigma_yy*sigma_zz + sigma_xy*sigma_xy +sigma_xz*sigma_xz +sigma_yz*sigma_yz;
     1400                        d = sigma_xx*sigma_yy*sigma_zz - sigma_xx*sigma_yz*sigma_yz -sigma_yy*sigma_xz*sigma_xz - sigma_zz*sigma_xy*sigma_xy + 2.*sigma_xy*sigma_xz*sigma_yz;
     1401                }
     1402
     1403                /*Get roots of polynomials*/
     1404                cubic(a,b,c,d,x,&numroots);
     1405
     1406                /*Initialize maximum eigne value*/
     1407                if(numroots>0){
     1408                        max = x[0];
     1409                }
     1410                else{
     1411                        _error_("No eigen value found");
     1412                }
     1413
     1414                /*Get max*/
     1415                for(int i=1;i<numroots;i++){
     1416                        if(x[i]>max) max = x[i];
     1417                }
     1418
     1419                maxprincipal[iv]=max;
     1420        }
     1421
     1422        /*Create input*/
     1423        this->AddInput(StressMaxPrincipalEnum,maxprincipal,P1Enum);
     1424
     1425        /*Clean up and return*/
     1426        xDelete<IssmDouble>(maxprincipal);
     1427        delete gauss;
     1428}
     1429/*}}}*/
    13401430void       Element::ViscousHeatingCreateInput(void){/*{{{*/
    13411431

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

    r18363 r18385  
    149149                void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array);
    150150                void       ViscousHeatingCreateInput(void);
     151                void       StressMaxPrincipalCreateInput(void);
    151152                void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    152153                void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);

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

    r18308 r18385  
    245245void       Penta::ComputeStressTensor(){/*{{{*/
    246246
    247         IssmDouble      xyz_list[NUMVERTICES][3];
    248         IssmDouble      pressure,viscosity;
    249         IssmDouble      epsilon[6]; /* epsilon=[exx,eyy,exy];*/
    250         IssmDouble      sigma_xx[NUMVERTICES];
    251         IssmDouble              sigma_yy[NUMVERTICES];
    252         IssmDouble              sigma_zz[NUMVERTICES];
    253         IssmDouble      sigma_xy[NUMVERTICES];
    254         IssmDouble              sigma_xz[NUMVERTICES];
    255         IssmDouble              sigma_yz[NUMVERTICES];
     247        IssmDouble  xyz_list[NUMVERTICES][3];
     248        IssmDouble  pressure,viscosity;
     249        IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
     250        IssmDouble  sigma_xx[NUMVERTICES];
     251        IssmDouble      sigma_yy[NUMVERTICES];
     252        IssmDouble      sigma_zz[NUMVERTICES];
     253        IssmDouble  sigma_xy[NUMVERTICES];
     254        IssmDouble      sigma_xz[NUMVERTICES];
     255        IssmDouble      sigma_yz[NUMVERTICES];
    256256        GaussPenta* gauss=NULL;
    257257

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r18383 r18385  
    563563        StressTensoryzEnum,
    564564        StressTensorzzEnum,
     565        StressMaxPrincipalEnum,
    565566        DeviatoricStressEnum,
    566567        DeviatoricStressxxEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r18383 r18385  
    554554                case StressTensoryzEnum : return "StressTensoryz";
    555555                case StressTensorzzEnum : return "StressTensorzz";
     556                case StressMaxPrincipalEnum : return "StressMaxPrincipal";
    556557                case DeviatoricStressEnum : return "DeviatoricStress";
    557558                case DeviatoricStressxxEnum : return "DeviatoricStressxx";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r18383 r18385  
    566566              else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
    567567              else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
     568              else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum;
    568569              else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
    569570              else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
    … …  
    628629              else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
    629630              else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
    630               else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"MinVel")==0) return MinVelEnum;
     634              if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
     635              else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    635636              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
    636637              else if (strcmp(name,"MinVx")==0) return MinVxEnum;

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r18383 r18385  
    546546def StressTensoryzEnum(): return StringToEnum("StressTensoryz")[0]
    547547def StressTensorzzEnum(): return StringToEnum("StressTensorzz")[0]
     548def StressMaxPrincipalEnum(): return StringToEnum("StressMaxPrincipal")[0]
    548549def DeviatoricStressEnum(): return StringToEnum("DeviatoricStress")[0]
    549550def DeviatoricStressxxEnum(): return StringToEnum("DeviatoricStressxx")[0]
Note: See TracChangeset for help on using the changeset viewer.