Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18910) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18911) @@ -155,4 +155,97 @@ } /*}}}*/ +void Tria::AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/ + + bool already = false; + int i,j; + int partition[NUMVERTICES]; + int offsetsid[NUMVERTICES]; + int offsetdof[NUMVERTICES]; + IssmDouble area; + IssmDouble mean; + + /*First, get the area: */ + area=this->GetArea(); + + /*Figure out the average for this element: */ + this->GetVerticesSidList(&offsetsid[0]); + this->GetVertexPidList(&offsetdof[0]); + mean=0; + for(i=0;i(qmu_part[offsetsid[i]]); + mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]]; + } + + /*Add contribution: */ + for(i=0;iSetValue(partition[i],mean*area,ADD_VAL); + partition_areas->SetValue(partition[i],area,ADD_VAL); + }; + } +} +/*}}}*/ +void Tria::CalvingRateLevermann(){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + GaussTria* gauss=NULL; + IssmDouble vx,vy,vel; + IssmDouble strainparallel; + IssmDouble propcoeff; + IssmDouble strainperpendicular; + IssmDouble calvingratex[NUMVERTICES]; + IssmDouble calvingratey[NUMVERTICES]; + IssmDouble calvingrate[NUMVERTICES]; + + + /* Get node coordinates and dof list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Retrieve all inputs and parameters we will need*/ + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum); _assert_(strainparallel_input); + Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input); + Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum); _assert_(levermanncoeff_input); + + /* Start looping on the number of vertices: */ + gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + /* Get the value we need*/ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel=vx*vx+vy*vy; + strainparallel_input->GetInputValue(&strainparallel,gauss); + strainperpendicular_input->GetInputValue(&strainperpendicular,gauss); + levermanncoeff_input->GetInputValue(&propcoeff,gauss); + + /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */ + calvingrate[iv]=propcoeff*strainparallel*strainperpendicular; + if(calvingrate[iv]<0){ + calvingrate[iv]=0; + } + calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6); + calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6); + } + + /*Add input*/ + this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum)); + this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum)); + this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; + +} +/*}}}*/ IssmDouble Tria::CharacteristicLength(void){/*{{{*/ @@ -162,4 +255,54 @@ void Tria::ComputeBasalStress(Vector* eps){/*{{{*/ _error_("Not Implemented yet"); +} +/*}}}*/ +void Tria::ComputeDeviatoricStressTensor(){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble viscosity; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble tau_xx[NUMVERTICES]; + IssmDouble tau_yy[NUMVERTICES]; + IssmDouble tau_zz[NUMVERTICES]={0,0,0}; + IssmDouble tau_xy[NUMVERTICES]; + IssmDouble tau_xz[NUMVERTICES]={0,0,0}; + IssmDouble tau_yz[NUMVERTICES]={0,0,0}; + GaussTria* gauss=NULL; + int domaintype,dim=2; + + /* Get node coordinates and dof list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Retrieve all inputs we will be needing: */ + this->FindParam(&domaintype,DomainTypeEnum); + if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of vertices: */ + gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + /*Compute strain rate and viscosity: */ + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input); + + /*Compute Stress*/ + tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps + tau_yy[iv]=2*viscosity*epsilon[1]; + tau_xy[iv]=2*viscosity*epsilon[2]; + } + + /*Add Stress tensor components into inputs*/ + this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum)); + this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum)); + this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum)); + this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum)); + this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum)); + this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; } /*}}}*/ @@ -276,201 +419,4 @@ } /*}}}*/ -void Tria::ComputeDeviatoricStressTensor(){/*{{{*/ - - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble viscosity; - IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - IssmDouble tau_xx[NUMVERTICES]; - IssmDouble tau_yy[NUMVERTICES]; - IssmDouble tau_zz[NUMVERTICES]={0,0,0}; - IssmDouble tau_xy[NUMVERTICES]; - IssmDouble tau_xz[NUMVERTICES]={0,0,0}; - IssmDouble tau_yz[NUMVERTICES]={0,0,0}; - GaussTria* gauss=NULL; - int domaintype,dim=2; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*Retrieve all inputs we will be needing: */ - this->FindParam(&domaintype,DomainTypeEnum); - if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /* Start looping on the number of vertices: */ - gauss=new GaussTria(); - for (int iv=0;ivGaussVertex(iv); - - /*Compute strain rate and viscosity: */ - this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input); - - /*Compute Stress*/ - tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps - tau_yy[iv]=2*viscosity*epsilon[1]; - tau_xy[iv]=2*viscosity*epsilon[2]; - } - - /*Add Stress tensor components into inputs*/ - this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum)); - this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum)); - this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum)); - this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum)); - this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum)); - this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum)); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ -void Tria::StrainRateparallel(){/*{{{*/ - - IssmDouble *xyz_list = NULL; - IssmDouble epsilon[3]; - GaussTria* gauss=NULL; - IssmDouble vx,vy,vel; - IssmDouble strainxx; - IssmDouble strainxy; - IssmDouble strainyy; - IssmDouble strainparallel[NUMVERTICES]; - - /* Get node coordinates and dof list: */ - this->GetVerticesCoordinates(&xyz_list); - - /*Retrieve all inputs we will need*/ - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /* Start looping on the number of vertices: */ - gauss=new GaussTria(); - for (int iv=0;ivGaussVertex(iv); - - /* Get the value we need*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vel=vx*vx+vy*vy; - - /*Compute strain rate viscosity and pressure: */ - this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - strainxx=epsilon[0]; - strainyy=epsilon[1]; - strainxy=epsilon[2]; - - /*strainparallel= Strain rate along the ice flow direction */ - strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6); - } - - /*Add input*/ - this->inputs->AddInput(new TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum)); - - /*Clean up and return*/ - delete gauss; - xDelete(xyz_list); -} -/*}}}*/ -void Tria::StrainRateperpendicular(){/*{{{*/ - - IssmDouble *xyz_list = NULL; - GaussTria* gauss=NULL; - IssmDouble epsilon[3]; - IssmDouble vx,vy,vel; - IssmDouble strainxx; - IssmDouble strainxy; - IssmDouble strainyy; - IssmDouble strainperpendicular[NUMVERTICES]; - - /* Get node coordinates and dof list: */ - this->GetVerticesCoordinates(&xyz_list); - - /*Retrieve all inputs we will need*/ - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - - /* Start looping on the number of vertices: */ - gauss=new GaussTria(); - for (int iv=0;ivGaussVertex(iv); - - /* Get the value we need*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vel=vx*vx+vy*vy; - - /*Compute strain rate viscosity and pressure: */ - this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - strainxx=epsilon[0]; - strainyy=epsilon[1]; - strainxy=epsilon[2]; - - /*strainperpendicular= Strain rate perpendicular to the ice flow direction */ - strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6); - } - - /*Add input*/ - this->inputs->AddInput(new TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum)); - - /*Clean up and return*/ - delete gauss; - xDelete(xyz_list); -} -/*}}}*/ -void Tria::CalvingRateLevermann(){/*{{{*/ - - IssmDouble xyz_list[NUMVERTICES][3]; - GaussTria* gauss=NULL; - IssmDouble vx,vy,vel; - IssmDouble strainparallel; - IssmDouble propcoeff; - IssmDouble strainperpendicular; - IssmDouble calvingratex[NUMVERTICES]; - IssmDouble calvingratey[NUMVERTICES]; - IssmDouble calvingrate[NUMVERTICES]; - - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*Retrieve all inputs and parameters we will need*/ - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum); _assert_(strainparallel_input); - Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input); - Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum); _assert_(levermanncoeff_input); - - /* Start looping on the number of vertices: */ - gauss=new GaussTria(); - for (int iv=0;ivGaussVertex(iv); - - /* Get the value we need*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vel=vx*vx+vy*vy; - strainparallel_input->GetInputValue(&strainparallel,gauss); - strainperpendicular_input->GetInputValue(&strainperpendicular,gauss); - levermanncoeff_input->GetInputValue(&propcoeff,gauss); - - /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */ - calvingrate[iv]=propcoeff*strainparallel*strainperpendicular; - if(calvingrate[iv]<0){ - calvingrate[iv]=0; - } - calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6); - calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6); - } - - /*Add input*/ - this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum)); - this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum)); - this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum)); - - /*Clean up and return*/ - delete gauss; - -} -/*}}}*/ void Tria::Configure(Elements* elementsin, Loads* loadsin,Nodes* nodesin,Vertices *verticesin,Materials* materialsin, Parameters* parametersin){/*{{{*/ @@ -508,4 +454,52 @@ } /*}}}*/ +void Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/ + + int vertexpidlist[NUMVERTICES]; + IssmDouble grad_list[NUMVERTICES]; + Input* grad_input=NULL; + + Input* input=inputs->GetInput(enum_type); + if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found"); + if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput"); + + GradientIndexing(&vertexpidlist[0],control_index); + for(int i=0;iSetGradient(grad_input); + +}/*}}}*/ +void Tria::ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum){/*{{{*/ + + Input* input=inputs->GetInput(control_enum); + if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found"); + if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput"); + + int sidlist[NUMVERTICES]; + int connectivity[NUMVERTICES]; + IssmPDouble values[NUMVERTICES]; + IssmPDouble gradients[NUMVERTICES]; + IssmDouble value,gradient; + + this->GetVerticesConnectivityList(&connectivity[0]); + this->GetVerticesSidList(&sidlist[0]); + + GaussTria* gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + ((ControlInput*)input)->GetInputValue(&value,gauss); + ((ControlInput*)input)->GetGradientValue(&gradient,gauss); + + values[iv] = reCast(value)/reCast(connectivity[iv]); + gradients[iv] = reCast(gradient)/reCast(connectivity[iv]); + } + delete gauss; + + vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL); + vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL); + +}/*}}}*/ void Tria::Delta18oParameterization(void){/*{{{*/ @@ -578,4 +572,105 @@ } /*}}}*/ +int Tria::EdgeOnBaseIndex(void){/*{{{*/ + + IssmDouble values[NUMVERTICES]; + int indices[3][2] = {{1,2},{2,0},{0,1}}; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&values[0],MeshVertexonbaseEnum); + + for(int i=0;i<3;i++){ + if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){ + return i; + } + } + + _printf_("list of vertices on bed: "<material->GetBbar(); + break; + + case VelEnum:{ + + /*Get input:*/ + IssmDouble vel; + Input* vel_input; + + vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input); + vel_input->GetInputAverage(&vel); + + /*Assign output pointers:*/ + *presponse=vel;} + break; + default: + _error_("Response type " << EnumToStringx(response_enum) << " not supported yet!"); + } + +} +/*}}}*/ void Tria::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){/*{{{*/ @@ -605,8 +700,88 @@ } /*}}}*/ -int Tria::ObjectEnum(void){/*{{{*/ - - return TriaEnum; - +void Tria::FSContactMigration(Vector* vertexgrounded,Vector* vertexfloating){/*{{{*/ + + if(!IsOnBase()) return; + + int approximation; + inputs->GetInputValue(&approximation,ApproximationEnum); + + if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){ + for(int i=0;iSetValue(vertices[i]->Pid(),+9999.,INS_VAL); + vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL); + } + } + else{ + /*Intermediaries*/ + IssmDouble* xyz_list = NULL; + IssmDouble* xyz_list_base = NULL; + IssmDouble pressure,water_pressure,sigma_nn,viscosity,bed,base; + IssmDouble bed_normal[2]; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble surface=0,value=0; + bool grounded; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list); + GetVerticesCoordinatesBase(&xyz_list_base); + + /*Retrieve all inputs we will be needing: */ + Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input); + Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input); + Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); + + /*Create gauss point in the middle of the basal edge*/ + Gauss* gauss=NewGaussBase(1); + gauss->GaussPoint(0); + + if(!IsFloating()){ + /*Check for basal force only if grounded and touching GL*/ + // if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){ + this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL); + pressure_input->GetInputValue(&pressure, gauss); + base_input->GetInputValue(&base, gauss); + + /*Compute Stress*/ + IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure; + IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure; + IssmDouble sigma_xy=2.*viscosity*epsilon[2]; + + /*Get normal vector to the bed */ + NormalBase(&bed_normal[0],xyz_list_base); + + /*basalforce*/ + sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1]; + + /*Compute water pressure*/ + IssmDouble rho_ice = matpar->GetRhoIce(); + IssmDouble rho_water = matpar->GetRhoWater(); + IssmDouble gravity = matpar->GetG(); + water_pressure=gravity*rho_water*base; + + /*Compare basal stress to water pressure and determine whether it should ground*/ + if (sigma_nnGetInputValue(&base, gauss); + bed_input->GetInputValue(&bed, gauss); + if(baseSetValue(vertices[i]->Pid(),+1.,INS_VAL); + else vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL); + } + + /*clean up*/ + delete gauss; + xDelete(xyz_list); + xDelete(xyz_list_base); + } } /*}}}*/ @@ -666,54 +841,4 @@ return this->element_type; -} -/*}}}*/ -void Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/ - - /*Computeportion of the element that has a positive levelset*/ - - bool negative=true; - int point; - const IssmPDouble epsilon= 1.e-15; - IssmDouble f1,f2; - - /*Be sure that values are not zero*/ - if(gl[0]==0.) gl[0]=gl[0]+epsilon; - if(gl[1]==0.) gl[1]=gl[1]+epsilon; - if(gl[2]==0.) gl[2]=gl[2]+epsilon; - - /*Check that not all nodes are positive or negative*/ - if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive - point=0; - f1=1.; - f2=1.; - } - else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative - point=0; - f1=0.; - f2=0.; - } - else{ - if(gl[0]*gl[1]*gl[2]<0) negative=false; - - if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 - point=2; - f1=gl[2]/(gl[2]-gl[0]); - f2=gl[2]/(gl[2]-gl[1]); - } - else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 - point=0; - f1=gl[0]/(gl[0]-gl[1]); - f2=gl[0]/(gl[0]-gl[2]); - } - else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 - point=1; - f1=gl[1]/(gl[1]-gl[2]); - f2=gl[1]/(gl[1]-gl[0]); - } - } - *point1=point; - *fraction1=f1; - *fraction2=f2; - *mainlynegative=negative; } /*}}}*/ @@ -889,147 +1014,4 @@ } /*}}}*/ -void Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/ - - int indices[2]; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Element XYZ list*/ - ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); - - /*Allocate Output*/ - IssmDouble* xyz_list_edge = xNew(2*3); - this->EdgeOnBaseIndices(&indices[0],&indices[1]); - for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j]; - - /*Assign output pointer*/ - *pxyz_list = xyz_list_edge; - -}/*}}}*/ -void Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/ - - int indices[2]; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Element XYZ list*/ - ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); - - /*Allocate Output*/ - IssmDouble* xyz_list_edge = xNew(2*3); - this->EdgeOnSurfaceIndices(&indices[0],&indices[1]); - for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j]; - - /*Assign output pointer*/ - *pxyz_list = xyz_list_edge; - -}/*}}}*/ -void Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/ - - /*Build unit outward pointing vector*/ - IssmDouble vector[2]; - IssmDouble norm; - - vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0]; - vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1]; - - norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]); - - normal[0]= + vector[1]/norm; - normal[1]= - vector[0]/norm; -} -/*}}}*/ -void Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ - - int normal_orientation=0; - IssmDouble s1,s2; - IssmDouble levelset[NUMVERTICES]; - - /*Recover parameters and values*/ - IssmDouble* xyz_zero = xNew(2*3); - GetInputListOnVertices(&levelset[0],levelsetenum); - - if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 - /*Portion of the segments*/ - s1=levelset[2]/(levelset[2]-levelset[1]); - s2=levelset[2]/(levelset[2]-levelset[0]); - - if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle depending on distribution of levelsetfunction - /*New point 1*/ - xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]); - xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]); - xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]); - - /*New point 0*/ - xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]); - xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]); - xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]); - } - else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 - /*Portion of the segments*/ - s1=levelset[0]/(levelset[0]-levelset[2]); - s2=levelset[0]/(levelset[0]-levelset[1]); - - if(levelset[0]<0.) normal_orientation=1; - /*New point 1*/ - xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]); - xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]); - xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]); - - /*New point 2*/ - xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]); - xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]); - xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]); - } - else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 - /*Portion of the segments*/ - s1=levelset[1]/(levelset[1]-levelset[0]); - s2=levelset[1]/(levelset[1]-levelset[2]); - - if(levelset[1]<0.) normal_orientation=1; - /*New point 0*/ - xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]); - xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]); - xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]); - - /*New point 2*/ - xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]); - xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]); - xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]); - } - else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1 - xyz_zero[3*0+0]=xyz_list[0*3+0]; - xyz_zero[3*0+1]=xyz_list[0*3+1]; - xyz_zero[3*0+2]=xyz_list[0*3+2]; - - /*New point 2*/ - xyz_zero[3*1+0]=xyz_list[1*3+0]; - xyz_zero[3*1+1]=xyz_list[1*3+1]; - xyz_zero[3*1+2]=xyz_list[1*3+2]; - } - else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1 - xyz_zero[3*0+0]=xyz_list[2*3+0]; - xyz_zero[3*0+1]=xyz_list[2*3+1]; - xyz_zero[3*0+2]=xyz_list[2*3+2]; - - /*New point 2*/ - xyz_zero[3*1+0]=xyz_list[0*3+0]; - xyz_zero[3*1+1]=xyz_list[0*3+1]; - xyz_zero[3*1+2]=xyz_list[0*3+2]; - } - else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1 - xyz_zero[3*0+0]=xyz_list[1*3+0]; - xyz_zero[3*0+1]=xyz_list[1*3+1]; - xyz_zero[3*0+2]=xyz_list[1*3+2]; - - /*New point 2*/ - xyz_zero[3*1+0]=xyz_list[2*3+0]; - xyz_zero[3*1+1]=xyz_list[2*3+1]; - xyz_zero[3*1+2]=xyz_list[2*3+2]; - } - else _error_("Case not covered"); - - /*Assign output pointer*/ - *pxyz_zero= xyz_zero; -} -/*}}}*/ void Tria::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ @@ -1072,4 +1054,16 @@ xDelete(indicesfront); }/*}}}*/ +void Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/ + + Input* input=inputs->GetInput(enumtype); + if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria"); + + GaussTria* gauss=new GaussTria(); + gauss->GaussVertex(this->GetNodeIndex(node)); + + input->GetInputValue(pvalue,gauss); + delete gauss; +} +/*}}}*/ void Tria::GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){/*{{{*/ @@ -1112,4 +1106,60 @@ xDelete(indicesfront); }/*}}}*/ +void Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/ + + /*Computeportion of the element that has a positive levelset*/ + + bool negative=true; + int point; + const IssmPDouble epsilon= 1.e-15; + IssmDouble f1,f2; + + /*Be sure that values are not zero*/ + if(gl[0]==0.) gl[0]=gl[0]+epsilon; + if(gl[1]==0.) gl[1]=gl[1]+epsilon; + if(gl[2]==0.) gl[2]=gl[2]+epsilon; + + /*Check that not all nodes are positive or negative*/ + if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive + point=0; + f1=1.; + f2=1.; + } + else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative + point=0; + f1=0.; + f2=0.; + } + else{ + if(gl[0]*gl[1]*gl[2]<0) negative=false; + + if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + point=2; + f1=gl[2]/(gl[2]-gl[0]); + f2=gl[2]/(gl[2]-gl[1]); + } + else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 + point=0; + f1=gl[0]/(gl[0]-gl[1]); + f2=gl[0]/(gl[0]-gl[2]); + } + else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 + point=1; + f1=gl[1]/(gl[1]-gl[2]); + f2=gl[1]/(gl[1]-gl[0]); + } + } + *point1=point; + *fraction1=f1; + *fraction2=f2; + *mainlynegative=negative; +} +/*}}}*/ +Node* Tria::GetNode(int node_number){/*{{{*/ + _assert_(node_number>=0); + _assert_(node_numberNumberofNodes(this->element_type)); + return this->nodes[node_number]; + +}/*}}}*/ int Tria::GetNodeIndex(Node* node){/*{{{*/ @@ -1135,22 +1185,219 @@ } /*}}}*/ -void Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/ - - Input* input=inputs->GetInput(enumtype); - if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria"); - +void Tria::GetSolutionFromInputsOneDof(Vector* solution, int enum_type){/*{{{*/ + + int *doflist = NULL; + IssmDouble value; + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(this->element_type); + + /*Fetch dof list and allocate solution vector*/ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numnodes); + + /*Get inputs*/ + Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input); + + /*Ok, we have the values, fill in the array: */ GaussTria* gauss=new GaussTria(); - gauss->GaussVertex(this->GetNodeIndex(node)); - - input->GetInputValue(pvalue,gauss); + for(int i=0;iGaussNode(this->element_type,i); + + enum_input->GetInputValue(&value,gauss); + values[i]=value; + } + + solution->SetValues(numnodes,doflist,values,INS_VAL); + + /*Free ressources:*/ + xDelete(doflist); + xDelete(values); delete gauss; } /*}}}*/ -Node* Tria::GetNode(int node_number){/*{{{*/ - _assert_(node_number>=0); - _assert_(node_numberNumberofNodes(this->element_type)); - return this->nodes[node_number]; +void Tria::GetVectorFromControlInputs(Vector* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/ + + int vertexidlist[NUMVERTICES]; + Input *input=NULL; + + /*Get out if this is not an element input*/ + if(!IsInput(control_enum)) return; + + /*Prepare index list*/ + GradientIndexing(&vertexidlist[0],control_index,onsid); + + /*Get input (either in element or material)*/ + input=(Input*)this->inputs->GetInput(control_enum); _assert_(input); + + /*Check that it is a ControlInput*/ + if (input->ObjectEnum()!=ControlInputEnum){ + _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); + } + + ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data); +} +/*}}}*/ +void Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/ + + int indices[2]; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Element XYZ list*/ + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + + /*Allocate Output*/ + IssmDouble* xyz_list_edge = xNew(2*3); + this->EdgeOnBaseIndices(&indices[0],&indices[1]); + for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j]; + + /*Assign output pointer*/ + *pxyz_list = xyz_list_edge; }/*}}}*/ +void Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/ + + int indices[2]; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Element XYZ list*/ + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + + /*Allocate Output*/ + IssmDouble* xyz_list_edge = xNew(2*3); + this->EdgeOnSurfaceIndices(&indices[0],&indices[1]); + for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j]; + + /*Assign output pointer*/ + *pxyz_list = xyz_list_edge; + +}/*}}}*/ +bool Tria::HasEdgeOnBase(){/*{{{*/ + + IssmDouble values[NUMVERTICES]; + IssmDouble sum; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&values[0],MeshVertexonbaseEnum); + sum = values[0]+values[1]+values[2]; + + _assert_(sum==0. || sum==1. || sum==2.); + + if(sum==3.) _error_("Two edges on bed not supported yet..."); + + if(sum>1.){ + return true; + } + else{ + return false; + } +} +/*}}}*/ +bool Tria::HasEdgeOnSurface(){/*{{{*/ + + IssmDouble values[NUMVERTICES]; + IssmDouble sum; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum); + sum = values[0]+values[1]+values[2]; + + _assert_(sum==0. || sum==1. || sum==2.); + + if(sum==3.) _error_("Two edges on surface not supported yet..."); + + if(sum>1.){ + return true; + } + else{ + return false; + } +} +/*}}}*/ +IssmDouble Tria::IceVolume(void){/*{{{*/ + + /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/ + IssmDouble base,surface,bed; + IssmDouble xyz_list[NUMVERTICES][3]; + + if(!IsIceInElement())return 0; + + /*First get back the area of the base*/ + base=this->GetArea(); + + /*Now get the average height*/ + Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); + Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input); + surface_input->GetInputAverage(&surface); + base_input->GetInputAverage(&bed); + + /*Return: */ + int domaintype; + parameters->FindParam(&domaintype,DomainTypeEnum); + if(domaintype==Domain2DverticalEnum){ + return base; + } + else{ + return base*(surface-bed); + } +} +/*}}}*/ +IssmDouble Tria::IceVolumeAboveFloatation(void){/*{{{*/ + + /*The volume above floatation: H + rho_water/rho_ice * bathymetry */ + IssmDouble rho_ice,rho_water; + IssmDouble base,surface,bed,bathymetry; + IssmDouble xyz_list[NUMVERTICES][3]; + + if(!IsIceInElement() || IsFloating())return 0; + + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*First calculate the area of the base (cross section triangle) + * http://en.wikipedia.org/wiki/Triangle + * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ + base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); + + /*Now get the average height and bathymetry*/ + Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); + Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input); + Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); + surface_input->GetInputAverage(&surface); + base_input->GetInputAverage(&bed); + bed_input->GetInputAverage(&bathymetry); + + /*Return: */ + return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.)); +} +/*}}}*/ +void Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/ + + /*Intermediary*/ + int num_controls; + int* control_type=NULL; + Input* input=NULL; + + /*retrieve some parameters: */ + this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + + for(int i=0;iinputs->GetInput(control_type[i]); _assert_(input); + if (input->ObjectEnum()!=ControlInputEnum){ + _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput"); + } + + ((ControlInput*)input)->UpdateValue(scalar); + ((ControlInput*)input)->Constrain(); + if (save_parameter) ((ControlInput*)input)->SaveValue(); + + } + + /*Clean up and return*/ + xDelete(control_type); +} +/*}}}*/ void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/ @@ -1371,4 +1618,57 @@ } /*}}}*/ +bool Tria::IsFaceOnBoundary(void){/*{{{*/ + + IssmDouble values[NUMVERTICES]; + IssmDouble sum; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum); + sum = values[0]+values[1]+values[2]; + + _assert_(sum==0. || sum==1. || sum==2.); + + if(sum==3.) _error_("Two edges on boundary not supported yet..."); + + if(sum>1.){ + return true; + } + else{ + return false; + } +}/*}}}*/ +bool Tria::IsIcefront(void){/*{{{*/ + + bool isicefront; + int i,nrice; + IssmDouble ls[NUMVERTICES]; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); + + /* If only one vertex has ice, there is an ice front here */ + isicefront=false; + if(IsIceInElement()){ + nrice=0; + for(i=0;iPid()]<0.){ + shelf=true; + break; + } + } + return shelf; +} +/*}}}*/ bool Tria::IsOnBase(){/*{{{*/ @@ -1397,4 +1697,23 @@ } /*}}}*/ +bool Tria::IsZeroLevelset(int levelset_enum){/*{{{*/ + + bool iszerols; + IssmDouble ls[NUMVERTICES]; + + /*Retrieve all inputs and parameters*/ + GetInputListOnVertices(&ls[0],levelset_enum); + + /*If the level set is awlays <0, there is no ice front here*/ + iszerols= false; + if(IsIceInElement()){ + if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){ + iszerols = true; + } + } + + return iszerols; +} +/*}}}*/ void Tria::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ @@ -1425,220 +1744,265 @@ } /*}}}*/ -bool Tria::HasEdgeOnBase(){/*{{{*/ - - IssmDouble values[NUMVERTICES]; - IssmDouble sum; - - /*Retrieve all inputs and parameters*/ - GetInputListOnVertices(&values[0],MeshVertexonbaseEnum); - sum = values[0]+values[1]+values[2]; - - _assert_(sum==0. || sum==1. || sum==2.); - - if(sum==3.) _error_("Two edges on bed not supported yet..."); - - if(sum>1.){ - return true; +IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/ + + + /*intermediary: */ + IssmDouble* values=NULL; + Input* thickness_input=NULL; + IssmDouble thickness; + IssmDouble weight; + IssmDouble Jdet; + IssmDouble volume; + IssmDouble rho_ice; + IssmDouble* xyz_list=NULL; + int point1; + IssmDouble fraction1,fraction2; + bool mainlynegative=true; + + /*Output:*/ + volume=0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list); + + /*Retrieve inputs required:*/ + thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); + + /*Retrieve material parameters: */ + rho_ice=matpar->GetRhoIce(); + + /*Retrieve values of the levelset defining the masscon: */ + values = xNew(NUMVERTICES); + for(int i=0;ivertices[i]->Sid()]; + } + + /*Ok, use the level set values to figure out where we put our gaussian points:*/ + this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values); + Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4); + + volume=0; + + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + this->JacobianDeterminant(&Jdet,xyz_list,gauss); + thickness_input->GetInputValue(&thickness, gauss); + + volume+=thickness*gauss->weight*Jdet; + } + + /* clean up and Return: */ + xDelete(xyz_list); + xDelete(values); + delete gauss; + return rho_ice*volume; +} +/*}}}*/ +IssmDouble Tria::MassFlux(IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/ + + int domaintype; + IssmDouble mass_flux=0.; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + IssmDouble length,rho_ice; + IssmDouble h1,h2; + IssmDouble vx1,vx2,vy1,vy2; + GaussTria* gauss_1=NULL; + GaussTria* gauss_2=NULL; + + /*Get material parameters :*/ + rho_ice=matpar->GetRhoIce(); + + /*First off, check that this segment belongs to this element: */ + if (segment_id!=this->id)_error_("error message: segment with id " << segment_id << " does not belong to element with id:" << this->id); + + /*Get xyz list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*get area coordinates of 0 and 1 locations: */ + gauss_1=new GaussTria(); + gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); + gauss_2=new GaussTria(); + gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); + + 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)); + + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + this->parameters->FindParam(&domaintype,DomainTypeEnum); + Input* vx_input=NULL; + Input* vy_input=NULL; + if(domaintype==Domain2DhorizontalEnum){ + vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); } else{ - return false; - } -} -/*}}}*/ -bool Tria::HasEdgeOnSurface(){/*{{{*/ - - IssmDouble values[NUMVERTICES]; - IssmDouble sum; - - /*Retrieve all inputs and parameters*/ - GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum); - sum = values[0]+values[1]+values[2]; - - _assert_(sum==0. || sum==1. || sum==2.); - - if(sum==3.) _error_("Two edges on surface not supported yet..."); - - if(sum>1.){ - return true; + vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input); + vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input); + } + + thickness_input->GetInputValue(&h1, gauss_1); + thickness_input->GetInputValue(&h2, gauss_2); + vx_input->GetInputValue(&vx1,gauss_1); + vx_input->GetInputValue(&vx2,gauss_2); + vy_input->GetInputValue(&vy1,gauss_1); + vy_input->GetInputValue(&vy2,gauss_2); + + mass_flux= rho_ice*length*( + (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+ + (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1] + ); + + /*clean up and return:*/ + delete gauss_1; + delete gauss_2; + return mass_flux; +} +/*}}}*/ +IssmDouble Tria::MassFlux(IssmDouble* segment){/*{{{*/ + + int domaintype; + IssmDouble mass_flux=0.; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + IssmDouble length,rho_ice; + IssmDouble x1,y1,x2,y2,h1,h2; + IssmDouble vx1,vx2,vy1,vy2; + GaussTria* gauss_1=NULL; + GaussTria* gauss_2=NULL; + + /*Get material parameters :*/ + rho_ice=matpar->GetRhoIce(); + + /*First off, check that this segment belongs to this element: */ + if (reCast(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast(*(segment+4)) << " does not belong to element with id:" << this->id); + + /*Recover segment node locations: */ + x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3); + + /*Get xyz list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*get area coordinates of 0 and 1 locations: */ + gauss_1=new GaussTria(); + gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); + gauss_2=new GaussTria(); + gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); + + 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)); + + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + this->parameters->FindParam(&domaintype,DomainTypeEnum); + Input* vx_input=NULL; + Input* vy_input=NULL; + if(domaintype==Domain2DhorizontalEnum){ + vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); } else{ - return false; - } -} -/*}}}*/ -void Tria::EdgeOnBaseIndices(int* pindex1,int* pindex2){/*{{{*/ - - IssmDouble values[NUMVERTICES]; - int indices[3][2] = {{1,2},{2,0},{0,1}}; - - /*Retrieve all inputs and parameters*/ - GetInputListOnVertices(&values[0],MeshVertexonbaseEnum); - - for(int i=0;i<3;i++){ - if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){ - *pindex1 = indices[i][0]; - *pindex2 = indices[i][1]; - return; - } - } - - _printf_("list of vertices on bed: "<* vertexgrounded,Vector* vertexfloating){/*{{{*/ - - if(!IsOnBase()) return; - - int approximation; - inputs->GetInputValue(&approximation,ApproximationEnum); - - if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){ - for(int i=0;iSetValue(vertices[i]->Pid(),+9999.,INS_VAL); - vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL); - } - } - else{ - /*Intermediaries*/ - IssmDouble* xyz_list = NULL; - IssmDouble* xyz_list_base = NULL; - IssmDouble pressure,water_pressure,sigma_nn,viscosity,bed,base; - IssmDouble bed_normal[2]; - IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - IssmDouble surface=0,value=0; - bool grounded; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list); - GetVerticesCoordinatesBase(&xyz_list_base); - - /*Retrieve all inputs we will be needing: */ - Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input); - Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input); - Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); - Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input); - - /*Create gauss point in the middle of the basal edge*/ - Gauss* gauss=NewGaussBase(1); - gauss->GaussPoint(0); - - if(!IsFloating()){ - /*Check for basal force only if grounded and touching GL*/ - // if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){ - this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL); - pressure_input->GetInputValue(&pressure, gauss); - base_input->GetInputValue(&base, gauss); - - /*Compute Stress*/ - IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure; - IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure; - IssmDouble sigma_xy=2.*viscosity*epsilon[2]; - - /*Get normal vector to the bed */ - NormalBase(&bed_normal[0],xyz_list_base); - - /*basalforce*/ - sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1]; - - /*Compute water pressure*/ - IssmDouble rho_ice = matpar->GetRhoIce(); - IssmDouble rho_water = matpar->GetRhoWater(); - IssmDouble gravity = matpar->GetG(); - water_pressure=gravity*rho_water*base; - - /*Compare basal stress to water pressure and determine whether it should ground*/ - if (sigma_nnGetInputValue(&base, gauss); - bed_input->GetInputValue(&bed, gauss); - if(baseSetValue(vertices[i]->Pid(),+1.,INS_VAL); - else vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL); - } - - /*clean up*/ - delete gauss; - xDelete(xyz_list); - xDelete(xyz_list_base); - } -} -/*}}}*/ -bool Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/ - - int i; - bool shelf=false; - - for(i=0;iPid()]<0.){ - shelf=true; - break; - } - } - return shelf; + vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input); + vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input); + } + + thickness_input->GetInputValue(&h1, gauss_1); + thickness_input->GetInputValue(&h2, gauss_2); + vx_input->GetInputValue(&vx1,gauss_1); + vx_input->GetInputValue(&vx2,gauss_2); + vy_input->GetInputValue(&vy1,gauss_1); + vy_input->GetInputValue(&vy2,gauss_2); + + mass_flux= rho_ice*length*( + (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+ + (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1] + ); + + /*clean up and return:*/ + delete gauss_1; + delete gauss_2; + return mass_flux; +} +/*}}}*/ +IssmDouble Tria::Misfit(int modelenum,int observationenum,int weightsenum){/*{{{*/ + + /*Intermediaries*/ + IssmDouble model,observation,weight; + IssmDouble Jdet; + IssmDouble Jelem = 0; + IssmDouble xyz_list[NUMVERTICES][3]; + GaussTria *gauss = NULL; + + /*If on water, return 0: */ + if(!IsIceInElement())return 0; + + /*Retrieve all inputs we will be needing: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* model_input=inputs->GetInput(modelenum); _assert_(model_input); + Input* observation_input=inputs->GetInput(observationenum);_assert_(observation_input); + Input* weights_input =inputs->GetInput(weightsenum); _assert_(weights_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Get parameters at gauss point*/ + model_input->GetInputValue(&model,gauss); + observation_input->GetInputValue(&observation,gauss); + weights_input->GetInputValue(&weight,gauss); + + /*compute misfit between model and observation */ + Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->weight; + } + + /* clean up and Return: */ + delete gauss; + return Jelem; +} +/*}}}*/ +IssmDouble Tria::MisfitArea(int weightsenum){/*{{{*/ + + /*Intermediaries*/ + IssmDouble weight; + IssmDouble Jdet; + IssmDouble Jelem = 0; + IssmDouble xyz_list[NUMVERTICES][3]; + GaussTria *gauss = NULL; + + /*If on water, return 0: */ + if(!IsIceInElement())return 0; + + /*Retrieve all inputs we will be needing: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* weights_input =inputs->GetInput(weightsenum); _assert_(weights_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Get parameters at gauss point*/ + weights_input->GetInputValue(&weight,gauss); + + /*compute misfit between model and observation */ + Jelem+=Jdet*weight*gauss->weight; + } + + /* clean up and Return: */ + delete gauss; + return Jelem; } /*}}}*/ @@ -1691,4 +2055,25 @@ } /*}}}*/ +void Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type); + +} +/*}}}*/ +void Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->VelocityInterpolation()); + +} +/*}}}*/ +void Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetNodalFunctions(basis,(GaussTria*)gauss,this->PressureInterpolation()); + +} +/*}}}*/ void Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/ @@ -1698,4 +2083,11 @@ } /*}}}*/ +void Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + + _assert_(gauss->Enum()==GaussTriaEnum); + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum); + +} +/*}}}*/ void Tria::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/ @@ -1705,22 +2097,8 @@ } /*}}}*/ -void Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ +void Tria::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/ _assert_(gauss->Enum()==GaussTriaEnum); - this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type); - -} -/*}}}*/ -void Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - - _assert_(gauss->Enum()==GaussTriaEnum); - this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum); - -} -/*}}}*/ -void Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - - _assert_(gauss->Enum()==GaussTriaEnum); - this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->VelocityInterpolation()); + this->GetNodalFunctions(basis,(GaussTria*)gauss,this->TensorInterpolation()); } @@ -1730,18 +2108,4 @@ _assert_(gauss->Enum()==GaussTriaEnum); this->GetNodalFunctions(basis,(GaussTria*)gauss,this->VelocityInterpolation()); - -} -/*}}}*/ -void Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/ - - _assert_(gauss->Enum()==GaussTriaEnum); - this->GetNodalFunctions(basis,(GaussTria*)gauss,this->PressureInterpolation()); - -} -/*}}}*/ -void Tria::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/ - - _assert_(gauss->Enum()==GaussTriaEnum); - this->GetNodalFunctions(basis,(GaussTria*)gauss,this->TensorInterpolation()); } @@ -1795,4 +2159,19 @@ } /*}}}*/ +void Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/ + + /*Build unit outward pointing vector*/ + IssmDouble vector[2]; + IssmDouble norm; + + vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0]; + vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1]; + + norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]); + + normal[0]= + vector[1]/norm; + normal[1]= - vector[0]/norm; +} +/*}}}*/ void Tria::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){/*{{{*/ @@ -1813,14 +2192,12 @@ } /*}}}*/ -int Tria::VelocityInterpolation(void){/*{{{*/ - return TriaRef::VelocityInterpolation(this->element_type); +int Tria::ObjectEnum(void){/*{{{*/ + + return TriaEnum; + } /*}}}*/ int Tria::PressureInterpolation(void){/*{{{*/ return TriaRef::PressureInterpolation(this->element_type); -} -/*}}}*/ -int Tria::TensorInterpolation(void){/*{{{*/ - return TriaRef::TensorInterpolation(this->element_type); } /*}}}*/ @@ -1901,4 +2278,31 @@ } /*}}}*/ +void Tria::PotentialUngrounding(Vector* potential_ungrounding){/*{{{*/ + + IssmDouble h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES]; + IssmDouble bed_hydro; + IssmDouble rho_water,rho_ice,density; + + /*material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + density=rho_ice/rho_water; + GetInputListOnVertices(&h[0],ThicknessEnum); + GetInputListOnVertices(&r[0],BedEnum); + GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum); + + /*go through vertices, and figure out which ones are grounded and want to unground: */ + for(int i=0;i0.){ + bed_hydro=-density*h[i]; + if(bed_hydro>r[i]){ + /*Vertex that could potentially unground, flag it*/ + potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL); + } + } + } +} +/*}}}*/ void Tria::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){/*{{{*/ @@ -2011,29 +2415,60 @@ } /*}}}*/ -Seg* Tria::SpawnSeg(int index1,int index2){/*{{{*/ - +void Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/ + + IssmDouble values[NUMVERTICES]; + int vertexpidlist[NUMVERTICES],control_init; + + + /*Get Domain type*/ + int domaintype; + parameters->FindParam(&domaintype,DomainTypeEnum); + + /*Specific case for depth averaged quantities*/ + control_init=control_enum; + if(domaintype==Domain2DverticalEnum){ + if(control_enum==MaterialsRheologyBbarEnum){ + control_enum=MaterialsRheologyBEnum; + if(!IsOnBase()) return; + } + if(control_enum==DamageDbarEnum){ + control_enum=DamageDEnum; + if(!IsOnBase()) return; + } + } + + /*Get out if this is not an element input*/ + if(!IsInput(control_enum)) return; + + /*Prepare index list*/ + GradientIndexing(&vertexpidlist[0],control_index); + + /*Get values on vertices*/ + for(int i=0;iinputs->GetInput(control_enum); _assert_(input); + if(input->ObjectEnum()!=ControlInputEnum){ + _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); + } + + ((ControlInput*)input)->SetInput(new_input); +} +/*}}}*/ +void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/ + + /*go into parameters and get the analysis_counter: */ int analysis_counter; - - /*go into parameters and get the analysis_counter: */ - this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum); - - /*Create Seg*/ - Seg* seg=new Seg(); - seg->id=this->id; - seg->inputs=(Inputs*)this->inputs->SpawnSegInputs(index1,index2); - seg->parameters=this->parameters; - seg->element_type=P1Enum; //Only P1 CG for now (TO BE CHANGED) - this->SpawnSegHook(xDynamicCast(seg),index1,index2); - - /*Spawn material*/ - seg->material=(Material*)this->material->copy2(seg); - - /*recover nodes, material and matpar: */ - seg->nodes = (Node**)seg->hnodes[analysis_counter]->deliverp(); - seg->vertices = (Vertex**)seg->hvertices->deliverp(); - seg->matpar = (Matpar*)seg->hmatpar->delivers(); - - /*Return new Seg*/ - return seg; + parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); + + /*Get Element type*/ + if(this->element_type_list) this->element_type=this->element_type_list[analysis_counter]; + + /*Pick up nodes*/ + if(this->hnodes && this->hnodes[analysis_counter]){ + this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); + } + } /*}}}*/ @@ -2056,4 +2491,31 @@ } /*}}}*/ +Seg* Tria::SpawnSeg(int index1,int index2){/*{{{*/ + + int analysis_counter; + + /*go into parameters and get the analysis_counter: */ + this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum); + + /*Create Seg*/ + Seg* seg=new Seg(); + seg->id=this->id; + seg->inputs=(Inputs*)this->inputs->SpawnSegInputs(index1,index2); + seg->parameters=this->parameters; + seg->element_type=P1Enum; //Only P1 CG for now (TO BE CHANGED) + this->SpawnSegHook(xDynamicCast(seg),index1,index2); + + /*Spawn material*/ + seg->material=(Material*)this->material->copy2(seg); + + /*recover nodes, material and matpar: */ + seg->nodes = (Node**)seg->hnodes[analysis_counter]->deliverp(); + seg->vertices = (Vertex**)seg->hvertices->deliverp(); + seg->matpar = (Matpar*)seg->hmatpar->delivers(); + + /*Return new Seg*/ + return seg; +} +/*}}}*/ Element* Tria::SpawnTopElement(void){/*{{{*/ @@ -2074,18 +2536,94 @@ } /*}}}*/ -void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/ - - /*go into parameters and get the analysis_counter: */ - int analysis_counter; - parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); - - /*Get Element type*/ - if(this->element_type_list) this->element_type=this->element_type_list[analysis_counter]; - - /*Pick up nodes*/ - if(this->hnodes && this->hnodes[analysis_counter]){ - this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); - } - +void Tria::StrainRateparallel(){/*{{{*/ + + IssmDouble *xyz_list = NULL; + IssmDouble epsilon[3]; + GaussTria* gauss=NULL; + IssmDouble vx,vy,vel; + IssmDouble strainxx; + IssmDouble strainxy; + IssmDouble strainyy; + IssmDouble strainparallel[NUMVERTICES]; + + /* Get node coordinates and dof list: */ + this->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will need*/ + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of vertices: */ + gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + /* Get the value we need*/ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel=vx*vx+vy*vy; + + /*Compute strain rate viscosity and pressure: */ + this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + strainxx=epsilon[0]; + strainyy=epsilon[1]; + strainxy=epsilon[2]; + + /*strainparallel= Strain rate along the ice flow direction */ + strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6); + } + + /*Add input*/ + this->inputs->AddInput(new TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); +} +/*}}}*/ +void Tria::StrainRateperpendicular(){/*{{{*/ + + IssmDouble *xyz_list = NULL; + GaussTria* gauss=NULL; + IssmDouble epsilon[3]; + IssmDouble vx,vy,vel; + IssmDouble strainxx; + IssmDouble strainxy; + IssmDouble strainyy; + IssmDouble strainperpendicular[NUMVERTICES]; + + /* Get node coordinates and dof list: */ + this->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will need*/ + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of vertices: */ + gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + /* Get the value we need*/ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel=vx*vx+vy*vy; + + /*Compute strain rate viscosity and pressure: */ + this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + strainxx=epsilon[0]; + strainyy=epsilon[1]; + strainxy=epsilon[2]; + + /*strainperpendicular= Strain rate perpendicular to the ice flow direction */ + strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6); + } + + /*Add input*/ + this->inputs->AddInput(new TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); } /*}}}*/ @@ -2115,4 +2653,8 @@ /*Return: */ return S; +} +/*}}}*/ +int Tria::TensorInterpolation(void){/*{{{*/ + return TriaRef::TensorInterpolation(this->element_type); } /*}}}*/ @@ -2156,4 +2698,32 @@ return dt; +} +/*}}}*/ +IssmDouble Tria::TotalSmb(void){/*{{{*/ + + /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/ + IssmDouble base,smb,rho_ice; + IssmDouble Total_Smb=0; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Get material parameters :*/ + rho_ice=matpar->GetRhoIce(); + + if(!IsIceInElement())return 0; + + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*First calculate the area of the base (cross section triangle) + * http://en.wikipedia.org/wiki/Triangle + * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ + base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2 + + /*Now get the average SMB over the element*/ + Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input); + smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1 + Total_Smb=rho_ice*base*smb; // smb on element in kg s-1 + + /*Return: */ + return Total_Smb; } /*}}}*/ @@ -2347,10 +2917,33 @@ } /*}}}*/ +int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/ + + int i; + int nflipped=0; + + /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */ + for(i=0;i<3;i++){ + if (reCast(vertices_potentially_ungrounding[vertices[i]->Pid()])){ + vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL); + + /*If node was not on ice shelf, we flipped*/ + if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){ + nflipped++; + } + } + } + return nflipped; +} +/*}}}*/ +void Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum); +} +/*}}}*/ void Tria::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/ TriaRef::GetInputValue(pvalue,values,gauss,P1Enum); } /*}}}*/ -void Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ - TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum); +int Tria::VelocityInterpolation(void){/*{{{*/ + return TriaRef::VelocityInterpolation(this->element_type); } /*}}}*/ @@ -2360,472 +2953,96 @@ } /*}}}*/ -bool Tria::IsZeroLevelset(int levelset_enum){/*{{{*/ - - bool iszerols; - IssmDouble ls[NUMVERTICES]; - - /*Retrieve all inputs and parameters*/ - GetInputListOnVertices(&ls[0],levelset_enum); - - /*If the level set is awlays <0, there is no ice front here*/ - iszerols= false; - if(IsIceInElement()){ - if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){ - iszerols = true; - } - } - - return iszerols; -} -/*}}}*/ -bool Tria::IsIcefront(void){/*{{{*/ - - bool isicefront; - int i,nrice; - IssmDouble ls[NUMVERTICES]; - - /*Retrieve all inputs and parameters*/ - GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); - - /* If only one vertex has ice, there is an ice front here */ - isicefront=false; - if(IsIceInElement()){ - nrice=0; - for(i=0;i1.){ - return true; - } - else{ - return false; - } -}/*}}}*/ -void Tria::AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/ - - bool already = false; - int i,j; - int partition[NUMVERTICES]; - int offsetsid[NUMVERTICES]; - int offsetdof[NUMVERTICES]; - IssmDouble area; - IssmDouble mean; - - /*First, get the area: */ - area=this->GetArea(); - - /*Figure out the average for this element: */ - this->GetVerticesSidList(&offsetsid[0]); - this->GetVertexPidList(&offsetdof[0]); - mean=0; - for(i=0;i(qmu_part[offsetsid[i]]); - mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]]; - } - - /*Add contribution: */ - for(i=0;iSetValue(partition[i],mean*area,ADD_VAL); - partition_areas->SetValue(partition[i],area,ADD_VAL); - }; - } -} -/*}}}*/ -IssmDouble Tria::IceVolume(void){/*{{{*/ - - /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/ - IssmDouble base,surface,bed; - IssmDouble xyz_list[NUMVERTICES][3]; - - if(!IsIceInElement())return 0; - - /*First get back the area of the base*/ - base=this->GetArea(); - - /*Now get the average height*/ - Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input); - surface_input->GetInputAverage(&surface); - base_input->GetInputAverage(&bed); - - /*Return: */ - int domaintype; - parameters->FindParam(&domaintype,DomainTypeEnum); - if(domaintype==Domain2DverticalEnum){ - return base; - } - else{ - return base*(surface-bed); - } -} -/*}}}*/ -IssmDouble Tria::IceVolumeAboveFloatation(void){/*{{{*/ - - /*The volume above floatation: H + rho_water/rho_ice * bathymetry */ - IssmDouble rho_ice,rho_water; - IssmDouble base,surface,bed,bathymetry; - IssmDouble xyz_list[NUMVERTICES][3]; - - if(!IsIceInElement() || IsFloating())return 0; - - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*First calculate the area of the base (cross section triangle) - * http://en.wikipedia.org/wiki/Triangle - * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ - base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); - - /*Now get the average height and bathymetry*/ - Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); - Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input); - Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); - surface_input->GetInputAverage(&surface); - base_input->GetInputAverage(&bed); - bed_input->GetInputAverage(&bathymetry); - - /*Return: */ - return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.)); -} -/*}}}*/ -IssmDouble Tria::MassFlux( IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/ - - int domaintype; - IssmDouble mass_flux=0.; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble length,rho_ice; - IssmDouble h1,h2; - IssmDouble vx1,vx2,vy1,vy2; - GaussTria* gauss_1=NULL; - GaussTria* gauss_2=NULL; - - /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); - - /*First off, check that this segment belongs to this element: */ - if (segment_id!=this->id)_error_("error message: segment with id " << segment_id << " does not belong to element with id:" << this->id); - - /*Get xyz list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*get area coordinates of 0 and 1 locations: */ - gauss_1=new GaussTria(); - gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); - gauss_2=new GaussTria(); - gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); - - 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)); - - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - this->parameters->FindParam(&domaintype,DomainTypeEnum); - Input* vx_input=NULL; - Input* vy_input=NULL; - if(domaintype==Domain2DhorizontalEnum){ - vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - } - else{ - vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input); - vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input); - } - - thickness_input->GetInputValue(&h1, gauss_1); - thickness_input->GetInputValue(&h2, gauss_2); - vx_input->GetInputValue(&vx1,gauss_1); - vx_input->GetInputValue(&vx2,gauss_2); - vy_input->GetInputValue(&vy1,gauss_1); - vy_input->GetInputValue(&vy2,gauss_2); - - mass_flux= rho_ice*length*( - (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+ - (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1] - ); - - /*clean up and return:*/ - delete gauss_1; - delete gauss_2; - return mass_flux; -} -/*}}}*/ -IssmDouble Tria::MassFlux( IssmDouble* segment){/*{{{*/ - - int domaintype; - IssmDouble mass_flux=0.; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble length,rho_ice; - IssmDouble x1,y1,x2,y2,h1,h2; - IssmDouble vx1,vx2,vy1,vy2; - GaussTria* gauss_1=NULL; - GaussTria* gauss_2=NULL; - - /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); - - /*First off, check that this segment belongs to this element: */ - if (reCast(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast(*(segment+4)) << " does not belong to element with id:" << this->id); - - /*Recover segment node locations: */ - x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3); - - /*Get xyz list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*get area coordinates of 0 and 1 locations: */ - gauss_1=new GaussTria(); - gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); - gauss_2=new GaussTria(); - gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); - - 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)); - - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - this->parameters->FindParam(&domaintype,DomainTypeEnum); - Input* vx_input=NULL; - Input* vy_input=NULL; - if(domaintype==Domain2DhorizontalEnum){ - vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - } - else{ - vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input); - vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input); - } - - thickness_input->GetInputValue(&h1, gauss_1); - thickness_input->GetInputValue(&h2, gauss_2); - vx_input->GetInputValue(&vx1,gauss_1); - vx_input->GetInputValue(&vx2,gauss_2); - vy_input->GetInputValue(&vy1,gauss_1); - vy_input->GetInputValue(&vy2,gauss_2); - - mass_flux= rho_ice*length*( - (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+ - (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1] - ); - - /*clean up and return:*/ - delete gauss_1; - delete gauss_2; - return mass_flux; -} -/*}}}*/ -void Tria::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/ - - switch(response_enum){ - case MaterialsRheologyBbarEnum: - *presponse=this->material->GetBbar(); - break; - - case VelEnum:{ - - /*Get input:*/ - IssmDouble vel; - Input* vel_input; - - vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input); - vel_input->GetInputAverage(&vel); - - /*Assign output pointers:*/ - *presponse=vel;} - break; - default: - _error_("Response type " << EnumToStringx(response_enum) << " not supported yet!"); - } - -} -/*}}}*/ -IssmDouble Tria::TotalSmb(void){/*{{{*/ - - /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/ - IssmDouble base,smb,rho_ice; - IssmDouble Total_Smb=0; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); - - if(!IsIceInElement())return 0; - - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*First calculate the area of the base (cross section triangle) - * http://en.wikipedia.org/wiki/Triangle - * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ - base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2 - - /*Now get the average SMB over the element*/ - Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input); - smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1 - Total_Smb=rho_ice*base*smb; // smb on element in kg s-1 - - /*Return: */ - return Total_Smb; -} -/*}}}*/ -IssmDouble Tria::MisfitArea(int weightsenum){/*{{{*/ - - /*Intermediaries*/ - IssmDouble weight; - IssmDouble Jdet; - IssmDouble Jelem = 0; - IssmDouble xyz_list[NUMVERTICES][3]; - GaussTria *gauss = NULL; - - /*If on water, return 0: */ - if(!IsIceInElement())return 0; - - /*Retrieve all inputs we will be needing: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* weights_input =inputs->GetInput(weightsenum); _assert_(weights_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get parameters at gauss point*/ - weights_input->GetInputValue(&weight,gauss); - - /*compute misfit between model and observation */ - Jelem+=Jdet*weight*gauss->weight; - } - - /* clean up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -IssmDouble Tria::Misfit(int modelenum,int observationenum,int weightsenum){/*{{{*/ - - /*Intermediaries*/ - IssmDouble model,observation,weight; - IssmDouble Jdet; - IssmDouble Jelem = 0; - IssmDouble xyz_list[NUMVERTICES][3]; - GaussTria *gauss = NULL; - - /*If on water, return 0: */ - if(!IsIceInElement())return 0; - - /*Retrieve all inputs we will be needing: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* model_input=inputs->GetInput(modelenum); _assert_(model_input); - Input* observation_input=inputs->GetInput(observationenum);_assert_(observation_input); - Input* weights_input =inputs->GetInput(weightsenum); _assert_(weights_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get parameters at gauss point*/ - model_input->GetInputValue(&model,gauss); - observation_input->GetInputValue(&observation,gauss); - weights_input->GetInputValue(&weight,gauss); - - /*compute misfit between model and observation */ - Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->weight; - } - - /* clean up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/ - - - /*intermediary: */ - IssmDouble* values=NULL; - Input* thickness_input=NULL; - IssmDouble thickness; - IssmDouble weight; - IssmDouble Jdet; - IssmDouble volume; - IssmDouble rho_ice; - IssmDouble* xyz_list=NULL; - int point1; - IssmDouble fraction1,fraction2; - bool mainlynegative=true; - - /*Output:*/ - volume=0; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list); - - /*Retrieve inputs required:*/ - thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); - - /*Retrieve material parameters: */ - rho_ice=matpar->GetRhoIce(); - - /*Retrieve values of the levelset defining the masscon: */ - values = xNew(NUMVERTICES); - for(int i=0;ivertices[i]->Sid()]; - } - - /*Ok, use the level set values to figure out where we put our gaussian points:*/ - this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values); - Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4); - - volume=0; - - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - this->JacobianDeterminant(&Jdet,xyz_list,gauss); - thickness_input->GetInputValue(&thickness, gauss); - - volume+=thickness*gauss->weight*Jdet; - } - - /* clean up and Return: */ - xDelete(xyz_list); - xDelete(values); - delete gauss; - return rho_ice*volume; +void Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ + + int normal_orientation=0; + IssmDouble s1,s2; + IssmDouble levelset[NUMVERTICES]; + + /*Recover parameters and values*/ + IssmDouble* xyz_zero = xNew(2*3); + GetInputListOnVertices(&levelset[0],levelsetenum); + + if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + /*Portion of the segments*/ + s1=levelset[2]/(levelset[2]-levelset[1]); + s2=levelset[2]/(levelset[2]-levelset[0]); + + if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle depending on distribution of levelsetfunction + /*New point 1*/ + xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]); + xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]); + xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]); + + /*New point 0*/ + xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]); + xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]); + xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]); + } + else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 + /*Portion of the segments*/ + s1=levelset[0]/(levelset[0]-levelset[2]); + s2=levelset[0]/(levelset[0]-levelset[1]); + + if(levelset[0]<0.) normal_orientation=1; + /*New point 1*/ + xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]); + xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]); + xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]); + + /*New point 2*/ + xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]); + xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]); + xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]); + } + else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 + /*Portion of the segments*/ + s1=levelset[1]/(levelset[1]-levelset[0]); + s2=levelset[1]/(levelset[1]-levelset[2]); + + if(levelset[1]<0.) normal_orientation=1; + /*New point 0*/ + xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]); + xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]); + xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]); + + /*New point 2*/ + xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]); + xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]); + xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]); + } + else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1 + xyz_zero[3*0+0]=xyz_list[0*3+0]; + xyz_zero[3*0+1]=xyz_list[0*3+1]; + xyz_zero[3*0+2]=xyz_list[0*3+2]; + + /*New point 2*/ + xyz_zero[3*1+0]=xyz_list[1*3+0]; + xyz_zero[3*1+1]=xyz_list[1*3+1]; + xyz_zero[3*1+2]=xyz_list[1*3+2]; + } + else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1 + xyz_zero[3*0+0]=xyz_list[2*3+0]; + xyz_zero[3*0+1]=xyz_list[2*3+1]; + xyz_zero[3*0+2]=xyz_list[2*3+2]; + + /*New point 2*/ + xyz_zero[3*1+0]=xyz_list[0*3+0]; + xyz_zero[3*1+1]=xyz_list[0*3+1]; + xyz_zero[3*1+2]=xyz_list[0*3+2]; + } + else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1 + xyz_zero[3*0+0]=xyz_list[1*3+0]; + xyz_zero[3*0+1]=xyz_list[1*3+1]; + xyz_zero[3*0+2]=xyz_list[1*3+2]; + + /*New point 2*/ + xyz_zero[3*1+0]=xyz_list[2*3+0]; + xyz_zero[3*1+1]=xyz_list[2*3+1]; + xyz_zero[3*1+2]=xyz_list[2*3+2]; + } + else _error_("Case not covered"); + + /*Assign output pointer*/ + *pxyz_zero= xyz_zero; } /*}}}*/ @@ -2959,177 +3176,43 @@ #endif -void Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/ - - /*Intermediary*/ - int num_controls; - int* control_type=NULL; - Input* input=NULL; - - /*retrieve some parameters: */ - this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); - this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); - - for(int i=0;iinputs->GetInput(control_type[i]); _assert_(input); - if (input->ObjectEnum()!=ControlInputEnum){ - _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput"); - } - - ((ControlInput*)input)->UpdateValue(scalar); - ((ControlInput*)input)->Constrain(); - if (save_parameter) ((ControlInput*)input)->SaveValue(); - - } - - /*Clean up and return*/ - xDelete(control_type); -} -/*}}}*/ -void Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/ - - int vertexpidlist[NUMVERTICES]; - IssmDouble grad_list[NUMVERTICES]; - Input* grad_input=NULL; - - Input* input=inputs->GetInput(enum_type); - if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found"); - if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput"); - - GradientIndexing(&vertexpidlist[0],control_index); - for(int i=0;iSetGradient(grad_input); - -}/*}}}*/ -void Tria::ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum){/*{{{*/ - - Input* input=inputs->GetInput(control_enum); - if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found"); - if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput"); - - int sidlist[NUMVERTICES]; - int connectivity[NUMVERTICES]; - IssmPDouble values[NUMVERTICES]; - IssmPDouble gradients[NUMVERTICES]; - IssmDouble value,gradient; - - this->GetVerticesConnectivityList(&connectivity[0]); - this->GetVerticesSidList(&sidlist[0]); - - GaussTria* gauss=new GaussTria(); - for (int iv=0;ivGaussVertex(iv); - - ((ControlInput*)input)->GetInputValue(&value,gauss); - ((ControlInput*)input)->GetGradientValue(&gradient,gauss); - - values[iv] = reCast(value)/reCast(connectivity[iv]); - gradients[iv] = reCast(gradient)/reCast(connectivity[iv]); - } - delete gauss; - - vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL); - vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL); - -}/*}}}*/ -void Tria::GetVectorFromControlInputs(Vector* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/ - - int vertexidlist[NUMVERTICES]; - Input *input=NULL; - - /*Get out if this is not an element input*/ - if(!IsInput(control_enum)) return; - - /*Prepare index list*/ - GradientIndexing(&vertexidlist[0],control_index,onsid); - - /*Get input (either in element or material)*/ - input=(Input*)this->inputs->GetInput(control_enum); _assert_(input); - - /*Check that it is a ControlInput*/ - if (input->ObjectEnum()!=ControlInputEnum){ - _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); - } - - ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data); -} -/*}}}*/ -void Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/ - - IssmDouble values[NUMVERTICES]; - int vertexpidlist[NUMVERTICES],control_init; - - - /*Get Domain type*/ - int domaintype; - parameters->FindParam(&domaintype,DomainTypeEnum); - - /*Specific case for depth averaged quantities*/ - control_init=control_enum; - if(domaintype==Domain2DverticalEnum){ - if(control_enum==MaterialsRheologyBbarEnum){ - control_enum=MaterialsRheologyBEnum; - if(!IsOnBase()) return; - } - if(control_enum==DamageDbarEnum){ - control_enum=DamageDEnum; - if(!IsOnBase()) return; - } - } - - /*Get out if this is not an element input*/ - if(!IsInput(control_enum)) return; - - /*Prepare index list*/ - GradientIndexing(&vertexpidlist[0],control_index); - - /*Get values on vertices*/ - for(int i=0;iinputs->GetInput(control_enum); _assert_(input); - if(input->ObjectEnum()!=ControlInputEnum){ - _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput"); - } - - ((ControlInput*)input)->SetInput(new_input); -} -/*}}}*/ -void Tria::GetSolutionFromInputsOneDof(Vector* solution, int enum_type){/*{{{*/ - - int *doflist = NULL; - IssmDouble value; - - /*Fetch number of nodes for this finite element*/ - int numnodes = this->NumberofNodes(this->element_type); - - /*Fetch dof list and allocate solution vector*/ - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - IssmDouble* values = xNew(numnodes); - - /*Get inputs*/ - Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input); - - /*Ok, we have the values, fill in the array: */ - GaussTria* gauss=new GaussTria(); - for(int i=0;iGaussNode(this->element_type,i); - - enum_input->GetInputValue(&value,gauss); - values[i]=value; - } - - solution->SetValues(numnodes,doflist,values,INS_VAL); - - /*Free ressources:*/ - xDelete(doflist); - xDelete(values); - delete gauss; -} -/*}}}*/ - #ifdef _HAVE_DAKOTA_ +void Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/ + + int i,t,row; + IssmDouble time; + TransientInput *transientinput = NULL; + IssmDouble values[3]; + + /*Check that name is an element input*/ + if (!IsInput(name)) return; + + switch(type){ + + case VertexEnum: + /*Create transient input: */ + for(t=0;tvertices[i]->Sid(); + values[i]=matrix[ncols*row+t]; + } + + /*time:*/ + time=matrix[(nrows-1)*ncols+t]; + + if(t==0) transientinput=new TransientInput(name); + transientinput->AddTimeInput(new TriaInput(name,values,P1Enum),time); + transientinput->Configure(parameters); + } + this->inputs->AddInput(transientinput); + break; + + default: + _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); + } + +} +/*}}}*/ void Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/ @@ -3223,88 +3306,3 @@ } /*}}}*/ -void Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/ - - int i,t,row; - IssmDouble time; - TransientInput *transientinput = NULL; - IssmDouble values[3]; - - /*Check that name is an element input*/ - if (!IsInput(name)) return; - - switch(type){ - - case VertexEnum: - /*Create transient input: */ - for(t=0;tvertices[i]->Sid(); - values[i]=matrix[ncols*row+t]; - } - - /*time:*/ - time=matrix[(nrows-1)*ncols+t]; - - if(t==0) transientinput=new TransientInput(name); - transientinput->AddTimeInput(new TriaInput(name,values,P1Enum),time); - transientinput->Configure(parameters); - } - this->inputs->AddInput(transientinput); - break; - - default: - _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); - } - -} -/*}}}*/ #endif - -void Tria::PotentialUngrounding(Vector* potential_ungrounding){/*{{{*/ - - IssmDouble h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES]; - IssmDouble bed_hydro; - IssmDouble rho_water,rho_ice,density; - - /*material parameters: */ - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - density=rho_ice/rho_water; - GetInputListOnVertices(&h[0],ThicknessEnum); - GetInputListOnVertices(&r[0],BedEnum); - GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum); - - /*go through vertices, and figure out which ones are grounded and want to unground: */ - for(int i=0;i0.){ - bed_hydro=-density*h[i]; - if(bed_hydro>r[i]){ - /*Vertex that could potentially unground, flag it*/ - potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL); - } - } - } -} -/*}}}*/ -int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/ - - int i; - int nflipped=0; - - /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */ - for(i=0;i<3;i++){ - if (reCast(vertices_potentially_ungrounding[vertices[i]->Pid()])){ - vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL); - - /*If node was not on ice shelf, we flipped*/ - if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){ - nflipped++; - } - } - } - return nflipped; -} -/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18910) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18911) @@ -42,103 +42,99 @@ /*}}}*/ /*Update virtual functions resolution: {{{*/ - void InputUpdateFromVector(IssmDouble* vector, int name, int type); #ifdef _HAVE_DAKOTA_ + void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type); void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type); - void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type); #endif void InputUpdateFromIoModel(int index, IoModel* iomodel); + void InputUpdateFromVector(IssmDouble* vector, int name, int type); /*}}}*/ /*Element virtual functions definitions: {{{*/ + void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); + void CalvingRateLevermann(); IssmDouble CharacteristicLength(void); void ComputeBasalStress(Vector* sigma_b); + void ComputeDeviatoricStressTensor(); void ComputeSigmaNN(); void ComputeStressTensor(); - void ComputeDeviatoricStressTensor(); - void StrainRateparallel(); - void StrainRateperpendicular(); void ComputeSurfaceNormalVelocity(); - void StressIntensityFactor(void){_error_("not implemented yet");}; - void CalvingRateLevermann(); void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); - void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); - void ResetHooks(); + void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); + void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); void Delta18oParameterization(void); + int EdgeOnBaseIndex(); + void EdgeOnBaseIndices(int* pindex1,int* pindex); + int EdgeOnSurfaceIndex(); + void EdgeOnSurfaceIndices(int* pindex1,int* pindex); + void ElementResponse(IssmDouble* presponse,int response_enum); void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); + int FiniteElement(void); void FSContactMigration(Vector* vertexgrounded,Vector* vertexfloating); - int FiniteElement(void); - Element* GetUpperElement(void){_error_("not implemented yet");}; Element* GetBasalElement(void){_error_("not implemented yet");}; void GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues); void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating); IssmDouble GetGroundedPortion(IssmDouble* xyz_list); + void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); + void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level); int GetNodeIndex(Node* node); int GetNumberOfNodes(void); int GetNumberOfNodes(int enum_type); int GetNumberOfVertices(void); - bool IsOnBase(); - bool IsOnSurface(); + void GetSolutionFromInputsOneDof(Vector* solution,int enum_type); + Element* GetUpperElement(void){_error_("not implemented yet");}; + void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,bool onsid); + void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); + void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); bool HasEdgeOnBase(); bool HasEdgeOnSurface(); - void EdgeOnSurfaceIndices(int* pindex1,int* pindex); - void EdgeOnBaseIndices(int* pindex1,int* pindex); - int EdgeOnBaseIndex(); - int EdgeOnSurfaceIndex(); - bool IsNodeOnShelfFromFlags(IssmDouble* flags); - int NumberofNodesVelocity(void); - int NumberofNodesPressure(void); - void GetSolutionFromInputsOneDof(Vector* solution,int enum_type); - void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); - void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); + IssmDouble IceVolume(void); + IssmDouble IceVolumeAboveFloatation(void); + void InputControlUpdate(IssmDouble scalar,bool save_parameter); void InputDepthAverageAtBase(int enum_type,int average_enum_type); void InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/}; void InputScale(int enum_type,IssmDouble scale_factor); + bool IsFaceOnBoundary(void); + bool IsIcefront(void); + bool IsNodeOnShelfFromFlags(IssmDouble* flags); + bool IsOnBase(); + bool IsOnSurface(); + bool IsZeroLevelset(int levelset_enum); + IssmDouble Masscon(IssmDouble* levelset); + IssmDouble MassFlux(IssmDouble* segment); + IssmDouble MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id); void MaterialUpdateFromTemperature(void){_error_("not implemented yet");}; + IssmDouble Misfit(int modelenum,int observationenum,int weightsenum); + IssmDouble MisfitArea(int weightsenum); int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum); + int NumberofNodesPressure(void); + int NumberofNodesVelocity(void); void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); + void PotentialUngrounding(Vector* potential_sheet_ungrounding); + int PressureInterpolation(); void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe); void ResetFSBasalBoundaryCondition(void); + void ResetHooks(); + void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index); + void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); Element* SpawnBasalElement(void); Element* SpawnTopElement(void); + void StrainRateparallel(); + void StrainRateperpendicular(); + void StressIntensityFactor(void){_error_("not implemented yet");}; + IssmDouble SurfaceArea(void); + int TensorInterpolation(); + IssmDouble TimeAdapt(); + IssmDouble TotalSmb(void); + void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement); + int UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf); + void ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss); + void ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss); int VelocityInterpolation(); - int PressureInterpolation(); - int TensorInterpolation(); - IssmDouble SurfaceArea(void); - void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement); - IssmDouble TimeAdapt(); - void ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss); - void ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss); int VertexConnectivity(int vertexindex); void VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");}; void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); - void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); - void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level); - bool IsZeroLevelset(int levelset_enum); - bool IsIcefront(void); - bool IsFaceOnBoundary(void); - - void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); - IssmDouble IceVolume(void); - IssmDouble IceVolumeAboveFloatation(void); - IssmDouble TotalSmb(void); - IssmDouble MassFlux(IssmDouble* segment); - IssmDouble MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id); - void ElementResponse(IssmDouble* presponse,int response_enum); - IssmDouble Masscon(IssmDouble* levelset); - IssmDouble Misfit(int modelenum,int observationenum,int weightsenum); - IssmDouble MisfitArea(int weightsenum); #ifdef _HAVE_GIA_ void GiaDeflection(Vector* wg,Vector* dwgdt,IssmDouble* x,IssmDouble* y); #endif - - void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,bool onsid); - void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index); - void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); - void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); - void InputControlUpdate(IssmDouble scalar,bool save_parameter); - - void PotentialUngrounding(Vector* potential_sheet_ungrounding); - int UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf); - /*}}}*/ /*Tria specific routines:{{{*/ @@ -148,7 +144,4 @@ void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); int GetElementType(void); - void NormalSection(IssmDouble* normal,IssmDouble* xyz_list); - void NormalTop(IssmDouble* normal,IssmDouble* xyz_list); - void NormalBase(IssmDouble* normal,IssmDouble* xyz_list); void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype); void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype); @@ -157,7 +150,7 @@ void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");}; void JacobianDeterminant(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); + void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); void JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); - void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); IssmDouble MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");}; @@ -171,21 +164,23 @@ Gauss* NewGaussTop(int order); void NodalFunctions(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); + void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); + void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; + void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss); void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); void NodalFunctionsP2(IssmDouble* basis,Gauss* gauss); - void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); - void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); - void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; - void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); + void NodalFunctionsTensor(IssmDouble* basis,Gauss* gauss); void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss); - void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss); - void NodalFunctionsTensor(IssmDouble* basis,Gauss* gauss); + void NormalBase(IssmDouble* normal,IssmDouble* xyz_list); + void NormalSection(IssmDouble* normal,IssmDouble* xyz_list); + void NormalTop(IssmDouble* normal,IssmDouble* xyz_list); void SetClone(int* minranks); void SetTemporaryElementType(int element_type_in){_error_("not implemented yet");}; Seg* SpawnSeg(int index1,int index2); IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; + void UpdateConstraintsExtrudeFromBase(void); + void UpdateConstraintsExtrudeFromTop(void); void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");}; - - void UpdateConstraintsExtrudeFromBase(void); - void UpdateConstraintsExtrudeFromTop(void); /*}}}*/