![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 5739
- Timestamp:
- 09/10/10 09:26:58 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 4 edited
-
Gauss/GaussTria.cpp (modified) (2 diffs)
-
Gauss/GaussTria.h (modified) (1 diff)
-
Loads/Numericalflux.cpp (modified) (9 diffs)
-
Loads/Numericalflux.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Gauss/GaussTria.cpp
r5734 r5739 56 56 weights=(double*)xmalloc(numgauss*sizeof(double)); 57 57 58 /*Reverse grid1 and 2 if necessary*/58 /*Reverse index1 and 2 if necessary*/ 59 59 if (index1>index2){ 60 60 index3=index1; index1=index2; index2=index3; … … 150 150 coord2=ONETHIRD; 151 151 coord3=ONETHIRD; 152 153 } 154 /*}}}*/ 155 /*FUNCTION GaussTria::GaussEdgeCenter{{{1*/ 156 void GaussTria::GaussEdgeCenter(int index1,int index2){ 157 158 int index3; 159 160 /*Reverse index1 and 2 if necessary*/ 161 if (index1>index2){ 162 index3=index1; index1=index2; index2=index3; 163 } 164 165 /*update static arrays*/ 166 if (index1==0 && index2==1){ 167 coord1=0.5; 168 coord2=0.5; 169 coord3=0.0; 170 } 171 else if (index1==0 && index2==2){ 172 coord1=0.5; 173 coord2=0.0; 174 coord3=0.5; 175 } 176 else if (index1==1 && index2==2){ 177 coord1=0.0; 178 coord2=0.5; 179 coord3=0.5; 180 } 181 else 182 ISSMERROR("The 2 indices provided are not supported yet (user provided %i and %i)",index1,index2); 152 183 153 184 } -
issm/trunk/src/c/objects/Gauss/GaussTria.h
r5734 r5739 42 42 void GaussVertex(int iv); 43 43 void GaussCenter(void); 44 void GaussEdgeCenter(int index1,int index2); 44 45 }; 45 46 #endif /* _GAUSSTRIA_H_ */ -
issm/trunk/src/c/objects/Loads/Numericalflux.cpp
r5738 r5739 20 20 21 21 /*Load macros*/ 22 #define NUMVERTICES 4 22 #define NUMVERTICES_INTERNAL 4 23 #define NUMVERTICES_BOUNDARY 2 23 24 #define NDOF1 1 24 25 … … 400 401 401 402 /* constants*/ 402 const int numdof=NDOF1*NUMVERTICES ;403 const int numdof=NDOF1*NUMVERTICES_INTERNAL; 403 404 404 405 /* Intermediaries*/ 405 406 int i,j,ig,index1,index2,analysis_type; 406 407 double DL1,DL2,Jdet,dt,vx,vy,UdotN; 407 double xyz_list[NUMVERTICES ][3];408 double xyz_list[NUMVERTICES_INTERNAL][3]; 408 409 double normal[2]; 409 410 double B[numdof]; 410 411 double Bprime[numdof]; 411 double Ke_g g1[numdof][numdof];412 double Ke_g g2[numdof][numdof];413 double Ke _gg[numdof][numdof] = {0.0};412 double Ke_g1[numdof][numdof]; 413 double Ke_g2[numdof][numdof]; 414 double Ke[numdof][numdof] = {0.0}; 414 415 int *doflist = NULL; 415 416 GaussTria *gauss; 416 417 417 GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES );418 GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES_INTERNAL); 418 419 GetDofList(&doflist); 419 420 … … 454 455 &DL1,1,1,0, 455 456 &Bprime[0],1,numdof,0, 456 &Ke_g g1[0][0],0);457 &Ke_g1[0][0],0); 457 458 TripleMultiply(&B[0],1,numdof,1, 458 459 &DL2,1,1,0, 459 460 &B[0],1,numdof,0, 460 &Ke_g g2[0][0],0);461 462 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke _gg[i][j]+=Ke_gg1[i][j];463 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke _gg[i][j]+=Ke_gg2[i][j];464 } 465 466 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke _gg,ADD_VALUES);461 &Ke_g2[0][0],0); 462 463 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g1[i][j]; 464 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g2[i][j]; 465 } 466 467 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES); 467 468 468 469 /*Free ressources:*/ … … 475 476 void Numericalflux::CreateKMatrixBoundary(Mat Kgg){ 476 477 477 /* local declarations */ 478 int i,j; 479 int analysis_type; 480 481 /* node data: */ 482 const int numgrids=2; 483 const int numdof=NDOF1*numgrids; 484 double xyz_list[numgrids][3]; 485 double normal[2]; 486 int* doflist=NULL; 487 488 /* gaussian points: */ 489 int num_gauss,ig; 490 double* gauss_coords =NULL; 491 double gauss_coord; 492 double* gauss_weights=NULL; 493 double gauss_weight; 494 495 /* matrices: */ 496 double L[numgrids]; 497 double DL; 498 double Ke_gg[numdof][numdof]={0.0}; 499 double Ke_gg_gaussian[numdof][numdof]; 500 double Jdet; 501 502 /*input parameters for structural analysis (diagnostic): */ 503 double vx,vy; 504 double mean_vx; 505 double mean_vy; 506 double UdotN; 507 double dt; 508 int dofs[1]={0}; 509 int dim; 510 511 /*dynamic objects pointed to by hooks: */ 512 Input* vxaverage_input=NULL; 513 Input* vyaverage_input=NULL; 514 515 /*Retrieve parameters: */ 478 /* constants*/ 479 const int numdof=NDOF1*NUMVERTICES_BOUNDARY; 480 481 /* Intermediaries*/ 482 int i,j,ig,index1,index2,analysis_type; 483 double DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN; 484 double xyz_list[NUMVERTICES_BOUNDARY][3]; 485 double normal[2]; 486 double L[numdof]; 487 double Ke_g[numdof][numdof]; 488 double Ke[numdof][numdof] = {0.0}; 489 int *doflist = NULL; 490 GaussTria *gauss; 491 492 GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY); 493 494 /*Retrieve all inputs and parameters we will be needing: */ 495 Tria* tria=(Tria*)element; 496 Input* vxaverage_input=tria->inputs->GetInput(VxEnum); 497 Input* vyaverage_input=tria->inputs->GetInput(VyEnum); 516 498 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); 517 parameters->FindParam(&dim,DimEnum);518 519 /*recover objects from hooks: */520 Tria* tria=(Tria*)element;521 522 /*recover parameters: */523 if (analysis_type==PrognosticAnalysisEnum){524 parameters->FindParam(&dt,DtEnum);525 }526 else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){527 /*No transient term is involved*/528 dt=1;529 }530 else{531 ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type));532 }533 534 /* Get node coordinates and normal vector: */535 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);536 499 GetNormal(&normal[0],xyz_list); 537 538 /*Retrieve all inputs we will be needing: */ 539 if(dim==2){ 540 vxaverage_input=tria->inputs->GetInput(VxEnum); 541 vyaverage_input=tria->inputs->GetInput(VyEnum); 500 switch(analysis_type){ 501 case PrognosticAnalysisEnum: 502 parameters->FindParam(&dt,DtEnum); break; 503 case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum: 504 dt=1; break;/*No transient term is involved*/ 505 default: 506 ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type)); 542 507 } 543 508 544 509 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ 545 this->GetParameterValue(&mean_vx,0.,vxaverage_input); 546 this->GetParameterValue(&mean_vy,0.,vyaverage_input); 510 index1=tria->GetNodeIndex(nodes[0]); 511 index2=tria->GetNodeIndex(nodes[1]); 512 513 gauss=new GaussTria(); 514 gauss->GaussEdgeCenter(index1,index2); 515 vxaverage_input->GetParameterValue(&mean_vx,gauss); 516 vyaverage_input->GetParameterValue(&mean_vy,gauss); 517 delete gauss; 518 547 519 UdotN=mean_vx*normal[0]+mean_vy*normal[1]; 548 520 if (UdotN<=0){ … … 551 523 } 552 524 553 /*Get dof list*/554 525 GetDofList(&doflist); 555 526 556 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */557 num_gauss=2;558 gaussSegment(&gauss_coords, &gauss_weights,num_gauss);559 560 527 /* Start looping on the number of gaussian points: */ 561 for (ig=0; ig<num_gauss; ig++){ 562 563 /*Pick up the gaussian point: */ 564 gauss_weight=gauss_weights[ig]; 565 gauss_coord =gauss_coords[ig]; 566 567 /* Get Jacobian determinant: */ 568 GetJacobianDeterminant(&Jdet,xyz_list,gauss_coord); 569 570 //Get vx, vy and v.n 571 this->GetParameterValue(&vx,gauss_coord,vxaverage_input); 572 this->GetParameterValue(&vy,gauss_coord,vyaverage_input); 573 528 gauss=new GaussTria(index1,index2,2); 529 for(ig=gauss->begin();ig<gauss->end();ig++){ 530 531 gauss->GaussPoint(ig); 532 533 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2); 534 535 vxaverage_input->GetParameterValue(&vx,gauss); 536 vyaverage_input->GetParameterValue(&vy,gauss); 574 537 UdotN=vx*normal[0]+vy*normal[1]; 575 576 /*Get L*/ 577 GetL(&L[0],gauss_coord,numdof); 578 579 /*Compute DLs*/ 580 DL=gauss_weight*Jdet*dt*UdotN; 581 582 /*Do triple product*/ 538 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); 539 DL=gauss->weight*Jdet*dt*UdotN; 540 583 541 TripleMultiply(&L[0],1,numdof,1, 584 542 &DL,1,1,0, 585 543 &L[0],1,numdof,0, 586 &Ke_gg_gaussian[0][0],0); 587 588 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ 589 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 590 591 } // for (ig=0; ig<num_gauss; ig++) 592 593 /*Add Ke_gg to global matrix Kgg: */ 594 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); 595 596 xfree((void**)&gauss_coords); 597 xfree((void**)&gauss_weights); 544 &Ke_g[0][0],0); 545 546 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g[i][j]; 547 } 548 549 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES); 598 550 599 551 /*Free ressources:*/ 552 delete gauss; 600 553 xfree((void**)&doflist); 601 602 554 } 603 555 /*}}}*/ … … 613 565 void Numericalflux::CreatePVectorBoundary(Vec pg){ 614 566 615 /* local declarations */ 616 int i,j; 617 int analysis_type; 618 619 /* node data: */ 620 const int numgrids=2; 621 const int numdof=NDOF1*numgrids; 622 double xyz_list[numgrids][3]; 623 double normal[2]; 624 int* doflist=NULL; 625 626 /* gaussian points: */ 627 int num_gauss,ig; 628 double* gauss_coords =NULL; 629 double gauss_coord; 630 double* gauss_weights=NULL; 631 double gauss_weight; 632 633 /* matrices: */ 634 double L[numgrids]; 635 double DL; 636 double pe_g[numdof]={0.0}; 637 double Jdet; 638 639 /*input parameters for structural analysis (diagnostic): */ 640 double vx,vy; 641 double mean_vx; 642 double mean_vy; 643 double UdotN; 644 double thickness; 645 double dt; 646 int dofs[1]={0}; 647 int dim; 648 649 /*dynamic objects pointed to by hooks: */ 650 Input* vxaverage_input=NULL; 651 Input* vyaverage_input=NULL; 652 Input* thickness_input=NULL; 653 654 /*recover objects from hooks: */ 655 Tria* tria=(Tria*)element; 656 657 /*Retrieve parameters: */ 567 /* constants*/ 568 const int numdof=NDOF1*NUMVERTICES_BOUNDARY; 569 570 /* Intermediaries*/ 571 int i,j,ig,index1,index2,analysis_type; 572 double DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness; 573 double xyz_list[NUMVERTICES_BOUNDARY][3]; 574 double normal[2]; 575 double L[numdof]; 576 double pe[numdof] = {0.0}; 577 int *doflist = NULL; 578 GaussTria *gauss; 579 580 GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY); 581 582 /*Retrieve all inputs and parameters we will be needing: */ 583 Tria* tria=(Tria*)element; 584 Input* vxaverage_input=tria->inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); 585 Input* vyaverage_input=tria->inputs->GetInput(VyEnum); ISSMASSERT(vyaverage_input); 586 Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum); 587 588 /*Here, as it is a forcing, we have H=Hobs by default (for control methods)*/ 589 if (!thickness_input){ 590 thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 591 } 592 658 593 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); 659 parameters->FindParam(&dim,DimEnum);660 661 /*recover parameters: */662 if (analysis_type==PrognosticAnalysisEnum){663 parameters->FindParam(&dt,DtEnum);664 }665 else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){666 /*No transient term is involved*/667 dt=1;668 }669 else{670 ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type));671 }672 673 /*Retrieve all inputs we will be needing: */674 if(dim==2){675 vxaverage_input=tria->inputs->GetInput(VxEnum);676 vyaverage_input=tria->inputs->GetInput(VyEnum);677 }678 /*Here, as it is a forcing, we have H=Hobs by default (for control methods)*/679 thickness_input=tria->inputs->GetInput(ThicknessObsEnum);680 if (!thickness_input) thickness_input=tria->inputs->GetInput(ThicknessEnum);681 682 /* Get node coordinates and normal vector: */683 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);684 594 GetNormal(&normal[0],xyz_list); 595 switch(analysis_type){ 596 case PrognosticAnalysisEnum: 597 parameters->FindParam(&dt,DtEnum); break; 598 case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum: 599 dt=1; break;/*No transient term is involved*/ 600 default: 601 ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type)); 602 } 685 603 686 604 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ 687 this->GetParameterValue(&mean_vx,0.,vxaverage_input); 688 this->GetParameterValue(&mean_vy,0.,vyaverage_input); 605 index1=tria->GetNodeIndex(nodes[0]); 606 index2=tria->GetNodeIndex(nodes[1]); 607 608 gauss=new GaussTria(); 609 gauss->GaussEdgeCenter(index1,index2); 610 vxaverage_input->GetParameterValue(&mean_vx,gauss); 611 vyaverage_input->GetParameterValue(&mean_vy,gauss); 612 delete gauss; 689 613 UdotN=mean_vx*normal[0]+mean_vy*normal[1]; 690 614 if (UdotN>0){ … … 693 617 } 694 618 695 /*Get dof list*/696 619 GetDofList(&doflist); 697 620 698 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */699 num_gauss=2;700 gaussSegment(&gauss_coords, &gauss_weights,num_gauss);701 702 621 /* Start looping on the number of gaussian points: */ 703 for (ig=0; ig<num_gauss; ig++){ 704 705 /*Pick up the gaussian point: */ 706 gauss_weight=gauss_weights[ig]; 707 gauss_coord =gauss_coords[ig]; 708 709 /* Get Jacobian determinant: */ 710 GetJacobianDeterminant(&Jdet,xyz_list,gauss_coord); 711 712 //Get vx, vy and v.n 713 this->GetParameterValue(&vx,gauss_coord,vxaverage_input); 714 this->GetParameterValue(&vy,gauss_coord,vyaverage_input); 715 this->GetParameterValue(&thickness,gauss_coord,thickness_input); 716 622 gauss=new GaussTria(index1,index2,2); 623 for(ig=gauss->begin();ig<gauss->end();ig++){ 624 625 gauss->GaussPoint(ig); 626 627 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2); 628 629 vxaverage_input->GetParameterValue(&vx,gauss); 630 vyaverage_input->GetParameterValue(&vy,gauss); 631 thickness_input->GetParameterValue(&thickness,gauss); 717 632 UdotN=vx*normal[0]+vy*normal[1]; 718 719 /*Get L*/ 720 GetL(&L[0],gauss_coord,numdof); 721 722 /*Compute DL*/ 723 DL= - gauss_weight*Jdet*dt*UdotN*thickness; 724 725 /* Add value into pe_g: */ 726 for( i=0; i<numdof; i++) pe_g[i] += DL* L[i]; 727 728 } // for (ig=0; ig<num_gauss; ig++) 633 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); 634 DL= - gauss->weight*Jdet*dt*UdotN*thickness; 635 636 for(i=0;i<numdof;i++) pe[i] += DL*L[i]; 637 } 729 638 730 639 /*Add pe_g to global matrix Kgg: */ 731 VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); 732 733 xfree((void**)&gauss_coords); 734 xfree((void**)&gauss_weights); 640 VecSetValues(pg,numdof,doflist,(const double*)pe,ADD_VALUES); 735 641 736 642 /*Free ressources:*/ 643 delete gauss; 737 644 xfree((void**)&doflist); 738 739 }740 /*}}}*/741 /*FUNCTION Numericalflux::GetB {{{1*/742 void Numericalflux::GetB(double* B, double gauss_coord){743 744 const int numgrids=4;745 double l1l4[numgrids];746 747 /*Get nodal functions: */748 GetNodalFunctions(&l1l4[0],gauss_coord);749 750 /*Build B: */751 B[0] = +l1l4[0];752 B[1] = +l1l4[1];753 B[2] = -l1l4[0];754 B[3] = -l1l4[1];755 645 } 756 646 /*}}}*/ … … 797 687 } 798 688 /*}}}*/ 799 /*FUNCTION Numericalflux::GetJacobianDeterminant{{{1*/800 void Numericalflux::GetJacobianDeterminant(double* pJdet,double xyz_list[4][3], double gauss_coord){801 802 double Jdet,length;803 804 length=sqrt(pow(xyz_list[1][0] - xyz_list[0][0],2.0) + pow(xyz_list[1][1] - xyz_list[0][1],2.0));805 Jdet=1.0/2.0*length;806 807 if(Jdet<0){808 ISSMERROR(" negative jacobian determinant!");809 }810 811 *pJdet=Jdet;812 }813 /*}}}*/814 /*FUNCTION Numericalflux::GetL {{{1*/815 void Numericalflux::GetL(double* L, double gauss_coord, int numdof){816 817 const int numgrids=4;818 double l1l4[numgrids];819 820 /*Check numdof*/821 ISSMASSERT(numdof==2 || numdof==4);822 823 /*Get nodal functions: */824 GetNodalFunctions(&l1l4[0],gauss_coord);825 826 /*Luild L: */827 L[0] = l1l4[0];828 L[1] = l1l4[1];829 if (numdof==4){830 L[2] = l1l4[2];831 L[3] = l1l4[3];832 }833 }834 /*}}}*/835 /*FUNCTION Numericalflux::GetNodalFunctions{{{1*/836 void Numericalflux::GetNodalFunctions(double* l1l4, double gauss_coord){837 838 /*This routine returns the values of the nodal functions at the gaussian point.*/839 840 l1l4[0]=-0.5*gauss_coord+0.5;841 l1l4[1]=+0.5*gauss_coord+0.5;842 l1l4[2]=-0.5*gauss_coord+0.5;843 l1l4[3]=+0.5*gauss_coord+0.5;844 845 }846 /*}}}*/847 689 /*FUNCTION Numericalflux::GetNormal {{{1*/ 848 690 void Numericalflux:: GetNormal(double* normal,double xyz_list[4][3]){ … … 862 704 } 863 705 /*}}}*/ 864 /*FUNCTION Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,int enumtype) {{{1*/865 void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,int enumtype){866 867 /*output: */868 double value;869 870 /*recover objects from hooks: */871 Tria* tria=(Tria*)element;872 873 /*Get value on Element 1*/874 tria->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,enumtype);875 876 /*Assign output pointers:*/877 *pvalue=value;878 }879 /*}}}*/880 /*FUNCTION Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in) {{{1*/881 void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){882 883 /*output: */884 double value;885 886 /*recover objects from hooks: */887 Tria* tria=(Tria*)element;888 889 /*Get value on Element 1*/890 tria->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,input_in);891 892 /*Assign output pointers:*/893 *pvalue=value;894 }895 /*}}}*/ -
issm/trunk/src/c/objects/Loads/Numericalflux.h
r5737 r5739 72 72 /*}}}*/ 73 73 /*Numericalflux management:{{{1*/ 74 void GetJacobianDeterminant(double* pJdet,double xyz_list[4][3], double gauss_coord);75 void GetNodalFunctions(double* l1l4, double gauss_coord);76 void GetB(double* B, double gauss_coord);77 void GetL(double* L, double gauss_coord,int numdof);78 74 void GetDofList(int** pdoflist); 79 75 void GetNormal(double* normal,double xyz_list[4][3]); 80 void GetParameterValue(double* pvalue, double gauss_coord,int enumtype);81 void GetParameterValue(double* pvalue, double gauss_coord,Input* input_in);82 76 void CreateKMatrixInternal(Mat Kgg); 83 77 void CreateKMatrixBoundary(Mat Kgg);
Note:
See TracChangeset
for help on using the changeset viewer.
