Changeset 15713


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

CHG: extending new framework of multiple nodes to other solutions

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15712 r15713  
    407407ElementVector* Tria::CreatePVectorSlope(void){
    408408
    409         /*Constants*/
    410         const int    numdof=NDOF1*NUMVERTICES;
    411 
    412409        /*Intermediaries */
    413         int        i;
    414         int        analysis_type;
     410        int        i,analysis_type;
    415411        IssmDouble Jdet;
    416412        IssmDouble xyz_list[NUMVERTICES][3];
    417413        IssmDouble slope[2];
    418         IssmDouble basis[3];
    419         GaussTria* gauss=NULL;
     414
     415        /*Fetch number of nodes and dof for this finite element*/
     416        int numnodes = this->NumberofNodes();
    420417
    421418        /*Initialize Element vector*/
    422         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     419        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     420        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    423421
    424422        /*Retrieve all inputs and parameters*/
    … …  
    434432
    435433        /* Start  looping on the number of gaussian points: */
    436         gauss=new GaussTria(2);
     434        GaussTria* gauss=new GaussTria(2);
    437435        for(int ig=gauss->begin();ig<gauss->end();ig++){
    438436
    … …  
    440438
    441439                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    442                 GetNodalFunctions(basis, gauss);
     440                GetNodalFunctions(basis,gauss);
    443441
    444442                slope_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    445443
    446                 if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){
    447                         for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[0]*basis[i];
    448                 }
    449                 if ( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){
    450                         for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[1]*basis[i];
     444                if( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){
     445                        for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*slope[0]*basis[i];
     446                }
     447                if( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){
     448                        for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*slope[1]*basis[i];
    451449                }
    452450        }
    453451
    454452        /*Clean up and return*/
     453        xDelete<IssmDouble>(basis);
    455454        delete gauss;
    456455        return pe;
    … …  
    31573156        /*Initialize Element vector and vectors*/
    31583157        ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,SSAApproximationEnum);
    3159         GaussTria*     gauss  = new GaussTria(2);
    31603158        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    31613159
    … …  
    31683166
    31693167        /* Start  looping on the number of gaussian points: */
     3168        GaussTria*     gauss  = new GaussTria(2);
    31703169        for(int ig=gauss->begin();ig<gauss->end();ig++){
    31713170
    … …  
    49174916ElementVector* Tria::CreatePVectorAdjointBalancethickness(void){
    49184917
    4919         /*Constants*/
    4920         const int    numdof=1*NUMVERTICES;
    4921 
    49224918        /*Intermediaries */
    4923         int         i,resp;
    4924         IssmDouble  Jdet;
    4925         IssmDouble  thickness,thicknessobs,weight;
    4926         int         num_responses;
    4927         IssmDouble  xyz_list[NUMVERTICES][3];
    4928         IssmDouble  basis[3];
    4929         IssmDouble  dbasis[NDOF2][NUMVERTICES];
    4930         IssmDouble  dH[2];
    4931         IssmDouble  vx,vy,vel;
    4932         GaussTria *gauss     = NULL;
     4919        int        i,resp;
     4920        IssmDouble Jdet;
     4921        IssmDouble thickness,thicknessobs,weight;
     4922        int        num_responses;
     4923        IssmDouble xyz_list[NUMVERTICES][3];
     4924        IssmDouble dH[2];
     4925        IssmDouble vx,vy,vel;
    49334926        int       *responses = NULL;
    49344927
    4935         /*Initialize Element vector*/
    4936         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     4928        /*Fetch number of nodes and dof for this finite element*/
     4929        int numnodes = this->NumberofNodes();
     4930
     4931        /*Initialize Element vector and vectors*/
     4932        ElementVector* pe     = new ElementVector(nodes,numnodes,this->parameters);
     4933        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     4934        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    49374935
    49384936        /*Retrieve all inputs and parameters*/
    … …  
    49474945
    49484946        /* Start  looping on the number of gaussian points: */
    4949         gauss=new GaussTria(2);
     4947        GaussTria* gauss=new GaussTria(2);
    49504948        for(int ig=gauss->begin();ig<gauss->end();ig++){
    49514949
    … …  
    49544952                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    49554953                GetNodalFunctions(basis, gauss);
    4956                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
     4954                GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
    49574955
    49584956                thickness_input->GetInputValue(&thickness, gauss);
    … …  
    49654963                        case ThicknessAbsMisfitEnum:
    49664964                                weights_input->GetInputValue(&weight, gauss,resp);
    4967                                 for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
     4965                                for(i=0;i<numnodes;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
    49684966                                break;
    49694967                        case ThicknessAbsGradientEnum:
    49704968                                weights_input->GetInputValue(&weight, gauss,resp);
    4971                                 for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[0]*dbasis[0][i]*Jdet*gauss->weight;
    4972                                 for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[1]*dbasis[1][i]*Jdet*gauss->weight;
     4969                                for(i=0;i<numnodes;i++) pe->values[i]+= - weight*dH[0]*dbasis[0*numnodes+i]*Jdet*gauss->weight;
     4970                                for(i=0;i<numnodes;i++) pe->values[i]+= - weight*dH[1]*dbasis[1*numnodes+i]*Jdet*gauss->weight;
    49734971                                break;
    49744972                        case ThicknessAlongGradientEnum:
    … …  
    49794977                                vx  = vx/(vel+1.e-9);
    49804978                                vy  = vy/(vel+1.e-9);
    4981                                 for(i=0;i<numdof;i++) pe->values[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0][i]*vx+dbasis[1][i]*vy)*Jdet*gauss->weight;
     4979                                for(i=0;i<numnodes;i++) pe->values[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0*numnodes+i]*vx+dbasis[1*numnodes+i]*vy)*Jdet*gauss->weight;
    49824980                                break;
    49834981                        case ThicknessAcrossGradientEnum:
    … …  
    49884986                                vx  = vx/(vel+1.e-9);
    49894987                                vy  = vy/(vel+1.e-9);
    4990                                 for(i=0;i<numdof;i++) pe->values[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0][i]*(-vy)+dbasis[1][i]*vx)*Jdet*gauss->weight;
     4988                                for(i=0;i<numnodes;i++) pe->values[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0*numnodes+i]*(-vy)+dbasis[1*numnodes+i]*vx)*Jdet*gauss->weight;
    49914989                                break;
    49924990                        default:
    … …  
    49964994
    49974995        /*Clean up and return*/
     4996        xDelete<IssmDouble>(basis);
     4997        xDelete<IssmDouble>(dbasis);
     4998        xDelete<int>(responses);
    49984999        delete gauss;
    4999         xDelete<int>(responses);
    50005000        return pe;
    50015001}
    … …  
    58745874        /*Fetch number of nodes and dof for this finite element*/
    58755875        int numnodes = this->NumberofNodes();
    5876         int numdof   = numnodes*NDOF2;
     5876        int numdof   = numnodes*NDOF1;
    58775877
    58785878        /*Initialize Element matrix and vectors*/
    … …  
    59285928ElementMatrix* Tria::CreateKMatrixHydrologyDCEfficient(void){
    59295929
    5930         /*constants: */
    5931         const int    numdof=NDOF1*NUMVERTICES;
    5932 
    59335930        /* Intermediaries */
    59345931        IssmDouble  D_scalar,Jdet;
    … …  
    59365933        IssmDouble  epl_storing;
    59375934        IssmDouble  xyz_list[NUMVERTICES][3];
    5938         IssmDouble  B[2][numdof];
    5939         IssmDouble  basis[NUMVERTICES];
    5940         IssmDouble  D[2][2];
    5941         GaussTria   *gauss = NULL;
    59425935
    59435936        /*Check that all nodes are active, else return empty matrix*/
    … …  
    59465939        }
    59475940
    5948         /*Initialize Element matrix*/
    5949         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     5941        /*Fetch number of nodes and dof for this finite element*/
     5942        int numnodes = this->NumberofNodes();
     5943        int numdof   = numnodes*NDOF1;
     5944
     5945        /*Initialize Element matrix and vectors*/
     5946        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
     5947        IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
     5948        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     5949        IssmDouble     D[2][2];
    59505950
    59515951        /*Retrieve all inputs and parameters*/
    … …  
    59565956
    59575957        /* Start looping on the number of gaussian points: */
    5958         gauss=new GaussTria(2);
     5958        GaussTria* gauss=new GaussTria(2);
    59595959        for(int ig=gauss->begin();ig<gauss->end();ig++){
    59605960
    … …  
    59675967                D[0][0]=D_scalar; D[0][1]=0.;
    59685968                D[1][0]=0.;       D[1][1]=D_scalar;
    5969                 GetBHydro(&B[0][0],&xyz_list[0][0],gauss);
    5970                 TripleMultiply(&B[0][0],2,numdof,1,
     5969                GetBHydro(B,&xyz_list[0][0],gauss);
     5970                TripleMultiply(B,2,numdof,1,
    59715971                                        &D[0][0],2,2,0,
    5972                                         &B[0][0],2,numdof,0,
     5972                                        B,2,numdof,0,
    59735973                                        &Ke->values[0],1);
    59745974
    59755975                /*Transient*/
    59765976                if(reCast<bool,IssmDouble>(dt)){
    5977                         GetNodalFunctions(&basis[0],gauss);
     5977                        GetNodalFunctions(basis,gauss);
    59785978                        D_scalar=epl_storing*gauss->weight*Jdet;
    59795979
    5980                         TripleMultiply(&basis[0],numdof,1,0,
     5980                        TripleMultiply(basis,numdof,1,0,
    59815981                                                &D_scalar,1,1,0,
    5982                                                 &basis[0],1,numdof,0,
     5982                                                basis,1,numdof,0,
    59835983                                                &Ke->values[0],1);
    59845984                }
    59855985        }
     5986
    59865987        /*Clean up and return*/
     5988        xDelete<IssmDouble>(basis);
     5989        xDelete<IssmDouble>(B);
    59875990        delete gauss;
    59885991        return Ke;
    … …  
    67346737
    67356738        /*Clean up and return*/
     6739        xDelete<IssmDouble>(basis);
    67366740        delete gauss;
    6737         xDelete<IssmDouble>(basis);
    67386741        return pe;
    67396742}
    … …  
    69446947ElementMatrix* Tria::CreateKMatrixBalancethickness_CG(void){
    69456948
    6946         /*Constants*/
    6947         const int    numdof=NDOF1*NUMVERTICES;
    6948 
    69496949        /*Intermediaries */
    69506950        int        stabilization;
    69516951        int        i,j,dim;
    69526952        IssmDouble Jdettria,vx,vy,dvxdx,dvydy,vel,h;
     6953        IssmDouble D_scalar;
    69536954        IssmDouble dvx[2],dvy[2];
    69546955        IssmDouble xyz_list[NUMVERTICES][3];
    6955         IssmDouble L[NUMVERTICES];
    6956         IssmDouble B[2][NUMVERTICES];
    6957         IssmDouble Bprime[2][NUMVERTICES];
    6958         IssmDouble K[2][2]                          = {0.0};
    6959         IssmDouble KDL[2][2]                        = {0.0};
    6960         IssmDouble DL[2][2]                         = {0.0};
    6961         IssmDouble DLprime[2][2]                    = {0.0};
    6962         IssmDouble DL_scalar;
    6963         GaussTria *gauss                            = NULL;
    6964 
    6965         /*Initialize Element matrix*/
    6966         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     6956
     6957        /*Fetch number of nodes and dof for this finite element*/
     6958        int numnodes = this->NumberofNodes();
     6959        int numdof   = numnodes*NDOF1;
     6960
     6961        /*Initialize Element matrix and vectors*/
     6962        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     6963        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     6964        IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
     6965        IssmDouble*    Bprime = xNew<IssmDouble>(2*numdof);
     6966        IssmDouble     D[2][2];
    69676967
    69686968        /*Retrieve all Inputs and parameters: */
    … …  
    69836983
    69846984        /*Start looping on the number of gaussian points:*/
    6985         gauss=new GaussTria(2);
     6985        GaussTria* gauss=new GaussTria(2);
    69866986        for(int ig=gauss->begin();ig<gauss->end();ig++){
    69876987
    … …  
    69896989
    69906990                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    6991                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
    6992                 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     6991                GetBPrognostic(B,&xyz_list[0][0],gauss);
     6992                GetBprimePrognostic(Bprime,&xyz_list[0][0],gauss);
    69936993
    69946994                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    69996999                dvxdx=dvx[0];
    70007000                dvydy=dvy[1];
    7001                 DL_scalar=gauss->weight*Jdettria;
    7002 
    7003                 DL[0][0]=DL_scalar*dvxdx;
    7004                 DL[1][1]=DL_scalar*dvydy;
    7005 
    7006                 DLprime[0][0]=DL_scalar*vx;
    7007                 DLprime[1][1]=DL_scalar*vy;
    7008 
    7009                 TripleMultiply( &B[0][0],2,numdof,1,
    7010                                         &DL[0][0],2,2,0,
    7011                                         &B[0][0],2,numdof,0,
     7001                D_scalar=gauss->weight*Jdettria;
     7002
     7003                D[0][0]=D_scalar*dvxdx;
     7004                D[0][1]=0.;
     7005                D[1][1]=D_scalar*dvydy;
     7006                D[1][1]=0.;
     7007                TripleMultiply(B,2,numdof,1,
     7008                                        &D[0][0],2,2,0,
     7009                                        B,2,numdof,0,
    70127010                                        &Ke->values[0],1);
    70137011
    7014                 TripleMultiply( &B[0][0],2,numdof,1,
    7015                                         &DLprime[0][0],2,2,0,
    7016                                         &Bprime[0][0],2,numdof,0,
     7012                D[0][0]=D_scalar*vx;
     7013                D[1][1]=D_scalar*vy;
     7014                TripleMultiply(B,2,numdof,1,
     7015                                        &D[0][0],2,2,0,
     7016                                        Bprime,2,numdof,0,
    70177017                                        &Ke->values[0],1);
    70187018
    … …  
    70207020                        /*Streamline upwinding*/
    70217021                        vel=sqrt(vx*vx+vy*vy);
    7022                         K[0][0]=h/(2*vel)*vx*vx;
    7023                         K[1][0]=h/(2*vel)*vy*vx;
    7024                         K[0][1]=h/(2*vel)*vx*vy;
    7025                         K[1][1]=h/(2*vel)*vy*vy;
     7022                        D[0][0]=h/(2*vel)*vx*vx;
     7023                        D[1][0]=h/(2*vel)*vy*vx;
     7024                        D[0][1]=h/(2*vel)*vx*vy;
     7025                        D[1][1]=h/(2*vel)*vy*vy;
    70267026                }
    70277027                else if(stabilization==2){
    … …  
    70297029                        vxaverage_input->GetInputAverage(&vx);
    70307030                        vyaverage_input->GetInputAverage(&vy);
    7031                         K[0][0]=h/2.0*fabs(vx);
    7032                         K[0][1]=0.;
    7033                         K[1][0]=0.;
    7034                         K[1][1]=h/2.0*fabs(vy);
     7031                        D[0][0]=h/2.0*fabs(vx);
     7032                        D[0][1]=0.;
     7033                        D[1][0]=0.;
     7034                        D[1][1]=h/2.0*fabs(vy);
    70357035                }
    70367036                if(stabilization==1 || stabilization==2){
    7037                         KDL[0][0]=DL_scalar*K[0][0];
    7038                         KDL[1][0]=DL_scalar*K[1][0];
    7039                         KDL[0][1]=DL_scalar*K[0][1];
    7040                         KDL[1][1]=DL_scalar*K[1][1];
    7041                         TripleMultiply( &Bprime[0][0],2,numdof,1,
    7042                                                 &KDL[0][0],2,2,0,
    7043                                                 &Bprime[0][0],2,numdof,0,
     7037                        D[0][0]=D_scalar*D[0][0];
     7038                        D[1][0]=D_scalar*D[1][0];
     7039                        D[0][1]=D_scalar*D[0][1];
     7040                        D[1][1]=D_scalar*D[1][1];
     7041                        TripleMultiply(Bprime,2,numdof,1,
     7042                                                &D[0][0],2,2,0,
     7043                                                Bprime,2,numdof,0,
    70447044                                                &Ke->values[0],1);
    70457045                }
    … …  
    70477047
    70487048        /*Clean up and return*/
     7049        xDelete<IssmDouble>(basis);
     7050        xDelete<IssmDouble>(B);
     7051        xDelete<IssmDouble>(Bprime);
    70497052        delete gauss;
    70507053        return Ke;
    … …  
    71227125ElementVector* Tria::CreatePVectorBalancethickness_CG(void){
    71237126
    7124         /*Constants*/
    7125         const int    numdof=NDOF1*NUMVERTICES;
    7126 
    71277127        /*Intermediaries */
    7128         int        i,j;
    71297128        IssmDouble xyz_list[NUMVERTICES][3];
    71307129        IssmDouble dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria;
    7131         IssmDouble basis[NUMVERTICES];
    7132         GaussTria* gauss=NULL;
    7133 
    7134         /*Initialize Element vector*/
    7135         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     7130
     7131        /*Fetch number of nodes and dof for this finite element*/
     7132        int numnodes = this->NumberofNodes();
     7133        int numdof   = numnodes*NDOF1;
     7134
     7135        /*Initialize Element vector and other vectors*/
     7136        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     7137        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    71367138
    71377139        /*Retrieve all inputs and parameters*/
    … …  
    71427144
    71437145        /* Start  looping on the number of gaussian points: */
    7144         gauss=new GaussTria(2);
     7146        GaussTria* gauss=new GaussTria(2);
    71457147        for(int ig=gauss->begin();ig<gauss->end();ig++){
    71467148
    … …  
    71527154
    71537155                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
    7154                 GetNodalFunctions(&basis[0],gauss);
    7155 
    7156                 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
     7156                GetNodalFunctions(basis,gauss);
     7157
     7158                for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
    71577159        }
    71587160
    71597161        /*Clean up and return*/
     7162        xDelete<IssmDouble>(basis);
    71607163        delete gauss;
    71617164        return pe;
Note: See TracChangeset for help on using the changeset viewer.