How to fill a planar 3D path in asymptote?

Question

Is there, in Asymptote an equivalent of the fill command that works for planar 3D paths?

I thought surface did this: I expected the commands

fill(project(p))

and

draw(surface(p))

to always produce identical results, for any planar path3 p. However, this does not work for the path I am working on (called tile(identity4) in the code below): basically it is some nonconvex shape, and the 3D surface contains parts of its convex hull not enclosed by the path.

I did not manage to produce a simple MWE. For what it's worth, here is my code; I'm afraid it is a bit complicated to understand in detail! I removed every definition irrelevant to the problem; but when I tried to simplify it further, it suddently stopped compiling in a reasonable time, so I dropped the matter.

You might ask why I do not simply use the former command. The problem is that I have other data in my picture, and it seems that when I draw a 2D picture, all previously drawn 3D objects get erased.

unitsize(100);
import three;
import graph3;
import solids;

real tmax = 1.15; //contraction strength
real translation_speed = 0.2; //ratio of translation part over contraction part of an infinitesimal generator
real height = 1; //height of the cutting plane
real radius = 5; //radius of the "projective sphere
real width = 2; //how far do the "infinite" surfaces actually go
int subdiv = 1; //how much we subdivide curves when projecting to the sphere

transform3 squish(real t) {
  return scale(exp(-t), 1, exp(t)) * shift(t*translation_speed*Y);
}
transform3 halfgplus = squish(tmax);
transform3 halfgminus = squish(-tmax);
transform3 gplus = squish(2*tmax);
transform3 gminus = squish(-2*tmax);
transform3 conj = rotate(90,X+Z);
transform3 halfhplus = conj * halfgplus * inverse(conj);
transform3 halfhminus = conj * halfgminus * inverse(conj);
transform3 hplus = conj * gplus * inverse(conj);
transform3 hminus = conj * gminus * inverse(conj);

triple c1 = (1,0,1);
triple c3 = (1,0,-1);
triple c5 = (-1,0,-1);
triple c7 = (-1,0,1);

path3 tennisball = scale3(width) *
(arc((1,0,0),c1,c3,-X) &
arc((0,0,-1),c3,c5,-Z) &
arc((-1,0,0),c5,c7,X) &
arc((0,0,1),c7,c1,Z));

path3[] cut(path3 p, triple center) {
//Cuts the path p by the plane x+z = height, then projects the upper part onto the plane.
//Returns {flattened upper part, lower part}.
//Assumes exactly 2 intersections.
  path p_projected = project(p, orthographic(1,0,-1));
  real [] cuttimes = times(p_projected, (0, height));
  real t1 = cuttimes[0];
  real t2 = cuttimes[1];
  pair test_point = point(p_projected, (t1+t2)/2);
  path3 upper_part = (test_point.y >= height) ? subpath(p, t1, t2) :
                                                subpath(p, t2, length(p)) & subpath(p, 0, t1);
  path3 lower_part = (test_point.y >= height) ? subpath(p, t2, length(p)) & subpath(p, 0, t1) :
                                                subpath(p, t1, t2);
  triple projection_target = planeproject(X+Z, sqrt(2)*height*Z) * center;
  path plane_upper_part = project(upper_part, perspective(center, target=projection_target));
  path3 flat_upper_part = invert(plane_upper_part, X+Z, sqrt(2)*height*Z, perspective(center, target=projection_target));
  return new path3[] {flat_upper_part, lower_part};
}

triple sph_proj(triple x, triple proj_cen) {
//Projects a point from some center (assumed to be inside the sphere) onto the sphere of center O and given radius.
  triple v = unit(x - proj_cen); //Unit vector generating the half-line.
  triple q = proj_cen - dot(v, proj_cen) * v; //Point of the line closest to O.
  return q + (radius +-+ abs(q)) * v;
}

path3 sph_proj_path(path3 p, triple proj_cen) {
//Projects a smooth, cyclic, non null path from some center (assumed to be inside the sphere) onto the sphere of center O and given radius.
  int n=length(p);
  guide3 g=sph_proj(point(p,0), proj_cen);
  for(int i=1; i < n*subdiv; ++i)
    g=g..sph_proj(point(p,i/subdiv), proj_cen);
  return g..cycle;
}

path3 side(transform3 gamma) {
  triple cen = gamma * O;
  return cut(sph_proj_path((gamma * tennisball), cen), cen)[0];
}

path3 tile(transform3 gamma) {
  path3 side_gplus = side(gamma * halfgplus);
  path3 side_gminus = side(gamma * halfgminus);
  path3 side_hplus = side(gamma * halfhplus * conj);
  path3 side_hminus = side(gamma * halfhminus * conj);
  path3 tile_contour = reverse(side_gplus) -- reverse(side_hplus) -- side_gminus -- side_hminus -- cycle;
  return planeproject(X+Z, sqrt(2)*height*Z) * tile_contour;
}

draw(surface(tile(identity4)), gray+opacity(0.8));

Update: I also tried

draw(surface(project(p, orthographic(n)), plane=myplane))

where myplane is the plane where I would like to draw the surface and n is its normal. This produces yet a different result, but still not the one I expect.

Answer 1

With Asymptore 3D surfaces are Bézier patches. A planar path3 does not define such a surface. The only solution I know is to use bezulate routines which transform a 2D path (cyclic or union of cyclic path) into a Bézier triangulation (2D). Then you can draw it with surface(bezulate(your 2D path))). It was developed to have 3D label.

I had some problems with your code (projections choice) and in particular for halfgminus and halfhminus the bezulate routine returns an error message. Here a part of your picture

path[]  aa=  project(side(identity4*halfgplus)--cycle ^^reverse(side(identity4*halfhplus*conj)--cycle),
               orthographic(X+Z));
path[] gg=bezulate(aa);
draw(planeproject(X+Z)*surface(gg),opacity(0.4));