![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 98
- Timestamp:
- 04/28/09 15:25:26 (16 years ago)
- File:
-
- 1 edited
-
issm/trunk/src/c/objects/Tria.cpp (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Tria.cpp
r77 r98 307 307 double B[3][numdofs]; 308 308 double Bprime[3][numdofs]; 309 double L[2][numdofs];310 309 double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}}; // material matrix, simple scalar matrix. 311 310 double D_scalar; 312 double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag313 double DL_scalar;314 311 315 312 /* local element matrices: */ 316 313 double Ke_gg[numdofs][numdofs]; //local element stiffness matrix 317 314 double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point. 318 double Ke_gg_drag_gaussian[numdofs][numdofs]; //stiffness matrix contribution from drag319 315 320 316 double Jdet; 321 317 322 /*slope: */323 double slope[2]={0.0,0.0};324 double slope_magnitude;325 326 318 /*input parameters for structural analysis (diagnostic): */ 327 319 double* velocity_param=NULL; … … 330 322 double oldvxvy_list[numgrids][2]; 331 323 double thickness; 332 333 /*friction: */334 double alpha2_list[numgrids]={0.0,0.0,0.0};335 double alpha2;336 337 double MAXSLOPE=.06; // 6 %338 double MOUNTAINKEXPONENT=12;339 340 324 341 325 /*recover extra inputs from users, at current convergence iteration: */ … … 390 374 391 375 392 /*Build alpha2_list used by drag stiffness matrix*/393 if (!shelf && (friction_type==2)){394 395 /*Allocate friction object: */396 Friction* friction=NewFriction();397 398 /*Initialize all fields: */399 friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));400 strcpy(friction->element_type,"2d");401 402 friction->gravity=matpar->GetG();403 friction->rho_ice=matpar->GetRhoIce();404 friction->rho_water=matpar->GetRhoWater();405 friction->K=&k[0];406 friction->bed=&b[0];407 friction->thickness=&h[0];408 friction->velocities=&vxvy_list[0][0];409 friction->p=p;410 friction->q=q;411 412 /*Compute alpha2_list: */413 FrictionGetAlpha2(&alpha2_list[0],friction);414 415 /*Erase friction object: */416 DeleteFriction(&friction);417 }418 419 #ifdef _DEBUGELEMENTS_420 if(my_rank==RANK && id==ELID){421 printf(" alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);422 }423 #endif424 425 376 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 426 377 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); … … 446 397 /*Compute thickness at gaussian point: */ 447 398 GetParameterValue(&thickness, &h[0],gauss_l1l2l3); 448 449 // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case,450 //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */451 if(!shelf){452 GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);453 slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));454 455 if (slope_magnitude>MAXSLOPE){456 alpha2_list[0]=pow(10,MOUNTAINKEXPONENT);457 alpha2_list[1]=pow(10,MOUNTAINKEXPONENT);458 alpha2_list[2]=pow(10,MOUNTAINKEXPONENT);459 }460 }461 399 462 400 /*Get strain rate from velocity: */ … … 501 439 GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); 502 440 503 /*Get L matrix if viscous basal drag present: */504 if((friction_type==2) && (!shelf)){505 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);506 }507 508 441 /* Do the triple product tB*D*Bprime: */ 509 442 TripleMultiply( &B[0][0],3,numdofs,1, … … 515 448 for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 516 449 517 /*Now, take care of the basal friction if there is any: */518 if((!shelf) && (friction_type==2)){519 520 GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);521 522 if (velocity_param){523 DL_scalar=alpha2*gauss_weight*Jdet;524 }525 else{526 DL_scalar=0;527 }528 529 for (i=0;i<2;i++){530 DL[i][i]=DL_scalar;531 }532 533 /* Do the triple producte tL*D*L: */534 TripleMultiply( &L[0][0],2,numdofs,1,535 &DL[0][0],2,2,0,536 &L[0][0],2,numdofs,0,537 &Ke_gg_drag_gaussian[0][0],0);538 539 for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_drag_gaussian[i][j];540 }541 542 543 450 #ifdef _DEBUGELEMENTS_ 544 451 if(my_rank==RANK && id==ELID){ 545 printf(" alpha2 %g\n",alpha2);546 452 printf(" B:\n"); 547 453 for(i=0;i<3;i++){ … … 558 464 printf("\n"); 559 465 } 560 printf(" L:\n");561 for(i=0;i<2;i++){562 for(j=0;j<numdofs;j++){563 printf("%g ",L[i][j]);564 }565 printf("\n");566 }567 568 466 } 569 467 #endif … … 572 470 /*Add Ke_gg to global matrix Kgg: */ 573 471 MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES); 472 473 474 /*Do not forget to include friction: */ 475 if(!shelf){ 476 CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type); 477 } 478 574 479 575 480 #ifdef _DEBUGELEMENTS_ … … 597 502 598 503 } 599 504 505 #undef __FUNCT__ 506 #define __FUNCT__ "Tria::CreateKMatrixDiagnosticHorizFriction" 507 void Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,ParameterInputs* inputs,int analysis_type){ 508 509 510 /* local declarations */ 511 int i,j; 512 513 /* node data: */ 514 const int numgrids=3; 515 const int numdofs=2*numgrids; 516 double xyz_list[numgrids][3]; 517 int doflist[numdofs]; 518 int numberofdofspernode; 519 520 /* gaussian points: */ 521 int num_gauss,ig; 522 double* first_gauss_area_coord = NULL; 523 double* second_gauss_area_coord = NULL; 524 double* third_gauss_area_coord = NULL; 525 double* gauss_weights = NULL; 526 double gauss_weight; 527 double gauss_l1l2l3[3]; 528 529 /* matrices: */ 530 double L[2][numdofs]; 531 double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag 532 double DL_scalar; 533 534 /* local element matrices: */ 535 double Ke_gg[numdofs][numdofs]; //local element stiffness matrix 536 double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix contribution from drag 537 538 double Jdet; 539 540 /*slope: */ 541 double slope[2]={0.0,0.0}; 542 double slope_magnitude; 543 544 /*input parameters for structural analysis (diagnostic): */ 545 double* velocity_param=NULL; 546 double vxvy_list[numgrids][2]; 547 548 /*friction: */ 549 double alpha2_list[numgrids]={0.0,0.0,0.0}; 550 double alpha2; 551 552 double MAXSLOPE=.06; // 6 % 553 double MOUNTAINKEXPONENT=10; 554 555 /*recover extra inputs from users, at current convergence iteration: */ 556 velocity_param=ParameterInputsRecover(inputs,"velocity"); 557 558 /* Get node coordinates and dof list: */ 559 GetElementNodeData( &xyz_list[0][0], nodes, numgrids); 560 GetDofList(&doflist[0],&numberofdofspernode); 561 562 /* Set Ke_gg to 0: */ 563 for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0; 564 565 if (shelf){ 566 /*no friction, do nothing*/ 567 return; 568 } 569 570 if (friction_type!=2)throw ErrorException(__FUNCT__," non-viscous friction not supported yet!"); 571 572 /*Initialize vxvy_list: */ 573 for(i=0;i<numgrids;i++){ 574 if(velocity_param){ 575 vxvy_list[i][0]=velocity_param[doflist[i*numberofdofspernode+0]]; 576 vxvy_list[i][1]=velocity_param[doflist[i*numberofdofspernode+1]]; 577 } 578 else{ 579 vxvy_list[i][0]=0; 580 vxvy_list[i][1]=0; 581 } 582 } 583 584 /*Build alpha2_list used by drag stiffness matrix*/ 585 Friction* friction=NewFriction(); 586 587 /*Initialize all fields: */ 588 friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); 589 strcpy(friction->element_type,"2d"); 590 591 friction->gravity=matpar->GetG(); 592 friction->rho_ice=matpar->GetRhoIce(); 593 friction->rho_water=matpar->GetRhoWater(); 594 friction->K=&k[0]; 595 friction->bed=&b[0]; 596 friction->thickness=&h[0]; 597 friction->velocities=&vxvy_list[0][0]; 598 friction->p=p; 599 friction->q=q; 600 601 /*Compute alpha2_list: */ 602 FrictionGetAlpha2(&alpha2_list[0],friction); 603 604 /*Erase friction object: */ 605 DeleteFriction(&friction); 606 607 #ifdef _DEBUGELEMENTS_ 608 if(my_rank==RANK && id==ELID){ 609 printf(" alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); 610 } 611 #endif 612 613 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 614 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); 615 616 #ifdef _DEBUGELEMENTS_ 617 if(my_rank==RANK && id==ELID){ 618 printf(" gaussian points: \n"); 619 for(i=0;i<num_gauss;i++){ 620 printf(" %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]); 621 } 622 } 623 #endif 624 625 /* Start looping on the number of gaussian points: */ 626 for (ig=0; ig<num_gauss; ig++){ 627 /*Pick up the gaussian point: */ 628 gauss_weight=*(gauss_weights+ig); 629 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 630 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 631 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 632 633 634 // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, 635 //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ 636 GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3); 637 slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); 638 639 if (slope_magnitude>MAXSLOPE){ 640 alpha2_list[0]=pow(10,MOUNTAINKEXPONENT); 641 alpha2_list[1]=pow(10,MOUNTAINKEXPONENT); 642 alpha2_list[2]=pow(10,MOUNTAINKEXPONENT); 643 } 644 645 /* Get Jacobian determinant: */ 646 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); 647 648 /*Get L matrix: */ 649 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); 650 651 /*Now, take care of the basal friction if there is any: */ 652 GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3); 653 654 if (velocity_param){ 655 DL_scalar=alpha2*gauss_weight*Jdet; 656 } 657 else{ 658 DL_scalar=0; 659 } 660 661 for (i=0;i<2;i++){ 662 DL[i][i]=DL_scalar; 663 } 664 665 /* Do the triple producte tL*D*L: */ 666 TripleMultiply( &L[0][0],2,numdofs,1, 667 &DL[0][0],2,2,0, 668 &L[0][0],2,numdofs,0, 669 &Ke_gg_gaussian[0][0],0); 670 671 for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 672 673 } // for (ig=0; ig<num_gauss; ig++) 674 675 /*Add Ke_gg to global matrix Kgg: */ 676 MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES); 677 678 cleanup_and_return: 679 xfree((void**)&first_gauss_area_coord); 680 xfree((void**)&second_gauss_area_coord); 681 xfree((void**)&third_gauss_area_coord); 682 xfree((void**)&gauss_weights); 683 684 } 600 685 601 686 #undef __FUNCT__ … … 1928 2013 return Jelem; 1929 2014 } 2015 2016 #undef __FUNCT__ 2017 #define __FUNCT__ "Tria::NodeConfiguration" 2018 void Tria::NodeConfiguration(int* tria_node_ids,Node* tria_nodes[3],int* tria_node_offsets){ 2019 2020 int i; 2021 for(i=0;i<3;i++){ 2022 node_ids[i]=tria_node_ids[i]; 2023 nodes[i]=tria_nodes[i]; 2024 node_offsets[i]=tria_node_offsets[i]; 2025 } 2026 2027 } 2028 #undef __FUNCT__ 2029 #define __FUNCT__ "Tria::MaticeConfiguration" 2030 void Tria::MaticeConfiguration(Matice* tria_matice,int tria_matice_offset){ 2031 matice=tria_matice; 2032 matice_offset=tria_matice_offset; 2033 } 2034 2035 #undef __FUNCT__ 2036 #define __FUNCT__ "Tria::MaticeConfiguration" 2037 void Tria::MatparConfiguration(Matpar* tria_matpar,int tria_matpar_offset){ 2038 2039 matpar=tria_matpar; 2040 matpar_offset=tria_matpar_offset; 2041 2042 } 2043
Note:
See TracChangeset
for help on using the changeset viewer.
