Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h =================================================================== --- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 25262) +++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 25263) @@ -29,8 +29,6 @@ ElementMatrix* CreateJacobianMatrix(Element* element); ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); - 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); void GetSolutionFromInputs(Vector* solution,Element* element); void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); void InputUpdateFromSolution(IssmDouble* solution,Element* element); Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 25262) +++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 25263) @@ -453,14 +453,12 @@ /*Initialize Element vector*/ ElementMatrix* Ke = element->NewElementMatrix(); IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(dim*numnodes); - IssmDouble* Bprime = xNew(dim*numnodes); + IssmDouble* dbasis = xNew(dim*numnodes); IssmDouble* D = xNewZeroInit(dim*dim); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); element->FindParam(&dt,TimesteppingTimeStepEnum); - //printf("dt %f\n", dt); element->FindParam(&stabilization,DamageStabilizationEnum); Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input); Input2* vy_input = element->GetInput2(VyEnum); _assert_(vy_input); @@ -478,6 +476,7 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); vx_input->GetInputValue(&vx,gauss); vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); @@ -489,38 +488,36 @@ vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); } + /*Transient term*/ D_scalar=gauss->weight*Jdet; - TripleMultiply(basis,1,numnodes,1, - &D_scalar,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); + for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[i]*basis[j]; - GetB(B,element,dim,xyz_list,gauss); - GetBprime(Bprime,element,dim,xyz_list,gauss); - dvxdx=dvx[0]; dvydy=dvy[1]; - if(dim==3) dvzdz=dvz[2]; D_scalar=dt*gauss->weight*Jdet; + if(dim==2){ + for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy); + /*\phi_i v\cdot\nabla\phi_j*/ + Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]); + } + } + } + else{/*3D*/ + _assert_(dim==3); + dvzdz=dvz[2]; + for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy+dvzdz); + /*\phi_i v\cdot\nabla\phi_j*/ + Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j] + vz*dbasis[2*numnodes+j]); + } + } + } - D[0*dim+0]=D_scalar*dvxdx; - D[1*dim+1]=D_scalar*dvydy; - if(dim==3) D[2*dim+2]=D_scalar*dvzdz; - - TripleMultiply(B,dim,numnodes,1, - D,dim,dim,0, - B,dim,numnodes,0, - &Ke->values[0],1); - - D[0*dim+0]=D_scalar*vx; - D[1*dim+1]=D_scalar*vy; - if(dim==3) D[2*dim+2]=D_scalar*vz; - - TripleMultiply(B,dim,numnodes,1, - D,dim,dim,0, - Bprime,dim,numnodes,0, - &Ke->values[0],1); - if(stabilization==2){ if(dim==3){ vel=sqrt(vx*vx+vy*vy+vz*vz)+1.e-8; @@ -565,6 +562,14 @@ D[1*dim+0]=D_scalar*D[1*dim+0]; D[0*dim+1]=D_scalar*D[0*dim+1]; D[1*dim+1]=D_scalar*D[1*dim+1]; + for(int i=0;ivalues[i*numnodes+j] += ( + dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) + + dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) + ); + } + } } else if(dim==3){ D[0*dim+0]=D_scalar*D[0*dim+0]; @@ -576,20 +581,23 @@ D[0*dim+2]=D_scalar*D[0*dim+2]; D[1*dim+2]=D_scalar*D[1*dim+2]; D[2*dim+2]=D_scalar*D[2*dim+2]; + for(int i=0;ivalues[i*numnodes+j] += ( + dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j] + D[0*dim+2]*dbasis[2*numnodes+j]) + + dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j] + D[1*dim+2]*dbasis[2*numnodes+j]) + + dbasis[2*numnodes+i] *(D[2*dim+0]*dbasis[0*numnodes+j] + D[2*dim+1]*dbasis[1*numnodes+j] + D[2*dim+2]*dbasis[2*numnodes+j]) + ); + } + } } - TripleMultiply(Bprime,dim,numnodes,1, - D,dim,dim,0, - Bprime,dim,numnodes,0, - &Ke->values[0],1); } - } /*Clean up and return*/ xDelete(xyz_list); xDelete(basis); - xDelete(B); - xDelete(Bprime); + xDelete(dbasis); xDelete(D); delete gauss; return Ke; @@ -665,63 +673,6 @@ delete gauss; return pe; }/*}}}*/ -void DamageEvolutionAnalysis::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 DamageEvolutionAnalysis::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); - -}/*}}}*/ void DamageEvolutionAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ element->GetSolutionFromInputsOneDof(solution,DamageDbarEnum); }/*}}}*/ @@ -787,9 +738,8 @@ /*Initialize Element vector and other vectors*/ ElementMatrix* Ke = element->NewElementMatrix(); - IssmDouble* B = xNew(dim*numnodes); - IssmDouble* Bprime = xNew(dim*numnodes); - IssmDouble* D = xNewZeroInit(dim*dim); + IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(dim*numnodes); /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); @@ -802,26 +752,24 @@ gauss->GaussPoint(ig); 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); - + for(int i=0;ivalues[i*numnodes+j] += -gauss->weight*Jdet*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]); + } + } } /*Clean up and return*/ xDelete(xyz_list); - xDelete(B); - xDelete(Bprime); - xDelete(D); + xDelete(basis); + xDelete(dbasis); delete gauss; return Ke; }/*}}}*/ @@ -853,10 +801,7 @@ element->NodalFunctions(basis,gauss); D=gauss->weight*Jdet; - TripleMultiply(basis,1,numnodes,1, - &D,1,1,0, - basis,1,numnodes,0, - &Me->values[0],1); + for(int i=0;ivalues[i*numnodes+j] += D*basis[i]*basis[j]; } /*Clean up and return*/ Index: ../trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 25262) +++ ../trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 25263) @@ -107,11 +107,7 @@ basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); basalelement->NodalFunctions(basis,gauss); D=gauss->weight*Jdet; - - TripleMultiply(basis,1,numnodes,1, - &D,1,1,0, - basis,1,numnodes,0, - &Ke->values[0],1); + for(int i=0;ivalues[i*numnodes+j] += D*basis[i]*basis[j]; } /*Clean up and return*/ Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 25262) +++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 25263) @@ -336,11 +336,7 @@ /*Transient term*/ D_scalar=gauss->weight*Jdet; - for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[i]*basis[j]; - } - } + for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[i]*basis[j]; /*Advection terms*/ vxaverage_input->GetInputValue(&vx,gauss);