Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 24428) +++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 24429) @@ -116,10 +116,9 @@ /*Initialize Element vector and other vectors*/ - ElementMatrix* Ke = workelement->NewElementMatrix(); - IssmDouble* B = xNew(dim*numnodes); - IssmDouble* Bprime = xNew(dim*numnodes); - IssmDouble* dlsf = xNew(dim); - IssmDouble* normal = xNew(dim); - IssmDouble* D = xNewZeroInit(dim*dim); + ElementMatrix *Ke = workelement->NewElementMatrix(); + IssmDouble *basis = xNew(numnodes); + IssmDouble *dbasis = xNew(dim*numnodes); + IssmDouble dlsf[3]; + IssmDouble normal[3]; /*Retrieve all inputs and parameters*/ @@ -143,18 +142,15 @@ workelement->JacobianDeterminant(&Jdet,xyz_list,gauss); - GetB(B,workelement,xyz_list,gauss,dim); - GetBprime(Bprime,workelement,xyz_list,gauss,dim); + workelement->NodalFunctions(basis,gauss); + workelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); D_scalar=gauss->weight*Jdet; if(extrapolatebydiffusion){ - - /* diffuse values outward only along the xy-plane*/ - for(int i=0;i<2;i++) D[i*dim+i] = D_scalar; - - TripleMultiply(Bprime,dim,numnodes,1, - D,dim,dim,0, - Bprime,dim,numnodes,0, - &Ke->values[0],1); + for(int i=0;ivalues[i*numnodes+j] += D_scalar*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]); + } + } } else{ @@ -175,14 +171,9 @@ for(i=0;ivalues[0],1); + for(int i=0;ivalues[i*numnodes+j] += D_scalar*(normal[0]*dbasis[0*numnodes+j]*basis[i] + normal[1]*dbasis[1*numnodes+j]*basis[i]); + } + } /* stabilization */ @@ -195,14 +186,10 @@ h=sqrt(pow(hx*normal[0],2) + pow(hy*normal[1],2)); kappa=h/2.+1.e-14; - for(row=0;rowvalues[0],1); + + for(int i=0;ivalues[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]); + } + } } } @@ -211,9 +198,6 @@ /*Clean up and return*/ xDelete(xyz_list); - xDelete(B); - xDelete(Bprime); - xDelete(D); - xDelete(dlsf); - xDelete(normal); + xDelete(basis); + xDelete(dbasis); delete gauss; if(extrapolationcase==0){workelement->DeleteMaterials(); delete workelement;}; Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 24428) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 24429) @@ -521,14 +521,9 @@ basalelement->ElementSizes(&hx,&hy,&hz); h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) ); - IssmDouble D[9]; - for(row=0;rowvalues[0],1); + for(int i=0;ivalues[i*numnodes+j] += D_scalar*h/(2.*vel)*w[i]*(w[j]*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + w[j]*dbasis[1*numnodes+j]*dbasis[1*numnodes+i]); + } + } } break; @@ -595,59 +590,4 @@ return pe; -}/*}}}*/ -void LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. - * 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(NDOF1*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 LevelsetAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. - * 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(NDOF2*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); - }/*}}}*/ IssmDouble LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 24428) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 24429) @@ -26,6 +26,4 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); - void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); - void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss); IssmDouble GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1); void GetSolutionFromInputs(Vector* solution,Element* element); Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 24428) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 24429) @@ -425,4 +425,5 @@ if(stabilization==2){ /*Streamline upwind*/ + _assert_(dim==2); for(int i=0;iNewElementMatrix(); IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]; + IssmDouble* dbasis = xNew(3*numnodes); /*Retrieve all inputs and parameters*/ @@ -517,4 +516,5 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); vxaverage_input->GetInputValue(&vx,gauss); @@ -522,23 +522,17 @@ D_scalar=gauss->weight*Jdet; - TripleMultiply(basis,1,numnodes,1, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); - - /*WARNING: B and Bprime are inverted compared to CG*/ - GetB(Bprime,element,2,xyz_list,gauss); - GetBprime(B,element,2,xyz_list,gauss); - + for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[i]*basis[j]; + } + } + + /*WARNING: basis and dbasis are inverted compared to CG*/ D_scalar = - dt*gauss->weight*Jdet; - D[0][0] = D_scalar*vx; - D[0][1] = 0.; - 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*(vx*dbasis[0*numnodes+i]*basis[j] + vy*dbasis[1*numnodes+i]*basis[j]); + } + } } @@ -546,6 +540,5 @@ xDelete(xyz_list); xDelete(basis); - xDelete(B); - xDelete(Bprime); + xDelete(dbasis); delete gauss; return Ke; @@ -989,5 +982,5 @@ /*Intermediaries */ - IssmDouble Jdet; + IssmDouble Jdet,D_scalar; IssmDouble vx,vy; IssmDouble* xyz_list = NULL; @@ -999,6 +992,6 @@ /*Initialize Element vector and other vectors*/ ElementMatrix* Ke = element->NewElementMatrix(); - IssmDouble* B = xNew(dim*numnodes); - IssmDouble* Bprime = xNew(dim*numnodes); + IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(dim*numnodes); IssmDouble* D = xNewZeroInit(dim*dim); @@ -1014,24 +1007,21 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); - GetB(B,element,dim,xyz_list,gauss); - GetBprime(Bprime,element,dim,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); vxaverage_input->GetInputValue(&vx,gauss); vyaverage_input->GetInputValue(&vy,gauss); - D[0*dim+0] = -gauss->weight*vx*Jdet; - D[1*dim+1] = -gauss->weight*vy*Jdet; - - TripleMultiply(B,dim,numnodes,1, - D,dim,dim,0, - Bprime,dim,numnodes,0, - &Ke->values[0],1); - + D_scalar = gauss->weight*Jdet; + for(int i=0;ivalues[i*numnodes+j] += -D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]); + } + } } /*Clean up and return*/ xDelete(xyz_list); - xDelete(B); - xDelete(Bprime); - xDelete(D); + xDelete(basis); + xDelete(dbasis); delete gauss; return Ke; Index: /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp (revision 24428) +++ /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp (revision 24429) @@ -124,9 +124,9 @@ /*Initialize Element matrix and vectors*/ - ElementVector* pe = element->NewElementVector(); - IssmDouble* basis = xNew(numnodes); - IssmDouble* dvx = xNew(dim); - IssmDouble* dvy = xNew(dim); - IssmDouble* dvz = xNew(dim); + ElementVector *pe = element->NewElementVector(); + IssmDouble *basis = xNew(numnodes); + IssmDouble dvx[3]; + IssmDouble dvy[3]; + IssmDouble dvz[3]; Input2* vx_input=element->GetInput2(VxEnum); _assert_(vx_input); @@ -159,7 +159,4 @@ xDelete(xyz_list); xDelete(basis); - xDelete(dvx); - xDelete(dvy); - xDelete(dvz); return pe; }/*}}}*/