linear_reorder

[fwd_model] = linear_reorder(fwd_model,ccw)
Function to reorder local nodes (counter)clockwise per element
Input:  - fwd_model structure
        - ccw = -1 (default) - counter clockwise OR 1 - clockwise   
Output: - fwd_model structure (only .elems changes)
NOTE:Function only for linear triangles, since in this case, identity:
         No. of nodes/element = No. spatial dimensions + 1

0001 function [fwd_model] = linear_reorder(fwd_model,ccw)
0002 %[fwd_model] = linear_reorder(fwd_model,ccw)
0003 %Function to reorder local nodes (counter)clockwise per element
0004 %Input:  - fwd_model structure
0005 %        - ccw = -1 (default) - counter clockwise OR 1 - clockwise
0006 %Output: - fwd_model structure (only .elems changes)
0007 %NOTE:Function only for linear triangles, since in this case, identity:
0008 %         No. of nodes/element = No. spatial dimensions + 1
0009 
0010 % (C) 2011 Michael Crabb. License: GPL version 2 or version 3
0011 % $Id: linear_reorder.m 5112 2015-06-14 13:00:41Z aadler $
0012 
0013 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0014 
0015 if (nargin==1) 
0016     ccw=-1; %Default specify counter-clockwise nodes
0017 end
0018 
0019    fwd_model =     do_reorder(fwd_model, ccw);
0020 %  fwd_model = old_do_reorder(fwd_model, ccw);
0021 end
0022 
0023 
0024 function fmdl = do_reorder(fmdl, ccw)
0025    dim=size( fmdl.nodes,2);
0026    els=num_elems( fmdl );
0027 
0028    if dim==3 && elem_dim( fmdl ) == 2
0029       % For a surface in 3D, we need to walk across elements
0030       fmdl = boundary_align( fmdl );
0031       return
0032    end
0033 
0034      xx= reshape( fmdl.nodes( fmdl.elems, 1 ), els, dim+1);
0035      xx= xx - xx(:,1)*ones(1,dim+1);
0036      yy= reshape( fmdl.nodes( fmdl.elems, 2 ), els, dim+1);
0037      yy= yy - yy(:,1)*ones(1,dim+1);
0038    if dim==2;
0039      vol = xx(:,2).*yy(:,3) - xx(:,3).*yy(:,2);
0040    elseif dim==3
0041      zz= reshape( fmdl.nodes( fmdl.elems, 3 ), els, dim+1);
0042      zz= zz - zz(:,1)*ones(1,dim+1);
0043 
0044      vol = zz(:,4).*( xx(:,2).*yy(:,3) - xx(:,3).*yy(:,2) ) ...
0045          - zz(:,3).*( xx(:,2).*yy(:,4) - xx(:,4).*yy(:,2) ) ...
0046          + zz(:,2).*( xx(:,3).*yy(:,4) - xx(:,4).*yy(:,3) );
0047 
0048 % Looks more elegant, but slower
0049 %    vol = sum(zz(:,[4,3,2]).*xx(:,[2,4,3]).*yy(:,[3,2,4]) ,2) ...
0050 %        - sum(zz(:,[4,3,2]).*xx(:,[3,2,4]).*yy(:,[2,4,3]) ,2);
0051 
0052    else
0053      error('reorder for 2 or 3 dimensions');
0054    end
0055 
0056    ff = find( sign(vol) == ccw );
0057    % reverse first two nodes
0058    fmdl.elems(ff,1:2) = fmdl.elems(ff,[2,1]);
0059 end
0060 
0061 function mdl = boundary_align( mdl );
0062    Ne = num_elems(mdl);
0063    checked = ( 1:Ne )';
0064    tocheck = zeros(Ne,1); tocheck(1) = 1; maxptr = 2;
0065    edges = mdl.elem2edge;
0066 
0067    for i=1:num_elems(mdl)
0068       curr = tocheck(i);
0069       curredges = edges(curr,:);
0070       for j=curredges(:)'
0071          ff = any( j == edges,2 );
0072          ff(curr) = 0;
0073          ff = find(ff); % Should only have 1
0074          if length(ff)>1;
0075             error('Only one elem should share this edge. unexpected error.');
0076          end
0077          if checked(ff) == 0; % we've already done it
0078             continue;
0079          end
0080          nodes_j = mdl.edges(j,:); % Nodes in this edge
0081          a1= find(mdl.elems(curr,:) == nodes_j(1));
0082          a2= find(mdl.elems(curr,:) == nodes_j(2));
0083          curr_o = rem(3+a1-a2,3);
0084          b1= find(mdl.elems(ff  ,:) == nodes_j(1));
0085          b2= find(mdl.elems(ff  ,:) == nodes_j(2));
0086          ff_o = rem(3+b1-b2,3);
0087          if (ff_o ~= curr_o); 
0088             mdl.elems(ff,[b1,b2]) = mdl.elems(ff,[b2,b1]);
0089          end
0090          % Fix orientation
0091          tocheck(maxptr) = ff; maxptr= maxptr+1;     
0092          checked(ff)= 0; % set as done
0093       end
0094    end
0095       
0096 end 
0097 
0098 function fwd_model = old_do_reorder(fwd_model, ccw)
0099    nodecoords = fwd_model.nodes; %Cache coorindates of nodes [nnodesxnodedim]
0100    elementnodes = fwd_model.elems; %Cache matrix of elements [eletotalxelenode]
0101 
0102    eletotal = size(elementnodes,1); %No. of elements
0103    elenode = size(elementnodes,2); %No. of nodes per element
0104    nodedim = size(fwd_model.nodes,2);
0105    midpoint = mean(fwd_model.nodes(unique(fwd_model.elems),:));
0106 
0107    for e=1:eletotal; %Loop over all elements
0108        %Row vector of global nodes [1xelenode]
0109        enodes = elementnodes(e,:); 
0110        %Matrix of nodal positions [elenodexdim] (Linear dimension==elenode-1)
0111        nd = nodecoords(enodes,:);
0112        
0113        % surface meshes need tweaking. Use the midpoint to fit the 3D formula.
0114        % This will not work for non simply-connected surfaces.
0115        if elenode == 3 && nodedim == 3
0116           nd = [nd; midpoint]; 
0117        end
0118        %Calculate area of triangle/volume defined by the elements nodes
0119        %In 2D this is area and in 3D this is volume
0120        area= det([ones(length(nd),1),nd]);
0121        areasign=sign(area); 
0122        
0123        %If sign is (pos) neg swap two nodes (last two will suffice..)
0124        if(areasign == ccw) %Swap last two entries of enodes
0125            temp=enodes(elenode-1);
0126            enodes(elenode-1)=enodes(elenode);
0127            enodes(elenode) = temp;
0128            %elementnodes(e,:)=enodes; %Put back into elementnodes matrix
0129        end
0130        elementnodes(e,:)=enodes; %Put enodes back into elementnodes matrix
0131    end
0132    fwd_model.elems=elementnodes; %Reassign fwd_model.elems
0133 end
0134 
0135 function do_unit_test
0136    for i = 1:10
0137      clear imdl;
0138      switch i,
0139        case 1; imdl = mk_common_model('n3r2',[16,2]);
0140        case 2; imdl = mk_common_model('a2c2',8);
0141        case 3; imdl = mk_common_model('d3cr',[16,2]);
0142        case 4; imdl = mk_common_model('f3cr',[16,2]);
0143      end
0144      if ~exist('imdl'); continue ; end
0145 
0146      fmdl = imdl.fwd_model;
0147      vol = test_linear_reorder( fmdl ); ok = std(sign(vol))==0; % not ok before
0148      t = cputime;
0149      fm0 = linear_reorder(fmdl, -1);
0150      fm1 = linear_reorder(fmdl, 1);
0151      t = cputime - t;
0152      vol0= test_linear_reorder( fm0 );
0153      vol1= test_linear_reorder( fm1 );
0154 
0155      str= sprintf('test%02d(t=%4.2f): OK=%d=>(%d,%d)',i, t, ...
0156              ok, all(vol0>0), all(vol1<0));
0157      unit_test_cmp(str, all(vol0>0) & all(vol1<0), 1);
0158    end
0159 end
0160    
0161 
0162 function vol = test_linear_reorder(fwd_model)
0163 
0164    dim=size(fwd_model.nodes,2); elee=size(fwd_model.elems,1);
0165 
0166    for e=1:elee
0167        b=fwd_model.elems(e,:);  [v]=fwd_model.nodes(b,:);
0168            for i=1:dim
0169                vv1(i,:)=v(i+1,:)-v(1,:);
0170            end
0171        vol(e)=det([vv1]);
0172    end
0173 end
0174