Changeset 3873

Timestamp:

05/21/10 11:04:58 (15 years ago)

Author:

seroussi

Message:

fixed vertical velocity

Location:

issm/trunk/src/m/solutions/jpl

Files:

• ModelUpdateInputsFromVector.m (modified) (2 diffs)
• diagnostic_core.m (modified) (5 diffs)
• diagnostic_core_linear.m (modified) (2 diffs)

issm/trunk/src/m/solutions/jpl/ModelUpdateInputsFromVector.m

@@ -1 @@
-function ModelUpdateInputsFromVector(models, vector, enum, typeenum);
+function models=ModelUpdateInputsFromVector(models, vector, enum, typeenum);
 %MODELUPDATEINPUTSFROMVECTOR - Update inputs using a vector
 %
@@ -26 +26 @@
         if isstruct(model), %there is an analysis_type model
                 [model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters] = UpdateInputsFromVector(model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters,vector,enum, typeenum);
+                models.(field)=model;
         end
 

issm/trunk/src/m/solutions/jpl/diagnostic_core.m

@@ -6 +6 @@
 %
 
-%recover fem models
-m_dh=model.dh;
-m_dv=model.dv;
-m_ds=model.ds;
-m_dhu=model.dhu;
-m_sl=model.sl;
-
 %recover parameters common to all solutions
-verbose=m_dhu.parameters.Verbose;
-dim=m_dh.parameters.Dim;
-ishutter=m_dhu.parameters.IsHutter;
-ismacayealpattyn=m_dh.parameters.IsMacAyealPattyn;
-isstokes=m_ds.parameters.IsStokes;
-numrifts=m_dhu.parameters.NumRifts;
-qmu_analysis=m_dh.parameters.QmuAnalysis;
+verbose=model.dhu.parameters.Verbose;
+dim=model.dh.parameters.Dim;
+ishutter=model.dhu.parameters.IsHutter;
+ismacayealpattyn=model.dh.parameters.IsMacAyealPattyn;
+isstokes=model.ds.parameters.IsStokes;
+numrifts=model.dhu.parameters.NumRifts;
+qmu_analysis=model.dh.parameters.QmuAnalysis;
 
         
 %for qmu analysis, be sure the velocity input we are starting from  is the one in the parameters: 
 if qmu_analysis,
-        ModelUpdateInputsFromVector(model,m_dh.vx,VxEnum,VertexEnum);
-        ModelUpdateInputsFromVector(model,m_dh.vy,VyEnum,VertexEnum);
-        ModelUpdateInputsFromVector(model,m_dh.vz,VzEnum,VertexEnum);
+        model=ModelUpdateInputsFromVector(model,m_dh.vx,VxEnum,VertexEnum);
+        model=ModelUpdateInputsFromVector(model,m_dh.vy,VyEnum,VertexEnum);
+        model=ModelUpdateInputsFromVector(model,m_dh.vz,VzEnum,VertexEnum);
 end
 
 %Compute slopes: 
-[surfaceslopex,surfaceslopey]=slope_core(m_sl,SurfaceAnalysisEnum);
-[bedslopex,bedslopey]=slope_core(m_sl,BedAnalysisEnum);
+[surfaceslopex,surfaceslopey]=slope_core(model.sl,SurfaceAnalysisEnum);
+[bedslopex,bedslopey]=slope_core(model.sl,BedAnalysisEnum);
 
 %Update:
-ModelUpdateInputsFromVector(model,surfaceslopex,SurfaceSlopexEnum,VertexEnum);
-ModelUpdateInputsFromVector(model,surfaceslopey,SurfaceSlopeyEnum,VertexEnum);
-ModelUpdateInputsFromVector(model,bedslopex,BedSlopexEnum,VertexEnum);
-ModelUpdateInputsFromVector(model,bedslopey,BedSlopeyEnum,VertexEnum);
+model=ModelUpdateInputsFromVector(model,surfaceslopex,SurfaceSlopexEnum,VertexEnum);
+model=ModelUpdateInputsFromVector(model,surfaceslopey,SurfaceSlopeyEnum,VertexEnum);
+model=ModelUpdateInputsFromVector(model,bedslopex,BedSlopexEnum,VertexEnum);
+model=ModelUpdateInputsFromVector(model,bedslopey,BedSlopeyEnum,VertexEnum);
 
 if ishutter,
 
         displaystring(verbose,'\n%s',['computing hutter velocities...']);
-        u_g=diagnostic_core_linear(m_dhu,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
+        u_g=diagnostic_core_linear(model.dhu,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
 
         displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']);
-        p_g=ComputePressure(m_dhu.elements,m_dhu.nodes,m_dhu.vertices,m_dhu.loads,m_dhu.materials,m_dhu.parameters,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
+        p_g=ComputePressure(model.dhu.elements,model.dhu.nodes,model.dhu.vertices,model.dhu.loads,model.dhu.materials,model.dhu.parameters,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
 
         displaystring(verbose,'\n%s',['update boundary conditions for macyeal pattyn using hutter results...']);
         if ismacayealpattyn,
-                m_dh.y_g=u_g;
-                [m_dh.ys m_dh.ys0]=Reducevectorgtos(m_dh.y_g,m_dh.nodesets);
+                model.dh.y_g=u_g;
+                [model.dh.ys model.dh.ys0]=Reducevectorgtos(model.dh.y_g,model.dh.nodesets);
         end
 
@@ -59 +52 @@
 
         displaystring(verbose,'\n%s',['computing horizontal velocities...']);
-        [u_g m_dh.loads]=diagnostic_core_nonlinear(m_dh,m_dh.loads,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+        [u_g model.dh.loads]=diagnostic_core_nonlinear(model.dh,model.dh.loads,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 
         if dim==2,
                 displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']);
-                p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+                p_g=ComputePressure(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
         end
 end
@@ -70 +63 @@
 
         displaystring(verbose,'\n%s',['extruding horizontal velocities...']);
-        u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,u_g,'velocity',1);
+        u_g_horiz=FieldExtrude(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,u_g,'velocity',1);
 
-        [vx,vy]=SplitSolutionVector(u_g_horiz,m_dh.parameters.NumberOfNodes,m_dh.parameters.NumberOfDofsPerNode);
-        ModelUpdateInputsFromVector(model,vx,VxEnum,VertexEnum);
-        ModelUpdateInputsFromVector(model,vy,VyEnum,VertexEnum);
+        [vx,vy]=SplitSolutionVector(u_g_horiz,model.dh.parameters.NumberOfNodes,model.dh.parameters.NumberOfDofsPerNode);
+        model=ModelUpdateInputsFromVector(model,vx,VxEnum,VertexEnum);
+        model=ModelUpdateInputsFromVector(model,vy,VyEnum,VertexEnum);
                 
         displaystring(verbose,'\n%s',['computing vertical velocities...']);
-        u_g_vert=diagnostic_core_linear(m_dv,DiagnosticAnalysisEnum(),VertAnalysisEnum());
-        ModelUpdateInputsFromVector(model,u_g_vert,VzEnum,VertexEnum);
+        u_g_vert=diagnostic_core_linear(model.dv,DiagnosticAnalysisEnum(),VertAnalysisEnum());
+        model=ModelUpdateInputsFromVector(model,u_g_vert,VzEnum,VertexEnum);
 
         displaystring(verbose,'\n%s',['computing pressure according to Pattyn...']);
-        p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
-        ModelUpdateInputsFromVector(model,p_g,PressureEnum,VertexEnum);
+        p_g=ComputePressure(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+        model=ModelUpdateInputsFromVector(model,p_g,PressureEnum,VertexEnum);
+        u_g=zeros(model.dh.parameters.NumberOfNodes*3,1); %%%%%%%%%%%%%%%%%%%%%%% NEED TO BE CLEANED 
+        u_g(1:3:end)=vx;
+        u_g(2:3:end)=vy;
+        u_g(3:3:end)=u_g_vert;
         
         if isstokes,
@@ -90 +87 @@
 
                 displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']);
-                m_ds.y_g=zeros(m_ds.nodesets.gsize,1);
-                m_ds.y_g(dofsetgen([1,2],4,m_ds.nodesets.gsize))=u_g;
-                m_ds.y_g(dofsetgen([3],4,m_ds.nodesets.gsize))=u_g_vert;
-                [m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets);
+                model.ds.y_g=zeros(model.ds.nodesets.gsize,1);
+                model.ds.y_g(dofsetgen([1,2],4,model.ds.nodesets.gsize))=u_g;
+                model.ds.y_g(dofsetgen([3],4,model.ds.nodesets.gsize))=u_g_vert;
+                [model.ds.ys model.ds.ys0]=Reducevectorgtos(model.ds.y_g,model.ds.nodesets);
 
                 displaystring(verbose,'\n%s',['computing stokes velocities and pressure ...']);
-                u_g=diagnostic_core_nonlinear(m_ds,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
+                u_g=diagnostic_core_nonlinear(model.ds,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
         
                 %"decondition" pressure
-                p_g=u_g(4:4:end)*m_dh.parameters.stokesreconditioning;
+                p_g=u_g(4:4:end)*model.dh.parameters.stokesreconditioning;
         end
 end
@@ -109 +106 @@
 
 if numrifts,
-        results.riftproperties=OutputRifts(m_dh.loads,m_dh.parameters);
+        results.riftproperties=OutputRifts(model.dh.loads,model.dh.parameters);
 end

issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m

@@ -8 +8 @@
         m.parameters.Kflag=1; m.parameters.Pflag=1;
 
-        %Update inputs in datasets
-        [m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters);
-        
         %system matrices
         [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
@@ -21 +18 @@
         %Reduce load from g size to f size
         [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
-        
+
         %Solve  
         [u_f]=Solver(K_ff,p_f,[],m.parameters);