#  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 ;