Changeset 17754

Timestamp:

04/16/14 16:45:24 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: better split in tetras

File:

• issm/trunk-jpl/src/m/classes/model.m (modified) (1 diff)

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

@@ -1144 +1144 @@
                         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)];
+                        %Subdivision from Philipp Furnstahl (http://studierstube.icg.tugraz.at/thesis/fuernstahl_thesis.pdf)
+                        steiner  = 0;
+                        nbv      = md.mesh.numberofvertices;
+                        nbt      = 3*md.mesh.numberofelements;
+                        elements = zeros(nbt,4);
+                        for i=1:md.mesh.numberofelements
+                                v1=md.mesh.elements(i,1); v2=md.mesh.elements(i,2); v3=md.mesh.elements(i,3);
+                                v4=md.mesh.elements(i,4); v5=md.mesh.elements(i,5); v6=md.mesh.elements(i,6);
+                                if(min(v2,v4)<min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)<min(v2,v6)),
+                                        steiner = steiner+1; nbv = nbv+1; nbt = nbt+5; v7 = nbv;
+                                        md.mesh.x=[md.mesh.x; mean(md.mesh.x(md.mesh.elements(i,:)))];
+                                        md.mesh.y=[md.mesh.y; mean(md.mesh.y(md.mesh.elements(i,:)))];
+                                        md.mesh.z=[md.mesh.z; mean(md.mesh.z(md.mesh.elements(i,:)))];
+                                        elements(3*(i-1)+1,:) = [v1 v2 v3 v7];
+                                        elements(3*(i-1)+2,:) = [v1 v2 v4 v7];
+                                        elements(3*(i-1)+3,:) = [v2 v4 v5 v7];
+                                        elements(end+1,:) = [v2 v3 v5 v7];
+                                        elements(end+1,:) = [v3 v5 v6 v7];
+                                        elements(end+1,:) = [v1 v3 v6 v7];
+                                        elements(end+1,:) = [v1 v4 v6 v7];
+                                        elements(end+1,:) = [v4 v5 v6 v7];
+                                elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)>min(v2,v6)),
+                                        elements(3*(i-1)+1,:) = [v1 v2 v4 v6];
+                                        elements(3*(i-1)+2,:) = [v2 v4 v5 v6];
+                                        elements(3*(i-1)+3,:) = [v1 v2 v3 v6];
+                                elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)<min(v2,v6)),
+                                        elements(3*(i-1)+1,:) = [v1 v2 v3 v4];
+                                        elements(3*(i-1)+2,:) = [v2 v3 v4 v5];
+                                        elements(3*(i-1)+3,:) = [v3 v4 v5 v6];
+                                elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)>min(v2,v6)),
+                                        elements(3*(i-1)+1,:) = [v1 v2 v3 v4];
+                                        elements(3*(i-1)+2,:) = [v2 v4 v5 v6];
+                                        elements(3*(i-1)+3,:) = [v2 v3 v4 v6];
+                                elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)>min(v2,v6)),
+                                        elements(3*(i-1)+1,:) = [v1 v4 v5 v6];
+                                        elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
+                                        elements(3*(i-1)+3,:) = [v1 v3 v5 v6];
+                                elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)>min(v2,v6)),
+                                        elements(3*(i-1)+1,:) = [v1 v4 v5 v6];
+                                        elements(3*(i-1)+2,:) = [v1 v2 v5 v6];
+                                        elements(3*(i-1)+3,:) = [v1 v2 v3 v6];
+                                elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)<min(v2,v6)),
+                                        elements(3*(i-1)+1,:) = [v1 v3 v4 v5];
+                                        elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
+                                        elements(3*(i-1)+3,:) = [v3 v4 v5 v6];
+                                elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)<min(v2,v6)),
+                                        %NEW
+                                        elements(3*(i-1)+1,:) = [v1 v5 v6 v4];
+                                        elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
+                                        elements(3*(i-1)+3,:) = [v5 v6 v3 v1];
+                                elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)>min(v2,v6)),
+                                        steiner = steiner+1; nbv = nbv+1; nbt = nbt+5; v7 = nbv;
+                                        md.mesh.x=[md.mesh.x; mean(md.mesh.x(md.mesh.elements(i,:)))];
+                                        md.mesh.y=[md.mesh.y; mean(md.mesh.y(md.mesh.elements(i,:)))];
+                                        md.mesh.z=[md.mesh.z; mean(md.mesh.z(md.mesh.elements(i,:)))];
+                                        elements(3*(i-1)+1,:) = [v1 v2 v3 v7];
+                                        elements(3*(i-1)+2,:) = [v1 v4 v5 v7];
+                                        elements(3*(i-1)+3,:) = [v1 v2 v5 v7];
+                                        elements(end+1,:) = [v2 v5 v6 v7];
+                                        elements(end+1,:) = [v2 v3 v6 v7];
+                                        elements(end+1,:) = [v3 v4 v6 v7];
+                                        elements(end+1,:) = [v1 v3 v4 v7];
+                                        elements(end+1,:) = [v4 v5 v6 v7];
+                                else
+                                        error('Case not supported'); %not supposed to happen!
+                                end
+                        end
+                        %%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)];
+                        disp([num2str(steiner) ' Steiner points had to be included'])
+
                         pos_elements = repmat([1:md.mesh.numberofelements]',3,1);