Changeset 7070


Ignore:
Timestamp:
01/13/11 16:03:17 (14 years ago)
Author:
seroussi
Message:

added coupling macayeal/stokes friction

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

Legend:

Unmodified
Added
Removed

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

    r7027 r7070  
    704704        /*compute all stiffness matrices for this element*/
    705705        ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealStokesViscous();
    706         ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealPattynFriction();
     706        ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealStokesFriction();
    707707        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    708708
    … …  
    804804}
    805805/*}}}*/
    806 /*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokesFriction{{{1*/
     806/*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokesFriction {{{1*/
    807807ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokesFriction(void){
    808808
    809         /*Initialize Element matrix and return if necessary*/
     809        /*Constants*/
     810        const int numdof=NUMVERTICES*NDOF4;
     811        const int numdofm=NUMVERTICES*NDOF2;
     812        const int numdof2d=NUMVERTICES2D*NDOF4;
     813        const int numdof2dm=NUMVERTICES2D*NDOF2;
     814        const int numdoftot=numdof+numdofm;
     815
     816        /*Intermediaries */
     817        int        i,j,ig;
     818        int        analysis_type,approximation;
     819        double     stokesreconditioning;
     820        double     viscosity,alpha2_gauss,Jdet2d;
     821        double    bed_normal[3];
     822        double     epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     823        double     xyz_list[NUMVERTICES][3];
     824        double    xyz_list_tria[NUMVERTICES2D][3];
     825        double     LMacAyealStokes[8][numdof2dm];
     826        double     LprimeMacAyealStokes[8][numdof2d];
     827        double     DLMacAyealStokes[8][8]={0.0};
     828        double     LStokesMacAyeal[4][numdof2d];
     829        double     LprimeStokesMacAyeal[4][numdof2dm];
     830        double     DLStokesMacAyeal[4][4]={0.0};
     831        double     Ke_drag_gaussian[numdof2dm][numdof2d];
     832        double     Ke_drag_gaussian2[numdof2d][numdof2dm];
     833        Friction*  friction=NULL;
     834        GaussPenta *gauss=NULL;
     835
     836        /*If on water or not Stokes, skip stiffness: */
     837        inputs->GetParameterValue(&approximation,ApproximationEnum);
    810838        if(IsOnWater() || IsOnShelf() || !IsOnBed()) return NULL;
    811 
    812         Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    813         ElementMatrix* Ke=tria->CreateKMatrixCouplingMacAyealPattynFriction();
    814         delete tria->matice; delete tria;
    815 
     839        ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
     840        ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     841        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
     842        delete Ke1; delete Ke2;
     843
     844        /*Retrieve all inputs and parameters*/
     845        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     846        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     847        parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
     848        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     849        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     850        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     851        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     852
     853        /*build friction object, used later on: */
     854        friction=new Friction("3d",inputs,matpar,analysis_type);
     855
     856        /* Start  looping on the number of gaussian points: */
     857        gauss=new GaussPenta(0,1,2,2);
     858        for (ig=gauss->begin();ig<gauss->end();ig++){
     859
     860                gauss->GaussPoint(ig);
     861
     862                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
     863                GetLMacAyealStokes(&LMacAyealStokes[0][0], gauss);
     864                GetLprimeMacAyealStokes(&LprimeMacAyealStokes[0][0], &xyz_list[0][0], gauss);
     865                GetLStokesMacAyeal(&LStokesMacAyeal[0][0], gauss);
     866                GetLprimeStokesMacAyeal(&LprimeStokesMacAyeal[0][0], &xyz_list[0][0], gauss);
     867
     868                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     869                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     870
     871                BedNormal(&bed_normal[0],xyz_list_tria);
     872                friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
     873
     874                DLMacAyealStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
     875                DLMacAyealStokes[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
     876                DLMacAyealStokes[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
     877                DLMacAyealStokes[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
     878                DLMacAyealStokes[4][4]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
     879                DLMacAyealStokes[5][5]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
     880                DLMacAyealStokes[6][6]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[0];
     881                DLMacAyealStokes[7][7]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[1];
     882
     883                DLStokesMacAyeal[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
     884                DLStokesMacAyeal[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
     885                DLStokesMacAyeal[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
     886                DLStokesMacAyeal[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
     887               
     888                TripleMultiply( &LMacAyealStokes[0][0],8,numdof2dm,1,
     889                                        &DLMacAyealStokes[0][0],8,8,0,
     890                                        &LprimeMacAyealStokes[0][0],8,numdof2d,0,
     891                                        &Ke_drag_gaussian[0][0],0);
     892
     893                TripleMultiply( &LStokesMacAyeal[0][0],4,numdof2d,1,
     894                                        &DLStokesMacAyeal[0][0],4,4,0,
     895                                        &LprimeStokesMacAyeal[0][0],4,numdof2dm,0,
     896                                        &Ke_drag_gaussian2[0][0],0);
     897
     898                for(i=0;i<numdof2dm;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdoftot+j+numdofm]+=Ke_drag_gaussian[i][j];
     899                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2dm;j++) Ke->values[(i+numdofm)*numdoftot+j]+=Ke_drag_gaussian2[i][j];
     900        }
     901
     902        /*Clean up and return*/
     903        delete gauss;
     904        delete friction;
    816905        return Ke;
    817906}

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

    r7000 r7070  
    713713}
    714714/*}}}*/
     715/*FUNCTION PentaRef::GetLMacAyealStokes {{{1*/
     716void PentaRef::GetLMacAyealStokes(double* LStokes, GaussPenta* gauss){
     717        /*
     718         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     719         * For grid i, Li can be expressed in the actual coordinate system
     720         * by:
     721         *       Li=[ h    0 ]
     722         *                    [ 0    h ]
     723         *                    [ h    0 ]
     724         *                    [ 0    h ]
     725         *                    [ h    0 ]
     726         *                    [ 0    h ]
     727         *                    [ h    0 ]
     728         *                    [ 0    h ]
     729         * where h is the interpolation function for grid i.
     730         */
     731
     732        int i;
     733        int num_dof=2;
     734
     735        double l1l2l3[NUMNODESP1_2d];
     736
     737
     738        /*Get l1l2l3 in actual coordinate system: */
     739        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     740        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     741        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     742
     743        /*Build LStokes: */
     744        for (i=0;i<3;i++){
     745                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
     746                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
     747                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
     748                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
     749                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
     750                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
     751                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
     752                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
     753                *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=l1l2l3[i];
     754                *(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
     755                *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
     756                *(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=l1l2l3[i];
     757                *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=l1l2l3[i];
     758                *(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
     759                *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
     760                *(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=l1l2l3[i];
     761
     762        }
     763}
     764/*}}}*/
     765/*FUNCTION PentaRef::GetLprimeMacAyealStokes {{{1*/
     766void PentaRef::GetLprimeMacAyealStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){
     767
     768        /*
     769         * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
     770         * For grid i, Lpi can be expressed in the actual coordinate system
     771         * by:
     772         *       Lpi=[ h    0    0   0]
     773         *                     [ 0    h    0   0]
     774         *                     [ 0    0    h   0]
     775         *                     [ 0    0    h   0]
     776         *                     [ 0    0  dh/dz 0]
     777         *                     [ 0    0  dh/dz 0]
     778         *           [ 0    0    0   h]
     779         *           [ 0    0    0   h]
     780         * where h is the interpolation function for grid i.
     781         */
     782        int i;
     783        int num_dof=4;
     784
     785        double l1l2l3[NUMNODESP1_2d];
     786        double dh1dh6[3][NUMNODESP1];
     787
     788        /*Get l1l2l3 in actual coordinate system: */
     789        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     790        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     791        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     792
     793        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     794
     795        /*Build LprimeStokes: */
     796        for (i=0;i<3;i++){
     797                *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
     798                *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
     799                *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
     800                *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
     801                *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
     802                *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
     803                *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
     804                *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
     805                *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
     806                *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
     807                *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
     808                *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
     809                *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
     810                *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
     811                *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
     812                *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
     813                *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
     814                *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
     815                *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=dh1dh6[2][i];
     816                *(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
     817                *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
     818                *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
     819                *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=dh1dh6[2][i];
     820                *(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
     821                *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
     822                *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
     823                *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=0;
     824                *(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=l1l2l3[i];
     825                *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
     826                *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
     827                *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=0;
     828                *(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=l1l2l3[i];
     829        }
     830}
     831/*}}}*/
     832/*FUNCTION PentaRef::GetLStokesMacAyeal {{{1*/
     833void PentaRef::GetLStokesMacAyeal(double* LStokes, GaussPenta* gauss){
     834        /*
     835         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     836         * For grid i, Li can be expressed in the actual coordinate system
     837         * by:
     838         *       Li=[ h    0    0   0]
     839         *                    [ 0    h    0   0]
     840         *                    [ 0    0    h   0]
     841         *                    [ 0    0    h   0]
     842         * where h is the interpolation function for grid i.
     843         */
     844
     845        int i;
     846        int num_dof=4;
     847
     848        double l1l2l3[NUMNODESP1_2d];
     849
     850
     851        /*Get l1l2l3 in actual coordinate system: */
     852        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     853        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     854        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     855
     856        /*Build LStokes: */
     857        for (i=0;i<3;i++){
     858                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
     859                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
     860                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
     861                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
     862                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
     863                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
     864                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
     865                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
     866                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
     867                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
     868                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
     869                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
     870                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
     871                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
     872                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
     873                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
     874
     875        }
     876}
     877/*}}}*/
     878/*FUNCTION PentaRef::GetLprimeStokesMacAyeal {{{1*/
     879void PentaRef::GetLprimeStokesMacAyeal(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){
     880
     881        /*
     882         * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
     883         * For grid i, Lpi can be expressed in the actual coordinate system
     884         * by:
     885         *       Lpi=[ h    0 ]
     886         *                     [ 0    h ]
     887         *                     [ h    0 ]
     888         *                     [ 0    h ]
     889         * where h is the interpolation function for grid i.
     890         */
     891        int i;
     892        int num_dof=2;
     893
     894        double l1l2l3[NUMNODESP1_2d];
     895        double dh1dh6[3][NUMNODESP1];
     896
     897        /*Get l1l2l3 in actual coordinate system: */
     898        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     899        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     900        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     901
     902        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     903
     904        /*Build LprimeStokes: */
     905        for (i=0;i<3;i++){
     906                *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
     907                *(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
     908                *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
     909                *(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
     910                *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
     911                *(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
     912                *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
     913                *(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
     914        }
     915}
     916/*}}}*/
    715917/*FUNCTION PentaRef::GetJacobian {{{1*/
    716918void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){

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

    r6987 r7070  
    5151                void GetLStokes(double* LStokes, GaussPenta* gauss);
    5252                void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss);
     53                void GetLMacAyealStokes(double* LMacAyealStokes, GaussPenta* gauss);
     54                void GetLprimeMacAyealStokes(double* LprimeMacAyealStokes, double* xyz_list, GaussPenta* gauss);
     55                void GetLStokesMacAyeal(double* LStokesMacAyeal, GaussPenta* gauss);
     56                void GetLprimeStokesMacAyeal(double* LprimeStokesMacAyeal, double* xyz_list, GaussPenta* gauss);
    5357                void GetParameterValue(double* pvalue,double* plist, GaussPenta* gauss);
    5458                void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){_error_("only PentaGauss are supported");};
Note: See TracChangeset for help on using the changeset viewer.