Changeset 5719
- Timestamp:
- 09/09/10 12:00:59 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Element.h
r5660 r5719 32 32 virtual void GetSolutionFromInputs(Vec solution)=0; 33 33 virtual void GetNodes(void** nodes)=0; 34 virtual int GetNodeIndex(Node* node)=0; 34 35 virtual bool GetShelf()=0; 35 36 virtual bool GetOnBed()=0; -
issm/trunk/src/c/objects/Elements/Penta.cpp
r5707 r5719 916 916 pnodes[i]=nodes[i]; 917 917 } 918 } 919 /*}}}*/ 920 /*FUNCTION Penta::GetNodeIndex {{{1*/ 921 int Penta::GetNodeIndex(Node* node){ 922 923 int i; 924 925 for(i=0;i<NUMVERTICES;i++){ 926 if(node==nodes[i]) 927 return i; 928 } 929 ISSMERROR("Node provided not found among element nodes"); 930 918 931 } 919 932 /*}}}*/ -
issm/trunk/src/c/objects/Elements/Penta.h
r5660 r5719 79 79 void DeleteResults(void); 80 80 void GetNodes(void** nodes); 81 int GetNodeIndex(Node* node); 81 82 bool GetOnBed(); 82 83 bool GetShelf(); -
issm/trunk/src/c/objects/Elements/Tria.cpp
r5669 r5719 820 820 } 821 821 /*}}}*/ 822 /*FUNCTION Tria::GetNodeIndex {{{1*/ 823 int Tria::GetNodeIndex(Node* node){ 824 825 int i; 826 827 for(i=0;i<NUMVERTICES;i++){ 828 if(node==nodes[i]) 829 return i; 830 } 831 ISSMERROR("Node provided not found among element nodes"); 832 } 833 /*}}}*/ 822 834 /*FUNCTION Tria::GetOnBed {{{1*/ 823 835 bool Tria::GetOnBed(){ -
issm/trunk/src/c/objects/Elements/Tria.h
r5669 r5719 75 75 void CreatePVector(Vec pg); 76 76 void GetNodes(void** nodes); 77 int GetNodeIndex(Node* node); 77 78 bool GetOnBed(); 78 79 bool GetShelf(); -
issm/trunk/src/c/objects/Elements/TriaRef.cpp
r5638 r5719 432 432 *(J+NDOF2*0+1)=0.5*(y2-y1); 433 433 *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2); 434 } 435 /*}}}*/ 436 /*FUNCTION TriaRef::GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss) {{{1*/ 437 void TriaRef::GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss){ 438 /*The Jacobian is constant over the element, discard the gaussian points. 439 * J is assumed to have been allocated.*/ 440 441 double x1,y1,x2,y2; 442 443 x1=*(xyz_list+3*0+0); 444 y1=*(xyz_list+3*0+1); 445 x2=*(xyz_list+3*1+0); 446 y2=*(xyz_list+3*1+1); 447 448 *J=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.)); 449 } 450 /*}}}*/ 451 /*FUNCTION TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss) {{{1*/ 452 void TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss){ 453 /*The Jacobian determinant is constant over the element, discard the gaussian points. 454 * J is assumed to have been allocated*/ 455 456 GetSegmentJacobian(Jdet,xyz_list,gauss); 457 if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); 458 434 459 } 435 460 /*}}}*/ … … 557 582 } 558 583 /*}}}*/ 584 /*FUNCTION TriaRef::GetSegmentNodalFunctions(double* l1l2,GaussTria* gauss) {{{1*/ 585 void TriaRef::GetSegmentNodalFunctions(double* l1l2,GaussTria* gauss,int index1,int index2){ 586 /*This routine returns the values of the nodal functions at the gaussian point.*/ 587 588 double BasisFunctions[3]; 589 590 GetNodalFunctions(&BasisFunctions[0],gauss); 591 592 ISSMASSERT(index1>=0 && index1<3); 593 ISSMASSERT(index2>=0 && index2<3); 594 l1l2[0]=BasisFunctions[index1]; 595 l1l2[1]=BasisFunctions[index2]; 596 } 597 /*}}}*/ 559 598 /*FUNCTION TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss) {{{1*/ 560 599 void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){ -
issm/trunk/src/c/objects/Elements/TriaRef.h
r5638 r5719 37 37 void GetJacobian(double* J, double* xyz_list,double* gauss); 38 38 void GetJacobian(double* J, double* xyz_list,GaussTria* gauss); 39 void GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss); 40 void GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss); 39 41 void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss); 40 42 void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss); … … 45 47 void GetNodalFunctions(double* l1l2l3, double* gauss); 46 48 void GetNodalFunctions(double* l1l2l3,GaussTria* gauss); 49 void GetSegmentNodalFunctions(double* l1l2l3,GaussTria* gauss, int index1,int index2); 47 50 void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss); 48 51 void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, GaussTria* gauss); -
issm/trunk/src/c/objects/Gauss/GaussTria.cpp
r5641 r5719 36 36 coord3=UNDEF; 37 37 38 } 39 /*}}}*/ 40 /*FUNCTION GaussTria::GaussTria(int index1,int index2,int order) {{{1*/ 41 GaussTria::GaussTria(int index1,int index2,int order){ 42 43 /*Intermediaties*/ 44 double *seg_coords = NULL; 45 double *seg_weights = NULL; 46 int i,index3; 47 48 /*Get Segment gauss points*/ 49 numgauss=order; 50 GaussLegendreLinear(&seg_coords,&seg_weights,numgauss); 51 52 /*Allocate GaussTria fields*/ 53 coords1=(double*)xmalloc(numgauss*sizeof(double)); 54 coords2=(double*)xmalloc(numgauss*sizeof(double)); 55 coords3=(double*)xmalloc(numgauss*sizeof(double)); 56 weights=(double*)xmalloc(numgauss*sizeof(double)); 57 58 /*Reverse grid1 and 2 if necessary*/ 59 if (index1>index2){ 60 index3=index1; index1=index2; index2=index3; 61 for(i=0;i<numgauss;i++) seg_coords[i]=-seg_coords[i]; 62 } 63 64 /*Build Triangle Gauss point*/ 65 if (index1==0 && index2==1){ 66 for(i=0;i<numgauss;i++) coords1[i]= 0.5*(1-seg_coords[i]); 67 for(i=0;i<numgauss;i++) coords2[i]=1-0.5*(1.-seg_coords[i]); 68 for(i=0;i<numgauss;i++) coords3[i]=0; 69 for(i=0;i<numgauss;i++) weights[i]=seg_weights[i]; 70 } 71 else if (index1==0 && index2==2){ 72 for(i=0;i<numgauss;i++) coords1[i]= 0.5*(1-seg_coords[i]); 73 for(i=0;i<numgauss;i++) coords2[i]=0; 74 for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]); 75 for(i=0;i<numgauss;i++) weights[i]=seg_weights[i]; 76 } 77 else if (index1==1 && index2==2){ 78 for(i=0;i<numgauss;i++) coords1[i]=0; 79 for(i=0;i<numgauss;i++) coords2[i]= 0.5*(1-seg_coords[i]); 80 for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]); 81 for(i=0;i<numgauss;i++) weights[i]=seg_weights[i]; 82 } 83 else 84 ISSMERROR("The 2 indices provided are not supported yet (user provided %i and %i)",index1,index2); 85 86 /*Initialize static fields as undefined*/ 87 weight=UNDEF; 88 coord1=UNDEF; 89 coord2=UNDEF; 90 coord3=UNDEF; 91 92 /*clean up*/ 93 xfree((void**)&seg_coords); 94 xfree((void**)&seg_weights); 38 95 } 39 96 /*}}}*/ -
issm/trunk/src/c/objects/Gauss/GaussTria.h
r5641 r5719 31 31 GaussTria(); 32 32 GaussTria(int order); 33 GaussTria(int index1,int index2,int order); 33 34 ~GaussTria(); 34 35 -
issm/trunk/src/c/objects/Loads/Icefront.cpp
r5715 r5719 425 425 void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg){ 426 426 427 int i;427 /*Constants*/ 428 428 const int numnodes= NUMVERTICESSEG; 429 429 const int numdofs = numnodes *NDOF2; 430 int* doflist=NULL; 431 double xyz_list[numnodes][3]; 432 433 /*Objects: */ 434 double pe_g[numdofs] = {0.0}; 435 436 /*Segment*/ 437 double normal[2]; 438 int num_gauss; 439 double* segment_gauss_coord=NULL; 440 double* gauss_weights=NULL; 441 double gauss_weight; 442 double gauss_coord; 443 int ig; 444 double Jdet; 445 double L[2]; 446 double thickness,bed,pressure; 447 double ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity; 448 double surface_under_water,base_under_water; 449 int fill; 450 451 /*Inputs*/ 452 Input* thickness_input=NULL; 453 Input* bed_input=NULL; 454 Tria* tria=NULL; 455 456 BeamRef* beam=NULL; 457 458 /*Recover inputs: */ 459 thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum); 460 bed_input =((Tria*)element)->inputs->GetInput(BedEnum); 461 inputs->GetParameterValue(&fill,FillEnum); 430 431 /*Intermediary*/ 432 int ig,index1,index2,fill; 433 double Jdet; 434 double thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity; 435 double water_pressure,air_pressure,surface_under_water,base_under_water; 436 int *doflist = NULL; 437 double xyz_list[numnodes][3]; 438 double pe_g[numdofs] = {0.0}; 439 double normal[2]; 440 double L[2]; 441 GaussTria *gauss; 442 443 Tria* tria=((Tria*)element); 462 444 463 445 /* Get dof list and node coordinates: */ 464 446 GetDofList(&doflist,MacAyealApproximationEnum); 465 447 GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes); 466 467 /*Get Matpar parameters*/ 448 449 /*Recover inputs and parameters */ 450 Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 451 Input* bed_input =tria->inputs->GetInput(BedEnum); ISSMASSERT(bed_input); 452 inputs->GetParameterValue(&fill,FillEnum); 453 468 454 rho_water=matpar->GetRhoWater(); 469 455 rho_ice =matpar->GetRhoIce(); 470 456 gravity =matpar->GetG(); 471 457 472 /*Get normal*/473 458 GetNormal(&normal[0],xyz_list); 474 459 475 //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions 476 num_gauss=3; 477 gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss); 478 479 for(ig=0;ig<num_gauss;ig++){ 480 481 /*Pick up the gaussian point: */ 482 gauss_weight=gauss_weights[ig]; 483 gauss_coord=segment_gauss_coord[ig]; 484 485 //compute Jacobian of segment 486 beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord); 487 488 /*Get thickness and bed*/ 489 this->GetParameterValue(&thickness,gauss_coord,thickness_input); 490 this->GetParameterValue(&bed, gauss_coord,bed_input); 491 492 /*Compute total pressure applied to ice front*/ 493 if (fill==WaterEnum){ 494 //icefront ends in water: 495 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 496 air_pressure=0; 497 498 //Now deal with water pressure 499 surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level 500 base_under_water=min(0,bed); // 0 if the bottom of the glacier is above water level 501 water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); 460 index1=tria->GetNodeIndex(nodes[0]); 461 index2=tria->GetNodeIndex(nodes[1]); 462 gauss=new GaussTria(index1,index2,3); 463 464 for(ig=gauss->begin();ig<gauss->end();ig++){ 465 466 gauss->GaussPoint(ig); 467 468 thickness_input->GetParameterValue(&thickness,gauss); 469 bed_input->GetParameterValue(&bed,gauss); 470 471 switch(fill){ 472 case WaterEnum: 473 surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level 474 base_under_water=min(0,bed); // 0 if the bottom of the glacier is above water level 475 water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); 476 break; 477 case AirEnum: 478 water_pressure=0; 479 break; 480 default: 481 ISSMERROR("fill type %s not supported yet",EnumToString(fill)); 502 482 } 503 else if (fill==AirEnum){ 504 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 505 air_pressure=0; 506 water_pressure=0; 483 ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); 484 air_pressure=0; 485 pressure = ice_pressure + water_pressure + air_pressure; 486 487 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); 488 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2); 489 490 for (int i=0;i<numnodes;i++){ 491 pe_g[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*L[i]; 492 pe_g[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*L[i]; 507 493 } 508 else{ 509 ISSMERROR("fill type %s not supported yet",EnumToString(fill)); 510 } 511 pressure = ice_pressure + water_pressure + air_pressure; 512 513 /*Get L*/ 514 beam->GetNodalFunctions(&L[0],gauss_coord); 515 516 for (i=0;i<numnodes;i++){ 517 pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i]; 518 pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i]; 519 } 520 } //for(ig=0;ig<num_gauss;ig++) 521 522 523 /*Plug pe_g into vector: */ 494 } 495 524 496 VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); 525 497 526 498 /*Free ressources:*/ 527 xfree((void**)&segment_gauss_coord); 528 xfree((void**)&gauss_weights); 499 delete gauss; 529 500 xfree((void**)&doflist); 530 531 501 } 532 502 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.