From 3ef209bad0cd7c1ae83a65702633e43f5a91c97f Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Sun, 3 Apr 2016 17:01:06 +0200 Subject: Compute fake conductivity (for PML). --- src/physics.c | 140 ++++++++++++++++++++++++++++++++++++--------------------- src/physics.h | 5 ++- src/renderer.c | 42 ++++++++--------- src/ui.c | 12 ++--- 4 files changed, 120 insertions(+), 79 deletions(-) (limited to 'src') diff --git a/src/physics.c b/src/physics.c index a56ad97..bbe2499 100644 --- a/src/physics.c +++ b/src/physics.c @@ -1,4 +1,4 @@ -#include "common.h" +#include "common.h" #include "itc.h" #include "physics.h" @@ -16,24 +16,46 @@ void phy_sim_destroy(phy_sim *sim) } } +static float phy_calc_fake_cond(float *fake_cond, size_t *x, size_t *widths, + size_t *dims, float dt) +{ + size_t i; + + for (i = 0; i < 3; i++) { + float f; + + if (x[i] <= widths[i]) + f = 1.0f - (float)(x[i] - 1) / widths[i]; + else if (x[i] >= dims[i] - widths[i] - 1) + f = (float)(x[i] + widths[i] - dims[i] + 2) / widths[i]; + else + continue; + + fake_cond[i] = EPS_ZERO / (2 * dt) * pow(f, 3); + } +} + void phy_sim_reset_aux(phy_sim *sim) { - size_t x, y, z; + size_t x[3], i; phy_field_info *fi = &sim->field_info; - for (z = 0; z < fi->depth; z++) - for (y = 0; y < fi->height; y++) - for (x = 0; x < fi->width; x++) { - phy_field_aux_point *point; + 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]++) { + phy_field_aux_point *aux; + + aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 + + x[0] * fi->xstr1; - point = sim->aux + z * fi->zstr1 + y * fi->ystr1 + - x * fi->xstr1; + aux->eps = 1.0f / EPS_ZERO; + aux->mu = -1.0f / MU_ZERO; - point->eps = 1.0f / EPS_ZERO; - point->mu = -1.0f / MU_ZERO; + v3_set(aux->int_curl_E, 0, 0, 0); + v3_set(aux->int_curl_H, 0, 0, 0); - v3_set(point->int_curl_E, 0, 0, 0); - v3_set(point->int_curl_H, 0, 0, 0); + phy_calc_fake_cond(aux->fake_cond, x, sim->pml_widths, + fi->dims, sim->time_delta); } } @@ -65,24 +87,28 @@ int phy_sim_create(phy_sim *sim) fi = &sim->field_info; - fi->width = 100; - fi->height = 100; - fi->depth = 100; - fi->spacing = 10e+9f / max3(fi->width, fi->height, fi->depth); + fi->dims[0] = 50; + fi->dims[1] = 50; + fi->dims[2] = 50; + fi->spacing = 10e+9f / max3(fi->dims[0], fi->dims[1], fi->dims[2]); sim->time_delta = 0.1; printf("phy_sim_create: Courant number is %f\n", 3 * LIGHT_SPEED * sim->time_delta / fi->spacing); + sim->pml_widths[0] = fi->dims[0] / 5; + sim->pml_widths[1] = fi->dims[1] / 5; + sim->pml_widths[2] = fi->dims[2] / 5; + fi->xstr = 3; - fi->ystr = fi->xstr * fi->width; - fi->zstr = fi->ystr * fi->height; - fi->size = fi->zstr * fi->depth; + fi->ystr = fi->xstr * fi->dims[0]; + fi->zstr = fi->ystr * fi->dims[1]; + fi->size = fi->zstr * fi->dims[2]; fi->xstr1 = 1; - fi->ystr1 = fi->xstr1 * fi->width; - fi->zstr1 = fi->ystr1 * fi->height; - fi->size1 = fi->zstr1 * fi->depth; + fi->ystr1 = fi->xstr1 * fi->dims[0]; + fi->zstr1 = fi->ystr1 * fi->dims[1]; + fi->size1 = fi->zstr1 * fi->dims[2]; for (i = 0; i < 3; i++) { sim->fields[i].E = malloc(fi->size * sizeof(float)); @@ -118,9 +144,9 @@ void phy_sim_update_aux_H(phy_sim *sim) phy_field_info *fi = &sim->field_info; size_t x, y, z; - for (z = 1; z < fi->depth - 1; z++) - for (y = 1; y < fi->height - 1; y++) - for (x = 1; x < fi->width - 1; x++) { + for (z = 1; z < fi->dims[2] - 1; z++) + for (y = 1; y < fi->dims[1] - 1; y++) + for (x = 1; x < fi->dims[0] - 1; x++) { size_t offs, offs_px, offs_py, offs_pz; phy_field_aux_point *aux; float *E = sim->fields[1].E, *H = sim->fields[1].H; @@ -150,9 +176,9 @@ void phy_sim_update_aux_E(phy_sim *sim) phy_field_info *fi = &sim->field_info; size_t x, y, z; - for (z = 1; z < fi->depth - 1; z++) - for (y = 1; y < fi->height - 1; y++) - for (x = 1; x < fi->width - 1; x++) { + for (z = 1; z < fi->dims[2] - 1; z++) + for (y = 1; y < fi->dims[1] - 1; y++) + for (x = 1; x < fi->dims[0] - 1; x++) { size_t offs, offs_nx, offs_ny, offs_nz; phy_field_aux_point *aux; float *E = sim->fields[1].E, *H = sim->fields[1].H; @@ -187,19 +213,19 @@ void phy_sim_update_em_H(phy_sim *sim) Hp = sim->fields[0].H; dt = sim->time_delta; - for (z = 1; z < fi->depth - 1; z++) - for (y = 1; y < fi->height - 1; y++) - for (x = 1; x < fi->width - 1; x++) { + for (z = 1; z < fi->dims[2] - 1; z++) + for (y = 1; y < fi->dims[1] - 1; y++) + for (x = 1; x < fi->dims[0] - 1; x++) { size_t offs; phy_field_aux_point *aux; offs = x * fi->xstr + y * fi->ystr + z * fi->zstr; aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + - z * fi->zstr1; + z * fi->zstr1; - Hn[offs+0] = Hp[offs+0] + dt * aux->eps * aux->curl_E[0]; - Hn[offs+1] = Hp[offs+1] + dt * aux->eps * aux->curl_E[1]; - Hn[offs+2] = Hp[offs+2] + dt * aux->eps * aux->curl_E[2]; + Hn[offs+0] = Hp[offs+0] + dt * aux->mu * aux->curl_E[0]; + Hn[offs+1] = Hp[offs+1] + dt * aux->mu * aux->curl_E[1]; + Hn[offs+2] = Hp[offs+2] + dt * aux->mu * aux->curl_E[2]; } } @@ -213,9 +239,9 @@ void phy_sim_update_em_E(phy_sim *sim) Ep = sim->fields[0].E; dt = sim->time_delta; - for (z = 1; z < fi->depth - 1; z++) - for (y = 1; y < fi->height - 1; y++) - for (x = 1; x < fi->width - 1; x++) { + for (z = 1; z < fi->dims[2] - 1; z++) + for (y = 1; y < fi->dims[1] - 1; y++) + for (x = 1; x < fi->dims[0] - 1; x++) { size_t offs; phy_field_aux_point *aux; @@ -223,9 +249,9 @@ void phy_sim_update_em_E(phy_sim *sim) aux = sim->aux + x * fi->xstr1 + y * fi->ystr1 + z * fi->zstr1; - En[offs+0] = Ep[offs+0] + dt * aux->mu * aux->curl_H[0]; - En[offs+1] = Ep[offs+1] + dt * aux->mu * aux->curl_H[1]; - En[offs+2] = Ep[offs+2] + dt * aux->mu * aux->curl_H[2]; + En[offs+0] = Ep[offs+0] + dt * aux->eps * aux->curl_H[0]; + En[offs+1] = Ep[offs+1] + dt * aux->eps * aux->curl_H[1]; + En[offs+2] = Ep[offs+2] + dt * aux->eps * aux->curl_H[2]; } } @@ -235,34 +261,44 @@ void phy_sim_add_sources_E(phy_sim *sim) float *E; size_t x, y, z; - x = 0; - for (y = 1; y < fi->height - 1; y++) - for (z = 1; z < fi->depth - 1; z++) { + + x = fi->dims[0]/2; + y = fi->dims[1]/2; + z = fi->dims[2]/2; + E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + + x * fi->xstr; + + E[0] += exp(-40.0f * pow(sim->time - 3, 2)) * 50; + +/* + x = 10; + for (y = 1; y < fi->dims[1] - 1; y++) + for (z = 1; z < fi->dims[2] - 1; z++) { E = sim->fields[2].E + z * fi->zstr + y * fi->ystr + x * fi->xstr; - E[0] += sin(sim->time) * 100; + E[0] += sin(sim->time * 3) * 1; E[1] += 0; E[2] += 0; - } + }*/ } void phy_sim_add_sources_H(phy_sim *sim) { - phy_field_info *fi = &sim->field_info; + /*phy_field_info *fi = &sim->field_info; float *H; size_t x, y, z; x = 0; - for (y = 1; y < fi->height - 1; y++) - for (z = 1; z < fi->depth - 1; z++) { + for (y = 1; y < fi->dims[1] - 1; y++) + for (z = 1; z < fi->dims[2] - 1; z++) { H = sim->fields[2].H + z * fi->zstr + y * fi->ystr + x * fi->xstr; H[0] += 0; H[1] += 0; H[2] += sin(sim->time + sim->time_delta / 2) / IMPEDANCE_ZERO * 100; - } + }*/ } @@ -310,12 +346,14 @@ void phy_sim_debug(phy_sim *sim, size_t *coords) coords[2], offs1); printf("CE = [%E %E %E]\n", aux->curl_E[0], aux->curl_E[1], aux->curl_E[2]); - printf("CH = [%E %E Ef]\n", aux->curl_H[0], aux->curl_H[1], + printf("CH = [%E %E %E]\n", aux->curl_H[0], aux->curl_H[1], aux->curl_H[2]); printf("ICE = [%E %E %E]\n", aux->int_curl_E[0], aux->int_curl_E[1], aux->int_curl_E[2]); printf("ICH = [%E %E %E]\n", aux->int_curl_H[0], aux->int_curl_H[1], aux->int_curl_H[2]); + printf("FC = [%E %E %E]\n", aux->fake_cond[0], aux->fake_cond[1], + aux->fake_cond[2]); printf("=======================\n"); } diff --git a/src/physics.h b/src/physics.h index 2557d3c..06c5d4f 100644 --- a/src/physics.h +++ b/src/physics.h @@ -10,7 +10,7 @@ #define IMPEDANCE_ZERO (sqrt(MU_ZERO/EPS_ZERO)) typedef struct { - size_t width, height, depth; + size_t dims[3]; size_t zstr, ystr, xstr, size; // strides in no. of points (not bytes) size_t zstr1, ystr1, xstr1, size1; // same as above but for aux float spacing; // in meters @@ -24,6 +24,8 @@ typedef struct { typedef struct { float eps, mu; // permittivity and permeability + float fake_cond[3]; + vec3_t curl_E, curl_H; vec3_t int_E, int_H; @@ -43,6 +45,7 @@ typedef struct { phy_field_info field_info; phy_field_em fields[3]; phy_field_aux_point *aux; + size_t pml_widths[3]; float time, time_delta; // UI stuff diff --git a/src/renderer.c b/src/renderer.c index 9cf8db3..94f0333 100644 --- a/src/renderer.c +++ b/src/renderer.c @@ -376,14 +376,14 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, return 0; // nothing's changed if (flags & XSECTION_XY) { - width = fi->width; - height = fi->height; + width = fi->dims[0]; + height = fi->dims[1]; } else if (flags & XSECTION_XZ) { - width = fi->width; - height = fi->depth; + width = fi->dims[0]; + height = fi->dims[2]; } else { - width = fi->height; - height = fi->depth; + width = fi->dims[1]; + height = fi->dims[2]; } if (flags & XSECTION_E) { @@ -408,12 +408,12 @@ 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->depth; - if (z >= fi->depth) - z = fi->depth - 1; + z = frac * fi->dims[2]; + if (z >= fi->dims[2]) + z = fi->dims[2] - 1; - for (y = 0; y < fi->height; y++) - for (x = 0; x < fi->width; x++) { + for (y = 0; y < fi->dims[1]; y++) + for (x = 0; x < fi->dims[0]; x++) { uint8_t *pixel; float *point; @@ -426,12 +426,12 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, pixel[2] = r_float_to_u8(point[2] * scale); } } else if (flags & XSECTION_XZ) { - y = frac * fi->height; - if (y >= fi->height) - y = fi->height - 1; + y = frac * fi->dims[1]; + if (y >= fi->dims[1]) + y = fi->dims[1] - 1; - for (z = 0; z < fi->depth; z++) - for (x = 0; x < fi->width; x++) { + for (z = 0; z < fi->dims[2]; z++) + for (x = 0; x < fi->dims[0]; x++) { uint8_t *pixel; float *point; @@ -444,12 +444,12 @@ int r_xsection_update(r_window *rw, r_xsection *xsection, phy_sim *sim, pixel[2] = r_float_to_u8(point[2] * scale); } } else { - x = frac * fi->width; - if (x >= fi->width) - x = fi->width - 1; + x = frac * fi->dims[0]; + if (x >= fi->dims[0]) + x = fi->dims[0] - 1; - for (z = 0; z < fi->depth; z++) - for (y = 0; y < fi->height; y++) { + for (z = 0; z < fi->dims[2]; z++) + for (y = 0; y < fi->dims[1]; y++) { uint8_t *pixel; float *point; diff --git a/src/ui.c b/src/ui.c index f485cfe..4245603 100644 --- a/src/ui.c +++ b/src/ui.c @@ -268,9 +268,9 @@ void ui_event_window(SDL_Event *event, ui_window *uiw) coords = calloc(3, sizeof(size_t)); assert(coords); - coords[0] = sv->select[0] * fi->width; - coords[1] = sv->select[1] * fi->height; - coords[2] = sv->select[2] * fi->depth; + coords[0] = sv->select[0] * fi->dims[0]; + coords[1] = sv->select[1] * fi->dims[1]; + coords[2] = sv->select[2] * fi->dims[2]; itc_chan_push(&sv->sim->ctl, PHY_CMD_DEBUG, coords); @@ -464,9 +464,9 @@ void ui_draw_simview_bars(ui_window *uiw, float dt) r_draw_rect(uiw->rw, 0, uiw->h - 4 * theme.font_size, uiw->w, 3 * theme.font_size, theme.color_main); - x = sv->select[0] * fi->width; - y = sv->select[1] * fi->height; - z = sv->select[2] * fi->depth; + x = sv->select[0] * fi->dims[0]; + y = sv->select[1] * fi->dims[1]; + z = sv->select[2] * fi->dims[2]; offs = z * fi->zstr + y * fi->ystr + x * fi->xstr; r_draw_text(uiw->rw, 0, uiw->h - 4 * theme.font_size, -- cgit