Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 18469) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 18470) @@ -75,8 +75,5 @@ iomodel->Constant(&finiteelement,DamageElementinterpEnum); - - if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(1,MeshVertexonbaseEnum); ::CreateNodes(nodes,iomodel,DamageEvolutionAnalysisEnum,finiteelement); - iomodel->DeleteData(1,MeshVertexonbaseEnum); }/*}}}*/ void DamageEvolutionAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ @@ -109,33 +106,24 @@ /* Check if ice in element */ if(!element->IsIceInElement()) return NULL; - /*Intermediaries*/ - Element* basalelement; int domaintype,dim; int stabilization; - IssmDouble Jdet,dt,D_scalar,h; - IssmDouble vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2]; + IssmDouble Jdet,dt,D_scalar,h,hx,hy,hz; + IssmDouble vel,vx,vy,vz,dvxdx,dvydy,dvzdz,dvx[3],dvy[3],dvz[3]; IssmDouble *xyz_list = NULL; - /*Get problem dimension and basal element*/ + /*Get problem dimension*/ element->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ - case Domain2DhorizontalEnum: - basalelement = element; - dim = 2; - break; - case Domain3DEnum: - if(!element->IsOnBase()) return NULL; - basalelement = element->SpawnBasalElement(); - dim = 2; - break; - default: _error_("mesh "<GetNumberOfNodes(); + int numnodes = element->GetNumberOfNodes(); /*Initialize Element vector*/ - ElementMatrix* Ke = basalelement->NewElementMatrix(); + ElementMatrix* Ke = element->NewElementMatrix(); IssmDouble* basis = xNew(numnodes); IssmDouble* B = xNew(dim*numnodes); @@ -144,37 +132,32 @@ /*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(domaintype==Domain2DhorizontalEnum){ - vxaverage_input=basalelement->GetInput(VxEnum); _assert_(vxaverage_input); - vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input); - } - else{ - if(dim==1){ - vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input); - } - if(dim==2){ - vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input); - vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input); - } - } - h=basalelement->CharacteristicLength(); + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + element->FindParam(&stabilization,DamageStabilizationEnum); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = NULL; + if(dim==3){ + vz_input=element->GetInput(VzEnum); _assert_(vz_input); + } + + if(dim==2) h=element->CharacteristicLength(); /* Start looping on the number of gaussian points: */ - Gauss* gauss=basalelement->NewGauss(2); + Gauss* gauss=element->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); - basalelement->NodalFunctions(basis,gauss); + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); - vxaverage_input->GetInputValue(&vx,gauss); - vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - if(dim==2){ - vyaverage_input->GetInputValue(&vy,gauss); - vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputValue(&vy,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + + if(dim==3){ + vz_input->GetInputValue(&vz,gauss); + vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); } @@ -185,13 +168,15 @@ &Ke->values[0],1); - GetB(B,basalelement,dim,xyz_list,gauss); - GetBprime(Bprime,basalelement,dim,xyz_list,gauss); + GetB(B,element,dim,xyz_list,gauss); + GetBprime(Bprime,element,dim,xyz_list,gauss); dvxdx=dvx[0]; - if(dim==2) dvydy=dvy[1]; + dvydy=dvy[1]; + if(dim==3) dvzdz=dvz[2]; D_scalar=dt*gauss->weight*Jdet; D[0*dim+0]=D_scalar*dvxdx; - if(dim==2) D[1*dim+1]=D_scalar*dvydy; + D[1*dim+1]=D_scalar*dvydy; + if(dim==3) D[2*dim+2]=D_scalar*dvzdz; TripleMultiply(B,dim,numnodes,1, @@ -201,5 +186,6 @@ D[0*dim+0]=D_scalar*vx; - if(dim==2) D[1*dim+1]=D_scalar*vy; + D[1*dim+1]=D_scalar*vy; + if(dim==3) D[2*dim+2]=D_scalar*vz; TripleMultiply(B,dim,numnodes,1, @@ -209,7 +195,15 @@ if(stabilization==2){ - if(dim==1){ - vel=fabs(vx)+1.e-8; - D[0]=h/(2.0*vel)*vx*vx; + if(dim==3){ + vel=sqrt(vx*vx+vy*vy+vz*vz)+1.e-8; + D[0*dim+0]=h/(2.0*vel)*vx*vx; + D[1*dim+0]=h/(2.0*vel)*vy*vx; + D[2*dim+0]=h/(2.0*vel)*vz*vx; + D[0*dim+1]=h/(2.0*vel)*vx*vy; + D[1*dim+1]=h/(2.0*vel)*vy*vy; + D[2*dim+1]=h/(2.0*vel)*vy*vz; + D[0*dim+2]=h/(2.0*vel)*vx*vz; + D[1*dim+2]=h/(2.0*vel)*vy*vz; + D[2*dim+2]=h/(2.0*vel)*vz*vz; } else{ @@ -223,12 +217,21 @@ } else if(stabilization==1){ - vxaverage_input->GetInputAverage(&vx); - if(dim==2) vyaverage_input->GetInputAverage(&vy); - D[0*dim+0]=h/2.0*fabs(vx); - if(dim==2) D[1*dim+1]=h/2.0*fabs(vy); + if(dim==2){ + vx_input->GetInputAverage(&vx); + vy_input->GetInputAverage(&vy); + D[0*dim+0]=h/2.0*fabs(vx); + D[1*dim+1]=h/2.0*fabs(vy); + } + else if(dim==3){ + element->ElementSizes(&hx,&hy,&hz); + vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14; + h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2)); + D[0*dim+0]=h/(2.*vel)*fabs(vx*vx); D[0*dim+1]=h/(2.*vel)*fabs(vx*vy); D[0*dim+2]=h/(2.*vel)*fabs(vx*vz); + D[1*dim+0]=h/(2.*vel)*fabs(vy*vx); D[1*dim+1]=h/(2.*vel)*fabs(vy*vy); D[1*dim+2]=h/(2.*vel)*fabs(vy*vz); + D[2*dim+0]=h/(2.*vel)*fabs(vz*vx); D[2*dim+1]=h/(2.*vel)*fabs(vz*vy); D[2*dim+2]=h/(2.*vel)*fabs(vz*vz); + } } if(stabilization==1 || stabilization==2){ - if(dim==1) D[0]=D_scalar*D[0]; - else{ + if(dim==2){ D[0*dim+0]=D_scalar*D[0*dim+0]; D[1*dim+0]=D_scalar*D[1*dim+0]; @@ -236,5 +239,15 @@ D[1*dim+1]=D_scalar*D[1*dim+1]; } - + else if(dim==3){ + D[0*dim+0]=D_scalar*D[0*dim+0]; + D[1*dim+0]=D_scalar*D[1*dim+0]; + D[2*dim+0]=D_scalar*D[2*dim+0]; + D[0*dim+1]=D_scalar*D[0*dim+1]; + D[1*dim+1]=D_scalar*D[1*dim+1]; + D[2*dim+1]=D_scalar*D[2*dim+1]; + 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]; + } TripleMultiply(Bprime,dim,numnodes,1, D,dim,dim,0, @@ -252,5 +265,4 @@ xDelete(D); delete gauss; - if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; return Ke; }/*}}}*/ @@ -262,51 +274,52 @@ /*Intermediaries*/ int domaintype,damagelaw; - Element* basalelement; IssmDouble Jdet,dt; IssmDouble f,damage; IssmDouble* xyz_list = NULL; - - /*Get basal element*/ + /*Get element*/ element->FindParam(&domaintype,DomainTypeEnum); - switch(domaintype){ - case Domain2DhorizontalEnum: - basalelement = element; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = element->NewElementVector(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + element->FindParam(&damagelaw,DamageLawEnum); + switch(damagelaw){ + case 1: + this->CreateDamageFInputPralong(element); break; - case Domain3DEnum: - if(!element->IsOnBase()) return NULL; - basalelement = element->SpawnBasalElement(); + case 2: + this->CreateDamageFInputPralong(element); break; - default: _error_("mesh "<GetNumberOfNodes(); - - /*Initialize Element vector and other vectors*/ - ElementVector* pe = basalelement->NewElementVector(); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - basalelement->GetVerticesCoordinates(&xyz_list); - basalelement->FindParam(&dt,TimesteppingTimeStepEnum); - basalelement->FindParam(&damagelaw,DamageLawEnum); - if(damagelaw==1 | damagelaw==2){ - this->CreateDamageFInputPralong(basalelement); - } - else if(damagelaw==3){ - this->CreateDamageFInputExp(basalelement); - } - Input* damaged_input = basalelement->GetInput(DamageDEnum); _assert_(damaged_input); - Input* damagef_input = basalelement->GetInput(DamageFEnum); _assert_(damagef_input); + case 3: + this->CreateDamageFInputExp(element); + break; + default: + _error_("not implemented yet"); + } + + Input* damaged_input = NULL; + Input* damagef_input = element->GetInput(DamageFEnum); _assert_(damagef_input); + if(domaintype==Domain2DhorizontalEnum){ + damaged_input = element->GetInput(DamageDbarEnum); _assert_(damaged_input); + } + else{ + damaged_input = element->GetInput(DamageDEnum); _assert_(damaged_input); + } /* Start looping on the number of gaussian points: */ - Gauss* gauss=basalelement->NewGauss(2); + Gauss* gauss=element->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); - basalelement->NodalFunctions(basis,gauss); + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); damaged_input->GetInputValue(&damage,gauss); @@ -317,9 +330,7 @@ } } - /*Clean up and return*/ xDelete(xyz_list); xDelete(basis); - if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; delete gauss; return pe; @@ -393,23 +404,16 @@ IssmDouble max_damage; int *doflist = NULL; - Element* basalelement=NULL; element->FindParam(&domaintype,DomainTypeEnum); - if(domaintype!=Domain2DhorizontalEnum){ - if(!element->IsOnBase()) return; - basalelement=element->SpawnBasalElement(); - } - else{ - basalelement = element; - } + /*Fetch number of nodes and dof for this finite element*/ - int numnodes = basalelement->GetNumberOfNodes(); + int numnodes = element->GetNumberOfNodes(); /*Fetch dof list and allocate solution vector*/ - basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); IssmDouble* newdamage = xNew(numnodes); /*Get user-supplied max_damage: */ - basalelement->FindParam(&max_damage,DamageMaxDamageEnum); + element->FindParam(&max_damage,DamageMaxDamageEnum); /*Use the dof list to index into the solution vector: */ @@ -425,8 +429,8 @@ /*Get all inputs and parameters*/ if(domaintype==Domain2DhorizontalEnum){ - basalelement->AddInput(DamageDbarEnum,newdamage,element->GetElementType()); + element->AddInput(DamageDbarEnum,newdamage,element->GetElementType()); } else{ - basalelement->AddBasalInput(DamageDEnum,newdamage,element->GetElementType()); + element->AddInput(DamageDEnum,newdamage,element->GetElementType()); } @@ -434,5 +438,4 @@ xDelete(newdamage); xDelete(doflist); - if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ @@ -446,8 +449,8 @@ /*Intermediaries */ IssmDouble c1,c2,c3,healing,stress_threshold; - IssmDouble s_xx,s_xy,s_yy,s1,s2,stmp; + IssmDouble s_xx,s_xy,s_xz,s_yy,s_yz,s_zz,s1,s2,s3,stmp; IssmDouble J2s,Chi,Psi,PosPsi,NegPsi; - IssmDouble damage,tau_xx,tau_xy,tau_yy; - int equivstress,domaintype,damagelaw; + IssmDouble damage,tau_xx,tau_xy,tau_xz,tau_yy,tau_yz,tau_zz,stressMaxPrincipal; + int equivstress,domaintype,damagelaw,dim; /*Fetch number of vertices and allocate output*/ @@ -464,11 +467,30 @@ element->FindParam(&damagelaw,DamageLawEnum); - /*Compute stress tensor: */ + /*Get problem dimension*/ + switch(domaintype){ + case Domain2DhorizontalEnum: dim = 2; break; + case Domain3DEnum: dim = 3; break; + default: _error_("not implemented"); + } + /*Compute stress tensor and Stress Max Principal: */ element->ComputeDeviatoricStressTensor(); - + if(dim==3){ + /*Only works in 3d because the pressure is defined*/ + element->StressMaxPrincipalCreateInput(); + } /*retrieve what we need: */ Input* tau_xx_input = element->GetInput(DeviatoricStressxxEnum); _assert_(tau_xx_input); Input* tau_xy_input = element->GetInput(DeviatoricStressxyEnum); _assert_(tau_xy_input); Input* tau_yy_input = element->GetInput(DeviatoricStressyyEnum); _assert_(tau_yy_input); + Input* tau_xz_input = NULL; + Input* tau_yz_input = NULL; + Input* tau_zz_input = NULL; + Input* stressMaxPrincipal_input = NULL; + if(dim==3){ + tau_xz_input = element->GetInput(DeviatoricStressxzEnum); _assert_(tau_xz_input); + tau_yz_input = element->GetInput(DeviatoricStressyzEnum); _assert_(tau_yz_input); + tau_zz_input = element->GetInput(DeviatoricStresszzEnum); _assert_(tau_zz_input); + stressMaxPrincipal_input = element->GetInput(StressMaxPrincipalEnum); _assert_(stressMaxPrincipal_input); + } Input* damage_input = NULL; if(domaintype==Domain2DhorizontalEnum){ @@ -491,35 +513,64 @@ tau_xy_input->GetInputValue(&tau_xy,gauss); tau_yy_input->GetInputValue(&tau_yy,gauss); - + if(dim==3){ + tau_xz_input->GetInputValue(&tau_xz,gauss); + tau_yz_input->GetInputValue(&tau_yz,gauss); + tau_zz_input->GetInputValue(&tau_zz,gauss); + } /*Calculate effective stress components*/ s_xx=tau_xx/(1.-damage); s_xy=tau_xy/(1.-damage); s_yy=tau_yy/(1.-damage); - + if(dim==3){ + s_xz=tau_xz/(1.-damage); + s_yz=tau_yz/(1.-damage); + s_zz=tau_zz/(1.-damage); + } /*Calculate principal effective stresses*/ - s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2)); - s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2)); - if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;} - - if(equivstress==0){ /* von Mises */ - Chi=sqrt(s1*s1-s1*s2+s2*s2); - } - else if(equivstress==1){ /* max principal stress */ - Chi=s1; - } - Psi=Chi-stress_threshold; - NegPsi=max(-Chi,0.); /* healing only for compressive stresses */ - - if(damagelaw==1){ - PosPsi=max(Psi,0.); - f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3); - } - else if(damagelaw==2){ - PosPsi=max(Psi,1.); - f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3); - } - else _error_("damage law not supported"); - } - + if(dim==2){ + s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2)); + s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2)); + if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;} + + if(equivstress==0){ /* von Mises */ + Chi=sqrt(s1*s1-s1*s2+s2*s2); + } + else if(equivstress==1){ /* max principal stress */ + Chi=s1; + } + Psi=Chi-stress_threshold; + NegPsi=max(-Chi,0.); /* healing only for compressive stresses */ + + if(damagelaw==1){ + PosPsi=max(Psi,0.); + f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3); + } + else if(damagelaw==2){ + PosPsi=max(Psi,1.); + f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3); + } + else _error_("damage law not supported"); + } + else{ + if(equivstress==1){/* max principal stress */ + stressMaxPrincipal_input->GetInputValue(&stressMaxPrincipal,gauss); + Chi=stressMaxPrincipal/(1.-damage); + } + else if(equivstress==0){/* von Mises */ + Chi=sqrt(((s_xx-s_yy)*(s_xx-s_yy)+(s_yy-s_zz)*(s_yy-s_zz)+(s_zz-s_xx)*(s_zz-s_xx)+6.*(s_xy*s_xy+s_yz*s_yz+s_xz*s_xz))/2.); + } + Psi=Chi-stress_threshold; + NegPsi=max(-Chi,0.); /* healing only for compressive stresses */ + if(damagelaw==1){ + PosPsi=max(Psi,0.); + f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3); + } + else if(damagelaw==2){ + PosPsi=max(Psi,1.); + f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3); + } + else _error_("damage law not supported"); + } + } /*Add input*/ element->AddInput(DamageFEnum,f,element->GetElementType()); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18469) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18470) @@ -1408,5 +1408,5 @@ /*Initialize maximum eigne value*/ if(numroots>0){ - max = x[0]; + max = fabs(x[0]); } else{ @@ -1416,5 +1416,5 @@ /*Get max*/ for(int i=1;imax) max = x[i]; + if(fabs(x[i])>max) max = fabs(x[i]); } Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 18469) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 18470) @@ -94,5 +94,4 @@ analysis_enum==BalancethicknessAnalysisEnum || analysis_enum==HydrologyDCInefficientAnalysisEnum || - analysis_enum==DamageEvolutionAnalysisEnum || analysis_enum==HydrologyDCEfficientAnalysisEnum || analysis_enum==LevelsetAnalysisEnum || Index: /issm/trunk-jpl/src/c/cores/damage_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/damage_core.cpp (revision 18469) +++ /issm/trunk-jpl/src/c/cores/damage_core.cpp (revision 18470) @@ -34,5 +34,5 @@ femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); } - + /*Free resources:*/ if(numoutputs){