Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17329) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17330) @@ -2899,4 +2899,5 @@ if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base); if(migration_style==SubelementMigration2Enum){ + if(meshtype==Mesh2DverticalEnum) _error_("Subelement Migration 2 not implemented yet for Flowline"); gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17329) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17330) @@ -347,4 +347,21 @@ for(int iv=0;ivGaussNode(this->FiniteElement(),iv); + input->GetInputValue(&pvalue[iv],gauss); + } + delete gauss; +} +/*}}}*/ +void Element::GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype){/*{{{*/ + + _assert_(pvalue); + + int numnodes = this->GetNumberOfNodesVelocity(); + Input *input = this->GetInput(enumtype); + if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); + + /* Start looping on the number of vertices: */ + Gauss* gauss=this->NewGauss(); + for(int iv=0;ivGaussNode(this->VelocityInterpolation(),iv); input->GetInputValue(&pvalue[iv],gauss); } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17329) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17330) @@ -66,4 +66,5 @@ void GetInputListOnNodes(IssmDouble* pvalue,int enumtype); void GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); + void GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17329) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17330) @@ -573,10 +573,12 @@ bool mainlyfloating = true; + int meshtype,index1,index2; const IssmPDouble epsilon = 1.e-15; - IssmDouble phi,s1,s2,area_init,area_grounded; - IssmDouble gl[NUMVERTICES]; - IssmDouble xyz_bis[3][3]; + IssmDouble phi,s1,s2,area_init,area_grounded; + IssmDouble gl[NUMVERTICES]; + IssmDouble xyz_bis[3][3]; /*Recover parameters and values*/ + parameters->FindParam(&meshtype,MeshTypeEnum); GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum); @@ -586,82 +588,97 @@ if(gl[2]==0.) gl[2]=gl[2]+epsilon; - /*Check that not all nodes are grounded or floating*/ - if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded - phi=1; - } - else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating - phi=0; - } - else{ - /*Figure out if two nodes are floating or grounded*/ - if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=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 - /*Coordinates of point 2: same as initial point 2*/ - xyz_bis[2][0]=*(xyz_list+3*2+0); - xyz_bis[2][1]=*(xyz_list+3*2+1); - xyz_bis[2][2]=*(xyz_list+3*2+2); - - /*Portion of the segments*/ - s1=gl[2]/(gl[2]-gl[1]); - s2=gl[2]/(gl[2]-gl[0]); - - /*New point 1*/ - xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0)); - xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1)); - xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2)); - - /*New point 0*/ - xyz_bis[0][0]=*(xyz_list+3*2+0)+s2*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0)); - xyz_bis[0][1]=*(xyz_list+3*2+1)+s2*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1)); - xyz_bis[0][2]=*(xyz_list+3*2+2)+s2*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2)); - } - 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 - /*Coordinates of point 0: same as initial point 2*/ - xyz_bis[0][0]=*(xyz_list+3*0+0); - xyz_bis[0][1]=*(xyz_list+3*0+1); - xyz_bis[0][2]=*(xyz_list+3*0+2); - - /*Portion of the segments*/ - s1=gl[0]/(gl[0]-gl[1]); - s2=gl[0]/(gl[0]-gl[2]); - - /*New point 1*/ - xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0)); - xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1)); - xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2)); - - /*New point 2*/ - xyz_bis[2][0]=*(xyz_list+3*0+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0)); - xyz_bis[2][1]=*(xyz_list+3*0+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1)); - xyz_bis[2][2]=*(xyz_list+3*0+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*0+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 - /*Coordinates of point 1: same as initial point 2*/ - xyz_bis[1][0]=*(xyz_list+3*1+0); - xyz_bis[1][1]=*(xyz_list+3*1+1); - xyz_bis[1][2]=*(xyz_list+3*1+2); - - /*Portion of the segments*/ - s1=gl[1]/(gl[1]-gl[0]); - s2=gl[1]/(gl[1]-gl[2]); - - /*New point 0*/ - xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0)); - xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1)); - xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2)); - - /*New point 2*/ - xyz_bis[2][0]=*(xyz_list+3*1+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0)); - xyz_bis[2][1]=*(xyz_list+3*1+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1)); - xyz_bis[2][2]=*(xyz_list+3*1+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2)); - } - - /*Compute fraction of grounded element*/ - GetJacobianDeterminant(&area_init, xyz_list,NULL); - GetJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL); - if(mainlyfloating==true) area_grounded=area_init-area_grounded; - phi=area_grounded/area_init; - } + if(meshtype==Mesh2DverticalEnum){ + this->EdgeOnBedIndices(&index1,&index2); + if(gl[index1]>0 && gl[index2]>0) phi=1; // All grounded + else if(gl[index1]<0 && gl[index2]<0) phi=0; // All floating + else if(gl[index1]<0 && gl[index2]>0){ //index2 grounded + phi=1./(1.-gl[index1]/gl[index2]); + } + else if(gl[index2]<0 && gl[index1]>0){ //index1 grounded + phi=1./(1.-gl[index2]/gl[index1]); + } + + } + else if(meshtype==Mesh3DEnum){ + /*Check that not all nodes are grounded or floating*/ + if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded + phi=1; + } + else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating + phi=0; + } + else{ + /*Figure out if two nodes are floating or grounded*/ + if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=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 + /*Coordinates of point 2: same as initial point 2*/ + xyz_bis[2][0]=*(xyz_list+3*2+0); + xyz_bis[2][1]=*(xyz_list+3*2+1); + xyz_bis[2][2]=*(xyz_list+3*2+2); + + /*Portion of the segments*/ + s1=gl[2]/(gl[2]-gl[1]); + s2=gl[2]/(gl[2]-gl[0]); + + /*New point 1*/ + xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0)); + xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1)); + xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2)); + + /*New point 0*/ + xyz_bis[0][0]=*(xyz_list+3*2+0)+s2*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0)); + xyz_bis[0][1]=*(xyz_list+3*2+1)+s2*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1)); + xyz_bis[0][2]=*(xyz_list+3*2+2)+s2*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2)); + } + 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 + /*Coordinates of point 0: same as initial point 2*/ + xyz_bis[0][0]=*(xyz_list+3*0+0); + xyz_bis[0][1]=*(xyz_list+3*0+1); + xyz_bis[0][2]=*(xyz_list+3*0+2); + + /*Portion of the segments*/ + s1=gl[0]/(gl[0]-gl[1]); + s2=gl[0]/(gl[0]-gl[2]); + + /*New point 1*/ + xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0)); + xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1)); + xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2)); + + /*New point 2*/ + xyz_bis[2][0]=*(xyz_list+3*0+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0)); + xyz_bis[2][1]=*(xyz_list+3*0+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1)); + xyz_bis[2][2]=*(xyz_list+3*0+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*0+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 + /*Coordinates of point 1: same as initial point 2*/ + xyz_bis[1][0]=*(xyz_list+3*1+0); + xyz_bis[1][1]=*(xyz_list+3*1+1); + xyz_bis[1][2]=*(xyz_list+3*1+2); + + /*Portion of the segments*/ + s1=gl[1]/(gl[1]-gl[0]); + s2=gl[1]/(gl[1]-gl[2]); + + /*New point 0*/ + xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0)); + xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1)); + xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2)); + + /*New point 2*/ + xyz_bis[2][0]=*(xyz_list+3*1+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0)); + xyz_bis[2][1]=*(xyz_list+3*1+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1)); + xyz_bis[2][2]=*(xyz_list+3*1+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2)); + } + + /*Compute fraction of grounded element*/ + GetJacobianDeterminant(&area_init, xyz_list,NULL); + GetJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL); + if(mainlyfloating==true) area_grounded=area_init-area_grounded; + phi=area_grounded/area_init; + } + } + else _error_("mesh type "<1 || phi<0) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1"); @@ -1939,39 +1956,53 @@ void Tria::ResetFSBasalBoundaryCondition(void){ - int approximation; - int numindices; - int *indices = NULL; - IssmDouble slope; - IssmDouble xz_plane[6]; + int numnodes = this->GetNumberOfNodes(); + + int approximation; + IssmDouble* vertexonbed= NULL; + IssmDouble slope,groundedice; + IssmDouble xz_plane[6]; /*For FS only: we want the CS to be tangential to the bedrock*/ inputs->GetInputValue(&approximation,ApproximationEnum); - if(IsFloating() || !HasEdgeOnBed() || approximation!=FSApproximationEnum) return; - - /*Get number of nodes for velocity only and base*/ - int index = this->EdgeOnBedIndex(); - NodeOnEdgeIndices(&numindices,&indices,index,this->VelocityInterpolation()); - + if(IsFloating() || !HasNodeOnBed() || approximation!=FSApproximationEnum) return; + + //printf("element number %i \n",this->id); /*Get inputs*/ - Input* slope_input=inputs->GetInput(BedSlopeXEnum); _assert_(slope_input); + Input* slope_input=inputs->GetInput(BedSlopeXEnum); _assert_(slope_input); + Input* groundedicelevelset_input=inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(groundedicelevelset_input); + vertexonbed = xNew(numnodes); + this->GetInputListOnNodesVelocity(&vertexonbed[0],MeshVertexonbedEnum); /*Loop over basal nodes and update their CS*/ GaussTria* gauss = new GaussTria(); - for(int i=0;iGaussNode(this->VelocityInterpolation(),indices[i]); - slope_input->GetInputValue(&slope,gauss); - IssmDouble theta = atan(slope); - - /*New X axis New Z axis*/ - xz_plane[0]=cos(theta); xz_plane[3]=0.; - xz_plane[1]=sin(theta); xz_plane[4]=0.; - xz_plane[2]=0.; xz_plane[5]=1.; - - XZvectorsToCoordinateSystem(&this->nodes[indices[i]]->coord_system[0][0],&xz_plane[0]); + for(int i=0;iNumberofNodesVelocity();i++){ + + if(vertexonbed[i]==1){ + gauss->GaussNode(this->VelocityInterpolation(),i); + slope_input->GetInputValue(&slope,gauss); + groundedicelevelset_input->GetInputValue(&groundedice,gauss); + IssmDouble theta = atan(slope); + + if(groundedice>0){ + /*New X axis New Z axis*/ + xz_plane[0]=cos(theta); xz_plane[3]=0.; + xz_plane[1]=sin(theta); xz_plane[4]=0.; + xz_plane[2]=0.; xz_plane[5]=1.; + this->nodes[i]->DofInSSet(1); //vy + } + else{ + /*New X axis New Z axis*/ + xz_plane[0]=1.; xz_plane[3]=0; + xz_plane[1]=0.; xz_plane[4]=0; + xz_plane[2]=0; xz_plane[5]=1.; + this->nodes[i]->DofInFSet(1); //vy + } + + XZvectorsToCoordinateSystem(&this->nodes[i]->coord_system[0][0],&xz_plane[0]); + } } /*cleanup*/ - xDelete(indices); + xDelete(vertexonbed); delete gauss; }