Changeset 15757


Ignore:
Timestamp:
08/08/13 15:53:27 (12 years ago)
Author:
seroussi
Message:

BUG: coupling SSA/FS now running with errors

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

Legend:

Unmodified
Added
Removed

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

    r15753 r15757  
    68766876
    68776877        /*Constants*/
    6878         const int numnodes    = 2 *NUMVERTICES;
    6879         const int numdofm     = NDOF2 *NUMVERTICES2D;
    6880         const int numdofs     = NDOF4 *NUMVERTICES;
    6881         const int numdoftotal = 2 *numdofm+numdofs;
     6878        const int numdofm      = NDOF2 *NUMVERTICES2D;
     6879        const int numdofs      = NDOF4 *NUMVERTICES;
     6880        const int numdofstotal = NDOF4 *NUMVERTICES + NDOF3;
     6881        const int numdoftotal  = 2 *numdofm+numdofstotal;
    68826882
    68836883        /*Intermediaries */
    6884         int         i,j;
     6884        int        i,j;
    68856885        IssmDouble Jdet;
    68866886        IssmDouble viscosity,FSreconditioning; //viscosity
    68876887        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    68886888        IssmDouble xyz_list[NUMVERTICES][3];
    6889         IssmDouble B[4][numdofs+3];
     6889        IssmDouble B[4][numdofstotal];
    68906890        IssmDouble Bprime[4][numdofm];
    68916891        IssmDouble B2[3][numdofm];
    6892         IssmDouble Bprime2[3][numdofs+3];
     6892        IssmDouble Bprime2[3][numdofstotal];
    68936893        IssmDouble D[4][4]={0.0};            // material matrix, simple scalar matrix.
    68946894        IssmDouble D2[3][3]={0.0};            // material matrix, simple scalar matrix.
    … …  
    68966896        IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix
    68976897        IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix
    6898         IssmDouble Ke_gg_gaussian[numdofs+3][numdofm]; //stiffness matrix evaluated at the gaussian point.
    6899         IssmDouble Ke_gg_gaussian2[numdofm][numdofs+3]; //stiffness matrix evaluated at the gaussian point.
     6898        IssmDouble Ke_gg_gaussian[numdofstotal][numdofm]; //stiffness matrix evaluated at the gaussian point.
     6899        IssmDouble Ke_gg_gaussian2[numdofm][numdofstotal]; //stiffness matrix evaluated at the gaussian point.
    69006900        GaussPenta *gauss=NULL;
    69016901        GaussTria  *gauss_tria=NULL;
    6902         Node       *node_list[numnodes];
    6903         int         cs_list[numnodes];
     6902        Node       *node_list[20];
    69046903
    69056904        /*Find penta on bed as FS must be coupled to the dofs on the bed: */
    … …  
    69076906        Tria* tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
    69086907
     6908        int vnumnodes = this->NumberofNodesVelocity();
     6909        int pnumnodes = this->NumberofNodesPressure();
     6910        int numnodes  = 2*vnumnodes-1+pnumnodes;
     6911
    69096912        /*Prepare node list*/
    6910         for(i=0;i<NUMVERTICES;i++){
    6911                 node_list[i+0*NUMVERTICES] = pentabase->nodes[i];
    6912                 node_list[i+1*NUMVERTICES] = this->nodes[i];
    6913                 cs_list[i+0*NUMVERTICES] = XYEnum;
    6914                 cs_list[i+1*NUMVERTICES] = XYZEnum;
    6915         }
     6913        int* cs_list = xNew<int>(2*vnumnodes+pnumnodes);
     6914        for(i=0;i<vnumnodes-1;i++){
     6915                node_list[i] = pentabase->nodes[i];
     6916                cs_list[i] = XYEnum;
     6917        }
     6918        for(i=0;i<vnumnodes;i++){
     6919                node_list[i+vnumnodes-1] = this->nodes[i];
     6920                cs_list[i+vnumnodes-1] = XYZEnum;
     6921        }
     6922        for(i=0;i<pnumnodes;i++){
     6923                node_list[2*vnumnodes-1+i] = this->nodes[vnumnodes+i];
     6924                cs_list[2*vnumnodes-1+i] = PressureEnum;
     6925        }
     6926
    69166927
    69176928        /*Initialize Element matrix and return if necessary*/
    6918         ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,SSAApproximationEnum);
    6919         ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,FSvelocityEnum);
    6920         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
     6929        ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,    this->parameters,SSAApproximationEnum);
     6930        ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,2*NUMVERTICES+1,this->parameters,FSvelocityEnum);
     6931        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    69216932        delete Ke1; delete Ke2;
    69226933
    … …  
    69506961                for (i=0;i<3;i++) D2[i][i]=D_scalar;
    69516962
    6952                 TripleMultiply( &B[0][0],4,numdofs+3,1,
     6963                TripleMultiply( &B[0][0],4,numdofstotal,1,
    69536964                                        &D[0][0],4,4,0,
    69546965                                        &Bprime[0][0],4,numdofm,0,
    … …  
    69576968                TripleMultiply( &B2[0][0],3,numdofm,1,
    69586969                                        &D2[0][0],3,3,0,
    6959                                         &Bprime2[0][0],3,numdofs+3,0,
     6970                                        &Bprime2[0][0],3,numdofstotal,0,
    69606971                                        &Ke_gg_gaussian2[0][0],0);
    69616972
    6962                 for( i=0; i<numdofs; i++) for(j=0;j<numdofm; j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    6963                 for( i=0; i<numdofm; i++) for(j=0;j<numdofs; j++) Ke_gg2[i][j]+=Ke_gg_gaussian2[i][j];
     6973                for( i=0; i<numdofstotal; i++) for(j=0;j<numdofm; j++)      Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     6974                for( i=0; i<numdofm; i++)      for(j=0;j<numdofstotal; j++) Ke_gg2[i][j]+=Ke_gg_gaussian2[i][j];
    69646975        }
    69656976        for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j];
    … …  
    71117122        ElementMatrix* Ke2=new ElementMatrix(this->nodes,2*NUMVERTICES+1,this->parameters,FSvelocityEnum);
    71127123        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    7113         delete Ke1;
    7114         delete Ke2;
     7124        delete Ke1; delete Ke2;
    71157125       
    71167126        /*Compute HO Matrix with P1 element type\n");*/
    … …  
    71327142        }
    71337143
    7134         /*Transform Coordinate System*/ //Do not transform, already sone in the matrices
     7144        /*Transform Coordinate System*/ //Do not transform, already done in the matrices
    71357145        //TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
    71367146
    … …  
    81168126        inputs->GetInputValue(&approximation,ApproximationEnum);
    81178127        if(approximation!=SSAFSApproximationEnum) return NULL;
    8118         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSvelocityEnum);
     8128        int vnumnodes = this->NumberofNodesVelocity();
     8129        int pnumnodes = this->NumberofNodesPressure();
     8130
     8131        /*Prepare coordinate system list*/
     8132        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     8133        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     8134        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     8135        ElementVector* pe=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    81198136
    81208137        /*Retrieve all inputs and parameters*/
    … …  
    81428159
    81438160                for(i=0;i<NUMVERTICES;i++){
    8144                         pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
    8145                         pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
    8146                         pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
    8147                         pe->values[i*NDOF4+3]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
     8161                        pe->values[i*NDOF3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
     8162                        pe->values[i*NDOF3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
     8163                        pe->values[i*NDOF3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
     8164                        pe->values[NDOF3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
    81488165                }
    81498166        }
    81508167
    81518168        /*Transform coordinate system*/
    8152         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum);
     8169        TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
    81538170
    81548171        /*Clean up and return*/
     8172        xDelete<int>(cs_list);
    81558173        delete gauss;
    81568174        return pe;

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

    r15741 r15757  
    107107
    108108        /*Build B: */
    109         for(i=0;i<NUMNODESP1b;i++){
    110                 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = dbasismini[0][i];
    111                 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = 0.;
    112                 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
    113                 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = 0.;
    114                 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = dbasismini[1][i];
    115                 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
    116                 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = 0.5*dbasismini[1][i];
    117                 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = 0.5*dbasismini[0][i];
    118                 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = 0.;
    119                 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+0] = 0.;
    120                 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1] = 0.;
    121                 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2] = 0.;
     109        for(i=0;i<NUMNODESP1;i++){
     110                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+0] = dbasismini[0][i];
     111                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+1] = 0.;
     112                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+2] = 0.;
     113                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+0] = 0.;
     114                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+1] = dbasismini[1][i];
     115                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+2] = 0.;
     116                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+0] = 0.5*dbasismini[1][i];
     117                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+1] = 0.5*dbasismini[0][i];
     118                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+2] = 0.;
     119                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*i+0] = 0.;
     120                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*i+1] = 0.;
     121                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*i+2] = 0.;
     122        }
     123        for(i=0;i<1;i++){
     124                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+0] = 0.;
     125                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+1] = 0.;
     126                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+2] = 0.;
     127                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+0] = 0.;
     128                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+1] = 0.;
     129                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+2] = 0.;
     130                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+0] = 0.;
     131                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+1] = 0.;
     132                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+2] = 0.;
     133                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*(NUMNODESP1+i)+0] = 0.;
     134                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*(NUMNODESP1+i)+1] = 0.;
     135                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*(NUMNODESP1+i)+2] = 0.;
    122136        }
    123137
    124138        for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
    125                 B[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0;
    126                 B[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0;
    127                 B[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0;
    128                 B[(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3] = basis[i];
     139                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NUMNODESP1b*NDOF3+i] = 0;
     140                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NUMNODESP1b*NDOF3+i] = 0;
     141                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NUMNODESP1b*NDOF3+i] = 0;
     142                B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NUMNODESP1b*NDOF3+i] = basis[i];
    129143        }
    130144}
    … …  
    230244
    231245        /*Build Bprime: */
    232         for(i=0;i<NUMNODESP1b;i++){
    233                 Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+0] = 2.*dbasismini[0][i];
    234                 Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1] = dbasismini[1][i];
    235                 Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2] = 0.;
    236                 Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+0] = dbasismini[0][i];
    237                 Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1] = 2.*dbasismini[1][i];
    238                 Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2] = 0.;
    239                 Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+0] = dbasismini[1][i];
    240                 Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1] = dbasismini[0][i];
    241                 Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2] = 0.;
     246        for(i=0;i<NUMNODESP1;i++){
     247                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+0] = 2.*dbasismini[0][i];
     248                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+1] = dbasismini[1][i];
     249                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+2] = 0.;
     250                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+0] = dbasismini[0][i];
     251                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+1] = 2.*dbasismini[1][i];
     252                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+2] = 0.;
     253                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+0] = dbasismini[1][i];
     254                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+1] = dbasismini[0][i];
     255                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+2] = 0.;
     256        }
     257
     258        for(i=0;i<1;i++){ //Add zeros for the bubble function
     259                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+0] = 0.;
     260                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+1] = 0.;
     261                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+2] = 0.;
     262                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+0] = 0.;
     263                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+1] = 0.;
     264                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+2] = 0.;
     265                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+0] = 0.;
     266                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+1] = 0.;
     267                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+2] = 0.;
    242268        }
    243269
    244270        for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
    245                 Bprime[(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3] = 0.;
    246                 Bprime[(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3] = 0.;
    247                 Bprime[(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3] = 0.;
     271                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NUMNODESP1b*NDOF3+i] = 0.;
     272                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NUMNODESP1b*NDOF3+i] = 0.;
     273                Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NUMNODESP1b*NDOF3+i] = 0.;
    248274        }
    249275
Note: See TracChangeset for help on using the changeset viewer.