Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 25261) +++ /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 25262) @@ -264,6 +264,6 @@ /*Initialize Element vector and other vectors*/ ElementMatrix* Ke = element->NewElementMatrix(); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); + IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(2*numnodes); IssmDouble D[2][2]; @@ -288,4 +288,6 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); vxaverage_input->GetInputValue(&vx,gauss); vyaverage_input->GetInputValue(&vy,gauss); @@ -294,8 +296,5 @@ D_scalar=gauss->weight*Jdet; - /*WARNING: B and Bprime are inverted compared to CG*/ - GetB(Bprime,element,xyz_list,gauss); - GetBprime(B,element,xyz_list,gauss); - + /*WARNING: inverted compared to CG*/ D_scalar = - gauss->weight*Jdet; D[0][0] = D_scalar*vx; @@ -303,15 +302,16 @@ D[1][0] = 0.; D[1][1] = D_scalar*vy; - TripleMultiply(B,2,numnodes,1, - &D[0][0],2,2,0, - Bprime,2,numnodes,0, - &Ke->values[0],1); - + for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy); + Ke->values[i*numnodes+j] += D_scalar*basis[j]*(vx*dbasis[0*numnodes+i] + vy*dbasis[1*numnodes+i]); + } + } } /*Clean up and return*/ xDelete(xyz_list); - xDelete(B); - xDelete(Bprime); + xDelete(basis); + xDelete(dbasis); delete gauss; return Ke; @@ -420,59 +420,4 @@ delete gauss; return pe; -}/*}}}*/ -void BalancethicknessAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*2. - * For node i, Bi can be expressed in the actual coordinate system - * by: - * Bi=[ N ] - * [ N ] - * where N is the finiteelement function for node i. - * - * We assume B_prog has been allocated already, of size: 2x(1*numnodes) - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Get nodal functions*/ - IssmDouble* basis=xNew(numnodes); - element->NodalFunctions(basis,gauss); - - /*Build B: */ - for(int i=0;i(basis); -}/*}}}*/ -void BalancethicknessAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. - * For node i, Bi' can be expressed in the actual coordinate system - * by: - * Bi_prime=[ dN/dx ] - * [ dN/dy ] - * where N is the finiteelement function for node i. - * - * We assume B' has been allocated already, of size: 3x(2*numnodes) - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Get nodal functions derivatives*/ - IssmDouble* dbasis=xNew(2*numnodes); - element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - - /*Build B': */ - for(int i=0;i(dbasis); - }/*}}}*/ void BalancethicknessAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.h (revision 25261) +++ /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.h (revision 25262) @@ -30,6 +30,4 @@ ElementVector* CreatePVectorCG(Element* element); ElementVector* CreatePVectorDG(Element* element); - void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); void GetSolutionFromInputs(Vector* solution,Element* element); void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); Index: /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (revision 25261) +++ /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (revision 25262) @@ -109,5 +109,4 @@ basalelement->NodalFunctions(basis,gauss); basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - GetBprime(Bprime,basalelement,dim,xyz_list,gauss); /*Get velocity*/ @@ -161,8 +160,12 @@ D[0][0]=D_scalar*h/(2.*vel)*fabs(vx*vx); D[0][1]=D_scalar*h/(2.*vel)*fabs(vx*vy); D[1][0]=D_scalar*h/(2.*vel)*fabs(vy*vx); D[1][1]=D_scalar*h/(2.*vel)*fabs(vy*vy); - TripleMultiply(Bprime,dim,numnodes,1, - &D[0][0],2,2,0, - Bprime,dim,numnodes,0, - &Ke->values[0],1); + for(int i=0;ivalues[i*numnodes+j] += ( + dbasis[0*numnodes+i] *(D[0][0]*dbasis[0*numnodes+j] + D[0][1]*dbasis[1*numnodes+j]) + + dbasis[1*numnodes+i] *(D[1][0]*dbasis[0*numnodes+j] + D[1][1]*dbasis[1*numnodes+j]) + ); + } + } } } @@ -172,5 +175,4 @@ xDelete(basis); xDelete(dbasis); - xDelete(Bprime); delete gauss; if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; @@ -235,60 +237,2 @@ return; }/*}}}*/ - -void GLheightadvectionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*2. - * For node i, Bi can be expressed in the actual coordinate system - * by: - * Bi=[ N ] - * [ N ] - * where N is the finiteelement function for node i. - * - * We assume B_prog has been allocated already, of size: 2x(1*numnodes) - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Get nodal functions*/ - IssmDouble* basis=xNew(numnodes); - element->NodalFunctions(basis,gauss); - - /*Build B: */ - for(int i=0;i(basis); -}/*}}}*/ -void GLheightadvectionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. - * For node i, Bi' can be expressed in the actual coordinate system - * by: - * Bi_prime=[ dN/dx ] - * [ dN/dy ] - * where N is the finiteelement function for node i. - * - * We assume B' has been allocated already, of size: 3x(2*numnodes) - */ - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Get nodal functions derivatives*/ - IssmDouble* dbasis=xNew(dim*numnodes); - element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - - /*Build B': */ - for(int i=0;i(dbasis); - -}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.h (revision 25261) +++ /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.h (revision 25262) @@ -30,7 +30,4 @@ void InputUpdateFromSolution(IssmDouble* solution,Element* element); void UpdateConstraints(FemModel* femmodel); - - void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); - void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); }; #endif