Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16035) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16036) @@ -3263,8 +3263,8 @@ for(int i=0;inodes[i]->Sid()]){ + if(!flags[this->nodes[i]->Lid()]){ /*flag current node so that no other element processes it*/ - flags[this->nodes[i]->Sid()]=true; + flags[this->nodes[i]->Lid()]=true; /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16035) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16036) @@ -131,8 +131,8 @@ for(int i=0;inodes[i]->Sid()]){ + if(!flags[this->nodes[i]->Lid()]){ /*flag current node so that no other element processes it*/ - flags[this->nodes[i]->Sid()]=true; + flags[this->nodes[i]->Lid()]=true; /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/ @@ -7230,7 +7230,4 @@ IssmDouble h,gamma,thickness; IssmDouble hnx,hny,dhnx[2],dhny[2]; - IssmDouble vel,vx,vy,dvxdx,dvydy; - IssmDouble D[2][2]; - int stabilization=0; /*Fetch number of nodes and dof for this finite element*/ @@ -7250,23 +7247,16 @@ /*Retrieve all Inputs and parameters: */ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* H_input=inputs->GetInput(ThicknessEnum); _assert_(H_input); - - /*Get vector N for all nodes*/ + Input* H_input =inputs->GetInput(ThicknessEnum); _assert_(H_input); + h=sqrt(2.*this->GetArea()); + + /*Get vector N for all nodes and build HNx and HNy*/ GetInputListOnNodes(Nx,SurfaceSlopeXEnum); GetInputListOnNodes(Ny,SurfaceSlopeYEnum); + GetInputListOnNodes(H,ThicknessEnum); for(int i=0;iGetArea()); + HNx[i] = -H[i]*Nx[i]/norm; + HNy[i] = -H[i]*Ny[i]/norm; + } /*Start looping on the number of gaussian points:*/ @@ -7296,33 +7286,4 @@ ); } - } - - vx=hnx; - vy=hny; - if(stabilization==2){ - /*Streamline upwinding*/ - vel=sqrt(vx*vx+vy*vy)+1.e-8; - D[0][0]=h/(2*vel)*vx*vx; - D[1][0]=h/(2*vel)*vy*vx; - D[0][1]=h/(2*vel)*vx*vy; - D[1][1]=h/(2*vel)*vy*vy; - } - else if(stabilization==1){ - /*MacAyeal*/ - D[0][0]=h/2.0*fabs(vx); - D[0][1]=0.; - D[1][0]=0.; - D[1][1]=h/2.0*fabs(vy); - } - if(stabilization==1 || stabilization==2){ - GetBprimeMasstransport(B,&xyz_list[0][0],gauss); - D[0][0]=gauss->weight*Jdet*D[0][0]; - D[1][0]=gauss->weight*Jdet*D[1][0]; - D[0][1]=gauss->weight*Jdet*D[0][1]; - D[1][1]=gauss->weight*Jdet*D[1][1]; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - B,2,numnodes,0, - &Ke->values[0],1); } } @@ -7348,5 +7309,5 @@ IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble D[2][2]; - IssmDouble l=12.; + IssmDouble l=8.; /*Fetch number of nodes and dof for this finite element*/ @@ -7518,22 +7479,15 @@ Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); Input* H_input=inputs->GetInput(ThicknessEnum); _assert_(H_input); + h=sqrt(2.*this->GetArea()); /*Get vector N for all nodes*/ GetInputListOnNodes(Nx,SurfaceSlopeXEnum); GetInputListOnNodes(Ny,SurfaceSlopeYEnum); + GetInputListOnNodes(H,ThicknessEnum); for(int i=0;iGetArea()); + Nx[i] = -H[i]*Nx[i]/norm; + Ny[i] = -H[i]*Ny[i]/norm; + } /* Start looping on the number of gaussian points: */ @@ -7586,5 +7540,5 @@ IssmDouble Jdet; IssmDouble thickness,slope[2]; - IssmDouble taud_x; + IssmDouble taud_x,norms,normv,vx,vy; /*Fetch number of nodes and dof for this finite element*/ @@ -7594,4 +7548,6 @@ ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); IssmDouble* basis = xNew(numnodes); + IssmDouble* Vx = xNew(numnodes); + IssmDouble* Vy = xNew(numnodes); /*Retrieve all inputs and parameters*/ @@ -7599,4 +7555,6 @@ Input* H_input = inputs->GetInput(ThicknessEnum); _assert_(H_input); Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); + Input* vx_input = inputs->GetInput(VxEnum); + Input* vy_input = inputs->GetInput(VyEnum); /* Start looping on the number of gaussian points: */ @@ -7608,6 +7566,12 @@ H_input->GetInputValue(&thickness,gauss); surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); + if(vx_input && vy_input){ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]); + normv = sqrt(vx*vx + vy*vy); + if(normv>15./(365.*24.*3600.)) slope[0] = -vx/normv*norms; + } taud_x = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[0]; - //taud_x = slope[0]; GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); @@ -7630,5 +7594,5 @@ IssmDouble Jdet; IssmDouble thickness,slope[2]; - IssmDouble taud_y; + IssmDouble taud_y,norms,normv,vx,vy; /*Fetch number of nodes and dof for this finite element*/ @@ -7637,5 +7601,5 @@ /*Initialize Element vector and other vectors*/ ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); - IssmDouble* basis = xNew(numnodes); + IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ @@ -7643,4 +7607,6 @@ Input* H_input = inputs->GetInput(ThicknessEnum); _assert_(H_input); Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); + Input* vx_input = inputs->GetInput(VxEnum); + Input* vy_input = inputs->GetInput(VyEnum); /* Start looping on the number of gaussian points: */ @@ -7652,6 +7618,12 @@ H_input->GetInputValue(&thickness,gauss); surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); + if(vx_input && vy_input){ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]); + normv = sqrt(vx*vx + vy*vy); + if(normv>15./(365.*24.*3600.)) slope[1] = -vy/normv*norms; + } taud_y = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[1]; - //taud_y = slope[1]; GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);