Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 5718) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 5719) @@ -32,4 +32,5 @@ virtual void GetSolutionFromInputs(Vec solution)=0; virtual void GetNodes(void** nodes)=0; + virtual int GetNodeIndex(Node* node)=0; virtual bool GetShelf()=0; virtual bool GetOnBed()=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5718) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5719) @@ -916,4 +916,17 @@ pnodes[i]=nodes[i]; } +} +/*}}}*/ +/*FUNCTION Penta::GetNodeIndex {{{1*/ +int Penta::GetNodeIndex(Node* node){ + + int i; + + for(i=0;i=0 && index1<3); + ISSMASSERT(index2>=0 && index2<3); + l1l2[0]=BasisFunctions[index1]; + l1l2[1]=BasisFunctions[index2]; +} +/*}}}*/ /*FUNCTION TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss) {{{1*/ void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){ Index: /issm/trunk/src/c/objects/Elements/TriaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaRef.h (revision 5718) +++ /issm/trunk/src/c/objects/Elements/TriaRef.h (revision 5719) @@ -37,4 +37,6 @@ void GetJacobian(double* J, double* xyz_list,double* gauss); void GetJacobian(double* J, double* xyz_list,GaussTria* gauss); + void GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss); + void GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss); void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss); void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss); @@ -45,4 +47,5 @@ void GetNodalFunctions(double* l1l2l3, double* gauss); void GetNodalFunctions(double* l1l2l3,GaussTria* gauss); + void GetSegmentNodalFunctions(double* l1l2l3,GaussTria* gauss, int index1,int index2); void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss); void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, GaussTria* gauss); Index: /issm/trunk/src/c/objects/Gauss/GaussTria.cpp =================================================================== --- /issm/trunk/src/c/objects/Gauss/GaussTria.cpp (revision 5718) +++ /issm/trunk/src/c/objects/Gauss/GaussTria.cpp (revision 5719) @@ -36,4 +36,61 @@ coord3=UNDEF; +} +/*}}}*/ +/*FUNCTION GaussTria::GaussTria(int index1,int index2,int order) {{{1*/ +GaussTria::GaussTria(int index1,int index2,int order){ + + /*Intermediaties*/ + double *seg_coords = NULL; + double *seg_weights = NULL; + int i,index3; + + /*Get Segment gauss points*/ + numgauss=order; + GaussLegendreLinear(&seg_coords,&seg_weights,numgauss); + + /*Allocate GaussTria fields*/ + coords1=(double*)xmalloc(numgauss*sizeof(double)); + coords2=(double*)xmalloc(numgauss*sizeof(double)); + coords3=(double*)xmalloc(numgauss*sizeof(double)); + weights=(double*)xmalloc(numgauss*sizeof(double)); + + /*Reverse grid1 and 2 if necessary*/ + if (index1>index2){ + index3=index1; index1=index2; index2=index3; + for(i=0;iinputs->GetInput(ThicknessEnum); - bed_input =((Tria*)element)->inputs->GetInput(BedEnum); - inputs->GetParameterValue(&fill,FillEnum); + + /*Intermediary*/ + int ig,index1,index2,fill; + double Jdet; + double thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity; + double water_pressure,air_pressure,surface_under_water,base_under_water; + int *doflist = NULL; + double xyz_list[numnodes][3]; + double pe_g[numdofs] = {0.0}; + double normal[2]; + double L[2]; + GaussTria *gauss; + + Tria* tria=((Tria*)element); /* Get dof list and node coordinates: */ GetDofList(&doflist,MacAyealApproximationEnum); GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes); - - /*Get Matpar parameters*/ + + /*Recover inputs and parameters */ + Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + Input* bed_input =tria->inputs->GetInput(BedEnum); ISSMASSERT(bed_input); + inputs->GetParameterValue(&fill,FillEnum); + rho_water=matpar->GetRhoWater(); rho_ice =matpar->GetRhoIce(); gravity =matpar->GetG(); - /*Get normal*/ GetNormal(&normal[0],xyz_list); - //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions - num_gauss=3; - gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss); - - for(ig=0;igGetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord); - - /*Get thickness and bed*/ - this->GetParameterValue(&thickness,gauss_coord,thickness_input); - this->GetParameterValue(&bed, gauss_coord,bed_input); - - /*Compute total pressure applied to ice front*/ - if (fill==WaterEnum){ - //icefront ends in water: - ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); - air_pressure=0; - - //Now deal with water pressure - surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level - base_under_water=min(0,bed); // 0 if the bottom of the glacier is above water level - water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); + index1=tria->GetNodeIndex(nodes[0]); + index2=tria->GetNodeIndex(nodes[1]); + gauss=new GaussTria(index1,index2,3); + + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + thickness_input->GetParameterValue(&thickness,gauss); + bed_input->GetParameterValue(&bed,gauss); + + switch(fill){ + case WaterEnum: + surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level + base_under_water=min(0,bed); // 0 if the bottom of the glacier is above water level + water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2)); + break; + case AirEnum: + water_pressure=0; + break; + default: + ISSMERROR("fill type %s not supported yet",EnumToString(fill)); } - else if (fill==AirEnum){ - ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); - air_pressure=0; - water_pressure=0; + ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); + air_pressure=0; + pressure = ice_pressure + water_pressure + air_pressure; + + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2); + + for (int i=0;iweight*normal[0]*L[i]; + pe_g[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*L[i]; } - else{ - ISSMERROR("fill type %s not supported yet",EnumToString(fill)); - } - pressure = ice_pressure + water_pressure + air_pressure; - - /*Get L*/ - beam->GetNodalFunctions(&L[0],gauss_coord); - - for (i=0;i