Changeset 8409


Ignore:
Timestamp:
05/24/11 11:05:39 (14 years ago)
Author:
seroussi
Message:

moved CreateKMatrixThermalShelf from Tria to Penta

Location:
issm/trunk/src/c/objects/Elements
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r8399 r8409  
    531531
    532532}/*}}}*/
     533/*FUNCTION Penta::CreateDVectorDiagnosticHoriz {{{1*/
     534ElementVector* Penta::CreateDVectorDiagnosticHoriz(void){
     535
     536        int approximation;
     537        inputs->GetParameterValue(&approximation,ApproximationEnum);
     538
     539        switch(approximation){
     540                case StokesApproximationEnum:
     541                        return CreateDVectorDiagnosticStokes();
     542                default:
     543                        return NULL; //no need for doftypes outside of stokes approximation
     544        }
     545}
     546/*}}}*/
     547/*FUNCTION Penta::CreateDVectorDiagnosticStokes{{{1*/
     548ElementVector* Penta::CreateDVectorDiagnosticStokes(void){
     549
     550        /*output: */
     551        ElementVector* De=NULL;
     552        /*intermediary: */
     553        int approximation;
     554        int i;
     555
     556        /*Initialize Element vector and return if necessary*/
     557        inputs->GetParameterValue(&approximation,ApproximationEnum);
     558        if(approximation!=StokesApproximationEnum) return NULL;
     559
     560        De=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     561
     562        for (i=0;i<NUMVERTICES;i++){
     563                De->values[i*4+0]=VelocityEnum;
     564                De->values[i*4+1]=VelocityEnum;
     565                De->values[i*4+2]=VelocityEnum;
     566                De->values[i*4+3]=PressureEnum;
     567        }
     568
     569        return De;
     570}
     571/*}}}*/
    533572/*FUNCTION Penta::CreateKMatrix {{{1*/
    534573void  Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg, Vec df){
    … …  
    10151054}
    10161055/*}}}*/
    1017 /*FUNCTION Penta::CreateDVectorDiagnosticHoriz {{{1*/
    1018 ElementVector* Penta::CreateDVectorDiagnosticHoriz(void){
    1019 
    1020         int approximation;
    1021         inputs->GetParameterValue(&approximation,ApproximationEnum);
    1022 
    1023         switch(approximation){
    1024                 case StokesApproximationEnum:
    1025                         return CreateDVectorDiagnosticStokes();
    1026                 default:
    1027                         return NULL; //no need for doftypes outside of stokes approximation
    1028         }
    1029 }
    1030 /*}}}*/
    10311056/*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{1*/
    10321057ElementMatrix* Penta::CreateKMatrixDiagnosticHutter(void){
    … …  
    13881413        delete Ke2;
    13891414        return Ke;
    1390 }
    1391 /*}}}*/
    1392 /*FUNCTION Penta::CreateDVectorDiagnosticStokes{{{1*/
    1393 ElementVector* Penta::CreateDVectorDiagnosticStokes(void){
    1394 
    1395         /*output: */
    1396         ElementVector* De=NULL;
    1397         /*intermediary: */
    1398         int approximation;
    1399         int i;
    1400 
    1401         /*Initialize Element vector and return if necessary*/
    1402         inputs->GetParameterValue(&approximation,ApproximationEnum);
    1403         if(approximation!=StokesApproximationEnum) return NULL;
    1404 
    1405         De=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
    1406 
    1407         for (i=0;i<NUMVERTICES;i++){
    1408                 De->values[i*4+0]=VelocityEnum;
    1409                 De->values[i*4+1]=VelocityEnum;
    1410                 De->values[i*4+2]=VelocityEnum;
    1411                 De->values[i*4+3]=PressureEnum;
    1412         }
    1413 
    1414         return De;
    14151415}
    14161416/*}}}*/
    … …  
    18661866ElementMatrix* Penta::CreateKMatrixThermalShelf(void){
    18671867
     1868
     1869        /*Constants*/
     1870        const int    numdof=NDOF1*NUMVERTICES;
     1871
     1872        /*Intermediaries */
     1873        int       i,j,ig;
     1874        double    mixed_layer_capacity,thermal_exchange_velocity;
     1875        double    rho_ice,rho_water,heatcapacity;
     1876        double    Jdet2d,dt;
     1877        double    xyz_list[NUMVERTICES][3];
     1878        double   xyz_list_tria[NUMVERTICES2D][3];
     1879        double    l1l6[NUMVERTICES];
     1880        double    D_scalar;
     1881        double    Ke_gaussian[numdof][numdof]={0.0};
     1882        GaussPenta *gauss=NULL;
     1883
     1884        /*Initialize Element matrix and return if necessary*/
    18681885        if (!IsOnBed() || !IsOnShelf()) return NULL;
    1869 
    1870         /*Call Tria function*/
    1871         Tria* tria=(Tria*)SpawnTria(0,1,2);
    1872         ElementMatrix* Ke=tria->CreateKMatrixThermal();
    1873         delete tria->matice; delete tria;
    1874 
     1886        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     1887
     1888        /*Retrieve all inputs and parameters*/
     1889        this->parameters->FindParam(&dt,DtEnum);
     1890        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
     1891        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
     1892        rho_water=matpar->GetRhoWater();
     1893        rho_ice=matpar->GetRhoIce();
     1894        heatcapacity=matpar->GetHeatCapacity();
     1895        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     1896        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     1897
     1898        /* Start looping on the number of gauss (nodes on the bedrock) */
     1899        gauss=new GaussPenta(0,1,2,2);
     1900        for (ig=gauss->begin();ig<gauss->end();ig++){
     1901
     1902                gauss->GaussPoint(ig);
     1903               
     1904                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     1905                GetNodalFunctionsP1(&l1l6[0], gauss);
     1906                               
     1907                D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);
     1908                if(dt) D_scalar=dt*D_scalar;
     1909
     1910                TripleMultiply(&l1l6[0],numdof,1,0,
     1911                                        &D_scalar,1,1,0,
     1912                                        &l1l6[0],1,numdof,0,
     1913                                        &Ke_gaussian[0][0],0);
     1914
     1915                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];
     1916        }
     1917       
     1918        /*Clean up and return*/
     1919        delete gauss;
    18751920        return Ke;
    18761921}
    … …  
    26872732                water_pressure=gravity*rho_water*bed;
    26882733
    2689                 for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) pe->values[i*NDOF4+j]+=(water_pressure+damper)*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j];
     2734                for(i=0;i<NUMVERTICES;i++) for(j=0;j<3;j++) pe->values[i*NDOF4+j]+=(water_pressure+damper)*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j];
    26902735        }
    26912736

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r8399 r8409  
    15221522        }
    15231523
    1524         /*Clean up and return*/
    1525         delete gauss;
    1526         return Ke;
    1527 }
    1528 /*}}}*/
    1529 /*FUNCTION Tria::CreateKMatrixThermal {{{1*/
    1530 ElementMatrix* Tria::CreateKMatrixThermal(void){
    1531 
    1532         /*Constants*/
    1533         const int    numdof=NDOF1*NUMVERTICES;
    1534 
    1535         /*Intermediaries */
    1536         int       i,j,ig;
    1537         double    mixed_layer_capacity,thermal_exchange_velocity;
    1538         double    rho_ice,rho_water,heatcapacity;
    1539         double    Jdet,dt;
    1540         double    xyz_list[NUMVERTICES][3];
    1541         double    l1l2l3[NUMVERTICES];
    1542         double    D_scalar;
    1543         double    Ke_gaussian[numdof][numdof]={0.0};
    1544         GaussTria *gauss=NULL;
    1545 
    1546         /*Initialize Element matrix and return if necessary*/
    1547         if(!IsOnShelf()) return NULL;
    1548         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    1549 
    1550         /*Retrieve all inputs and parameters*/
    1551         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1552         this->parameters->FindParam(&dt,DtEnum);
    1553         mixed_layer_capacity=matpar->GetMixedLayerCapacity();
    1554         thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
    1555         rho_water=matpar->GetRhoWater();
    1556         rho_ice=matpar->GetRhoIce();
    1557         heatcapacity=matpar->GetHeatCapacity();
    1558 
    1559         /* Start looping on the number of gauss (nodes on the bedrock) */
    1560         gauss=new GaussTria(2);
    1561         for (ig=gauss->begin();ig<gauss->end();ig++){
    1562 
    1563                 gauss->GaussPoint(ig);
    1564                
    1565                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss);
    1566                 GetNodalFunctions(&l1l2l3[0], gauss);
    1567                                
    1568                 D_scalar=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);
    1569                 if(dt) D_scalar=dt*D_scalar;
    1570 
    1571                 TripleMultiply(&l1l2l3[0],numdof,1,0,
    1572                                         &D_scalar,1,1,0,
    1573                                         &l1l2l3[0],1,numdof,0,
    1574                                         &Ke_gaussian[0][0],0);
    1575 
    1576                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];
    1577         }
    1578        
    15791524        /*Clean up and return*/
    15801525        delete gauss;

  • issm/trunk/src/c/objects/Elements/Tria.h

    r8365 r8409  
    160160                ElementMatrix* CreateKMatrixPrognostic_DG(void);
    161161                ElementMatrix* CreateKMatrixSlope(void);
    162                 ElementMatrix* CreateKMatrixThermal(void);
    163162                ElementVector* CreatePVectorBalancethickness(void);
    164163                ElementVector* CreatePVectorBalancethickness_DG(void);
Note: See TracChangeset for help on using the changeset viewer.