The code shown is for packing identical circles inside a Cassini oval. Results are shown for that curve and also for a cardioid. For small numbers of circles, best result were obtained using NMaximize. For larger numbers of circles, best results were obtained using ParametricIPOPTMinimize. The code to implement that has been previously posted.

 Packing Code

 formula for Cassini Oval

In[21]:= oval[{x_, y_}] = ((x - 1)^2 + y^2) ((x + 1)^2 + y^2) - (21/20)^4;

 plot of Cassini Oval

In[22]:= ovalplot = 
  ContourPlot[Evaluate[oval[{x, y}] == 0], {x, -1.5, 1.5}, {y, -1.5, 1.5}, 
   ImageSize -> Small];

 signed distance from a point {xc,yc} to the oval

In[23]:= dist[{xc_, yc_}] := -N @Sign[oval[{xc, yc}]]* 
  Sqrt[MinValue[{(x - xc)^2 + (y - yc)^2, oval[{x, y}] == 0}, {x, y}]]

 evaluate the distance at an array of points

In[24]:= AbsoluteTiming[
 Quiet[disttbl = 
    Flatten[Table[{xc, yc, dist[{xc, yc}]}, {xc, -2, 2, 1/10}, {yc, -1, 1, 1/
       10}], 1]];]

Out[24]= {14.3459, Null}

 create an Interpolation function for the distance

In[25]:= disti = Interpolation[disttbl];

 constraints that two circles be disjoint

In[26]:= disjointcircs[i_Integer, 
  j_Integer] := (xc[i] - xc[j])^2 + (yc[i] - yc[j])^2 >= 4 radius^2

 maximize the radii of the circles, constrained to be in the oval, using NMaximize

In[27]:= pack[n_Integer] := 
 Quiet[NMaximize[{radius, 
    Sequence @@ Table[disti[xc[i], yc[i]] >= radius, {i, n}] , 
    Sequence @@ 
     Flatten[Table[disjointcircs[i, j], {i, n - 1}, {j, i + 1, n}]]}, 
   Join[{radius}, Table[xc[i], {i, n}], Table[yc[i], {i, n}]]]]

plot the result 

In[28]:= plt[sln_, n_Integer] := 
 Show[ovalplot, 
  Graphics @ Table[Circle[{xc[i], yc[i]}, radius] /. sln[[2]], {i, n}]]

 perform the optimization using Ipopt

In[29]:= iMinPack[n_Integer, nrands_Integer, seed_Integer] := 
 iMin[-radius, {Sequence @@ Table[disti[xc[i], yc[i]] >= radius, {i, n}] , 
   Sequence @@ 
    Flatten[Table[disjointcircs[i, j], {i, n - 1}, {j, i + 1, n}]]}, 
  Join[{{radius, 0, 1}}, Table[{xc[i], -2, 2}, {i, n}], 
   Table[{yc[i], -1, 1}, {i, n}]], nrands, seed, 
  "IPOPTOptions" -> {"tol" -> 10^-6}]