Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22378) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22379) @@ -200,4 +200,5 @@ virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0; virtual void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum)=0; + virtual void CreateDistanceInputFromSegmentlist(IssmDouble* segments,int numsegments,int distanceenum){_error_("not implemented yet");}; virtual void ElementResponse(IssmDouble* presponse,int response_enum)=0; virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; @@ -303,4 +304,5 @@ virtual void VerticalSegmentIndicesBase(int** pindices,int* pnumseg)=0; virtual void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; + virtual void WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue){_error_("not implemented yet");}; virtual void WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");}; virtual void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 22378) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 22379) @@ -2194,5 +2194,4 @@ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); - InputDuplicate(MaskIceLevelsetEnum,PressureEnum); /*Get distance from list of segments and reset ls*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22378) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22379) @@ -525,5 +525,5 @@ } /*}}}*/ -void Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/ +void Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/ IssmDouble xyz_list[NUMVERTICES][3]; @@ -762,4 +762,67 @@ }/*}}}*/ +void Tria::CreateDistanceInputFromSegmentlist(IssmDouble* segments,int numsegments,int distanceenum){/*{{{*/ + + /*Intermediaries*/ + IssmDouble d,xn,yn; + + /*Get current field and vertex coordinates*/ + IssmDouble ls[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + GetInputListOnVertices(&ls[0],distanceenum); + + /*Get distance from list of segments and reset ls*/ + for(int j=0;j 1.0){ + // Beyond the 'w' end of the segment + d = (x-segments[4*i+2])*(x-segments[4*i+2])+(y-segments[4*i+3])*(y-segments[4*i+3]); + } + else{ + // Projection falls on the segment + xn = segments[4*i+0] + t * (segments[4*i+2] - segments[4*i+0]); + yn = segments[4*i+1] + t * (segments[4*i+3] - segments[4*i+1]); + d = (x-xn)*(x-xn)+(y-yn)*(y-yn); + } + + if(d10000) dmin=10000; + if(ls[j]>0){ + ls[j] = dmin; + } + else{ + ls[j] = - dmin; + } + } + + /*Update Levelset*/ + this->inputs->AddInput(new TriaInput(distanceenum,&ls[0],P1Enum)); +} +/*}}}*/ int Tria::EdgeOnBaseIndex(void){/*{{{*/ @@ -1402,5 +1465,5 @@ }/*}}}*/ -void Tria::GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/ +void Tria::GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/ /* GetLevelsetIntersection computes: @@ -2821,5 +2884,4 @@ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); - InputDuplicate(MaskIceLevelsetEnum,PressureEnum); /*Get distance from list of segments and reset ls*/ @@ -3515,4 +3577,30 @@ } /*}}}*/ +void Tria::WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue){/*{{{*/ + + _assert_(fieldvalue==0.); //field value != 0 not implemented yet + + /*1. check that we do cross fieldvalue in this element*/ + Input* input = inputs->GetInput(fieldenum); + if(!input) _error_("Cannot calculate distance to "<Min(); + IssmDouble maxvalue = input->Max(); + if(minvalue>fieldvalue || maxvaluevertices,NUMVERTICES); + this->ZeroLevelsetCoordinates(&xyz_list_zero,&xyz_list[0][0],fieldenum); + if(xyz_list_zero){ + IssmDouble x[2],y[2]; + x[0] = xyz_list_zero[0*3 + 0]; x[1] = xyz_list_zero[1*3 + 0]; + y[0] = xyz_list_zero[0*3 + 1]; y[1] = xyz_list_zero[1*3 + 1]; + segments->AddObject(new Contour(segments->Size()+1,2,&x[0],&y[0],false)); + } + xDelete(xyz_list_zero); +} +/*}}}*/ void Tria::WriteLevelsetSegment(DataSet* segments){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22378) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22379) @@ -65,4 +65,5 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); + void CreateDistanceInputFromSegmentlist(IssmDouble* segments,int numsegments,int distanceenum); int EdgeOnBaseIndex(); void EdgeOnBaseIndices(int* pindex1,int* pindex); @@ -136,4 +137,5 @@ void VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");}; void VerticalSegmentIndicesBase(int** pindices,int* pnumseg){_error_("not implemented yet");}; + void WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue); void WriteLevelsetSegment(DataSet* segments); void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22378) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22379) @@ -992,4 +992,75 @@ } /*}}}*/ +void FemModel::DistanceToFieldValue(int fieldenum,IssmDouble fieldvalue,int distanceenum){/*{{{*/ + + /*recover my_rank:*/ + int my_rank = IssmComm::GetRank(); + int num_procs = IssmComm::GetSize(); + + /*Get domain type (2d or 3d)*/ + int domaintype; + this->parameters->FindParam(&domaintype,DomainTypeEnum); + + /*1: go throug all elements of this partition and figure out how many + * segments we have (corresopnding to field = value)*/ + DataSet* segments=new DataSet(); + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + if(!element->IsOnBase()) continue; + Element* basalelement = element->SpawnBasalElement(); + basalelement->WriteFieldIsovalueSegment(segments,fieldenum,fieldvalue); + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + } + + /*2: now get the segments from all partitions*/ + int segcount=segments->Size(); + int* allsegcount=xNew(num_procs); + ISSM_MPI_Gather(&segcount,1,ISSM_MPI_INT,allsegcount,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(allsegcount,num_procs,ISSM_MPI_INT,0,IssmComm::GetComm()); + + /* Every cpu should start its own dof count at the end of the dofcount from cpu-1*/ + int numseg_offset=0; + int numseg=0; + for(int i=0;i(4*numseg); + IssmDouble* allsegmentlist = xNewZeroInit(4*numseg); + for(int i=0;iSize();i++){ + Contour* segment=(Contour*)segments->GetObjectByOffset(i); + _assert_(segment->nods == 2); + segmentlist[(numseg_offset+i)*4 + 0] = segment->x[0]; + segmentlist[(numseg_offset+i)*4 + 1] = segment->y[0]; + segmentlist[(numseg_offset+i)*4 + 2] = segment->x[1]; + segmentlist[(numseg_offset+i)*4 + 3] = segment->y[1]; + } + ISSM_MPI_Allreduce((void*)segmentlist,(void*)allsegmentlist,4*numseg,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + delete segments; + xDelete(segmentlist); + xDelete(allsegcount); + + /*3: Add distance input to all elements*/ + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + if(!element->IsOnBase()) continue; + element->CreateDistanceInputFromSegmentlist(allsegmentlist,numseg,distanceenum); + } + + + /*Extrude if necessary*/ + int elementtype; + this->parameters->FindParam(&elementtype,MeshElementtypeEnum); + if(elementtype==PentaEnum){ + InputExtrudex(this,distanceenum,-1); + } + else if(elementtype==TriaEnum){ + /*no need to extrude*/ + } + else{ + _error_("not implemented yet"); + } + + /*Clean up and return*/ + xDelete(allsegmentlist); +}/*}}}*/ void FemModel::Divergencex(IssmDouble* pdiv){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22378) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22379) @@ -107,4 +107,5 @@ void MinVyx(IssmDouble* presponse); void MinVzx(IssmDouble* presponse); + void DistanceToFieldValue(int fieldenum,IssmDouble fieldvalue,int distanceenum); void ResetLevelset(); void StrainRateparallelx(); Index: /issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryx.cpp (revision 22378) +++ /issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryx.cpp (revision 22379) @@ -20,9 +20,9 @@ /*initialize thread parameters: */ DistanceToMaskBoundaryxThreadStruct gate; - gate.distance = distance; - gate.x = x; - gate.y = y; - gate.mask = mask; - gate.nods = nods; + gate.distance = distance; + gate.x = x; + gate.y = y; + gate.mask = mask; + gate.nods = nods; /*launch the thread manager with DistanceToMaskBoundaryxt as a core: */ Index: /issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp (revision 22378) +++ /issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp (revision 22379) @@ -37,6 +37,5 @@ for(i=i0;i