From 6dd784058ae8e45b306311bd303d09a32113f50e Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Sun, 3 Apr 2016 21:56:12 +0200 Subject: PML. --- src/common.h | 1 + src/physics.c | 128 +++++++++++++++++++++++++++++++++++++++++++++++----------- src/physics.h | 7 ++-- src/ui.c | 58 ++++++++++++++++++-------- 4 files changed, 151 insertions(+), 43 deletions(-) (limited to 'src') diff --git a/src/common.h b/src/common.h index 743dea4..9dc85e9 100644 --- a/src/common.h +++ b/src/common.h @@ -18,6 +18,7 @@ #define max3(a, b, c) (max2(max2(a, b), c)) #define min2(a, b) ((a) < (b) ? (a) : (b)) #define min3(a, b, c) (min2(min2(a, b), c)) +#define sq(x) ((x)*(x)) char *va(const char *fmt, ...); int64_t get_time(void); diff --git a/src/physics.c b/src/physics.c index bbe2499..952e53a 100644 --- a/src/physics.c +++ b/src/physics.c @@ -16,8 +16,8 @@ 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) +static void phy_calc_fc(float *fc, size_t *x, size_t *widths, + size_t *dims, float dt) { size_t i; @@ -31,31 +31,88 @@ static float phy_calc_fake_cond(float *fake_cond, size_t *x, size_t *widths, else continue; - fake_cond[i] = EPS_ZERO / (2 * dt) * pow(f, 3); + fc[i] = EPS_ZERO / (2 * dt) * pow(f, 3); } } +static void phy_calc_uec_H(float *uec, float dt, float fc0, float fc1, + float fc2, float mu) +{ + float tmp; + + tmp = 1.0f / dt + (fc1 + fc2) / (2 * EPS_ZERO) + + fc1 * fc2 * dt / sq(2.0f * EPS_ZERO); + uec[0] = 1.0f / tmp * (1.0f / dt - (fc1 + fc2) / (2 * EPS_ZERO) + - fc1 * fc2 * dt / sq(2.0f * EPS_ZERO)); + uec[1] = -1.0f / tmp * LIGHT_SPEED / mu; + uec[2] = -1.0f / tmp * LIGHT_SPEED * dt / EPS_ZERO * fc0 / mu; + uec[3] = -1.0f / tmp * dt / sq(EPS_ZERO) * fc1 * fc2; +} + +static void phy_calc_uec_E(float *uec, float dt, float fc0, float fc1, + float fc2, float eps) +{ + float tmp; + + tmp = 1.0f / dt + (fc1 + fc2) / (2 * EPS_ZERO) + + fc1 * fc2 * dt / sq(2.0f * EPS_ZERO); + uec[0] = 1.0f / tmp * (1.0f / dt - (fc1 + fc2) / (2 * EPS_ZERO) + - fc1 * fc2 * dt / sq(2.0f * EPS_ZERO)); + uec[1] = LIGHT_SPEED / tmp; + uec[2] = 1.0f / tmp * LIGHT_SPEED * dt * fc0 / EPS_ZERO; + uec[3] = 1.0f / tmp * dt / sq(EPS_ZERO) * fc1 * fc2; + + // equations above are for the D field and E = 1/eps * D + /* + uec[0] *= 1.0f / eps; + uec[1] *= 1.0f / eps; + uec[2] *= 1.0f / eps; + uec[3] *= 1.0f / eps;*/ +} + void phy_sim_reset_aux(phy_sim *sim) { - size_t x[3], i; + size_t x[3]; phy_field_info *fi = &sim->field_info; + float dt = sim->time_delta; + + memset(sim->aux, 0, fi->size1 * sizeof(phy_field_aux_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; + vec3_t fc = {0, 0, 0}; aux = sim->aux + x[2] * fi->zstr1 + x[1] * fi->ystr1 + x[0] * fi->xstr1; - aux->eps = 1.0f / EPS_ZERO; - aux->mu = -1.0f / MU_ZERO; + aux->eps = EPS_ZERO; + aux->mu = MU_ZERO; v3_set(aux->int_curl_E, 0, 0, 0); v3_set(aux->int_curl_H, 0, 0, 0); - phy_calc_fake_cond(aux->fake_cond, x, sim->pml_widths, - fi->dims, sim->time_delta); + phy_calc_fc(fc, x, sim->pml_widths, fi->dims, dt); + /* + for (size_t i = 0; i < 3; i++) { + aux->uec_H[i][0] = 1; + aux->uec_H[i][1] = -dt / MU_ZERO; + aux->uec_H[i][2] = 0; + aux->uec_H[i][3] = 0; + aux->uec_E[i][0] = 1; + aux->uec_E[i][1] = +dt / EPS_ZERO; + aux->uec_E[i][2] = 0; + aux->uec_E[i][3] = 0; + }*/ + + phy_calc_uec_H(aux->uec_H[0], dt, fc[0], fc[1], fc[2], 1); + phy_calc_uec_H(aux->uec_H[1], dt, fc[1], fc[0], fc[2], 1); + phy_calc_uec_H(aux->uec_H[2], dt, fc[2], fc[0], fc[1], 1); + + phy_calc_uec_E(aux->uec_E[0], dt, fc[0], fc[1], fc[2], 1); + phy_calc_uec_E(aux->uec_E[1], dt, fc[1], fc[0], fc[2], 1); + phy_calc_uec_E(aux->uec_E[2], dt, fc[2], fc[0], fc[1], 1); } } @@ -96,9 +153,15 @@ int phy_sim_create(phy_sim *sim) printf("phy_sim_create: Courant number is %f\n", 3 * LIGHT_SPEED * sim->time_delta / fi->spacing); + #if 1 sim->pml_widths[0] = fi->dims[0] / 5; sim->pml_widths[1] = fi->dims[1] / 5; sim->pml_widths[2] = fi->dims[2] / 5; + #else + sim->pml_widths[0] = 0; + sim->pml_widths[1] = 0; + sim->pml_widths[2] = 0; + #endif fi->xstr = 3; fi->ystr = fi->xstr * fi->dims[0]; @@ -207,25 +270,30 @@ void phy_sim_update_em_H(phy_sim *sim) { phy_field_info *fi = &sim->field_info; size_t x, y, z; - float *Hn, *Hp, dt; + float *Hn, *Hp; Hn = sim->fields[2].H; Hp = sim->fields[0].H; - dt = sim->time_delta; 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; + size_t offs, i; 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; - 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]; + for (i = 0; i < 3; i++) + Hn[offs+i] = aux->uec_H[i][0] * Hp[offs+i] + + aux->uec_H[i][1] * aux->curl_E[i] + + aux->uec_H[i][2] * aux->int_curl_E[i] + + aux->uec_H[i][3] * aux->int_H[i]; + + //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]; } } @@ -233,25 +301,30 @@ void phy_sim_update_em_E(phy_sim *sim) { phy_field_info *fi = &sim->field_info; size_t x, y, z; - float *En, *Ep, dt; + float *En, *Ep; En = sim->fields[2].E; Ep = sim->fields[0].E; - dt = sim->time_delta; 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; + size_t offs, i; 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; - 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]; + for (i = 0; i < 3; i++) + En[offs+i] = aux->uec_E[i][0] * Ep[offs+i] + + aux->uec_E[i][1] * aux->curl_H[i] + + aux->uec_E[i][2] * aux->int_curl_H[i] + + aux->uec_E[i][3] * aux->int_E[i]; + + //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]; } } @@ -342,8 +415,9 @@ void phy_sim_debug(phy_sim *sim, size_t *coords) aux = sim->aux + offs1; printf("===[ phy_sim_debug ]===\n"); - printf("x = [%zu %zu %zu], offs1=%zu\n", coords[0], coords[1], - coords[2], offs1); + printf("x = [%zu %zu %zu], offs1=%zu, aux=%p\n", coords[0], coords[1], + coords[2], offs1, aux); + printf("eps = %f, mu = %f\n", aux->eps, aux->mu); printf("CE = [%E %E %E]\n", aux->curl_E[0], aux->curl_E[1], aux->curl_E[2]); printf("CH = [%E %E %E]\n", aux->curl_H[0], aux->curl_H[1], @@ -352,8 +426,14 @@ void phy_sim_debug(phy_sim *sim, size_t *coords) 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("UEC_E = [[%E %E %E %E] [%E %E %E %E] [%E %E %E %E]]\n", + aux->uec_E[0][0], aux->uec_E[0][1], aux->uec_E[0][2], aux->uec_E[0][3], + aux->uec_E[1][0], aux->uec_E[1][1], aux->uec_E[1][2], aux->uec_E[1][3], + aux->uec_E[2][0], aux->uec_E[2][1], aux->uec_E[2][2], aux->uec_E[2][3]); + printf("UEC_H = [[%E %E %E %E] [%E %E %E %E] [%E %E %E %E]]\n", + aux->uec_H[0][0], aux->uec_H[0][1], aux->uec_H[0][2], aux->uec_H[0][3], + aux->uec_H[1][0], aux->uec_H[1][1], aux->uec_H[1][2], aux->uec_H[1][3], + aux->uec_H[2][0], aux->uec_H[2][1], aux->uec_H[2][2], aux->uec_H[2][3]); printf("=======================\n"); } diff --git a/src/physics.h b/src/physics.h index 06c5d4f..09ee23e 100644 --- a/src/physics.h +++ b/src/physics.h @@ -5,7 +5,7 @@ #include "itc.h" #define LIGHT_SPEED (299792458.0f) -#define MU_ZERO (4.0f * M_PI * 10e-7f) +#define MU_ZERO (4.0f * M_PI * 1e-7f) #define EPS_ZERO (1.0f / (MU_ZERO * LIGHT_SPEED * LIGHT_SPEED)) #define IMPEDANCE_ZERO (sqrt(MU_ZERO/EPS_ZERO)) @@ -22,9 +22,10 @@ typedef struct { } phy_field_em; typedef struct { - float eps, mu; // permittivity and permeability + float eps, mu; // for visualisation, not needed for stepping - float fake_cond[3]; + float uec_H[3][4]; + float uec_E[3][4]; vec3_t curl_E, curl_H; diff --git a/src/ui.c b/src/ui.c index 4245603..3f2214a 100644 --- a/src/ui.c +++ b/src/ui.c @@ -24,7 +24,7 @@ struct { typedef struct { phy_sim *sim; - bool dragging, selecting; + bool dragging, dragging_frac, selecting; int xsection_flags; float xsection_frac; bool select_valid; @@ -167,31 +167,54 @@ void ui_event_window(SDL_Event *event, ui_window *uiw) case SDL_MOUSEBUTTONDOWN: if (event->button.y > sv->margin_top && event->button.y < uiw->h - sv->margin_bottom) { - if (event->button.button == SDL_BUTTON_LEFT) + switch (event->button.button) { + case SDL_BUTTON_LEFT: sv->dragging = true; - - if (event->button.button == SDL_BUTTON_RIGHT) + break; + case SDL_BUTTON_RIGHT: + sv->dragging_frac = true; + break; + case SDL_BUTTON_MIDDLE: sv->selecting = true; + break; + } } break; case SDL_MOUSEBUTTONUP: - if (event->button.button == SDL_BUTTON_LEFT) + switch (event->button.button) { + case SDL_BUTTON_LEFT: sv->dragging = false; - else if (event->button.button == SDL_BUTTON_RIGHT) + break; + case SDL_BUTTON_RIGHT: + sv->dragging_frac = false; + break; + case SDL_BUTTON_MIDDLE: sv->selecting = false; - break; + break; + } case SDL_MOUSEMOTION: uiw->mouse[0] = event->motion.x; uiw->mouse[1] = event->motion.y; - if (sv->dragging) { + if (sv->dragging || sv->dragging_frac) { vec2_t delta_v; v2_set(delta_v, event->motion.xrel, event->motion.yrel); v2_div_mst2_nt(delta_v, delta_v, sv->tf_v2s); - v2_add(sv->origin, sv->origin, delta_v); + + if (sv->dragging) + v2_add(sv->origin, sv->origin, delta_v); + + if (sv->dragging_frac) { + sv->xsection_frac += delta_v[1]; + + if (sv->xsection_frac > 1.0f) + sv->xsection_frac = 1.0f; + else if (sv->xsection_frac < 0.0f) + sv->xsection_frac = 0.0f; + } } break; @@ -407,14 +430,17 @@ void ui_draw_simview_selection(ui_window *uiw) char *ui_draw_simview_top_bar_text(ui_simview *sv) { int flags = sv->xsection_flags; + char *rv; + + rv = va("Przekrój %s płaszczyzną %s=%f", + (flags & XSECTION_E ? "E" : + "H"), + (flags & XSECTION_XY ? "Z" : + flags & XSECTION_XZ ? "Y" : + "X"), + sv->xsection_frac); - va("Przekrój %s płaszczyzną %s=%f", - (flags & XSECTION_E ? "E" : - "H"), - (flags & XSECTION_XY ? "Z" : - flags & XSECTION_XZ ? "Y" : - "X"), - sv->xsection_frac); + return rv; } void ui_draw_simview_bars(ui_window *uiw, float dt) -- cgit