Changeset 26161

Timestamp:

03/30/21 13:43:45 (4 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added super cool function to refine 4x a mesh automatically

File:

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

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

@@ -818 +818 @@
                         md2.mesh.extractedvertices=pos_node;
                         md2.mesh.extractedelements=pos_elem;
+                end % }}}
+                function md2 = refine(md) % {{{
+                        %refine - split all triangles into 3 to refine the mesh everywhere
+                        %
+                        %   This function only works for 2d triangle meshes
+                        %
+                        %   Usage:
+                        %      md2=refine(md);
+                        %
+                        %   See also: EXTRUDE, COLLAPSE, EXTRACT
+
+                        %Check incoming 
+                        if ~strcmp(elementtype(md.mesh),'Tria')
+                                error('not supported for 3d meshes');
+                        end
+
+                        %copy model
+                        md2=md;
+
+                        disp('Getting edges');
+                        %initialization of some variables
+                        nbe   = md.mesh.numberofelements;
+                        nbv   = md.mesh.numberofvertices;
+                        index = md.mesh.elements;
+                        elementslist=1:nbe;
+                        %1: list of edges
+                        edges=[index(:,[1,2]); index(:,[2,3]); index(:,[3,1])];
+                        %2: find unique edges
+                        [edges,I,J]=unique(sort(edges,2),'rows');
+                        %3: unique edge numbers
+                        vec=J;
+                        %4: unique edges numbers in each triangle (2 triangles sharing the same edge will have the same edge number)
+                        edges_tria=[vec(elementslist+nbe) vec(elementslist+2*nbe) vec(elementslist)];
+
+                        % We divide each element as follows
+                        %
+                        %                   e2
+                        %    n1 ------------+------------ n3
+                        %       \          / \          /
+                        %        \    1   /   \   3    /
+                        %         \      /     \      /
+                        %          \    /   2   \    /
+                        %           \  /         \  /
+                        %         e3 +____________\/ e1
+                        %             \           /
+                        %              \         /
+                        %               \   4   /
+                        %                \     /
+                        %                 \   /
+                        %                   n2
+
+                        %Create new coordinates
+                        disp('Remeshing...');
+                        x_edges = 0.5*(md.mesh.x(edges(:,1)) + md.mesh.x(edges(:,2)));
+                        y_edges = 0.5*(md.mesh.y(edges(:,1)) + md.mesh.y(edges(:,2)));
+                        md2.mesh.x = [md2.mesh.x;x_edges];
+                        md2.mesh.y = [md2.mesh.y;y_edges];
+                        md2.mesh.elements = [...
+                                index(:,1)          nbv+edges_tria(:,3) nbv+edges_tria(:,2);...
+                                nbv+edges_tria(:,2) nbv+edges_tria(:,3) nbv+edges_tria(:,1);...
+                                nbv+edges_tria(:,2) nbv+edges_tria(:,1) index(:,3);...
+                                nbv+edges_tria(:,3) index(:,2)          nbv+edges_tria(:,1)];
+                        md2.mesh.numberofelements = 4*nbe;
+                        md2.mesh.numberofvertices = nbv + size(edges,1);
+                        disp(['   Old number of elements: ' num2str(nbe)]);
+                        disp(['   New number of elements: ' num2str(4*nbe)]);
+
+                        disp('Interpolate all fields');
+                        numberofvertices1 = md.mesh.numberofvertices;
+                        numberofelements1 = md.mesh.numberofelements;
+                        nbv2 = md2.mesh.numberofvertices;
+                        %Create transformation vectors
+                        nbedges = size(edges,1);
+                        Pelem = sparse(1:4*nbe,repmat([1:nbe],1,4),ones(4*nbe,1),4*nbe,nbe);
+                        Pnode = sparse([1:nbv,repmat([nbv+1:nbv+nbedges],1,2)],[1:nbv edges(:)'],[ones(nbv,1);1/2*ones(2*nbedges,1)],md2.mesh.numberofvertices,nbv);
+                        %loop over model fields
+                        model_fields=setxor(fields(md),{'mesh'});
+                        %remove mesh from this field
+                        for i=1:length(model_fields),
+                                %get field
+                                field=md.(model_fields{i});
+                                fieldsize=size(field);
+                                if isobject(field), %recursive call
+                                        object_fields=fields(md.(model_fields{i}));
+                                        for j=1:length(object_fields),
+                                                %get field
+                                                field=md.(model_fields{i}).(object_fields{j});
+                                                fieldsize=size(field);
+                                                %size = number of nodes * n
+                                                if fieldsize(1)==numberofvertices1
+                                                        md2.(model_fields{i}).(object_fields{j})=Pnode*field;
+                                                elseif (fieldsize(1)==numberofvertices1+1)
+                                                        md2.(model_fields{i}).(object_fields{j})=[Pnode*field(1:end-1,:); field(end,:)];
+                                                        %size = number of elements * n
+                                                elseif fieldsize(1)==numberofelements1
+                                                        md2.(model_fields{i}).(object_fields{j})=Pelem*field;
+                                                elseif (fieldsize(1)==numberofelements1+1)
+                                                        md2.(model_fields{i}).(object_fields{j})=[Pelem*field(1:end-1,:); field(end,:)];
+                                                end
+                                        end
+                                else
+                                        %size = number of nodes * n
+                                        if fieldsize(1)==numberofvertices1
+                                                md2.(model_fields{i})=Pnode*field;
+                                        elseif (fieldsize(1)==numberofvertices1+1)
+                                                md2.(model_fields{i})=[Pnode*field(1:end-1,:); field(end,:)];
+                                                %size = number of elements * n
+                                        elseif fieldsize(1)==numberofelements1
+                                                md2.(model_fields{i})=Pelem*field;
+                                        elseif (fieldsize(1)==numberofelements1+1)
+                                                md2.(model_fields{i})=[Pelem*field(1:end-1,:); field(end,:)];
+                                        end
+                                end
+                        end
+
                 end % }}}
                 function md = extrude(md,varargin) % {{{