Changeset 16036


Ignore:
Timestamp:
08/30/13 10:31:53 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: balance velocitis now use observed directions if provided

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

Legend:

Unmodified
Added
Removed

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

    r16034 r16036  
    32633263        for(int i=0;i<numnodes;i++){
    32643264
    3265                 if(!flags[this->nodes[i]->Sid()]){
     3265                if(!flags[this->nodes[i]->Lid()]){
    32663266
    32673267                        /*flag current node so that no other element processes it*/
    3268                         flags[this->nodes[i]->Sid()]=true;
     3268                        flags[this->nodes[i]->Lid()]=true;
    32693269
    32703270                        /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/

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

    r16021 r16036  
    131131        for(int i=0;i<numnodes;i++){
    132132
    133                 if(!flags[this->nodes[i]->Sid()]){
     133                if(!flags[this->nodes[i]->Lid()]){
    134134
    135135                        /*flag current node so that no other element processes it*/
    136                         flags[this->nodes[i]->Sid()]=true;
     136                        flags[this->nodes[i]->Lid()]=true;
    137137
    138138                        /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
    … …  
    72307230        IssmDouble h,gamma,thickness;
    72317231        IssmDouble hnx,hny,dhnx[2],dhny[2];
    7232         IssmDouble vel,vx,vy,dvxdx,dvydy;
    7233         IssmDouble D[2][2];
    7234         int        stabilization=0;
    72357232
    72367233        /*Fetch number of nodes and dof for this finite element*/
    … …  
    72507247        /*Retrieve all Inputs and parameters: */
    72517248        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7252         Input* H_input=inputs->GetInput(ThicknessEnum); _assert_(H_input);
    7253 
    7254         /*Get vector N for all nodes*/
     7249        Input* H_input =inputs->GetInput(ThicknessEnum); _assert_(H_input);
     7250        h=sqrt(2.*this->GetArea());
     7251
     7252        /*Get vector N for all nodes and build HNx and HNy*/
    72557253        GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
    72567254        GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
     7255        GetInputListOnNodes(H,ThicknessEnum);
    72577256        for(int i=0;i<numnodes;i++){
    72587257                IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
    7259                 Nx[i] = -Nx[i]/norm;
    7260                 Ny[i] = -Ny[i]/norm;
    7261         }
    7262 
    7263         /*Build HNx and HNy*/
    7264         GetInputListOnNodes(H,ThicknessEnum);
    7265         for(int i=0;i<numnodes;i++){
    7266                 HNx[i]=H[i]*Nx[i];
    7267                 HNy[i]=H[i]*Ny[i];
    7268         }
    7269 
    7270         h=sqrt(2.*this->GetArea());
     7258                HNx[i] = -H[i]*Nx[i]/norm;
     7259                HNy[i] = -H[i]*Ny[i]/norm;
     7260        }
    72717261
    72727262        /*Start looping on the number of gaussian points:*/
    … …  
    72967286                                                        );
    72977287                        }
    7298                 }
    7299 
    7300                 vx=hnx;
    7301                 vy=hny;
    7302                 if(stabilization==2){
    7303                         /*Streamline upwinding*/
    7304                         vel=sqrt(vx*vx+vy*vy)+1.e-8;
    7305                         D[0][0]=h/(2*vel)*vx*vx;
    7306                         D[1][0]=h/(2*vel)*vy*vx;
    7307                         D[0][1]=h/(2*vel)*vx*vy;
    7308                         D[1][1]=h/(2*vel)*vy*vy;
    7309                 }
    7310                 else if(stabilization==1){
    7311                         /*MacAyeal*/
    7312                         D[0][0]=h/2.0*fabs(vx);
    7313                         D[0][1]=0.;
    7314                         D[1][0]=0.;
    7315                         D[1][1]=h/2.0*fabs(vy);
    7316                 }
    7317                 if(stabilization==1 || stabilization==2){
    7318                         GetBprimeMasstransport(B,&xyz_list[0][0],gauss);
    7319                         D[0][0]=gauss->weight*Jdet*D[0][0];
    7320                         D[1][0]=gauss->weight*Jdet*D[1][0];
    7321                         D[0][1]=gauss->weight*Jdet*D[0][1];
    7322                         D[1][1]=gauss->weight*Jdet*D[1][1];
    7323                         TripleMultiply(B,2,numnodes,1,
    7324                                                 &D[0][0],2,2,0,
    7325                                                 B,2,numnodes,0,
    7326                                                 &Ke->values[0],1);
    73277288                }
    73287289        }
    … …  
    73487309        IssmDouble xyz_list[NUMVERTICES][3];
    73497310        IssmDouble D[2][2];
    7350         IssmDouble l=12.;
     7311        IssmDouble l=8.;
    73517312
    73527313        /*Fetch number of nodes and dof for this finite element*/
    … …  
    75187479        Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
    75197480        Input* H_input=inputs->GetInput(ThicknessEnum); _assert_(H_input);
     7481        h=sqrt(2.*this->GetArea());
    75207482
    75217483        /*Get vector N for all nodes*/
    75227484        GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
    75237485        GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
     7486        GetInputListOnNodes(H,ThicknessEnum);
    75247487        for(int i=0;i<numnodes;i++){
    75257488                IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
    7526                 Nx[i] = -Nx[i]/norm;
    7527                 Ny[i] = -Ny[i]/norm;
    7528         }
    7529 
    7530         /*Build HNx and HNy*/
    7531         GetInputListOnNodes(H,ThicknessEnum);
    7532         for(int i=0;i<numnodes;i++){
    7533                 HNx[i]=H[i]*Nx[i];
    7534                 HNy[i]=H[i]*Ny[i];
    7535         }
    7536 
    7537         h=sqrt(2.*this->GetArea());
     7489                Nx[i] = -H[i]*Nx[i]/norm;
     7490                Ny[i] = -H[i]*Ny[i]/norm;
     7491        }
    75387492
    75397493        /* Start  looping on the number of gaussian points: */
    … …  
    75867540        IssmDouble Jdet;
    75877541        IssmDouble thickness,slope[2];
    7588         IssmDouble taud_x;
     7542        IssmDouble taud_x,norms,normv,vx,vy;
    75897543
    75907544        /*Fetch number of nodes and dof for this finite element*/
    … …  
    75947548        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    75957549        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     7550        IssmDouble*    Vx     = xNew<IssmDouble>(numnodes);
     7551        IssmDouble*    Vy     = xNew<IssmDouble>(numnodes);
    75967552
    75977553        /*Retrieve all inputs and parameters*/
    … …  
    75997555        Input* H_input       = inputs->GetInput(ThicknessEnum); _assert_(H_input);
    76007556        Input* surface_input = inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
     7557        Input* vx_input      = inputs->GetInput(VxEnum);
     7558        Input* vy_input      = inputs->GetInput(VyEnum);
    76017559
    76027560        /* Start  looping on the number of gaussian points: */
    … …  
    76087566                H_input->GetInputValue(&thickness,gauss);
    76097567                surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
     7568                if(vx_input && vy_input){
     7569                        vx_input->GetInputValue(&vx,gauss);
     7570                        vy_input->GetInputValue(&vy,gauss);
     7571                        norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
     7572                        normv = sqrt(vx*vx + vy*vy);
     7573                        if(normv>15./(365.*24.*3600.)) slope[0] = -vx/normv*norms;
     7574                }
    76107575                taud_x = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[0];
    7611                 //taud_x = slope[0];
    76127576
    76137577                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    … …  
    76307594        IssmDouble Jdet;
    76317595        IssmDouble thickness,slope[2];
    7632         IssmDouble taud_y;
     7596        IssmDouble taud_y,norms,normv,vx,vy;
    76337597
    76347598        /*Fetch number of nodes and dof for this finite element*/
    … …  
    76377601        /*Initialize Element vector and other vectors*/
    76387602        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    7639         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     7603        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    76407604
    76417605        /*Retrieve all inputs and parameters*/
    … …  
    76437607        Input* H_input       = inputs->GetInput(ThicknessEnum); _assert_(H_input);
    76447608        Input* surface_input = inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
     7609        Input* vx_input      = inputs->GetInput(VxEnum);
     7610        Input* vy_input      = inputs->GetInput(VyEnum);
    76457611
    76467612        /* Start  looping on the number of gaussian points: */
    … …  
    76527618                H_input->GetInputValue(&thickness,gauss);
    76537619                surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
     7620                if(vx_input && vy_input){
     7621                        vx_input->GetInputValue(&vx,gauss);
     7622                        vy_input->GetInputValue(&vy,gauss);
     7623                        norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
     7624                        normv = sqrt(vx*vx + vy*vy);
     7625                        if(normv>15./(365.*24.*3600.)) slope[1] = -vy/normv*norms;
     7626                }
    76547627                taud_y = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[1];
    7655                 //taud_y = slope[1];
    76567628
    76577629                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
Note: See TracChangeset for help on using the changeset viewer.