Changeset 15712


Ignore:
Timestamp:
08/05/13 15:56:17 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: extending new framework of multiple nodes to other solutions

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

Legend:

Unmodified
Added
Removed

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

    r15711 r15712  
    84418441
    84428442        /*Clean up and return*/
     8443        xDelete<IssmDouble>(basis);
    84438444        delete gauss;
    8444         xDelete<IssmDouble>(basis);
    84458445        return pe;
    84468446}

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

    r15711 r15712  
    31443144
    31453145        /*Intermediaries */
    3146         int            i,j;
    3147         IssmDouble     driving_stress_baseline,thickness;
    3148         IssmDouble     Jdet;
    3149         IssmDouble     xyz_list[NUMVERTICES][3];
    3150         IssmDouble     slope[2];
    3151         IssmDouble     icefrontlevel[3];
     3146        int        i;
     3147        IssmDouble driving_stress_baseline,thickness;
     3148        IssmDouble Jdet;
     3149        IssmDouble xyz_list[NUMVERTICES][3];
     3150        IssmDouble slope[2];
     3151        IssmDouble icefrontlevel[3];
    31523152
    31533153        /*Fetch number of nodes and dof for this finite element*/
    … …  
    50085008
    50095009        /*Intermediaries */
    5010         int        i,resp;
    5011         int       *responses=NULL;
    5012         int        num_responses;
    5013         IssmDouble     Jdet;
    5014         IssmDouble     obs_velocity_mag,velocity_mag;
    5015         IssmDouble     dux,duy;
    5016         IssmDouble     epsvel=2.220446049250313e-16;
    5017         IssmDouble     meanvel=3.170979198376458e-05; /*1000 m/yr*/
    5018         IssmDouble     scalex=0,scaley=0,scale=0,S=0;
    5019         IssmDouble     vx,vy,vxobs,vyobs,weight;
    5020         IssmDouble     xyz_list[NUMVERTICES][3];
    5021         IssmDouble     basis[3];
     5010        int         i,resp;
     5011        int        *responses=NULL;
     5012        int         num_responses;
     5013        IssmDouble Jdet;
     5014        IssmDouble obs_velocity_mag,velocity_mag;
     5015        IssmDouble dux,duy;
     5016        IssmDouble epsvel=2.220446049250313e-16;
     5017        IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
     5018        IssmDouble scalex=0,scaley=0,scale=0,S=0;
     5019        IssmDouble vx,vy,vxobs,vyobs,weight;
     5020        IssmDouble xyz_list[NUMVERTICES][3];
     5021        IssmDouble basis[3];
    50225022        GaussTria* gauss=NULL;
    50235023
    … …  
    56395639ElementMatrix* Tria::CreateKMatrixMelting(void){
    56405640
    5641         /*Constants*/
    5642         const int  numdof=NUMVERTICES*NDOF1;
    5643 
    56445641        /*Intermediaries */
    56455642        IssmDouble heatcapacity,latentheat;
    56465643        IssmDouble Jdet,D_scalar;
    56475644        IssmDouble xyz_list[NUMVERTICES][3];
    5648         IssmDouble basis[3];
    5649         GaussTria *gauss=NULL;
     5645
     5646        /*Fetch number of nodes and dof for this finite element*/
     5647        int numnodes = this->NumberofNodes();
     5648        int numdof   = numnodes*NDOF1;
     5649
     5650        /*Initialize Element vector and vectors*/
     5651        ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     5652        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    56505653
    56515654        /*Initialize Element matrix*/
    5652         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    56535655
    56545656        /*Retrieve all inputs and parameters*/
    56555657        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5656         latentheat=matpar->GetLatentHeat();
    5657         heatcapacity=matpar->GetHeatCapacity();
     5658        latentheat   = matpar->GetLatentHeat();
     5659        heatcapacity = matpar->GetHeatCapacity();
    56585660
    56595661        /* Start looping on the number of gauss  (nodes on the bedrock) */
    5660         gauss=new GaussTria(2);
     5662        GaussTria* gauss=new GaussTria(2);
    56615663        for(int ig=gauss->begin();ig<gauss->end();ig++){
    56625664
    … …  
    56685670                D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
    56695671
    5670                 TripleMultiply(&basis[0],numdof,1,0,
     5672                TripleMultiply(basis,numdof,1,0,
    56715673                                        &D_scalar,1,1,0,
    5672                                         &basis[0],1,numdof,0,
     5674                                        basis,1,numdof,0,
    56735675                                        &Ke->values[0],1);
    56745676        }
    56755677
    56765678        /*Clean up and return*/
     5679        xDelete<IssmDouble>(basis);
    56775680        delete gauss;
    56785681        return Ke;
    … …  
    57645767ElementMatrix* Tria::CreateKMatrixHydrologyShreve(void){
    57655768
    5766         /*Constants*/
    5767         const int  numdof=NDOF1*NUMVERTICES;
    5768 
    5769 /*Intermediaries */
     5769        /*Intermediaries */
    57705770        IssmDouble diffusivity;
    5771         IssmDouble Jdettria,DL_scalar,dt,h;
     5771        IssmDouble Jdettria,D_scalar,dt,h;
    57725772        IssmDouble vx,vy,vel,dvxdx,dvydy;
    57735773        IssmDouble dvx[2],dvy[2];
    5774         IssmDouble v_gauss[2]={0.0};
    57755774        IssmDouble xyz_list[NUMVERTICES][3];
    5776         IssmDouble basis[NUMVERTICES];
    5777         IssmDouble B[2][NUMVERTICES];
    5778         IssmDouble Bprime[2][NUMVERTICES];
    5779         IssmDouble K[2][2]                        ={0.0};
    5780         IssmDouble KDL[2][2]                      ={0.0};
    5781         IssmDouble DL[2][2]                        ={0.0};
    5782         IssmDouble DLprime[2][2]                   ={0.0};
    5783         GaussTria *gauss=NULL;
    5784 
    5785 /*Skip if water or ice shelf element*/
     5775
     5776        /*Skip if water or ice shelf element*/
    57865777        if(IsOnWater() | IsFloating()) return NULL;
    57875778
    5788 /*Initialize Element matrix*/
    5789         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    5790 
    5791 /*Create water velocity vx and vy from current inputs*/
     5779        /*Fetch number of nodes and dof for this finite element*/
     5780        int numnodes = this->NumberofNodes();
     5781        int numdof   = numnodes*NDOF1;
     5782
     5783        /*Initialize Element matrix and vectors*/
     5784        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     5785        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     5786        IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
     5787        IssmDouble*    Bprime = xNew<IssmDouble>(2*numdof);
     5788        IssmDouble     D[2][2];
     5789
     5790        /*Create water velocity vx and vy from current inputs*/
    57925791        CreateHydrologyWaterVelocityInput();
    57935792
    … …  
    58005799        h=sqrt(2*this->GetArea());
    58015800
    5802 /* Start  looping on the number of gaussian points: */
    5803         gauss=new GaussTria(2);
     5801        /* Start  looping on the number of gaussian points: */
     5802        GaussTria* gauss=new GaussTria(2);
    58045803        for(int ig=gauss->begin();ig<gauss->end();ig++){
    58055804
    … …  
    58075806
    58085807                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    5809                 GetNodalFunctions(&basis[0],gauss);
     5808                GetNodalFunctions(basis,gauss);
    58105809
    58115810                vx_input->GetInputValue(&vx,gauss);
    … …  
    58145813                vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    58155814
    5816                 DL_scalar=gauss->weight*Jdettria;
    5817 
    5818                 TripleMultiply(&basis[0],1,numdof,1,
    5819                                         &DL_scalar,1,1,0,
    5820                                         &basis[0],1,numdof,0,
    5821                                         &Ke->values[0],1);
    5822 
    5823                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
    5824                 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     5815                D_scalar=gauss->weight*Jdettria;
     5816
     5817                TripleMultiply(basis,1,numdof,1,
     5818                                        &D_scalar,1,1,0,
     5819                                        basis,1,numdof,0,
     5820                                        Ke->values,1);
     5821
     5822                GetBPrognostic(B,&xyz_list[0][0], gauss);
     5823                GetBprimePrognostic(Bprime,&xyz_list[0][0], gauss);
    58255824
    58265825                dvxdx=dvx[0];
    58275826                dvydy=dvy[1];
    5828                 DL_scalar=dt*gauss->weight*Jdettria;
    5829 
    5830                 DL[0][0]=DL_scalar*dvxdx;
    5831                 DL[1][1]=DL_scalar*dvydy;
    5832                 DLprime[0][0]=DL_scalar*vx;
    5833                 DLprime[1][1]=DL_scalar*vy;
    5834 
    5835                 TripleMultiply( &B[0][0],2,numdof,1,
    5836                                                                                 &DL[0][0],2,2,0,
    5837                                                                                 &B[0][0],2,numdof,0,
    5838                                                                                 &Ke->values[0],1);
    5839 
    5840                 TripleMultiply( &B[0][0],2,numdof,1,
    5841                                                                                 &DLprime[0][0],2,2,0,
    5842                                                                                 &Bprime[0][0],2,numdof,0,
    5843                                                                                 &Ke->values[0],1);
     5827                D_scalar=dt*gauss->weight*Jdettria;
     5828
     5829                D[0][0]=D_scalar*dvxdx;
     5830                D[0][1]=0.;
     5831                D[1][1]=D_scalar*dvydy;
     5832                D[1][1]=0.;
     5833                TripleMultiply(B,2,numdof,1,
     5834                                        &D[0][0],2,2,0,
     5835                                        B,2,numdof,0,
     5836                                        &Ke->values[0],1);
     5837
     5838                D[0][0]=D_scalar*vx;
     5839                D[1][1]=D_scalar*vy;
     5840                TripleMultiply(B,2,numdof,1,
     5841                                        &D[0][0],2,2,0,
     5842                                        Bprime,2,numdof,0,
     5843                                        &Ke->values[0],1);
    58445844
    58455845                /*Artificial diffusivity*/
    58465846                vel=sqrt(vx*vx+vy*vy);
    5847                 K[0][0]=diffusivity*h/(2*vel)*vx*vx;
    5848                 K[1][0]=diffusivity*h/(2*vel)*vy*vx;
    5849                 K[0][1]=diffusivity*h/(2*vel)*vx*vy;
    5850                 K[1][1]=diffusivity*h/(2*vel)*vy*vy;
    5851                 KDL[0][0]=DL_scalar*K[0][0];
    5852                 KDL[1][0]=DL_scalar*K[1][0];
    5853                 KDL[0][1]=DL_scalar*K[0][1];
    5854                 KDL[1][1]=DL_scalar*K[1][1];
    5855 
    5856                 TripleMultiply( &Bprime[0][0],2,numdof,1,
    5857                                                                                 &KDL[0][0],2,2,0,
    5858                                                                                 &Bprime[0][0],2,numdof,0,
    5859                                                                                 &Ke->values[0],1);
    5860         }
    5861 
    5862 /*Clean up and return*/
     5847                D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx;
     5848                D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx;
     5849                D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;
     5850                D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;
     5851                TripleMultiply(Bprime,2,numdof,1,
     5852                                        &D[0][0],2,2,0,
     5853                                        Bprime,2,numdof,0,
     5854                                        &Ke->values[0],1);
     5855        }
     5856
     5857        /*Clean up and return*/
     5858        xDelete<IssmDouble>(basis);
     5859        xDelete<IssmDouble>(B);
     5860        xDelete<IssmDouble>(Bprime);
    58635861        delete gauss;
    58645862        return Ke;
    … …  
    58675865/*FUNCTION Tria::CreateKMatrixHydrologyDCInefficient{{{*/
    58685866ElementMatrix* Tria::CreateKMatrixHydrologyDCInefficient(void){
    5869 
    5870         /*constants: */
    5871         const int    numdof=NDOF1*NUMVERTICES;
    58725867
    58735868        /* Intermediaries */
    … …  
    58765871        IssmDouble  sediment_storing;
    58775872        IssmDouble  xyz_list[NUMVERTICES][3];
    5878         IssmDouble  B[2][numdof];
    5879         IssmDouble  basis[NUMVERTICES];
    5880         IssmDouble  D[2][2];
    5881         GaussTria   *gauss = NULL;
    5882 
    5883         /*Initialize Element matrix*/
    5884         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     5873
     5874        /*Fetch number of nodes and dof for this finite element*/
     5875        int numnodes = this->NumberofNodes();
     5876        int numdof   = numnodes*NDOF2;
     5877
     5878        /*Initialize Element matrix and vectors*/
     5879        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
     5880        IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
     5881        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     5882        IssmDouble     D[2][2];
    58855883
    58865884        /*Retrieve all inputs and parameters*/
    … …  
    58915889
    58925890        /* Start looping on the number of gaussian points: */
    5893         gauss=new GaussTria(2);
     5891        GaussTria* gauss=new GaussTria(2);
    58945892        for(int ig=gauss->begin();ig<gauss->end();ig++){
    58955893
    … …  
    59035901                D[0][0]=D_scalar; D[0][1]=0.;
    59045902                D[1][0]=0.;       D[1][1]=D_scalar;
    5905                 GetBHydro(&B[0][0],&xyz_list[0][0],gauss);
    5906                 TripleMultiply(&B[0][0],2,numdof,1,
     5903                GetBHydro(B,&xyz_list[0][0],gauss);
     5904                TripleMultiply(B,2,numdof,1,
    59075905                                        &D[0][0],2,2,0,
    5908                                         &B[0][0],2,numdof,0,
     5906                                        B,2,numdof,0,
    59095907                                        &Ke->values[0],1);
    59105908
    … …  
    59145912                        D_scalar=sediment_storing*gauss->weight*Jdet;
    59155913
    5916                         TripleMultiply(&basis[0],numdof,1,0,
     5914                        TripleMultiply(basis,numdof,1,0,
    59175915                                                &D_scalar,1,1,0,
    5918                                                 &basis[0],1,numdof,0,
     5916                                                basis,1,numdof,0,
    59195917                                                &Ke->values[0],1);
    59205918                }
    59215919        }
    59225920        /*Clean up and return*/
     5921        xDelete<IssmDouble>(basis);
     5922        xDelete<IssmDouble>(B);
    59235923        delete gauss;
    59245924        return Ke;
    … …  
    59695969                GetBHydro(&B[0][0],&xyz_list[0][0],gauss);
    59705970                TripleMultiply(&B[0][0],2,numdof,1,
    5971                                                                          &D[0][0],2,2,0,
    5972                                                                          &B[0][0],2,numdof,0,
    5973                                                                          &Ke->values[0],1);
     5971                                        &D[0][0],2,2,0,
     5972                                        &B[0][0],2,numdof,0,
     5973                                        &Ke->values[0],1);
    59745974
    59755975                /*Transient*/
    … …  
    59925992ElementVector* Tria::CreatePVectorHydrologyShreve(void){
    59935993
    5994         /*Constants*/
    5995         const int    numdof=NDOF1*NUMVERTICES;
    5996 
    59975994        /*Intermediaries */
    5998         int        i,j;
     5995        int        i;
    59995996        IssmDouble Jdettria,dt;
    60005997        IssmDouble basal_melting_g;
    60015998        IssmDouble old_watercolumn_g;
    60025999        IssmDouble xyz_list[NUMVERTICES][3];
    6003         IssmDouble basis[numdof];
    6004         GaussTria* gauss=NULL;
    60056000
    60066001        /*Skip if water or ice shelf element*/
    60076002        if(IsOnWater() | IsFloating()) return NULL;
    60086003
     6004        /*Fetch number of nodes and dof for this finite element*/
     6005        int numnodes = this->NumberofNodes();
     6006        int numdof   = numnodes*NDOF1;
     6007
    60096008        /*Initialize Element vector*/
    6010         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     6009        ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters);
     6010        IssmDouble*    basis = xNewZeroInit<IssmDouble>(numnodes);
    60116011
    60126012        /*Retrieve all inputs and parameters*/
    … …  
    60186018        /*Initialize basal_melting_correction_g to 0, do not forget!:*/
    60196019        /* Start  looping on the number of gaussian points: */
    6020         gauss=new GaussTria(2);
     6020        GaussTria* gauss=new GaussTria(2);
    60216021        for(int ig=gauss->begin();ig<gauss->end();ig++){
    60226022
    … …  
    60296029                old_watercolumn_input->GetInputValue(&old_watercolumn_g,gauss);
    60306030
    6031                 if(reCast<int,IssmDouble>(dt))for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
    6032                 else  for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
     6031                if(reCast<int,IssmDouble>(dt)){
     6032                        for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
     6033                }
     6034                else{
     6035                        for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
     6036                }
    60336037        }
    60346038
    60356039        /*Clean up and return*/
     6040        xDelete<IssmDouble>(basis);
    60366041        delete gauss;
    60376042        return pe;
    … …  
    64846489        IssmDouble vel,vx,vy,dvxdx,dvydy;
    64856490        IssmDouble dvx[2],dvy[2];
    6486         IssmDouble v_gauss[2]={0.0};
    64876491        IssmDouble xyz_list[NUMVERTICES][3];
    64886492
    … …  
    65896593
    65906594        /*Clean up and return*/
    6591         delete gauss;
    65926595        xDelete<IssmDouble>(basis);
    65936596        xDelete<IssmDouble>(B);
    65946597        xDelete<IssmDouble>(Bprime);
     6598        delete gauss;
    65956599        return Ke;
    65966600}
Note: See TracChangeset for help on using the changeset viewer.