Changeset 21402

Timestamp:

11/18/16 13:56:00 (8 years ago)

Author:

youngmc3

Message:

CHG: calving 3D

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Penta.cpp (modified) (3 diffs)
• Elements/Penta.h (modified) (2 diffs)
• FemModel.cpp (modified) (3 diffs)

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

@@ -174 +174 @@
 void       Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
         _error_("Not supported yet!");
+}
+/*}}}*/
+void       Penta::CalvingRateDev(){/*{{{*/
+
+        if(!this->IsOnBase()) return;
+
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+        IssmDouble  calvingratex[NUMVERTICES];
+        IssmDouble  calvingratey[NUMVERTICES];
+        IssmDouble  calvingrate[NUMVERTICES];
+        IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
+        IssmDouble  sigma_vm,sigma_max,epse_2,groundedice;
+
+        /* Get node coordinates and dof list: */
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+        /*Depth average B for stress calculation*/
+        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+        /*Retrieve all inputs and parameters we will need*/
+        Input* vx_input = inputs->GetInput(VxAverageEnum); _assert_(vx_input);
+        Input* vy_input = inputs->GetInput(VyAverageEnum); _assert_(vy_input);
+        Input* gr_input = inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input);
+        IssmDouble  B   = this->GetMaterialParameter(MaterialsRheologyBbarEnum);
+        IssmDouble  n   = this->GetMaterialParameter(MaterialsRheologyNEnum);
+
+        /* Start looping on the number of vertices: */
+        GaussPenta* gauss=new GaussPenta();
+        for(int iv=0;iv<3;iv++){
+                gauss->GaussVertex(iv);
+
+                /*Get velocity components and thickness*/
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                gr_input->GetInputValue(&groundedice,gauss);
+                vel=sqrt(vx*vx+vy*vy)+1.e-14;
+
+                /*Compute strain rate and viscosity: */
+                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+
+                /*Get Eigen values*/
+                Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
+                _assert_(!xIsNan<IssmDouble>(lambda1));
+                _assert_(!xIsNan<IssmDouble>(lambda2));
+
+                /*Process Eigen values (only account for extension)*/
+                lambda1 = max(lambda1,0.);
+                lambda2 = max(lambda2,0.);
+
+                /*Calculate sigma_vm*/
+                epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
+                sigma_vm  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
+                sigma_max = 1000.e+3;
+
+                if(groundedice<0) sigma_max=200.e+3;
+
+                /*Assign values*/
+                calvingratex[iv]=vx*sigma_vm/sigma_max;
+                calvingratey[iv]=vy*sigma_vm/sigma_max;
+                calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
+        }
+
+        /*Add input*/
+        this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
+        this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
+        this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+
+        this->InputExtrude(CalvingratexEnum,-1);
+        this->InputExtrude(CalvingrateyEnum,-1);
+        this->InputExtrude(CalvingCalvingrateEnum,-1);
+
+        /*Clean up and return*/
+        delete gauss;
 }
 /*}}}*/
@@ -641 +717 @@
 
                 /*Create gauss point in the middle of the basal edge*/
-                Gauss* gauss=NewGaussBase(1);
+                Gauss* gauss=NewGauss();
                 gauss->GaussPoint(0);
 
@@ -2213 +2289 @@
 }
 /*}}}*/
+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);
+        InputDuplicate(MaskIceLevelsetEnum,PressureEnum);
+
+        /*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 PentaInput(MaskIceLevelsetEnum,&ls[0],P1Enum));
+}
+/*}}}*/
 void       Penta::SetClone(int* minranks){/*{{{*/
 
         _error_("not implemented yet");
 }
+
 /*}}}*/
 void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/

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

@@ -48 +48 @@
                 void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
                 void           AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
+                void           CalvingRateDev();
                 void           CalvingRateLevermann();
                 IssmDouble     CharacteristicLength(void){_error_("not implemented yet");};
@@ -146 +147 @@
                 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);

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

@@ -1734 +1734 @@
         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)*/
@@ -1739 +1743 @@
         for(int i=0;i<elements->Size();i++){
                 Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
-                element->WriteLevelsetSegment(segments);
+                if(!element->IsOnBase()) continue;
+                Element* basalelement = element->SpawnBasalElement();
+                basalelement->WriteLevelsetSegment(segments);
+                if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
         }
 
@@ -1771 +1778 @@
         for(int i=0;i<elements->Size();i++){
                 Element* element=dynamic_cast<Element*>(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");
         }