![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 17391
- Timestamp:
- 03/07/14 17:42:04 (11 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 4 edited
-
analyses/DamageEvolutionAnalysis.cpp (modified) (10 diffs)
-
analyses/DamageEvolutionAnalysis.h (modified) (1 diff)
-
classes/Loads/Pengrid.cpp (modified) (1 diff)
-
cores/damage_core.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
r17369 r17391 117 117 ElementMatrix* DamageEvolutionAnalysis::CreateKMatrix(Element* element){/*{{{*/ 118 118 119 /* Check if ice in element */ 120 if(!element->IsIceInElement()) return NULL; 121 122 /*Intermediaries*/ 123 Element* basalelement; 124 int meshtype,dim; 125 int stabilization; 126 IssmDouble Jdet,dt,D_scalar,h; 127 IssmDouble vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2]; 128 IssmDouble *xyz_list = NULL; 129 130 /*Get problem dimension and basal element*/ 131 element->FindParam(&meshtype,MeshTypeEnum); 132 switch(meshtype){ 133 case Mesh2DhorizontalEnum: 134 basalelement = element; 135 dim = 2; 136 break; 137 case Mesh3DEnum: 138 if(!element->IsOnBed()) return NULL; 139 basalelement = element->SpawnBasalElement(); 140 dim = 2; 141 break; 142 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet"); 143 } 144 145 /*Fetch number of nodes and dof for this finite element*/ 146 int numnodes = basalelement->GetNumberOfNodes(); 147 148 /*Initialize Element vector*/ 149 ElementMatrix* Ke = basalelement->NewElementMatrix(); 150 IssmDouble* basis = xNew<IssmDouble>(numnodes); 151 IssmDouble* B = xNew<IssmDouble>(dim*numnodes); 152 IssmDouble* Bprime = xNew<IssmDouble>(dim*numnodes); 153 IssmDouble* D = xNewZeroInit<IssmDouble>(dim*dim); 154 155 /*Retrieve all inputs and parameters*/ 156 basalelement->GetVerticesCoordinates(&xyz_list); 157 basalelement->FindParam(&dt,TimesteppingTimeStepEnum); 158 basalelement->FindParam(&stabilization,DamageStabilizationEnum); 159 Input* vxaverage_input=NULL; 160 Input* vyaverage_input=NULL; 161 if(meshtype==Mesh2DhorizontalEnum){ 162 vxaverage_input=basalelement->GetInput(VxEnum); _assert_(vxaverage_input); 163 vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input); 164 } 165 else{ 166 if(dim==1){ 167 vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input); 168 } 169 if(dim==2){ 170 vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input); 171 vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input); 172 } 173 } 174 h=basalelement->CharacteristicLength(); 175 176 /* Start looping on the number of gaussian points: */ 177 Gauss* gauss=basalelement->NewGauss(2); 178 for(int ig=gauss->begin();ig<gauss->end();ig++){ 179 gauss->GaussPoint(ig); 180 181 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); 182 basalelement->NodalFunctions(basis,gauss); 183 184 vxaverage_input->GetInputValue(&vx,gauss); 185 vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); 186 if(dim==2){ 187 vyaverage_input->GetInputValue(&vy,gauss); 188 vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); 189 } 190 191 D_scalar=gauss->weight*Jdet; 192 TripleMultiply(basis,1,numnodes,1, 193 &D_scalar,1,1,0, 194 basis,1,numnodes,0, 195 &Ke->values[0],1); 196 197 GetB(B,basalelement,dim,xyz_list,gauss); 198 GetBprime(Bprime,basalelement,dim,xyz_list,gauss); 199 200 dvxdx=dvx[0]; 201 if(dim==2) dvydy=dvy[1]; 202 D_scalar=dt*gauss->weight*Jdet; 203 204 D[0*dim+0]=D_scalar*dvxdx; 205 if(dim==2) D[1*dim+1]=D_scalar*dvydy; 206 207 TripleMultiply(B,dim,numnodes,1, 208 D,dim,dim,0, 209 B,dim,numnodes,0, 210 &Ke->values[0],1); 211 212 D[0*dim+0]=D_scalar*vx; 213 if(dim==2) D[1*dim+1]=D_scalar*vy; 214 215 TripleMultiply(B,dim,numnodes,1, 216 D,dim,dim,0, 217 Bprime,dim,numnodes,0, 218 &Ke->values[0],1); 219 220 if(stabilization==2){ 221 if(dim==1){ 222 vel=fabs(vx)+1.e-8; 223 D[0]=h/(2*vel)*vx*vx; 224 } 225 else{ 226 /*Streamline upwinding*/ 227 vel=sqrt(vx*vx+vy*vy)+1.e-8; 228 D[0*dim+0]=h/(2*vel)*vx*vx; 229 D[1*dim+0]=h/(2*vel)*vy*vx; 230 D[0*dim+1]=h/(2*vel)*vx*vy; 231 D[1*dim+1]=h/(2*vel)*vy*vy; 232 } 233 } 234 else if(stabilization==1){ 235 vxaverage_input->GetInputAverage(&vx); 236 if(dim==2) vyaverage_input->GetInputAverage(&vy); 237 D[0*dim+0]=h/2.0*fabs(vx); 238 if(dim==2) D[1*dim+1]=h/2.0*fabs(vy); 239 } 240 if(stabilization==1 || stabilization==2){ 241 if(dim==1) D[0]=D_scalar*D[0]; 242 else{ 243 D[0*dim+0]=D_scalar*D[0*dim+0]; 244 D[1*dim+0]=D_scalar*D[1*dim+0]; 245 D[0*dim+1]=D_scalar*D[0*dim+1]; 246 D[1*dim+1]=D_scalar*D[1*dim+1]; 247 } 248 249 TripleMultiply(Bprime,dim,numnodes,1, 250 D,dim,dim,0, 251 Bprime,dim,numnodes,0, 252 &Ke->values[0],1); 253 } 254 255 } 256 257 /*Clean up and return*/ 258 xDelete<IssmDouble>(xyz_list); 259 xDelete<IssmDouble>(basis); 260 xDelete<IssmDouble>(B); 261 xDelete<IssmDouble>(Bprime); 262 xDelete<IssmDouble>(D); 263 delete gauss; 264 if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; 265 return Ke; 266 }/*}}}*/ 267 ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/ 268 269 /* Check if ice in element */ 270 if(!element->IsIceInElement()) return NULL; 271 119 272 /*Intermediaries*/ 120 273 int meshtype; 121 274 Element* basalelement; 275 IssmDouble Jdet,dt; 276 IssmDouble f,damage; 277 IssmDouble* xyz_list = NULL; 122 278 123 279 /*Get basal element*/ … … 134 290 } 135 291 136 /*Intermediaries */137 int stabilization;138 IssmDouble Jdet,dt,D_scalar;139 IssmDouble vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2];140 IssmDouble *xyz_list = NULL;141 142 292 /*Fetch number of nodes and dof for this finite element*/ 143 293 int numnodes = basalelement->GetNumberOfNodes(); 144 294 145 /*Initialize Element vector*/ 146 ElementMatrix* Ke = basalelement->NewElementMatrix(NoneApproximationEnum); 147 IssmDouble* basis = xNew<IssmDouble>(numnodes); 148 IssmDouble* B = xNew<IssmDouble>(2*numnodes); 149 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes); 150 IssmDouble D[2][2]={0.}; 295 /*Initialize Element vector and other vectors*/ 296 ElementVector* pe = basalelement->NewElementVector(); 297 IssmDouble* basis = xNew<IssmDouble>(numnodes); 151 298 152 299 /*Retrieve all inputs and parameters*/ 153 300 basalelement->GetVerticesCoordinates(&xyz_list); 154 301 basalelement->FindParam(&dt,TimesteppingTimeStepEnum); 155 basalelement->FindParam(&stabilization,DamageStabilizationEnum); 156 Input* vxaverage_input=NULL; 157 Input* vyaverage_input=NULL; 158 if(meshtype==Mesh2DhorizontalEnum){ 159 vxaverage_input=basalelement->GetInput(VxEnum); _assert_(vxaverage_input); 160 vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input); 161 } 162 else{ 163 vxaverage_input=basalelement->GetInput(VxAverageEnum); _assert_(vxaverage_input); 164 vyaverage_input=basalelement->GetInput(VyAverageEnum); _assert_(vyaverage_input); 165 } 166 IssmDouble h=basalelement->CharacteristicLength(); 302 this->CreateDamageFInput(basalelement); 303 Input* damaged_input = basalelement->GetInput(DamageDEnum); _assert_(damaged_input); 304 Input* damagef_input = basalelement->GetInput(DamageFEnum); _assert_(damagef_input); 305 167 306 168 307 /* Start looping on the number of gaussian points: */ … … 173 312 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); 174 313 basalelement->NodalFunctions(basis,gauss); 175 D_scalar=gauss->weight*Jdet;176 177 vxaverage_input->GetInputValue(&vx,gauss);178 vyaverage_input->GetInputValue(&vy,gauss);179 vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);180 vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);181 182 TripleMultiply(basis,1,numnodes,1,183 &D_scalar,1,1,0,184 basis,1,numnodes,0,185 &Ke->values[0],1);186 187 GetB(B,element,xyz_list,gauss);188 GetBprime(Bprime,element,xyz_list,gauss);189 190 dvxdx=dvx[0];191 dvydy=dvy[1];192 D_scalar=dt*gauss->weight*Jdet;193 194 D[0][0]=D_scalar*dvxdx;195 D[1][1]=D_scalar*dvydy;196 TripleMultiply(B,2,numnodes,1,197 &D[0][0],2,2,0,198 B,2,numnodes,0,199 &Ke->values[0],1);200 201 D[0][0]=D_scalar*vx;202 D[1][1]=D_scalar*vy;203 TripleMultiply(B,2,numnodes,1,204 &D[0][0],2,2,0,205 Bprime,2,numnodes,0,206 &Ke->values[0],1);207 208 if(stabilization==2){209 /*Streamline upwinding*/210 vel=sqrt(vx*vx+vy*vy)+1.e-8;211 D[0][0]=h/(2*vel)*vx*vx;212 D[1][0]=h/(2*vel)*vy*vx;213 D[0][1]=h/(2*vel)*vx*vy;214 D[1][1]=h/(2*vel)*vy*vy;215 }216 else if(stabilization==1){217 /*SSA*/218 vxaverage_input->GetInputAverage(&vx);219 vyaverage_input->GetInputAverage(&vy);220 D[0][0]=h/2.0*fabs(vx);221 D[0][1]=0.;222 D[1][0]=0.;223 D[1][1]=h/2.0*fabs(vy);224 }225 if(stabilization==1 || stabilization==2){226 D[0][0]=D_scalar*D[0][0];227 D[1][0]=D_scalar*D[1][0];228 D[0][1]=D_scalar*D[0][1];229 D[1][1]=D_scalar*D[1][1];230 TripleMultiply(Bprime,2,numnodes,1,231 &D[0][0],2,2,0,232 Bprime,2,numnodes,0,233 &Ke->values[0],1);234 }235 236 }237 238 /*Clean up and return*/239 xDelete<IssmDouble>(xyz_list);240 xDelete<IssmDouble>(basis);241 xDelete<IssmDouble>(B);242 xDelete<IssmDouble>(Bprime);243 delete gauss;244 if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};245 return Ke;246 }/*}}}*/247 ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/248 249 /*Intermediaries*/250 int meshtype;251 Element* basalelement;252 253 /*Get basal element*/254 element->FindParam(&meshtype,MeshTypeEnum);255 switch(meshtype){256 case Mesh2DhorizontalEnum:257 basalelement = element;258 break;259 case Mesh3DEnum:260 if(!element->IsOnBed()) return NULL;261 basalelement = element->SpawnBasalElement();262 break;263 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");264 }265 266 /*Intermediaries */267 IssmDouble Jdet,dt;268 IssmDouble f,damage;269 IssmDouble* xyz_list = NULL;270 271 /*Fetch number of nodes and dof for this finite element*/272 int numnodes = element->GetNumberOfNodes();273 274 /*Initialize Element vector and other vectors*/275 ElementVector* pe = element->NewElementVector();276 IssmDouble* basis = xNew<IssmDouble>(numnodes);277 278 /*Retrieve all inputs and parameters*/279 element->GetVerticesCoordinates(&xyz_list);280 element->FindParam(&dt,TimesteppingTimeStepEnum);281 this->CreateDamageFInput(element);282 Input* damaged_input = element->GetInput(DamageDbarEnum); _assert_(damaged_input);283 Input* damagef_input = element->GetInput(DamageFEnum); _assert_(damagef_input);284 285 /* Start looping on the number of gaussian points: */286 Gauss* gauss=element->NewGauss(2);287 for(int ig=gauss->begin();ig<gauss->end();ig++){288 gauss->GaussPoint(ig);289 290 element->JacobianDeterminant(&Jdet,xyz_list,gauss);291 element->NodalFunctions(basis,gauss);292 314 293 315 damaged_input->GetInputValue(&damage,gauss); 294 316 damagef_input->GetInputValue(&f,gauss); 295 317 296 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i]; 318 for(int i=0;i<numnodes;i++){ 319 pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i]; 320 } 297 321 } 298 322 … … 304 328 return pe; 305 329 }/*}}}*/ 306 void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/330 void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 307 331 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 308 332 * For node i, Bi can be expressed in the actual coordinate system … … 324 348 /*Build B: */ 325 349 for(int i=0;i<numnodes;i++){ 326 B[numnodes*0+i] = basis[i]; 327 B[numnodes*1+i] = basis[i]; 350 for(int j=0;j<dim;j++){ 351 B[numnodes*j+i] = basis[i]; 352 } 328 353 } 329 354 … … 331 356 xDelete<IssmDouble>(basis); 332 357 }/*}}}*/ 333 void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/358 void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ 334 359 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 335 360 * For node i, Bi' can be expressed in the actual coordinate system … … 346 371 347 372 /*Get nodal functions derivatives*/ 348 IssmDouble* dbasis=xNew<IssmDouble>( 2*numnodes);373 IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes); 349 374 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); 350 375 351 376 /*Build B': */ 352 377 for(int i=0;i<numnodes;i++){ 353 Bprime[numnodes*0+i] = dbasis[0*numnodes+i]; 354 Bprime[numnodes*1+i] = dbasis[1*numnodes+i]; 378 for(int j=0;j<dim;j++){ 379 Bprime[numnodes*j+i] = dbasis[j*numnodes+i]; 380 } 355 381 } 356 382 … … 364 390 void DamageEvolutionAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ 365 391 392 int meshtype; 366 393 IssmDouble max_damage; 367 394 int *doflist = NULL; 368 395 Element* basalelement=NULL; 396 397 element->FindParam(&meshtype,MeshTypeEnum); 398 if(meshtype!=Mesh2DhorizontalEnum){ 399 if(!element->IsOnBed()) return; 400 basalelement=element->SpawnBasalElement(); 401 } 402 else{ 403 basalelement = element; 404 } 369 405 /*Fetch number of nodes and dof for this finite element*/ 370 int numnodes = element->GetNumberOfNodes();406 int numnodes = basalelement->GetNumberOfNodes(); 371 407 372 408 /*Fetch dof list and allocate solution vector*/ 373 element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);374 IssmDouble* values= xNew<IssmDouble>(numnodes);409 basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); 410 IssmDouble* newdamage = xNew<IssmDouble>(numnodes); 375 411 376 412 /*Get user-supplied max_damage: */ 377 element->FindParam(&max_damage,DamageMaxDamageEnum);413 basalelement->FindParam(&max_damage,DamageMaxDamageEnum); 378 414 379 415 /*Use the dof list to index into the solution vector: */ 380 416 for(int i=0;i<numnodes;i++){ 381 values[i]=solution[doflist[i]];417 newdamage[i]=solution[doflist[i]]; 382 418 /*Check solution*/ 383 if(xIsNan<IssmDouble>( values[i])) _error_("NaN found in solution vector");419 if(xIsNan<IssmDouble>(newdamage[i])) _error_("NaN found in solution vector"); 384 420 /*Enforce D < max_damage and D > 0 */ 385 if( values[i]>max_damage) values[i]=max_damage;386 else if( values[i]<0.) values[i]=0.;421 if(newdamage[i]>max_damage) newdamage[i]=max_damage; 422 else if(newdamage[i]<0.) newdamage[i]=0.; 387 423 } 388 424 389 425 /*Get all inputs and parameters*/ 390 element->AddInput(DamageDbarEnum,values,P1Enum);426 basalelement->AddBasalInput(DamageDEnum,newdamage,P1Enum); 391 427 392 428 /*Free ressources:*/ 393 xDelete<IssmDouble>( values);429 xDelete<IssmDouble>(newdamage); 394 430 xDelete<int>(doflist); 431 if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; 395 432 }/*}}}*/ 396 433 void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ … … 427 464 Input* sigma_xy_input = element->GetInput(StressTensorxyEnum); _assert_(sigma_xy_input); 428 465 Input* sigma_yy_input = element->GetInput(StressTensoryyEnum); _assert_(sigma_yy_input); 429 Input* damage_input = element->GetInput(DamageD barEnum); _assert_(damage_input);466 Input* damage_input = element->GetInput(DamageDEnum); _assert_(damage_input); 430 467 431 468 /*retrieve the desired type of equivalent stress*/ … … 446 483 s_yy=sigma_yy/(1.-damage); 447 484 448 if(equivstress==1){ /* von Mises */ 449 J2s=1./sqrt(2.)*sqrt(s_xx*s_xx + s_yy*s_yy + s_xy*s_xy); 450 Chi=sqrt(3.0)*J2s; 485 if(equivstress==0){ /* von Mises */ 486 Chi=sqrt(s_xx*s_xx - s_xx*s_yy + s_yy*s_yy + 3*s_xy*s_xy); 451 487 } 452 488 Psi=Chi-stress_threshold; 453 489 PosPsi=max(Psi,0.); 454 NegPsi=max(- Psi,0.);490 NegPsi=max(-Chi,0.); /* healing only for compressive stresses */ 455 491 456 492 f[iv]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1.-damage),-c3); -
issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h
r17212 r17391 26 26 ElementMatrix* CreateKMatrix(Element* element); 27 27 ElementVector* CreatePVector(Element* element); 28 void GetB(IssmDouble* B,Element* element, IssmDouble* xyz_list,Gauss* gauss);29 void GetBprime(IssmDouble* B,Element* element, IssmDouble* xyz_list,Gauss* gauss);28 void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); 29 void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); 30 30 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element); 31 31 void InputUpdateFromSolution(IssmDouble* solution,Element* element); -
issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
r17372 r17391 713 713 714 714 //First recover damage using the element: */ 715 element->GetInputValue(&damage,node,DamageD barEnum);715 element->GetInputValue(&damage,node,DamageDEnum); 716 716 717 717 //Recover our data: -
issm/trunk-jpl/src/c/cores/damage_core.cpp
r17369 r17391 19 19 char **requested_outputs = NULL; 20 20 21 if(VerboseSolution()) _printf0_(" computing damage\n");22 21 23 22 //first recover parameters common to all solutions … … 33 32 } 34 33 34 if(VerboseSolution()) _printf0_(" computing damage\n"); 35 35 femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum); 36 36 solutionsequence_linear(femmodel);
Note:
See TracChangeset
for help on using the changeset viewer.
