Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 16874) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 16875) @@ -96,7 +96,4 @@ /*Finite Element Analysis*/ ElementMatrix* DamageEvolutionAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); -}/*}}}*/ -ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/ /*Intermediaries*/ @@ -118,4 +115,134 @@ /*Intermediaries */ + int stabilization; + IssmDouble Jdet,dt,D_scalar; + IssmDouble vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2]; + IssmDouble *xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = basalelement->GetNumberOfNodes(); + + /*Initialize Element vector*/ + ElementMatrix* Ke = basalelement->NewElementMatrix(NoneApproximationEnum); + IssmDouble* basis = xNew(numnodes); + IssmDouble* B = xNew(2*numnodes); + IssmDouble* Bprime = xNew(2*numnodes); + IssmDouble D[2][2]={0.}; + + /*Retrieve all inputs and parameters*/ + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->FindParam(&dt,TimesteppingTimeStepEnum); + basalelement->FindParam(&stabilization,DamageStabilizationEnum); + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; + if(meshtype==Mesh2DhorizontalEnum){ + vxaverage_input=basalelement->GetInput(VxEnum); _assert_(vxaverage_input); + vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input); + } + else{ + vxaverage_input=basalelement->GetInput(VxAverageEnum); _assert_(vxaverage_input); + vyaverage_input=basalelement->GetInput(VyAverageEnum); _assert_(vyaverage_input); + } + IssmDouble h=basalelement->CharacteristicLength(); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctions(basis,gauss); + D_scalar=gauss->weight*Jdet; + + vxaverage_input->GetInputValue(&vx,gauss); + vyaverage_input->GetInputValue(&vy,gauss); + vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + + TripleMultiply(basis,1,numnodes,1, + &D_scalar,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + + GetB(B,element,xyz_list,gauss); + GetBprime(Bprime,element,xyz_list,gauss); + + dvxdx=dvx[0]; + dvydy=dvy[1]; + D_scalar=dt*gauss->weight*Jdet; + + D[0][0]=D_scalar*dvxdx; + D[1][1]=D_scalar*dvydy; + TripleMultiply(B,2,numnodes,1, + &D[0][0],2,2,0, + B,2,numnodes,0, + &Ke->values[0],1); + + D[0][0]=D_scalar*vx; + 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); + + if(stabilization==2){ + /*Streamline upwinding*/ + vel=sqrt(vx*vx+vy*vy)+1.e-8; + 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){ + /*SSA*/ + vxaverage_input->GetInputAverage(&vx); + vyaverage_input->GetInputAverage(&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){ + 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,numnodes,1, + &D[0][0],2,2,0, + Bprime,2,numnodes,0, + &Ke->values[0],1); + } + + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(B); + xDelete(Bprime); + delete gauss; + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + return Ke; +}/*}}}*/ +ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/ + + /*Intermediaries*/ + int meshtype; + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + basalelement = element; + break; + case Mesh3DEnum: + if(!element->IsOnBed()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<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,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); + }/*}}}*/ void DamageEvolutionAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 16874) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 16875) @@ -23,4 +23,6 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(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 InputUpdateFromSolution(IssmDouble* solution,Element* element);