Dan Wood

Best I can think of is something like this - set all the border points to have a targetmeshsize that roughly matches the local length of the longest polygon edge, and all the interior points to something smaller, then use Adaptive edge lengths mode, Gradation of 1.0, and switch on the Mesh Size Attribute control attrib.

This won't be a guarantee that it'll avoid splitting border edges, but you should be able to "massage" the values to have the effect you want.

(Edit: found doing a 2-pass thing with some inward value spreading to yield a bit cleaner result while letting you push the inner target size a little lower)

(Edit 2: border_preserve_remesh_03.hiplc - or... one with a guaranteed border-preserving tidy-up step)

I'm far from sure what the limits of the FLIP solver's capabilities are, but I presume the answer is going to be something along the lines of "it's not built to model it in a physically-correct way, but there are ways to fake it".

You can simply set f@density to 1 on air particles, 1000 on water particles, fill an entire container with particles, and switch on variable density... it'll sort-of work - the broad motion will look correct (at least for a couple of seconds), but any cells that end up with a mix of particle densities will be unable to separate those particles correctly (those voxels will just be seen by the pressure solve as a homogeneous fluid with a density the average of the particles it contains), and so droplets of water that should be treated as ballistic flying through the air will end up floating instead, and the water will quickly just "diffuse" into the air.

I've had some success in the past running a secondary smoke-solve alongside a FLIP solve, using the FLIP water surface/velocity as a collider in the smoke solve, and then using the smoke solve's velocity to advect near-surface/droplet particles, and especially white-water secondary spray particles (at least I think that's how I did it, from memory)... that's possibly the most practical way to do it using the tools available. Doesn't help for sub-surface incompressible pockets, but it adds life to large splashy motion.

I wonder if there's a way to run a custom separate pass that will work out the local cell fraction at the fluid surface, and push the different density particles to the correct side of the boundary. I doubt that would be enough on its own though.
I get the feeling others out there have probably gone further down this rabbit-hole (and/or actually understand the maths behind two-phase fluids). Would love to know more if anyone's got it working in a remotely physically-correct way.

...or we just wait for SideFX to implement that "Adaptive Phase-Field-FLIP for Very Large Scale Two-Phase Fluid Simulation" paper from Siggraph. I'm hoping it'll come along eventually (nudge nudge)

Ah, no, the big initial explosions are intentional - I left some gaps in the fluid for it to crash into, so I could stress test it with some reasonably violent motion.

I'm not using any of Houdini's default feedback here - the FLIP and Bullet solvers are completely decoupled (while still running in parallel in the same DOP loop), and I built my own feedback system to calculate the buoyancy, which is now using accurate physical values - the bodies are set to 50 density in that test, and it's using the FLIP output pressure field 1:1 to generate per-face forces, and then uses a script solver to edit the uniform force attached to each Bullet body. On the FLIP side, I'm importing the rigid body transform data and using it to generate the collision/collisionvel fields directly from high res meshes (while the bullet bodies themselves are using compound capsule colliders).