Changeset 21417

Timestamp:

11/23/16 09:04:15 (8 years ago)

Author:

Mathieu Morlighem

Message:

CHG: make sure that 3d mesh follows geometry, and added collapse of some other quantities

Location:

issm/trunk-jpl/src/m/classes

Files:

• mesh3dprisms.m (modified) (1 diff)
• model.m (modified) (3 diffs)

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

@@ -104 +104 @@
                         md = checkfield(md,'fieldname','mesh.vertexonbase','size',[md.mesh.numberofvertices 1],'values',[0 1]);
                         md = checkfield(md,'fieldname','mesh.vertexonsurface','size',[md.mesh.numberofvertices 1],'values',[0 1]);
+                        md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',24,'message','''mesh.average_vertex_connectivity'' should be at least 24 in 3d');
+
+                        %Check that mesh follows the geometry
                         md = checkfield(md,'fieldname','mesh.z','>=',md.geometry.base-10^-10,'message','''mesh.z'' lower than bedrock');
                         md = checkfield(md,'fieldname','mesh.z','<=',md.geometry.surface+10^-10,'message','''mesh.z'' higher than surface elevation');
-                        md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',24,'message','''mesh.average_vertex_connectivity'' should be at least 24 in 3d');
+                        if any(max(abs(project2d(md,md.mesh.z,1)-project2d(md,md.geometry.base,1)))>1e-11),
+                                md = checkmessage(md,'md.mesh.z is not consistent with md.geometry.base, you changed the geometry after extrusion');
+                        end
+                        if any(max(abs(project2d(md,md.mesh.z,md.mesh.numberoflayers)-project2d(md,md.geometry.surface,1)))>1e-11),
+                                md = checkmessage(md,'md.mesh.z is not consistent with md.geometry.surface, you changed the geometry after extrusion !!');
+                        end
+                        if any(max(abs(project2d(md,md.mesh.z,md.mesh.numberoflayers)-project2d(md,md.mesh.z,1) - project2d(md,md.geometry.thickness,1)))>1e-11),
+                                md = checkmessage(md,'md.mesh.z is not consistent with md.geometry.thickness, you changed the geometry after extrusion !!');
+                        end
                 end % }}}
                 function disp(self) % {{{

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

@@ -223 +223 @@
                         if ~isnan(md.initialization.epl_head),md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);end;
                         if ~isnan(md.initialization.epl_thickness),md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);end;
+                        if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project2d(md,md.initialization.waterfraction,1);end;
+                        if ~isnan(md.initialization.watercolumn),md.initialization.watercolumn=project2d(md,md.initialization.watercolumn,1);end;
                         %gia
                         if ~isnan(md.gia.mantle_viscosity), md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1); end
@@ -244 +246 @@
                         md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
                         if ~isnan(md.damage.spcdamage), md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers); end
+                        if numel(md.levelset.spclevelset)>1, md.levelset.spclevelset=project2d(md,md.levelset.spclevelset,md.mesh.numberoflayers); end
                         md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
 
@@ -270 +273 @@
                                 md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1); 
                         end
-                        if ~isnan(md.basalforcings.floatingice_melting_rate),
+                        if isprop(md.basalforcings,'floatingice_melting_rate') & ~isnan(md.basalforcings.floatingice_melting_rate),
                                 md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1); 
                         end
                         md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux
+
+                        if isprop(md.calving,'coeff') & ~isnan(md.calving.coeff),
+                                md.calving.coeff=project2d(md,md.calving.coeff,1); 
+                        end
+                        if isprop(md.calving,'meltingrate') & ~isnan(md.calving.meltingrate),
+                                md.calving.meltingrate=project2d(md,md.calving.meltingrate,1); 
+                        end
 
                         %update of connectivity matrix