# How to fill a planar 3D path in asymptote?

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 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;
real t2 = cuttimes;
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);
}

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.

``````path[]  aa=  project(side(identity4*halfgplus)--cycle ^^reverse(side(identity4*halfhplus*conj)--cycle),