Changeset 7000


Ignore:
Timestamp:
01/07/11 16:34:33 (14 years ago)
Author:
seroussi
Message:

finished coupling viscous matrix MacAyeal/Stokes

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

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r6988 r7000  
    719719        const int    numdofm=NDOF2*NUMVERTICES2D;
    720720        const int    numdofs=NDOF4*NUMVERTICES;
    721         const int    numdoftotal=numdofm+numdofs;
     721        const int    numdoftotal=2*numdofm+numdofs;
    722722
    723723        /*Intermediaries */
    724724        int         i,j,ig;
    725725        double      Jdet;
    726         double      viscosity; //viscosity
    727         double      epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     726        double      viscosity,stokesreconditioning; //viscosity
     727        double      epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    728728        double      xyz_list[NUMVERTICES][3];
    729         double      B[4][numdofs];
     729        double      B[4][numdofs+3];
    730730        double      Bprime[4][numdofm];
    731731        double      B2[3][numdofm];
    732         double      Bprime2[3][numdofs];
     732        double      Bprime2[3][numdofs+3];
    733733        double      D[4][4]={0.0};            // material matrix, simple scalar matrix.
    734734        double      D2[3][3]={0.0};            // material matrix, simple scalar matrix.
     
    736736        double      Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix
    737737        double      Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix
    738         double      Ke_gg_gaussian[numdofs][numdofm]; //stiffness matrix evaluated at the gaussian point.
    739         double      Ke_gg_gaussian2[numdofm][numdofs]; //stiffness matrix evaluated at the gaussian point.
     738        double      Ke_gg_gaussian[numdofs+3][numdofm]; //stiffness matrix evaluated at the gaussian point.
     739        double      Ke_gg_gaussian2[numdofm][numdofs+3]; //stiffness matrix evaluated at the gaussian point.
    740740        GaussPenta *gauss=NULL;
    741741        GaussTria  *gauss_tria=NULL;
     
    754754        /* Get node coordinates and dof list: */
    755755        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     756        parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    756757        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    757758        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
     
    777778                D_scalar=2*viscosity*gauss->weight*Jdet;
    778779                for (i=0;i<3;i++) D[i][i]=D_scalar;
    779                 D[4][4]=(-1)*gauss->weight*Jdet;
     780                D[3][3]=-gauss->weight*Jdet*stokesreconditioning;
    780781                for (i=0;i<3;i++) D2[i][i]=D_scalar;
    781782
    782                 TripleMultiply( &B[0][0],4,numdofs,1,
     783                TripleMultiply( &B[0][0],4,numdofs+3,1,
    783784                                        &D[0][0],4,4,0,
    784785                                        &Bprime[0][0],4,numdofm,0,
     
    787788                TripleMultiply( &B2[0][0],3,numdofm,1,
    788789                                        &D2[0][0],3,3,0,
    789                                         &Bprime2[0][0],3,numdofs,0,
     790                                        &Bprime2[0][0],3,numdofs+3,0,
    790791                                        &Ke_gg_gaussian2[0][0],0);
    791792
     
    794795        }
    795796        for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j];
    796         for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[j][i];
    797         _error_("coupling not finished yet");
     797        for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j];
    798798
    799799        /*Clean-up and return*/
     
    990990        double      viscosity, oldviscosity, newviscosity, viscosity_overshoot;
    991991        double      epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     992        double      epsilons[6]; //6 for stokes
    992993        double      xyz_list[NUMVERTICES][3];
    993994        double      B[3][numdof2d];
     
    10401041                }
    10411042                else if (approximation==MacAyealStokesApproximationEnum){
    1042                         this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    1043                         matice->GetViscosity3dStokes(&newviscosity,&epsilon[0]);
     1043                        this->GetStrainRate3d(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     1044                        matice->GetViscosity3dStokes(&newviscosity,&epsilons[0]);
    10441045                }
    10451046                else _error_("approximation %i (%s) not supported yet",approximation,EnumToString(approximation));
  • issm/trunk/src/c/objects/Elements/PentaRef.cpp

    r6987 r7000  
    134134                *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=l1l6[i];
    135135        }
    136 
    137136}
    138137/*}}}*/
     
    245244                *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=dh1dh7[0][i];
    246245                *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
    247                 *(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=0;
    248                 *(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=0;
    249                 *(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
    250246        }
    251247
  • issm/trunk/src/c/objects/Node.cpp

    r6412 r7000  
    104104                        if (iomodel->vertices_type[io_index]==MacAyealStokesApproximationEnum && iomodel->gridonmacayeal[io_index]){
    105105                                if(!iomodel->gridonbed[io_index]){
    106                                         for(k=1;k<=gsize;k++) this->FreezeDof(k);
     106                                        for(k=1;k<=2;k++) this->FreezeDof(k);
    107107                                }
    108108                        }
Note: See TracChangeset for help on using the changeset viewer.