Changeset 14529

Timestamp:

04/10/13 15:52:49 (12 years ago)

Author:

seroussi

Message:

NEW: added variable forcing load for stokes

Location:

issm/trunk-jpl/src

Files:

• c/EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• c/classes/objects/Elements/Penta.cpp (modified) (3 diffs)
• c/modules/EnumToStringx/EnumToStringx.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (modified) (2 diffs)
• c/modules/StringToEnumx/StringToEnumx.cpp (modified) (5 diffs)
• m/boundaryconditions/SetIceSheetBC.m (modified) (1 diff)
• m/boundaryconditions/SetIceSheetBC.py (modified) (1 diff)
• m/boundaryconditions/SetIceShelfBC.m (modified) (1 diff)
• m/boundaryconditions/SetIceShelfBC.py (modified) (1 diff)
• m/boundaryconditions/SetMarineIceSheetBC.m (modified) (1 diff)
• m/boundaryconditions/SetMarineIceSheetBC.py (modified) (1 diff)
• m/classes/diagnostic.m (modified) (4 diffs)
• m/classes/diagnostic.py (modified) (4 diffs)
• m/classes/model/model.m (modified) (3 diffs)
• m/enum/EnumDefinitions.py (modified) (2 diffs)
• m/enum/LoadingforceXEnum.m (added)
• m/enum/LoadingforceYEnum.m (added)
• m/enum/LoadingforceZEnum.m (added)
• m/enum/MaximumNumberOfEnums.m (modified) (1 diff)

issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h

@@ -63 +63 @@
         DiagnosticVertexPairingEnum,
         DiagnosticViscosityOvershootEnum,
+        LoadingforceXEnum,
+        LoadingforceYEnum,
+        LoadingforceZEnum,
         FlowequationBordermacayealEnum,
         FlowequationBorderpattynEnum,

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

@@ -7665 +7665 @@
         IssmDouble Jdet,viscosity;
         IssmDouble gravity,rho_ice,stokesreconditioning;
+        IssmDouble forcex,forcey,forcez;
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
@@ -7690 +7691 @@
         Input* vy_input=inputs->GetInput(VyEnum);   _assert_(vy_input);
         Input* vz_input=inputs->GetInput(VzEnum);   _assert_(vz_input);
+        Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
+        Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
+        Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
 
         /* Start  looping on the number of gaussian points: */
@@ -7705 +7709 @@
                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
 
+                loadingforcex_input->GetInputValue(&forcex, gauss);
+                loadingforcey_input->GetInputValue(&forcey, gauss);
+                loadingforcez_input->GetInputValue(&forcez, gauss);
+
                 for(i=0;i<NUMVERTICES+1;i++){
                         Pe_gaussian[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l7[i];
+                        Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*l1l7[i];
+                        Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*l1l7[i];
+                        Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*l1l7[i];
                 }
 

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

@@ -68 +68 @@
                 case DiagnosticVertexPairingEnum : return "DiagnosticVertexPairing";
                 case DiagnosticViscosityOvershootEnum : return "DiagnosticViscosityOvershoot";
+                case LoadingforceXEnum : return "LoadingforceX";
+                case LoadingforceYEnum : return "LoadingforceY";
+                case LoadingforceZEnum : return "LoadingforceZ";
                 case FlowequationBordermacayealEnum : return "FlowequationBordermacayeal";
                 case FlowequationBorderpattynEnum : return "FlowequationBorderpattyn";

issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp

@@ -69 +69 @@
         iomodel->FetchDataToInput(elements,VxEnum);
         iomodel->FetchDataToInput(elements,VyEnum);
+        iomodel->FetchDataToInput(elements,LoadingforceXEnum);
+        iomodel->FetchDataToInput(elements,LoadingforceYEnum);
         if(materials_type==MatdamageiceEnum){
                 iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum);
@@ -79 +81 @@
                 iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
                 iomodel->FetchDataToInput(elements,FlowequationBorderstokesEnum);
+                iomodel->FetchDataToInput(elements,LoadingforceZEnum);
         }
 

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

@@ -69 +69 @@
               else if (strcmp(name,"DiagnosticVertexPairing")==0) return DiagnosticVertexPairingEnum;
               else if (strcmp(name,"DiagnosticViscosityOvershoot")==0) return DiagnosticViscosityOvershootEnum;
+              else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
+              else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
+              else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum;
               else if (strcmp(name,"FlowequationBordermacayeal")==0) return FlowequationBordermacayealEnum;
               else if (strcmp(name,"FlowequationBorderpattyn")==0) return FlowequationBorderpattynEnum;
@@ -135 +138 @@
               else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
               else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum;
-              else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum;
-              else if (strcmp(name,"MeshDimension")==0) return MeshDimensionEnum;
-              else if (strcmp(name,"MeshEdges")==0) return MeshEdgesEnum;
          else stage=2;
    }
    if(stage==2){
-              if (strcmp(name,"MeshElementconnectivity")==0) return MeshElementconnectivityEnum;
+              if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum;
+              else if (strcmp(name,"MeshDimension")==0) return MeshDimensionEnum;
+              else if (strcmp(name,"MeshEdges")==0) return MeshEdgesEnum;
+              else if (strcmp(name,"MeshElementconnectivity")==0) return MeshElementconnectivityEnum;
               else if (strcmp(name,"MeshElementonbed")==0) return MeshElementonbedEnum;
               else if (strcmp(name,"MeshElementonsurface")==0) return MeshElementonsurfaceEnum;
@@ -258 +261 @@
               else if (strcmp(name,"PrognosticSolution")==0) return PrognosticSolutionEnum;
               else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
-              else if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum;
-              else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
-              else if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"SurfaceSlopeYAnalysis")==0) return SurfaceSlopeYAnalysisEnum;
+              if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum;
+              else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
+              else if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum;
+              else if (strcmp(name,"SurfaceSlopeYAnalysis")==0) return SurfaceSlopeYAnalysisEnum;
               else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
               else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
@@ -381 +384 @@
               else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
               else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
-              else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
-              else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
-              else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
+              if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
+              else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
+              else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
+              else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
               else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
               else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
@@ -504 +507 @@
               else if (strcmp(name,"PetscOptionsAnalyses")==0) return PetscOptionsAnalysesEnum;
               else if (strcmp(name,"PetscOptionsStrings")==0) return PetscOptionsStringsEnum;
-              else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
-              else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
-              else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"Regular")==0) return RegularEnum;
+              if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
+              else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
+              else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
+              else if (strcmp(name,"Regular")==0) return RegularEnum;
               else if (strcmp(name,"Scaled")==0) return ScaledEnum;
               else if (strcmp(name,"Separate")==0) return SeparateEnum;

issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m

@@ -16 +16 @@
 md.diagnostic.spcvz(pos)=0;
 md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6);
+md.diagnostic.loadingforce=0*ones(md.mesh.numberofvertices,3);
 
 %Dirichlet Values

issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py

@@ -22 +22 @@
         md.diagnostic.spcvz[pos]=0
         md.diagnostic.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
+        md.diagnostic.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
 
         #Dirichlet Values

issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m

@@ -35 +35 @@
 md.diagnostic.spcvz(pos)=0;
 md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6);
+md.diagnostic.loadingforce=0*ones(md.mesh.numberofvertices,3);
 
 %Dirichlet Values

issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py

@@ -39 +39 @@
         md.diagnostic.spcvz[pos]=0
         md.diagnostic.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
+        md.diagnostic.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
 
         #Dirichlet Values

issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m

@@ -42 +42 @@
 md.diagnostic.spcvz(pos)=0;
 md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6);
+md.diagnostic.loadingforce=0*ones(md.mesh.numberofvertices,3);
 
 %Dirichlet Values

issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py

@@ -47 +47 @@
         md.diagnostic.spcvz[pos]=0
         md.diagnostic.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
+        md.diagnostic.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
 
         #Dirichlet Values

issm/trunk-jpl/src/m/classes/diagnostic.m

@@ -23 +23 @@
                 rift_penalty_threshold   = 0;
                 referential              = NaN;
+                loadingforce             = NaN;
                 requested_outputs        = NaN;
         end
@@ -87 +88 @@
                         md = checkfield(md,'diagnostic.maxiter','size',[1 1],'>=',1);
                         md = checkfield(md,'diagnostic.referential','size',[md.mesh.numberofvertices 6]);
+                        md = checkfield(md,'diagnostic.loadingforce','size',[md.mesh.numberofvertices 3]);
                         if ~isempty(md.diagnostic.requested_outputs),
                                 md = checkfield(md,'diagnostic.requested_outputs','size',[NaN 1]);
@@ -144 +146 @@
                         fielddisplay(obj,'stokesreconditioning','multiplier for incompressibility equation. Only for Stokes model');
                         fielddisplay(obj,'referential','local referential');
+                        fielddisplay(obj,'loadingforce','loading force applied on each point');
                         fielddisplay(obj,'requested_outputs','additional outputs requested');
 
@@ -164 +167 @@
                         WriteData(fid,'object',obj,'fieldname','rift_penalty_threshold','format','Integer');
                         WriteData(fid,'object',obj,'fieldname','referential','format','DoubleMat','mattype',1);
+                        %WriteData(fid,'object',obj,'fieldname','loadingforce(:,1)','format','DoubleMat','mattype',1,'enum',LoadingforceXEnum);
+                        %WriteData(fid,'object',obj,'fieldname','loadingforce(:,2)','format','DoubleMat','mattype',1,'enum',LoadingforceYEnum);
+                        %WriteData(fid,'object',obj,'fieldname','loadingforce(:,3)','format','DoubleMat','mattype',1,'enum',LoadingforceZEnum);
+                        WriteData(fid,'data',obj.loadingforce(:,1),'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum);
+                        WriteData(fid,'data',obj.loadingforce(:,2),'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum);
+                        WriteData(fid,'data',obj.loadingforce(:,3),'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum);
                         WriteData(fid,'object',obj,'fieldname','requested_outputs','format','DoubleMat','mattype',3);
 

issm/trunk-jpl/src/m/classes/diagnostic.py

@@ -36 +36 @@
                 self.rift_penalty_threshold   = 0
                 self.referential              = float('NaN')
+                self.loadingforce             = float('NaN')
                 self.requested_outputs        = float('NaN')
 
@@ -74 +75 @@
                 string="%s\n%s"%(string,fielddisplay(self,'stokesreconditioning','multiplier for incompressibility equation. Only for Stokes model'))
                 string="%s\n%s"%(string,fielddisplay(self,'referential','local referential'))
+                string="%s\n%s"%(string,fielddisplay(self,'loadingforce','loading force applied on each point'))
                 string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
 
@@ -136 +138 @@
                 md = checkfield(md,'diagnostic.maxiter','size',[1],'>=',1)
                 md = checkfield(md,'diagnostic.referential','size',[md.mesh.numberofvertices,6])
+                md = checkfield(md,'diagnostic.loadingforce','size',[md.mesh.loadingforce,3])
                 if not md.diagnostic.requested_outputs:
                         md = checkfield(md,'diagnostic.requested_outputs','size',[float('NaN'),1])
@@ -184 +187 @@
                 WriteData(fid,'object',self,'fieldname','rift_penalty_threshold','format','Integer')
                 WriteData(fid,'object',self,'fieldname','referential','format','DoubleMat','mattype',1)
+                WriteData(fid,'data',self.loadingforce[:,0],'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum())
+                WriteData(fid,'data',self.loadingforce[:,1],'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum())
+                WriteData(fid,'data',self.loadingforce[:,2],'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum())
                 WriteData(fid,'object',self,'fieldname','requested_outputs','format','DoubleMat','mattype',3)
 

issm/trunk-jpl/src/m/classes/model/model.m

@@ -172 +172 @@
                         md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers);
                         md.diagnostic.referential=project2d(md,md.diagnostic.referential,md.mesh.numberoflayers);
+                        md.diagnostic.loadingforce=project2d(md,md.diagnostic.loadingforce,md.mesh.numberoflayers);
                         md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers);
                         md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
@@ -718 +719 @@
                         md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node');
                         md.diagnostic.referential=project3d(md,'vector',md.diagnostic.referential,'type','node');
+                        md.diagnostic.loadingforce=project3d(md,'vector',md.diagnostic.loadingforce,'type','node');
 
                         %in 3d, pressureload: [node1 node2 node3 node4 element]
@@ -1055 +1057 @@
                         end
 
+                        if ~isfield(structmd,'loadingforce');
+                                md.diagnostic.loadingforce=0*ones(md.mesh.numberofvertices,3);
+                        end
                 end% }}}
                 function md = setdefaultparameters(md) % {{{

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

@@ -519 +519 @@
         return StringToEnum('DiagnosticViscosityOvershoot')[0]
 
+def LoadingforceXEnum():
+        """
+        LOADINGFORCEXENUM - Enum of LoadingforceX
+
+           Usage:
+              macro=LoadingforceXEnum()
+        """
+
+        return StringToEnum('LoadingforceX')[0]
+
+def LoadingforceYEnum():
+        """
+        LOADINGFORCEYENUM - Enum of LoadingforceY
+
+           Usage:
+              macro=LoadingforceYEnum()
+        """
+
+        return StringToEnum('LoadingforceY')[0]
+
+def LoadingforceZEnum():
+        """
+        LOADINGFORCEZENUM - Enum of LoadingforceZ
+
+           Usage:
+              macro=LoadingforceZEnum()
+        """
+
+        return StringToEnum('LoadingforceZ')[0]
+
 def FlowequationBordermacayealEnum():
         """
@@ -5007 +5037 @@
         """
 
-        return 499
-
+        return 502
+

issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m

@@ -9 +9 @@
 %      macro=MaximumNumberOfEnums()
 
-macro=499;
+macro=502;