Changeset 17774

Timestamp:

04/18/14 13:22:29 (11 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixing split of prismatic mesh into tets

File:

• issm/trunk-jpl/src/m/classes/model.m (modified) (4 diffs)

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

@@ -1210 +1210 @@
                                         error('Case not supported'); %not supposed to happen!
                                 end
+                                %Reorder elements to make sure they are direct
+                                for j=1:3
+                                        element = elements(3*(i-1)+j,:);
+                                        matrix = [md.mesh.x(element), md.mesh.y(element), md.mesh.z(element), ones(4,1)];
+                                        if det(matrix)>0,
+                                                elements(3*(i-1)+j,1)=element(2);
+                                                elements(3*(i-1)+j,2)=element(1);
+                                        end
+                                end
                         end
                         %%Split in 3 tetras
@@ -1216 +1225 @@
                         %subelement3 = [5 6 3 1];
                         %elements=[md.mesh.elements(:,subelement1);md.mesh.elements(:,subelement2);md.mesh.elements(:,subelement3)];
-                        disp([num2str(steiner) ' Steiner points had to be included'])
+                        if steiner==0,
+                                disp('No Steiner point required to split prismatic mesh into tets');
+                        else
+                                disp([num2str(steiner) ' Steiner points had to be included'])
+                                error('Steiner point not supported yet');
+                        end
 
                         pos_elements = repmat([1:md.mesh.numberofelements]',3,1);
@@ -1224 +1238 @@
 
                         %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);
+                        if ~isnan(md.friction.p),
+                                md.friction.p=md.friction.p(pos_elements);
+                                md.friction.q=md.friction.q(pos_elements);
+                        end
 
                         %elementstype
@@ -1237 +1253 @@
 
                         %materials
-                        md.materials.rheology_n=md.materials.rheology_n(pos_elements);
+                        if ~isnan(md.materials.rheology_n),
+                                md.materials.rheology_n=md.materials.rheology_n(pos_elements);
+                        end
 
                         %increase connectivity if less than 25: