Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h =================================================================== --- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 17390) +++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 17391) @@ -25,8 +25,8 @@ ElementMatrix* CreateJacobianMatrix(Element* element); 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 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 InputUpdateFromSolution(IssmDouble* solution,Element* element); void UpdateConstraints(FemModel* femmodel); Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17390) +++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17391) @@ -116,38 +116,41 @@ }/*}}}*/ ElementMatrix* DamageEvolutionAnalysis::CreateKMatrix(Element* element){/*{{{*/ + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + /*Intermediaries*/ - int meshtype; - Element* basalelement; + Element* basalelement; + int meshtype,dim; + int stabilization; + IssmDouble Jdet,dt,D_scalar,h; + IssmDouble vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2]; + IssmDouble *xyz_list = NULL; - /*Get basal element*/ + /*Get problem dimension and basal element*/ element->FindParam(&meshtype,MeshTypeEnum); switch(meshtype){ case Mesh2DhorizontalEnum: basalelement = element; + dim = 2; break; case Mesh3DEnum: if(!element->IsOnBed()) return NULL; basalelement = element->SpawnBasalElement(); + dim = 2; break; default: _error_("mesh "<GetNumberOfNodes(); /*Initialize Element vector*/ - ElementMatrix* Ke = basalelement->NewElementMatrix(NoneApproximationEnum); + ElementMatrix* Ke = basalelement->NewElementMatrix(); IssmDouble* basis = xNew(numnodes); - IssmDouble* B = xNew(2*numnodes); - IssmDouble* Bprime = xNew(2*numnodes); - IssmDouble D[2][2]={0.}; + IssmDouble* B = xNew(dim*numnodes); + IssmDouble* Bprime = xNew(dim*numnodes); + IssmDouble* D = xNewZeroInit(dim*dim); /*Retrieve all inputs and parameters*/ basalelement->GetVerticesCoordinates(&xyz_list); @@ -160,10 +163,15 @@ 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); + 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); + } } - IssmDouble h=basalelement->CharacteristicLength(); + h=basalelement->CharacteristicLength(); /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(2); @@ -172,64 +180,75 @@ 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); + if(dim==2){ + vyaverage_input->GetInputValue(&vy,gauss); + vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + } + D_scalar=gauss->weight*Jdet; 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); + GetB(B,basalelement,dim,xyz_list,gauss); + GetBprime(Bprime,basalelement,dim,xyz_list,gauss); dvxdx=dvx[0]; - dvydy=dvy[1]; + if(dim==2) 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, + D[0*dim+0]=D_scalar*dvxdx; + if(dim==2) D[1*dim+1]=D_scalar*dvydy; + + TripleMultiply(B,dim,numnodes,1, + D,dim,dim,0, + B,dim,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, + D[0*dim+0]=D_scalar*vx; + if(dim==2) D[1*dim+1]=D_scalar*vy; + + TripleMultiply(B,dim,numnodes,1, + D,dim,dim,0, + Bprime,dim,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; + if(dim==1){ + vel=fabs(vx)+1.e-8; + D[0]=h/(2*vel)*vx*vx; + } + else{ + /*Streamline upwinding*/ + vel=sqrt(vx*vx+vy*vy)+1.e-8; + D[0*dim+0]=h/(2*vel)*vx*vx; + D[1*dim+0]=h/(2*vel)*vy*vx; + D[0*dim+1]=h/(2*vel)*vx*vy; + D[1*dim+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(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(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, + if(dim==1) D[0]=D_scalar*D[0]; + else{ + D[0*dim+0]=D_scalar*D[0*dim+0]; + 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]; + } + + TripleMultiply(Bprime,dim,numnodes,1, + D,dim,dim,0, + Bprime,dim,numnodes,0, &Ke->values[0],1); } @@ -240,15 +259,22 @@ xDelete(basis); xDelete(B); xDelete(Bprime); + xDelete(D); delete gauss; if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; return Ke; }/*}}}*/ ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/ + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + /*Intermediaries*/ int meshtype; Element* basalelement; + IssmDouble Jdet,dt; + IssmDouble f,damage; + IssmDouble* xyz_list = NULL; /*Get basal element*/ element->FindParam(&meshtype,MeshTypeEnum); @@ -263,37 +289,35 @@ default: _error_("mesh "<GetNumberOfNodes(); + int numnodes = basalelement->GetNumberOfNodes(); /*Initialize Element vector and other vectors*/ - ElementVector* pe = element->NewElementVector(); + ElementVector* pe = basalelement->NewElementVector(); IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - element->FindParam(&dt,TimesteppingTimeStepEnum); - this->CreateDamageFInput(element); - Input* damaged_input = element->GetInput(DamageDbarEnum); _assert_(damaged_input); - Input* damagef_input = element->GetInput(DamageFEnum); _assert_(damagef_input); + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->FindParam(&dt,TimesteppingTimeStepEnum); + this->CreateDamageFInput(basalelement); + Input* damaged_input = basalelement->GetInput(DamageDEnum); _assert_(damaged_input); + Input* damagef_input = basalelement->GetInput(DamageFEnum); _assert_(damagef_input); + /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(2); + Gauss* gauss=basalelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - element->NodalFunctions(basis,gauss); + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctions(basis,gauss); damaged_input->GetInputValue(&damage,gauss); damagef_input->GetInputValue(&f,gauss); - for(int i=0;ivalues[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i]; + for(int i=0;ivalues[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i]; + } } /*Clean up and return*/ @@ -303,7 +327,7 @@ delete gauss; return pe; }/*}}}*/ -void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +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*NDOF2. * For node i, Bi can be expressed in the actual coordinate system * by: @@ -323,14 +347,15 @@ /*Build B: */ for(int i=0;i(basis); }/*}}}*/ -void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +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*NDOF2. * For node i, Bi' can be expressed in the actual coordinate system * by: @@ -345,13 +370,14 @@ int numnodes = element->GetNumberOfNodes(); /*Get nodal functions derivatives*/ - IssmDouble* dbasis=xNew(2*numnodes); + IssmDouble* dbasis=xNew(dim*numnodes); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B': */ for(int i=0;iFindParam(&meshtype,MeshTypeEnum); + if(meshtype!=Mesh2DhorizontalEnum){ + if(!element->IsOnBed()) return; + basalelement=element->SpawnBasalElement(); + } + else{ + basalelement = element; + } /*Fetch number of nodes and dof for this finite element*/ - int numnodes = element->GetNumberOfNodes(); + int numnodes = basalelement->GetNumberOfNodes(); /*Fetch dof list and allocate solution vector*/ - element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - IssmDouble* values = xNew(numnodes); + basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* newdamage = xNew(numnodes); /*Get user-supplied max_damage: */ - element->FindParam(&max_damage,DamageMaxDamageEnum); + basalelement->FindParam(&max_damage,DamageMaxDamageEnum); /*Use the dof list to index into the solution vector: */ for(int i=0;i(values[i])) _error_("NaN found in solution vector"); + if(xIsNan(newdamage[i])) _error_("NaN found in solution vector"); /*Enforce D < max_damage and D > 0 */ - if(values[i]>max_damage) values[i]=max_damage; - else if(values[i]<0.) values[i]=0.; + if(newdamage[i]>max_damage) newdamage[i]=max_damage; + else if(newdamage[i]<0.) newdamage[i]=0.; } /*Get all inputs and parameters*/ - element->AddInput(DamageDbarEnum,values,P1Enum); + basalelement->AddBasalInput(DamageDEnum,newdamage,P1Enum); /*Free ressources:*/ - xDelete(values); + xDelete(newdamage); xDelete(doflist); + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ @@ -426,7 +463,7 @@ Input* sigma_xx_input = element->GetInput(StressTensorxxEnum); _assert_(sigma_xx_input); Input* sigma_xy_input = element->GetInput(StressTensorxyEnum); _assert_(sigma_xy_input); Input* sigma_yy_input = element->GetInput(StressTensoryyEnum); _assert_(sigma_yy_input); - Input* damage_input = element->GetInput(DamageDbarEnum); _assert_(damage_input); + Input* damage_input = element->GetInput(DamageDEnum); _assert_(damage_input); /*retrieve the desired type of equivalent stress*/ element->FindParam(&equivstress,DamageEquivStressEnum); @@ -445,13 +482,12 @@ s_xy=sigma_xy/(1.-damage); s_yy=sigma_yy/(1.-damage); - if(equivstress==1){ /* von Mises */ - J2s=1./sqrt(2.)*sqrt(s_xx*s_xx + s_yy*s_yy + s_xy*s_xy); - Chi=sqrt(3.0)*J2s; + if(equivstress==0){ /* von Mises */ + Chi=sqrt(s_xx*s_xx - s_xx*s_yy + s_yy*s_yy + 3*s_xy*s_xy); } Psi=Chi-stress_threshold; PosPsi=max(Psi,0.); - NegPsi=max(-Psi,0.); + NegPsi=max(-Chi,0.); /* healing only for compressive stresses */ f[iv]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1.-damage),-c3); } Index: ../trunk-jpl/src/c/cores/damage_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/damage_core.cpp (revision 17390) +++ ../trunk-jpl/src/c/cores/damage_core.cpp (revision 17391) @@ -18,7 +18,6 @@ int numoutputs = 0; char **requested_outputs = NULL; - if(VerboseSolution()) _printf0_(" computing damage\n"); //first recover parameters common to all solutions femmodel->parameters->FindParam(&save_results,SaveResultsEnum); @@ -32,6 +31,7 @@ ResetConstraintsx(femmodel); } + if(VerboseSolution()) _printf0_(" computing damage\n"); femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum); solutionsequence_linear(femmodel); Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17390) +++ ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17391) @@ -712,7 +712,7 @@ } //First recover damage using the element: */ - element->GetInputValue(&damage,node,DamageDbarEnum); + element->GetInputValue(&damage,node,DamageDEnum); //Recover our data: parameters->FindParam(&penalty_lock,DamagePenaltyLockEnum);