Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15485) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15486) @@ -247,9 +247,9 @@ /*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*/ ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); - IssmDouble* basis = xNew(numdof); + IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ @@ -266,7 +266,7 @@ D=gauss->weight*Jdet; - TripleMultiply(basis,1,numnodes,1, + TripleMultiply(basis,1,numdof,1, &D,1,1,0, - basis,1,numnodes,0, + basis,1,numdof,0, &Ke->values[0],1); } @@ -6244,5 +6244,5 @@ switch(GetElementType()){ - case P1Enum: + case P1Enum: case P2Enum: return CreateKMatrixPrognostic_CG(); case P1DGEnum: @@ -6256,25 +6256,23 @@ ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){ - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - /*Intermediaries */ int stabilization; int dim; - IssmDouble Jdettria,DL_scalar,dt,h; + IssmDouble Jdettria,D_scalar,dt,h; IssmDouble vel,vx,vy,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}; - - /*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*/ @@ -6302,5 +6300,5 @@ GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); - GetNodalFunctions(&basis[0],gauss); + GetNodalFunctions(basis,gauss); vxaverage_input->GetInputValue(&vx,gauss); @@ -6309,31 +6307,32 @@ vyaverage_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, + D_scalar=gauss->weight*Jdettria; + + TripleMultiply(basis,1,numdof,1, + &D_scalar,1,1,0, + basis,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); + 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, + D_scalar=dt*gauss->weight*Jdettria; + + D[0][0]=D_scalar*dvxdx; + D[0][1]=0.; + D[1][1]=D_scalar*dvydy; + D[1][0]=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); @@ -6341,8 +6340,8 @@ /*Streamline upwinding*/ vel=sqrt(vx*vx+vy*vy)+1.e-8; - 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==1){ @@ -6350,17 +6349,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); } @@ -6369,4 +6368,7 @@ /*Clean up and return*/ delete gauss; + xDelete(basis); + xDelete(B); + xDelete(Bprime); return Ke; } @@ -6451,5 +6453,5 @@ switch(GetElementType()){ - case P1Enum: + case P1Enum: case P2Enum: return CreatePVectorPrognostic_CG(); case P1DGEnum: @@ -6463,28 +6465,28 @@ ElementVector* Tria::CreatePVectorPrognostic_CG(void){ - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - /*Intermediaries */ IssmDouble Jdettria,dt; IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g; IssmDouble xyz_list[NUMVERTICES][3]; - 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*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input); - Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); - Input* basal_melting_correction_input=inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* surface_mass_balance_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input); + Input* basal_melting_input = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); + Input* basal_melting_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum); + Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); /*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++){ @@ -6492,5 +6494,5 @@ GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); - GetNodalFunctions(&basis[0],gauss); + GetNodalFunctions(basis,gauss); surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss); @@ -6507,4 +6509,5 @@ /*Clean up and return*/ delete gauss; + xDelete(basis); return pe; }