//////////////////////////////////////////////////////////////// // Copyright G. Allaire, A. Karrman, October 2009 // Last revision, January 2012, with G. Michailidis // // A Scilab toolbox for 2-d structural optimization // by the level set method. // // Based on the work of G. Allaire, F. Jouve, A.-M. Toader, // J. Comp. Phys. Vol 194/1, pp.363-393 (2004). //////////////////////////////////////////////////////////////// // // This file contains the subroutines: mesh0, mesh00, shift2n, // minmod, g, FEdensity, FE, FEtop, lk, c, compliance, volume, perimeter, // curv, solvelvlset, solvelvlset_top, regularize. // // The routines FE and lk (finite element solver) are courtesy of // O. Sigmund (A 99 line topology optimization code written in Matlab, // Struct. Multidisc. Optim. 21, pp.120-127, 2001) whose kind // permission is greatefully acknowledged. // //////////////////////////////////////////////////////////////// // // INITIALIZE THE STRUCTURE // This function will just distribute holes // uniformly throughout our mesh. hx is the number // of horizontal holes while hy is number of vertical // holes. r is a variable that lets you change the // size of the holes; a default size would be .5; .1 // would give you very small holes and .9 would give // you very large holes. (0 < r < 1) function [phi0] = mesh0(hx,hy,r) if isdef('ntop') then phi0 = -ones(nely+1,nelx+1) ; // Full-domain initialisation else phi0 = zeros(nely+1,nelx+1) ; // First create an empty matrix //In order to make our holes, we use a 3d function //that maps x and y coordinates to z = cos(x)*cos(y) //and then we flatten holes to .1 and the structure //part to -.1 phi0 = -cos((hy+1)*(nodey*%pi)/yheight)'*cos((hx+1)*(nodex*%pi)/xlength)+r-1 ; phi0 = .2*ceil(max(phi0,0))-.1 ; end endfunction // REINITIALIZE THE LEVEL SET FUNCTION // After intialization as well as periodically while // solving the transport equation, we reinitialize // the level set function, which means we solve // d(phi)/dt + sign(phi0)(|grad(phi)|-1) = 0 with // phi(t=0,x) = phi_0(x) // We do this to make sure the level set function // does not get too flat or steep at the zero // level set, which is where our structure boundary // is. phi is the level set funtion we are // reinitializing and RIiter is the number of times // we're solving the reinitialization equation. function phi00 = mesh00(phi,RIiter) cfl = 0.5 ; // CFL CONDITION // Using our CFL condition, we define the time step as dt0 = min(dx,dy)*cfl ; for n = 1 : RIiter // For our first derivatives: phin = shift2n('n',phi,'n') ; phis = shift2n('s',phi,'n') ; phie = shift2n('e',phi,'n') ; phiw = shift2n('w',phi,'n') ; // Our scheme is second order in space, so // we use these for our second derivatives: phinn = shift2n('n',phin,'n') ; phiss = shift2n('s',phis,'n') ; phiee = shift2n('e',phie,'n') ; phiww = shift2n('w',phiw,'n') ; // The first derivatives: dxm = (phi-phie)/dx ; dxp = (phiw-phi)/dx ; dym = (phi-phis)/dy ; dyp = (phin-phi)/dy ; // The second derivatives (our scheme is second // order so we use several different ones): dxmxm = (phi - 2*phie + phiee)/(dx^2) ; dxpxp = (phiww - 2*phiw + phi)/(dx^2) ; dxpxm = (phie - 2*phi + phiw)/(dx^2) ; dymym = (phi - 2*phis + phiss)/(dy^2) ; dypyp = (phinn - 2*phin + phi)/(dy^2) ; dypym = (phin - 2*phi + phis)/(dy^2) ; // From Sethian, our scheme requires four // main parts, as defined here: partA = dxm + .5*dx*minmod(dxmxm,dxpxm) ; partB = dxp - .5*dx*minmod(dxpxp,dxpxm) ; partC = dym + .5*dy*minmod(dymym,dypym) ; partD = dyp - .5*dy*minmod(dypyp,dypym) ; // Then we use these parts along with the // flux function to get: delp2 = g(partA,partB,partC,partD) ; delm2 = g(partB,partA,partD,partC) ; // sphi is the sign of phi and we added a small // value proportional to the mesh size in the // denominator to be sure we never divide by zero. nabla = 0.5*(dxm.^2 + dxp.^2 + dym.^2 + dyp.^2) ; sphi = phi./sqrt(phi.*phi+sqrt(dx^2+dy^2)*nabla/10) ; sphip = max(sphi,0) ; sphim = min(sphi,0) ; phi = phi - dt0*(sphip.*delp2 + sphim.*delm2 - sphi) ; end // Once our iterations are finished, we set // our output to the new level set function: phi00 = phi ; endfunction // SPACE SHIFT FUNCTION // Using Neumann, or Neumann for the directional gradient boundary conditions, we just shift // our matrix over one index in a certain direction // in order to take derivatives. function phishift = shift2n(direction,phi,conditions) // SHIFTS LEVEL SET FUNCTION WITH NEUMANN OR NEUMANN FOR THE GRADIENT CONDITIONS select direction case 'w' then // SHIFT WEST [m,n] = size(phi) ; phishift(1:m,1:n-1) = phi(1:m,2:n) ; select conditions case 'n' then phishift(1:m,n) = phi(1:m,n) ; // NEUMANN CONDITIONS case 'ng' then phishift(1:m,n) = 2*phi(1:m,n)-phi(1:m,n-1) ; // NEUMANN FOR THE DIRECTIONAL GRADIENT CONDITIONS end case 'e' then // SHIFT EAST [m,n] = size(phi) ; phishift(1:m,2:n) = phi(1:m,1:n-1) ; select conditions case 'n' then phishift(1:m,1) = phi(1:m,1) ; // NEUMANN CONDITIONS case 'ng' then phishift(1:m,1) = 2*phi(1:m,1)-phi(1:m,2) ; // NEUMANN FOR THE DIRECTIONAL GRADIENT CONDITIONS end case 'n' then // SHIFT NORTH [m,n] = size(phi) ; phishift(1:m-1,1:n) = phi(2:m,1:n) ; select conditions case 'n' then phishift(m,1:n) = phi(m,1:n) ; // NEUMANN CONDITIONS case 'ng' then phishift(m,1:n) = 2*phi(m,1:n)-phi(m-1,1:n) ; // NEUMANN FOR THE DIRECTIONAL GRADIENT CONDITIONS end case 's' then // SHIFT SOUTH [m,n] = size(phi) ; phishift(2:m,1:n) = phi(1:m-1,1:n) ; select conditions case 'n' then phishift(1,1:n) = phi(1,1:n) ; // NEUMANN CONDITIONS case 'ng' then phishift(1,1:n) = 2*phi(1,1:n)-phi(2,1:n) ; // NEUMANN FOR THE DIRECTIONAL GRADIENT CONDITIONS end else error('SHIFT N,S,E, OR W?') end endfunction // MINMOD FUNCTION // This function is essential for our second order // scheme used to solve the level set equation // (as well as for for reinitialization). function mout = minmod(phi1,phi2) sphi1=sign(phi1) ; sphi2=sign(phi2) ; mout = max(0,sphi1.*sphi2).*sphi1.*min(abs(phi1),abs(phi2)) ; endfunction // FLUX FUNCTION // This is the numerical flux of our scheme for Hamilton Jacobi equations function gout = g(u1,u2,v1,v2) gout = sqrt( max(u1,0.).^2 + min(u2,0.).^2... + max(v1,0.).^2 + min(v2,0.).^2 ) ; endfunction // FINITE ELEMENT MESH DENSITY FUNCTION // This function converts our level set function // into a density function defined for each rectangular // finite element. phi is the level set, eps is the // density of the very light ersatz material // that fills the holes. function FEtheta = FEdensity(phi,eps) // We originally begin with a matrix with empty elements // and then using phi we will compute material densities for each. FEtheta = zeros(nely,nelx) ; // In order to get rid of any nodes in phi with value 0 // (our algorithm does not deal well with this) we bump // the value above or below 0 (with 1/2 probability for each) e = .000001 ; indices = find(phi==0) ; for m = 1 : length(indices) phi(indices(m)) = (2*round(rand())-1)*e ; end for i = 1 : nely for j = 1 : nelx // We go through each element's phi values for the surrounding // nodes. If all are negative, the element has full (1) density // and if all are positive, then the element has 'eps' densitiy. phis = [phi(i,j) phi(i,j+1) phi(i+1,j+1) phi(i+1,j)] ; select sum(sign(phis)) case 4 then FEtheta(i,j) = eps ; case -4 then FEtheta(i,j) = 1 ; else // Otherwise, we split the element into four // triangles by making cuts across the elements // two diagonals. We then sum the surrounding nodes // to get a 5th value for the center value of the element. // We don't want it to be 0, so we bump it positive // or negative according to equal probabilities. phis($+1) = sum(phis(1:4)) ; if phis(5) == 0 then phis(5) = (2*round(rand())-1)*e ; end // The triangle matrix has 4 rows corresponding // to each triangle (the first is the top triangle, // followed by the other three in the clockwise direction. triMat = [phis(5) phis(1) phis(2) phis(5) phis(2) phis(3) phis(5) phis(3) phis(4) phis(5) phis(4) phis(1)] ; // We go through each triangle and add a certain amount of // density to the element's density value. for k = 1 : 4 select sum(sign(triMat(k,:))) // If the triangle is completely part of the structure: case -3 then FEtheta(i,j) = FEtheta(i,j)+.25 ; // If it is not at all part of the structure: case 3 then FEtheta(i,j) = FEtheta(i,j)+.25*eps ; else // Otherwise, we use linear interpolation to get an // accurate estimate of how much density the triangle // should contribute to this element. n = find(sign(triMat(k,:))~=sign(sum(triMat(k,:)))) ; phicc = triMat(k,modulo(n+1,3)+1) ; phimid = triMat(k,n) ; phic = triMat(k,modulo(n,3)+1) ; f1 = phimid/(phimid-phicc) ; f2 = phimid/(phimid-phic) ; if sum(sign(triMat(k,:))) == 1 then FEtheta(i,j) = FEtheta(i,j)+.25*(1-f1*f2)*eps+.25*f1*f2 ; else FEtheta(i,j) = FEtheta(i,j)+.25*(1-f1*f2)+.25*f1*f2*eps ; end end end end end end endfunction // FE ANALYSIS AND CALCULATION OF THE SHAPE GRADIENT // From Ole Sigmund's Matlab code, we use the // density function along with the key finite element // vectors and run the analysis in order to get // the elastic displacement which is used to // compute the compliance and the velocity field // for the transport level set equation. function [lvlAeueu, FEAeueu]=FE(FEtheta,KE,K,F,U) // ANALYSIS for elx = 1:nelx for ely = 1:nely n1 = (nely+1)*(elx-1)+ely; n2 = (nely+1)* elx +ely; edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2]; K(edof,edof) = K(edof,edof) + FEtheta(ely,elx)*KE; end end // SOLVING // Using Cholevsky factorization, we solve for U. If // we have multiple forces, we must solve for each // situation separately. sK = sparse(K(freedofs,freedofs)) ; sKcho = chfact(sK) ; for i = 1 : forceCount U(freedofs,i) = chsolve(sKcho,full(F(freedofs,i))) ; end U(fixeddofs,:)= 0; // COMPUTE THE VELOCITY FIELD (THE ELASTIC ENERGY DENSITY) FEAeueu = zeros(nely,nelx) ; for ely = 1:nely for elx = 1:nelx n1 = (nely+1)*(elx-1)+ely; n2 = (nely+1)* elx +ely; for i = 1 : forceCount Ue = U([2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2],i); FEAeueu(ely,elx) = FEAeueu(ely,elx)+FEtheta(ely,elx)*Ue'*KE*Ue; end end end // The elastic energy density has just been computed // in each finite element but we need it at each node // for further use in the transport level set equation. // So we just average the elastic energy density // for each finite element surrounding each node. lvlAeueu = zeros(nely+1,nelx+1) ; for i = 1:(nely+1) for j = 1:(nelx+1) bl = FEAeueu(max(1,i-1),max(1,j-1)) ; br = FEAeueu(min(nely,i),max(1,j-1)) ; tr = FEAeueu(min(nely,i),min(nelx,j)) ; tl = FEAeueu(max(1,i-1),min(nelx,j)) ; lvlAeueu(i,j) = (bl+br+tr+tl)/4 ; end end endfunction // FE ANALYSIS AND CALCULATION OF THE TOPOLOGICAL GRADIENT // From Ole Sigmund's Matlab code, we use the // density function along with the key finite element // vectors and run the analysis in order to get // the elastic displacement which is used to // compute the compliance and the topological gradient // for the change of the level set function. function [FEAeueu,FEAeueucomp]=FEtop(FEtheta,KE,K,F,U) E = 1.; // Young modulus nu = 0.3; // Poisson's coefficient lambda=E*nu/(1+nu)/(1-2*nu); // Lame coefficients mi=E/(1+nu)/2.; D=E/(1-nu^2)*[1 nu 0;nu 1 0;0 0 1-nu]; // Hooke's tensor // ANALYSIS for elx = 1:nelx for ely = 1:nely n1 = (nely+1)*(elx-1)+ely; n2 = (nely+1)* elx +ely; edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2]; K(edof,edof) = K(edof,edof) + FEtheta(ely,elx)*KE; end end // SOLVING // Using Cholevsky factorization, we solve for U. If // we have multiple forces, we must solve for each // situation separately. sK = sparse(K(freedofs,freedofs)) ; sKcho = chfact(sK) ; for i = 1 : forceCount U(freedofs,i) = chsolve(sKcho,full(F(freedofs,i))) ; end U(fixeddofs,:)= 0; // COMPUTE THE VELOCITY FIELD (TOPOLOGICAL GRADIENT) FEAeueu = zeros(nely+1,nelx+1) ; FEAeueucomp = zeros(nely,nelx); FEAeueu=FEAeueu(:); for ely = 1:nely for elx = 1:nelx n1 = (nely+1)*(elx-1)+ely; n2 = (nely+1)* elx +ely; for i = 1 : forceCount Ue = U([2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2],i); e(1,1)=(Ue(3)-Ue(1))/dx; e(1,2)=(Ue(8)-Ue(2))/dy; e(1,3)=0.5*((Ue(7)-Ue(1))/dy+(Ue(4)-Ue(2))/dx); e(2,1)=(Ue(5)-Ue(7))/dx; e(2,2)=(Ue(8)-Ue(2))/dy; e(2,3)=0.5*((Ue(7)-Ue(1))/dy+(Ue(6)-Ue(8))/dx); e(3,1)=(Ue(5)-Ue(7))/dx; e(3,2)=(Ue(6)-Ue(4))/dy; e(3,3)=0.5*((Ue(5)-Ue(3))/dy+(Ue(6)-Ue(8))/dx); e(4,1)=(Ue(3)-Ue(1))/dx; e(4,2)=(Ue(6)-Ue(4))/dy; e(4,3)=0.5*((Ue(5)-Ue(3))/dy+(Ue(4)-Ue(2))/dx); s=zeros(3,4); for j=1:4 s(:,j)=FEtheta(ely,elx)*D*(e(j,:)'); end s=s'; FEAeueu(n1)=FEAeueu(n1)+(%pi)*(lambda+2*mi)/(2*mi)/(lambda+mi)*(4*mi*(e(1,1)*s(1,1)+e(1,2)*s(1,2)+2*e(1,3)*s(1,3))+(lambda-mi)*(e(1,1)+e(1,2))*(s(1,1)+s(1,2))); FEAeueu(n1+1)=FEAeueu(n1+1)+(%pi)*(lambda+2*mi)/(2*mi)/(lambda+mi)*(4*mi*(e(2,1)*s(2,1)+e(2,2)*s(2,2)+2*e(2,3)*s(2,3))+(lambda-mi)*(e(2,1)+e(2,2))*(s(2,1)+s(2,2))); FEAeueu(n2)=FEAeueu(n2)+(%pi)*(lambda+2*mi)/(2*mi)/(lambda+mi)*(4*mi*(e(4,1)*s(4,1)+e(4,2)*s(4,2)+2*e(4,3)*s(4,3))+(lambda-mi)*(e(4,1)+e(4,2))*(s(4,1)+s(4,2))); FEAeueu(n2+1)=FEAeueu(n2+1)+(%pi)*(lambda+2*mi)/(2*mi)/(lambda+mi)*(4*mi*(e(3,1)*s(3,1)+e(3,2)*s(3,2)+2*e(3,3)*s(3,3))+(lambda-mi)*(e(3,1)+e(3,2))*(s(3,1)+s(3,2))); FEAeueucomp(ely,elx) = FEAeueucomp(ely,elx)+FEtheta(ely,elx)*Ue'*KE*Ue; end end end FEAeueu=matrix(FEAeueu,nely+1,nelx+1)/4; for i=1:nely+1 for j=1:nelx+1 if i==1 | i==nely+1 then FEAeueu(i,j)=FEAeueu(i,j)*2; end if j==1 | j==nelx+1 then FEAeueu(i,j)=FEAeueu(i,j)*2; end end end endfunction // STIFFNESS MATRIX // This matrix, found using symbolic manipulations, // is key for analyzing our elastic structure using // the FE routine. (Again, thanks to Ole Sigmund) function [KE]=lk() E = 1.; nu = 0.3; k=[ 1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ... -1/4+nu/12 -1/8-nu/8 nu/6 1/8-3*nu/8]; KE = E/(1-nu^2)*[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8) k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3) k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2) k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5) k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4) k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7) k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6) k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)]; endfunction // MAP A CONDITION LOCATION TO THE FE MATRIX // This function was meant to make placing boundary // conditions (thus the 'c') as well as forces // on the structure much easier and intuitive. // fV is a vector of 3 components. The first // two components indicate the fraction along the // x and y axis we are (with the top left corner // as the origin). The direction variable is binary // and means in the y direction if it is 1 and in the // x direction if it is 0. For forces, this direction // variable literally means what direction the force // is in whereas for the conditions we set, it indicates // in which direction a constraint is fixed. function elnum = c(fV) xfrac = fV(1) ; yfrac = fV(2); direction = fV(3) ; X = zeros(nely+1,nelx+1) ; xvector = 1:2:(2*(nely+1)-1) ; // From the FE vectors, we grab the x components of // the vector fields and then lay the index in the vector // on a matrix that represents our discretized working domain: for i = 1 : (nelx+1) X(:,i) = 2*(nely+1)*(i-1)+xvector' ; end // The y values are just one above the x values // because in the FE vectors, they're arranged in // pairs. So: Y = X+1 ; i = floor(yfrac*nely)+1 ; j = floor(xfrac*nelx)+1 ; // Here we just extract whichever index of the // vector that we need using the matrices we created: if direction == 0 then elnum = X(i,j) ; else elnum = Y(i,j) ; end endfunction // COMPUTE THE COMPLIANCE // The compliance is the main part of our objective function. // It is the integral of the elastic energy density, equal to // the total work done by the applied forces. function totcomp = compliance(FEAeueu) totcomp = sum(FEAeueu) ; endfunction // THE VOLUME OF THE STRUCTURE // To find the total volume of our structure, we just // integrate the density*dV: function totvol = volume(FEtheta) dV = dx*dy ; totvol = sum(FEtheta)*dV ; endfunction // THE PERIMETER // In order to roughly calculate the perimeter // we find the norm of the gradient of the sign // of phi and integrate it, then divide by 2: function totperim = perimeter(phi) // To smooth sign(phi): epsperim = min(dx,dy)/20 ; sx = phi./sqrt(phi.^2+epsperim^2) ; sxn = shift2n('n',sx,'ng') ; sxs = shift2n('s',sx,'ng') ; sxe = shift2n('e',sx,'ng') ; sxw = shift2n('w',sx,'ng') ; // We now calculate d(phi)/dx and d(phi)/dy: dsxx = (sxw-sxe)/(2*dx) ; dsxy = (sxn-sxs)/(2*dy) ; dV = dx*dy ; // And then integrate: totperim = .5*sum(sqrt(dsxx.^2+dsxy.^2))*dV ; endfunction // THE CURVATURE // The shape gradient of the perimeter is the // mean curvature which we compute by // div(grad(phi)/|grad(phi)|) // where phi is the level set function. function H = curv(phi) // When finding the normal vector, we need // to divide by the norm of the gradient of phi; // we use this small value to make sure that // the gradient of phi never goes to 0: epscurv = min(dx,dy)/20 ; // Here are the first derivatives for finding // the gradient of phi: phin = shift2n('n',phi,'ng') ; phis = shift2n('s',phi,'ng') ; phie = shift2n('e',phi,'ng') ; phiw = shift2n('w',phi,'ng') ; dphix = (phiw-phie)/(2*dx) ; dphiy = (phin-phis)/(2*dy) ; // Here is |grad(phi)| and then the x and y // components of the normal vector field: mag = sqrt(dphix.^2+dphiy.^2+epscurv^2) ; nx = dphix./mag ; ny = dphiy./mag ; // Now to find the divergence, just take the // partials with repsect to x and y of the // x and y components of our normal vector field // (respectively) and then add these together // to get our end mean curvature, which is a // function across our working domain. nxe = shift2n('e',nx,'ng') ; nxw = shift2n('w',nx,'ng') ; nyn = shift2n('n',ny,'ng') ; nys = shift2n('s',ny,'ng') ; divnx = (nxw-nxe)/(2*dx) ; divny = (nyn-nys)/(2*dy) ; H = divnx+divny ; endfunction // SOLVE THE TRANSPORT LEVEL SET EQUATION // This function solves an Hamilton Jacobi equation // to transport the level set function: // d(phi)/dt + V |grad(phi)| = 0 with // where V is the normal velocity (the shape gradient // of the objective function). // The scheme is second order in space and explicit // first order in time. The time step is dt and the // number of time steps (or iterations) is HJiter. function solvedphi = solvelvlset(phi,V,dt,HJiter,lagP) for i = 1 : HJiter // We reinitialize the level set function for every RIfreq // time steps while solving the transport level set equation: if modulo(i,RIfreq) == 0 then phi = mesh00(phi,RIiter) ; end // Our scheme takes into account whether our // front is moving forward or backward: Vp = max(V,0) ; Vm = min(V,0) ; // We use these for our first derivatives: phin = shift2n('n',phi,'ng') ; phis = shift2n('s',phi,'ng') ; phie = shift2n('e',phi,'ng') ; phiw = shift2n('w',phi,'ng') ; // We use these for our second derivatives: phinn = shift2n('n',phin,'ng') ; phiss = shift2n('s',phis,'ng') ; phiee = shift2n('e',phie,'ng') ; phiww = shift2n('w',phiw,'ng') ; // The first derivatives: dxm = (phi-phie)/dx ; dxp = (phiw-phi)/dx ; dym = (phi-phis)/dy ; dyp = (phin-phi)/dy ; // Because our scheme is second order in space, // we use the second derivatives: dxmxm = (phi - 2*phie + phiee)/(dx^2) ; dxpxp = (phiww - 2*phiw + phi)/(dx^2) ; dxpxm = (phie - 2*phi + phiw)/(dx^2) ; dymym = (phi - 2*phis + phiss)/(dy^2) ; dypyp = (phinn - 2*phin + phi)/(dy^2) ; dypym = (phin - 2*phi + phis)/(dy^2) ; // Our scheme uses certain parts to which // we apply the flux function partA = dxm + .5*dx*minmod(dxmxm,dxpxm) ; partB = dxp - .5*dx*minmod(dxpxp,dxpxm) ; partC = dym + .5*dy*minmod(dymym,dypym) ; partD = dyp - .5*dy*minmod(dypyp,dypym) ; delp2 = g(partA,partB,partC,partD) ; delm2 = g(partB,partA,partD,partC) ; // Here we compute the |grad(phi)|. epscurv = min(dx,dy)/20 ; dphix = (phiw-phie)/(2*dx) ; dphiy = (phin-phis)/(2*dy) ; mag = sqrt(dphix.^2+dphiy.^2+epscurv^2) ; // We update our old level set: phi = phi - dt*(delp2.*Vp + delm2.*Vm)+dt*lagP*curv(phi).*mag ; end // We let our solved level set be the output: solvedphi = phi ; endfunction // CHANGE OF THE LEVEL-SET FUNCTION IN THE TOPOLOGICAL GRADIENT METHOD function phiTest = solvelvlset_top(phi,V,percentage,FEtheta) // To update the shape in a topological gradient step, we remove a certain percentage of volume // at the areas where the topological gradient takes the most negative values. totvolin = volume(FEtheta); minV=10^30.; phiTest=phi; // We find the most negative value of the topological gradient. for i=1:nely+1 for j=1:nelx+1 if phi(i,j)<0. then if V(i,j)<0. then minV=min(minV,V(i,j)); end end end end // Then, we change the sign of the level-set function at the points // where the value of the velocity is lower than perv*minV. Perv is // chosen arbitrarily in the beginning and will be adjjusted in the // sequel in order to remove the desired percentage of volume. if abs(minV) ~= 0. then perv=0.95; totvoltar=(1-percentage)*totvolin; totvolTest=0.; rit=0.; for i=1:nely+1 for j=1:nelx+1 if phi(i,j)<0. then if V(i,j)<0. then if abs(V(i,j))>perv*abs(minV) then phiTest(i,j) = -phi(i,j); end end end end end // The new volume is calculated. phiTest = mesh00(phiTest,100); FEthetaTest = FEdensity(phiTest,eps) ; FEthetaTest = passive(FEthetaTest) ; totvolTest=volume(FEthetaTest); // We try to determine the appropriate value of perv // using a dichotomy algorithm. if (abs(totvolTest-totvoltar)/totvolin)>0.005 then if totvolTest>totvoltar then pervmax=perv; pervmin=perv; while (totvolTest>totvoltar) & (pervmin>0.) pervmax=pervmin; pervmin=max(pervmin-0.01,0.); phiTest=phi; for i=1:nely+1 for j=1:nelx+1 if phi(i,j)<0. then if V(i,j)<0. then if abs(V(i,j))>pervmin*abs(minV) then phiTest(i,j) = -phi(i,j); end end end end end phiTest=mesh00(phiTest,100); FEthetaTest = FEdensity(phiTest,eps) ; FEthetaTest = passive(FEthetaTest) ; totvolTest=volume(FEthetaTest); end else pervmax=perv; pervmin=perv; while (totvolTest<=totvoltar) & (pervmax<1.) pervmin=pervmax; pervmax=min(pervmax+0.01,1.); phiTest=phi; for i=1:nely+1 for j=1:nelx+1 if phi(i,j)<0. then if V(i,j)<0. then if abs(V(i,j))>pervmax*abs(minV) then phiTest(i,j) = -phi(i,j); end end end end end phiTest=mesh00(phiTest,100); FEthetaTest = FEdensity(phiTest,eps) ; FEthetaTest = passive(FEthetaTest) ; totvolTest=volume(FEthetaTest); end end end while ((abs(totvolTest-totvoltar)/totvolin)>0.005) & (rit<100) & (pervmax<1.) & (pervmin>0.) rit=rit+1.; perv=(pervmin+pervmax)/2.; phiTest=phi; for i=1:nely+1 for j=1:nelx+1 if phi(i,j)<0. then if V(i,j)<0. then if abs(V(i,j))>perv*abs(minV) then phiTest(i,j) = -phi(i,j); end end end end end phiTest=mesh00(phiTest,100); FEthetaTest = FEdensity(phiTest,eps) ; FEthetaTest = passive(FEthetaTest) ; totvolTest=volume(FEthetaTest); if totvolTest>totvoltar then pervmax=perv; else pervmin=perv; end end end phiTest=mesh00(phiTest,100); endfunction // REGULARIZATION OF THE ADVECTION VELOCITY // Here we regularize the velocity field by considering the H1 product instead of the L2 one. function v = regularize(phi,V,e2,M,lagP) k1 = e2/(dx^2) ; // Coefficients of the regularization matrix. k2 = e2/(dy^2) ; // Now we calculate the surface Dirac function. epsperim = min(dx,dy)/20 ; sx = phi./sqrt(phi.^2+epsperim^2) ; sxn = shift2n('n',sx,'ng') ; sxs = shift2n('s',sx,'ng') ; sxe = shift2n('e',sx,'ng') ; sxw = shift2n('w',sx,'ng') ; //We now calculate d(phi)/dx and d(phi)/dy: dsxx = (sxw-sxe)/(2*dx) ; dsxy = (sxn-sxs)/(2*dy) ; delta = 0.5*sqrt(dsxx.^2+dsxy.^2); // Surface Dirac function b = -V.*delta; // The right hand side is defined on the boundary so we use the delta function. b=b(:); // We make the right hand side a vector, by taking the matrix column by column. // We form the regularization matrix with Neumann conditions. for i = 1:nely+1 for j = 1:nelx+1 k = (j-1)*(nely+1)+i; M(k,k) = 1.+2*k1+2*k2; if (i>1) then M(k,k-1) = -k2; else M(k,k)=M(k,k)-k2; end if (i1) then M(k,k-(nely+1)) = -k1; else M(k,k)=M(k,k)-k1; end if (j