Changeset 22379

Timestamp:

01/29/18 14:16:04 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on automating distance to GL or IF

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Element.h (modified) (2 diffs)
• classes/Elements/Penta.cpp (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (5 diffs)
• classes/Elements/Tria.h (modified) (2 diffs)
• classes/FemModel.cpp (modified) (1 diff)
• classes/FemModel.h (modified) (1 diff)
• modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryx.cpp (modified) (1 diff)
• modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -200 +200 @@
                 virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;
                 virtual void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* 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 +304 @@
                 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;

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -2194 +2194 @@
         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
         GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
-        InputDuplicate(MaskIceLevelsetEnum,PressureEnum);
 
         /*Get distance from list of segments and reset ls*/

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -525 +525 @@
 }
 /*}}}*/
-void                    Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/
+void                        Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/
 
         IssmDouble  xyz_list[NUMVERTICES][3];
@@ -762 +762 @@
 
 }/*}}}*/
+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<NUMVERTICES;j++){
+                IssmDouble dmin = 1.e+50;
+                for(int i=0;i<numsegments;i++){
+                        IssmDouble x = xyz_list[j][0];
+                        IssmDouble y = xyz_list[j][1];
+                        IssmDouble l2 = (segments[4*i+2]-segments[4*i+0])*(segments[4*i+2]-segments[4*i+0]) + (segments[4*i+3]-segments[4*i+1])*(segments[4*i+3]-segments[4*i+1]);
+
+                        /*Segment has a length of 0*/
+                        if(l2==0.){
+                                d = (x-segments[4*i+0])*(x-segments[4*i+0])+(y-segments[4*i+1])*(y-segments[4*i+1]);
+                                if(d<dmin) dmin = d;
+                                continue;
+                        }
+
+                        /*Consider the line extending the segment, parameterized as v + t (w - v).
+                         *We find projection of point p onto the line.
+                         *It falls where t = [(p-v) . (w-v)] / |w-v|^2*/
+                        IssmDouble t = ((x-segments[4*i+0])*(segments[4*i+2]-segments[4*i+0]) + (y-segments[4*i+1])*(segments[4*i+3]-segments[4*i+1]))/l2;
+                        if(t < 0.0){
+                                // Beyond the 'v' end of the segment
+                                d = (x-segments[4*i+0])*(x-segments[4*i+0])+(y-segments[4*i+1])*(y-segments[4*i+1]);
+                        }
+                        else if (t > 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(d<dmin) dmin = d;
+                }
+
+                /*Update signed distance*/
+                dmin = sqrt(dmin);
+                if(dmin>10000) 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 +1465 @@
 
 }/*}}}*/
-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 +2884 @@
         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
         GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
-        InputDuplicate(MaskIceLevelsetEnum,PressureEnum);
 
         /*Get distance from list of segments and reset ls*/
@@ -3515 +3577 @@
 }
 /*}}}*/
+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 "<<EnumToStringx(fieldenum)<<", input not found");
+        IssmDouble minvalue = input->Min();
+        IssmDouble maxvalue = input->Max();
+        if(minvalue>fieldvalue || maxvalue<fieldvalue) return;
+
+
+        /*2. Write segments*/
+        IssmDouble* xyz_list_zero = NULL;
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,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<IssmDouble>(segments->Size()+1,2,&x[0],&y[0],false));
+        }
+        xDelete<IssmDouble>(xyz_list_zero);
+}
+/*}}}*/
 void       Tria::WriteLevelsetSegment(DataSet* segments){/*{{{*/
 

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -65 +65 @@
                 void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
                 void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
+                void        CreateDistanceInputFromSegmentlist(IssmDouble* segments,int numsegments,int distanceenum);
                 int         EdgeOnBaseIndex();
                 void        EdgeOnBaseIndices(int* pindex1,int* pindex);
@@ -136 +137 @@
                 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);

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -992 +992 @@
 }
 /*}}}*/
+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;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(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<int>(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<my_rank;  i++) numseg_offset+=allsegcount[i];
+        for(int i=0;i<num_procs;i++) numseg+=allsegcount[i];
+        IssmDouble* segmentlist    = xNewZeroInit<IssmDouble>(4*numseg);
+        IssmDouble* allsegmentlist = xNewZeroInit<IssmDouble>(4*numseg);
+        for(int i=0;i<segments->Size();i++){
+                Contour<IssmDouble>* segment=(Contour<IssmDouble>*)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<IssmDouble>(segmentlist);
+        xDelete<int>(allsegcount);
+
+        /*3: Add distance input to all elements*/
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(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<IssmDouble>(allsegmentlist);
+}/*}}}*/
 void FemModel::Divergencex(IssmDouble* pdiv){/*{{{*/
 

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -107 +107 @@
                 void MinVyx(IssmDouble* presponse);
                 void MinVzx(IssmDouble* presponse);
+                void DistanceToFieldValue(int fieldenum,IssmDouble fieldvalue,int distanceenum);
                 void ResetLevelset();
                 void StrainRateparallelx();

issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryx.cpp

@@ -20 +20 @@
         /*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: */

issm/trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp

@@ -37 +37 @@
         for(i=i0;i<i1;i++){
 
-                IssmDouble d0=pow(10.,10);
-
+                IssmDouble d0=1.e+10;
                 IssmDouble xi,yi;