# obstacle.mod # Calculate the position of a membrane pushed up through a (rectangular) # hole in a rigid plate; in addition, there are rigid obstacle(s) inside # the hole (perhaps not at the same level as the plate). These # obstacles can constrain the membrane from above or from below. The # correct position of the membrane (the position of minimum energy) is # determined by the minimization of a quadratic function of the membrane # position, subject to the constraints imposed by the hole and the # obstacle. The MCP below arises from the optimality conditions for # this QP. # Reference: Ciarlet, Philippe G., "The Finite Element Method for # Elliptic Problems", North-Holland, 1978. param M integer >= 1, default 50; # of interior grid pts in Y direction param N integer >= 1, default 50; # of interior grid pts in X direction set Y := 0 .. M+1 ; set X := 0 .. N+1 ; param xlo := 0; param xhi > xlo, := 1; param ylo := 0; param yhi > ylo, := 1; param dy := (yhi - ylo) / (M + 1) ; param dx := (xhi - xlo) / (N + 1) ; param c := 1; /* force constant */ param ub {i in Y, j in X} := if 1 <= i <= M and 1 <= j <= N then (sin(9.2*(xlo+dx*i))*sin(9.3*(ylo+j*dy)))^2 + 0.2; param lb {i in Y, j in X} := if 1 <= i <= M and 1 <= j <= N then (sin(9.2*(xlo+dx*i))*sin(9.3*(ylo+j*dy)))^3; # height of membrane var v {i in Y, j in X} >= lb[i,j] <= ub[i,j] := max(0,lb[i,j]); s.t. dv {i in 1..M, j in 1..N}: lb[i,j] <= v[i,j] <= ub[i,j] complements (dy/dx) * (2*v[i,j] - v[i+1,j] - v[i-1,j]) + (dx/dy) * (2*v[i,j] - v[i,j+1] - v[i,j-1]) - c * dx * dy ;