Changeset 6723


Ignore:
Timestamp:
12/10/10 11:19:50 (14 years ago)
Author:
seroussi
Message:

started to prepare Penta for MacAyealStokes coupling

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

Legend:

Unmodified
Added
Removed

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

    r6437 r6723  
    699699}
    700700/*}}}*/
     701/*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokes{{{1*/
     702ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokes(void){
     703
     704        /*compute all stiffness matrices for this element*/
     705        ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
     706        ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     707        _error_("KMatrix coupling MacAyealStokes not terminated yet");
     708        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
     709        delete Ke1;
     710        delete Ke2;
     711        Ke1=CreateKMatrixDiagnosticMacAyeal2d();
     712        Ke2=CreateKMatrixDiagnosticStokes();
     713
     714
     715        /*Constants*/
     716        const int    numdofm=NDOF2*NUMVERTICES;
     717        const int    numdofs=NDOF4*NUMVERTICES;
     718        const int    numdoftotal=(NDOF2+NDOF4)*NUMVERTICES;
     719        int          i,j;
     720
     721        for(i=0;i<numdofs;i++) for(j=0;j<NUMVERTICES;j++){
     722                Ke->values[(i+numdofm)*numdoftotal+NDOF2*j+0]+=Ke2->values[i*numdofs+NDOF4*j+0];
     723                Ke->values[(i+numdofm)*numdoftotal+NDOF2*j+1]+=Ke2->values[i*numdofs+NDOF4*j+1];
     724        }
     725        for(i=0;i<numdofm;i++) for(j=0;j<NUMVERTICES;j++){
     726                Ke->values[i*numdoftotal+numdofm+NDOF4*j+0]+=Ke1->values[i*numdofm+NDOF2*j+0];
     727                Ke->values[i*numdoftotal+numdofm+NDOF4*j+1]+=Ke1->values[i*numdofm+NDOF2*j+1];
     728        }
     729
     730        /*clean-up and return*/
     731        delete Ke1;
     732        delete Ke2;
     733        return Ke;
     734}
     735/*}}}*/
    701736/*FUNCTION Penta::CreateKMatrixCouplingPattynStokes{{{1*/
    702737ElementMatrix* Penta::CreateKMatrixCouplingPattynStokes(void){
    … …  
    751786                case MacAyealPattynApproximationEnum:
    752787                        return CreateKMatrixDiagnosticMacAyealPattyn();
     788                case MacAyealStokesApproximationEnum:
     789                        return CreateKMatrixDiagnosticMacAyealStokes();
    753790                case PattynStokesApproximationEnum:
    754791                        return CreateKMatrixDiagnosticPattynStokes();
    … …  
    9671004}
    9681005/*}}}*/
     1006/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyealStokes{{{1*/
     1007ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyealStokes(void){
     1008
     1009        /*compute all stiffness matrices for this element*/
     1010        ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyeal2d();
     1011        ElementMatrix* Ke2=CreateKMatrixDiagnosticStokes();
     1012        ElementMatrix* Ke3=CreateKMatrixCouplingMacAyealStokes();
     1013        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
     1014
     1015        /*clean-up and return*/
     1016        delete Ke1;
     1017        delete Ke2;
     1018        delete Ke3;
     1019        return Ke;
     1020}
     1021/*}}}*/
    9691022/*FUNCTION Penta::CreateKMatrixDiagnosticPattyn{{{1*/
    9701023ElementMatrix* Penta::CreateKMatrixDiagnosticPattyn(void){
    … …  
    11121165        /*If on water or not Stokes, skip stiffness: */
    11131166        inputs->GetParameterValue(&approximation,ApproximationEnum);
    1114         if(IsOnWater() || (approximation!=StokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL;
     1167        if(IsOnWater() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL;
    11151168        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
    11161169
    … …  
    11801233        /*If on water or not Stokes, skip stiffness: */
    11811234        inputs->GetParameterValue(&approximation,ApproximationEnum);
    1182         if(IsOnWater() || IsOnShelf() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL;
     1235        if(IsOnWater() || IsOnShelf() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum)) return NULL;
    11831236        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
    11841237
    … …  
    16961749}
    16971750/*}}}*/
    1698 /*FUNCTION Penta::CreatePVectorCouplingPattynStokes {{{1*/
    1699 ElementVector* Penta::CreatePVectorCouplingPattynStokes(void){
    1700 
     1751/*FUNCTION Penta::CreatePVectorCouplingMacAyealStokes {{{1*/
     1752ElementVector* Penta::CreatePVectorCouplingMacAyealStokes(void){
     1753
     1754        _error_("coupling MacAyeal Stokes not implemented yet");
    17011755        /*compute all load vectors for this element*/
    1702         ElementVector* pe1=CreatePVectorCouplingPattynStokesViscous();
    1703         ElementVector* pe2=CreatePVectorCouplingPattynStokesFriction();
     1756        ElementVector* pe1=CreatePVectorCouplingMacAyealStokesViscous();
     1757        ElementVector* pe2=CreatePVectorCouplingMacAyealStokesFriction();
    17041758        ElementVector* pe =new ElementVector(pe1,pe2);
    17051759
    … …  
    17101764}
    17111765/*}}}*/
    1712 /*FUNCTION Penta::CreatePVectorCouplingPattynStokesViscous {{{1*/
    1713 ElementVector* Penta::CreatePVectorCouplingPattynStokesViscous(void){
     1766/*FUNCTION Penta::CreatePVectorCouplingMacAyealStokesViscous {{{1*/
     1767ElementVector* Penta::CreatePVectorCouplingMacAyealStokesViscous(void){
    17141768
    17151769        /*Constants*/
    … …  
    17311785        if(IsOnWater()) return NULL;
    17321786        inputs->GetParameterValue(&approximation,ApproximationEnum);
    1733         if(approximation!=PattynStokesApproximationEnum) return NULL;
     1787        if(approximation!=MacAyealStokesApproximationEnum) return NULL;
    17341788        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
    17351789
    … …  
    17401794        Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    17411795        Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
    1742         Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   _assert_(vzpattyn_input);
     1796        Input* vzpattyn_input=inputs->GetInput(VzMacAyealEnum);   _assert_(vzpattyn_input);
    17431797
    17441798        /* Start  looping on the number of gaussian points: */
    … …  
    17701824}
    17711825/*}}}*/
    1772 /*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{1*/
    1773 ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){
     1826/*FUNCTION Penta::CreatePVectorCouplingMacAyealStokesFriction{{{1*/
     1827ElementVector* Penta::CreatePVectorCouplingMacAyealStokesFriction(void){
    17741828
    17751829        /*Constants*/
    … …  
    17951849        if(IsOnWater() || !IsOnBed() || IsOnShelf()) return NULL;
    17961850        inputs->GetParameterValue(&approximation,ApproximationEnum);
     1851        if(approximation!=MacAyealStokesApproximationEnum) return NULL;
     1852        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     1853
     1854        /*Retrieve all inputs and parameters*/
     1855        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     1856        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     1857        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
     1858        Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
     1859        Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
     1860        Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
     1861        Input* vzpattyn_input=inputs->GetInput(VzMacAyealEnum);   _assert_(vzpattyn_input);
     1862
     1863        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     1864
     1865        /*build friction object, used later on: */
     1866        friction=new Friction("3d",inputs,matpar,analysis_type);
     1867
     1868        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     1869        gauss=new GaussPenta(0,1,2,2);
     1870        for(ig=gauss->begin();ig<gauss->end();ig++){
     1871
     1872                gauss->GaussPoint(ig);
     1873
     1874                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     1875                GetNodalFunctionsP1(l1l6, gauss);
     1876
     1877                vzpattyn_input->GetParameterValue(&w, gauss);
     1878                vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
     1879
     1880                BedNormal(&bed_normal[0],xyz_list_tria);
     1881                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     1882                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     1883                friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
     1884
     1885                for(i=0;i<NUMVERTICES2D;i++){
     1886                        pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*l1l6[i];
     1887                        pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*l1l6[i];
     1888                        pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i];
     1889                }
     1890        }
     1891
     1892        /*Clean up and return*/
     1893        delete gauss;
     1894        return pe;
     1895}
     1896/*}}}*/
     1897/*FUNCTION Penta::CreatePVectorCouplingPattynStokes {{{1*/
     1898ElementVector* Penta::CreatePVectorCouplingPattynStokes(void){
     1899
     1900        /*compute all load vectors for this element*/
     1901        ElementVector* pe1=CreatePVectorCouplingPattynStokesViscous();
     1902        ElementVector* pe2=CreatePVectorCouplingPattynStokesFriction();
     1903        ElementVector* pe =new ElementVector(pe1,pe2);
     1904
     1905        /*clean-up and return*/
     1906        delete pe1;
     1907        delete pe2;
     1908        return pe;
     1909}
     1910/*}}}*/
     1911/*FUNCTION Penta::CreatePVectorCouplingPattynStokesViscous {{{1*/
     1912ElementVector* Penta::CreatePVectorCouplingPattynStokesViscous(void){
     1913
     1914        /*Constants*/
     1915        const int   numdof=NUMVERTICES*NDOF4;
     1916
     1917        /*Intermediaries */
     1918        int         i,j,ig;
     1919        int         approximation;
     1920        double      viscosity,Jdet;
     1921        double      stokesreconditioning;
     1922        double      epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     1923        double      dw[3];
     1924        double      xyz_list[NUMVERTICES][3];
     1925        double      l1l6[6]; //for the six nodes of the penta
     1926        double      dh1dh6[3][6]; //for the six nodes of the penta
     1927        GaussPenta *gauss=NULL;
     1928
     1929        /*Initialize Element vector and return if necessary*/
     1930        if(IsOnWater()) return NULL;
     1931        inputs->GetParameterValue(&approximation,ApproximationEnum);
     1932        if(approximation!=PattynStokesApproximationEnum) return NULL;
     1933        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     1934
     1935        /*Retrieve all inputs and parameters*/
     1936        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     1937        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
     1938        Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
     1939        Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
     1940        Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
     1941        Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   _assert_(vzpattyn_input);
     1942
     1943        /* Start  looping on the number of gaussian points: */
     1944        gauss=new GaussPenta(5,5);
     1945        for (ig=gauss->begin();ig<gauss->end();ig++){
     1946
     1947                gauss->GaussPoint(ig);
     1948
     1949                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     1950                GetNodalFunctionsP1(&l1l6[0], gauss);
     1951                GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
     1952
     1953                vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
     1954
     1955                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     1956                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     1957
     1958                for(i=0;i<NUMVERTICES;i++){
     1959                        pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i];
     1960                        pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i];
     1961                        pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+2*dw[2]*dh1dh6[2][i]);
     1962                        pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i];
     1963                }
     1964        }
     1965
     1966        /*Clean up and return*/
     1967        delete gauss;
     1968        return pe;
     1969}
     1970/*}}}*/
     1971/*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{1*/
     1972ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){
     1973
     1974        /*Constants*/
     1975        const int numdof=NUMVERTICES*NDOF4;
     1976
     1977        /*Intermediaries*/
     1978        int         i,j,ig;
     1979        int         approximation,analysis_type;
     1980        double      Jdet,Jdet2d;
     1981        double      stokesreconditioning;
     1982        double     bed_normal[3];
     1983        double      epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     1984        double      viscosity, w, alpha2_gauss;
     1985        double      dw[3];
     1986        double     xyz_list_tria[NUMVERTICES2D][3];
     1987        double      xyz_list[NUMVERTICES][3];
     1988        double      l1l6[6]; //for the six nodes of the penta
     1989        Tria*       tria=NULL;
     1990        Friction*   friction=NULL;
     1991        GaussPenta  *gauss=NULL;
     1992
     1993        /*Initialize Element vector and return if necessary*/
     1994        if(IsOnWater() || !IsOnBed() || IsOnShelf()) return NULL;
     1995        inputs->GetParameterValue(&approximation,ApproximationEnum);
    17971996        if(approximation!=PattynStokesApproximationEnum) return NULL;
    17981997        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
    … …  
    18602059                case MacAyealPattynApproximationEnum:
    18612060                        return CreatePVectorDiagnosticMacAyealPattyn();
     2061                case MacAyealStokesApproximationEnum:
     2062                        return CreatePVectorDiagnosticMacAyealStokes();
    18622063                case PattynStokesApproximationEnum:
    18632064                        return CreatePVectorDiagnosticPattynStokes();
    … …  
    18782079        delete pe1;
    18792080        delete pe2;
     2081        return pe;
     2082}
     2083/*}}}*/
     2084/*FUNCTION Penta::CreatePVectorDiagnosticMacAyealStokes{{{1*/
     2085ElementVector* Penta::CreatePVectorDiagnosticMacAyealStokes(void){
     2086
     2087        /*compute all load vectors for this element*/
     2088        ElementVector* pe1=CreatePVectorDiagnosticMacAyeal();
     2089        ElementVector* pe2=CreatePVectorDiagnosticStokes();
     2090        ElementVector* pe3=CreatePVectorCouplingMacAyealStokes();
     2091        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     2092
     2093        /*clean-up and return*/
     2094        delete pe1;
     2095        delete pe2;
     2096        delete pe3;
    18802097        return pe;
    18812098}
    … …  
    20872304        if(IsOnWater()) return NULL;
    20882305        inputs->GetParameterValue(&approximation,ApproximationEnum);
    2089         if(approximation!=StokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
     2306        if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
    20902307        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
    20912308
    … …  
    21572374        if(IsOnWater() || !IsOnBed() || !IsOnShelf()) return NULL;
    21582375        inputs->GetParameterValue(&approximation,ApproximationEnum);
    2159         if(approximation!=StokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
     2376        if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
    21602377        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
    21612378
    … …  
    36033820                else if (*(iomodel->elements_type+index)==StokesApproximationEnum){
    36043821                        this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum));
     3822                }
     3823                else if (*(iomodel->elements_type+index)==MacAyealStokesApproximationEnum){
     3824                        this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealStokesApproximationEnum));
    36053825                }
    36063826                else if (*(iomodel->elements_type+index)==PattynStokesApproximationEnum){

  • TabularUnified issm/trunk/src/c/objects/Elements/Penta.h

    r6268 r6723  
    131131                ElementMatrix* CreateKMatrixCouplingMacAyealPattynViscous(void);
    132132                ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void);
     133                ElementMatrix* CreateKMatrixCouplingMacAyealStokes(void);
    133134                ElementMatrix* CreateKMatrixCouplingPattynStokes(void);
    134135                ElementMatrix* CreateKMatrixDiagnosticHoriz(void);
    … …  
    139140                ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dFriction(void);
    140141                ElementMatrix* CreateKMatrixDiagnosticMacAyealPattyn(void);
     142                ElementMatrix* CreateKMatrixDiagnosticMacAyealStokes(void);
    141143                ElementMatrix* CreateKMatrixDiagnosticPattyn(void);
    142144                ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void);
    … …  
    160162                ElementVector* CreatePVectorAdjointMacAyeal(void);
    161163                ElementVector* CreatePVectorAdjointPattyn(void);
     164                ElementVector* CreatePVectorCouplingMacAyealStokes(void);
     165                ElementVector* CreatePVectorCouplingMacAyealStokesViscous(void);
     166                ElementVector* CreatePVectorCouplingMacAyealStokesFriction(void);
    162167                ElementVector* CreatePVectorCouplingPattynStokes(void);
    163168                ElementVector* CreatePVectorCouplingPattynStokesViscous(void);
    … …  
    167172                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
    168173                ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void);
     174                ElementVector* CreatePVectorDiagnosticMacAyealStokes(void);
    169175                ElementVector* CreatePVectorDiagnosticPattyn(void);
    170176                ElementVector* CreatePVectorDiagnosticPattynStokes(void);
Note: See TracChangeset for help on using the changeset viewer.