Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15710) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15711) @@ -6981,58 +6981,58 @@ ElementMatrix* Penta::CreateKMatrixDiagnosticSIA(void){ - /*Constants*/ - const int numdof=NDOF2*NUMVERTICES; - /*Intermediaries*/ - int connectivity[2]; - int i,i0,i1,j0,j1; - IssmDouble one0,one1; + IssmDouble connectivity[2]; + IssmDouble one0,one1; + int i,i0,i1,j0,j1; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); _assert_(numnodes==6); + int numdof = numnodes*NDOF2; /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Spawn 3 beam elements: */ + ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); + + /*3 vertical edges*/ for(i=0;i<3;i++){ + /*2 dofs of first node*/ - i0=2*i; - i1=2*i+1; + i0=2*i; i1=2*i+1; /*2 dofs of second node*/ - j0=2*(i+3); - j1=2*(i+3)+1; + j0=2*(i+3); j1=2*(i+3)+1; /*Find connectivity for the two nodes*/ - connectivity[0]=vertices[i]->Connectivity(); - connectivity[1]=vertices[i+3]->Connectivity(); - one0=1/(IssmDouble)connectivity[0]; - one1=1/(IssmDouble)connectivity[1]; + connectivity[0]=(IssmDouble)vertices[i]->Connectivity(); + connectivity[1]=(IssmDouble)vertices[i+3]->Connectivity(); + one0=1./connectivity[0]; + one1=1./connectivity[1]; /*Create matrix for these two nodes*/ if (IsOnBed() && IsOnSurface()){ - Ke->values[i0*numdof+i0]=one0; - Ke->values[i1*numdof+i1]=one0; - Ke->values[j0*numdof+i0]=-one1; - Ke->values[j0*numdof+j0]=one1; - Ke->values[j1*numdof+i1]=-one1; - Ke->values[j1*numdof+j1]=one1; + Ke->values[i0*numdof+i0] = +one0; + Ke->values[i1*numdof+i1] = +one0; + Ke->values[j0*numdof+i0] = -one1; + Ke->values[j0*numdof+j0] = +one1; + Ke->values[j1*numdof+i1] = -one1; + Ke->values[j1*numdof+j1] = +one1; } else if (IsOnBed()){ - Ke->values[i0*numdof+i0]=one0; - Ke->values[i1*numdof+i1]=one0; - Ke->values[j0*numdof+i0]=-2*one1; - Ke->values[j0*numdof+j0]=2*one1; - Ke->values[j1*numdof+i1]=-2*one1; - Ke->values[j1*numdof+j1]=2*one1; + Ke->values[i0*numdof+i0] = one0; + Ke->values[i1*numdof+i1] = one0; + Ke->values[j0*numdof+i0] = -2.*one1; + Ke->values[j0*numdof+j0] = +2.*one1; + Ke->values[j1*numdof+i1] = -2.*one1; + Ke->values[j1*numdof+j1] = +2.*one1; } else if (IsOnSurface()){ - Ke->values[j0*numdof+i0]=-one1; - Ke->values[j0*numdof+j0]=one1; - Ke->values[j1*numdof+i1]=-one1; - Ke->values[j1*numdof+j1]=one1; + Ke->values[j0*numdof+i0] = -one1; + Ke->values[j0*numdof+j0] = +one1; + Ke->values[j1*numdof+i1] = -one1; + Ke->values[j1*numdof+j1] = +one1; } else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity - Ke->values[j0*numdof+i0]=-2*one1; - Ke->values[j0*numdof+j0]=2*one1; - Ke->values[j1*numdof+i1]=-2*one1; - Ke->values[j1*numdof+j1]=2*one1; + Ke->values[j0*numdof+i0] = -2.*one1; + Ke->values[j0*numdof+j0] = +2.*one1; + Ke->values[j1*numdof+i1] = -2.*one1; + Ke->values[j1*numdof+j1] = +2.*one1; } } @@ -8264,19 +8264,23 @@ /*Intermediaries*/ - int i,j; - int node0,node1; - int connectivity[2]; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_segment[2][3]; - IssmDouble z_list[NUMVERTICES]; - IssmDouble slope[2]; - IssmDouble slope2,constant_part; - IssmDouble rho_ice,gravity,n,B; - IssmDouble ub,vb,z_g,surface,thickness; - GaussPenta* gauss=NULL; + int i,j; + int node0,node1; + IssmDouble connectivity[2]; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_segment[2][3]; + IssmDouble z_list[NUMVERTICES]; + IssmDouble slope[2]; + IssmDouble slope2,constant_part; + IssmDouble rho_ice,gravity,n,B; + IssmDouble ub,vb,z_g,surface,thickness; + GaussPenta* gauss=NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); _assert_(numnodes==6); + int numdof = numnodes*NDOF2; /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters); /*Retrieve all inputs and parameters*/ @@ -8294,6 +8298,6 @@ /*Loop on the three segments*/ for(i=0;i<3;i++){ - node0=i; - node1=i+3; + + node0=i; node1=i+3; for(j=0;j<3;j++){ @@ -8302,6 +8306,6 @@ } - connectivity[0]=vertices[node0]->Connectivity(); - connectivity[1]=vertices[node1]->Connectivity(); + connectivity[0]=(IssmDouble)vertices[node0]->Connectivity(); + connectivity[1]=(IssmDouble)vertices[node1]->Connectivity(); /*Loop on the Gauss points: */ @@ -8316,26 +8320,28 @@ slope2=pow(slope[0],2)+pow(slope[1],2); - constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2)); + constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.)); PentaRef::GetInputValue(&z_g,&z_list[0],gauss); GetSegmentJacobianDeterminant(&Jdet,&xyz_list_segment[0][0],gauss); - if (IsOnSurface()){ - for(j=0;jvalues[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight/(IssmDouble)connectivity[1]; + if(IsOnSurface()){ + pe->values[2*node1+0]+=constant_part*pow((surface-z_g)/B,n)*slope[0]*Jdet*gauss->weight/connectivity[1]; + pe->values[2*node1+1]+=constant_part*pow((surface-z_g)/B,n)*slope[1]*Jdet*gauss->weight/connectivity[1]; } - else{//connectivity is too large, should take only half on it - for(j=0;jvalues[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight*2/(IssmDouble)connectivity[1]; + else{/*connectivity is too large, should take only half on it*/ + pe->values[2*node1+0]+=constant_part*pow((surface-z_g)/B,n)*slope[0]*Jdet*gauss->weight*2./connectivity[1]; + pe->values[2*node1+1]+=constant_part*pow((surface-z_g)/B,n)*slope[1]*Jdet*gauss->weight*2./connectivity[1]; } } delete gauss; - //Deal with lower surface - if (IsOnBed()){ - constant_part=-1.58*pow((IssmDouble)10.0,-(IssmDouble)10.0)*rho_ice*gravity*thickness; + /*Deal with lower surface*/ + if(IsOnBed()){ + constant_part=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness; ub=constant_part*slope[0]; vb=constant_part*slope[1]; - pe->values[2*node0]+=ub/(IssmDouble)connectivity[0]; - pe->values[2*node0+1]+=vb/(IssmDouble)connectivity[0]; + pe->values[2*node0+0]+=ub/connectivity[0]; + pe->values[2*node0+1]+=vb/connectivity[0]; } } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15710) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15711) @@ -3106,15 +3106,18 @@ /*Intermediaries*/ - const int numdof=NUMVERTICES*NDOF2; - int i,connectivity; + IssmDouble connectivity; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); _assert_(numnodes==3); + int numdof = numnodes*NDOF2; /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); + ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); /*Create Element matrix*/ - for(i=0;iConnectivity(); - Ke->values[(2*i)*numdof +(2*i) ]=1/(IssmDouble)connectivity; - Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(IssmDouble)connectivity; + for(int i=0;iConnectivity(); + Ke->values[(2*i+0)*numdof+(2*i+0)]=1./connectivity; + Ke->values[(2*i+1)*numdof+(2*i+1)]=1./connectivity; } @@ -3276,13 +3279,15 @@ /*Intermediaries */ - int i,connectivity; - IssmDouble constant_part,ub,vb; - IssmDouble rho_ice,gravity,n,B; - IssmDouble slope2,thickness; - IssmDouble slope[2]; - GaussTria* gauss=NULL; + IssmDouble constant_part,ub,vb; + IssmDouble rho_ice,gravity,n,B; + IssmDouble slope2,thickness,connectivity; + IssmDouble slope[2]; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); _assert_(numnodes==3); + int numdof = numnodes*NDOF2; /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters); /*Retrieve all inputs and parameters*/ @@ -3296,10 +3301,10 @@ /*Spawn 3 sing elements: */ - gauss=new GaussTria(); - for(i=0;iGaussVertex(i); - connectivity=vertices[i]->Connectivity(); + connectivity=(IssmDouble)vertices[i]->Connectivity(); thickness_input->GetInputValue(&thickness,gauss); @@ -3308,11 +3313,11 @@ slope2=pow(slope[0],2)+pow(slope[1],2); - constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2)); - - ub=-1.58*pow((IssmDouble)10.0,(IssmDouble)-10.0)*rho_ice*gravity*thickness*slope[0]; - vb=-1.58*pow((IssmDouble)10.0,(IssmDouble)-10.0)*rho_ice*gravity*thickness*slope[1]; - - 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; - 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; + constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.)); + + ub=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness*slope[0]; + vb=-1.58*pow(10.,-10.)*rho_ice*gravity*thickness*slope[1]; + + 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; + 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; }