Level set constraints

[latexpage]

Constrained minimization is a common Morpho task. Two broad types of constraint can be imposed: Global constraints affect the entire body such as fixing the area or volume enclosed (we use this in the loop minimization example for instance). Local constraints on the other hand are enforced at every individual mesh point or element. One of the key ideas of Morpho is that any energy functional can equally well be used as either a global or local constraint, which provides a lot of flexibility.

A level set of a function is the locus of points within the domain for which the function returns the same value. For example, the contours on a topographic map form level sets of the height function:

Similarly, level sets of the function $f(\mathbf{x})=\mathbf{x} \cdot \mathbf{x}$ in 3D space are concentric spheres.

We can use level sets to prevent bodies in Morpho from entering particular regions of space. In the rest of this tutorial, we’ll minimize the shape of a loop at constant area subject to two solid boundaries above and below the loop.

Let’s begin by setting up the loop and defining the energy functional. This is very similar to loop minimization example so feel free to refer back to that sequence of tutorials. In this first section, we define some constants for later use and create an initially circular loop with 50 points:

/* Level set constraints */

np=50 // Number of points
A=1.0 // Initial radius of loop 
Nsteps=100 // Number of steps to take
y0 = 0.9 // Initial position of level set constraints
sl = 0.25 // Scale limit

// Create the manifold
m=manifoldline({A*cos(t), A*sin(t),0}, {t,0,2*pi-pi/np/2, 2*pi/np}, closed=true)

show(m)

If you run this, the following initial loop will be displayed:

Now we will add in the energy functionals, including linecurvsq and linetension terms that penalize the curvature squared and length respectively. We’ll also establish a global constraint, namely fixing the area enclosed by the loop.

// Line curvature squared
lc=new(linecurvsq)
lc.setprefactor(0.01)
m.addenergy(lc)

// Add line tension
lt=new(linetension)
m.addenergy(lt)

// Keep the area enclosed constant
la=new(enclosedarea)
m.addconstraint(la)

//Promote equal sized segments
le=new(equilength);

Notice here we’re also creating a functional equilength that we’ll use to regularize the mesh.

Now for the level set constraints. Let’s define the first step by step. First, create a levelset object:

llower=new(levelset);

We’ll next set up the coordinate system for the object so that we can specify the level set using these symbols

llower.setcoordinates({x,y,z});

Note that this only affects the levelset object and not the global coordinate system. It’s also important to know that, for now, Morpho only understands Cartesian coordinates (we expect to add other coordinate systems in future versions).

We now specify the level set function definition, here $y+y_0=0$. Only the left hand side of this equation needs to be given, the right hand side is assumes to be zero. This particular level set is a horizontal plane at $y=-y_0$. The gradient of the level set function also has to be specified, which is a vector in 3D space:

llower.setexpression(y+y0);
llower.setgradient({0,1,0});

Finally, we are going to make this constraint one sided, which means that it will only act to prevent mesh points from going below the $y=-y0$ plane.

llower.setonesided(true);

Notice that in the allowed region, the constraint function $y+y_0$ is positive while in the forbidden region it’s negative; this means you can always exchange the allowed and forbidden regions by changing the sign of the constraint function.

Finally, we tell the body to use this as a local constraint.

m.addlocalconstraint(llower)

Exercise 1. Make a second level set constraint, lupper, that prevents the body from going above the $y=y_0$ plane. [Hint: The constraint function should be positive for $y<y_0$]

[code lang=”js” gutter=”false” collapse=”true” title=”Solution”] lupper=new(levelset); lupper.setcoordinates({x,y,z}); lupper.setexpression(-y+y0); lupper.setgradient({0,-1,0}); lupper.setonesided(true); m.addlocalconstraint(lupper) [/code]

Now that we’re all done with constraints, we can create a function that performs a specified number of minimization and regularization steps and allows us to track the values of each functional:

// Minimize function
minimize(n) = {
	do(
		m.linesearch(scale=0.1, scalelimit=sl);
		m.linesearch(le, scale=0.1, scalelimit=sl);
		print({i, m.totalenergy(), m.evaluatefunctional(lc), m.evaluatefunctional(lt), m.evaluatefunctional(la), m.evaluatefunctional(le)});
	, {i, n})
};

And that’s it! If you run this code, you can now perform the minimization just by calling minimize(). We have initially set $y_0$ to be 0.9, so when we run the minimization, e.g. by typing

minimize(Nsteps)

you should see the loop deform into something like the following

You can modify $y_0$ by typing

y0=0.7

and then after reminimizing you’ll see

Exercise 2. Try minimizing with a few different values of $y_0$. What happens to the shape? What happens if you start with a fresh loop and go straight to a very small value of $y_0$ like 0.1?

You should find that if you start with $y_0$ too small, the minimization becomes more challenging. This is because the initial loop violates the level set constraints (it has points in the forbidden region). Morpho tries to correct for this by moving mesh points into the allowed region, but the further they are from the boundary the more challenging this is. Generally, it’s a good idea to do something like the above example, starting with the level set constraint just perturbing the body and then changing its position slowly until you get the desired final configuration.

That’s it for this tutorial. You should now be able to create and impose level constraints in your own code—and we’ll be using them in several later tutorials.