Changeset 15711


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

NEW: SIA now available for different interpolation

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

Legend:

Unmodified
Added
Removed

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

    r15707 r15711  
    69816981ElementMatrix* Penta::CreateKMatrixDiagnosticSIA(void){
    69826982
    6983         /*Constants*/
    6984         const int numdof=NDOF2*NUMVERTICES;
    6985 
    69866983        /*Intermediaries*/
    6987         int         connectivity[2];
    6988         int         i,i0,i1,j0,j1;
    6989         IssmDouble  one0,one1;
     6984        IssmDouble connectivity[2];
     6985        IssmDouble one0,one1;
     6986        int        i,i0,i1,j0,j1;
     6987
     6988        /*Fetch number of nodes and dof for this finite element*/
     6989        int numnodes = this->NumberofNodes(); _assert_(numnodes==6);
     6990        int numdof   = numnodes*NDOF2;
    69906991
    69916992        /*Initialize Element matrix*/
    6992         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    6993 
    6994         /*Spawn 3 beam elements: */
     6993        ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     6994
     6995        /*3 vertical edges*/
    69956996        for(i=0;i<3;i++){
     6997
    69966998                /*2 dofs of first node*/
    6997                 i0=2*i;
    6998                 i1=2*i+1;
     6999                i0=2*i;     i1=2*i+1;
    69997000                /*2 dofs of second node*/
    7000                 j0=2*(i+3);
    7001                 j1=2*(i+3)+1;
     7001                j0=2*(i+3); j1=2*(i+3)+1;
    70027002
    70037003                /*Find connectivity for the two nodes*/
    7004                 connectivity[0]=vertices[i]->Connectivity();
    7005                 connectivity[1]=vertices[i+3]->Connectivity();
    7006                 one0=1/(IssmDouble)connectivity[0];
    7007                 one1=1/(IssmDouble)connectivity[1];
     7004                connectivity[0]=(IssmDouble)vertices[i]->Connectivity();
     7005                connectivity[1]=(IssmDouble)vertices[i+3]->Connectivity();
     7006                one0=1./connectivity[0];
     7007                one1=1./connectivity[1];
    70087008
    70097009                /*Create matrix for these two nodes*/
    70107010                if (IsOnBed() && IsOnSurface()){
    7011                         Ke->values[i0*numdof+i0]=one0;
    7012                         Ke->values[i1*numdof+i1]=one0;
    7013                         Ke->values[j0*numdof+i0]=-one1;
    7014                         Ke->values[j0*numdof+j0]=one1;
    7015                         Ke->values[j1*numdof+i1]=-one1;
    7016                         Ke->values[j1*numdof+j1]=one1;
     7011                        Ke->values[i0*numdof+i0] = +one0;
     7012                        Ke->values[i1*numdof+i1] = +one0;
     7013                        Ke->values[j0*numdof+i0] = -one1;
     7014                        Ke->values[j0*numdof+j0] = +one1;
     7015                        Ke->values[j1*numdof+i1] = -one1;
     7016                        Ke->values[j1*numdof+j1] = +one1;
    70177017                }
    70187018                else if (IsOnBed()){
    7019                         Ke->values[i0*numdof+i0]=one0;
    7020                         Ke->values[i1*numdof+i1]=one0;
    7021                         Ke->values[j0*numdof+i0]=-2*one1;
    7022                         Ke->values[j0*numdof+j0]=2*one1;
    7023                         Ke->values[j1*numdof+i1]=-2*one1;
    7024                         Ke->values[j1*numdof+j1]=2*one1;
     7019                        Ke->values[i0*numdof+i0] = one0;
     7020                        Ke->values[i1*numdof+i1] = one0;
     7021                        Ke->values[j0*numdof+i0] = -2.*one1;
     7022                        Ke->values[j0*numdof+j0] = +2.*one1;
     7023                        Ke->values[j1*numdof+i1] = -2.*one1;
     7024                        Ke->values[j1*numdof+j1] = +2.*one1;
    70257025                }
    70267026                else if (IsOnSurface()){
    7027                         Ke->values[j0*numdof+i0]=-one1;
    7028                         Ke->values[j0*numdof+j0]=one1;
    7029                         Ke->values[j1*numdof+i1]=-one1;
    7030                         Ke->values[j1*numdof+j1]=one1;
     7027                        Ke->values[j0*numdof+i0] = -one1;
     7028                        Ke->values[j0*numdof+j0] = +one1;
     7029                        Ke->values[j1*numdof+i1] = -one1;
     7030                        Ke->values[j1*numdof+j1] = +one1;
    70317031                }
    70327032                else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity
    7033                         Ke->values[j0*numdof+i0]=-2*one1;
    7034                         Ke->values[j0*numdof+j0]=2*one1;
    7035                         Ke->values[j1*numdof+i1]=-2*one1;
    7036                         Ke->values[j1*numdof+j1]=2*one1;
     7033                        Ke->values[j0*numdof+i0] = -2.*one1;
     7034                        Ke->values[j0*numdof+j0] = +2.*one1;
     7035                        Ke->values[j1*numdof+i1] = -2.*one1;
     7036                        Ke->values[j1*numdof+j1] = +2.*one1;
    70377037                }
    70387038        }
    … …  
    82648264
    82658265        /*Intermediaries*/
    8266         int          i,j;
    8267         int          node0,node1;
    8268         int          connectivity[2];
    8269         IssmDouble   Jdet;
    8270         IssmDouble   xyz_list[NUMVERTICES][3];
    8271         IssmDouble   xyz_list_segment[2][3];
    8272         IssmDouble   z_list[NUMVERTICES];
    8273         IssmDouble   slope[2];
    8274         IssmDouble   slope2,constant_part;
    8275         IssmDouble   rho_ice,gravity,n,B;
    8276         IssmDouble   ub,vb,z_g,surface,thickness;
    8277         GaussPenta*  gauss=NULL;
     8266        int         i,j;
     8267        int         node0,node1;
     8268        IssmDouble  connectivity[2];
     8269        IssmDouble  Jdet;
     8270        IssmDouble  xyz_list[NUMVERTICES][3];
     8271        IssmDouble  xyz_list_segment[2][3];
     8272        IssmDouble  z_list[NUMVERTICES];
     8273        IssmDouble  slope[2];
     8274        IssmDouble  slope2,constant_part;
     8275        IssmDouble  rho_ice,gravity,n,B;
     8276        IssmDouble  ub,vb,z_g,surface,thickness;
     8277        GaussPenta* gauss=NULL;
     8278
     8279        /*Fetch number of nodes and dof for this finite element*/
     8280        int numnodes = this->NumberofNodes(); _assert_(numnodes==6);
     8281        int numdof   = numnodes*NDOF2;
    82788282
    82798283        /*Initialize Element vector*/
    8280         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     8284        ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters);
    82818285
    82828286        /*Retrieve all inputs and parameters*/
    … …  
    82948298        /*Loop on the three segments*/
    82958299        for(i=0;i<3;i++){
    8296                 node0=i;
    8297                 node1=i+3;
     8300
     8301                node0=i; node1=i+3;
    82988302
    82998303                for(j=0;j<3;j++){
    … …  
    83028306                }
    83038307
    8304                 connectivity[0]=vertices[node0]->Connectivity();
    8305                 connectivity[1]=vertices[node1]->Connectivity();
     8308                connectivity[0]=(IssmDouble)vertices[node0]->Connectivity();
     8309                connectivity[1]=(IssmDouble)vertices[node1]->Connectivity();
    83068310
    83078311                /*Loop on the Gauss points: */
    … …  
    83168320
    83178321                        slope2=pow(slope[0],2)+pow(slope[1],2);
    8318                         constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
     8322                        constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.));
    83198323
    83208324                        PentaRef::GetInputValue(&z_g,&z_list[0],gauss);
    83218325                        GetSegmentJacobianDeterminant(&Jdet,&xyz_list_segment[0][0],gauss);
    83228326
    8323                         if (IsOnSurface()){
    8324                                 for(j=0;j<NDOF2;j++) pe->values[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight/(IssmDouble)connectivity[1];
     8327                        if(IsOnSurface()){
     8328                                pe->values[2*node1+0]+=constant_part*pow((surface-z_g)/B,n)*slope[0]*Jdet*gauss->weight/connectivity[1];
     8329                                pe->values[2*node1+1]+=constant_part*pow((surface-z_g)/B,n)*slope[1]*Jdet*gauss->weight/connectivity[1];
    83258330                        }
    8326                         else{//connectivity is too large, should take only half on it
    8327                                 for(j=0;j<NDOF2;j++) pe->values[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight*2/(IssmDouble)connectivity[1];
     8331                        else{/*connectivity is too large, should take only half on it*/
     8332                                pe->values[2*node1+0]+=constant_part*pow((surface-z_g)/B,n)*slope[0]*Jdet*gauss->weight*2./connectivity[1];
     8333                                pe->values[2*node1+1]+=constant_part*pow((surface-z_g)/B,n)*slope[1]*Jdet*gauss->weight*2./connectivity[1];
    83288334                        }
    83298335                }
    83308336                delete gauss;
    83318337
    8332                 //Deal with lower surface
    8333                 if (IsOnBed()){
    8334                         constant_part=-1.58*pow((IssmDouble)10.0,-(IssmDouble)10.0)*rho_ice*gravity*thickness;
     8338                /*Deal with lower surface*/
     8339                if(IsOnBed()){
     8340                        constant_part=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness;
    83358341                        ub=constant_part*slope[0];
    83368342                        vb=constant_part*slope[1];
    83378343
    8338                         pe->values[2*node0]+=ub/(IssmDouble)connectivity[0];
    8339                         pe->values[2*node0+1]+=vb/(IssmDouble)connectivity[0];
     8344                        pe->values[2*node0+0]+=ub/connectivity[0];
     8345                        pe->values[2*node0+1]+=vb/connectivity[0];
    83408346                }
    83418347        }

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

    r15704 r15711  
    31063106
    31073107        /*Intermediaries*/
    3108         const int numdof=NUMVERTICES*NDOF2;
    3109         int    i,connectivity;
     3108        IssmDouble connectivity;
     3109
     3110        /*Fetch number of nodes and dof for this finite element*/
     3111        int numnodes = this->NumberofNodes(); _assert_(numnodes==3);
     3112        int numdof   = numnodes*NDOF2;
    31103113
    31113114        /*Initialize Element matrix*/
    3112         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     3115        ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    31133116
    31143117        /*Create Element matrix*/
    3115         for(i=0;i<NUMVERTICES;i++){
    3116                 connectivity=vertices[i]->Connectivity();
    3117                 Ke->values[(2*i)*numdof  +(2*i)  ]=1/(IssmDouble)connectivity;
    3118                 Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(IssmDouble)connectivity;
     3118        for(int i=0;i<numnodes;i++){
     3119                connectivity=(IssmDouble)vertices[i]->Connectivity();
     3120                Ke->values[(2*i+0)*numdof+(2*i+0)]=1./connectivity;
     3121                Ke->values[(2*i+1)*numdof+(2*i+1)]=1./connectivity;
    31193122        }
    31203123
    … …  
    32763279
    32773280        /*Intermediaries */
    3278         int        i,connectivity;
    3279         IssmDouble     constant_part,ub,vb;
    3280         IssmDouble     rho_ice,gravity,n,B;
    3281         IssmDouble     slope2,thickness;
    3282         IssmDouble     slope[2];
    3283         GaussTria* gauss=NULL;
     3281        IssmDouble constant_part,ub,vb;
     3282        IssmDouble rho_ice,gravity,n,B;
     3283        IssmDouble slope2,thickness,connectivity;
     3284        IssmDouble slope[2];
     3285
     3286        /*Fetch number of nodes and dof for this finite element*/
     3287        int numnodes = this->NumberofNodes(); _assert_(numnodes==3);
     3288        int numdof   = numnodes*NDOF2;
    32843289
    32853290        /*Initialize Element vector*/
    3286         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3291        ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters);
    32873292
    32883293        /*Retrieve all inputs and parameters*/
    … …  
    32963301
    32973302        /*Spawn 3 sing elements: */
    3298         gauss=new GaussTria();
    3299         for(i=0;i<NUMVERTICES;i++){
     3303        GaussTria* gauss=new GaussTria();
     3304        for(int i=0;i<numnodes;i++){
    33003305
    33013306                gauss->GaussVertex(i);
    33023307
    3303                 connectivity=vertices[i]->Connectivity();
     3308                connectivity=(IssmDouble)vertices[i]->Connectivity();
    33043309
    33053310                thickness_input->GetInputValue(&thickness,gauss);
    … …  
    33083313                slope2=pow(slope[0],2)+pow(slope[1],2);
    33093314
    3310                 constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
    3311 
    3312                 ub=-1.58*pow((IssmDouble)10.0,(IssmDouble)-10.0)*rho_ice*gravity*thickness*slope[0];
    3313                 vb=-1.58*pow((IssmDouble)10.0,(IssmDouble)-10.0)*rho_ice*gravity*thickness*slope[1];
    3314 
    3315                 pe->values[2*i]  =(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(IssmDouble)connectivity;
    3316                 pe->values[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(IssmDouble)connectivity;
     3315                constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.));
     3316
     3317                ub=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness*slope[0];
     3318                vb=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness*slope[1];
     3319
     3320                pe->values[2*i+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity;
     3321                pe->values[2*i+1]=(vb-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/connectivity;
    33173322        }
    33183323
Note: See TracChangeset for help on using the changeset viewer.