# Minimal surface with obstacle problem
# Michael Merritt - Summer 2000
# Version 2.0 - October 2000

model;

param nx > 0, integer; # grid points in 1st direction
param ny > 0, integer; # grid points in 2nd direction

param hx := 1/(nx+1);     # grid spacing
param hy := 1/(ny+1);     # grid spacing
param area := 0.5*hx*hy;  # area of triangle

# v defines the finite element approximation.

var v {i in 0..nx+1,0..ny+1};

minimize surf:
  (sum {i in 0..nx,j in 0..ny}
    area*(1+((v[i+1,j]-v[i,j])/hx)^2+((v[i,j+1]-v[i,j])/hy)^2)^(1/2)) +
  (sum {i in 1..nx+1,j in 1..ny+1}
    area*(1+((v[i-1,j]-v[i,j])/hx)^2+((v[i,j-1]-v[i,j])/hy)^2)^(1/2));

subject to b1 {j in 0..ny+1}: v[0,j] = 0;
subject to b2 {j in 0..ny+1}: v[nx+1,j] = 0;
subject to b3 {i in 0..nx+1}: v[i,0] = 1 - (2*i*hx-1)^2;
subject to b4 {i in 0..nx+1}: v[i,ny+1] = 1 - (2*i*hx-1)^2;
subject to b5 {i in 0..nx+1,j in 0..ny+1}: v[i,j] >= 0;
subject to b6 {i in floor(0.25/hx)..ceil(0.75/hx),
               j in floor(0.25/hy)..ceil(0.75/hy)}: v[i,j] >= 1;

# Minimal surface with obstacle problem
# Michael Merritt - Summer 2000
# Version 2.0 - October 2000

data;

# Set the number of discretization points and design parameters

param nx := 200;
param ny := 200;

# Set the starting point.

let {i in 0..nx+1,j in 0..ny+1} v[i,j]:= 1 - (2*i*hx-1)^2;