Changeset 19476

Timestamp:

08/03/15 13:03:29 (10 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added ResetLevelset function

Location:

issm/trunk-jpl/src/c

Files:

• analyses/LevelsetAnalysis.cpp (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (5 diffs)
• classes/Elements/Tria.h (modified) (1 diff)
• classes/FemModel.cpp (modified) (1 diff)
• classes/FemModel.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

@@ -143 +143 @@
         basalelement->FindParam(&calvinglaw,CalvingLawEnum);
         basalelement->FindParam(&stabilization,LevelsetStabilizationEnum);
+        calvinglaw=DefaultCalvingEnum;
         switch(domaintype){
                 case Domain2DverticalEnum:   dim = 1; break;

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

@@ -178 +178 @@
                 virtual void       CalvingRatePi(void)=0;
                 virtual void       CalvingRateDev(void){_error_("not implemented yet");};
+                virtual void       WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");};
+                virtual void       ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");};
                 virtual IssmDouble CharacteristicLength(void)=0;
                 virtual void       ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;

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

@@ -350 +350 @@
         IssmDouble  calvingrate[NUMVERTICES];
         IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
-        IssmDouble  sigma_vm,sigma_max;
+        IssmDouble  sigma_vm,sigma_max,epse_2;
 
         /* Get node coordinates and dof list: */
@@ -384 +384 @@
 
                 /*Calculate sigma_vm*/
-                sigma_vm  = sqrt(3./2.) * B * pow(lambda1*lambda1 + lambda2*lambda2,1./(2.*n));
+                epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
+                sigma_vm  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
+                //sigma_max = 125.e+3;
                 sigma_max = 350.e+3;
+                sigma_max = 450.e+3;
+                sigma_max = 800.e+3; //too much
+                sigma_max = 700.e+3;
+                sigma_max = 670.e+3;
+                //sigma_max = 550.e+3;
 
                 /*Assign values*/
@@ -400 +407 @@
         /*Clean up and return*/
         delete gauss;
+}
+/*}}}*/
+void       Tria::WriteLevelsetSegment(DataSet* segments){/*{{{*/
+
+        if(!this->IsZeroLevelset(MaskIceLevelsetEnum)) return;
+
+        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], MaskIceLevelsetEnum);
+        if(xyz_list_zero){
+                IssmDouble x[2];
+                IssmDouble 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);
+
+//      IssmDouble ls[NUMVERTICES];
+//      IssmDouble  xyz_list[NUMVERTICES][3];
+//
+//      if(IsIceInElement()){
+//
+//              /*Retrieve all inputs and parameters*/
+//              GetInputListOnVertices(&ls[0],levelset_enum);
+//
+//              /*If the level set is awlays <0, there is no ice front here*/
+//              bool iszerols= false;
+//              if(IsIceInElement()){
+//                      if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
+//                              iszerols = true;
+//                      }
+//              }
+//
+//              if(iszerols){
+//                      /*OK we have one segment!*/
+//                      ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+//                      int count = 0;
+//                      IssmDouble x[2];
+//                      IssmDouble y[2];
+//
+//                      for(int i=0;i<NUMVERTICES,i++){
+//                              int index1 = i;
+//                              int index1 = (i+1)%3;
+//                              if(ls[index1]<=0 && ls[index2]>=0){
+//
+//                              }
+//
+//                      }
+//                      Contour* segment = new Contour<IssmDouble>(segment->Size()+1,2,x,y,false);
+//              }
+//
+//      }
+//
+//      _error_("STOP");
+}
+/*}}}*/
+void       Tria::ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){/*{{{*/
+
+        /*Intermediaries*/
+        double 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);
+        InputDuplicate(MaskIceLevelsetEnum,PressureEnum);
+
+        /*Get distance from list of segments and reset ls*/
+        for(int j=0;j<NUMVERTICES;j++){
+                double 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*/
+                        double 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(MaskIceLevelsetEnum,&ls[0],P1Enum));
 }
 /*}}}*/
@@ -3144 +3271 @@
         IssmDouble* xyz_zero=NULL;
 
-        this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);
+        this->GetLevelsetIntersection(&indices, &numiceverts, s,levelsetenum,0.);
         
         //TODO: check if for 2 iceverts front segment is oriented in CCW way
@@ -3165 +3292 @@
         else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
                 IssmDouble lsf[NUMVERTICES];
-                this->GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+                this->GetInputListOnVertices(&lsf[0],levelsetenum);
                 counter=0;
                 for(i=0;i<NUMVERTICES;i++){

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

@@ -55 +55 @@
                 void                    CalvingRatePi();
                 void                    CalvingRateDev();
+                void                    WriteLevelsetSegment(DataSet* segments);
+                void        ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);
                 IssmDouble  CharacteristicLength(void);
                 void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);

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

@@ -842 +842 @@
                 (element->*function)();
         }
+
+}
+/*}}}*/
+void FemModel::ResetLevelset(void){/*{{{*/
+
+        /*recover my_rank:*/
+        int my_rank   = IssmComm::GetRank();
+        int num_procs = IssmComm::GetSize();
+
+        /*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;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+                element->WriteLevelsetSegment(segments);
+        }
+
+        /*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: update levelset for all elements*/
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+                element->ResetLevelsetFromSegmentlist(allsegmentlist,numseg);
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(allsegmentlist);
 
 }

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

@@ -94 +94 @@
                 void CalvingRatePix();
                 void CalvingRateDevx();
+                void ResetLevelset();
                 #ifdef  _HAVE_DAKOTA_
                 void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);