Changeset 15486


Ignore:
Timestamp:
07/11/13 16:30:33 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: mass matrix is now dynamically allocated for P2 elements

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15485 r15486  
    247247        /*Fetch number of nodes and dof for this finite element*/
    248248        int numnodes = this->NumberofNodes();
    249         int numdof   = numnodes*NDOF2;
     249        int numdof   = numnodes*NDOF1;
    250250
    251251        /*Initialize Element matrix and vectors*/
    252252        ElementMatrix* Ke    = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    253         IssmDouble*    basis = xNew<IssmDouble>(numdof);
     253        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    254254
    255255        /*Retrieve all inputs and parameters*/
    … …  
    266266                D=gauss->weight*Jdet;
    267267
    268                 TripleMultiply(basis,1,numnodes,1,
     268                TripleMultiply(basis,1,numdof,1,
    269269                                        &D,1,1,0,
    270                                         basis,1,numnodes,0,
     270                                        basis,1,numdof,0,
    271271                                        &Ke->values[0],1);
    272272        }
    … …  
    62446244
    62456245        switch(GetElementType()){
    6246                 case P1Enum:
     6246                case P1Enum: case P2Enum:
    62476247                        return CreateKMatrixPrognostic_CG();
    62486248                case P1DGEnum:
    … …  
    62566256ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){
    62576257
    6258         /*Constants*/
    6259         const int    numdof=NDOF1*NUMVERTICES;
    6260 
    62616258        /*Intermediaries */
    62626259        int        stabilization;
    62636260        int        dim;
    6264         IssmDouble Jdettria,DL_scalar,dt,h;
     6261        IssmDouble Jdettria,D_scalar,dt,h;
    62656262        IssmDouble vel,vx,vy,dvxdx,dvydy;
    62666263        IssmDouble dvx[2],dvy[2];
    62676264        IssmDouble v_gauss[2]={0.0};
    62686265        IssmDouble xyz_list[NUMVERTICES][3];
    6269         IssmDouble basis[NUMVERTICES];
    6270         IssmDouble B[2][NUMVERTICES];
    6271         IssmDouble Bprime[2][NUMVERTICES];
    6272         IssmDouble K[2][2]                        ={0.0};
    6273         IssmDouble KDL[2][2]                      ={0.0};
    6274         IssmDouble DL[2][2]                        ={0.0};
    6275         IssmDouble DLprime[2][2]                   ={0.0};
    6276 
    6277         /*Initialize Element matrix*/
    6278         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     6266
     6267        /*Fetch number of nodes and dof for this finite element*/
     6268        int numnodes = this->NumberofNodes();
     6269        int numdof   = numnodes*NDOF1;
     6270
     6271        /*Initialize Element matrix and vectors*/
     6272        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     6273        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     6274        IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
     6275        IssmDouble*    Bprime = xNew<IssmDouble>(2*numdof);
     6276        IssmDouble     D[2][2];
    62796277
    62806278        /*Retrieve all inputs and parameters*/
    … …  
    63026300
    63036301                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    6304                 GetNodalFunctions(&basis[0],gauss);
     6302                GetNodalFunctions(basis,gauss);
    63056303
    63066304                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    63096307                vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    63106308
    6311                 DL_scalar=gauss->weight*Jdettria;
    6312 
    6313                 TripleMultiply(&basis[0],1,numdof,1,
    6314                                         &DL_scalar,1,1,0,
    6315                                         &basis[0],1,numdof,0,
     6309                D_scalar=gauss->weight*Jdettria;
     6310
     6311                TripleMultiply(basis,1,numdof,1,
     6312                                        &D_scalar,1,1,0,
     6313                                        basis,1,numdof,0,
    63166314                                        &Ke->values[0],1);
    63176315
    6318                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
    6319                 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     6316                GetBPrognostic(B,&xyz_list[0][0],gauss);
     6317                GetBprimePrognostic(Bprime,&xyz_list[0][0],gauss);
    63206318
    63216319                dvxdx=dvx[0];
    63226320                dvydy=dvy[1];
    6323                 DL_scalar=dt*gauss->weight*Jdettria;
    6324 
    6325                 DL[0][0]=DL_scalar*dvxdx;
    6326                 DL[1][1]=DL_scalar*dvydy;
    6327                 DLprime[0][0]=DL_scalar*vx;
    6328                 DLprime[1][1]=DL_scalar*vy;
    6329 
    6330                 TripleMultiply( &B[0][0],2,numdof,1,
    6331                                         &DL[0][0],2,2,0,
    6332                                         &B[0][0],2,numdof,0,
     6321                D_scalar=dt*gauss->weight*Jdettria;
     6322
     6323                D[0][0]=D_scalar*dvxdx;
     6324                D[0][1]=0.;
     6325                D[1][1]=D_scalar*dvydy;
     6326                D[1][0]=0.;
     6327                TripleMultiply(B,2,numdof,1,
     6328                                        &D[0][0],2,2,0,
     6329                                        B,2,numdof,0,
    63336330                                        &Ke->values[0],1);
    63346331
    6335                 TripleMultiply( &B[0][0],2,numdof,1,
    6336                                         &DLprime[0][0],2,2,0,
    6337                                         &Bprime[0][0],2,numdof,0,
     6332                D[0][0]=D_scalar*vx;
     6333                D[1][1]=D_scalar*vy;
     6334                TripleMultiply(B,2,numdof,1,
     6335                                        &D[0][0],2,2,0,
     6336                                        Bprime,2,numdof,0,
    63386337                                        &Ke->values[0],1);
    63396338
    … …  
    63416340                        /*Streamline upwinding*/
    63426341                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
    6343                         K[0][0]=h/(2*vel)*vx*vx;
    6344                         K[1][0]=h/(2*vel)*vy*vx;
    6345                         K[0][1]=h/(2*vel)*vx*vy;
    6346                         K[1][1]=h/(2*vel)*vy*vy;
     6342                        D[0][0]=h/(2*vel)*vx*vx;
     6343                        D[1][0]=h/(2*vel)*vy*vx;
     6344                        D[0][1]=h/(2*vel)*vx*vy;
     6345                        D[1][1]=h/(2*vel)*vy*vy;
    63476346                }
    63486347                else if(stabilization==1){
    … …  
    63506349                        vxaverage_input->GetInputAverage(&vx);
    63516350                        vyaverage_input->GetInputAverage(&vy);
    6352                         K[0][0]=h/2.0*fabs(vx);
    6353                         K[0][1]=0.;
    6354                         K[1][0]=0.;
    6355                         K[1][1]=h/2.0*fabs(vy);
     6351                        D[0][0]=h/2.0*fabs(vx);
     6352                        D[0][1]=0.;
     6353                        D[1][0]=0.;
     6354                        D[1][1]=h/2.0*fabs(vy);
    63566355                }
    63576356                if(stabilization==1 || stabilization==2){
    6358                         KDL[0][0]=DL_scalar*K[0][0];
    6359                         KDL[1][0]=DL_scalar*K[1][0];
    6360                         KDL[0][1]=DL_scalar*K[0][1];
    6361                         KDL[1][1]=DL_scalar*K[1][1];
    6362                         TripleMultiply( &Bprime[0][0],2,numdof,1,
    6363                                                 &KDL[0][0],2,2,0,
    6364                                                 &Bprime[0][0],2,numdof,0,
     6357                        D[0][0]=D_scalar*D[0][0];
     6358                        D[1][0]=D_scalar*D[1][0];
     6359                        D[0][1]=D_scalar*D[0][1];
     6360                        D[1][1]=D_scalar*D[1][1];
     6361                        TripleMultiply(Bprime,2,numdof,1,
     6362                                                &D[0][0],2,2,0,
     6363                                                Bprime,2,numdof,0,
    63656364                                                &Ke->values[0],1);
    63666365                }
    … …  
    63696368        /*Clean up and return*/
    63706369        delete gauss;
     6370        xDelete<IssmDouble>(basis);
     6371        xDelete<IssmDouble>(B);
     6372        xDelete<IssmDouble>(Bprime);
    63716373        return Ke;
    63726374}
    … …  
    64516453
    64526454        switch(GetElementType()){
    6453                 case P1Enum:
     6455                case P1Enum: case P2Enum:
    64546456                        return CreatePVectorPrognostic_CG();
    64556457                case P1DGEnum:
    … …  
    64636465ElementVector* Tria::CreatePVectorPrognostic_CG(void){
    64646466
    6465         /*Constants*/
    6466         const int    numdof=NDOF1*NUMVERTICES;
    6467 
    64686467        /*Intermediaries */
    64696468        IssmDouble Jdettria,dt;
    64706469        IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g;
    64716470        IssmDouble xyz_list[NUMVERTICES][3];
    6472         IssmDouble basis[NUMVERTICES];
    6473         GaussTria* gauss=NULL;
    6474 
    6475         /*Initialize Element vector*/
    6476         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     6471
     6472        /*Fetch number of nodes and dof for this finite element*/
     6473        int numnodes = this->NumberofNodes();
     6474        int numdof   = numnodes*NDOF1;
     6475
     6476        /*Initialize Element vector and other vectors*/
     6477        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     6478        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    64776479
    64786480        /*Retrieve all inputs and parameters*/
    64796481        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    64806482        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6481         Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
    6482         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
    6483         Input* basal_melting_correction_input=inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
    6484         Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);
     6483        Input* surface_mass_balance_input     = inputs->GetInput(SurfaceforcingsMassBalanceEnum);        _assert_(surface_mass_balance_input);
     6484        Input* basal_melting_input            = inputs->GetInput(BasalforcingsMeltingRateEnum);           _assert_(basal_melting_input);
     6485        Input* basal_melting_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
     6486        Input* thickness_input                = inputs->GetInput(ThicknessEnum);                          _assert_(thickness_input);
    64856487
    64866488        /*Initialize basal_melting_correction_g to 0, do not forget!:*/
    64876489        /* Start  looping on the number of gaussian points: */
    6488         gauss=new GaussTria(2);
     6490        GaussTria* gauss=new GaussTria(2);
    64896491        for(int ig=gauss->begin();ig<gauss->end();ig++){
    64906492
    … …  
    64926494
    64936495                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    6494                 GetNodalFunctions(&basis[0],gauss);
     6496                GetNodalFunctions(basis,gauss);
    64956497
    64966498                surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
    … …  
    65076509        /*Clean up and return*/
    65086510        delete gauss;
     6511        xDelete<IssmDouble>(basis);
    65096512        return pe;
    65106513}
Note: See TracChangeset for help on using the changeset viewer.