![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 16993 for issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
- Timestamp:
- 12/03/13 10:36:26 (11 years ago)
- File:
-
- 1 edited
-
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r16982 r16993 191 191 } 192 192 /*}}}*/ 193 /*FUNCTION Tria::CreateKMatrix(void) {{{*/194 ElementMatrix* Tria::CreateKMatrix(void){195 196 /*retreive parameters: */197 int analysis_type;198 int meshtype;199 parameters->FindParam(&analysis_type,AnalysisTypeEnum);200 201 /*Checks in debugging mode{{{*/202 _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs);203 /*}}}*/204 205 /*Skip if water element*/206 if(NoIceInElement()) return NULL;207 208 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */209 switch(analysis_type){210 #ifdef _HAVE_STRESSBALANCE_211 case StressbalanceAnalysisEnum:212 int approximation;213 inputs->GetInputValue(&approximation,ApproximationEnum);214 switch(approximation){215 case SSAApproximationEnum:216 return CreateKMatrixStressbalanceSSA();217 case FSApproximationEnum:218 return CreateKMatrixStressbalanceFS();219 case SIAApproximationEnum:220 return NULL;221 case NoneApproximationEnum:222 return NULL;223 default:224 _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");225 }226 break;227 case StressbalanceSIAAnalysisEnum:228 return CreateKMatrixStressbalanceSIA();229 break;230 #endif231 #ifdef _HAVE_DAMAGE_232 case DamageEvolutionAnalysisEnum:233 return CreateKMatrixDamageEvolution();234 break;235 #endif236 case L2ProjectionBaseAnalysisEnum:237 parameters->FindParam(&meshtype,MeshTypeEnum);238 if(meshtype==Mesh2DverticalEnum){239 return CreateBasalMassMatrix();240 }241 else{242 return CreateMassMatrix();243 }244 break;245 #ifdef _HAVE_MASSTRANSPORT_246 case MasstransportAnalysisEnum:247 return CreateKMatrixMasstransport();248 break;249 case FreeSurfaceTopAnalysisEnum:250 parameters->FindParam(&meshtype,MeshTypeEnum);251 switch(meshtype){252 case Mesh2DverticalEnum:253 return CreateKMatrixFreeSurfaceTop1D();254 case Mesh3DEnum:255 return CreateKMatrixFreeSurfaceTop();256 default:257 _error_("Mesh not supported yet");258 }259 break;260 case FreeSurfaceBaseAnalysisEnum:261 parameters->FindParam(&meshtype,MeshTypeEnum);262 switch(meshtype){263 case Mesh2DverticalEnum:264 return CreateKMatrixFreeSurfaceBase1D();265 case Mesh3DEnum:266 return CreateKMatrixFreeSurfaceBase();267 default:268 _error_("Mesh not supported yet");269 }270 break;271 case ExtrudeFromBaseAnalysisEnum: case ExtrudeFromTopAnalysisEnum:272 return CreateKMatrixExtrusion();273 #endif274 #ifdef _HAVE_HYDROLOGY_275 case HydrologyShreveAnalysisEnum:276 return CreateKMatrixHydrologyShreve();277 break;278 case HydrologyDCInefficientAnalysisEnum:279 return CreateKMatrixHydrologyDCInefficient();280 break;281 case HydrologyDCEfficientAnalysisEnum:282 return CreateKMatrixHydrologyDCEfficient();283 break;284 case L2ProjectionEPLAnalysisEnum:285 return CreateEPLDomainMassMatrix();286 break;287 #endif288 #ifdef _HAVE_BALANCED_289 case BalancethicknessAnalysisEnum:290 return CreateKMatrixBalancethickness();291 break;292 case BalancevelocityAnalysisEnum:293 return CreateKMatrixBalancevelocity();294 break;295 case SmoothedSurfaceSlopeXAnalysisEnum: case SmoothedSurfaceSlopeYAnalysisEnum:296 return CreateKMatrixSmoothedSlope();297 break;298 #endif299 #ifdef _HAVE_CONTROL_300 case AdjointBalancethicknessAnalysisEnum:301 return CreateKMatrixAdjointBalancethickness();302 break;303 case AdjointHorizAnalysisEnum:304 return CreateKMatrixAdjointSSA();305 break;306 #endif307 default:308 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");309 }310 311 /*Make compiler happy*/312 return NULL;313 }314 /*}}}*/315 /*FUNCTION Tria::CreateMassMatrix {{{*/316 ElementMatrix* Tria::CreateMassMatrix(void){317 318 /* Intermediaries */319 IssmDouble D,Jdet;320 IssmDouble xyz_list[NUMVERTICES][3];321 322 /*Fetch number of nodes and dof for this finite element*/323 int numnodes = this->NumberofNodes();324 325 /*Initialize Element matrix and vectors*/326 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);327 IssmDouble* basis = xNew<IssmDouble>(numnodes);328 329 /*Retrieve all inputs and parameters*/330 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);331 332 /* Start looping on the number of gaussian points: */333 GaussTria* gauss=new GaussTria(2);334 for(int ig=gauss->begin();ig<gauss->end();ig++){335 336 gauss->GaussPoint(ig);337 338 GetNodalFunctions(basis,gauss);339 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);340 D=gauss->weight*Jdet;341 342 TripleMultiply(basis,1,numnodes,1,343 &D,1,1,0,344 basis,1,numnodes,0,345 &Ke->values[0],1);346 }347 348 /*Clean up and return*/349 delete gauss;350 xDelete<IssmDouble>(basis);351 return Ke;352 }353 /*}}}*/354 /*FUNCTION Tria::CreateBasalMassMatrix{{{*/355 ElementMatrix* Tria::CreateBasalMassMatrix(void){356 357 if(!HasEdgeOnBed()) return NULL;358 359 int index1,index2;360 this->EdgeOnBedIndices(&index1,&index2);361 362 Seg* seg=(Seg*)SpawnSeg(index1,index2);363 ElementMatrix* Ke=seg->CreateMassMatrix();364 delete seg->material; delete seg;365 366 /*clean up and return*/367 return Ke;368 }369 /*}}}*/370 193 /*FUNCTION Tria::CreateDVector {{{*/ 371 194 void Tria::CreateDVector(Vector<IssmDouble>* df){ 372 195 373 196 /*Nothing done yet*/ 374 }375 /*}}}*/376 /*FUNCTION Tria::CreateJacobianMatrix{{{*/377 void Tria::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){378 379 /*retrieve parameters: */380 ElementMatrix* Ke=NULL;381 int analysis_type;382 parameters->FindParam(&analysis_type,AnalysisTypeEnum);383 384 /*Checks in debugging {{{*/385 _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs);386 /*}}}*/387 388 /*Skip if water element*/389 if(NoIceInElement()) return;390 391 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */392 switch(analysis_type){393 #ifdef _HAVE_STRESSBALANCE_394 case StressbalanceAnalysisEnum:395 Ke=CreateJacobianStressbalanceSSA();396 break;397 #endif398 default:399 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");400 }401 402 /*Add to global matrix*/403 if(Ke){404 Ke->AddToGlobal(Jff);405 delete Ke;406 }407 197 } 408 198 /*}}}*/ … … 3273 3063 3274 3064 #ifdef _HAVE_STRESSBALANCE_ 3275 /*FUNCTION Tria::CreateKMatrixStressbalanceFS{{{*/3276 ElementMatrix* Tria::CreateKMatrixStressbalanceFS(void){3277 3278 ElementMatrix* Ke1 = NULL;3279 ElementMatrix* Ke2 = NULL;3280 ElementMatrix* Ke = NULL;3281 3282 /*compute all stiffness matrices for this element*/3283 Ke1=CreateKMatrixStressbalanceFSViscous();3284 Ke2=CreateKMatrixStressbalanceFSFriction();3285 Ke =new ElementMatrix(Ke1,Ke2);3286 3287 /*clean-up and return*/3288 delete Ke1;3289 delete Ke2;3290 return Ke;3291 3292 }3293 /*}}}*/3294 /*FUNCTION Tria::CreateKMatrixStressbalanceFSViscous {{{*/3295 ElementMatrix* Tria::CreateKMatrixStressbalanceFSViscous(void){3296 3297 /*Intermediaries */3298 int i,approximation;3299 IssmDouble Jdet,viscosity,FSreconditioning,D_scalar;3300 IssmDouble xyz_list[NUMVERTICES][3];3301 IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/3302 GaussTria *gauss=NULL;3303 3304 /*Fetch number of nodes and dof for this finite element*/3305 int vnumnodes = this->NumberofNodesVelocity();3306 int pnumnodes = this->NumberofNodesPressure();3307 int numdof = vnumnodes*NDOF2 + pnumnodes*NDOF1;3308 3309 /*Prepare coordinate system list*/3310 int* cs_list = xNew<int>(vnumnodes+pnumnodes);3311 for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;3312 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;3313 3314 /*Initialize Element matrix and vectors*/3315 ElementMatrix* Ke = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);3316 IssmDouble* B = xNew<IssmDouble>(5*numdof);3317 IssmDouble* Bprime = xNew<IssmDouble>(5*numdof);3318 IssmDouble* D = xNewZeroInit<IssmDouble>(5*5);3319 3320 /*Retrieve all inputs and parameters*/3321 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);3322 parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);3323 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);3324 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);3325 3326 /* Start looping on the number of gaussian points: */3327 gauss=new GaussTria(5);3328 for(int ig=gauss->begin();ig<gauss->end();ig++){3329 3330 gauss->GaussPoint(ig);3331 3332 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);3333 GetBFS(B,&xyz_list[0][0],gauss);3334 GetBprimeFS(Bprime,&xyz_list[0][0],gauss);3335 3336 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);3337 material->GetViscosity2dvertical(&viscosity,&epsilon[0]);3338 3339 D_scalar=gauss->weight*Jdet;3340 for(i=0;i<3;i++) D[i*5+i] = +D_scalar*2.*viscosity;3341 for(i=3;i<5;i++) D[i*5+i] = -D_scalar*FSreconditioning;3342 3343 TripleMultiply(B,5,numdof,1,3344 D,5,5,0,3345 Bprime,5,numdof,0,3346 Ke->values,1);3347 }3348 3349 /*Transform Coordinate System*/3350 ::TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);3351 3352 /*Clean up and return*/3353 delete gauss;3354 xDelete<int>(cs_list);3355 xDelete<IssmDouble>(B);3356 xDelete<IssmDouble>(Bprime);3357 xDelete<IssmDouble>(D);3358 return Ke;3359 }3360 /*}}}*/3361 /*FUNCTION Tria::CreateKMatrixStressbalanceFSFriction{{{*/3362 ElementMatrix* Tria::CreateKMatrixStressbalanceFSFriction(void){3363 3364 /*Intermediaries */3365 int i,j;3366 int analysis_type;3367 int indices[2];3368 IssmDouble alpha2,Jdet;3369 IssmDouble xyz_list[NUMVERTICES][3];3370 IssmDouble xyz_list_seg[NUMVERTICES1D][3];3371 Friction *friction = NULL;3372 3373 /*Initialize Element matrix and return if necessary*/3374 if(IsFloating() || !HasEdgeOnBed()) return NULL;3375 3376 /*Fetch number of nodes and dof for this finite element*/3377 int vnumnodes = this->NumberofNodesVelocity();3378 int pnumnodes = this->NumberofNodesPressure();3379 int numdof = vnumnodes*NDOF2 + pnumnodes*NDOF1;3380 3381 /*Prepare coordinate system list*/3382 int* cs_list = xNew<int>(vnumnodes+pnumnodes);3383 for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;3384 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;3385 3386 /*Initialize Element matrix and vectors*/3387 ElementMatrix* Ke = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);3388 IssmDouble* BFriction = xNew<IssmDouble>(1*numdof);3389 IssmDouble D;3390 3391 /*Retrieve all inputs and parameters*/3392 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);3393 parameters->FindParam(&analysis_type,AnalysisTypeEnum);3394 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);3395 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);3396 3397 /*Get vertex indices that lie on bed*/3398 this->EdgeOnBedIndices(&indices[0],&indices[1]);3399 for(i=0;i<NUMVERTICES1D;i++) for(j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];3400 3401 /*build friction object, used later on: */3402 friction=new Friction(this,1);3403 3404 /* Start looping on the number of gauss 1d (nodes on the bedrock) */3405 GaussTria* gauss=new GaussTria(indices[0],indices[1],3);3406 for(int ig=gauss->begin();ig<gauss->end();ig++){3407 3408 gauss->GaussPoint(ig);3409 3410 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);3411 GetLFS(BFriction,gauss);3412 3413 friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,NULL);3414 D = +alpha2*gauss->weight*Jdet; //taub_x = -alpha2 vx3415 3416 TripleMultiply(BFriction,1,numdof,1,3417 &D,1,1,0,3418 BFriction,1,numdof,0,3419 Ke->values,1);3420 }3421 3422 /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/3423 //::TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);3424 3425 /*Clean up and return*/3426 delete gauss;3427 delete friction;3428 xDelete<IssmDouble>(BFriction);3429 xDelete<int>(cs_list);3430 return Ke;3431 }3432 /*}}}*/3433 /*FUNCTION Tria::CreateKMatrixStressbalanceSSA {{{*/3434 ElementMatrix* Tria::CreateKMatrixStressbalanceSSA(void){3435 3436 /*compute all stiffness matrices for this element*/3437 ElementMatrix* Ke1=CreateKMatrixStressbalanceSSAViscous();3438 ElementMatrix* Ke2=CreateKMatrixStressbalanceSSAFriction();3439 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);3440 3441 /*clean-up and return*/3442 delete Ke1;3443 delete Ke2;3444 return Ke;3445 }3446 /*}}}*/3447 /*FUNCTION Tria::CreateKMatrixStressbalanceSSAViscous{{{*/3448 ElementMatrix* Tria::CreateKMatrixStressbalanceSSAViscous(void){3449 3450 /*Intermediaries*/3451 IssmDouble xyz_list[NUMVERTICES][3];3452 IssmDouble viscosity,newviscosity,oldviscosity;3453 IssmDouble viscosity_overshoot,thickness,Jdet;3454 IssmDouble epsilon[3],oldepsilon[3]; /* epsilon=[exx,eyy,exy]; */3455 IssmDouble D_scalar;3456 3457 /*Fetch number of nodes and dof for this finite element*/3458 int numnodes = this->NumberofNodes();3459 int numdof = numnodes*NDOF2;3460 3461 /*Initialize Element matrix and vectors*/3462 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,SSAApproximationEnum);3463 IssmDouble* B = xNew<IssmDouble>(3*numdof);3464 IssmDouble* Bprime = xNew<IssmDouble>(3*numdof);3465 IssmDouble* D = xNewZeroInit<IssmDouble>(3*3);3466 3467 /*Retrieve all inputs and parameters*/3468 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);3469 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);3470 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);3471 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);3472 Input* vxold_input=inputs->GetInput(VxPicardEnum); _assert_(vxold_input);3473 Input* vyold_input=inputs->GetInput(VyPicardEnum); _assert_(vyold_input);3474 this->parameters->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);3475 3476 /* Start looping on the number of gaussian points: */3477 GaussTria* gauss = new GaussTria(2);3478 for(int ig=gauss->begin();ig<gauss->end();ig++){3479 3480 gauss->GaussPoint(ig);3481 3482 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);3483 GetBSSA(&B[0], &xyz_list[0][0], gauss);3484 GetBprimeSSA(&Bprime[0], &xyz_list[0][0], gauss);3485 3486 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);3487 this->StrainRateSSA(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);3488 material->GetViscosity2d(&viscosity, &epsilon[0]);3489 material->GetViscosity2d(&oldviscosity, &oldepsilon[0]);3490 thickness_input->GetInputValue(&thickness, gauss);3491 3492 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);3493 D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet;3494 for(int i=0;i<3;i++) D[i*3+i]=D_scalar;3495 3496 TripleMultiply(B,3,numdof,1,3497 D,3,3,0,3498 Bprime,3,numdof,0,3499 &Ke->values[0],1);3500 }3501 3502 /*Transform Coordinate System*/3503 ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);3504 3505 /*Clean up and return*/3506 delete gauss;3507 xDelete<IssmDouble>(D);3508 xDelete<IssmDouble>(Bprime);3509 xDelete<IssmDouble>(B);3510 return Ke;3511 }3512 /*}}}*/3513 /*FUNCTION Tria::CreateKMatrixStressbalanceSSAFriction {{{*/3514 ElementMatrix* Tria::CreateKMatrixStressbalanceSSAFriction(void){3515 3516 /*Intermediaries*/3517 bool mainlyfloating;3518 int analysis_type,migration_style;3519 int point1;3520 IssmDouble alpha2,Jdet;3521 IssmDouble fraction1,fraction2;3522 IssmDouble phi=1.0;3523 IssmDouble D_scalar;3524 IssmDouble gllevelset;3525 IssmDouble xyz_list[NUMVERTICES][3];3526 Friction *friction = NULL;3527 GaussTria *gauss = NULL;3528 3529 /*Return if element is inactive*/3530 if(IsFloating()) return NULL;3531 3532 /*Fetch number of nodes and dof for this finite element*/3533 int numnodes = this->NumberofNodes();3534 int numdof = numnodes*NDOF2;3535 3536 /*Initialize Element matrix and vectors*/3537 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,SSAApproximationEnum);3538 IssmDouble* B = xNew<IssmDouble>(2*numdof);3539 IssmDouble* D = xNewZeroInit<IssmDouble>(2*2);3540 3541 /*Retrieve all inputs and parameters*/3542 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);3543 parameters->FindParam(&migration_style,GroundinglineMigrationEnum);3544 Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);3545 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);3546 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);3547 Input* gllevelset_input=NULL;3548 parameters->FindParam(&analysis_type,AnalysisTypeEnum);3549 3550 /*build friction object, used later on: */3551 friction=new Friction(this,2);3552 3553 /*Recover portion of element that is grounded*/3554 if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list[0][0]);3555 if(migration_style==SubelementMigration2Enum){3556 gllevelset_input=inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);3557 this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);3558 gauss=new GaussTria(point1,fraction1,fraction2,mainlyfloating,2);3559 }3560 else{3561 gauss=new GaussTria(2);3562 }3563 3564 /* Start looping on the number of gaussian points: */3565 for(int ig=gauss->begin();ig<gauss->end();ig++){3566 3567 gauss->GaussPoint(ig);3568 3569 friction->GetAlpha2(&alpha2, gauss,vx_input,vy_input,NULL);3570 if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;3571 if(migration_style==SubelementMigration2Enum){3572 gllevelset_input->GetInputValue(&gllevelset, gauss);3573 if(gllevelset<0) alpha2=0;3574 }3575 3576 GetBSSAFriction(&B[0], &xyz_list[0][0], gauss);3577 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);3578 D_scalar=alpha2*gauss->weight*Jdet;3579 for(int i=0;i<2;i++) D[i*2+i]=D_scalar;3580 3581 TripleMultiply(B,2,numdof,1,3582 D,2,2,0,3583 B,2,numdof,0,3584 &Ke->values[0],1);3585 }3586 3587 /*Transform Coordinate System*/3588 ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);3589 3590 /*Clean up and return*/3591 delete gauss;3592 delete friction;3593 xDelete<IssmDouble>(D);3594 xDelete<IssmDouble>(B);3595 return Ke;3596 }3597 /*}}}*/3598 /*FUNCTION Tria::CreateKMatrixStressbalanceSIA{{{*/3599 ElementMatrix* Tria::CreateKMatrixStressbalanceSIA(void){3600 3601 /*Intermediaries*/3602 IssmDouble connectivity;3603 3604 /*Fetch number of nodes and dof for this finite element*/3605 int numnodes = this->NumberofNodes(); _assert_(numnodes==3);3606 int numdof = numnodes*NDOF2;3607 3608 /*Initialize Element matrix*/3609 ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);3610 3611 /*Create Element matrix*/3612 for(int i=0;i<numnodes;i++){3613 connectivity=(IssmDouble)vertices[i]->Connectivity();3614 Ke->values[(2*i+0)*numdof+(2*i+0)]=1./connectivity;3615 Ke->values[(2*i+1)*numdof+(2*i+1)]=1./connectivity;3616 }3617 3618 /*Clean up and return*/3619 return Ke;3620 }3621 /*}}}*/3622 /*FUNCTION Tria::CreateJacobianStressbalanceSSA{{{*/3623 ElementMatrix* Tria::CreateJacobianStressbalanceSSA(void){3624 3625 /*Intermediaries */3626 int i,j;3627 IssmDouble xyz_list[NUMVERTICES][3];3628 IssmDouble Jdet,thickness;3629 IssmDouble eps1dotdphii,eps1dotdphij;3630 IssmDouble eps2dotdphii,eps2dotdphij;3631 IssmDouble mu_prime;3632 IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/3633 IssmDouble eps1[2],eps2[2];3634 GaussTria* gauss = NULL;3635 3636 /*Fetch number of nodes and dof for this finite element*/3637 int numnodes = this->NumberofNodes();3638 3639 /*Initialize Element matrix, vectors and Gaussian points*/3640 ElementMatrix* Ke=CreateKMatrixStressbalanceSSA(); //Initialize Jacobian with regular SSA (first part of the Gateau derivative)3641 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);3642 3643 /*Retrieve all inputs and parameters*/3644 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);3645 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);3646 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);3647 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);3648 3649 /* Start looping on the number of gaussian points: */3650 gauss=new GaussTria(2);3651 for(int ig=gauss->begin();ig<gauss->end();ig++){3652 3653 gauss->GaussPoint(ig);3654 3655 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);3656 GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);3657 3658 thickness_input->GetInputValue(&thickness, gauss);3659 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);3660 material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);3661 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];3662 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];3663 3664 for(i=0;i<numnodes;i++){3665 for(j=0;j<numnodes;j++){3666 eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];3667 eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];3668 eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];3669 eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];3670 3671 Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;3672 Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;3673 Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;3674 Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;3675 }3676 }3677 }3678 3679 /*Transform Coordinate System*/3680 ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);3681 3682 /*Clean up and return*/3683 xDelete<IssmDouble>(dbasis);3684 delete gauss;3685 return Ke;3686 }3687 /*}}}*/3688 3065 /*FUNCTION Tria::GetYcoord {{{*/ 3689 3066 IssmDouble Tria::GetYcoord(GaussTria* gauss){ … … 5037 4414 } 5038 4415 /*}}}*/ 5039 /*FUNCTION Tria::CreateKMatrixAdjointBalancethickness {{{*/5040 ElementMatrix* Tria::CreateKMatrixAdjointBalancethickness(void){5041 5042 ElementMatrix* Ke=NULL;5043 5044 /*Get Element Matrix of the forward model*/5045 switch(GetElementType()){5046 case P1Enum:5047 Ke=CreateKMatrixBalancethickness_CG();5048 break;5049 case P1DGEnum:5050 Ke=CreateKMatrixBalancethickness_DG();5051 break;5052 default:5053 _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");5054 }5055 5056 /*Transpose and return Ke*/5057 Ke->Transpose();5058 return Ke;5059 }5060 /*}}}*/5061 /*FUNCTION Tria::CreateKMatrixAdjointSSA{{{*/5062 ElementMatrix* Tria::CreateKMatrixAdjointSSA(void){5063 5064 /*Intermediaries */5065 int i,j;5066 bool incomplete_adjoint;5067 IssmDouble xyz_list[NUMVERTICES][3];5068 IssmDouble Jdet,thickness;5069 IssmDouble eps1dotdphii,eps1dotdphij;5070 IssmDouble eps2dotdphii,eps2dotdphij;5071 IssmDouble mu_prime;5072 IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/5073 IssmDouble eps1[2],eps2[2];5074 GaussTria *gauss=NULL;5075 5076 /*Fetch number of nodes and dof for this finite element*/5077 int numnodes = this->NumberofNodes();5078 5079 /*Initialize Jacobian with regular SSA (first part of the Gateau derivative)*/5080 parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);5081 ElementMatrix* Ke=CreateKMatrixStressbalanceSSA();5082 if(incomplete_adjoint) return Ke;5083 5084 /*Retrieve all inputs and parameters*/5085 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);5086 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);5087 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);5088 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);5089 5090 /*Allocate dbasis*/5091 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);5092 5093 /* Start looping on the number of gaussian points: */5094 gauss=new GaussTria(2);5095 for(int ig=gauss->begin();ig<gauss->end();ig++){5096 5097 gauss->GaussPoint(ig);5098 5099 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);5100 GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);5101 5102 thickness_input->GetInputValue(&thickness, gauss);5103 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);5104 material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);5105 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];5106 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];5107 5108 for(i=0;i<numnodes;i++){5109 for(j=0;j<numnodes;j++){5110 eps1dotdphii=eps1[0]*dbasis[0*numnodes+i]+eps1[1]*dbasis[1*numnodes+i];5111 eps1dotdphij=eps1[0]*dbasis[0*numnodes+j]+eps1[1]*dbasis[1*numnodes+j];5112 eps2dotdphii=eps2[0]*dbasis[0*numnodes+i]+eps2[1]*dbasis[1*numnodes+i];5113 eps2dotdphij=eps2[0]*dbasis[0*numnodes+j]+eps2[1]*dbasis[1*numnodes+j];5114 5115 Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;5116 Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;5117 Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;5118 Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;5119 }5120 }5121 }5122 5123 /*Transform Coordinate System*/5124 ::TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum);5125 5126 /*Clean up and return*/5127 delete gauss;5128 xDelete<IssmDouble>(dbasis);5129 //Ke->Transpose();5130 return Ke;5131 }5132 /*}}}*/5133 4416 /*FUNCTION Tria::GetVectorFromControlInputs{{{*/ 5134 4417 void Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){ … … 5191 4474 5192 4475 ((ControlInput*)input)->SetInput(new_input); 5193 }5194 /*}}}*/5195 #endif5196 5197 #ifdef _HAVE_THERMAL_5198 /*FUNCTION Tria::CreateKMatrixMelting {{{*/5199 ElementMatrix* Tria::CreateKMatrixMelting(void){5200 5201 /*Intermediaries */5202 IssmDouble heatcapacity,latentheat;5203 IssmDouble Jdet,D_scalar;5204 IssmDouble xyz_list[NUMVERTICES][3];5205 5206 /*Fetch number of nodes and dof for this finite element*/5207 int numnodes = this->NumberofNodes();5208 5209 /*Initialize Element vector and vectors*/5210 ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);5211 IssmDouble* basis = xNew<IssmDouble>(numnodes);5212 5213 /*Retrieve all inputs and parameters*/5214 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);5215 latentheat = matpar->GetLatentHeat();5216 heatcapacity = matpar->GetHeatCapacity();5217 5218 /* Start looping on the number of gauss (nodes on the bedrock) */5219 GaussTria* gauss=new GaussTria(2);5220 for(int ig=gauss->begin();ig<gauss->end();ig++){5221 5222 gauss->GaussPoint(ig);5223 5224 GetNodalFunctions(&basis[0],gauss);5225 GetJacobianDeterminant(&Jdet, &xyz_list[0][0], gauss);5226 5227 D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;5228 5229 TripleMultiply(basis,numnodes,1,0,5230 &D_scalar,1,1,0,5231 basis,1,numnodes,0,5232 &Ke->values[0],1);5233 }5234 5235 /*Clean up and return*/5236 xDelete<IssmDouble>(basis);5237 delete gauss;5238 return Ke;5239 4476 } 5240 4477 /*}}}*/ … … 5294 4531 this->inputs->AddInput(new TriaInput(HydrologyWaterVxEnum,vx,P1Enum)); 5295 4532 this->inputs->AddInput(new TriaInput(HydrologyWaterVyEnum,vy,P1Enum)); 5296 }5297 /*}}}*/5298 /*FUNCTION Tria::CreateKMatrixHydrologyShreve{{{*/5299 ElementMatrix* Tria::CreateKMatrixHydrologyShreve(void){5300 5301 /*Intermediaries */5302 IssmDouble diffusivity;5303 IssmDouble Jdet,D_scalar,dt,h;5304 IssmDouble vx,vy,vel,dvxdx,dvydy;5305 IssmDouble dvx[2],dvy[2];5306 IssmDouble xyz_list[NUMVERTICES][3];5307 5308 /*Skip if water or ice shelf element*/5309 if(NoIceInElement() || IsFloating()) return NULL;5310 5311 /*Fetch number of nodes and dof for this finite element*/5312 int numnodes = this->NumberofNodes();5313 5314 /*Initialize Element matrix and vectors*/5315 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);5316 IssmDouble* basis = xNew<IssmDouble>(numnodes);5317 IssmDouble* B = xNew<IssmDouble>(2*numnodes);5318 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);5319 IssmDouble D[2][2];5320 5321 /*Create water velocity vx and vy from current inputs*/5322 CreateHydrologyWaterVelocityInput();5323 5324 /*Retrieve all inputs and parameters*/5325 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);5326 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);5327 this->parameters->FindParam(&diffusivity,HydrologyshreveStabilizationEnum);5328 Input* vx_input=inputs->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);5329 Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);5330 h=sqrt(2*this->GetArea());5331 5332 /* Start looping on the number of gaussian points: */5333 GaussTria* gauss=new GaussTria(2);5334 for(int ig=gauss->begin();ig<gauss->end();ig++){5335 5336 gauss->GaussPoint(ig);5337 5338 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);5339 GetNodalFunctions(basis,gauss);5340 5341 vx_input->GetInputValue(&vx,gauss);5342 vy_input->GetInputValue(&vy,gauss);5343 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);5344 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);5345 5346 D_scalar=gauss->weight*Jdet;5347 5348 TripleMultiply(basis,1,numnodes,1,5349 &D_scalar,1,1,0,5350 basis,1,numnodes,0,5351 Ke->values,1);5352 5353 GetBMasstransport(B,&xyz_list[0][0], gauss);5354 GetBprimeMasstransport(Bprime,&xyz_list[0][0], gauss);5355 5356 dvxdx=dvx[0];5357 dvydy=dvy[1];5358 D_scalar=dt*gauss->weight*Jdet;5359 5360 D[0][0]=D_scalar*dvxdx;5361 D[0][1]=0.;5362 D[1][0]=0.;5363 D[1][1]=D_scalar*dvydy;5364 TripleMultiply(B,2,numnodes,1,5365 &D[0][0],2,2,0,5366 B,2,numnodes,0,5367 &Ke->values[0],1);5368 5369 D[0][0]=D_scalar*vx;5370 D[1][1]=D_scalar*vy;5371 TripleMultiply(B,2,numnodes,1,5372 &D[0][0],2,2,0,5373 Bprime,2,numnodes,0,5374 &Ke->values[0],1);5375 5376 /*Artificial diffusivity*/5377 vel=sqrt(vx*vx+vy*vy);5378 D[0][0]=D_scalar*diffusivity*h/(2*vel)*vx*vx;5379 D[1][0]=D_scalar*diffusivity*h/(2*vel)*vy*vx;5380 D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;5381 D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;5382 TripleMultiply(Bprime,2,numnodes,1,5383 &D[0][0],2,2,0,5384 Bprime,2,numnodes,0,5385 &Ke->values[0],1);5386 }5387 5388 /*Clean up and return*/5389 xDelete<IssmDouble>(basis);5390 xDelete<IssmDouble>(B);5391 xDelete<IssmDouble>(Bprime);5392 delete gauss;5393 return Ke;5394 }5395 /*}}}*/5396 /*FUNCTION Tria::CreateKMatrixHydrologyDCInefficient{{{*/5397 ElementMatrix* Tria::CreateKMatrixHydrologyDCInefficient(void){5398 5399 /* Intermediaries */5400 IssmDouble D_scalar,Jdet;5401 IssmDouble sediment_transmitivity,dt;5402 IssmDouble sediment_storing;5403 IssmDouble xyz_list[NUMVERTICES][3];5404 5405 /*Fetch number of nodes and dof for this finite element*/5406 int numnodes = this->NumberofNodes();5407 5408 /*Initialize Element matrix and vectors*/5409 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters);5410 IssmDouble* B = xNew<IssmDouble>(5*numnodes);5411 IssmDouble* basis = xNew<IssmDouble>(numnodes);5412 IssmDouble D[2][2];5413 5414 /*Retrieve all inputs and parameters*/5415 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);5416 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);5417 sediment_storing = matpar->GetSedimentStoring();5418 sediment_transmitivity = matpar->GetSedimentTransmitivity();5419 5420 /* Start looping on the number of gaussian points: */5421 GaussTria* gauss=new GaussTria(2);5422 for(int ig=gauss->begin();ig<gauss->end();ig++){5423 5424 gauss->GaussPoint(ig);5425 5426 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);5427 5428 /*Diffusivity*/5429 D_scalar=sediment_transmitivity*gauss->weight*Jdet;5430 if(reCast<bool,IssmDouble>(dt)) D_scalar=D_scalar*dt;5431 D[0][0]=D_scalar; D[0][1]=0.;5432 D[1][0]=0.; D[1][1]=D_scalar;5433 GetBHydro(B,&xyz_list[0][0],gauss);5434 TripleMultiply(B,2,numnodes,1,5435 &D[0][0],2,2,0,5436 B,2,numnodes,0,5437 &Ke->values[0],1);5438 5439 /*Transient*/5440 if(reCast<bool,IssmDouble>(dt)){5441 GetNodalFunctions(&basis[0],gauss);5442 D_scalar=sediment_storing*gauss->weight*Jdet;5443 5444 TripleMultiply(basis,numnodes,1,0,5445 &D_scalar,1,1,0,5446 basis,1,numnodes,0,5447 &Ke->values[0],1);5448 }5449 }5450 /*Clean up and return*/5451 xDelete<IssmDouble>(basis);5452 xDelete<IssmDouble>(B);5453 delete gauss;5454 return Ke;5455 }5456 /*}}}*/5457 /*FUNCTION Tria::CreateKMatrixHydrologyDCEfficient{{{*/5458 ElementMatrix* Tria::CreateKMatrixHydrologyDCEfficient(void){5459 5460 /* Intermediaries */5461 IssmDouble D_scalar,Jdet,dt;5462 IssmDouble epl_thickness;5463 IssmDouble epl_conductivity;5464 IssmDouble epl_specificstoring;5465 IssmDouble xyz_list[NUMVERTICES][3];5466 5467 /*Check that all nodes are active, else return empty matrix*/5468 if(!this->AllActive()){5469 return NULL;5470 }5471 5472 /*Fetch number of nodes and dof for this finite element*/5473 int numnodes = this->NumberofNodes();5474 5475 /*Initialize Element matrix and vectors*/5476 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters);5477 IssmDouble* B = xNew<IssmDouble>(5*numnodes);5478 IssmDouble* basis = xNew<IssmDouble>(numnodes);5479 IssmDouble D[2][2];5480 5481 /*Retrieve all inputs and parameters*/5482 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);5483 Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input);5484 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);5485 epl_specificstoring = matpar->GetEplSpecificStoring();5486 epl_conductivity = matpar->GetEplConductivity();5487 5488 5489 /* Start looping on the number of gaussian points: */5490 GaussTria* gauss=new GaussTria(2);5491 for(int ig=gauss->begin();ig<gauss->end();ig++){5492 5493 5494 gauss->GaussPoint(ig);5495 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);5496 thickness_input->GetInputValue(&epl_thickness,gauss);5497 5498 /*Diffusivity*/5499 D_scalar=epl_conductivity*epl_thickness*gauss->weight*Jdet;5500 if(reCast<bool,IssmDouble>(dt)) D_scalar=D_scalar*dt;5501 D[0][0]=D_scalar; D[0][1]=0.;5502 D[1][0]=0.; D[1][1]=D_scalar;5503 GetBHydro(B,&xyz_list[0][0],gauss);5504 TripleMultiply(B,2,numnodes,1,5505 &D[0][0],2,2,0,5506 B,2,numnodes,0,5507 &Ke->values[0],1);5508 5509 /*Transient*/5510 if(reCast<bool,IssmDouble>(dt)){5511 GetNodalFunctions(basis,gauss);5512 D_scalar=epl_specificstoring*epl_thickness*gauss->weight*Jdet;5513 5514 TripleMultiply(basis,numnodes,1,0,5515 &D_scalar,1,1,0,5516 basis,1,numnodes,0,5517 &Ke->values[0],1);5518 }5519 }5520 5521 /*Clean up and return*/5522 xDelete<IssmDouble>(basis);5523 xDelete<IssmDouble>(B);5524 delete gauss;5525 return Ke;5526 4533 } 5527 4534 /*}}}*/ … … 5898 4905 #endif 5899 4906 5900 #ifdef _HAVE_MASSTRANSPORT_5901 /*FUNCTION Tria::CreateKMatrixExtrusion {{{*/5902 ElementMatrix* Tria::CreateKMatrixExtrusion(void){5903 5904 /*compute all stiffness matrices for this element*/5905 ElementMatrix* Ke1=CreateKMatrixExtrusionVolume();5906 ElementMatrix* Ke2=CreateKMatrixExtrusionSurface();5907 ElementMatrix* Ke3=CreateKMatrixExtrusionBed();5908 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);5909 5910 /*clean-up and return*/5911 delete Ke1;5912 delete Ke2;5913 delete Ke3;5914 return Ke;5915 5916 }5917 /*}}}*/5918 /*FUNCTION Tria::CreateKMatrixExtrusionVolume {{{*/5919 ElementMatrix* Tria::CreateKMatrixExtrusionVolume(void){5920 5921 /*Intermediaries */5922 IssmDouble Jdet;5923 IssmDouble xyz_list[NUMVERTICES][3];5924 IssmDouble B[NDOF1][NUMVERTICES];5925 IssmDouble Bprime[NDOF1][NUMVERTICES];5926 IssmDouble DL_scalar;5927 GaussTria *gauss=NULL;5928 5929 /*Initialize Element matrix*/5930 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);5931 5932 /*Retrieve all inputs and parameters*/5933 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);5934 5935 /* Start looping on the number of gaussian points: */5936 gauss=new GaussTria(2);5937 for(int ig=gauss->begin();ig<gauss->end();ig++){5938 5939 gauss->GaussPoint(ig);5940 5941 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);5942 GetBExtrusion(&B[0][0], &xyz_list[0][0], gauss);5943 GetNodalFunctions(&Bprime[0][0],gauss);5944 5945 DL_scalar=gauss->weight*Jdet;5946 5947 TripleMultiply(&B[0][0],1,NUMVERTICES,1,5948 &DL_scalar,1,1,0,5949 &Bprime[0][0],1,NUMVERTICES,0,5950 &Ke->values[0],1);5951 }5952 5953 /*Clean up and return*/5954 delete gauss;5955 return Ke;5956 }5957 /*}}}*/5958 /*FUNCTION Tria::CreateKMatrixExtrusionSurface {{{*/5959 ElementMatrix* Tria::CreateKMatrixExtrusionSurface(void){5960 5961 if (!HasEdgeOnSurface()) return NULL;5962 5963 /*Constants*/5964 const int numdof=NDOF1*NUMVERTICES;5965 5966 /*Intermediaries */5967 int indices[2];5968 IssmDouble xyz_list[NUMVERTICES][3];5969 IssmDouble xyz_list_seg[NUMVERTICES1D][3];5970 IssmDouble normal[3];5971 IssmDouble Jdet,DL_scalar;5972 IssmDouble basis[NUMVERTICES];5973 GaussTria *gauss=NULL;5974 5975 /*Initialize Element matrix*/5976 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);5977 5978 /*Retrieve all inputs and parameters*/5979 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);5980 5981 /*Get vertex indices that lie on bed*/5982 this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);5983 for(int i=0;i<NUMVERTICES1D;i++) for(int j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];5984 NormalSection(&normal[0],&xyz_list_seg[0][0]);5985 5986 /* Start looping on the number of gaussian points: */5987 gauss=new GaussTria(indices[0],indices[1],2);5988 for(int ig=gauss->begin();ig<gauss->end();ig++){5989 5990 gauss->GaussPoint(ig);5991 5992 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);5993 GetNodalFunctions(&basis[0], gauss);5994 5995 DL_scalar= - gauss->weight*Jdet*normal[1];5996 5997 TripleMultiply( basis,1,numdof,1,5998 &DL_scalar,1,1,0,5999 basis,1,numdof,0,6000 &Ke->values[0],1);6001 }6002 6003 /*Clean up and return*/6004 delete gauss;6005 return Ke;6006 }6007 /*}}}*/6008 /*FUNCTION Tria::CreateKMatrixExtrusionBed {{{*/6009 ElementMatrix* Tria::CreateKMatrixExtrusionBed(void){6010 6011 if (!HasEdgeOnBed()) return NULL;6012 6013 /*Constants*/6014 const int numdof=NDOF1*NUMVERTICES;6015 6016 /*Intermediaries */6017 int indices[2];6018 IssmDouble xyz_list[NUMVERTICES][3];6019 IssmDouble xyz_list_seg[NUMVERTICES1D][3];6020 IssmDouble normal[3];6021 IssmDouble Jdet,DL_scalar;6022 IssmDouble basis[NUMVERTICES];6023 GaussTria *gauss=NULL;6024 6025 /*Initialize Element matrix*/6026 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);6027 6028 /*Retrieve all inputs and parameters*/6029 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);6030 6031 /*Get vertex indices that lie on bed*/6032 this->EdgeOnBedIndices(&indices[0],&indices[1]);6033 for(int i=0;i<NUMVERTICES1D;i++) for(int j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];6034 NormalSection(&normal[0],&xyz_list_seg[0][0]);6035 6036 /* Start looping on the number of gaussian points: */6037 gauss=new GaussTria(indices[0],indices[1],2);6038 for(int ig=gauss->begin();ig<gauss->end();ig++){6039 6040 gauss->GaussPoint(ig);6041 6042 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);6043 GetNodalFunctions(&basis[0], gauss);6044 6045 DL_scalar= - gauss->weight*Jdet*normal[1];6046 6047 TripleMultiply( basis,1,numdof,1,6048 &DL_scalar,1,1,0,6049 basis,1,numdof,0,6050 &Ke->values[0],1);6051 }6052 6053 /*Clean up and return*/6054 delete gauss;6055 return Ke;6056 }6057 /*}}}*/6058 /*FUNCTION Tria::CreateKMatrixMasstransport {{{*/6059 ElementMatrix* Tria::CreateKMatrixMasstransport(void){6060 6061 switch(GetElementType()){6062 case P1Enum: case P2Enum:6063 return CreateKMatrixMasstransport_CG();6064 case P1DGEnum:6065 return CreateKMatrixMasstransport_DG();6066 default:6067 _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");6068 }6069 }6070 /*}}}*/6071 /*FUNCTION Tria::CreateKMatrixMasstransport_CG {{{*/6072 ElementMatrix* Tria::CreateKMatrixMasstransport_CG(void){6073 6074 /*Intermediaries */6075 int stabilization;6076 int meshtype;6077 IssmDouble Jdet,D_scalar,dt,h;6078 IssmDouble vel,vx,vy,dvxdx,dvydy;6079 IssmDouble dvx[2],dvy[2];6080 IssmDouble xyz_list[NUMVERTICES][3];6081 6082 /*Fetch number of nodes for this finite element*/6083 int numnodes = this->NumberofNodes();6084 6085 /*Initialize Element matrix and vectors*/6086 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);6087 IssmDouble* basis = xNew<IssmDouble>(numnodes);6088 IssmDouble* B = xNew<IssmDouble>(2*numnodes);6089 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);6090 IssmDouble D[2][2];6091 6092 /*Retrieve all inputs and parameters*/6093 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);6094 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);6095 this->parameters->FindParam(&meshtype,MeshTypeEnum);6096 this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);6097 Input* vxaverage_input=NULL;6098 Input* vyaverage_input=NULL;6099 if(meshtype==Mesh2DhorizontalEnum){6100 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);6101 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);6102 }6103 else{6104 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);6105 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);6106 }6107 h=sqrt(2*this->GetArea());6108 6109 /* Start looping on the number of gaussian points: */6110 GaussTria *gauss=new GaussTria(2);6111 for(int ig=gauss->begin();ig<gauss->end();ig++){6112 6113 gauss->GaussPoint(ig);6114 6115 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);6116 GetNodalFunctions(basis,gauss);6117 6118 vxaverage_input->GetInputValue(&vx,gauss);6119 vyaverage_input->GetInputValue(&vy,gauss);6120 vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);6121 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);6122 6123 D_scalar=gauss->weight*Jdet;6124 6125 TripleMultiply(basis,1,numnodes,1,6126 &D_scalar,1,1,0,6127 basis,1,numnodes,0,6128 &Ke->values[0],1);6129 6130 GetBMasstransport(B,&xyz_list[0][0],gauss);6131 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);6132 6133 dvxdx=dvx[0];6134 dvydy=dvy[1];6135 D_scalar=dt*gauss->weight*Jdet;6136 6137 D[0][0]=D_scalar*dvxdx;6138 D[0][1]=0.;6139 D[1][0]=0.;6140 D[1][1]=D_scalar*dvydy;6141 TripleMultiply(B,2,numnodes,1,6142 &D[0][0],2,2,0,6143 B,2,numnodes,0,6144 &Ke->values[0],1);6145 6146 D[0][0]=D_scalar*vx;6147 D[1][1]=D_scalar*vy;6148 TripleMultiply(B,2,numnodes,1,6149 &D[0][0],2,2,0,6150 Bprime,2,numnodes,0,6151 &Ke->values[0],1);6152 6153 if(stabilization==2){6154 /*Streamline upwinding*/6155 vel=sqrt(vx*vx+vy*vy)+1.e-8;6156 D[0][0]=h/(2*vel)*vx*vx;6157 D[1][0]=h/(2*vel)*vy*vx;6158 D[0][1]=h/(2*vel)*vx*vy;6159 D[1][1]=h/(2*vel)*vy*vy;6160 }6161 else if(stabilization==1){6162 /*SSA*/6163 vxaverage_input->GetInputAverage(&vx);6164 vyaverage_input->GetInputAverage(&vy);6165 D[0][0]=h/2.0*fabs(vx);6166 D[0][1]=0.;6167 D[1][0]=0.;6168 D[1][1]=h/2.0*fabs(vy);6169 }6170 if(stabilization==1 || stabilization==2){6171 D[0][0]=D_scalar*D[0][0];6172 D[1][0]=D_scalar*D[1][0];6173 D[0][1]=D_scalar*D[0][1];6174 D[1][1]=D_scalar*D[1][1];6175 TripleMultiply(Bprime,2,numnodes,1,6176 &D[0][0],2,2,0,6177 Bprime,2,numnodes,0,6178 &Ke->values[0],1);6179 }6180 }6181 6182 /*Clean up and return*/6183 xDelete<IssmDouble>(basis);6184 xDelete<IssmDouble>(B);6185 xDelete<IssmDouble>(Bprime);6186 delete gauss;6187 return Ke;6188 }6189 /*}}}*/6190 /*FUNCTION Tria::CreateKMatrixMasstransport_DG {{{*/6191 ElementMatrix* Tria::CreateKMatrixMasstransport_DG(void){6192 6193 /*Intermediaries */6194 int meshtype;6195 IssmDouble xyz_list[NUMVERTICES][3];6196 IssmDouble Jdet,D_scalar,dt,vx,vy;6197 6198 /*Fetch number of nodes for this finite element*/6199 int numnodes = this->NumberofNodes();6200 6201 /*Initialize Element matrix and vectors*/6202 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);6203 IssmDouble* basis = xNew<IssmDouble>(numnodes);6204 IssmDouble* B = xNew<IssmDouble>(2*numnodes);6205 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);6206 IssmDouble D[2][2];6207 6208 /*Retrieve all inputs and parameters*/6209 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);6210 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);6211 this->parameters->FindParam(&meshtype,MeshTypeEnum);6212 Input* vxaverage_input=NULL;6213 Input* vyaverage_input=NULL;6214 if(meshtype==Mesh2DhorizontalEnum){6215 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);6216 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);6217 }6218 else{6219 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);6220 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);6221 }6222 6223 /* Start looping on the number of gaussian points: */6224 GaussTria* gauss=new GaussTria(2);6225 for(int ig=gauss->begin();ig<gauss->end();ig++){6226 6227 gauss->GaussPoint(ig);6228 6229 vxaverage_input->GetInputValue(&vx,gauss);6230 vyaverage_input->GetInputValue(&vy,gauss);6231 6232 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);6233 GetNodalFunctions(basis,gauss);6234 6235 D_scalar=gauss->weight*Jdet;6236 6237 TripleMultiply(basis,1,numnodes,1,6238 &D_scalar,1,1,0,6239 basis,1,numnodes,0,6240 &Ke->values[0],1);6241 6242 /*WARNING: B and Bprime are inverted compared to usual masstransport!!!!*/6243 GetBMasstransport(Bprime, &xyz_list[0][0], gauss);6244 GetBprimeMasstransport(B, &xyz_list[0][0], gauss);6245 6246 D_scalar=-dt*gauss->weight*Jdet;6247 D[0][0]=D_scalar*vx;6248 D[0][1]=0.;6249 D[1][0]=0.;6250 D[1][1]=D_scalar*vy;6251 6252 TripleMultiply(B,2,numnodes,1,6253 &D[0][0],2,2,0,6254 Bprime,2,numnodes,0,6255 &Ke->values[0],1);6256 }6257 6258 /*Clean up and return*/6259 xDelete<IssmDouble>(basis);6260 xDelete<IssmDouble>(B);6261 xDelete<IssmDouble>(Bprime);6262 delete gauss;6263 return Ke;6264 }6265 /*}}}*/6266 /*FUNCTION Tria::CreateKMatrixFreeSurfaceTop {{{*/6267 ElementMatrix* Tria::CreateKMatrixFreeSurfaceTop(void){6268 6269 /*Intermediaries */6270 int stabilization;6271 IssmDouble Jdet,D_scalar,dt,h;6272 IssmDouble vel,vx,vy;6273 IssmDouble xyz_list[NUMVERTICES][3];6274 6275 /*Fetch number of nodes for this finite element*/6276 int numnodes = this->NumberofNodes();6277 6278 /*Initialize Element matrix and vectors*/6279 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);6280 IssmDouble* basis = xNew<IssmDouble>(numnodes);6281 IssmDouble* B = xNew<IssmDouble>(2*numnodes);6282 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);6283 IssmDouble D[2][2];6284 6285 /*Retrieve all inputs and parameters*/6286 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);6287 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);6288 this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);6289 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);6290 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);6291 h=sqrt(2*this->GetArea());6292 6293 /* Start looping on the number of gaussian points: */6294 GaussTria *gauss=new GaussTria(2);6295 for(int ig=gauss->begin();ig<gauss->end();ig++){6296 6297 gauss->GaussPoint(ig);6298 6299 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);6300 GetNodalFunctions(basis,gauss);6301 6302 vx_input->GetInputValue(&vx,gauss);6303 vy_input->GetInputValue(&vy,gauss);6304 6305 D_scalar=gauss->weight*Jdet;6306 6307 TripleMultiply(basis,1,numnodes,1,6308 &D_scalar,1,1,0,6309 basis,1,numnodes,0,6310 &Ke->values[0],1);6311 6312 GetBMasstransport(B,&xyz_list[0][0],gauss);6313 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);6314 6315 D_scalar=dt*gauss->weight*Jdet;6316 6317 D[0][0]=D_scalar*vx;6318 D[0][1]=0.;6319 D[1][0]=0.;6320 D[1][1]=D_scalar*vy;6321 TripleMultiply(B,2,numnodes,1,6322 &D[0][0],2,2,0,6323 Bprime,2,numnodes,0,6324 &Ke->values[0],1);6325 6326 if(stabilization==2){6327 /*Streamline upwinding*/6328 vel=sqrt(vx*vx+vy*vy)+1.e-8;6329 D[0][0]=h/(2*vel)*vx*vx;6330 D[1][0]=h/(2*vel)*vy*vx;6331 D[0][1]=h/(2*vel)*vx*vy;6332 D[1][1]=h/(2*vel)*vy*vy;6333 }6334 else if(stabilization==1){6335 /*SSA*/6336 vx_input->GetInputAverage(&vx);6337 vy_input->GetInputAverage(&vy);6338 D[0][0]=h/2.0*fabs(vx);6339 D[0][1]=0.;6340 D[1][0]=0.;6341 D[1][1]=h/2.0*fabs(vy);6342 }6343 if(stabilization==1 || stabilization==2){6344 D[0][0]=D_scalar*D[0][0];6345 D[1][0]=D_scalar*D[1][0];6346 D[0][1]=D_scalar*D[0][1];6347 D[1][1]=D_scalar*D[1][1];6348 TripleMultiply(Bprime,2,numnodes,1,6349 &D[0][0],2,2,0,6350 Bprime,2,numnodes,0,6351 &Ke->values[0],1);6352 }6353 }6354 6355 /*Clean up and return*/6356 xDelete<IssmDouble>(basis);6357 xDelete<IssmDouble>(B);6358 xDelete<IssmDouble>(Bprime);6359 delete gauss;6360 return Ke;6361 }6362 /*}}}*/6363 /*FUNCTION Tria::CreateKMatrixFreeSurfaceTop1D {{{*/6364 ElementMatrix* Tria::CreateKMatrixFreeSurfaceTop1D(void){6365 6366 if(!HasEdgeOnSurface()) return NULL;6367 6368 int index1,index2;6369 this->EdgeOnSurfaceIndices(&index1,&index2);6370 6371 Seg* seg=(Seg*)SpawnSeg(index1,index2);6372 ElementMatrix* Ke=seg->CreateKMatrixFreeSurfaceTop();6373 delete seg->material; delete seg;6374 6375 /*clean up and return*/6376 return Ke;6377 }6378 /*}}}*/6379 /*FUNCTION Tria::CreateKMatrixFreeSurfaceBase {{{*/6380 ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase(void){6381 6382 /*Intermediaries */6383 int stabilization;6384 IssmDouble Jdet,D_scalar,dt,h;6385 IssmDouble vel,vx,vy;6386 IssmDouble xyz_list[NUMVERTICES][3];6387 6388 /*Fetch number of nodes for this finite element*/6389 int numnodes = this->NumberofNodes();6390 6391 /*Initialize Element matrix and vectors*/6392 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);6393 IssmDouble* basis = xNew<IssmDouble>(numnodes);6394 IssmDouble* B = xNew<IssmDouble>(2*numnodes);6395 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);6396 IssmDouble D[2][2];6397 6398 /*Retrieve all inputs and parameters*/6399 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);6400 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);6401 this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);6402 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);6403 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);6404 h=sqrt(2*this->GetArea());6405 6406 /* Start looping on the number of gaussian points: */6407 GaussTria *gauss=new GaussTria(2);6408 for(int ig=gauss->begin();ig<gauss->end();ig++){6409 6410 gauss->GaussPoint(ig);6411 6412 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);6413 GetNodalFunctions(basis,gauss);6414 6415 vx_input->GetInputValue(&vx,gauss);6416 vy_input->GetInputValue(&vy,gauss);6417 6418 D_scalar=gauss->weight*Jdet;6419 6420 TripleMultiply(basis,1,numnodes,1,6421 &D_scalar,1,1,0,6422 basis,1,numnodes,0,6423 &Ke->values[0],1);6424 6425 GetBMasstransport(B,&xyz_list[0][0],gauss);6426 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);6427 6428 D_scalar=dt*gauss->weight*Jdet;6429 6430 D[0][0]=D_scalar*vx;6431 D[0][1]=0.;6432 D[1][0]=0.;6433 D[1][1]=D_scalar*vy;6434 TripleMultiply(B,2,numnodes,1,6435 &D[0][0],2,2,0,6436 Bprime,2,numnodes,0,6437 &Ke->values[0],1);6438 6439 if(stabilization==2){6440 /*Streamline upwinding*/6441 vel=sqrt(vx*vx+vy*vy)+1.e-8;6442 D[0][0]=h/(2*vel)*vx*vx;6443 D[1][0]=h/(2*vel)*vy*vx;6444 D[0][1]=h/(2*vel)*vx*vy;6445 D[1][1]=h/(2*vel)*vy*vy;6446 }6447 else if(stabilization==1){6448 /*SSA*/6449 vx_input->GetInputAverage(&vx);6450 vy_input->GetInputAverage(&vy);6451 D[0][0]=h/2.0*fabs(vx);6452 D[0][1]=0.;6453 D[1][0]=0.;6454 D[1][1]=h/2.0*fabs(vy);6455 }6456 if(stabilization==1 || stabilization==2){6457 D[0][0]=D_scalar*D[0][0];6458 D[1][0]=D_scalar*D[1][0];6459 D[0][1]=D_scalar*D[0][1];6460 D[1][1]=D_scalar*D[1][1];6461 TripleMultiply(Bprime,2,numnodes,1,6462 &D[0][0],2,2,0,6463 Bprime,2,numnodes,0,6464 &Ke->values[0],1);6465 }6466 }6467 6468 /*Clean up and return*/6469 xDelete<IssmDouble>(basis);6470 xDelete<IssmDouble>(B);6471 xDelete<IssmDouble>(Bprime);6472 delete gauss;6473 return Ke;6474 }6475 /*}}}*/6476 /*FUNCTION Tria::CreateKMatrixFreeSurfaceBase1D {{{*/6477 ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase1D(void){6478 6479 if(!HasEdgeOnBed()) return NULL;6480 6481 int index1,index2;6482 this->EdgeOnBedIndices(&index1,&index2);6483 6484 Seg* seg=(Seg*)SpawnSeg(index1,index2);6485 ElementMatrix* Ke=seg->CreateKMatrixFreeSurfaceBase();6486 delete seg->material; delete seg;6487 6488 /*clean up and return*/6489 return Ke;6490 }6491 /*}}}*/6492 #endif6493 6494 4907 #ifdef _HAVE_DAMAGE_ 6495 /*FUNCTION Tria::CreateKMatrixDamageEvolution {{{*/6496 ElementMatrix* Tria::CreateKMatrixDamageEvolution(void){6497 6498 switch(GetElementType()){6499 case P1Enum: case P2Enum:6500 return CreateKMatrixDamageEvolution_CG();6501 case P1DGEnum:6502 _error_("DG not implemented yet!");break;6503 default:6504 _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");6505 }6506 }6507 /*}}}*/6508 /*FUNCTION Tria::CreateKMatrixDamageEvolution_CG {{{*/6509 ElementMatrix* Tria::CreateKMatrixDamageEvolution_CG(void){6510 6511 /*Intermediaries */6512 int stabilization;6513 int meshtype;6514 IssmDouble Jdet,D_scalar,dt,h;6515 IssmDouble vel,vx,vy,dvxdx,dvydy;6516 IssmDouble dvx[2],dvy[2];6517 IssmDouble xyz_list[NUMVERTICES][3];6518 6519 /*Fetch number of nodes for this finite element*/6520 int numnodes = this->NumberofNodes();6521 6522 /*Initialize Element matrix and vectors*/6523 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);6524 IssmDouble* basis = xNew<IssmDouble>(numnodes);6525 IssmDouble* B = xNew<IssmDouble>(2*numnodes);6526 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);6527 IssmDouble D[2][2];6528 6529 /*Retrieve all inputs and parameters*/6530 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);6531 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);6532 this->parameters->FindParam(&meshtype,MeshTypeEnum);6533 this->parameters->FindParam(&stabilization,DamageStabilizationEnum);6534 Input* vxaverage_input=NULL;6535 Input* vyaverage_input=NULL;6536 if(meshtype==Mesh2DhorizontalEnum){6537 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);6538 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);6539 }6540 else{6541 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);6542 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);6543 }6544 h=sqrt(2*this->GetArea());6545 6546 /* Start looping on the number of gaussian points: */6547 GaussTria *gauss=new GaussTria(2);6548 for(int ig=gauss->begin();ig<gauss->end();ig++){6549 6550 gauss->GaussPoint(ig);6551 6552 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);6553 GetNodalFunctions(basis,gauss);6554 6555 vxaverage_input->GetInputValue(&vx,gauss);6556 vyaverage_input->GetInputValue(&vy,gauss);6557 vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);6558 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);6559 6560 D_scalar=gauss->weight*Jdet;6561 6562 TripleMultiply(basis,1,numnodes,1,6563 &D_scalar,1,1,0,6564 basis,1,numnodes,0,6565 &Ke->values[0],1);6566 GetBMasstransport(B,&xyz_list[0][0],gauss);6567 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);6568 6569 dvxdx=dvx[0];6570 dvydy=dvy[1];6571 D_scalar=dt*gauss->weight*Jdet;6572 6573 D[0][0]=D_scalar*dvxdx;6574 D[0][1]=0.;6575 D[1][0]=0.;6576 D[1][1]=D_scalar*dvydy;6577 TripleMultiply(B,2,numnodes,1,6578 &D[0][0],2,2,0,6579 B,2,numnodes,0,6580 &Ke->values[0],1);6581 6582 D[0][0]=D_scalar*vx;6583 D[1][1]=D_scalar*vy;6584 TripleMultiply(B,2,numnodes,1,6585 &D[0][0],2,2,0,6586 Bprime,2,numnodes,0,6587 &Ke->values[0],1);6588 6589 if(stabilization==2){6590 /*Streamline upwinding*/6591 vel=sqrt(vx*vx+vy*vy)+1.e-8;6592 D[0][0]=h/(2*vel)*vx*vx;6593 D[1][0]=h/(2*vel)*vy*vx;6594 D[0][1]=h/(2*vel)*vx*vy;6595 D[1][1]=h/(2*vel)*vy*vy;6596 }6597 else if(stabilization==1){6598 /*SSA*/6599 vxaverage_input->GetInputAverage(&vx);6600 vyaverage_input->GetInputAverage(&vy);6601 D[0][0]=h/2.0*fabs(vx);6602 D[0][1]=0.;6603 D[1][0]=0.;6604 D[1][1]=h/2.0*fabs(vy);6605 }6606 if(stabilization==1 || stabilization==2){6607 D[0][0]=D_scalar*D[0][0];6608 D[1][0]=D_scalar*D[1][0];6609 D[0][1]=D_scalar*D[0][1];6610 D[1][1]=D_scalar*D[1][1];6611 TripleMultiply(Bprime,2,numnodes,1,6612 &D[0][0],2,2,0,6613 Bprime,2,numnodes,0,6614 &Ke->values[0],1);6615 }6616 }6617 6618 /*Clean up and return*/6619 xDelete<IssmDouble>(basis);6620 xDelete<IssmDouble>(B);6621 xDelete<IssmDouble>(Bprime);6622 delete gauss;6623 return Ke;6624 }6625 /*}}}*/6626 4908 /*FUNCTION Tria::DamageEvolutionF{{{*/ 6627 4909 void Tria::DamageEvolutionF(IssmDouble* f){ … … 6819 5101 } 6820 5102 6821 }6822 /*}}}*/6823 #endif6824 6825 #ifdef _HAVE_BALANCED_6826 /*FUNCTION Tria::CreateKMatrixBalancethickness {{{*/6827 ElementMatrix* Tria::CreateKMatrixBalancethickness(void){6828 6829 switch(GetElementType()){6830 case P1Enum:6831 return CreateKMatrixBalancethickness_CG();6832 case P1DGEnum:6833 return CreateKMatrixBalancethickness_DG();6834 default:6835 _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");6836 }6837 6838 }6839 /*}}}*/6840 /*FUNCTION Tria::CreateKMatrixBalancethickness_CG {{{*/6841 ElementMatrix* Tria::CreateKMatrixBalancethickness_CG(void){6842 6843 /*Intermediaries */6844 int stabilization,meshtype;6845 IssmDouble Jdet,vx,vy,dvxdx,dvydy,vel,h;6846 IssmDouble D_scalar;6847 IssmDouble dvx[2],dvy[2];6848 IssmDouble xyz_list[NUMVERTICES][3];6849 6850 /*Fetch number of nodes and dof for this finite element*/6851 int numnodes = this->NumberofNodes();6852 6853 /*Initialize Element matrix and vectors*/6854 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);6855 IssmDouble* basis = xNew<IssmDouble>(numnodes);6856 IssmDouble* B = xNew<IssmDouble>(2*numnodes);6857 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);6858 IssmDouble D[2][2];6859 6860 /*Retrieve all Inputs and parameters: */6861 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);6862 this->parameters->FindParam(&stabilization,BalancethicknessStabilizationEnum);6863 this->parameters->FindParam(&meshtype,MeshTypeEnum);6864 Input* vxaverage_input=NULL;6865 Input* vyaverage_input=NULL;6866 if(meshtype==Mesh2DhorizontalEnum){6867 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);6868 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);6869 }6870 else{6871 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);6872 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);6873 }6874 h=sqrt(2.*this->GetArea());6875 6876 /*Start looping on the number of gaussian points:*/6877 GaussTria* gauss=new GaussTria(2);6878 for(int ig=gauss->begin();ig<gauss->end();ig++){6879 6880 gauss->GaussPoint(ig);6881 6882 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);6883 GetBMasstransport(B,&xyz_list[0][0],gauss);6884 GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);6885 6886 vxaverage_input->GetInputValue(&vx,gauss);6887 vyaverage_input->GetInputValue(&vy,gauss);6888 vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);6889 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);6890 6891 dvxdx=dvx[0];6892 dvydy=dvy[1];6893 D_scalar=gauss->weight*Jdet;6894 6895 D[0][0]=D_scalar*dvxdx;6896 D[0][1]=0.;6897 D[1][0]=0.;6898 D[1][1]=D_scalar*dvydy;6899 TripleMultiply(B,2,numnodes,1,6900 &D[0][0],2,2,0,6901 B,2,numnodes,0,6902 &Ke->values[0],1);6903 6904 D[0][0]=D_scalar*vx;6905 D[1][1]=D_scalar*vy;6906 TripleMultiply(B,2,numnodes,1,6907 &D[0][0],2,2,0,6908 Bprime,2,numnodes,0,6909 &Ke->values[0],1);6910 6911 if(stabilization==1){6912 /*Streamline upwinding*/6913 vel=sqrt(vx*vx+vy*vy);6914 D[0][0]=h/(2*vel)*vx*vx;6915 D[1][0]=h/(2*vel)*vy*vx;6916 D[0][1]=h/(2*vel)*vx*vy;6917 D[1][1]=h/(2*vel)*vy*vy;6918 }6919 else if(stabilization==2){6920 /*SSA*/6921 vxaverage_input->GetInputAverage(&vx);6922 vyaverage_input->GetInputAverage(&vy);6923 D[0][0]=h/2.0*fabs(vx);6924 D[0][1]=0.;6925 D[1][0]=0.;6926 D[1][1]=h/2.0*fabs(vy);6927 }6928 if(stabilization==1 || stabilization==2){6929 D[0][0]=D_scalar*D[0][0];6930 D[1][0]=D_scalar*D[1][0];6931 D[0][1]=D_scalar*D[0][1];6932 D[1][1]=D_scalar*D[1][1];6933 TripleMultiply(Bprime,2,numnodes,1,6934 &D[0][0],2,2,0,6935 Bprime,2,numnodes,0,6936 &Ke->values[0],1);6937 }6938 }6939 6940 /*Clean up and return*/6941 xDelete<IssmDouble>(basis);6942 xDelete<IssmDouble>(B);6943 xDelete<IssmDouble>(Bprime);6944 delete gauss;6945 return Ke;6946 }6947 /*}}}*/6948 /*FUNCTION Tria::CreateKMatrixBalancethickness_DG {{{*/6949 ElementMatrix* Tria::CreateKMatrixBalancethickness_DG(void){6950 6951 /*Intermediaries*/6952 IssmDouble vx,vy,D_scalar,Jdet;6953 IssmDouble xyz_list[NUMVERTICES][3];6954 6955 /*Fetch number of nodes for this finite element*/6956 int numnodes = this->NumberofNodes();6957 6958 /*Initialize Element matrix and vectors*/6959 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);6960 IssmDouble* basis = xNew<IssmDouble>(numnodes);6961 IssmDouble* B = xNew<IssmDouble>(2*numnodes);6962 IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);6963 IssmDouble D[2][2];6964 6965 /*Retrieve all inputs and parameters*/6966 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);6967 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);6968 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);6969 6970 /*Start looping on the number of gaussian points:*/6971 GaussTria* gauss=new GaussTria(2);6972 for(int ig=gauss->begin();ig<gauss->end();ig++){6973 6974 gauss->GaussPoint(ig);6975 6976 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);6977 /*WARNING: B and Bprime are inverted compared to usual masstransport!!!!*/6978 GetBMasstransport(Bprime,&xyz_list[0][0],gauss);6979 GetBprimeMasstransport(B,&xyz_list[0][0],gauss);6980 6981 vx_input->GetInputValue(&vx,gauss);6982 vy_input->GetInputValue(&vy,gauss);6983 6984 D_scalar=-gauss->weight*Jdet;6985 D[0][0]=D_scalar*vx;6986 D[0][1]=0.;6987 D[1][0]=0.;6988 D[1][1]=D_scalar*vy;6989 6990 TripleMultiply(B,2,numnodes,1,6991 &D[0][0],2,2,0,6992 Bprime,2,numnodes,0,6993 &Ke->values[0],1);6994 }6995 6996 /*Clean up and return*/6997 xDelete<IssmDouble>(basis);6998 xDelete<IssmDouble>(B);6999 xDelete<IssmDouble>(Bprime);7000 delete gauss;7001 return Ke;7002 }7003 /*}}}*/7004 /*FUNCTION Tria::CreateKMatrixBalancevelocity{{{*/7005 ElementMatrix* Tria::CreateKMatrixBalancevelocity(void){7006 7007 /*Intermediaries */7008 IssmDouble xyz_list[NUMVERTICES][3];7009 IssmDouble dhdt_g,mb_g,ms_g,Jdet;7010 IssmDouble h,gamma,thickness;7011 IssmDouble hnx,hny,dhnx[2],dhny[2];7012 7013 /*Fetch number of nodes and dof for this finite element*/7014 int numnodes = this->NumberofNodes();7015 7016 /*Initialize Element matrix and vectors*/7017 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);7018 IssmDouble* B = xNew<IssmDouble>(2*numnodes);7019 IssmDouble* basis = xNew<IssmDouble>(numnodes);7020 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);7021 IssmDouble* HNx = xNew<IssmDouble>(numnodes);7022 IssmDouble* HNy = xNew<IssmDouble>(numnodes);7023 IssmDouble* H = xNew<IssmDouble>(numnodes);7024 IssmDouble* Nx = xNew<IssmDouble>(numnodes);7025 IssmDouble* Ny = xNew<IssmDouble>(numnodes);7026 7027 /*Retrieve all Inputs and parameters: */7028 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);7029 Input* H_input =inputs->GetInput(ThicknessEnum); _assert_(H_input);7030 h=sqrt(2.*this->GetArea());7031 7032 /*Get vector N for all nodes and build HNx and HNy*/7033 GetInputListOnNodes(Nx,SurfaceSlopeXEnum);7034 GetInputListOnNodes(Ny,SurfaceSlopeYEnum);7035 GetInputListOnNodes(H,ThicknessEnum);7036 for(int i=0;i<numnodes;i++){7037 IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);7038 HNx[i] = -H[i]*Nx[i]/norm;7039 HNy[i] = -H[i]*Ny[i]/norm;7040 }7041 7042 /*Start looping on the number of gaussian points:*/7043 GaussTria* gauss=new GaussTria(2);7044 for(int ig=gauss->begin();ig<gauss->end();ig++){7045 7046 gauss->GaussPoint(ig);7047 7048 H_input->GetInputValue(&thickness,gauss);7049 if(thickness<50.) thickness=50.;7050 TriaRef::GetInputDerivativeValue(&dhnx[0],HNx,&xyz_list[0][0],gauss);7051 TriaRef::GetInputDerivativeValue(&dhny[0],HNy,&xyz_list[0][0],gauss);7052 TriaRef::GetInputValue(&hnx,HNx,gauss);7053 TriaRef::GetInputValue(&hny,HNy,gauss);7054 7055 gamma=h/(2.*thickness+1.e-10);7056 7057 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);7058 GetNodalFunctions(basis,gauss);7059 GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);7060 7061 for(int i=0;i<numnodes;i++){7062 for(int j=0;j<numnodes;j++){7063 Ke->values[i*numnodes+j] += gauss->weight*Jdet*(7064 (basis[i]+gamma*(basis[i]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+i]*hnx + dbasis[1*numnodes+i]*hny))*7065 (basis[j]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)7066 );7067 }7068 }7069 }7070 7071 /*Clean up and return*/7072 xDelete<IssmDouble>(basis);7073 xDelete<IssmDouble>(dbasis);7074 xDelete<IssmDouble>(H);7075 xDelete<IssmDouble>(Nx);7076 xDelete<IssmDouble>(Ny);7077 xDelete<IssmDouble>(HNx);7078 xDelete<IssmDouble>(HNy);7079 xDelete<IssmDouble>(B);7080 delete gauss;7081 return Ke;7082 }7083 /*}}}*/7084 /*FUNCTION Tria::CreateKMatrixSmoothedSlope {{{*/7085 ElementMatrix* Tria::CreateKMatrixSmoothedSlope(void){7086 7087 /* Intermediaries */7088 IssmDouble D_scalar,Jdet,thickness;7089 IssmDouble xyz_list[NUMVERTICES][3];7090 IssmDouble l=8.;7091 7092 /*Fetch number of nodes and dof for this finite element*/7093 int numnodes = this->NumberofNodes();7094 7095 /*Initialize Element matrix and vectors*/7096 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters);7097 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes);7098 IssmDouble* basis = xNew<IssmDouble>(numnodes);7099 7100 /*Retrieve all inputs and parameters*/7101 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);7102 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);7103 7104 /* Start looping on the number of gaussian points: */7105 GaussTria* gauss=new GaussTria(2);7106 for(int ig=gauss->begin();ig<gauss->end();ig++){7107 7108 gauss->GaussPoint(ig);7109 7110 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);7111 thickness_input->GetInputValue(&thickness,gauss);7112 if(thickness<50.) thickness=50.;7113 7114 GetNodalFunctions(basis,gauss);7115 GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);7116 7117 for(int i=0;i<numnodes;i++){7118 for(int j=0;j<numnodes;j++){7119 Ke->values[i*numnodes+j] += gauss->weight*Jdet*(7120 basis[i]*basis[j]7121 +(l*thickness)*(l*thickness)*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])7122 );7123 }7124 }7125 }7126 7127 /*Clean up and return*/7128 delete gauss;7129 xDelete<IssmDouble>(dbasis);7130 xDelete<IssmDouble>(basis);7131 return Ke;7132 5103 } 7133 5104 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.
