Changeset 22989


Ignore:
Timestamp:
07/20/18 16:15:01 (7 years ago)
Author:
Mathieu Morlighem
Message:

CHG: simplifying levelset distance

Location:
issm/trunk-jpl/src/c
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r22974 r22989  
    319319                virtual void       ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
    320320                virtual void       WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue){_error_("not implemented yet");};
    321                 virtual void       WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");};
    322321                virtual void       ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0;
    323322

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

    r22988 r22989  
    40664066        if(!input) _error_("Cannot calculate distance to "<<EnumToStringx(fieldenum)<<", input not found");
    40674067        IssmDouble minvalue = input->Min();
     4068        if(minvalue>fieldvalue) return;
    40684069        IssmDouble maxvalue = input->Max();
    4069         if(minvalue>fieldvalue || maxvalue<fieldvalue) return;
    4070 
    4071         /* check #2: If only one vertex is on fieldvalue, there is no segment here */
    4072         IssmDouble lsf[NUMVERTICES];
    4073         this->GetInputListOnVertices(&lsf[0],fieldenum);
    4074         int nrice=0;       
    4075         for(int i=0;i<NUMVERTICES;i++) if(lsf[i]==fieldvalue) nrice++;
    4076         if(nrice==1) return;
    4077 
    4078         /*2. Write segments*/
    4079         IssmDouble* xyz_list_zero = NULL;
     4070        if(maxvalue<fieldvalue) return;
     4071
     4072        /*2. Find coordinates of where levelset crosses 0*/
     4073        int         numiceverts;
     4074        IssmDouble  s[2],x[2],y[2];
     4075        int        *indices = NULL;
     4076        this->GetLevelsetIntersection(&indices, &numiceverts,&s[0],fieldenum,fieldvalue);
     4077        _assert_(numiceverts);
     4078
     4079        /*3 Write coordinates*/
    40804080        IssmDouble  xyz_list[NUMVERTICES][3];
    40814081        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
    4082         this->ZeroLevelsetCoordinates(&xyz_list_zero,&xyz_list[0][0],fieldenum);
    4083         if(xyz_list_zero){
    4084                 IssmDouble x[2],y[2];
    4085                 x[0] = xyz_list_zero[0*3 + 0]; x[1] = xyz_list_zero[1*3 + 0];
    4086                 y[0] = xyz_list_zero[0*3 + 1]; y[1] = xyz_list_zero[1*3 + 1];
    4087                 segments->AddObject(new Contour<IssmDouble>(segments->Size()+1,2,&x[0],&y[0],false));
    4088         }
    4089         xDelete<IssmDouble>(xyz_list_zero);
    4090 }
    4091 /*}}}*/
    4092 void       Tria::WriteLevelsetSegment(DataSet* segments){/*{{{*/
    4093 
    4094         if(!this->IsZeroLevelset(MaskIceLevelsetEnum)) return;
    4095 
    4096         IssmDouble* xyz_list_zero = NULL;
    4097         IssmDouble  xyz_list[NUMVERTICES][3];
    4098         ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
    4099         this->ZeroLevelsetCoordinates(&xyz_list_zero,&xyz_list[0][0], MaskIceLevelsetEnum);
    4100         if(xyz_list_zero){
    4101                 IssmDouble x[2];
    4102                 IssmDouble y[2];
    4103                 x[0] = xyz_list_zero[0*3 + 0]; x[1] = xyz_list_zero[1*3 + 0];
    4104                 y[0] = xyz_list_zero[0*3 + 1]; y[1] = xyz_list_zero[1*3 + 1];
    4105                 segments->AddObject(new Contour<IssmDouble>(segments->Size()+1,2,&x[0],&y[0],false));
    4106         }
    4107         xDelete<IssmDouble>(xyz_list_zero);
    4108 
    4109 //      IssmDouble ls[NUMVERTICES];
    4110 //      IssmDouble  xyz_list[NUMVERTICES][3];
    4111 //
    4112 //      if(IsIceInElement()){
    4113 //
    4114 //              /*Retrieve all inputs and parameters*/
    4115 //              GetInputListOnVertices(&ls[0],levelset_enum);
    4116 //
    4117 //              /*If the level set is awlays <0, there is no ice front here*/
    4118 //              bool iszerols= false;
    4119 //              if(IsIceInElement()){
    4120 //                      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.)){
    4121 //                              iszerols = true;
    4122 //                      }
    4123 //              }
    4124 //
    4125 //              if(iszerols){
    4126 //                      /*OK we have one segment!*/
    4127 //                      ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4128 //                      int count = 0;
    4129 //                      IssmDouble x[2];
    4130 //                      IssmDouble y[2];
    4131 //
    4132 //                      for(int i=0;i<NUMVERTICES,i++){
    4133 //                              int index1 = i;
    4134 //                              int index1 = (i+1)%3;
    4135 //                              if(ls[index1]<=0 && ls[index2]>=0){
    4136 //
    4137 //                              }
    4138 //
    4139 //                      }
    4140 //                      Contour* segment = new Contour<IssmDouble>(segment->Size()+1,2,x,y,false);
    4141 //              }
    4142 //
    4143 //      }
    4144 //
    4145 //      _error_("STOP");
     4082        int counter = 0;
     4083        if((numiceverts>0) && (numiceverts<NUMVERTICES)){
     4084                for(int i=0;i<numiceverts;i++){
     4085                        for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
     4086                                x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
     4087                                y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
     4088                                counter++;
     4089                        }
     4090                }
     4091        }
     4092        else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
     4093                IssmDouble lsf[NUMVERTICES];
     4094                this->GetInputListOnVertices(&lsf[0],fieldenum);
     4095                for(int i=0;i<NUMVERTICES;i++){
     4096                        if(lsf[indices[i]]==0.){
     4097                                x[counter]=xyz_list[indices[i]][0];
     4098                                y[counter]=xyz_list[indices[i]][1];
     4099                                counter++;
     4100                        }
     4101                        if(counter==2) break;
     4102                }
     4103                if(counter==1){
     4104                        /*We actually have only 1 vertex on levelset, write a single point as a segment*/
     4105                        _error_("TODO");
     4106                        counter==2;
     4107                }
     4108        }
     4109        else{
     4110                _error_("not sure what's going on here...");
     4111        }
     4112
     4113        /*4. Write segment*/
     4114        _assert_(counter==2);
     4115        segments->AddObject(new Contour<IssmDouble>(segments->Size()+1,2,&x[0],&y[0],false));
     4116
     4117        /*Cleanup and return*/
     4118        xDelete<int>(indices);
    41464119}
    41474120/*}}}*/
    … …  
    41514124
    41524125        /*Intermediaries*/
    4153         const int dim=3;
    4154         int numiceverts;
    4155         int i, n, e, counter;
    4156         IssmDouble s[2];
    4157         int* indices=NULL;
    4158         IssmDouble* xyz_zero=NULL;
     4126        int         numiceverts;
     4127        int         i,n,e,counter;
     4128        IssmDouble  s[2];
     4129        int        *indices  = NULL;
     4130        IssmDouble *xyz_zero = NULL;
    41594131
    41604132        this->GetLevelsetIntersection(&indices, &numiceverts, s,levelsetenum,0.);
    41614133       
    4162         //TODO: check if for 2 iceverts front segment is oriented in CCW way
    4163        
    4164         if(numiceverts>0) xyz_zero=xNew<IssmDouble>(2*dim);
    4165         if((numiceverts>0)&&(numiceverts<NUMVERTICES)){
     4134        printf("numiceverts = %i\n",numiceverts);
     4135        if(numiceverts>0) xyz_zero=xNew<IssmDouble>(2*3);
     4136
     4137        if((numiceverts>0) && (numiceverts<NUMVERTICES)){
    41664138                counter=0;
    4167                 for(i=0;i<numiceverts;i++){     // iterate over ice vertices
     4139                for(i=0;i<numiceverts;i++){
    41684140                        for(n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
    4169                                 for(e=0;e<dim;e++){ // spatial direction
    4170                                         int ind_ice             =dim*indices[i]+e;
    4171                                         int ind_noice   =dim*indices[n]+e;
    4172                                         int ind                 =dim*counter+e;
     4141                                for(e=0;e<3;e++){ // spatial direction
     4142                                        int ind_ice   = 3 *indices[i]+e;
     4143                                        int ind_noice = 3 *indices[n]+e;
     4144                                        int ind       = 3 *counter+e;
    41734145                                        xyz_zero[ind]=xyz_list[ind_ice]+s[counter]*(xyz_list[ind_noice]-xyz_list[ind_ice]);
    41744146                                }
    … …  
    41834155                for(i=0;i<NUMVERTICES;i++){
    41844156                        if(lsf[indices[i]]==0.){
    4185                                 for(e=0;e<dim;e++)      xyz_zero[dim*counter+e]=xyz_list[dim*indices[i]+e];
     4157                                for(e=0;e<3;e++)        xyz_zero[3*counter+e]=xyz_list[3*indices[i]+e];
    41864158                                counter++;
    41874159                        }

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

    r22974 r22989  
    144144                void        VerticalSegmentIndicesBase(int** pindices,int* pnumseg){_error_("not implemented yet");};
    145145                void                    WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue);
    146                 void                    WriteLevelsetSegment(DataSet* segments);
    147146                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    148147

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

    r22988 r22989  
    20302030                if(!element->IsOnBase()) continue;
    20312031                Element* basalelement = element->SpawnBasalElement();
    2032                 basalelement->WriteLevelsetSegment(segments);
     2032                basalelement->WriteFieldIsovalueSegment(segments,MaskIceLevelsetEnum,0.);
    20332033                if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    20342034        }

  • TabularUnified issm/trunk-jpl/src/c/cores/movingfront_core.cpp

    r22987 r22989  
    8585
    8686        /*Reset levelset if needed*/
    87         printf("-------------- file: movingfront_core.cpp line: %i\n",__LINE__);
    8887        if(reinit_frequency && (step%reinit_frequency==0)){
    8988                if(VerboseSolution()) _printf0_("   reinitializing level set\n");
Note: See TracChangeset for help on using the changeset viewer.