Changeset 10480


Ignore:
Timestamp:
11/04/11 16:41:35 (13 years ago)
Author:
Mathieu Morlighem
Message:

Working on better friction for FS

Location:
issm/trunk/src/c
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r10444 r10480  
    59865986        ElementMatrix* Ke1=CreateKMatrixDiagnosticStokesViscous();
    59875987        ElementMatrix* Ke2=CreateKMatrixDiagnosticStokesFriction();
     5988        //ElementMatrix* Ke2=CreateKMatrixDiagnosticStokesFriction2();
    59885989        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    59895990
    … …  
    61366137
    61376138        /*Transform Coordinate System*/
     6139        if(id==1) printarray(&bed_normal[0],1,3);
     6140        if(id==1) printarray(&nodes[0]->coord_system[0][0],3,3);
     6141        if(id==1) printf("Before rotation\n");
     6142        if(id==1) printarray(Ke->values,4,Ke->ncols);
    61386143        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4);
     6144        if(id==1) printf("After rotation\n");
     6145        if(id==1) printarray(Ke->values,4,Ke->ncols);
     6146        if(id==1) printf("New:\n");
     6147        if(id==1) Ke=CreateKMatrixDiagnosticStokesFriction2();
     6148
     6149        //if(id==56) printarray(Ke->values,1,9);
     6150        //if(id==56) printarray(&nodes[0]->coord_system[0][0],3,3);
     6151        //if(id==56) _error_("STOP");
    61396152
    61406153        /*Clean up and return*/
    61416154        delete gauss;
    61426155        delete friction;
     6156        return Ke;
     6157}
     6158/*}}}*/
     6159/*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction2{{{1*/
     6160ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction2(void){
     6161
     6162        /*Constants*/
     6163        const int numdof=NUMVERTICES*NDOF4;
     6164        const int numdof2d=NUMVERTICES2D*NDOF4;
     6165
     6166        /*Intermediaries */
     6167        int        i,j,ig;
     6168        int        analysis_type,approximation;
     6169        double     alpha2,Jdet2d;
     6170        double     stokesreconditioning,viscosity;
     6171        double     epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     6172        double     xyz_list[NUMVERTICES][3];
     6173        double    xyz_list_tria[NUMVERTICES2D][3];
     6174        double     LStokes2[2][numdof2d];
     6175        //double     LprimeStokes2[4][numdof2d];
     6176        double     DLStokes[2][2]={0.0};
     6177        double     Ke_drag_gaussian[numdof2d][numdof2d];
     6178        Friction*  friction=NULL;
     6179        GaussPenta *gauss=NULL;
     6180
     6181        /*If on water or not Stokes, skip stiffness: */
     6182        inputs->GetInputValue(&approximation,ApproximationEnum);
     6183        if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum)) return NULL;
     6184        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     6185
     6186        /*Retrieve all inputs and parameters*/
     6187        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     6188        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     6189        parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     6190        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     6191        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     6192        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     6193        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     6194
     6195        /*build friction object, used later on: */
     6196        friction=new Friction("3d",inputs,matpar,analysis_type);
     6197
     6198        /* Start  looping on the number of gaussian points: */
     6199        gauss=new GaussPenta(0,1,2,2);
     6200        for (ig=gauss->begin();ig<gauss->end();ig++){
     6201
     6202                gauss->GaussPoint(ig);
     6203
     6204                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
     6205                GetLStokes3(&LStokes2[0][0], gauss);
     6206                //GetLprimeStokes2(&LprimeStokes2[0][0], &xyz_list[0][0], gauss);
     6207
     6208                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     6209                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     6210
     6211                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
     6212
     6213                DLStokes[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
     6214                DLStokes[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy
     6215                //DLStokes[2][2] = +2*viscosity*gauss->weight*Jdet2d;
     6216                //DLStokes[3][3] = -stokesreconditioning*gauss->weight*Jdet2d;
     6217
     6218                //TripleMultiply( &LStokes2[0][0],4,numdof2d,1,
     6219                //                      &DLStokes[0][0],4,4,0,
     6220                //                      &LprimeStokes2[0][0],4,numdof2d,0,
     6221                //                      &Ke_drag_gaussian[0][0],0);
     6222                TripleMultiply( &LStokes2[0][0],2,numdof2d,1,
     6223                                        &DLStokes[0][0],2,2,0,
     6224                                        &LStokes2[0][0],2,numdof2d,0,
     6225                                        &Ke_drag_gaussian[0][0],0);
     6226
     6227                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof+j]+=Ke_drag_gaussian[i][j];
     6228        }
     6229
     6230        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
     6231        //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4);
     6232       
     6233        //printf("New:\n");
     6234        //printarray(&Ke_drag_gaussian[0][0],numdof2d,numdof2d);
     6235        //printf("Old:\n");
     6236        //Ke=CreateKMatrixDiagnosticStokesFriction();
     6237        //_error_("STOP");
     6238        //if(id==1) printarray(Ke->values,4,Ke->ncols);
     6239        //if(id==1) printarray(Ke->values,1,9);
     6240
     6241        /*Clean up and return*/
     6242        delete gauss;
     6243        delete friction;
     6244
     6245        //return CreateKMatrixDiagnosticStokesFriction();
    61436246        return Ke;
    61446247}

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

    r10440 r10480  
    223223                ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void);
    224224                ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void);
     225                ElementMatrix* CreateKMatrixDiagnosticStokesFriction2(void);
    225226                ElementMatrix* CreateKMatrixDiagnosticVert(void);
    226227                ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);

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

    r10135 r10480  
    651651}
    652652/*}}}*/
     653/*FUNCTION PentaRef::GetLStokes2{{{1*/
     654void PentaRef::GetLStokes2(double* LStokes, GaussPenta* gauss){
     655        /*
     656         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     657         * For node i, Li can be expressed in the actual coordinate system
     658         * by:
     659         *       Li=[ h 0 0 0]
     660         *                    [ 0 h 0 0]
     661         *                    [ 0 0 h h]
     662         *                    [ 0 0 0 0]
     663         * where h is the interpolation function for node i.
     664         */
     665
     666        int i;
     667        int num_dof=4;
     668
     669        double l1l2l3[NUMNODESP1_2d];
     670
     671
     672        /*Get l1l2l3 in actual coordinate system: */
     673        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     674        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     675        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     676
     677        /*Build LStokes: */
     678        for (i=0;i<3;i++){
     679                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
     680                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.;
     681                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.;
     682                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.;
     683
     684                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.;
     685                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
     686                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.;
     687                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.;
     688
     689                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+0)=0.;
     690                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0.;
     691                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
     692                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0.;
     693
     694                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+0)=0.;
     695                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0.;
     696                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
     697                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0.;
     698        }
     699}
     700/*}}}*/
     701/*FUNCTION PentaRef::GetLStokes3{{{1*/
     702void PentaRef::GetLStokes3(double* LStokes, GaussPenta* gauss){
     703        /*
     704         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     705         * For node i, Li can be expressed in the actual coordinate system
     706         * by:
     707         *       Li=[ h 0 ]
     708         *                    [ 0 h ]
     709         *                    [ 0 0 ]
     710         *                    [ 0 0 ]
     711         * where h is the interpolation function for node i.
     712         */
     713
     714        int i;
     715        int num_dof=4;
     716
     717        double l1l2l3[NUMNODESP1_2d];
     718
     719
     720        /*Get l1l2l3 in actual coordinate system: */
     721        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     722        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     723        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     724
     725        /*Build LStokes: */
     726        for (i=0;i<3;i++){
     727                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
     728                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.;
     729                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.;
     730                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.;
     731
     732                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.;
     733                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
     734                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.;
     735                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.;
     736        }
     737}
     738/*}}}*/
     739/*FUNCTION PentaRef::GetLprimeStokes2{{{1*/
     740void PentaRef::GetLprimeStokes2(double* LStokes,double* xyz_list,GaussPenta* gauss){
     741        /*
     742         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     743         * For node i, Li can be expressed in the actual coordinate system
     744         * by:
     745         *       Li=[ h 0 0 0]
     746         *                    [ 0 h 0 0]
     747         *                    [ 0 0 dh\dz 0]
     748         *                    [ 0 0 0 h]
     749         * where h is the interpolation function for node i.
     750         */
     751
     752        int i;
     753        int num_dof=4;
     754
     755        double l1l2l3[NUMNODESP1_2d];
     756        double dh1dh6[3][NUMNODESP1];
     757
     758
     759        /*Get l1l2l3 in actual coordinate system: */
     760        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     761        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     762        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     763
     764        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     765
     766        /*Build LStokes: */
     767        for (i=0;i<3;i++){
     768                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
     769                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.;
     770                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.;
     771                *(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.;
     772
     773                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.;
     774                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
     775                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.;
     776                *(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.;
     777
     778                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+0)=0.;
     779                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0.;
     780                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=dh1dh6[2][i];
     781                *(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0.;
     782
     783                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+0)=0.;
     784                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0.;
     785                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0.;
     786                *(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=l1l2l3[i];
     787        }
     788}
     789/*}}}*/
    653790/*FUNCTION PentaRef::GetLprimeStokes {{{1*/
    654791void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){
    … …  
    658795         * For node i, Lpi can be expressed in the actual coordinate system
    659796         * by:
    660          *       Lpi=[ h    0    0   0]
    661          *                     [ 0    h    0   0]
    662          *                     [ h    0    0   0]
    663          *                     [ 0    h    0   0]
    664          *                     [ 0    0    h   0]
    665          *                     [ 0    0    h   0]
    666          *                     [ 0    0  dh/dz 0]
    667          *                     [ 0    0  dh/dz 0]
    668          *                     [ 0    0  dh/dz 0]
    669          *                     [dh/dz 0  dh/dx 0]
    670          *                     [ 0 dh/dz dh/dy 0]
    671          *           [ 0    0    0   h]
    672          *           [ 0    0    0   h]
    673          *           [ 0    0    0   h]
     797         *       Lpi=[ h    0    0   0]1
     798         *                     [ 0    h    0   0]2
     799         *                     [ h    0    0   0]3
     800         *                     [ 0    h    0   0]4
     801         *                     [ 0    0    h   0]5
     802         *                     [ 0    0    h   0]6
     803         *                     [ 0    0  dh/dz 0]7
     804         *                     [ 0    0  dh/dz 0]8
     805         *                     [ 0    0  dh/dz 0]9
     806         *                     [dh/dz 0  dh/dx 0]0
     807         *                     [ 0 dh/dz dh/dy 0]1
     808         *           [ 0    0    0   h]2
     809         *           [ 0    0    0   h]3
     810         *           [ 0    0    0   h]4
     811         *
     812         *       Li=[ h    0    0   0]1
     813         *                    [ 0    h    0   0]2
     814         *                    [ 0    0    h   0]3
     815         *                    [ 0    0    h   0]4
     816         *                    [ h    0    0   0]5
     817         *                    [ 0    h    0   0]6
     818         *                    [ h    0    0   0]7
     819         *                    [ 0    h    0   0]8
     820         *                    [ 0    0    h   0]9
     821         *                    [ 0    0    h   0]0
     822         *                    [ 0    0    h   0]1
     823         *                    [ h    0    0   0]2
     824         *                    [ 0    h    0   0]3
     825         *                    [ 0    0    h   0]4
    674826         * where h is the interpolation function for node i.
    675827         */

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

    r10135 r10480  
    5151                void GetL(double* L, GaussPenta* gauss,int numdof);
    5252                void GetLStokes(double* LStokes, GaussPenta* gauss);
     53                void GetLStokes2(double* LStokes, GaussPenta* gauss);
     54                void GetLStokes3(double* LStokes, GaussPenta* gauss);
     55                void GetLprimeStokes2(double* LStokes,double* xyz_list,GaussPenta* gauss);
    5356                void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss);
    5457                void GetLMacAyealStokes(double* LMacAyealStokes, GaussPenta* gauss);

  • issm/trunk/src/c/shared/Elements/elements.h

    r10407 r10480  
    2828        for(int i=0;i<lines;i++){ 
    2929                printf("   [ ");
    30                 for(int j=0;j<cols;j++) printf(" %7.5g ",array[i*cols+j]);
     30                for(int j=0;j<cols;j++) printf(" %12.7g ",array[i*cols+j]);
    3131                printf(" ]\n");
    3232        } 
Note: See TracChangeset for help on using the changeset viewer.