Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 22986) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 22987) @@ -663,57 +663,4 @@ } }/*}}}*/ -void LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector* vec_signed_dist, Element* element){/*{{{*/ - - if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) - return; - - /* Intermediaries */ - const int dim=3; - int i,d; - int numvertices=element->GetNumberOfVertices(); - IssmDouble s0[dim], s1[dim], v[dim]; - IssmDouble dist,lsf_old; - - IssmDouble* lsf = xNew(numvertices); - IssmDouble* sign_lsf = xNew(numvertices); - IssmDouble* signed_dist = xNew(numvertices); - IssmDouble* xyz_list = NULL; - IssmDouble* xyz_list_zero = NULL; - - /* retrieve inputs and parameters */ - element->GetVerticesCoordinates(&xyz_list); - element->GetInputListOnVertices(lsf,MaskIceLevelsetEnum); - - /* get sign of levelset function */ - for(i=0;i=0.?1.:-1.); - - element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum); - for(d=0;dGetValue(&lsf_old, element->vertices[i]->Sid()); - if(xIsNan(lsf_old) || fabs(signed_dist[i])SetValue(element->vertices[i]->Sid(),signed_dist[i],INS_VAL); - } - - xDelete(lsf); - xDelete(sign_lsf); - xDelete(signed_dist); - xDelete(xyz_list); - xDelete(xyz_list_zero); -}/*}}}*/ void LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 22986) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h (revision 22987) @@ -32,5 +32,4 @@ void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); void InputUpdateFromSolution(IssmDouble* solution,Element* element); - void SetDistanceToZeroLevelsetElement(Vector* vec_signed_dist, Element* element); void UpdateConstraints(FemModel* femmodel); }; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22986) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22987) @@ -3301,60 +3301,19 @@ } /*}}}*/ -void Tria::ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){/*{{{*/ - - /*Intermediaries*/ - IssmDouble d,xn,yn; +void Tria::ResetLevelsetFromSegmentlist(IssmDouble* distances){/*{{{*/ + /*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; + ls[j] = distances[this->vertices->Pid()]; } else{ - ls[j] = - dmin; + ls[j] = - distances[this->vertices->Pid()]; } } Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22986) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22987) @@ -1050,24 +1050,54 @@ /*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"); - } + IssmDouble* distances = xNew(vertices->Size()); + IssmDouble d,xn,yn; + for(int i=0;iSize();i++){ + Vertex* vertex=dynamic_cast(this->vertices->GetObjectByOffset(i)); + IssmDouble dmin = 1.e+50; + + for(int i=0;ix; + IssmDouble y = vertex->y; + IssmDouble l2 = (allsegmentlist[4*i+2]-allsegmentlist[4*i+0])*(allsegmentlist[4*i+2]-allsegmentlist[4*i+0]) + (allsegmentlist[4*i+3]-allsegmentlist[4*i+1])*(allsegmentlist[4*i+3]-allsegmentlist[4*i+1]); + + /*Segment has a length of 0*/ + if(l2==0.){ + d = (x-allsegmentlist[4*i+0])*(x-allsegmentlist[4*i+0])+(y-allsegmentlist[4*i+1])*(y-allsegmentlist[4*i+1]); + if(d 1.0){ + // Beyond the 'w' end of the segment + d = (x-allsegmentlist[4*i+2])*(x-allsegmentlist[4*i+2])+(y-allsegmentlist[4*i+3])*(y-allsegmentlist[4*i+3]); + } + else{ + // Projection falls on the segment + xn = allsegmentlist[4*i+0] + t * (allsegmentlist[4*i+2] - allsegmentlist[4*i+0]); + yn = allsegmentlist[4*i+1] + t * (allsegmentlist[4*i+3] - allsegmentlist[4*i+1]); + d = (x-xn)*(x-xn)+(y-yn)*(y-yn); + } + + if(dpid] = sqrt(dmin); + } + + //element->CreateDistanceInputFromSegmentlist(allsegmentlist,numseg,distanceenum); + InputUpdateFromVectorx(this,distances,distanceenum,VertexPIdEnum); + printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); /*Clean up and return*/ + xDelete(distances); xDelete(allsegmentlist); }/*}}}*/ Index: /issm/trunk-jpl/src/c/cores/movingfront_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/movingfront_core.cpp (revision 22986) +++ /issm/trunk-jpl/src/c/cores/movingfront_core.cpp (revision 22987) @@ -85,4 +85,5 @@ /*Reset levelset if needed*/ + printf("-------------- file: movingfront_core.cpp line: %i\n",__LINE__); if(reinit_frequency && (step%reinit_frequency==0)){ if(VerboseSolution()) _printf0_(" reinitializing level set\n");