Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15711) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15712) @@ -8441,6 +8441,6 @@ /*Clean up and return*/ + xDelete(basis); delete gauss; - xDelete(basis); return pe; } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15711) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15712) @@ -3144,10 +3144,10 @@ /*Intermediaries */ - int i,j; - IssmDouble driving_stress_baseline,thickness; - IssmDouble Jdet; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble slope[2]; - IssmDouble icefrontlevel[3]; + int i; + IssmDouble driving_stress_baseline,thickness; + IssmDouble Jdet; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble slope[2]; + IssmDouble icefrontlevel[3]; /*Fetch number of nodes and dof for this finite element*/ @@ -5008,16 +5008,16 @@ /*Intermediaries */ - int i,resp; - int *responses=NULL; - int num_responses; - IssmDouble Jdet; - IssmDouble obs_velocity_mag,velocity_mag; - IssmDouble dux,duy; - IssmDouble epsvel=2.220446049250313e-16; - IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ - IssmDouble scalex=0,scaley=0,scale=0,S=0; - IssmDouble vx,vy,vxobs,vyobs,weight; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[3]; + int i,resp; + int *responses=NULL; + int num_responses; + IssmDouble Jdet; + IssmDouble obs_velocity_mag,velocity_mag; + IssmDouble dux,duy; + IssmDouble epsvel=2.220446049250313e-16; + IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ + IssmDouble scalex=0,scaley=0,scale=0,S=0; + IssmDouble vx,vy,vxobs,vyobs,weight; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[3]; GaussTria* gauss=NULL; @@ -5639,24 +5639,26 @@ ElementMatrix* Tria::CreateKMatrixMelting(void){ - /*Constants*/ - const int numdof=NUMVERTICES*NDOF1; - /*Intermediaries */ IssmDouble heatcapacity,latentheat; IssmDouble Jdet,D_scalar; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[3]; - GaussTria *gauss=NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF1; + + /*Initialize Element vector and vectors*/ + ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); + IssmDouble* basis = xNew(numnodes); /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - latentheat=matpar->GetLatentHeat(); - heatcapacity=matpar->GetHeatCapacity(); + latentheat = matpar->GetLatentHeat(); + heatcapacity = matpar->GetHeatCapacity(); /* Start looping on the number of gauss (nodes on the bedrock) */ - gauss=new GaussTria(2); + GaussTria* gauss=new GaussTria(2); for(int ig=gauss->begin();igend();ig++){ @@ -5668,11 +5670,12 @@ D_scalar=latentheat/heatcapacity*gauss->weight*Jdet; - TripleMultiply(&basis[0],numdof,1,0, + TripleMultiply(basis,numdof,1,0, &D_scalar,1,1,0, - &basis[0],1,numdof,0, + basis,1,numdof,0, &Ke->values[0],1); } /*Clean up and return*/ + xDelete(basis); delete gauss; return Ke; @@ -5764,30 +5767,26 @@ ElementMatrix* Tria::CreateKMatrixHydrologyShreve(void){ - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - -/*Intermediaries */ + /*Intermediaries */ IssmDouble diffusivity; - IssmDouble Jdettria,DL_scalar,dt,h; + IssmDouble Jdettria,D_scalar,dt,h; IssmDouble vx,vy,vel,dvxdx,dvydy; IssmDouble dvx[2],dvy[2]; - IssmDouble v_gauss[2]={0.0}; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[NUMVERTICES]; - IssmDouble B[2][NUMVERTICES]; - IssmDouble Bprime[2][NUMVERTICES]; - IssmDouble K[2][2] ={0.0}; - IssmDouble KDL[2][2] ={0.0}; - IssmDouble DL[2][2] ={0.0}; - IssmDouble DLprime[2][2] ={0.0}; - GaussTria *gauss=NULL; - -/*Skip if water or ice shelf element*/ + + /*Skip if water or ice shelf element*/ if(IsOnWater() | IsFloating()) return NULL; -/*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - -/*Create water velocity vx and vy from current inputs*/ + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF1; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); + IssmDouble* basis = xNew(numnodes); + IssmDouble* B = xNew(2*numdof); + IssmDouble* Bprime = xNew(2*numdof); + IssmDouble D[2][2]; + + /*Create water velocity vx and vy from current inputs*/ CreateHydrologyWaterVelocityInput(); @@ -5800,6 +5799,6 @@ h=sqrt(2*this->GetArea()); -/* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); + /* Start looping on the number of gaussian points: */ + GaussTria* gauss=new GaussTria(2); for(int ig=gauss->begin();igend();ig++){ @@ -5807,5 +5806,5 @@ GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); - GetNodalFunctions(&basis[0],gauss); + GetNodalFunctions(basis,gauss); vx_input->GetInputValue(&vx,gauss); @@ -5814,51 +5813,50 @@ vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); - DL_scalar=gauss->weight*Jdettria; - - TripleMultiply(&basis[0],1,numdof,1, - &DL_scalar,1,1,0, - &basis[0],1,numdof,0, - &Ke->values[0],1); - - GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss); - GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); + D_scalar=gauss->weight*Jdettria; + + TripleMultiply(basis,1,numdof,1, + &D_scalar,1,1,0, + basis,1,numdof,0, + Ke->values,1); + + GetBPrognostic(B,&xyz_list[0][0], gauss); + GetBprimePrognostic(Bprime,&xyz_list[0][0], gauss); dvxdx=dvx[0]; dvydy=dvy[1]; - DL_scalar=dt*gauss->weight*Jdettria; - - DL[0][0]=DL_scalar*dvxdx; - DL[1][1]=DL_scalar*dvydy; - DLprime[0][0]=DL_scalar*vx; - DLprime[1][1]=DL_scalar*vy; - - TripleMultiply( &B[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &B[0][0],2,numdof,0, - &Ke->values[0],1); - - TripleMultiply( &B[0][0],2,numdof,1, - &DLprime[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke->values[0],1); + D_scalar=dt*gauss->weight*Jdettria; + + D[0][0]=D_scalar*dvxdx; + D[0][1]=0.; + D[1][1]=D_scalar*dvydy; + D[1][1]=0.; + TripleMultiply(B,2,numdof,1, + &D[0][0],2,2,0, + B,2,numdof,0, + &Ke->values[0],1); + + D[0][0]=D_scalar*vx; + D[1][1]=D_scalar*vy; + TripleMultiply(B,2,numdof,1, + &D[0][0],2,2,0, + Bprime,2,numdof,0, + &Ke->values[0],1); /*Artificial diffusivity*/ vel=sqrt(vx*vx+vy*vy); - K[0][0]=diffusivity*h/(2*vel)*vx*vx; - K[1][0]=diffusivity*h/(2*vel)*vy*vx; - K[0][1]=diffusivity*h/(2*vel)*vx*vy; - K[1][1]=diffusivity*h/(2*vel)*vy*vy; - KDL[0][0]=DL_scalar*K[0][0]; - KDL[1][0]=DL_scalar*K[1][0]; - KDL[0][1]=DL_scalar*K[0][1]; - KDL[1][1]=DL_scalar*K[1][1]; - - TripleMultiply( &Bprime[0][0],2,numdof,1, - &KDL[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke->values[0],1); - } - -/*Clean up and return*/ + D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx; + D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx; + D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy; + D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy; + TripleMultiply(Bprime,2,numdof,1, + &D[0][0],2,2,0, + Bprime,2,numdof,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + xDelete(basis); + xDelete(B); + xDelete(Bprime); delete gauss; return Ke; @@ -5867,7 +5865,4 @@ /*FUNCTION Tria::CreateKMatrixHydrologyDCInefficient{{{*/ ElementMatrix* Tria::CreateKMatrixHydrologyDCInefficient(void){ - - /*constants: */ - const int numdof=NDOF1*NUMVERTICES; /* Intermediaries */ @@ -5876,11 +5871,14 @@ IssmDouble sediment_storing; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B[2][numdof]; - IssmDouble basis[NUMVERTICES]; - IssmDouble D[2][2]; - GaussTria *gauss = NULL; - - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters); + IssmDouble* B = xNew(5*numdof); + IssmDouble* basis = xNew(numnodes); + IssmDouble D[2][2]; /*Retrieve all inputs and parameters*/ @@ -5891,5 +5889,5 @@ /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); + GaussTria* gauss=new GaussTria(2); for(int ig=gauss->begin();igend();ig++){ @@ -5903,8 +5901,8 @@ D[0][0]=D_scalar; D[0][1]=0.; D[1][0]=0.; D[1][1]=D_scalar; - GetBHydro(&B[0][0],&xyz_list[0][0],gauss); - TripleMultiply(&B[0][0],2,numdof,1, + GetBHydro(B,&xyz_list[0][0],gauss); + TripleMultiply(B,2,numdof,1, &D[0][0],2,2,0, - &B[0][0],2,numdof,0, + B,2,numdof,0, &Ke->values[0],1); @@ -5914,11 +5912,13 @@ D_scalar=sediment_storing*gauss->weight*Jdet; - TripleMultiply(&basis[0],numdof,1,0, + TripleMultiply(basis,numdof,1,0, &D_scalar,1,1,0, - &basis[0],1,numdof,0, + basis,1,numdof,0, &Ke->values[0],1); } } /*Clean up and return*/ + xDelete(basis); + xDelete(B); delete gauss; return Ke; @@ -5969,7 +5969,7 @@ GetBHydro(&B[0][0],&xyz_list[0][0],gauss); TripleMultiply(&B[0][0],2,numdof,1, - &D[0][0],2,2,0, - &B[0][0],2,numdof,0, - &Ke->values[0],1); + &D[0][0],2,2,0, + &B[0][0],2,numdof,0, + &Ke->values[0],1); /*Transient*/ @@ -5992,21 +5992,21 @@ ElementVector* Tria::CreatePVectorHydrologyShreve(void){ - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - /*Intermediaries */ - int i,j; + int i; IssmDouble Jdettria,dt; IssmDouble basal_melting_g; IssmDouble old_watercolumn_g; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[numdof]; - GaussTria* gauss=NULL; /*Skip if water or ice shelf element*/ if(IsOnWater() | IsFloating()) return NULL; + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF1; + /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters); + IssmDouble* basis = xNewZeroInit(numnodes); /*Retrieve all inputs and parameters*/ @@ -6018,5 +6018,5 @@ /*Initialize basal_melting_correction_g to 0, do not forget!:*/ /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); + GaussTria* gauss=new GaussTria(2); for(int ig=gauss->begin();igend();ig++){ @@ -6029,9 +6029,14 @@ old_watercolumn_input->GetInputValue(&old_watercolumn_g,gauss); - if(reCast(dt))for(i=0;ivalues[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i]; - else for(i=0;ivalues[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i]; + if(reCast(dt)){ + for(i=0;ivalues[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i]; + } + else{ + for(i=0;ivalues[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i]; + } } /*Clean up and return*/ + xDelete(basis); delete gauss; return pe; @@ -6484,5 +6489,4 @@ IssmDouble vel,vx,vy,dvxdx,dvydy; IssmDouble dvx[2],dvy[2]; - IssmDouble v_gauss[2]={0.0}; IssmDouble xyz_list[NUMVERTICES][3]; @@ -6589,8 +6593,8 @@ /*Clean up and return*/ - delete gauss; xDelete(basis); xDelete(B); xDelete(Bprime); + delete gauss; return Ke; }