diff --git a/src/kernel.cpp b/src/kernel.cpp
index bca41ef..35dcb1a 100644
--- a/src/kernel.cpp
+++ b/src/kernel.cpp
@@ -2383,6 +2383,94 @@ string opencl_c_container() { return R( // ########################## begin of O
 	}
 } // voxelize_mesh()
 
+)+R(kernel void voxelize_sdf)+"("+R(global uchar* flags, const uchar flag, const global float* sdf_data, const global float* sdf_params)+") {"+R( // voxelize_sdf()
+	const uxx n = get_global_id(0);
+	if(n>=(uxx)def_N) return;
+
+	// Unpack parameters
+	const uint sdf_Nx = as_uint(sdf_params[0]);
+	const uint sdf_Ny = as_uint(sdf_params[1]);
+	const uint sdf_Nz = as_uint(sdf_params[2]);
+	const float sdf_cell_size = sdf_params[3];
+	const float3 sdf_origin = (float3)(sdf_params[4], sdf_params[5], sdf_params[6]);
+	const float3 center_lbm = (float3)(sdf_params[7], sdf_params[8], sdf_params[9]);
+	const float world_to_lbm = sdf_params[10];
+	// Inverse rotation matrix (row-major)
+	const float3 rot_row0 = (float3)(sdf_params[11], sdf_params[12], sdf_params[13]);
+	const float3 rot_row1 = (float3)(sdf_params[14], sdf_params[15], sdf_params[16]);
+	const float3 rot_row2 = (float3)(sdf_params[17], sdf_params[18], sdf_params[19]);
+
+	// Get LBM coordinates
+	const uint3 xyz = coordinates(n);
+	const float3 lbm_pos = (float3)(xyz.x + 0.5f, xyz.y + 0.5f, xyz.z + 0.5f);
+
+	// Transform LBM -> World -> SDF coordinates
+	// 1. Compute offset from LBM center
+	const float3 offset_lbm = (lbm_pos - center_lbm) / world_to_lbm;
+	// 2. Apply inverse rotation to get offset in SDF world space
+	const float3 offset_world = (float3)(dot(rot_row0, offset_lbm), dot(rot_row1, offset_lbm), dot(rot_row2, offset_lbm));
+	// 3. Add SDF world center to get world position
+	const float3 sdf_world_size = (float3)(sdf_Nx, sdf_Ny, sdf_Nz) * sdf_cell_size;
+	const float3 sdf_world_center = sdf_origin + 0.5f * sdf_world_size;
+	const float3 world_pos = offset_world + sdf_world_center;
+	// 4. Convert to SDF grid coordinates
+	const float3 sdf_coord = (world_pos - sdf_origin) / sdf_cell_size;
+
+	// Check bounds
+	if(sdf_coord.x < 0.0f || sdf_coord.x >= (float)(sdf_Nx-1u) ||
+	   sdf_coord.y < 0.0f || sdf_coord.y >= (float)(sdf_Ny-1u) ||
+	   sdf_coord.z < 0.0f || sdf_coord.z >= (float)(sdf_Nz-1u)) {
+		return;  // Outside SDF bounds
+	}
+
+	// Trilinear interpolation
+	const uint x0 = (uint)sdf_coord.x;
+	const uint y0 = (uint)sdf_coord.y;
+	const uint z0 = (uint)sdf_coord.z;
+	const uint x1 = min(x0 + 1u, sdf_Nx - 1u);
+	const uint y1 = min(y0 + 1u, sdf_Ny - 1u);
+	const uint z1 = min(z0 + 1u, sdf_Nz - 1u);
+
+	const float fx = sdf_coord.x - (float)x0;
+	const float fy = sdf_coord.y - (float)y0;
+	const float fz = sdf_coord.z - (float)z0;
+
+	// Sample SDF at 8 corners (SDF uses x + Nx*(y + Ny*z) indexing)
+	const ulong sdf_index_000 = (ulong)x0 + (ulong)sdf_Nx * ((ulong)y0 + (ulong)sdf_Ny * (ulong)z0);
+	const ulong sdf_index_001 = (ulong)x0 + (ulong)sdf_Nx * ((ulong)y0 + (ulong)sdf_Ny * (ulong)z1);
+	const ulong sdf_index_010 = (ulong)x0 + (ulong)sdf_Nx * ((ulong)y1 + (ulong)sdf_Ny * (ulong)z0);
+	const ulong sdf_index_011 = (ulong)x0 + (ulong)sdf_Nx * ((ulong)y1 + (ulong)sdf_Ny * (ulong)z1);
+	const ulong sdf_index_100 = (ulong)x1 + (ulong)sdf_Nx * ((ulong)y0 + (ulong)sdf_Ny * (ulong)z0);
+	const ulong sdf_index_101 = (ulong)x1 + (ulong)sdf_Nx * ((ulong)y0 + (ulong)sdf_Ny * (ulong)z1);
+	const ulong sdf_index_110 = (ulong)x1 + (ulong)sdf_Nx * ((ulong)y1 + (ulong)sdf_Ny * (ulong)z0);
+	const ulong sdf_index_111 = (ulong)x1 + (ulong)sdf_Nx * ((ulong)y1 + (ulong)sdf_Ny * (ulong)z1);
+
+	const float c000 = sdf_data[sdf_index_000];
+	const float c001 = sdf_data[sdf_index_001];
+	const float c010 = sdf_data[sdf_index_010];
+	const float c011 = sdf_data[sdf_index_011];
+	const float c100 = sdf_data[sdf_index_100];
+	const float c101 = sdf_data[sdf_index_101];
+	const float c110 = sdf_data[sdf_index_110];
+	const float c111 = sdf_data[sdf_index_111];
+
+	// Trilinear interpolation
+	const float c00 = c000 * (1.0f - fx) + c100 * fx;
+	const float c01 = c001 * (1.0f - fx) + c101 * fx;
+	const float c10 = c010 * (1.0f - fx) + c110 * fx;
+	const float c11 = c011 * (1.0f - fx) + c111 * fx;
+
+	const float c0 = c00 * (1.0f - fy) + c10 * fy;
+	const float c1 = c01 * (1.0f - fy) + c11 * fy;
+
+	const float sdf_value = c0 * (1.0f - fz) + c1 * fz;
+
+	// Mark as solid if inside (negative SDF)
+	if(sdf_value <= 0.0f) {
+		flags[n] = (flags[n] & ~TYPE_BO) | flag;
+	}
+} // voxelize_sdf()
+
 )+R(kernel void unvoxelize_mesh(global uchar* flags, const uchar flag, float x0, float y0, float z0, float x1, float y1, float z1) { // remove voxelized triangle mesh
 	const uxx n = get_global_id(0);
 	const float3 p = position(coordinates(n))+(float3)(0.5f*(float)((int)def_Nx+2*def_Ox)-0.5f, 0.5f*(float)((int)def_Ny+2*def_Oy)-0.5f, 0.5f*(float)((int)def_Nz+2*def_Oz)-0.5f);
diff --git a/src/lbm.cpp b/src/lbm.cpp
index 64fca45..5e32e3c 100644
--- a/src/lbm.cpp
+++ b/src/lbm.cpp
@@ -325,6 +325,43 @@ void LBM_Domain::voxelize_mesh_on_device(const Mesh* mesh, const uchar flag, con
 	bounding_box_and_velocity.write_to_device();
 	kernel_voxelize_mesh.run();
 }
+void LBM_Domain::voxelize_sdf_on_device(const SDF* sdf, const float3& center, const float3x3& rotation, const float scale, const uchar flag) { // voxelize SDF on GPU
+	// Upload SDF data to GPU
+	const ulong sdf_size = (ulong)sdf->Nx * (ulong)sdf->Ny * (ulong)sdf->Nz;
+	Memory<float> sdf_data(device, sdf_size, 1u, sdf->data);
+
+	// Compute inverse rotation matrix (to transform LBM coords back to SDF space)
+	// For orthogonal matrices, inverse = transpose
+	const float3x3 rotation_inv = transpose(rotation);
+
+	// Pack parameters into array (similar to voxelize_mesh approach)
+	Memory<float> sdf_params(device, 20u);
+	const float3 sdf_world_size = sdf->get_world_size();
+	const float sdf_max_world_dim = fmax(fmax(sdf_world_size.x, sdf_world_size.y), sdf_world_size.z);
+	const float world_to_lbm = scale / sdf_max_world_dim;
+
+	sdf_params[0] = as_float(sdf->Nx);      // pack uint as float
+	sdf_params[1] = as_float(sdf->Ny);
+	sdf_params[2] = as_float(sdf->Nz);
+	sdf_params[3] = sdf->cell_size;
+	sdf_params[4] = sdf->origin.x;          // SDF origin
+	sdf_params[5] = sdf->origin.y;
+	sdf_params[6] = sdf->origin.z;
+	sdf_params[7] = center.x;               // LBM center
+	sdf_params[8] = center.y;
+	sdf_params[9] = center.z;
+	sdf_params[10] = world_to_lbm;
+	// Inverse rotation matrix (row-major: row0, row1, row2)
+	sdf_params[11] = rotation_inv.xx; sdf_params[12] = rotation_inv.xy; sdf_params[13] = rotation_inv.xz;
+	sdf_params[14] = rotation_inv.yx; sdf_params[15] = rotation_inv.yy; sdf_params[16] = rotation_inv.yz;
+	sdf_params[17] = rotation_inv.zx; sdf_params[18] = rotation_inv.zy; sdf_params[19] = rotation_inv.zz;
+
+	// Create kernel and run
+	Kernel kernel_voxelize_sdf(device, get_N(), "voxelize_sdf", flags, flag, sdf_data, sdf_params);
+	sdf_data.write_to_device();
+	sdf_params.write_to_device();
+	kernel_voxelize_sdf.run();
+}
 void LBM_Domain::enqueue_unvoxelize_mesh_on_device(const Mesh* mesh, const uchar flag) { // remove voxelized triangle mesh from LBM grid
 	const float x0=mesh->pmin.x, y0=mesh->pmin.y, z0=mesh->pmin.z, x1=mesh->pmax.x, y1=mesh->pmax.y, z1=mesh->pmax.z; // remove all flags in bounding box of mesh
 	Kernel kernel_unvoxelize_mesh(device, get_N(), "unvoxelize_mesh", flags, flag, x0, y0, z0, x1, y1, z1);
@@ -1111,6 +1148,42 @@ void LBM::voxelize_stl(const string& path, const float size, const uchar flag) {
 	voxelize_stl(path, center(), float3x3(1.0f), size, flag);
 }
 
+void LBM::voxelize_sdf(const string& path, const float3& center, const float3x3& rotation, const float scale, const uchar flag) { // voxelize SDF
+	const SDF* sdf = read_sdf(path);
+	if(sdf == nullptr) return; // Error already printed by read_sdf
+
+	flags.write_to_device();
+
+	print_info("Voxelizing SDF on GPU...");
+	print_info("LBM grid: " + to_string(get_Nx()) + " x " + to_string(get_Ny()) + " x " + to_string(get_Nz()));
+
+	// Voxelize on GPU using kernel
+	for(uint d = 0u; d < get_D(); d++) {
+		lbm_domain[d]->voxelize_sdf_on_device(sdf, center, rotation, scale, flag);
+	}
+	flags.read_from_device();
+
+	delete sdf;
+
+	// Count total solid voxels
+	uint total_solid = 0;
+	for(ulong n = 0; n < get_N(); n++) {
+		if(flags[n] == flag) total_solid++;
+	}
+
+	flags.write_to_device(); // Push the modified flags back to GPU
+	print_info("SDF voxelization complete. Total solid voxels: " + to_string(total_solid));
+}
+void LBM::voxelize_sdf(const string& path, const float3x3& rotation, const float scale, const uchar flag) { // read and voxelize binary SDF file (place in box center)
+	voxelize_sdf(path, center(), rotation, scale, flag);
+}
+void LBM::voxelize_sdf(const string& path, const float3& center, const float scale, const uchar flag) { // read and voxelize binary SDF file (no rotation)
+	voxelize_sdf(path, center, float3x3(1.0f), scale, flag);
+}
+void LBM::voxelize_sdf(const string& path, const float scale, const uchar flag) { // read and voxelize binary SDF file (place in box center, no rotation)
+	voxelize_sdf(path, center(), float3x3(1.0f), scale, flag);
+}
+
 #ifdef GRAPHICS
 int* LBM::Graphics::draw_frame() {
 #ifndef UPDATE_FIELDS
diff --git a/src/lbm.hpp b/src/lbm.hpp
index a4c9efa..9ff87f2 100644
--- a/src/lbm.hpp
+++ b/src/lbm.hpp
@@ -144,6 +144,7 @@ class LBM_Domain {
 	void set_f(const float fx, const float fy, const float fz) { set_fx(fx); set_fy(fy); set_fz(fz); } // set global froce per volume
 
 	void voxelize_mesh_on_device(const Mesh* mesh, const uchar flag=TYPE_S, const float3& rotation_center=float3(0.0f), const float3& linear_velocity=float3(0.0f), const float3& rotational_velocity=float3(0.0f)); // voxelize mesh
+	void voxelize_sdf_on_device(const SDF* sdf, const float3& center, const float3x3& rotation, const float scale, const uchar flag); // voxelize SDF on GPU
 	void enqueue_unvoxelize_mesh_on_device(const Mesh* mesh, const uchar flag=TYPE_S); // remove voxelized triangle mesh from LBM grid
 
 #ifdef GRAPHICS
@@ -542,6 +543,10 @@ class LBM {
 	void voxelize_stl(const string& path, const float3x3& rotation, const float size=0.0f, const uchar flag=TYPE_S); // read and voxelize binary .stl file (place in box center)
 	void voxelize_stl(const string& path, const float3& center, const float size=0.0f, const uchar flag=TYPE_S); // read and voxelize binary .stl file (no rotation)
 	void voxelize_stl(const string& path, const float size=0.0f, const uchar flag=TYPE_S); // read and voxelize binary .stl file (place in box center, no rotation)
+	void voxelize_sdf(const string& path, const float3& center, const float3x3& rotation, const float scale=1.0f, const uchar flag=TYPE_S); // read and voxelize binary SDF file
+	void voxelize_sdf(const string& path, const float3x3& rotation, const float scale=1.0f, const uchar flag=TYPE_S); // read and voxelize binary SDF file (place in box center)
+	void voxelize_sdf(const string& path, const float3& center, const float scale=1.0f, const uchar flag=TYPE_S); // read and voxelize binary SDF file (no rotation)
+	void voxelize_sdf(const string& path, const float scale=1.0f, const uchar flag=TYPE_S); // read and voxelize binary SDF file (place in box center, no rotation)
 
 #ifdef GRAPHICS
 	class Graphics {
diff --git a/src/utilities.hpp b/src/utilities.hpp
index e579dbc..0ac0a43 100644
--- a/src/utilities.hpp
+++ b/src/utilities.hpp
@@ -1211,6 +1211,9 @@ inline float3 operator*(const float3& v, const float3x3& m) { // multiply vector
 inline float3 operator*(const float3x3& m, const float3& v) { // multiply matrix with vector
 		return float3(m.xx*v.x+m.xy*v.y+m.xz*v.z, m.yx*v.x+m.yy*v.y+m.yz*v.z, m.zx*v.x+m.zy*v.y+m.zz*v.z);
 }
+inline float3x3 transpose(const float3x3& m) { // transpose matrix
+	return float3x3(m.xx, m.yx, m.zx, m.xy, m.yy, m.zy, m.xz, m.yz, m.zz);
+}
 inline float3::float3(const float3x3& m) { // extract diagonal of matrix
 	this->x = m.xx;
 	this->y = m.yy;
@@ -4527,6 +4530,151 @@ struct Mesh { // triangle mesh
 		return fmin(fmin(box_size.x/(pmax.x-pmin.x), box_size.y/(pmax.y-pmin.y)), box_size.z/(pmax.z-pmin.z));
 	}
 };
+
+// ########################### SDF (Signed Distance Field) Support ###########################
+
+struct SDF { // 3D signed distance field
+	uint Nx = 0u, Ny = 0u, Nz = 0u; // Grid dimensions
+	float cell_size = 1.0f; // Size of each cell in world units
+	float3 origin; // World-space origin (lower corner) of the SDF grid
+	float* data = nullptr; // SDF values: negative=inside, positive=outside, zero=surface
+
+	inline SDF(const uint Nx, const uint Ny, const uint Nz, const float cell_size, const float3& origin) {
+		this->Nx = Nx;
+		this->Ny = Ny;
+		this->Nz = Nz;
+		this->cell_size = cell_size;
+		this->origin = origin;
+		this->data = new float[Nx*Ny*Nz];
+	}
+
+	inline ~SDF() {
+		if(data) delete[] data;
+	}
+
+	inline uint index(const uint x, const uint y, const uint z) const {
+		return x + Nx * (y + Ny * z);
+	}
+
+	inline float get(const uint x, const uint y, const uint z) const {
+		return data[index(x, y, z)];
+	}
+
+	inline void set(const uint x, const uint y, const uint z, const float value) {
+		data[index(x, y, z)] = value;
+	}
+
+	inline float3 get_world_size() const {
+		return float3((float)Nx * cell_size, (float)Ny * cell_size, (float)Nz * cell_size);
+	}
+
+	inline float3 get_bounding_box_size() const {
+		return get_world_size();
+	}
+
+	inline float3 get_bounding_box_center() const {
+		return origin + 0.5f * get_world_size();
+	}
+};
+
+/**
+ * @brief Read binary SDF file
+ * @param path Path to binary SDF file
+ * @return Pointer to SDF structure
+ *
+ * File format (binary, little-endian):
+ * - int32: Nx (grid dimension X)
+ * - int32: Ny (grid dimension Y)
+ * - int32: Nz (grid dimension Z)
+ * - float32: bounds_min_x (bounding box minimum X)
+ * - float32: bounds_min_y (bounding box minimum Y)
+ * - float32: bounds_min_z (bounding box minimum Z)
+ * - float32: bounds_max_x (bounding box maximum X)
+ * - float32: bounds_max_y (bounding box maximum Y)
+ * - float32: bounds_max_z (bounding box maximum Z)
+ * - float32[Nx*Ny*Nz]: SDF values (negative=inside, positive=outside, zero=surface)
+ *   Layout: data[x + Nx*(y + Ny*z)]
+ */
+inline SDF* read_sdf(const string& path) {
+	const string filename = create_file_extension(path, ".sdf");
+	std::ifstream file(filename, std::ios::in | std::ios::binary);
+
+	if(!file.is_open()) {
+		print_error("Error: SDF file \"" + filename + "\" could not be opened!");
+		return nullptr;
+	}
+
+	// Read header: grid dimensions + bounding box
+	int Nx, Ny, Nz;
+	float bounds_min_x, bounds_min_y, bounds_min_z;
+	float bounds_max_x, bounds_max_y, bounds_max_z;
+
+	file.read(reinterpret_cast<char*>(&Nx), sizeof(int));
+	file.read(reinterpret_cast<char*>(&Ny), sizeof(int));
+	file.read(reinterpret_cast<char*>(&Nz), sizeof(int));
+	file.read(reinterpret_cast<char*>(&bounds_min_x), sizeof(float));
+	file.read(reinterpret_cast<char*>(&bounds_min_y), sizeof(float));
+	file.read(reinterpret_cast<char*>(&bounds_min_z), sizeof(float));
+	file.read(reinterpret_cast<char*>(&bounds_max_x), sizeof(float));
+	file.read(reinterpret_cast<char*>(&bounds_max_y), sizeof(float));
+	file.read(reinterpret_cast<char*>(&bounds_max_z), sizeof(float));
+
+	if(file.fail()) {
+		print_error("Error: Failed to read SDF header from \"" + filename + "\"!");
+		file.close();
+		return nullptr;
+	}
+
+	// Calculate cell size and origin from bounding box
+	const float3 bounds_min(bounds_min_x, bounds_min_y, bounds_min_z);
+	const float3 bounds_max(bounds_max_x, bounds_max_y, bounds_max_z);
+	const float3 bounds_size = bounds_max - bounds_min;
+
+	// Cell size = bounding box size / (grid dimensions - 1) for grid corners
+	// Or = bounding box size / grid dimensions for grid centers
+	// Using grid centers interpretation:
+	const float cell_size_x = bounds_size.x / (float)Nx;
+	const float cell_size_y = bounds_size.y / (float)Ny;
+	const float cell_size_z = bounds_size.z / (float)Nz;
+	const float cell_size = (cell_size_x + cell_size_y + cell_size_z) / 3.0f; // Average
+
+	const float3 origin = bounds_min;
+	SDF* sdf = new SDF((uint)Nx, (uint)Ny, (uint)Nz, cell_size, origin);
+
+	// Read SDF data
+	const ulong data_size = (ulong)Nx * (ulong)Ny * (ulong)Nz;
+	float* temp_data = new float[data_size];
+	file.read(reinterpret_cast<char*>(temp_data), data_size * sizeof(float));
+
+	if(file.fail()) {
+		print_error("Error: Failed to read SDF data from \"" + filename + "\"!");
+		delete[] temp_data;
+		delete sdf;
+		file.close();
+		return nullptr;
+	}
+
+	file.close();
+
+	// Transpose data from file layout (k varies fastest) to our layout (x varies fastest)
+	// SDFGen writes: for(i) for(j) for(k) write(value) -> index: i*Ny*Nz + j*Nz + k
+	// We use: index(x,y,z) = x + Nx*(y + Ny*z) -> x varies fastest
+	for(uint i = 0; i < (uint)Nx; i++) {
+		for(uint j = 0; j < (uint)Ny; j++) {
+			for(uint k = 0; k < (uint)Nz; k++) {
+				const ulong file_index = (ulong)i * (ulong)Ny * (ulong)Nz + (ulong)j * (ulong)Nz + (ulong)k;
+				sdf->set(i, j, k, temp_data[file_index]);
+			}
+		}
+	}
+	delete[] temp_data;
+
+	print_info("Loaded SDF: " + to_string(Nx) + "x" + to_string(Ny) + "x" + to_string(Nz) +
+	           " cells, cell_size=" + to_string(cell_size) + "m");
+
+	return sdf;
+}
+
 inline Mesh* read_stl_raw(const string& path, const bool reposition, const float3& box_size, const float3& center, const float3x3& rotation, const float size) { // read binary .stl file
 	const string filename = create_file_extension(path, ".stl");
 	std::ifstream file(filename, std::ios::in|std::ios::binary);