Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15712) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15713) @@ -407,18 +407,16 @@ ElementVector* Tria::CreatePVectorSlope(void){ - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - /*Intermediaries */ - int i; - int analysis_type; + int i,analysis_type; IssmDouble Jdet; IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble slope[2]; - IssmDouble basis[3]; - GaussTria* gauss=NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); + IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ @@ -434,5 +432,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++){ @@ -440,17 +438,18 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis, gauss); + GetNodalFunctions(basis,gauss); slope_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); - if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){ - for(i=0;ivalues[i]+=Jdet*gauss->weight*slope[0]*basis[i]; - } - if ( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){ - for(i=0;ivalues[i]+=Jdet*gauss->weight*slope[1]*basis[i]; + if( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){ + for(i=0;ivalues[i]+=Jdet*gauss->weight*slope[0]*basis[i]; + } + if( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){ + for(i=0;ivalues[i]+=Jdet*gauss->weight*slope[1]*basis[i]; } } /*Clean up and return*/ + xDelete(basis); delete gauss; return pe; @@ -3157,5 +3156,4 @@ /*Initialize Element vector and vectors*/ ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,SSAApproximationEnum); - GaussTria* gauss = new GaussTria(2); IssmDouble* basis = xNew(numnodes); @@ -3168,4 +3166,5 @@ /* Start looping on the number of gaussian points: */ + GaussTria* gauss = new GaussTria(2); for(int ig=gauss->begin();igend();ig++){ @@ -4917,22 +4916,21 @@ ElementVector* Tria::CreatePVectorAdjointBalancethickness(void){ - /*Constants*/ - const int numdof=1*NUMVERTICES; - /*Intermediaries */ - int i,resp; - IssmDouble Jdet; - IssmDouble thickness,thicknessobs,weight; - int num_responses; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basis[3]; - IssmDouble dbasis[NDOF2][NUMVERTICES]; - IssmDouble dH[2]; - IssmDouble vx,vy,vel; - GaussTria *gauss = NULL; + int i,resp; + IssmDouble Jdet; + IssmDouble thickness,thicknessobs,weight; + int num_responses; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dH[2]; + IssmDouble vx,vy,vel; int *responses = NULL; - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Initialize Element vector and vectors*/ + ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); + IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(2*numnodes); /*Retrieve all inputs and parameters*/ @@ -4947,5 +4945,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++){ @@ -4954,5 +4952,5 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); GetNodalFunctions(basis, gauss); - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); + GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss); thickness_input->GetInputValue(&thickness, gauss); @@ -4965,10 +4963,10 @@ case ThicknessAbsMisfitEnum: weights_input->GetInputValue(&weight, gauss,resp); - for(i=0;ivalues[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i]; + for(i=0;ivalues[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i]; break; case ThicknessAbsGradientEnum: weights_input->GetInputValue(&weight, gauss,resp); - for(i=0;ivalues[i]+= - weight*dH[0]*dbasis[0][i]*Jdet*gauss->weight; - for(i=0;ivalues[i]+= - weight*dH[1]*dbasis[1][i]*Jdet*gauss->weight; + for(i=0;ivalues[i]+= - weight*dH[0]*dbasis[0*numnodes+i]*Jdet*gauss->weight; + for(i=0;ivalues[i]+= - weight*dH[1]*dbasis[1*numnodes+i]*Jdet*gauss->weight; break; case ThicknessAlongGradientEnum: @@ -4979,5 +4977,5 @@ vx = vx/(vel+1.e-9); vy = vy/(vel+1.e-9); - for(i=0;ivalues[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0][i]*vx+dbasis[1][i]*vy)*Jdet*gauss->weight; + for(i=0;ivalues[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0*numnodes+i]*vx+dbasis[1*numnodes+i]*vy)*Jdet*gauss->weight; break; case ThicknessAcrossGradientEnum: @@ -4988,5 +4986,5 @@ vx = vx/(vel+1.e-9); vy = vy/(vel+1.e-9); - for(i=0;ivalues[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0][i]*(-vy)+dbasis[1][i]*vx)*Jdet*gauss->weight; + for(i=0;ivalues[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0*numnodes+i]*(-vy)+dbasis[1*numnodes+i]*vx)*Jdet*gauss->weight; break; default: @@ -4996,6 +4994,8 @@ /*Clean up and return*/ + xDelete(basis); + xDelete(dbasis); + xDelete(responses); delete gauss; - xDelete(responses); return pe; } @@ -5874,5 +5874,5 @@ /*Fetch number of nodes and dof for this finite element*/ int numnodes = this->NumberofNodes(); - int numdof = numnodes*NDOF2; + int numdof = numnodes*NDOF1; /*Initialize Element matrix and vectors*/ @@ -5928,7 +5928,4 @@ ElementMatrix* Tria::CreateKMatrixHydrologyDCEfficient(void){ - /*constants: */ - const int numdof=NDOF1*NUMVERTICES; - /* Intermediaries */ IssmDouble D_scalar,Jdet; @@ -5936,8 +5933,4 @@ IssmDouble epl_storing; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble B[2][numdof]; - IssmDouble basis[NUMVERTICES]; - IssmDouble D[2][2]; - GaussTria *gauss = NULL; /*Check that all nodes are active, else return empty matrix*/ @@ -5946,6 +5939,13 @@ } - /*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*NDOF1; + + /*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*/ @@ -5956,5 +5956,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++){ @@ -5967,22 +5967,25 @@ 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); /*Transient*/ if(reCast(dt)){ - GetNodalFunctions(&basis[0],gauss); + GetNodalFunctions(basis,gauss); D_scalar=epl_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; @@ -6734,6 +6737,6 @@ /*Clean up and return*/ + xDelete(basis); delete gauss; - xDelete(basis); return pe; } @@ -6944,25 +6947,22 @@ ElementMatrix* Tria::CreateKMatrixBalancethickness_CG(void){ - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - /*Intermediaries */ int stabilization; int i,j,dim; IssmDouble Jdettria,vx,vy,dvxdx,dvydy,vel,h; + IssmDouble D_scalar; IssmDouble dvx[2],dvy[2]; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble L[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}; - IssmDouble DL_scalar; - 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*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]; /*Retrieve all Inputs and parameters: */ @@ -6983,5 +6983,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++){ @@ -6989,6 +6989,6 @@ GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); - GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss); - GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); + GetBPrognostic(B,&xyz_list[0][0],gauss); + GetBprimePrognostic(Bprime,&xyz_list[0][0],gauss); vxaverage_input->GetInputValue(&vx,gauss); @@ -6999,20 +6999,20 @@ dvxdx=dvx[0]; dvydy=dvy[1]; - DL_scalar=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, + D_scalar=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); - TripleMultiply( &B[0][0],2,numdof,1, - &DLprime[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, + 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); @@ -7020,8 +7020,8 @@ /*Streamline upwinding*/ vel=sqrt(vx*vx+vy*vy); - K[0][0]=h/(2*vel)*vx*vx; - K[1][0]=h/(2*vel)*vy*vx; - K[0][1]=h/(2*vel)*vx*vy; - K[1][1]=h/(2*vel)*vy*vy; + 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==2){ @@ -7029,17 +7029,17 @@ vxaverage_input->GetInputAverage(&vx); vyaverage_input->GetInputAverage(&vy); - K[0][0]=h/2.0*fabs(vx); - K[0][1]=0.; - K[1][0]=0.; - K[1][1]=h/2.0*fabs(vy); + 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){ - 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, + D[0][0]=D_scalar*D[0][0]; + D[1][0]=D_scalar*D[1][0]; + D[0][1]=D_scalar*D[0][1]; + D[1][1]=D_scalar*D[1][1]; + TripleMultiply(Bprime,2,numdof,1, + &D[0][0],2,2,0, + Bprime,2,numdof,0, &Ke->values[0],1); } @@ -7047,4 +7047,7 @@ /*Clean up and return*/ + xDelete(basis); + xDelete(B); + xDelete(Bprime); delete gauss; return Ke; @@ -7122,16 +7125,15 @@ ElementVector* Tria::CreatePVectorBalancethickness_CG(void){ - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - /*Intermediaries */ - int i,j; IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria; - IssmDouble basis[NUMVERTICES]; - GaussTria* gauss=NULL; - - /*Initialize Element vector*/ - ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF1; + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); + IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ @@ -7142,5 +7144,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++){ @@ -7152,10 +7154,11 @@ GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); - GetNodalFunctions(&basis[0],gauss); - - for(i=0;ivalues[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i]; + GetNodalFunctions(basis,gauss); + + for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i]; } /*Clean up and return*/ + xDelete(basis); delete gauss; return pe;