obstacle.gms : An Obstacle Problem

**Description**

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.

**Reference**

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

**Small Model of Type :** MCP

**Category :** GAMS Model library

**Main file :** obstacle.gms

```
$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
*
* On the major platforms, GMSPYTHONHOME gets set automatically, otherwise the
* user has to set it. This condition can also be removed, if one has set up its
* Python environment appropriately
$if not setEnv GMSPYTHONHOME
$if %gams.pySetup%==1 $abort.noError Embedded code Python not ready to be used
$if not setEnv GMSPYTHONHOME $log --- Using external Python
$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:
mi.sync_db[sym].copy_symbol(gams.db[sym])
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;
);
```