Changeset 15717


Ignore:
Timestamp:
08/06/13 08:40:49 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: extended arbitrary number of nodes to other solutions

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15715 r15717  
    66066606ElementMatrix* Tria::CreateKMatrixPrognostic_DG(void){
    66076607
    6608         /*Constants*/
    6609         const int    numnodes=NUMVERTICES;
    6610 
    66116608        /*Intermediaries */
    66126609        int        dim;
    66136610        IssmDouble xyz_list[NUMVERTICES][3];
    6614         IssmDouble Jdettria,dt,vx,vy;
    6615         IssmDouble basis[NUMVERTICES];
    6616         IssmDouble B[2][NUMVERTICES];
    6617         IssmDouble Bprime[2][NUMVERTICES];
    6618         IssmDouble DL[2][2]={0.0};
    6619         IssmDouble DLprime[2][2]={0.0};
    6620         IssmDouble DL_scalar;
    6621         GaussTria  *gauss=NULL;
    6622 
    6623         /*Initialize Element matrix*/
    6624         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     6611        IssmDouble Jdettria,D_scalar,dt,vx,vy;
     6612
     6613        /*Fetch number of nodes for this finite element*/
     6614        int numnodes = this->NumberofNodes();
     6615
     6616        /*Initialize Element matrix and vectors*/
     6617        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     6618        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     6619        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
     6620        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
     6621        IssmDouble     D[2][2];
    66256622
    66266623        /*Retrieve all inputs and parameters*/
    … …  
    66406637
    66416638        /* Start  looping on the number of gaussian points: */
    6642         gauss=new GaussTria(2);
     6639        GaussTria* gauss=new GaussTria(2);
    66436640        for(int ig=gauss->begin();ig<gauss->end();ig++){
    66446641
    … …  
    66496646
    66506647                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    6651                 GetNodalFunctions(&basis[0],gauss);
    6652 
    6653                 DL_scalar=gauss->weight*Jdettria;
    6654 
    6655                 TripleMultiply(&basis[0],1,numnodes,1,
    6656                                         &DL_scalar,1,1,0,
    6657                                         &basis[0],1,numnodes,0,
     6648                GetNodalFunctions(basis,gauss);
     6649
     6650                D_scalar=gauss->weight*Jdettria;
     6651
     6652                TripleMultiply(basis,1,numnodes,1,
     6653                                        &D_scalar,1,1,0,
     6654                                        basis,1,numnodes,0,
    66586655                                        &Ke->values[0],1);
    66596656
    66606657                /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
    6661                 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    6662                 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
    6663 
    6664                 DL_scalar=-dt*gauss->weight*Jdettria;
    6665 
    6666                 DLprime[0][0]=DL_scalar*vx;
    6667                 DLprime[1][1]=DL_scalar*vy;
    6668 
    6669                 TripleMultiply( &B[0][0],2,numnodes,1,
    6670                                         &DLprime[0][0],2,2,0,
    6671                                         &Bprime[0][0],2,numnodes,0,
     6658                GetBPrognostic(Bprime, &xyz_list[0][0], gauss);
     6659                GetBprimePrognostic(B, &xyz_list[0][0], gauss);
     6660
     6661                D_scalar=-dt*gauss->weight*Jdettria;
     6662                D[0][0]=D_scalar*vx;
     6663                D[0][1]=0.;
     6664                D[1][0]=0.;
     6665                D[1][1]=D_scalar*vy;
     6666
     6667                TripleMultiply(B,2,numnodes,1,
     6668                                        &D[0][0],2,2,0,
     6669                                        Bprime,2,numnodes,0,
    66726670                                        &Ke->values[0],1);
    66736671        }
    66746672
    66756673        /*Clean up and return*/
     6674        xDelete<IssmDouble>(basis);
     6675        xDelete<IssmDouble>(B);
     6676        xDelete<IssmDouble>(Bprime);
    66766677        delete gauss;
    66776678        return Ke;
    … …  
    67446745ElementVector* Tria::CreatePVectorPrognostic_DG(void){
    67456746
    6746         /*Constants*/
    6747         const int    numnodes=NUMVERTICES;
    6748 
    67496747        /*Intermediaries */
    67506748        IssmDouble Jdettria,dt;
    67516749        IssmDouble surface_mass_balance_g,basal_melting_g,thickness_g;
    67526750        IssmDouble xyz_list[NUMVERTICES][3];
    6753         IssmDouble basis[NUMVERTICES];
    6754         GaussTria* gauss=NULL;
    6755 
    6756         /*Initialize Element vector*/
    6757         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     6751
     6752        /*Fetch number of nodes and dof for this finite element*/
     6753        int numnodes = this->NumberofNodes();
     6754
     6755        /*Initialize Element vector and other vectors*/
     6756        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     6757        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    67586758
    67596759        /*Retrieve all inputs and parameters*/
    … …  
    67656765
    67666766        /* Start  looping on the number of gaussian points: */
    6767         gauss=new GaussTria(2);
     6767        GaussTria* gauss=new GaussTria(2);
    67686768        for(int ig=gauss->begin();ig<gauss->end();ig++){
    67696769
    … …  
    67716771
    67726772                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    6773                 GetNodalFunctions(&basis[0],gauss);
     6773                GetNodalFunctions(basis,gauss);
    67746774
    67756775                surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
    … …  
    67816781
    67826782        /*Clean up and return*/
     6783        xDelete<IssmDouble>(basis);
    67836784        delete gauss;
    67846785        return pe;
    … …  
    70557056ElementMatrix* Tria::CreateKMatrixBalancethickness_DG(void){
    70567057
    7057         /*Constants*/
    7058         const int  numnodes=NUMVERTICES;
    7059 
    70607058        /*Intermediaries*/
    7061         int        i,j,dim;
    7062         IssmDouble vx,vy,Jdettria;
     7059        int        dim;
     7060        IssmDouble vx,vy,D_scalar,Jdettria;
    70637061        IssmDouble xyz_list[NUMVERTICES][3];
    7064         IssmDouble B[2][NUMVERTICES];
    7065         IssmDouble Bprime[2][NUMVERTICES];
    7066         IssmDouble DL[2][2]={0.0};
    7067         IssmDouble DL_scalar;
    7068         GaussTria  *gauss=NULL;
    7069 
    7070         /*Initialize Element matrix*/
    7071         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     7062
     7063        /*Fetch number of nodes for this finite element*/
     7064        int numnodes = this->NumberofNodes();
     7065
     7066        /*Initialize Element matrix and vectors*/
     7067        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     7068        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     7069        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
     7070        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
     7071        IssmDouble     D[2][2];
    70727072
    70737073        /*Retrieve all inputs and parameters*/
    … …  
    70787078
    70797079        /*Start looping on the number of gaussian points:*/
    7080         gauss=new GaussTria(2);
     7080        GaussTria* gauss=new GaussTria(2);
    70817081        for(int ig=gauss->begin();ig<gauss->end();ig++){
    70827082
    … …  
    70857085                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    70867086                /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
    7087                 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    7088                 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
     7087                GetBPrognostic(Bprime,&xyz_list[0][0],gauss);
     7088                GetBprimePrognostic(B,&xyz_list[0][0],gauss);
    70897089
    70907090                vx_input->GetInputValue(&vx,gauss);
    70917091                vy_input->GetInputValue(&vy,gauss);
    70927092
    7093                 DL_scalar=-gauss->weight*Jdettria;
    7094                 DL[0][0]=DL_scalar*vx;
    7095                 DL[1][1]=DL_scalar*vy;
    7096 
    7097                 TripleMultiply( &B[0][0],2,numnodes,1,
    7098                                         &DL[0][0],2,2,0,
    7099                                         &Bprime[0][0],2,numnodes,0,
     7093                D_scalar=-gauss->weight*Jdettria;
     7094                D[0][0]=D_scalar*vx;
     7095                D[0][1]=0.;
     7096                D[1][0]=0.;
     7097                D[1][1]=D_scalar*vy;
     7098
     7099                TripleMultiply(B,2,numnodes,1,
     7100                                        &D[0][0],2,2,0,
     7101                                        Bprime,2,numnodes,0,
    71007102                                        &Ke->values[0],1);
    71017103        }
    71027104
    71037105        /*Clean up and return*/
     7106        xDelete<IssmDouble>(basis);
     7107        xDelete<IssmDouble>(B);
     7108        xDelete<IssmDouble>(Bprime);
    71047109        delete gauss;
    71057110        return Ke;
    … …  
    71657170ElementVector* Tria::CreatePVectorBalancethickness_DG(void){
    71667171
    7167         /*Constants*/
    7168         const int    numnodes=NUMVERTICES;
    7169 
    71707172        /*Intermediaries */
    7171         int        i,j;
    71727173        IssmDouble xyz_list[NUMVERTICES][3];
    71737174        IssmDouble basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria;
    7174         IssmDouble basis[NUMVERTICES];
    7175         GaussTria* gauss=NULL;
    7176 
    7177         /*Initialize Element vector*/
    7178         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     7175
     7176        /*Fetch number of nodes and dof for this finite element*/
     7177        int numnodes = this->NumberofNodes();
     7178
     7179        /*Initialize Element vector and other vectors*/
     7180        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     7181        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    71797182
    71807183        /*Retrieve all inputs and parameters*/
    … …  
    71857188
    71867189        /* Start  looping on the number of gaussian points: */
    7187         gauss=new GaussTria(2);
     7190        GaussTria* gauss=new GaussTria(2);
    71887191        for(int ig=gauss->begin();ig<gauss->end();ig++){
    71897192
    … …  
    71957198
    71967199                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    7197                 GetNodalFunctions(&basis[0],gauss);
    7198 
    7199                 for(i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
     7200                GetNodalFunctions(basis,gauss);
     7201
     7202                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
    72007203        }
    72017204
    72027205        /*Clean up and return*/
     7206        xDelete<IssmDouble>(basis);
    72037207        delete gauss;
    72047208        return pe;
Note: See TracChangeset for help on using the changeset viewer.