tic %Start timer % Inputs % Cantilever beam initialisation inputVariable = [97.61832484 149.7719077 365.1375156 223.4120554 53.6675781 358.4963287 190.4325727 993.5017445 458.1047723 104.1558069 1263.632702 157.1124904 1856.558758 484.0095117 86.65508219 2278.803746 110.8489292 2989.997407 494.4067854 99.17563938 2.903465949 77.84350147 1285.693456 94.20420511 175.5631069 1047.479593 432.1372272 1225.784109 178.5534327 54.16930538 1233.155207 39.12020394 2310.260023 77.5680831 112.5339576 1944.043104 445.809374 2190.869333 112.3500192 55.5377348 0.818624205 922.6232978 1323.079043 902.6336413 176.7152959 1047.523507 515.8261431 1272.838431 819.9036931 57.18174301 1195.373141 962.5920799 2247.169066 934.1637617 103.1597299 1895.849529 530.0174616 2160.376388 897.5392847 56.11519557 107.4683731 831.9824568 270.5886823 784.0666087 23.4602835 236.9764153 852.3365967 1055.797134 465.0637232 103.3865752 1306.811322 842.1776914 1950.598061 456.3307501 91.55734746 2229.68268 906.8339149 2993.322836 508.9093891 99.43492064]; N = 16; %Number of features % % MBB beam initialisation % inputVariable = [-3.965747491 109.6268401 1151.236636 99.98034054 215.3900746 1163.238955 36.83021514 2954.854123 50.89371887 88.68074753 90.3745719 896.3311845 962.4840866 181.2865674 150.3156322 47.85777149 931.9167069 2245.057279 882.9774955 180.8406982 1087.894813 181.1343394 1988.786778 833.8556109 137.4338256 2207.253945 822.7514084 2985.965954 51.21906711 126.1429372]; % N = 6; featureVars = 5;%Number of design variables for each feature % Constants % Regularised Heaviside parameters alpha = 0; epsilon = 0.3; % Feature construction parameters p = 6; %TDF penalty parameter. aggregationSelection = 'MAX'; % Design domain properties x0 = 0; y0 = 0; xDim = 3000; %Domain x dimension yDim = 1000; %Domain y dimension elSize = 5; %Mesh element size % Threshold Parameters t1Lim = 1; %Threshold for collinear features with endpoints in close proximity parentSnapThreshold = 0; %Either geometric (e.g. 150, 50) or parametric (e.g. 0.0, 0.1, 0.2) childEndpointSelection = 0.5; snapType = 'GEOM'; %Either 'GEOM' or 'PARAM' widthScale = 1; %Snap for close endpoint when within the a multiple of the parent feature width % Visualisation showFeatureSpines = true; showFeatureIds = false; % Create mesh [xElCount,yElCount,nodeCoords,elementNodes,x,y] = CreateMesh(x0,y0,xDim,yDim,elSize); % Variable initialisation updatedFeatureVariableValues = inputVariable(:); masterFeatureIdList = 1:N; N_update = N+1; while N < N_update % Variable initialisation deletedList = []; N_update = N; restartCommand = 0; attemptedConnectionInfo = cell(N,1); % records childConnectionInfoEndA = cell(N,1); childConnectionInfoEndB = cell(N,1); childEndpointACounter = ones(N,1); childEndpointBCounter = ones(N,1); parentConnectionInfoEndA = cell(N,1); parentConnectionInfoEndB = cell(N,1); parentConnectionInfoSpine = cell(N,1); parentEndpointACounter = ones(N,1); parentEndpointBCounter = ones(N,1); parentSpineCounter = ones(N,1); featureIdList = 1:N; newIndependentVariableList = updatedFeatureVariableValues(:); [~,~,~,H_features]=FeatureConstruction(updatedFeatureVariableValues,featureVars,N,nodeCoords,p,alpha,epsilon,aggregationSelection,xElCount,yElCount); [featureSize,uniqueSize,intersectingFeatures,intersectingFeatureSize] = FeatureComparison(H_features,alpha,N,elementNodes,elSize); [~,I] = sort((featureSize/max(featureSize)),'descend'); for k = 1:N parentFeature = featureIdList(I(k)); attemptedConnectionInfo{k,1} = parentFeature; parentFeatureVars = updatedFeatureVariableValues(featureVars*parentFeature-featureVars+1: featureVars*parentFeature); childFeatureList = intersectingFeatures{parentFeature}; [childFeatureList,J] = sort(childFeatureList,'descend'); childFeatureList = childFeatureList(J); if isempty(childFeatureList) attemptedConnectionInfo{k,2} = 'No intersecting features'; %This is unlikely! else for i = 1:length(childFeatureList) childFeature = childFeatureList(i); attemptedConnectionInfo{k,1+i} = childFeature; childFeatureVars = updatedFeatureVariableValues(featureVars*childFeature-featureVars+1: featureVars*childFeature); % Remove the parent feature from the child feature's list of intersecting features to prevent duplicate connections index = find(intersectingFeatures{childFeature} == parentFeature); if isempty(index) %log error. Child intersects parent but parent does not intersect child else intersectingFeatures{childFeature}(index) = []; intersectingFeatureSize{childFeature}(index) = []; end % Calculate intersection positions [tParent,tChild,tLenParent,parentLength,tLenChild,childLength] = IntersectionCalc(parentFeatureVars,childFeatureVars); if abs(tParent-0.5)>t1Lim %This should be updated to check for collinear %check end point proximity % 4 options. EndA-EndA, EndA-EndB, EndB-EndA, EndB-EndB dist = zeros(1,4); dist(1) = ((parentFeatureVars(1)-childFeatureVars(1))^2+(parentFeatureVars(2)-childFeatureVars(2))^2)^0.5; dist(2) = ((parentFeatureVars(1)-childFeatureVars(3))^2+(parentFeatureVars(2)-childFeatureVars(4))^2)^0.5; dist(3) = ((parentFeatureVars(3)-childFeatureVars(1))^2+(parentFeatureVars(4)-childFeatureVars(2))^2)^0.5; dist(4) = ((parentFeatureVars(3)-childFeatureVars(3))^2+(parentFeatureVars(4)-childFeatureVars(4))^2)^0.5; [dist,removedDesignVariableLocations] = min(dist); %loc records the endpoint connection option if dist < parentFeatureVars(5)*widthScale % use loc to select which points to snap if removedDesignVariableLocations == 1 && childEndpointACounter(childFeature) == 1 %snap EndA_EndA updatedFeatureVariableValues(featureVars*childFeature-featureVars+1: featureVars*childFeature-featureVars+2) = parentFeatureVars(1:2); newIndependentVariableList(featureVars*childFeature-featureVars+1: featureVars*childFeature-featureVars+2) = nan; childConnectionInfoEndA{childFeature,1} = [parentFeature,1]; parentConnectionInfoEndA{parentFeature,parentEndpointACounter(parentFeature)} = [childFeature,1]; parentEndpointACounter(parentFeature) = parentEndpointACounter(parentFeature)+1; childEndpointACounter(childFeature) = 2; elseif removedDesignVariableLocations == 2 && childEndpointBCounter(childFeature) == 1 %snap EndA_EndB updatedFeatureVariableValues(featureVars*childFeature-featureVars+3: featureVars*childFeature-featureVars+4) = parentFeatureVars(1:2); newIndependentVariableList(featureVars*childFeature-featureVars+3: featureVars*childFeature-featureVars+4) = nan; childConnectionInfoEndB{childFeature,1} = [parentFeature,1]; parentConnectionInfoEndA{parentFeature,parentEndpointACounter(parentFeature)} = [childFeature,2]; parentEndpointACounter(parentFeature) = parentEndpointACounter(parentFeature)+1; childEndpointBCounter(childFeature) = 2; elseif removedDesignVariableLocations == 3 && childEndpointACounter(childFeature) == 1 %snap EndB_EndA updatedFeatureVariableValues(featureVars*childFeature-featureVars+1: featureVars*childFeature-featureVars+2) = parentFeatureVars(3:4); newIndependentVariableList(featureVars*childFeature-featureVars+1: featureVars*childFeature-featureVars+2) = nan; childConnectionInfoEndA{childFeature,1} = [parentFeature,2]; parentConnectionInfoEndB{parentFeature,parentEndpointBCounter(parentFeature)} = [childFeature,1]; parentEndpointBCounter(parentFeature) = parentEndpointBCounter(parentFeature)+1; childEndpointACounter(childFeature) = 2; elseif removedDesignVariableLocations == 4 && childEndpointBCounter(childFeature) == 1 %snap EndB_EndB updatedFeatureVariableValues(featureVars*childFeature-featureVars+3: featureVars*childFeature-featureVars+4) = parentFeatureVars(3:4); newIndependentVariableList(featureVars*childFeature-featureVars+3: featureVars*childFeature-featureVars+4) = nan; childConnectionInfoEndB{childFeature,1} = [parentFeature,2]; childEndpointBCounter(childFeature) = 2; parentConnectionInfoEndB{parentFeature,parentEndpointBCounter(parentFeature)} = [childFeature,2]; parentEndpointBCounter(parentFeature) = parentEndpointBCounter(parentFeature)+1; end else if mod(removedDesignVariableLocations,2) == 0 childConnectionInfoEndB{childFeature,1} = 'No suitable connection'; childEndpointBCounter(childFeature) = 2; else childConnectionInfoEndA{childFeature,1} = 'No suitable connection'; childEndpointACounter(childFeature) = 2; end end % elseif abs(t1) < snap_proximity elseif abs(tLenParent) < parentSnapThreshold if tChild < childEndpointSelection && childEndpointACounter(childFeature) == 1 %snap end A to end A of master updatedFeatureVariableValues(featureVars*childFeature-featureVars+1: featureVars*childFeature-featureVars+2) = parentFeatureVars(1:2); newIndependentVariableList(featureVars*childFeature-featureVars+1: featureVars*childFeature-featureVars+2) = nan; childConnectionInfoEndA{childFeature,1} = [parentFeature,1]; childEndpointACounter(childFeature) = 2; parentConnectionInfoEndA{parentFeature,parentEndpointACounter(parentFeature)} = [childFeature,1]; parentEndpointACounter(parentFeature) = parentEndpointACounter(parentFeature)+1; elseif tChild > (1-childEndpointSelection) && childEndpointBCounter(childFeature) == 1 %snap end B to end A of master updatedFeatureVariableValues(featureVars*childFeature-featureVars+3: featureVars*childFeature-featureVars+4) = parentFeatureVars(1:2); newIndependentVariableList(featureVars*childFeature-featureVars+3: featureVars*childFeature-featureVars+4) = nan; childConnectionInfoEndB{childFeature,1} = [parentFeature,1]; childEndpointBCounter(childFeature) = 2; parentConnectionInfoEndA{parentFeature,parentEndpointACounter(parentFeature)} = [childFeature,2]; parentEndpointACounter(parentFeature) = parentEndpointACounter(parentFeature)+1; end % elseif abs(1-abs(t1)) < snap_proximity elseif abs(parentLength-tLenParent) < parentSnapThreshold if tChild < childEndpointSelection && childEndpointACounter(childFeature) == 1 %snap end A to end B of master updatedFeatureVariableValues(featureVars*childFeature-featureVars+1: featureVars*childFeature-featureVars+2) = parentFeatureVars(3:4); newIndependentVariableList(featureVars*childFeature-featureVars+1: featureVars*childFeature-featureVars+2) = nan; childConnectionInfoEndA{childFeature,1} = [parentFeature,2]; childEndpointACounter(childFeature) = 2; parentConnectionInfoEndB{parentFeature,parentEndpointBCounter(parentFeature)} = [childFeature,1]; parentEndpointBCounter(parentFeature) = parentEndpointBCounter(parentFeature)+1; elseif tChild > (1-childEndpointSelection) && childEndpointBCounter(childFeature) == 1 %snap end B to end B of master updatedFeatureVariableValues(featureVars*childFeature-featureVars+3: featureVars*childFeature-featureVars+4) = parentFeatureVars(3:4); newIndependentVariableList(featureVars*childFeature-featureVars+3: featureVars*childFeature-featureVars+4) = nan; childConnectionInfoEndB{childFeature,1} = [parentFeature,2]; childEndpointBCounter(childFeature) = 2; parentConnectionInfoEndB{parentFeature,parentEndpointBCounter(parentFeature)} = [childFeature,2]; parentEndpointBCounter(parentFeature) = parentEndpointBCounter(parentFeature)+1; end else %snap to spine if tChild < childEndpointSelection && childEndpointACounter(childFeature) == 1 % Snap child endpoint A to parent spine updatedFeatureVariableValues(featureVars*childFeature-featureVars+1: featureVars*childFeature-featureVars+2) = childFeatureVars(1:2) + (childFeatureVars(3:4)-childFeatureVars(1:2))*tChild; newIndependentVariableList(featureVars*childFeature-featureVars+1) = tLenParent; newIndependentVariableList(featureVars*childFeature-featureVars+2) = nan; childConnectionInfoEndA{childFeature,1} = [parentFeature,3,tLenParent]; childEndpointACounter(childFeature) = 2; parentConnectionInfoSpine{parentFeature,parentSpineCounter(parentFeature)} = [childFeature,1,tLenParent]; parentSpineCounter(parentFeature) = parentSpineCounter(parentFeature)+1; elseif tChild > (1-childEndpointSelection) && childEndpointBCounter(childFeature) == 1 % Snap child endpoint B to parent spine updatedFeatureVariableValues(featureVars*childFeature-featureVars+3: featureVars*childFeature-featureVars+4) = childFeatureVars(1:2) + (childFeatureVars(3:4)-childFeatureVars(1:2))*tChild; newIndependentVariableList(featureVars*childFeature-featureVars+3) = tLenParent; newIndependentVariableList(featureVars*childFeature-featureVars+4) = nan; childConnectionInfoEndB{childFeature,1} = [parentFeature,3,tLenParent]; childEndpointBCounter(childFeature) = 2; parentConnectionInfoSpine{parentFeature,parentSpineCounter(parentFeature)} = [childFeature,2,tLenParent]; parentSpineCounter(parentFeature) = parentSpineCounter(parentFeature)+1; end end % Check if the snapping has merged features (may happend when almost collinear) [~,~,~,H_connectedFeatures]=FeatureConstruction(updatedFeatureVariableValues,featureVars,N,nodeCoords,p,alpha,epsilon,aggregationSelection,xElCount,yElCount); [featureSize,uniqueSize,~,~] = FeatureComparison(H_connectedFeatures,alpha,N,elementNodes,elSize); uniqueRatioConnectedFeatures = uniqueSize./featureSize; if uniqueRatioConnectedFeatures(childFeature) == 0 deletedList(end+1) = featureIdList(childFeature); % Restart to prevent other features connecting to the tested feature. restartCommand = 1; end end end end if restartCommand == 1 for i = 1:length(deletedList) N = N - 1; updatedFeatureVariableValues(featureVars*deletedList(i)-featureVars+1: featureVars*deletedList(i)) = []; newIndependentVariableList(featureVars*deletedList(i)-featureVars+1: featureVars*deletedList(i)) = []; end % break end end combinedParentConnections = {parentConnectionInfoEndA;parentConnectionInfoEndB}; combinedChildConnections = {childConnectionInfoEndA;childConnectionInfoEndB}; combinedParentCounter = {parentEndpointACounter;parentEndpointBCounter}; % Find the global parent. A child feature may otherwise be connected to a parent that has a super parent for i = 1:N checkChild = childConnectionInfoEndA{i}; if isempty(checkChild) % Not necessary if the child has no parent elseif checkChild(2) == 3 % Not necessary for spine connections else checkIfParent = parentConnectionInfoEndA{i}; if isempty(checkIfParent) % The child end A is not used as a parent else for j = 1:size(parentConnectionInfoEndA,2) if isempty(parentConnectionInfoEndA{i,j}) % Do nothing, because, % The global cell width does not mean each cell has a value else % If the temp location in "connection" has %associated dependents in "sens_connection", these %should be moved to the same master var as in the %"connection" array. combinedParentConnections{checkChild(2),1}{checkChild(1),combinedParentCounter{checkChild(2),1}(checkChild(1))} = parentConnectionInfoEndA{i,j}; combinedParentCounter{checkChild(2),1}(checkChild(1)) = combinedParentCounter{checkChild(2),1}(checkChild(1)) + 1; update_con_info_loc = parentConnectionInfoEndA{i,j}; combinedChildConnections{update_con_info_loc(2),1}{update_con_info_loc(1),1} = checkChild; parentConnectionInfoEndA{i,j} = []; combinedParentConnections{1,1}{i,j} = []; end end end end checkChild = childConnectionInfoEndB{i}; if isempty(checkChild) % Not necessary if the child has no parent elseif checkChild(2) == 3 % Not necessary for spine connections else checkIfParent = parentConnectionInfoEndB{i}; if isempty(checkIfParent) % The child end B is not used as a parent else for j = 1:size(parentConnectionInfoEndB,2) if isempty(parentConnectionInfoEndB{i,j}) % Do nothing, because, % The global cell width does not mean each cell has a value else combinedParentConnections{checkChild(2),1}{checkChild(1),combinedParentCounter{checkChild(2),1}(checkChild(1))} = parentConnectionInfoEndB{i,j}; combinedParentCounter{checkChild(2),1}(checkChild(1)) = combinedParentCounter{checkChild(2),1}(checkChild(1)) + 1; update_con_info_loc = parentConnectionInfoEndB{i,j}; combinedChildConnections{update_con_info_loc(2),1}{update_con_info_loc(1),1} = checkChild; parentConnectionInfoEndB{i,j} = []; combinedParentConnections{2,1}{i,j} = []; end end end end end removedDesignVariableLocations = isnan(newIndependentVariableList); designVariablesList = newIndependentVariableList(~removedDesignVariableLocations)'; % without nan designVariablesTable = reshape(newIndependentVariableList,5,N)'; % reshape xy00 countDesignVariablesUpdated = sum(~removedDesignVariableLocations); tableOfDesignVaribleIds = newIndependentVariableList; tableOfDesignVaribleIds(~removedDesignVariableLocations) = 1:countDesignVariablesUpdated; tableOfDesignVaribleIds(removedDesignVariableLocations) = 0; tableOfDesignVaribleIds = reshape(tableOfDesignVaribleIds,5,N)'; % reshaped array of design variable numbers allFeatureVariables = updatedFeatureVariableValues'; % for use as master variable for constructing components featureWidthLocations = zeros(1,countDesignVariablesUpdated); featureWidthLocations(tableOfDesignVaribleIds(:,5)) = 1; %Visualise updated geometry with feature spine [Phi,~,~,~]=FeatureConstruction(updatedFeatureVariableValues,featureVars,N,nodeCoords,p,alpha,epsilon,aggregationSelection,xElCount,yElCount); FeaturePlot(Phi,N,alpha,epsilon,xDim,yDim,x,y,xElCount,yElCount,showFeatureSpines,allFeatureVariables,featureVars,showFeatureIds) toc %End timer function [tParent,tChild,tLenParent,parentLength,tLenChild,childLength] = IntersectionCalc(parentFeatureVars,childFeatureVars) D = (parentFeatureVars(3)-parentFeatureVars(1))/(parentFeatureVars(4)-parentFeatureVars(2)); D2 = (childFeatureVars(3)-childFeatureVars(1))/(childFeatureVars(4)-childFeatureVars(2)); if ~isinf(D) %Where the child feature is intersected by the parent feature tChild = (childFeatureVars(1)-parentFeatureVars(1) + D*(parentFeatureVars(2)-childFeatureVars(2)))/(D*(childFeatureVars(4)-childFeatureVars(2))-childFeatureVars(3)+childFeatureVars(1)); %Where the parent feature is intersected by the child feature tParent = (childFeatureVars(1)+(childFeatureVars(3)-childFeatureVars(1))*tChild-parentFeatureVars(1))/(parentFeatureVars(3)-parentFeatureVars(1)); elseif ~isinf(D2) %Where the parent feature is intersected by the child feature tParent = (parentFeatureVars(1)-childFeatureVars(1) + D2*(childFeatureVars(2)-parentFeatureVars(2)))/(D2*(parentFeatureVars(4)-parentFeatureVars(2))-parentFeatureVars(3)+parentFeatureVars(1)); %Where the child feature is intersected by the parent feature tChild = (parentFeatureVars(1)+(parentFeatureVars(3)-parentFeatureVars(1))*tParent-childFeatureVars(1))/(childFeatureVars(3)-childFeatureVars(1)); else tParent = Inf; tChild = Inf; end parentLength = ((parentFeatureVars(3)-parentFeatureVars(1))^2+(parentFeatureVars(4)-parentFeatureVars(2))^2)^0.5; childLength = ((childFeatureVars(3)-childFeatureVars(1))^2+(childFeatureVars(4)-childFeatureVars(2))^2)^0.5; tLenParent = tParent*parentLength; tLenChild = tChild*childLength; end