tic %Start timer % Inputs % Cantilever beam initialisation inputVariable = [-108.797979211583,65.4435901909376,488.286882509882,235.094790586772,97.1191736430205,353.358259669744,183.229340471321,1040.64752625714,483.499281766377,105.482993125990,1245.17000000000,143.880000000000,1914.36000000000,519.210000000000,85.8400000000000,2266.72458336935,103.599913896207,3075.49878607669,541.543891303967,98.6392171290896,-165.385770418896,34.6702323605109,1294.71690116466,82.4941426407934,225.427701893167,1013.96554469604,481.869437007546,1276.10572684622,104.548183536541,53.0824222272433,1164.72286094920,20.6376767011557,2306.01016629544,80.7801099250218,124.690171263730,1936.53000000000,456.560000000000,2205.78000000000,90.7800000000000,55.3800000000000,84.8554968244647,935.110381354288,1337.04740554466,898.808202407775,189.928655698751,1018.69363673509,475.766727602642,1321.87899384919,888.498219393463,56.3593086547422,1012.23199075083,980.208552595332,2260.61617008326,933.296479022833,116.346615518285,1816.20000000000,417.230000000000,2195.25000000000,948.590000000000,54.5600000000000,-53.8673576299499,959.525237989423,390.022754316555,826.811407004563,173.207446220761,280.210000000000,837.990000000000,1064.90000000000,455.680000000000,105.010000000000,1295.66000000000,849.570000000000,1948.78000000000,457.130000000000,91.7100000000000,2217.66548916056,913.984739997730,3077.34357000791,464.633438460650,98.2481879010957]; N=16; %Numer of features % % MBB Beam % inputVariable = [-4265.034285 42.82528376 1155.260374 144.9266098 173.4104155 -1243674.234 -12545.20966 23339899.54 235650.6125 108.1001602 -2071.071538 2728.667554 44518.39879 -36553.42301 148.4214507 -282.4563604 985.2175981 2221.396183 879.5440855 216.8443741 -3616.24827 -3314.856766 2102.056271 925.4278817 133.0192464 2196.50903 834.0281359 3055.965485 -18.87213008 125.5687262]; % N=6; %Numer of features 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'; calculationSelection = 'NodeBased'; % Either 'NodeBased' or 'ElementBased' % Design domain properties x0 = 0; y0 = 0; xDim = 3000; %Domain x dimension yDim = 1000; %Domain y dimension elSize = 5; %Mesh element size % Threshold Parameters angleLimitThreshold = 0; %for identifying and enlarging collinear components at the end of the program (e.g. 0, 5) voidPenalty = 10; solidPenalty = 10; % Visualisation showFeatureSpines = false; showFeatureIds = false; % Variable initialisation compIdList = 1:N; xy00=inputVariable(:); enlargeInfo = cell(N,6); % records % Create mesh [xElCount,yElCount,nodeCoords,elementNodes,x,y] = CreateMesh(x0,y0,xDim,yDim,elSize); % Create feature representation [Phi,~,H_master,H_features]=FeatureConstruction(xy00,featureVars,N,nodeCoords,p,alpha,epsilon,aggregationSelection,xElCount,yElCount); [featureSize,~,intersectingFeatures,~] = FeatureComparison(H_features,alpha,N,elementNodes,elSize); [~,I] = sort((featureSize/max(featureSize)),'descend'); for k = 1:N selectedFeature = compIdList(I(k)); enlargeInfo{k,1} = selectedFeature; intersectingFeatureList = intersectingFeatures{selectedFeature}; if isempty(intersectingFeatureList) enlargeInfo{k,2} = "No intersecting features"; else enlargeInfo{k,2} = intersectingFeatureList; for i = 1:length(intersectingFeatureList) x_init = xy00(featureVars*(selectedFeature)-featureVars+1: featureVars*(selectedFeature)); phiSubSelection=cell(2,1); phiSubSelection{1,1} = Phi{selectedFeature}; phiSubSelection{2,1} = Phi{intersectingFeatureList(i)}; aggregatedPhi = FeatureAggregation(phiSubSelection,size(phiSubSelection,1),aggregationSelection); H_featuresN_M = Heaviside(aggregatedPhi,alpha,xElCount,yElCount,epsilon); H_enlargingFeature = Heaviside(Phi{selectedFeature},alpha,xElCount,yElCount,epsilon); switch calculationSelection case 'NodeBased' f = @(result)EnlargeFeature(result,alpha,xElCount,yElCount,epsilon,nodeCoords(:,1),nodeCoords(:,2),p,H_featuresN_M,H_enlargingFeature,voidPenalty,solidPenalty); case 'ElementBased' featuresN_MArea = AreaCalc(H_featuresN_M,elementNodes,elSize); enlargingFeatureArea = AreaCalc(H_enlargingFeature,elementNodes,elSize); f = @(result)EnlargeFeature2(result,alpha,xElCount,yElCount,epsilon,nodeCoords(:,1),nodeCoords(:,2),p,featuresN_MArea,enlargingFeatureArea,voidPenalty,solidPenalty,elementNodes,elSize); otherwise warning('Unexpected calculation selection, reverted to NodeBased') f = @(result)EnlargeFeature(result,alpha,xElCount,yElCount,epsilon,nodeCoords(:,1),nodeCoords(:,2),p,H_featuresN_M,H_enlargingFeature,voidPenalty,solidPenalty); end % Optimisation options = optimoptions(@lsqnonlin,'MaxFunEvals',5000); result = lsqnonlin(f,x_init,[],[],options); xy00(featureVars*(selectedFeature)-featureVars+1: featureVars*(selectedFeature)) = result(1: featureVars); [Phi,~,~,~]=FeatureConstruction(xy00,featureVars,N,nodeCoords,p,alpha,epsilon,aggregationSelection,xElCount,yElCount); end [originalAreaSum,deltaAreaSum] = AreaDeviationFromInput(H_master,xy00,featureVars,N,nodeCoords,p,xElCount,yElCount,alpha,epsilon,aggregationSelection,elementNodes,elSize); enlargeInfo{k,3} = deltaAreaSum/originalAreaSum*100; for i = 1:N index = find(intersectingFeatures{i} == selectedFeature); %Remove the enlarged feature from the lists of intersecting features to ensure other features don't attempt to overlap it if ~isempty(index) intersectingFeatures{i}(index) = []; end end end end [Phi,~,~,H_features]=FeatureConstruction(xy00,featureVars,N,nodeCoords,p,alpha,epsilon,aggregationSelection,xElCount,yElCount); [featureSize,~,intersectingFeatures,~] = FeatureComparison(H_features,alpha,N,elementNodes,elSize); % Calculating feature angles featureCoords = reshape(xy00,featureVars,N)'; featureAngles = atan2(featureCoords(:,4)-featureCoords(:,2),featureCoords(:,3)-featureCoords(:,1)); featureAngles(logical(featureAngles(:,1)<0),1) = featureAngles(logical(featureAngles(:,1)<0),1) + pi(); %Increase the overlap of collinear components. for i = 1:N enlargeInfo{i,4} = i; if isempty(intersectingFeatures{i}) enlargeInfo{i,5} = 'No collinear comps to test'; else testedFeatures = []; interlist = intersectingFeatures{i}; for j = 1:length(interlist) if abs(featureAngles(interlist(j))-featureAngles(i)) < pi()/180*angleLimitThreshold || abs(pi()-featureAngles(interlist(j))-featureAngles(i)) < pi()/180*angleLimitThreshold testedFeatures(end+1) = interlist(j); x_init = xy00(featureVars*(i)-featureVars+1: featureVars*(i)); aggregatedPhi=Phi{i}; aggregatedPhi=max(aggregatedPhi,Phi{interlist(j)}); H_featuresN_M = Heaviside(aggregatedPhi,alpha,xElCount,yElCount,epsilon); H_enlargingFeature = Heaviside(Phi{i},alpha,xElCount,yElCount,epsilon); switch calculationSelection case 'NodeBased' f = @(result)EnlargeFeature(result,alpha,xElCount,yElCount,epsilon,nodeCoords(:,1),nodeCoords(:,2),p,H_featuresN_M,H_enlargingFeature,voidPenalty,solidPenalty); case 'ElementBased' featuresN_MArea = AreaCalc(H_featuresN_M,elementNodes,elSize); enlargingFeatureArea = AreaCalc(H_enlargingFeature,elementNodes,elSize); f = @(result)EnlargeFeature2(result,alpha,xElCount,yElCount,epsilon,nodeCoords(:,1),nodeCoords(:,2),p,featuresN_MArea,enlargingFeatureArea,voidPenalty,solidPenalty,elementNodes,elSize); otherwise warning('Unexpected calculation selection, reverted to NodeBased') f = @(result)EnlargeFeature(result,alpha,xElCount,yElCount,epsilon,nodeCoords(:,1),nodeCoords(:,2),p,H_featuresN_M,H_enlargingFeature,voidPenalty,solidPenalty); end % Optimisation options = optimoptions(@lsqnonlin,'MaxFunEvals',5000); result = lsqnonlin(f,x_init,[],[],options); xy00(featureVars*(i)-featureVars+1: featureVars*(i)) = result(1: featureVars); [Phi,~,~,~]=FeatureConstruction(xy00,featureVars,N,nodeCoords,p,alpha,epsilon,aggregationSelection,xElCount,yElCount); end end enlargeInfo{i,5} = testedFeatures; [originalAreaSum,deltaAreaSum] = AreaDeviationFromInput(H_master,xy00,featureVars,N,nodeCoords,p,xElCount,yElCount,alpha,epsilon,aggregationSelection,elementNodes,elSize); enlargeInfo{i,6} = deltaAreaSum/originalAreaSum*100; end end %Row vector output of updated feature design variables updatedDesignVariables = xy00'; %Generate updated feature representation [Phi,~,~,~]=FeatureConstruction(xy00,featureVars,N,nodeCoords,p,alpha,epsilon,aggregationSelection,xElCount,yElCount); FeaturePlot(Phi,N,alpha,epsilon,xDim,yDim,x,y,xElCount,yElCount,showFeatureSpines,xy00,featureVars,showFeatureIds) toc %End timer function F=EnlargeFeature(result,alpha,xElCount,yElCount,epsilon,nodeCoordsX,nodeCoordsY,p,H_featuresN_M,H_enlargingFeature,voidPenalty,solidPenalty) phi = tPhi(result,nodeCoordsX,nodeCoordsY,p); H_updatedFeatureN = Heaviside(phi,alpha,xElCount,yElCount,epsilon); voidElementsInverse = logical(H_featuresN_M); %Locations where the union of features n and m is non-void solidElements = logical(H_enlargingFeature==1); %Locations where the feature being enlarged (n) is solid (doesn't account for partially filled elements) F = (H_featuresN_M(:) - H_updatedFeatureN(:)); F(~voidElementsInverse) = F(~voidElementsInverse).*voidPenalty; %Locations where the union of features n and m is void F(solidElements) = F(solidElements).*solidPenalty; %Locations where feature n is solid end function F=EnlargeFeature2(result,alpha,xElCount,yElCount,epsilon,nodeCoordsX,nodeCoordsY,p,featuresN_MArea,enlargingFeatureArea,voidPenalty,solidPenalty,elementNodes,elSize) phi = tPhi(result,nodeCoordsX,nodeCoordsY,p); H_updatedFeatureN = Heaviside(phi,alpha,xElCount,yElCount,epsilon); featureNArea = AreaCalc(H_updatedFeatureN,elementNodes,elSize); voidElementsInverse = logical(featuresN_MArea); %Locations where the union of features n and m is non-void solidElements = logical(enlargingFeatureArea==1); %Locations where the feature being enlarged (n) is solid (doesn't account for partially filled elements) F = (featuresN_MArea(:) - featureNArea(:)); F(~voidElementsInverse) = F(~voidElementsInverse).*voidPenalty; %Locations where the union of features n and m is void F(solidElements) = F(solidElements).*solidPenalty; %Locations where feature n is solid end