Changeset 15779


Ignore:
Timestamp:
08/09/13 15:10:41 (12 years ago)
Author:
seroussi
Message:

BUG: fixed coupling SSA/FS sheet

Location:
issm/trunk-jpl/src/c/classes/Elements
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r15776 r15779  
    69246924        }
    69256925
    6926 
    69276926        /*Initialize Element matrix and return if necessary*/
    69286927        ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,    this->parameters,SSAApproximationEnum);
    … …  
    69886987/*FUNCTION Penta::CreateKMatrixCouplingSSAFSFriction {{{*/
    69896988ElementMatrix* Penta::CreateKMatrixCouplingSSAFSFriction(void){
    6990 
    69916989        /*Constants*/
    6992         const int numdof    = NUMVERTICES *NDOF4;
     6990        const int numdofs   = (NUMVERTICES+1)*NDOF3 + NUMVERTICES*NDOF1;
    69936991        const int numdofm   = NUMVERTICES *NDOF2;
    6994         const int numdof2d  = NUMVERTICES2D *NDOF4;
     6992        const int numdof2d  = NUMVERTICES2D *NDOF3;
    69956993        const int numdof2dm = NUMVERTICES2D *NDOF2;
    6996         const int numdoftot = numdof+numdofm;
     6994        const int numdoftot = NUMVERTICES*2 + (NUMVERTICES+1)*3 +NUMVERTICES; // HO + FS vel + FS Pressure
    69976995
    69986996        /*Intermediaries */
    … …  
    70067004        IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    70077005        IssmDouble LSSAFS[8][numdof2dm];
    7008         IssmDouble LprimeSSAFS[8][numdof2d];
     7006        IssmDouble LprimeSSAFS[8][numdofs];
    70097007        IssmDouble DLSSAFS[8][8]={0.0};
    70107008        IssmDouble LFSSSA[4][numdof2d];
    70117009        IssmDouble LprimeFSSSA[4][numdof2dm];
    70127010        IssmDouble DLFSSSA[4][4]={0.0};
    7013         IssmDouble Ke_drag_gaussian[numdof2dm][numdof2d];
     7011        IssmDouble Ke_drag_gaussian[numdof2dm][numdofs];
    70147012        IssmDouble Ke_drag_gaussian2[numdof2d][numdof2dm];
    70157013        Friction*  friction=NULL;
    70167014        GaussPenta *gauss=NULL;
    7017         Node *node_list[20];
     7015        Node       *node_list[20];
    70187016
    70197017        /*If on water or not FS, skip stiffness: */
    … …  
    70257023        int numnodes  = 2*vnumnodes-1+pnumnodes;
    70267024
    7027         int* cs_list = xNew<int>(numnodes);
    7028 
    7029         ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,SSAApproximationEnum);
    7030         ElementMatrix* Ke2=new ElementMatrix(this->nodes,numnodes,this->parameters,FSvelocityEnum);
     7025        /*Prepare node list*/
     7026        int* cs_list = xNew<int>(2*vnumnodes+pnumnodes);
     7027        for(i=0;i<vnumnodes-1;i++){
     7028                node_list[i] = this->nodes[i];
     7029                cs_list[i] = XYEnum;
     7030        }
     7031        for(i=0;i<vnumnodes;i++){
     7032                node_list[i+vnumnodes-1] = this->nodes[i];
     7033                cs_list[i+vnumnodes-1] = XYZEnum;
     7034        }
     7035        for(i=0;i<pnumnodes;i++){
     7036                node_list[2*vnumnodes-1+i] = this->nodes[vnumnodes+i];
     7037                cs_list[2*vnumnodes-1+i] = PressureEnum;
     7038        }
     7039
     7040        ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,        this->parameters,SSAApproximationEnum);
     7041        ElementMatrix* Ke2=new ElementMatrix(this->nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    70317042        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    70327043        delete Ke1; delete Ke2;
    7033 
    7034         /*Prepare node list*/
    7035         for(i=0;i<numnodes+NUMVERTICES;i++){
    7036                 node_list[i+0*NUMVERTICES] = this->nodes[i];
    7037                 node_list[i+1*NUMVERTICES] = this->nodes[i];
    7038                 cs_list[i+0*NUMVERTICES] = XYEnum;
    7039                 cs_list[i+1*NUMVERTICES] = XYZEnum;
    7040         }
    70417044
    70427045        /*Retrieve all inputs and parameters*/
    … …  
    70867089                TripleMultiply( &LSSAFS[0][0],8,numdof2dm,1,
    70877090                                        &DLSSAFS[0][0],8,8,0,
    7088                                         &LprimeSSAFS[0][0],8,numdof2d,0,
     7091                                        &LprimeSSAFS[0][0],8,numdofs,0,
    70897092                                        &Ke_drag_gaussian[0][0],0);
    70907093
    … …  
    70937096                                        &LprimeFSSSA[0][0],4,numdof2dm,0,
    70947097                                        &Ke_drag_gaussian2[0][0],0);
    7095 
    7096                 for(i=0;i<numdof2dm;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdoftot+j+numdofm]+=Ke_drag_gaussian[i][j];
     7098                for(i=0;i<numdof2dm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftot+j+numdofm]+=Ke_drag_gaussian[i][j];
    70977099                for(i=0;i<numdof2d;i++) for(j=0;j<numdof2dm;j++) Ke->values[(i+numdofm)*numdoftot+j]+=Ke_drag_gaussian2[i][j];
    70987100        }
    … …  
    71187120
    71197121        /*Intermediaries*/
    7120         int       i,j,init;
     7122        int   i,j,init;
    71217123        Node  *node_list[NUMVERTICES*3+1];
    71227124        int   cs_list[NUMVERTICES*3+1];
    … …  
    82208222        inputs->GetInputValue(&approximation,ApproximationEnum);
    82218223        if(approximation!=SSAFSApproximationEnum) return NULL;
    8222         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSvelocityEnum);
     8224        int vnumnodes = this->NumberofNodesVelocity();
     8225        int pnumnodes = this->NumberofNodesPressure();
     8226
     8227        /*Prepare coordinate system list*/
     8228        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     8229        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     8230        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     8231        ElementVector* pe=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    82238232
    82248233        /*Retrieve all inputs and parameters*/
    … …  
    82268235        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    82278236        this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    8228         Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    8229         Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    8230         Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
    8231         Input* vzSSA_input=inputs->GetInput(VzSSAEnum);   _assert_(vzSSA_input);
     8237        Input* vx_input=inputs->GetInput(VxEnum);        _assert_(vx_input);
     8238        Input* vy_input=inputs->GetInput(VyEnum);        _assert_(vy_input);
     8239        Input* vz_input=inputs->GetInput(VzEnum);        _assert_(vz_input);
     8240        Input* vzSSA_input=inputs->GetInput(VzSSAEnum);  _assert_(vzSSA_input);
    82328241
    82338242        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    … …  
    82548263
    82558264                for(i=0;i<NUMVERTICES2D;i++){
    8256                         pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
    8257                         pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
    8258                         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])*basis[i];
     8265                        pe->values[i*NDOF3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
     8266                        pe->values[i*NDOF3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
     8267                        pe->values[i*NDOF3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
    82598268                }
    82608269        }
    82618270
    82628271        /*Transform coordinate system*/
    8263         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum);
     8272        TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
    82648273
    82658274        /*Clean up and return*/

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

    r15757 r15779  
    981981         * where h is the interpolation function for node i.
    982982         */
    983         int num_dof=4;
     983        int num_dof=3;
     984        int num_dof_vel=3*NUMNODESP1b;
     985        int num_dof_total=3*NUMNODESP1b+1*NUMNODESP1;
    984986        IssmDouble L1L2l3[NUMNODESP1_2d];
    985987        IssmDouble dbasis[3][NUMNODESP1];
    … …  
    994996        /*Build LprimeFS: */
    995997        for(int i=0;i<3;i++){
    996                 LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
    997                 LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
    998                 LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
    999                 LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
    1000                 LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
    1001                 LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
    1002                 LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
    1003                 LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
    1004                 LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
    1005                 LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
    1006                 LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = L1L2l3[i];
    1007                 LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
    1008                 LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
    1009                 LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
    1010                 LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = L1L2l3[i];
    1011                 LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
    1012                 LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.;
    1013                 LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.;
    1014                 LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = dbasis[2][i];
    1015                 LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.;
    1016                 LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.;
    1017                 LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.;
    1018                 LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = dbasis[2][i];
    1019                 LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.;
    1020                 LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.;
    1021                 LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.;
    1022                 LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = 0.;
    1023                 LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = L1L2l3[i];
    1024                 LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.;
    1025                 LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.;
    1026                 LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = 0.;
    1027                 LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = L1L2l3[i];
     998                LprimeFS[num_dof_total*0+num_dof*i+0] = L1L2l3[i];
     999                LprimeFS[num_dof_total*0+num_dof*i+1] = 0.;
     1000                LprimeFS[num_dof_total*0+num_dof*i+2] = 0.;
     1001                LprimeFS[num_dof_total*1+num_dof*i+0] = 0.;
     1002                LprimeFS[num_dof_total*1+num_dof*i+1] = L1L2l3[i];
     1003                LprimeFS[num_dof_total*1+num_dof*i+2] = 0.;
     1004                LprimeFS[num_dof_total*2+num_dof*i+0] = 0.;
     1005                LprimeFS[num_dof_total*2+num_dof*i+1] = 0.;
     1006                LprimeFS[num_dof_total*2+num_dof*i+2] = L1L2l3[i];
     1007                LprimeFS[num_dof_total*3+num_dof*i+0] = 0.;
     1008                LprimeFS[num_dof_total*3+num_dof*i+1] = 0.;
     1009                LprimeFS[num_dof_total*3+num_dof*i+2] = L1L2l3[i];
     1010                LprimeFS[num_dof_total*4+num_dof*i+0] = 0.;
     1011                LprimeFS[num_dof_total*4+num_dof*i+1] = 0.;
     1012                LprimeFS[num_dof_total*4+num_dof*i+2] = dbasis[2][i];
     1013                LprimeFS[num_dof_total*5+num_dof*i+0] = 0.;
     1014                LprimeFS[num_dof_total*5+num_dof*i+1] = 0.;
     1015                LprimeFS[num_dof_total*5+num_dof*i+2] = dbasis[2][i];
     1016                LprimeFS[num_dof_total*6+num_dof*i+0] = 0.;
     1017                LprimeFS[num_dof_total*6+num_dof*i+1] = 0.;
     1018                LprimeFS[num_dof_total*6+num_dof*i+2] = 0.;
     1019                LprimeFS[num_dof_total*7+num_dof*i+0] = 0.;
     1020                LprimeFS[num_dof_total*7+num_dof*i+1] = 0.;
     1021                LprimeFS[num_dof_total*7+num_dof*i+2] = 0.;
     1022        }
     1023        for(int i=3;i<7;i++){
     1024                LprimeFS[num_dof_total*0+num_dof*i+0] = 0.;
     1025                LprimeFS[num_dof_total*0+num_dof*i+1] = 0.;
     1026                LprimeFS[num_dof_total*0+num_dof*i+2] = 0.;
     1027                LprimeFS[num_dof_total*1+num_dof*i+0] = 0.;
     1028                LprimeFS[num_dof_total*1+num_dof*i+1] = 0.;
     1029                LprimeFS[num_dof_total*1+num_dof*i+2] = 0.;
     1030                LprimeFS[num_dof_total*2+num_dof*i+0] = 0.;
     1031                LprimeFS[num_dof_total*2+num_dof*i+1] = 0.;
     1032                LprimeFS[num_dof_total*2+num_dof*i+2] = 0.;
     1033                LprimeFS[num_dof_total*3+num_dof*i+0] = 0.;
     1034                LprimeFS[num_dof_total*3+num_dof*i+1] = 0.;
     1035                LprimeFS[num_dof_total*3+num_dof*i+2] = 0.;
     1036                LprimeFS[num_dof_total*4+num_dof*i+0] = 0.;
     1037                LprimeFS[num_dof_total*4+num_dof*i+1] = 0.;
     1038                LprimeFS[num_dof_total*4+num_dof*i+2] = 0.;
     1039                LprimeFS[num_dof_total*5+num_dof*i+0] = 0.;
     1040                LprimeFS[num_dof_total*5+num_dof*i+1] = 0.;
     1041                LprimeFS[num_dof_total*5+num_dof*i+2] = 0.;
     1042                LprimeFS[num_dof_total*6+num_dof*i+0] = 0.;
     1043                LprimeFS[num_dof_total*6+num_dof*i+1] = 0.;
     1044                LprimeFS[num_dof_total*6+num_dof*i+2] = 0.;
     1045                LprimeFS[num_dof_total*7+num_dof*i+0] = 0.;
     1046                LprimeFS[num_dof_total*7+num_dof*i+1] = 0.;
     1047                LprimeFS[num_dof_total*7+num_dof*i+2] = 0.;
     1048        }
     1049        for(int i=0;i<3;i++){
     1050                LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.;
     1051                LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.;
     1052                LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.;
     1053                LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.;
     1054                LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.;
     1055                LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.;
     1056                LprimeFS[num_dof_total*6+num_dof_vel+i] = L1L2l3[i];
     1057                LprimeFS[num_dof_total*7+num_dof_vel+i] = L1L2l3[i];
     1058        }
     1059        for(int i=3;i<6;i++){
     1060                LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.;
     1061                LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.;
     1062                LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.;
     1063                LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.;
     1064                LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.;
     1065                LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.;
     1066                LprimeFS[num_dof_total*6+num_dof_vel+i] = 0.;
     1067                LprimeFS[num_dof_total*7+num_dof_vel+i] = 0.;
    10281068        }
    10291069}
    … …  
    10341074         * For node i, Li can be expressed in the actual coordinate system
    10351075         * by:
    1036          *       Li=[ h    0    0   0]
    1037          *                    [ 0    h    0   0]
    1038          *                    [ 0    0    h   0]
    1039          *                    [ 0    0    h   0]
     1076         *       Li=[ h    0    0 ]
     1077         *                    [ 0    h    0 ]
     1078         *                    [ 0    0    h ]
     1079         *                    [ 0    0    h ]
    10401080         * where h is the interpolation function for node i.
    10411081         */
    10421082
    1043         int num_dof=4;
     1083        int num_dof=3;
    10441084        IssmDouble L1L2l3[NUMNODESP1_2d];
    10451085
    … …  
    10541094                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
    10551095                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
    1056                 LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
    10571096                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
    10581097                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
    10591098                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
    1060                 LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
    10611099                LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
    10621100                LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
    10631101                LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = L1L2l3[i];
    1064                 LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
    10651102                LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
    10661103                LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
    10671104                LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = L1L2l3[i];
    1068                 LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
    10691105        }
    10701106}
Note: See TracChangeset for help on using the changeset viewer.