$title An Obstacle Problem

$onText
Program to 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.

The MCP is solved for multiple force constants c using GUSS and
GAMSModelInstance. This is the first model in the library of this
model type that is solved this way.

Ciarlet, P G, The Finite Element Method for Elliptic Problems.
Elsevier Science, Studies in Mathematics and its Applications, 1978.
$offText

Sets
   I       / 0 * 11 /
   J       / 0 * 11 /
   boundary(I,J)
   interior(I,J)
   ubnd(I,J) 'points with an upper bound';

boundary(I,J)$[(ord(J) eq card(J)) or (ord(J) eq 1) or
               (ord(I) eq card(I)) or (ord(I) eq 1)] = yes;
interior(I,J) = not boundary(I,J);

Scalars
   xlo     / 0 /
   xhi     / 1 /
   ylo     / 0 /
   yhi     / 1 /
   dx, dy
   c   'force constant' / 1 /
   M   'number of divisions in I'
   N   'number of divisions in J';

M = (card(I)-1);  
N = (card(J)-1);  
dx = (xhi - xlo) / M;
dy = (yhi - ylo) / N;

Parameters
   d(I,J)   'distance of point from center of rectangle'
   tmp(I,J)
   ub(I,J);

Variable v(I,J)  'height of membrane';

Equation dv(I,J);

dv(I,J) ..
  (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))
  =N=
  c * dx * dy;

Model obstacle / dv.v /;

* problem with a ball constraining from above
Scalar
  midI, midJ
  cheight  'height of center of ball' / .55 /
  ball_rad 'ball radius'              / 0.5 /
  ;
midI = 1 + M / 2;
midJ = (card(J) + 1) / 2;
d(I,J) = sqrt(sqr((ord(I)-midI)/M) + sqr((ord(J)-midJ)/N));

ubnd(I,J) = [d(I,J) <= ball_rad];
v.up(ubnd(I,J)) = cheight - sqrt(ball_rad-sqr(d(I,J)));

* the position v is fixed at zero on the boundary
v.fx(boundary) = 0;

v.l(interior(I,J)) = min(cheight, v.up(I,J));

Set s 'force scenarios' / s1*s5 /;
Parameter cval(s);

cval(s) = (ord(s)-1)*(1.5-0.5)/(card(s)-1) + 0.5;
display cval;

Set mattrib / system.GUSSModelAttributes /;

Parameter
   vs(s,I,J)
   srep(s, mattrib)  'model attibutes like modelstat etc'
   o(*)              'GUSS options' / SkipBaseCase 1 /;

Set obstdict / s.scenario.'', o.opt.srep, c.param.cval, v.level.vs /;

solve obstacle using mcp scenario obstdict;

option vs:4:1:1; display vs;

* Now do the same with a GAMSModelInstance

$log --- Using Python library %sysEnv.GMSPYTHONLIB%

$set solverlog
$if set useSolverLog $set solverlog output=sys.stdout

embeddedCode Python:
def solveMI(mi, symIn=[], symOut=[]):
  for sym in symIn:
    gams.db[sym].copy_symbol(mi.sync_db[sym])
  mi.solve(%solverlog%)
  for sym in symOut:
    try:
      gams.db[sym].clear() # Explicitly clear the symbol to ensure setting "writtenTo" flag for sym
      mi.sync_db[sym].copy_symbol(gams.db[sym])
    except:
      pass
gams.wsWorkingDir = "."    
pauseEmbeddedCode
abort$execerror 'Python error. Check the log';

$libInclude pyEmbMI miObst 'obstacle using mcp' -all_model_types=path c.Zero

option v:4:1:1;
loop(s,
   c = cval(s);
   continueEmbeddedCode:
   gams.db['c'].copy_symbol(miObst.sync_db['c'])
   solveMI(miObst,['c'],['v'])
   pauseEmbeddedCode v
   display v.l;
);