Changeset 18786

Timestamp:

11/14/14 16:42:09 (10 years ago)

Author:

Eric.Larour

Message:

NEW: introduced new Masscon output definition, relying on a levelset definition of the masscon shape.

Location:

issm/trunk-jpl/src

Files:

• c/Makefile.am (modified) (1 diff)
• c/classes/Elements/Element.h (modified) (2 diffs)
• c/classes/Elements/Penta.h (modified) (2 diffs)
• c/classes/Elements/Seg.h (modified) (2 diffs)
• c/classes/Elements/Tetra.h (modified) (2 diffs)
• c/classes/Elements/Tria.cpp (modified) (2 diffs)
• c/classes/Elements/Tria.h (modified) (2 diffs)
• c/classes/Masscon.h (added)
• c/classes/classes.h (modified) (1 diff)
• c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (4 diffs)
• m/classes/masscon.m (added)
• m/enum/EnumDefinitions.py (modified) (1 diff)
• m/enum/MassconEnum.m (added)
• m/enum/MassconLevelsetEnum.m (added)
• m/enum/MassconNameEnum.m (added)
• m/exp/exp_to_levelset.m (modified) (2 diffs)
• m/plot/applyoptions.m (modified) (1 diff)

issm/trunk-jpl/src/c/Makefile.am ¶

@@ -437 +437 @@
                                         ./classes/Massfluxatgate.h \
                                         ./classes/Misfit.h \
+                                        ./classes/Masscon.h \
                                         ./modules/ModelProcessorx/CreateOutputDefinitions.cpp\
                                         ./modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.h\

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

@@ -205 +205 @@
                 virtual bool   IsOnSurface()=0;
                 virtual void   GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0;
+                virtual void   GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues)=0;
                 virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0;
                 virtual void   GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
@@ -273 +274 @@
                 virtual IssmDouble TotalSmb(void)=0;
                 virtual IssmDouble Misfit(int modelenum,int observationenum,int weightsenum)=0;
+                virtual IssmDouble Masscon(IssmDouble* levelset)=0;
                 virtual IssmDouble MisfitArea(int weightsenum)=0;
                 virtual int    VertexConnectivity(int vertexindex)=0;

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

@@ -74 +74 @@
                 Element* GetUpperElement(void);
                 Element* GetBasalElement(void);
+                void        GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues){_error_("not implemented yet");};
                 void   GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
                 IssmDouble GetGroundedPortion(IssmDouble* xyz_list);
@@ -132 +133 @@
                 IssmDouble Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");};
                 IssmDouble MisfitArea(int weightsenum){_error_("not implemented yet");};
+                IssmDouble Masscon(IssmDouble* levelset){_error_("not implemented yet");};
 
                 void   MigrateGroundingLine(IssmDouble* sheet_ungrounding);

issm/trunk-jpl/src/c/classes/Elements/Seg.h ¶

@@ -109 +109 @@
                 int         VelocityInterpolation(void){_error_("not implemented yet");};
                 int         TensorInterpolation(void){_error_("not implemented yet");};
+                void        GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues){_error_("not implemented yet");};
                 void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
                 IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
@@ -157 +158 @@
                 IssmDouble Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");};
                 IssmDouble MisfitArea(int weightsenum){_error_("not implemented yet");};
+                IssmDouble Masscon(IssmDouble* levelset){_error_("not implemented yet");};
 
 #ifdef _HAVE_GIA_

issm/trunk-jpl/src/c/classes/Elements/Tetra.h ¶

@@ -114 +114 @@
                 int         VelocityInterpolation(void);
                 int         TensorInterpolation(void);
+                void        GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues){_error_("not implemented yet");};
                 void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
                 IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
@@ -161 +162 @@
                 IssmDouble Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");};
                 IssmDouble MisfitArea(int weightsenum){_error_("not implemented yet");};
+                IssmDouble Masscon(IssmDouble* levelset){_error_("not implemented yet");};
 
 #ifdef _HAVE_GIA_

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

@@ -666 +666 @@
         return this->element_type;
 
+}
+/*}}}*/
+void        Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
+        
+        /*Computeportion of the element that has a positive levelset*/ 
+
+        bool               negative=true;
+        int                point;
+        const IssmPDouble  epsilon= 1.e-15;
+        IssmDouble         f1,f2;
+
+        /*Be sure that values are not zero*/
+        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+        /*Check that not all nodes are positive or negative*/
+        if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive
+                point=0;
+                f1=1.;
+                f2=1.;
+        }
+        else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative
+                point=0;
+                f1=0.;
+                f2=0.;
+        }
+        else{
+                if(gl[0]*gl[1]*gl[2]<0) negative=false;
+
+                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+                        point=2;
+                        f1=gl[2]/(gl[2]-gl[0]);
+                        f2=gl[2]/(gl[2]-gl[1]);
+                }
+                else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+                        point=0;
+                        f1=gl[0]/(gl[0]-gl[1]);
+                        f2=gl[0]/(gl[0]-gl[2]);
+                }
+                else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+                        point=1;
+                        f1=gl[1]/(gl[1]-gl[2]);
+                        f2=gl[1]/(gl[1]-gl[0]);
+                }
+        }
+        *point1=point;
+        *fraction1=f1;
+        *fraction2=f2;
+        *mainlynegative=negative;
 }
 /*}}}*/
@@ -2716 +2766 @@
 }
 /*}}}*/
+IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/
+
+
+        /*intermediary: */
+        IssmDouble* values=NULL;
+        Input*      thickness_input=NULL;
+        IssmDouble  thickness;
+        IssmDouble  weight;
+        IssmDouble  Jdet;
+        IssmDouble  volume;
+        IssmDouble  rho_ice;
+        IssmDouble* xyz_list=NULL;
+        int         point1;
+        IssmDouble  fraction1,fraction2;
+        bool        mainlynegative=true;
+        
+        /*Output:*/
+        volume=0;
+
+        /* Get node coordinates and dof list: */
+        GetVerticesCoordinates(&xyz_list);
+
+        /*Retrieve inputs required:*/
+        thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
+        
+        /*Retrieve material parameters: */
+        rho_ice=matpar->GetRhoIce();
+
+        /*Retrieve values of the levelset defining the masscon: */
+        values = xNew<IssmDouble>(NUMVERTICES);
+        for(int i=0;i<NUMVERTICES;i++){
+                values[i]=levelset[this->vertices[i]->Sid()];
+        }
+                
+        /*Ok, use the level set values to figure out where we put our gaussian points:*/
+        this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
+        Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
+
+        volume=0;
+
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                this->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                thickness_input->GetInputValue(&thickness, gauss);
+
+                volume+=thickness*gauss->weight*Jdet;
+        }
+
+        /* clean up and Return: */
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(values);
+        delete gauss;
+        return rho_ice*volume;
+}
+/*}}}*/
 
 #ifdef _HAVE_GIA_

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

@@ -69 +69 @@
                 Element*    GetUpperElement(void){_error_("not implemented yet");};
                 Element*    GetBasalElement(void){_error_("not implemented yet");};
+                void        GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues);
                 void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
                 IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
@@ -123 +124 @@
                 IssmDouble MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id);
                 void       ElementResponse(IssmDouble* presponse,int response_enum);
+                IssmDouble Masscon(IssmDouble* levelset);
                 IssmDouble Misfit(int modelenum,int observationenum,int weightsenum);
                 IssmDouble MisfitArea(int weightsenum);

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

@@ -19 +19 @@
 #include "./Massfluxatgate.h"
 #include "./Misfit.h"
+#include "./Masscon.h"
 
 /*Constraints: */

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp ¶

@@ -114 +114 @@
                                 /*}}}*/
                         }
+                        else if (output_definition_enums[i]==MassconEnum){
+                                /*Deal with masscons: {{{*/
+                                
+                                /*masscon variables: */
+                                int          nummasscons;
+                                char**       masscon_name_s             = NULL;    
+                                IssmDouble** masscon_levelset_s           = NULL;
+                                int*         masscon_levelset_M_s    = NULL;
+                                int*         masscon_levelset_N_s    = NULL;
 
+                                /*Fetch name and levelset, etc ... (see src/m/classes/masscon.m): */
+                                iomodel->FetchMultipleData(&masscon_name_s,&nummasscons,MassconNameEnum);
+                                iomodel->FetchMultipleData(&masscon_levelset_s,&masscon_levelset_M_s,&masscon_levelset_N_s,&nummasscons,MassconLevelsetEnum);
+                                for(j=0;j<nummasscons;j++){
+
+                                        /*Create a masscon object: */
+                                        output_definitions->AddObject(new Masscon(masscon_name_s[j],masscon_levelset_s[j],masscon_levelset_M_s[j]));
+
+                                }
+
+                                /*Free ressources:*/
+                                for(j=0;j<nummasscons;j++){
+                                        char* string=NULL;
+                                        IssmDouble* matrix = NULL;
+
+                                        string = masscon_name_s[j];    xDelete<char>(string);
+                                        matrix = masscon_levelset_s[j]; xDelete<IssmDouble>(matrix);
+                                }
+                                xDelete<char*>(masscon_name_s);
+                                xDelete<IssmDouble*>(masscon_levelset_s);
+                                xDelete<int>(masscon_levelset_M_s);
+                                xDelete<int>(masscon_levelset_N_s);
+                                /*}}}*/
+                        }
                         else _error_("output definition enum " << output_definition_enums[i] << "not supported yet!");
                 }

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -491 +491 @@
         ProfilerEnum,
         MatrixParamEnum,
+        MassconEnum,
+        MassconNameEnum,
+        MassconLevelsetEnum,
         NodeSIdEnum,
         VectorParamEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp ¶

@@ -488 +488 @@
                 case ProfilerEnum : return "Profiler";
                 case MatrixParamEnum : return "MatrixParam";
+                case MassconEnum : return "Masscon";
+                case MassconNameEnum : return "MassconName";
+                case MassconLevelsetEnum : return "MassconLevelset";
                 case NodeSIdEnum : return "NodeSId";
                 case VectorParamEnum : return "VectorParam";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp ¶

@@ -497 +497 @@
               else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
               else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
+              else if (strcmp(name,"Masscon")==0) return MassconEnum;
+              else if (strcmp(name,"MassconName")==0) return MassconNameEnum;
+              else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
               else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
               else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
@@ -503 +506 @@
               else if (strcmp(name,"Segment")==0) return SegmentEnum;
               else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
-              else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
-              else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
-              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
+              if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
+              else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
+              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
+              else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
               else if (strcmp(name,"StringParam")==0) return StringParamEnum;
               else if (strcmp(name,"Seg")==0) return SegEnum;
@@ -626 +629 @@
               else if (strcmp(name,"P1")==0) return P1Enum;
               else if (strcmp(name,"P1DG")==0) return P1DGEnum;
-              else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
-              else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
-              else if (strcmp(name,"P2")==0) return P2Enum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
+              if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
+              else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
+              else if (strcmp(name,"P2")==0) return P2Enum;
+              else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
               else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
               else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
@@ -749 +752 @@
               else if (strcmp(name,"OptionCell")==0) return OptionCellEnum;
               else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum;
-              else if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
-              else if (strcmp(name,"Paterson")==0) return PatersonEnum;
-              else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
+              if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
+              else if (strcmp(name,"Paterson")==0) return PatersonEnum;
+              else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
+              else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
               else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum;
               else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;

issm/trunk-jpl/src/m/enum/EnumDefinitions.py ¶

@@ -480 +480 @@
 def ProfilerEnum(): return StringToEnum("Profiler")[0]
 def MatrixParamEnum(): return StringToEnum("MatrixParam")[0]
+def MassconEnum(): return StringToEnum("Masscon")[0]
+def MassconNameEnum(): return StringToEnum("MassconName")[0]
+def MassconLevelsetEnum(): return StringToEnum("MassconLevelset")[0]
 def NodeSIdEnum(): return StringToEnum("NodeSId")[0]
 def VectorParamEnum(): return StringToEnum("VectorParam")[0]

issm/trunk-jpl/src/m/exp/exp_to_levelset.m ¶

@@ -23 +23 @@
         cleanup=1;
         if cleanup,
-                flags=zeros(length(segments),1);
-                for j=1:length(segments),
+                flags=zeros(size(segments,1),1);
+                for j=1:size(segments,1),
                         segment=segments(j,:);
                         x1=segment(1); x2=segment(3); 
@@ -74 +74 @@
         %cleanup: remove 0 length segments: 
         if cleanup,
-                flags=zeros(length(segments),1);
-                for j=1:length(segments),
+                flags=zeros(size(segments,1),1);
+                for j=1:size(segments,1),
                         segment=segments(j,:);
                         x1=segment(1); x2=segment(3); 

issm/trunk-jpl/src/m/plot/applyoptions.m ¶

@@ -427 +427 @@
         end
         hold on,p=plot(x,y,'k.');
-        markersize=getfieldvalue(options,'markersize',5);
+        markersize=getfieldvalue(options,'markersize',12);
         color=getfieldvalue(options,'cloudcolor','k');
         set(p,'Color',color);