Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22990) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22991) @@ -290,5 +290,4 @@ virtual void ResetFSBasalBoundaryCondition()=0; virtual void ResetHooks()=0; - virtual void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");}; virtual void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M)=0; virtual void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 22990) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 22991) @@ -2286,67 +2286,4 @@ if(this->hneighbors) this->hneighbors->reset(); -} -/*}}}*/ -void Penta::ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){/*{{{*/ - - /*Intermediaries*/ - IssmDouble d,xn,yn; - - /*Get current levelset and vertex coordinates*/ - IssmDouble ls[NUMVERTICES]; - IssmDouble xyz_list[NUMVERTICES][3]; - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); - - /*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 PentaInput(MaskIceLevelsetEnum,&ls[0],P1Enum)); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22990) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22991) @@ -149,5 +149,4 @@ void ResetFSBasalBoundaryCondition(void); void ResetHooks(); - void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments); void SetClone(int* minranks); void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N,int M); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22990) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22991) @@ -873,5 +873,6 @@ /*Update signed distance*/ dmin = sqrt(dmin); - // if(dmin>10000) dmin=10000; + /*FIXME: do we really need this?*/ + if(distanceenum==MaskIceLevelsetEnum) if(dmin>10000) dmin=10000; if(ls[j]>0){ ls[j] = dmin; @@ -3301,67 +3302,4 @@ } /*}}}*/ -void Tria::ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){/*{{{*/ - - /*Intermediaries*/ - IssmDouble d,xn,yn; - - /*Get current levelset and vertex coordinates*/ - IssmDouble ls[NUMVERTICES]; - IssmDouble xyz_list[NUMVERTICES][3]; - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); - - /*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(MaskIceLevelsetEnum,&ls[0],P1Enum)); -} -/*}}}*/ void Tria::SetClone(int* minranks){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22990) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22991) @@ -2015,71 +2015,5 @@ void FemModel::ResetLevelset(void){/*{{{*/ - /*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 levelset = 0)*/ - 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,MaskIceLevelsetEnum,0.); - 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: update levelset for all elements*/ - for(int i=0;iSize();i++){ - Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); - if(!element->IsOnBase()) continue; - element->ResetLevelsetFromSegmentlist(allsegmentlist,numseg); - } - - /*Extrude if necessary*/ - int elementtype; - this->parameters->FindParam(&elementtype,MeshElementtypeEnum); - if(elementtype==PentaEnum){ - InputExtrudex(this,MaskIceLevelsetEnum,-1); - } - else if(elementtype==TriaEnum){ - /*no need to extrude*/ - } - else{ - _error_("not implemented yet"); - } - - /*Clean up and return*/ - xDelete(allsegmentlist); + this->DistanceToFieldValue(MaskIceLevelsetEnum,0.,MaskIceLevelsetEnum); }