Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.cpp (revision 4945) +++ /issm/trunk/src/c/objects/Loads/Icefront.cpp (revision 4946) @@ -397,91 +397,123 @@ void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg){ - int i,j; - - const int numgrids=2; - const int NDOF2=2; - int numberofdofspernode; - const int numdofs=numgrids*NDOF2; - int doflist[numdofs]; - double xyz_list[numgrids][3]; - - double x1,y1,x2,y2; - - /* input parameters: */ - double thickness_list_element[6]; - double thickness_list[2]; - double bed_list_element[6]; - double bed_list[2]; - - int grid1,grid2; - double normal[2]; - double length; - int element_type; + int i; + const int numgrids = 2; + const int NDOF2 = 2; + const int numdofs = numgrids *NDOF2; + int numberofdofspernode; + int doflist[numdofs]; + double xyz_list[numgrids][3]; /*Objects: */ - double pe_g[numdofs]={0.0}; - Matpar* matpar=NULL; - Node** element_nodes=NULL; - Node** nodes=NULL; - Element* element=NULL; + double pe_g[numdofs] = {0.0}; + Matpar *matpar = NULL; + Node **nodes = NULL; + Element *element = NULL; + + /*Segment*/ + double normal[2]; + int num_gauss; + double* segment_gauss_coord=NULL; + double* gauss_weights=NULL; + double gauss_weight; + double gauss_coord; + int ig; + double Jdet; + double L[2]; + double thickness,bed,pressure; + double ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity; + double surface_under_water,base_under_water; + int fill; + + /*Inputs*/ + Input* thickness_input=NULL; + Input* bed_input=NULL; + Tria* tria=NULL; + Penta* penta=NULL; /*Recover hook objects: */ - matpar=(Matpar*)hmatpar->delivers(); + matpar =(Matpar*)hmatpar->delivers(); element=(Element*)helement->delivers(); - nodes=(Node**)hnodes->deliverp(); - - element_type=element->Enum(); + nodes =(Node**)hnodes->deliverp(); + BeamRef* beam=NULL; //check that the element is onbed (collapsed formulation) otherwise:pe=0 - if(element_type==PentaEnum){ - if (!element->GetOnBed()){ - return; - } - } - - /*Identify which grids are comprised in the segment. */ - if(element_type==TriaEnum)element_nodes=(Node**)xmalloc(3*sizeof(Node*)); - if(element_type==PentaEnum)element_nodes=(Node**)xmalloc(6*sizeof(Node*)); - element->GetNodes((void**)element_nodes); - - /*go through first 3 grids (all grids for tria, bed grids for penta) and identify 1st and 2nd grid: */ - for(i=0;i<3;i++){ - if (nodes[0]==element_nodes[i])grid1=i; - if (nodes[1]==element_nodes[i])grid2=i; - } - + if(element->Enum()==PentaEnum){ + if(!element->GetOnBed()) return; + thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum); + bed_input =((Penta*)element)->inputs->GetInput(BedEnum); + } + else{ + thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum); + bed_input =((Tria*)element)->inputs->GetInput(BedEnum); + } + + /*Recover inputs: */ + inputs->GetParameterValue(&fill,FillEnum); + /* Get dof list and node coordinates: */ GetDofList(&doflist[0],&numberofdofspernode); - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - - /*Now build thickness_list_element and bed_list_element: */ - element->GetThicknessList(&thickness_list_element[0]); - element->GetBedList(&bed_list_element[0]); - - /*Build thickness_list and bed_list: */ - thickness_list[0]=thickness_list_element[grid1]; - thickness_list[1]=thickness_list_element[grid2]; - bed_list[0]=bed_list_element[grid1]; - bed_list[1]=bed_list_element[grid2]; - - //Recover grid coordinates - x1=xyz_list[0][0]; - y1=xyz_list[0][1]; - x2=xyz_list[1][0]; - y2=xyz_list[1][1]; - - //Compute length and normal of segment - normal[0]=cos(atan2(x1-x2,y2-y1)); - normal[1]=sin(atan2(x1-x2,y2-y1)); - length=sqrt(pow(x2-x1,2)+pow(y2-y1,2)); - - //Compute load contribution for this segment: - SegmentPressureLoad(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list,bed_list,normal,length); + GetVerticesCoordinates(&xyz_list[0][0],nodes,numgrids); + + /*Get Matpar parameters*/ + 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)); + } + else if (fill==AirEnum){ + ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2); + air_pressure=0; + water_pressure=0; + } + else{ + ISSMERROR("fill type %i not supported yet",fill); + } + pressure = ice_pressure + water_pressure + air_pressure; + + /*Get L*/ + beam->GetNodalFunctions(&L[0],gauss_coord); + + for (i=0;iGetParameterValue(&fill,FillEnum); - - nx=normal[0]; - ny=normal[1]; - - - //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); - - //recover thickness and bed at two grids - h1=thickness_list[0]; - h2=thickness_list[1]; - b1=bed_list[0]; - b2=bed_list[1]; - - //compute Jacobian of segment - Jdet=1./2.*length; - - for(ig=0;igdelivers(); + nodes=(Node**)hnodes->deliverp(); + + /*Get value on Element 1*/ + element->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,input_in); + + /*Assign output pointers:*/ + *pvalue=value; +} +/*}}}*/ +/*FUNCTION Icefront::GetNormal {{{1*/ +void Icefront:: GetNormal(double* normal,double xyz_list[2][3]){ + + /*Build unit outward pointing vector*/ + const int numgrids=2; + double vector[2]; + double norm; + + vector[0]=xyz_list[1][0] - xyz_list[0][0]; + vector[1]=xyz_list[1][1] - xyz_list[0][1]; + + norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0)); + + normal[0]= + vector[1]/norm; + normal[1]= - vector[0]/norm; +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Icefront.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.h (revision 4945) +++ /issm/trunk/src/c/objects/Loads/Icefront.h (revision 4946) @@ -73,9 +73,10 @@ void CreatePVectorDiagnosticStokes( Vec pg); void GetDofList(int* doflist,int* pnumberofdofs); - void SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length); void QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list); void QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list); + void GetParameterValue(double* pvalue, double gauss_coord,Input* input_in); + void GetNormal(double* normal,double xyz_list[2][3]); /*}}}*/ };