FE Verification using photoelasticity and constrained shape optimization

function fiberangles = surfacegen(perturbation,intpts) %AUTHOR: Veera Sundararaghavan %example: fiberangles = surfacegen(5*rand(5),intpts) %input perturbations to nurbs gcontrol points (5x5 matrix for z coordinates only) %input integration points of elements nintptsx2 matrix %output fiberangles at each integration point in degrees nintptsx1 vector %Code also displays the composite fiber distribution close all nrb = nrbtestsrf; % Default Surface res = [20,20]; [x,y,z] = nrb2xyz(nrb,res); [xn,yn,zn] = nrb2net(nrb);%5x5 control point grids zn = yn';%initialize with 90 degree fiber orientation %TODO can also perturb xn and yn for inner control points and xn or yn (not %both) for boundary control points. But also need to identify bounds on the %perturbation to avoid interpenetration. Use gnurbs to graphically play %with the control points and surfaces generated zn = zn + perturbation; nrb.coefs = net2nrb(xn,yn,zn); [x,y,z] = nrb2xyz(nrb,res); x = (x-5)*20.85./10; y = (y-5)*64.90./10; cc = contour3(x,y,z,30);axis equal;hold on; [xc,yc,zc] = sphere(10);f = 12.98/2;surf(xc.*f,yc.*f,(zc.*f)+100,'FaceColor','k'); view(2); fiberangles = contourmatch(intpts,cc'); function [x,y,z] = nrb2net(srf) % NRB2NET Get control net from NURBS structure x=squeeze(srf.coefs(1,:,:)); y=squeeze(srf.coefs(2,:,:)); z=squeeze(srf.coefs(3,:,:)); function nrb = net2nrb(x,y,z) % NET2NRB Get NURBS coefficients from control net [m,n]=size(x); x = reshape(x,1,m,n); y = reshape(y,1,m,n); z = reshape(z,1,m,n); one = ones(1,m,n); nrb = cat(1,x,y); nrb = cat(1,nrb,z); nrb = cat(1,nrb,one); function [x,y,z] = nrb2xyz(nrb,res) % NRB2XYZ Evaluate NURBS with a given resolution n = length(nrb.order); if n > 1 nu = res(1); nv = res(2); p = nrbeval(nrb,{linspace(0,1,nu),linspace(0,1,nv)}); else p = nrbeval(nrb,linspace(0,1,res(1))); end x = squeeze(p(1,:,:)); y = squeeze(p(2,:,:)); z = squeeze(p(3,:,:)); function fiberangles = contourmatch(points,cc) linenum = 1; index = 1; A = []; while linenum < length(cc) npts = cc(linenum,2); A = [A; [cc(linenum+1:linenum+npts,:) repmat(index,npts,1)]]; index = index+1; linenum = linenum+npts+1; end for i = 1:size(points,1) dist = A(:,1:2) - repmat(points(i,:),length(A),1); norms = sqrt(sum(dist.^2,2)); [~,I] = min(norms); curpt = A(I,1:2); index = A(I,3); I2 = find(A(:,3)==index); B = A(I2,1:2); dist2 = B - repmat(points(i,:),length(B),1); norms2 = sqrt(sum(dist2.^2,2)); [~,I] = sort(norms2); nextpt = B(I(2),1:2); if (nextpt(1)-curpt(1)) ~= 0 fiberangles(i) = atand((nextpt(2)-curpt(2))/(nextpt(1)-curpt(1))); else fiberangles(i) = 90; end %%uncomment for graphic display for point location(with angle) and contour selected % text(curpt(1),curpt(2),'A'); % text(nextpt(1),nextpt(2),'B') % text(points(i,1),points(i,2),num2str(fiberangles(i))); end