Changeset 17483

Timestamp:

03/19/14 19:15:02 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added tetras method in model that splits mesh3dprism into mesh3dtetra

Location:

issm/trunk-jpl/src/m/classes

Files:

• mesh3dtetras.m (added)
• model.m (modified) (1 diff)

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

@@ -1103 +1103 @@
                 end
                 %}}}
+                function md = tetras(md,varargin) % {{{
+                        %TETRAS - split 3d prismatic mesh into 3 tetrahedrons
+                        %
+                        %   Usage:
+                        %      md=tetra(md)
+
+                        if ~isa(md.mesh,'mesh3dprisms')
+                                error('mesh is not a 3d prismatic mesh');
+                        end
+
+                        %Initialize tetra mesh
+                        md.mesh=mesh3dtetras(md.mesh);
+
+                        %Split in 3 tetras
+                        subelement1 = [1 2 3 5];
+                        subelement2 = [4 6 5 1];
+                        subelement3 = [5 6 3 1];
+                        elements=[md.mesh.elements(:,subelement1);md.mesh.elements(:,subelement2);md.mesh.elements(:,subelement3)];
+                        pos_elements = repmat([1:md.mesh.numberofelements]',3,1);
+
+                        md.mesh.elements=elements;
+                        md.mesh.numberofelements=size(elements,1);
+
+                        %p and q (same deal, except for element that are on the bedrock: )
+                        md.friction.p=md.friction.p(pos_elements);
+                        md.friction.q=md.friction.q(pos_elements);
+
+                        %elementstype
+                        if ~isnan(md.flowequation.element_equation)
+                                oldelements_type=md.flowequation.element_equation;
+                                md.flowequation.element_equation=zeros(number_el3d,1);
+                                md.flowequation.element_equation=md.flowequation.element_equation(pos_elements);
+                        end
+
+                        %connectivity
+                        md.mesh.elementconnectivity=NaN;
+
+                        %materials
+                        md.materials.rheology_n=md.materials.rheology_n(pos_elements);
+
+                        %increase connectivity if less than 25:
+                        if md.mesh.average_vertex_connectivity<=25,
+                                md.mesh.average_vertex_connectivity=100;
+                        end
+                end % }}}
                 function disp(obj) % {{{
                         disp(sprintf('%19s: %-22s -- %s','mesh'            ,['[1x1 ' class(obj.mesh) ']'],'mesh properties'));