Voxillization creates errant geometry

[FCLC]

Hi Moritz,

Following up on the post from mastodon, I wanted to point out a potential bug. Specifically during the voxel creation phase when converting an STL mesh, the converter seems to be creating "fake" geometry.

Example of behaviour:

I'm running the latest version

The file in question is:

RB18_alligned.zip

And the code being run is:

setup.zip

#include "setup.hpp"
#ifndef BENCHMARK

void main_setup() { // F1 car
	// ######################################################### define simulation box size, viscosity and volume force ############################################################################
	const uint L = 672u; // 2152u on 8x MI200
	const float kmh = 100.0f;
	const float si_u = kmh/3.6f;
	const float si_x = 2.0f;
	const float si_rho = 1.225f;
	const float si_nu = 1.48E-5f;
	const float Re = units.si_Re(si_x, si_u, si_nu);
	print_info("Re = "+to_string(Re));
	const float u = 0.08f;
	const float size = 1.6f*(float)L;
	units.set_m_kg_s(size*2.0f/5.5f, u, 1.0f, si_x, si_u, si_rho);
	const float nu = units.nu(si_nu);
	print_info("1s = "+to_string(units.t(1.0f)));
	LBM lbm(0.6*L, 1.65*L, 0.35*L, nu); // (width, depth, height)
//	LBM lbm(L, 2*L, 3*L, nu); //work with L=256u, u=0.08f
	// #############################################################################################################################################################################################
	const float3 center = float3(lbm.center().x, lbm.center().y, lbm.center().z);
//	const float3 center = float3(lbm.center().x, 0.525f*size, 0.116f*size);
	lbm.voxelize_stl("/home/felix/cad/RB18_alligned.stl", center, size); // https://addons-redbullracing-com2020.redbull.com/b3993c8955aad441ec69/assets/cars/RB18/model/RB_18_v05.glb
	const ulong N=lbm.get_N(); const uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); for(ulong n=0ull; n<N; n++) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
		// ########################################################################### define geometry #############################################################################################
		//if(lbm.flags[n]!=TYPE_S) lbm.u.y[n] = u;
		if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==Nz-1u) lbm.flags[n] = TYPE_E;
		if(z==0u) lbm.flags[n] = TYPE_S;
		const float3 p = lbm.position(x, y, z);
		const float W = 1.05f*(0.312465f-0.179692f)*(float)Nx;
		const float R = 1.05f*0.5f*(0.361372f-0.255851f)*(float)Ny;
//rotating bodies

/*		const float3 FL = float3(0.247597f*(float)Nx, -0.308710f*(float)Ny, -0.260423f*(float)Nz);
		const float3 HL = float3(0.224739f*(float)Nx, 0.210758f*(float)Ny, -0.264461f*(float)Nz);
		const float3 FR = float3(-FL.x, FL.y, FL.z);
		const float3 HR = float3(-HL.x, HL.y, HL.z);
		if((lbm.flags[n]&TYPE_S) && cylinder(x, y, z, lbm.center()+FL, float3(W, 0.0f, 0.0f), R)) {
			const float3 uW = u/R*float3(0.0f, FL.z-p.z, p.y-FL.y);
			lbm.u.y[n] = uW.y;
			lbm.u.z[n] = uW.z;
		}
		if((lbm.flags[n]&TYPE_S) && cylinder(x, y, z, lbm.center()+HL, float3(W, 0.0f, 0.0f), R)) {
			const float3 uW = u/R*float3(0.0f, HL.z-p.z, p.y-HL.y);
			lbm.u.y[n] = uW.y;
			lbm.u.z[n] = uW.z;
		}
		if((lbm.flags[n]&TYPE_S) && cylinder(x, y, z, lbm.center()+FR, float3(W, 0.0f, 0.0f), R)) {
			const float3 uW = u/R*float3(0.0f, FR.z-p.z, p.y-FR.y);
			lbm.u.y[n] = uW.y;
			lbm.u.z[n] = uW.z;
		}
		if((lbm.flags[n]&TYPE_S) && cylinder(x, y, z, lbm.center()+HR, float3(W, 0.0f, 0.0f), R)) {
			const float3 uW = u/R*float3(0.0f, HR.z-p.z, p.y-HR.y);
			lbm.u.y[n] = uW.y;
			lbm.u.z[n] = uW.z;
		}
*/	}	// #########################################################################################################################################################################################
//	key_4 = true;
	//Clock clock;
	//lbm.run(0u);
	//while(lbm.get_t()<=units.t(1.0f)) {
	//	lbm.graphics.set_camera_free(float3(0.779346f*(float)Nx, -0.315650f*(float)Ny, 0.329444f*(float)Nz), -27.0f, 19.0f, 100.0f);
	//	lbm.graphics.write_frame_png(get_exe_path()+"export/a/");
	//	lbm.graphics.set_camera_free(float3(0.556877f*(float)Nx, 0.228191f*(float)Ny, 1.159613f*(float)Nz), 19.0f, 53.0f, 100.0f);
	//	lbm.graphics.write_frame_png(get_exe_path()+"export/b/");
	//	lbm.graphics.set_camera_free(float3(0.220650f*(float)Nx, -0.589529f*(float)Ny, 0.085407f*(float)Nz), -72.0f, 21.0f, 86.0f);
	//	lbm.graphics.write_frame_png(get_exe_path()+"export/c/");
	//	lbm.run(units.t(0.5f/600.0f)); // run LBM in parallel while CPU is voxelizing the next frame
	//}
	//write_file(get_exe_path()+"time.txt", print_time(clock.stop()));
	lbm.run();
} /**/



#endif // SURFACE
#ifdef TEMPERATURE



/*void main_setup() { // Rayleigh-Benard convection
	// ######################################################### define simulation box size, viscosity and volume force ############################################################################
	LBM lbm(256u, 256u, 64u, 0.02f, 0.0f, 0.0f, -0.001f, 0.0f, 1.0f, 1.0f);
	// #############################################################################################################################################################################################
	const ulong N=lbm.get_N(); const uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); for(ulong n=0ull; n<N; n++) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
		// ########################################################################### define geometry #############################################################################################
		lbm.u.x[n] = random_symmetric(0.015f);
		lbm.u.y[n] = random_symmetric(0.015f);
		lbm.u.z[n] = random_symmetric(0.015f);
		if(z==1u) {
			lbm.T[n] = 1.75f;
			lbm.flags[n] = TYPE_T;
		} else if(z==Nz-2u) {
			lbm.T[n] = 0.25f;
			lbm.flags[n] = TYPE_T;
		}
		//if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) lbm.flags[n] = TYPE_S; // all non periodic
		if(z==0u||z==Nz-1u) lbm.flags[n] = TYPE_S;
	}	// #########################################################################################################################################################################################
	lbm.run();
} /**/



/*void main_setup() { // TEMPERATURE test
	// ######################################################### define simulation box size, viscosity and volume force ############################################################################
	LBM lbm(32u, 196u, 60u, 1u, 1u, 1u, 0.02f, 0.0f, 0.0f, -0.001f, 0.0f, 1.0f, 1.0f);
	// #############################################################################################################################################################################################
	const ulong N=lbm.get_N(); const uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); for(ulong n=0ull; n<N; n++) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
		// ########################################################################### define geometry #############################################################################################
		if(y==1) {
			lbm.T[n] = 1.8f;
			lbm.flags[n] = TYPE_T;
		} else if(y==Ny-2) {
			lbm.T[n] = 0.3f;
			lbm.flags[n] = TYPE_T;
		}
		if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) lbm.flags[n] = TYPE_S; // all non periodic
	}	// #########################################################################################################################################################################################
	lbm.run();
	//lbm.run(1000u); lbm.u.read_from_device(); println(lbm.u.x[lbm.index(Nx/2u, Ny/2u, Nz/2u)]); wait(); // test for binary identity
} /**/


/*
#endif // TEMPERATURE
#else // BENCHMARK
#include "info.hpp"
void main_setup() { // benchmark
	uint mlups = 0u;
	{ // ######################################################## define simulation box size, viscosity and volume force ###########################################################################
		//LBM lbm( 32u,  32u,  32u, 1.0f);
		//LBM lbm( 48u,  48u,  48u, 1.0f);
		//LBM lbm( 64u,  64u,  64u, 1.0f);
		//LBM lbm( 96u,  96u,  96u, 1.0f);
		//LBM lbm(128u, 128u, 128u, 1.0f);
		//LBM lbm(192u, 192u, 192u, 1.0f);
		LBM lbm(256u, 256u, 256u, 1.0f);
		//LBM lbm(384u, 384u, 384u, 1.0f);
		//LBM lbm(464u, 464u, 464u, 1.0f);
		//LBM lbm(480u, 480u, 480u, 1.0f);
		//LBM lbm(512u, 512u, 512u, 1.0f);

		//const uint memory = 4096u; // in MB
		//const uint L = ((uint)cbrt(fmin((float)memory*1048576.0f/(19.0f*(float)sizeof(fpxx)+17.0f), (float)max_uint))/2u)*2u;
		//LBM lbm(1u*L, 1u*L, 1u*L, 1u, 1u, 1u, 1.0f); // 1 GPU
		//LBM lbm(2u*L, 1u*L, 1u*L, 2u, 1u, 1u, 1.0f); // 2 GPUs
		//LBM lbm(2u*L, 2u*L, 1u*L, 2u, 2u, 1u, 1.0f); // 4 GPUs
		//LBM lbm(2u*L, 2u*L, 2u*L, 2u, 2u, 2u, 1.0f); // 8 GPUs
		// #########################################################################################################################################################################################
		for(uint i=0u; i<1000u; i++) {
			lbm.run(10u);
			mlups = max(mlups, to_uint((double)lbm.get_N()*1E-6/info.dt_smooth));
		}
	} // make lbm object go out of scope to free its memory
	print_info("Peak MLUPs/s = "+to_string(mlups));
#if defined(_WIN32)
	wait();
#endif // Windows
}*/
#endif // BENCHMARK

[ibonito1]

Do you know whether the geometry is watertight? Pure speculation, but (especially) CAD models from the internet often contain (probably even microscopic) holes or missing triangles. The geometry might look good, and some software programs might be able to deal with that, but for something like that might throw off the voxelization process.

For “classical” meshers for FVM CFD, geometries usually must be completely watertight and often need a lot of preparation (if you don't yourself control the source CAD file) before the meshing works.

[FCLC]

Do you know whether the geometry is watertight? Pure speculation, but (especially) CAD models from the internet often contain (probably even microscopic) holes or missing triangles. The geometry might look good, and some software programs might be able to deal with that, but for something like that might throw off the voxelization process.

For “classical” meshers for FVM CFD, geometries usually must be completely watertight and often need a lot of preparation (if you don't yourself control the source CAD file) before the meshing works.

This model has been fine for me in an OpenFOAM setup, presumably it should be fine here? There's a Zip of the STL in question in the OP (or click here: https://github.com/ProjectPhysX/FluidX3D/files/10855164/RB18_alligned.zip)

GH doesn't accept raw STLs, but in this case it's the single STL in the ZIP

[ProjectPhysX]

@FCLC the stl mesh must be watertight for the voxelization algorithm to work properly. I've built in some failsaves for non-watertight meshes, but they aren't perfect.

If you can't get the mesh watertight, to at least reduce artifacts, in this case you can override the voxelization ray direction to z-direction by inserting "direction = 2u;" between lines 224 and 225 in src/lbm.cpp.

[FCLC]

I've since adjusted/fixed the model to be considered solid (confirmed in FreeCAD and Blender 3D print submodule)

Unfortunately the behaviour remains the same:

rb-18-watertight.zip

The issue persists

[ProjectPhysX]

Can you try with lower mesh resolution?

Could be an issue with very small triangles and floating-point accuracy during Möller-Trumbore.

[FCLC]

Can you try with lower mesh resolution?

Dropping from 672 to 256 yields

At u=196:

Down to u=128

Dropping to u=64

[FCLC]

Tried 671(1 below my highest attempt before) and got a very "clean" mesh:

I suspect this is to do with small geometry not "properly" landing in place during a voxelization of even number unit geo?

Edit:

Testing further, issue is "gone" at 680, 688, 692.

Seems to be specific points cause issues

[ProjectPhysX]

@FCLC yes specific single points can be problematic, when a ray crosses exactly on the boundary between two/multiple adjacent triangles. I've already implemented it mathematically clean as ">=" for lower and "<" for upper boundary on the intersection plane, but sometimes it still counts none/both triangles, and this will trigger an artifact column where inside/outside states are inverted.

Probability of such edge/corner hits is increased when the triangle mesh is very fine compared to grid resolution, so either reducing mesh resolution or increasing grid resolution should help.

I"ll think about how to further harden the algorithm to be more robust in such cases. Thanks a lot for your help!

[FCLC]

Sounds good!

I'll close here.

Also maybe worth mentioning:
In terms of usability, it may be "useful" to have an L=automatic option to use ~95% of available ram, using the params as supplied in the grid defining size.

For some usecases, especially on "consumer" cards that don't have much VRAM, being able to maximize based on your hardware could be useful.

Another small usability idea would be to allow zooming in and out using the + and - keys (either only through the numpad or also via the keys on the main part of the keyboard ). Idea is to be slightly less cryptic than the current , and . options.

[ProjectPhysX]

@FCLC great suggestions, thanks! Will implement +/-. The L=automatic is a bit more difficult, but I can auto-rescale the box to 95% mem with a warning print if the input is too large.

[FCLC]

@FCLC great suggestions, thanks! Will implement +/-. The L=automatic is a bit more difficult, but I can auto-rescale the box to 95% mem with a warning print if the input is too large.

Other idea: showing what interactive rendering modes are active along side the simulation statistics?