Classifying Constraints

There are several schemes of classifying constraints. Here, we will consider those that have the most implications for computational implementation.

Based on the geometric entity concerned, we can have point (isolated), distributed, and global constraints. For example, on one hand, the boundary conditions in a 1D problem are constraints at isolated points. On the other hand, a condition that should hold at every point is a distributed constraint. Global constraints specify some norm (usually an integral) of the solution. For example, specifying the length of a catenary cable or surface area of a soap film provides a global constraint.

Some people use the term pointwise constraint for a distributed constraint. We would like to make a clear distinction between that and what we are calling a point constraint here. A point constraint is enforced at a single point or a finite number of isolated points. This set of points has no length, area, or volume. Distributed constraints, however, hold at every point of a region. This can be every point of an edge, surface, or domain of a 3D object.

Another classification is equality constraints versus inequality constraints. A familiar inequality constraint in structural mechanics arises in contact mechanics. The gap between contacting objects in an assembly has to be nonnegative. In chemical reaction engineering, lower bounds on species concentrations are also inequalities.

These classifications overlap. For instance, we can have distributed inequality constraints and distributed equality constraints and so on. Inequality constraints are mathematically more challenging, so we will focus on equality constraints first and move on to the former a bit later in the series.

Brief Introduction to the Theory of Equality Constraints

The calculus problem

is solved by finding the stationary points of the augmented objective function

with respect to both the coordinate  and the Lagrange multiplier .

This same idea is used, albeit with appropriate extensions, in variational calculus. Consider the problem

Here, the distributed constraint  has to be satisfied at all points in our domain and not just at one point. Thus, each point will have its own Lagrange multiplier, making  a function and not just one number. Accordingly, the augmented functional is

(1)

Having turned the constrained variational problem for one field, , to an unconstrained problem in two fields,  and , we turn our attention to the optimality criteria. It is important to note here that we can make independent variation in both the solution field and the Lagrange multiplier field. That gives us the first-order optimality criteria

In terms of  and , we have

These are the equations we need to enter in the Weak Form PDE interface. We will quickly return to that. First, let’s derive the corresponding conditions for global (integral) and point (isolated) constraints.

When we have a global constraint

the augmented functional is

where  is one number and not a field.

The first-order optimality conditions are

Finally, using the properties of the Dirac delta function, the point (isolated) constraint  can be thought of as the global constraint

Dirac delta functions can be analogously utilized to derive constraints on edges in 2D and 3D, and surfaces in 3D.

Plugging the above result in the formulation for global constraints, we get the first-order optimality conditions

(2)

(3)

If we have more than one isolated point constraint, we will have a Lagrange multiplier for each point. The Lagrange multiplier will not be a field, but a finite set of scalars, one valid at each isolated point. If we have a distributed constraint that is not imposed on the whole domain but over parts of a domain, we can define the Lagrange multiplier only over that part. Here, we will have a Lagrange multiplier field.

Implementing Constraints in COMSOL Multiphysics®

Let us now see how to implement constraints in COMSOL Multiphysics. We will consider the same soap film problem from the previous blog post, but with the following boundary conditions.

The first boundary condition is something we could have specified using the Dirichlet Boundary Condition node, but for pedagogic reasons, we will use the more general constraint framework. The two boundary conditions above can be rewritten as

We need the partial derivatives of the constraint equations with respect to  and .

Finally, we will plug these as weak contributions at the corresponding points. The contribution from  remains the same as before.

Specifying point constraints using point weak contributions.

Note that the point contribution from (Eq. 3) and the second term in (Eq. 2) have been added together in the Weak Contribution 1. The logic here is that, since  and  are independent variations, setting the sum of terms containing those variations to zero is equivalent to setting terms containing each variation to zero.

Solution with radius 2 at the left end and zero slope at the right end.

In the theory section, we discussed different types of constraints. The table below summarizes the recommended places in the software where contributions from constraints are specified and unknowns (Lagrange multipliers) are defined based on the type of constraint.

Constraint Type	Examples	Constraint Contribution Containing	Constraint Contribution Containing	Where to Define
Distributed constraint	Constraints on material behavior Lower or upper bounds on solutions Boundary conditions in 2D and 3D problems	Weak contribution	Weak contribution	Auxiliary dependent variable Attribute under weak contribution
Global (integral) constraint	Length of an edge Weight of an object Average displacement (temperature) on a surface	Weak contribution	Global Equations node	Global Equationsnode
Isolated point constraint	Boundary conditions in 1D Constraints on an isolated internal point	Weak contribution	Weak contribution	Auxiliary dependent variable Attribute under weak contribution

Specifying Fluxes (Forces)

So far, we discussed specifying constraints. To this end, we introduced unknown Lagrange multipliers. In many physical problems, Lagrange multipliers are reaction forces or fluxes necessary to enforce a constraint. If instead of the constraint, you know the applied forces or fluxes, how do you specify that?

Go back to the functional to be minimized and reformulate, including the forces (fluxes). For example, for boundary loads in structural mechanics, this adds the virtual work due to the loads. The expression to enter in the COMSOL Multiphysics interface is similar to what we had to do for constraints. The known value of the force (flux) goes in place of the Lagrange multiplier, and we do not need auxiliary variables for Lagrange multipliers. As a result, the term containing the variation of the Lagrange multiplier disappears.

Find more information in our previous blog post on how to add point loads (sources) using weak contributions.

Adding Constraints to Points That Are Not in the Geometry Sequence

Occasionally, we want to add a constraint to a point in our domain that is not explicitly inserted (or created by intersection of curves) in the geometry sequence. We cannot select the point and associate it with the point weak contribution as shown in the above example. However, we can add a global weak contribution and use a domain point probe to refer to the solution and its variation.

Consider a catenary cable supported at two ends. For a cable with uniform weight per unit length the variational problem is the same as the axisymmetric soap film problem. Thus, we will use the same COMSOL model. Additionally, we want to constrain the height  to 1.95 at the center without adding a point node to the geometry sequence. First, we add a domain point probe for . Let the name of the probe be ppb1. The weak contribution corresponding to (Eq. 1) is lam_c*test(comp1.ppb1), where lam_c is a new Lagrange multiplier. We enter this contribution in the global Weak Contribution node as shown in the screenshot below. This node does not allow the creation of auxiliary variables, as opposed to weak contributions on explicitly defined geometric entities. However, we can add a Global Equation node where we can define the Lagrange multiplier and specify the constraint.

Adding a constraint on a point not in the geometry.

The plots below show the solution including that constraint. Circles are added at mesh nodes. The plot on the left shows the case where the central point does not have a mesh node associated with it. As such, the constraint is satisfied only approximately. On the right, we see the case of a solution on a finer mesh, where the mesher added a node where we want to impose a constraint. Here, the constraint is satisfied exactly.

Verifying internal constraints on points not in the geometry sequence. Aspect ratio is not preserved in plots.

A similar strategy can be used to add point loads (sources) for problems such as structural mechanics, heat transfer, or chemical transport. As mentioned earlier, such loads mathematically correspond to Lagrange multipliers. If we know the Lagrange multiplier, we do not add the Global Equation. Rather, we just have the weak contribution and, in place of lam_c above, we type the applied mechanical, thermal, chemical, or other type of load.

Stay Tuned!

So far in this blog series, we have shown how to solve variational problems using the Weak Form PDE interface and how to include equality constraints.

Constraints introduce additional complexity in the numerical solution. Mainly, they lead to saddle point problems, and a Lagrange multiplier implementation destroys the positive definiteness of the stiffness matrix. For nonlinear constraints, there is the additional danger of singular matrices during the nonlinear iteration. This can be especially problematic in global and distributed constraints.

In the next installment of this series, we will show numerical strategies for circumventing these issues. After we discuss these strategies using equality constraints, we will proceed to inequality constraints later in the series.