Changeset 23805

Timestamp:

03/20/19 15:32:22 (6 years ago)

Author:

tsantos

Message:

BUG: depth average for element fields (P0) is working now

Location:

issm/trunk-jpl/src/m/extrusion

Files:

• DepthAverage.m (modified) (1 diff)
• DepthAverage.py (modified) (1 diff)

issm/trunk-jpl/src/m/extrusion/DepthAverage.m ¶

@@ -25 +25 @@
         vector_average=zeros(md.mesh.numberofelements2d,1);
         for i=1:md.mesh.numberoflayers-1,
-                vector_average=vector_average+project2d(md,vector,i).*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i));
+                vertices_dz = (project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i));
+                elements_dz = mean(vertices_dz(md.mesh.elements2d),2);
+                vector_average = vector_average+project2d(md,vector,i).*elements_dz;
+                %vector_average=vector_average+project2d(md,vector,i).*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i));
         end
-        vector_average=vector_average./project2d(md,md.geometry.thickness,1);
+        vertices_thickness = project2d(md,md.geometry.thickness,1);
+        elements_thickness = mean(vertices_thickness(md.mesh.elements2d),2);
+        vector_average = vector_average./elements_thickness;
+        %vector_average=vector_average./project2d(md,md.geometry.thickness,1);
 
 else

issm/trunk-jpl/src/m/extrusion/DepthAverage.py ¶

@@ -40 +40 @@
         vector_average = np.zeros(md.mesh.numberofelements2d)
         for i in range(1, md.mesh.numberoflayers):
-            vector_average = vector_average + project2d(md, vector, i) * (project2d(md, md.mesh.z, i + 1) - project2d(md, md.mesh.z, i))
-        vector_average = vector_average / project2d(md, md.geometry.thickness, 1)
+            vertices_dz = (project2d(md, md.mesh.z, i + 1) - project2d(md, md.mesh.z, i))
+            elements_dz = vertices_dz.mean(1)
+            vector_average = vector_average + project2d(md, vector, i) * elements_dz
+            #vector_average = vector_average + project2d(md, vector, i) * (project2d(md, md.mesh.z, i + 1) - project2d(md, md.mesh.z, i))
+        vertices_thickness = project2d(md, md.geometry.thickness, 1)
+        elements_thickness = vertices_thickness.mean(1)
+        vector_average = vector_average / elements_thickness
+        #vector_average = vector_average / project2d(md, md.geometry.thickness, 1)
 
     else: