From 85422fbdaf5c56ca817f0a8d8f17f2d725dd1c3a Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Wed, 6 Apr 2016 02:55:02 +0200 Subject: Finish CSG, permittivity visualisation. --- src/common.h | 7 +++ src/console.c | 14 +++++ src/physics.c | 190 ++++++++++++++++++++++++++++++++++++++++++--------------- src/physics.h | 23 ++++++- src/renderer.c | 73 +++++++++------------- src/renderer.h | 2 + src/ui.c | 10 ++- 7 files changed, 223 insertions(+), 96 deletions(-) (limited to 'src') diff --git a/src/common.h b/src/common.h index 86d98ad..b6e90d2 100644 --- a/src/common.h +++ b/src/common.h @@ -178,4 +178,11 @@ static inline void mst2_dump(mst2_t M) printf("[%f\t%f\t%f]\n", 0.0, 0.0, 1.0); } +static inline float remap(float x, float a0, float b0, float a1, float b1) +{ + return (x - a0) * (b1 - a1) / (b0 - a0) + a1; +} + +#define DEG2RAD(x) ((x) / (180.0f / M_PI)) + #endif // _COMMON_H diff --git a/src/console.c b/src/console.c index d415500..9c7cebb 100644 --- a/src/console.c +++ b/src/console.c @@ -155,6 +155,18 @@ char *cmd_load(char *arg) return va("Załadowano plik %s", arg); } +char *cmd_csg_clear(char *arg) +{ + itc_chan_push(&con_sim->ctl, PHY_CMD_CSG_CLEAR, NULL); + return "Wyzerowano przenikalność elektryczną"; +} + +char *cmd_csg_push(char *arg) +{ + itc_chan_push(&con_sim->ctl, PHY_CMD_CSG_PUSH, NULL); + return "Wykonano operację na przenikalności elektrycznej"; +} + char *cmd_get(char *arg) { if (!arg) @@ -205,6 +217,8 @@ char *cmd_zero(char *arg) static const con_cmd con_cmds[] = { + {"csg-clear", cmd_csg_clear}, + {"csg-push", cmd_csg_push}, {"get", cmd_get}, {"getf", cmd_getf}, {"getz", cmd_getz}, diff --git a/src/physics.c b/src/physics.c index fedc4c4..c7cc477 100644 --- a/src/physics.c +++ b/src/physics.c @@ -55,9 +55,48 @@ static void phy_calc_uec_E(float *uec, float dt, float fc0, float fc1, uec[3] *= 1.0f / eps; } -static float remap(float x, float a0, float b0, float a1, float b1) + +static float phy_sim_csg_subsample(phy_csg_step *stack, float *x) +{ + phy_csg_step *step; + float rv = 1.0f; + + for (step = stack; step; step = step->next) { + switch (step->shape) { + case CSG_SPHERE: + if (sq(x[0] - step->x[0]) + sq(x[1] - step->x[1]) + + sq(x[2] - step->x[2]) > sq(step->d[0])) + goto next_step; + break; + } + + rv = step->value; + + next_step: + ; + } + + return rv; +} + +static float phy_sim_csg_sample(phy_csg_step *stack, size_t *x) { - return (x - a0) * (b1 - a1) / (b0 - a0) + a1; + vec3_t xf = {x[0], x[1], x[2]}; + size_t dx, dy, dz; + float total; + + for (dz = 0; dz < 2; dz++) + for (dy = 0; dy < 2; dy++) + for (dx = 0; dx < 2; dx++) { + vec3_t xf; + + v3_set(xf, x[0] + 0.5f * dx, x[1] + 0.5f * dy, + x[2] + 0.5f * dz); + + total += phy_sim_csg_subsample(stack, xf); + } + + return total / 8.0f; } static void phy_sim_zero_aux(phy_sim *sim) @@ -68,6 +107,8 @@ static void phy_sim_zero_aux(phy_sim *sim) memset(sim->aux, 0, fi->size1 * sizeof(phy_field_aux_point)); + sim->eps_min = sim->eps_max = 1.0f; + for (x[2] = 1; x[2] < fi->dims[2] - 1; x[2]++) for (x[1] = 1; x[1] < fi->dims[1] - 1; x[1]++) for (x[0] = 1; x[0] < fi->dims[0] - 1; x[0]++) { @@ -77,23 +118,14 @@ static void phy_sim_zero_aux(phy_sim *sim) aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 + x[0] * fi->xstr1; - aux->eps = 1; - aux->mu = 1; - - if (0){ - float cx, cy, cz, cr, d; + aux->eps = phy_sim_csg_sample(sim->csg_stack, x); - cx = 25, cy = 25, cz = 65; - cr = 10; + if (aux->eps > sim->eps_max) + sim->eps_max = aux->eps; + if (aux->eps < sim->eps_min) + sim->eps_min = aux->eps; - d = sqrt(sq(x[0] - cx) + sq(x[1] - cy) + sq(x[2] - cz)); - - if (d < cr) { - aux->eps = remap(d, cr, cr * 0.8, 1, 2); - if (aux->eps > 2) - aux->eps = 2; - } - } + aux->mu = 1; v3_set(aux->int_curl_E, 0, 0, 0); v3_set(aux->int_curl_H, 0, 0, 0); @@ -122,6 +154,10 @@ static void phy_sim_zero(phy_sim *sim) size_t i; phy_field_info *fi; + sim->pml_widths[0] = con_var_get_size("pml_x"); + sim->pml_widths[1] = con_var_get_size("pml_y"); + sim->pml_widths[2] = con_var_get_size("pml_z"); + fi = &sim->field_info; for (i = 0; i < 3; i++) { @@ -149,6 +185,18 @@ int phy_sim_create(phy_sim *sim) con_var_set("depth", "50"); con_var_set("scale", "7e+9"); con_var_set("time_delta", "0.1"); + con_var_set("pml_x", "0"); + con_var_set("pml_y", "0"); + con_var_set("pml_z", "0"); + + con_var_set("src_z", "10"); + con_var_set("src_angle", "0"); + + con_var_set("csg_x", NULL); + con_var_set("csg_y", NULL); + con_var_set("csg_z", NULL); + con_var_set("csg_r", NULL); + con_var_set("csg_v", NULL); phy_sim_create_fields(sim); @@ -157,16 +205,28 @@ int phy_sim_create(phy_sim *sim) void phy_sim_destroy(phy_sim *sim) { - size_t i; - + phy_sim_csg_clear(sim); + phy_sim_destroy_fields(sim); SDL_DestroyMutex(sim->rotate_lock); - itc_chan_destroy(&sim->ctl); +} + +static void phy_sim_configure(phy_sim *sim) +{ + phy_field_info *fi = &sim->field_info; + + fi->spacing = con_var_get_float("scale") / + min3(fi->dims[0], fi->dims[1], fi->dims[2]); + sim->time_delta = con_var_get_float("time_delta"); + + con_var_set("courant", va("%f", 3 * LIGHT_SPEED * sim->time_delta / + fi->spacing)); + + sim->source_xyplane.z = con_var_get_size("src_z"); + sim->source_xyplane.Px = cos(DEG2RAD(con_var_get_float("src_angle"))); + sim->source_xyplane.Py = sin(DEG2RAD(con_var_get_float("src_angle"))); + sim->source_xyplane.enabled = true; - for (i = 0; i < 3; i++) { - free(sim->fields[i].E); - free(sim->fields[i].H); - } } int phy_sim_create_fields(phy_sim *sim) @@ -179,29 +239,9 @@ int phy_sim_create_fields(phy_sim *sim) fi->dims[0] = con_var_get_size("width"); fi->dims[1] = con_var_get_size("height"); fi->dims[2] = con_var_get_size("depth"); - fi->spacing = con_var_get_float("scale") / - min3(fi->dims[0], fi->dims[1], fi->dims[2]); - sim->time_delta = con_var_get_float("time_delta"); - con_var_set("courant", va("%f", 3 * LIGHT_SPEED * sim->time_delta / - fi->spacing)); - - #if 0 - sim->pml_widths[0] = 0; - sim->pml_widths[1] = 0; - sim->pml_widths[2] = fi->dims[2] / 5 * 1; - #else - sim->pml_widths[0] = 0; - sim->pml_widths[1] = 0; - sim->pml_widths[2] = 0; - #endif - - sim->source_xyplane.z = fi->dims[2] / 2;//sim->pml_widths[2] + 5; - sim->source_xyplane.Px = 0; - sim->source_xyplane.Py = 1; - sim->source_xyplane.lambda = 20; - sim->source_xyplane.time = 2; - sim->source_xyplane.enabled = true; + printf("phy_sim_create_fields: %zux%zux%zu\n", + fi->dims[0], fi->dims[1], fi->dims[2]); fi->xstr = 3; fi->ystr = fi->xstr * fi->dims[0]; @@ -228,6 +268,7 @@ int phy_sim_create_fields(phy_sim *sim) goto error; } + phy_sim_configure(sim); phy_sim_zero(sim); sim->running = false; @@ -470,6 +511,44 @@ void phy_sim_step(phy_sim *sim) { SDL_mutexV(sim->rotate_lock); } +void phy_sim_csg_push(phy_sim *sim) +{ + char *shape_str; + phy_csg_shape shape; + phy_csg_step *step; + + // TODO shape_str = con_var_get("csg_shape"); + + shape = CSG_SPHERE; + + step = malloc(sizeof(phy_csg_step)); + assert(step); + + step->next = sim->csg_stack; + sim->csg_stack = step; + + step->value = con_var_get_float("csg_v"); + + switch (shape) { + case CSG_SPHERE: + step->x[0] = con_var_get_float("csg_x"); + step->x[1] = con_var_get_float("csg_y"); + step->x[2] = con_var_get_float("csg_z"); + step->d[0] = con_var_get_float("csg_r"); + break; + } +} + +void phy_sim_csg_clear(phy_sim *sim) +{ + phy_csg_step *step, *next; + + for (step = sim->csg_stack; step; step = next) { + next = step->next; + free(step); + } +} + void phy_sim_debug(phy_sim *sim, size_t *coords) { phy_field_info *fi = &sim->field_info; @@ -538,10 +617,25 @@ void phy_run_command(phy_sim *sim, int number, void *data) case PHY_CMD_RESET: SDL_LockMutex(sim->rotate_lock); - phy_sim_destroy_fields(sim); - phy_sim_create_fields(sim); + if (sim->field_info.dims[0] == con_var_get_size("width") && + sim->field_info.dims[1] == con_var_get_size("height") && + sim->field_info.dims[2] == con_var_get_size("depth")) { + phy_sim_configure(sim); + phy_sim_zero(sim); + } else { + phy_sim_destroy_fields(sim); + phy_sim_create_fields(sim); + } SDL_UnlockMutex(sim->rotate_lock); break; + + case PHY_CMD_CSG_PUSH: + phy_sim_csg_push(sim); + break; + + case PHY_CMD_CSG_CLEAR: + phy_sim_csg_clear(sim); + break; } } diff --git a/src/physics.h b/src/physics.h index 60db8bf..80f3fdf 100644 --- a/src/physics.h +++ b/src/physics.h @@ -38,8 +38,6 @@ typedef struct { size_t z; float Px, Py; - - float lambda, time; } phy_source_xyplane; enum { @@ -49,7 +47,22 @@ enum { PHY_CMD_STEP, PHY_CMD_ZERO, PHY_CMD_DEBUG, - PHY_CMD_RESET + PHY_CMD_RESET, + PHY_CMD_CSG_PUSH, + PHY_CMD_CSG_CLEAR +}; + +typedef enum { + CSG_SPHERE +} phy_csg_shape; + +typedef struct phy_csg_step_s phy_csg_step; +struct phy_csg_step_s { + phy_csg_shape shape; + vec3_t x; + vec3_t d; + float value; + phy_csg_step *next; }; typedef struct { @@ -63,6 +76,9 @@ typedef struct { phy_source_xyplane source_xyplane; + phy_csg_step *csg_stack; + float eps_min, eps_max; + // UI stuff bool running; itc_chan ctl; @@ -79,6 +95,7 @@ int phy_sim_create_fields(phy_sim *sim); void phy_sim_destroy_fields(phy_sim *sim); void phy_sim_compute_const_fields(phy_sim *sim); void phy_sim_step(phy_sim *sim); +void phy_sim_csg_clear(phy_sim *sim); int phy_thread(phy_sim *sim); diff --git a/src/renderer.c b/src/renderer.c index d68287c..6d159e7 100644 --- a/src/renderer.c +++ b/src/renderer.c @@ -373,7 +373,7 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, { phy_field_info *fi = &sim->field_info; float *field; - size_t width, height, x, y, z; + size_t width, height, x[3], dx, dy, d0; uint8_t *pixels; int pitch; float scale; @@ -387,12 +387,21 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, if (flags & XSECTION_XY) { width = fi->dims[0]; height = fi->dims[1]; + d0 = 2; + dx = 0; + dy = 1; } else if (flags & XSECTION_XZ) { width = fi->dims[0]; height = fi->dims[2]; + d0 = 1; + dx = 0; + dy = 2; } else { width = fi->dims[1]; height = fi->dims[2]; + d0 = 0; + dx = 1; + dy = 2; } if (flags & XSECTION_E) { @@ -419,61 +428,39 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, SDL_LockTexture(xsection->texture, NULL, (void**)&pixels, &pitch); - if (flags & XSECTION_XY) { - z = frac * fi->dims[2]; - if (z >= fi->dims[2]) - z = fi->dims[2] - 1; + x[d0] = frac * fi->dims[d0]; + if (x[d0] >= fi->dims[d0]) + x[d0] = fi->dims[d0] - 1; - for (y = 0; y < fi->dims[1]; y++) - for (x = 0; x < fi->dims[0]; x++) { + if (flags & XSECTION_EPS) + for (x[dy] = 0; x[dy] < fi->dims[dy]; x[dy]++) + for (x[dx] = 0; x[dx] < fi->dims[dx]; x[dx]++) { uint8_t *pixel; - float *point; + float gray; + phy_field_aux_point *aux; - pixel = pixels + (y * pitch + x * 4); - point = field + z * fi->zstr + y * fi->ystr + - x * fi->xstr; + pixel = pixels + (x[dy] * pitch + x[dx] * 4); + aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 + + x[0] * fi->xstr1; - pixel[0] = r_float_to_u8(point[0] * scale); - pixel[1] = r_float_to_u8(point[1] * scale); - pixel[2] = r_float_to_u8(point[2] * scale); + gray = remap(aux->eps, sim->eps_min, sim->eps_max, + 0.0f, 1.0f); + pixel[0] = pixel[1] = pixel[2] = r_float_to_u8(gray); } - } else if (flags & XSECTION_XZ) { - y = frac * fi->dims[1]; - if (y >= fi->dims[1]) - y = fi->dims[1] - 1; - - for (z = 0; z < fi->dims[2]; z++) - for (x = 0; x < fi->dims[0]; x++) { - uint8_t *pixel; - float *point; - - pixel = pixels + (z * pitch + x * 4); - point = field + z * fi->zstr + y * fi->ystr + - x * fi->xstr; - - pixel[0] = r_float_to_u8(point[0] * scale); - pixel[1] = r_float_to_u8(point[1] * scale); - pixel[2] = r_float_to_u8(point[2] * scale); - } - } else { - x = frac * fi->dims[0]; - if (x >= fi->dims[0]) - x = fi->dims[0] - 1; - - for (z = 0; z < fi->dims[2]; z++) - for (y = 0; y < fi->dims[1]; y++) { + else + for (x[dy] = 0; x[dy] < fi->dims[dy]; x[dy]++) + for (x[dx] = 0; x[dx] < fi->dims[dx]; x[dx]++) { uint8_t *pixel; float *point; - pixel = pixels + (z * pitch + y * 4); - point = field + z * fi->zstr + y * fi->ystr + - x * fi->xstr; + pixel = pixels + (x[dy] * pitch + x[dx] * 4); + point = field + x[2] * fi->zstr + x[1] * fi->ystr + + x[0] * fi->xstr; pixel[0] = r_float_to_u8(point[0] * scale); pixel[1] = r_float_to_u8(point[1] * scale); pixel[2] = r_float_to_u8(point[2] * scale); } - } SDL_UnlockTexture(xsection->texture); diff --git a/src/renderer.h b/src/renderer.h index fadcd68..3ac10b8 100644 --- a/src/renderer.h +++ b/src/renderer.h @@ -47,6 +47,8 @@ void r_draw_text(r_window *rw, float x, float y, float h, char *text, #define XSECTION_PLANES (XSECTION_XY|XSECTION_XZ|XSECTION_YZ) #define XSECTION_E 0x0008 #define XSECTION_H 0x0010 +#define XSECTION_EPS 0x0020 +#define XSECTION_FIELDS (XSECTION_E|XSECTION_H|XSECTION_EPS) typedef struct { diff --git a/src/ui.c b/src/ui.c index 5b58be7..868b6a1 100644 --- a/src/ui.c +++ b/src/ui.c @@ -338,16 +338,22 @@ void ui_event_window(SDL_Event *event, ui_window *uiw) case SDLK_4: ui_infof(uiw, "Przekrój: pola elektrycznego"); - sv->xsection_flags &= ~XSECTION_H; + sv->xsection_flags &= ~XSECTION_FIELDS; sv->xsection_flags |= XSECTION_E; break; case SDLK_5: ui_infof(uiw, "Przekrój: pola magnetycznego"); - sv->xsection_flags &= ~XSECTION_E; + sv->xsection_flags &= ~XSECTION_FIELDS; sv->xsection_flags |= XSECTION_H; break; + case SDLK_6: + ui_infof(uiw, "Przekrój: pola przenikalności"); + sv->xsection_flags &= ~XSECTION_FIELDS; + sv->xsection_flags |= XSECTION_EPS; + break; + case SDLK_z: if (sv->select_valid) { -- cgit