From e288e590211380dc937760038dcf7ed89f71068f Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 3 Jul 2016 16:45:38 -0700 Subject: [PATCH] PEC and PMC for all wave operators --- examples/test.py | 12 ++++++--- fdfd_tools/operators.py | 54 ++++++++++++++++++++++++++++++++--------- 2 files changed, 51 insertions(+), 15 deletions(-) diff --git a/examples/test.py b/examples/test.py index 69b726c..b41438a 100644 --- a/examples/test.py +++ b/examples/test.py @@ -177,12 +177,18 @@ def test1(): J = waveguide_mode.compute_source(**wg_args, **wg_results) H_overlap = waveguide_mode.compute_overlap_e(**wg_args, **wg_results) - pecg = gridlock.Grid(edge_coords, initial=0, num_grids=3) + pecg = gridlock.Grid(edge_coords, initial=0.0, num_grids=3) # pecg.draw_cuboid(center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1) - pecg.grids = [numpy.sign(r) for r in pecg.grids] # pecg.visualize_isosurface() - A = fdfd_tools.operators.e_full(omega, dxes, vec(grid.grids), pec=vec(pecg.grids)).tocsr() + pmcg = gridlock.Grid(edge_coords, initial=0.0, num_grids=3) + # pmcg.draw_cuboid(center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1) + # pmcg.visualize_isosurface() + + A = fdfd_tools.operators.e_full(omega, dxes, + epsilon=vec(grid.grids), + pec=vec(pecg.grids), + pmc=vec(pmcg.grids)).tocsr() b = -1j * omega * vec(J) x = solve_A(A, b) E = unvec(x, grid.shape) diff --git a/fdfd_tools/operators.py b/fdfd_tools/operators.py index f3d501a..8045aaa 100644 --- a/fdfd_tools/operators.py +++ b/fdfd_tools/operators.py @@ -47,6 +47,7 @@ def e_full(omega: complex, epsilon: vfield_t, mu: vfield_t = None, pec: vfield_t = None, + pmc: vfield_t = None, ) -> sparse.spmatrix: """ Wave operator del x (1/mu * del x) - omega**2 * epsilon, for use with E-field, @@ -61,6 +62,8 @@ def e_full(omega: complex, :param mu: Vectorized magnetic permeability (default 1 everywhere). :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted as containing a perfect electrical conductor (PEC). + :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted + as containing a perfect magnetic conductor (PMC). :return: Sparse matrix containing the wave operator """ ce = curl_e(dxes) @@ -68,10 +71,15 @@ def e_full(omega: complex, ev = epsilon if numpy.any(numpy.equal(pec, None)): + pe = sparse.eye(epsilon.size) + else: + pe = sparse.diags(numpy.where(pec, 0, 1)) # Set pe to (not PEC) + ev = numpy.where(pec, 1.0, ev) # Set epsilon to 1 at PEC + + if numpy.any(numpy.equal(pmc, None)): pm = sparse.eye(epsilon.size) else: - pm = sparse.diags(numpy.where(pec, 0, 1)) # Set pm to (not PEC) - ev = numpy.where(pec, 1.0, ev) # Set epsilon to 1 at PEC + pm = sparse.diags(numpy.where(pmc, 0, 1)) # set pm to (not PMC) e = sparse.diags(ev) if numpy.any(numpy.equal(mu, None)): @@ -79,7 +87,7 @@ def e_full(omega: complex, else: m_div = sparse.diags(1 / mu) - op = pm @ ch @ m_div @ ce @ pm - omega**2 * e + op = pe @ ch @ pm @ m_div @ ce @ pe - omega**2 * e return op @@ -110,6 +118,7 @@ def h_full(omega: complex, dxes: dx_lists_t, epsilon: vfield_t, mu: vfield_t = None, + pec: vfield_t = None, pmc: vfield_t = None, ) -> sparse.spmatrix: """ @@ -121,6 +130,8 @@ def h_full(omega: complex, :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header :param epsilon: Vectorized dielectric constant :param mu: Vectorized magnetic permeability (default 1 everywhere) + :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted + as containing a perfect electrical conductor (PEC). :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted as containing a perfect magnetic conductor (PMC). :return: Sparse matrix containing the wave operator @@ -134,19 +145,24 @@ def h_full(omega: complex, mv = mu if numpy.any(numpy.equal(pmc, None)): + pm = sparse.eye(epsilon.size) + else: + pm = sparse.diags(numpy.where(pmc, 0, 1)) # Set pe to (not PMC) + mv = numpy.where(pmc, 1.0, mv) # Set mu to 1 at PMC + + if numpy.any(numpy.equal(pec, None)): pe = sparse.eye(epsilon.size) else: - pe = sparse.diags(numpy.where(pmc, 0, 1)) # Set pe to (not PMC) - mv = numpy.where(pmc, 1.0, mv) # Set mu to 1 at PMC + pe = sparse.diags(numpy.where(pec, 0, 1)) # set pe to (not PEC) e_div = sparse.diags(1 / epsilon) m = sparse.diags(mv) - A = pe @ ec @ e_div @ hc @ pe - omega**2 * m + A = pm @ ec @ pe @ e_div @ hc @ pm - omega**2 * m return A -def eh_full(omega, dxes, epsilon, mu=None): +def eh_full(omega, dxes, epsilon, mu=None, pec=None, pmc=None): """ Wave operator for [E, H] field representation. This operator implements Maxwell's equations without cancelling out either E or H. The operator is @@ -159,18 +175,32 @@ def eh_full(omega, dxes, epsilon, mu=None): :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header :param epsilon: Vectorized dielectric constant :param mu: Vectorized magnetic permeability (default 1 everywhere) + :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted + as containing a perfect electrical conductor (PEC). + :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted + as containing a perfect magnetic conductor (PMC). :return: Sparse matrix containing the wave operator """ - A2 = curl_e(dxes) - A1 = curl_h(dxes) + if numpy.any(numpy.equal(pec, None)): + pe = sparse.eye(epsilon.size) + else: + pe = sparse.diags(numpy.where(pec, 0, 1)) # set pe to (not PEC) - iwe = 1j * omega * sparse.diags(epsilon) - iwm = 1j * omega + if numpy.any(numpy.equal(pmc, None)): + pm = sparse.eye(epsilon.size) + else: + pm = sparse.diags(numpy.where(pmc, 0, 1)) # set pm to (not PMC) + + iwe = pe @ (1j * omega * sparse.diags(epsilon)) + iwm = pm * 1j * omega if not numpy.any(numpy.equal(mu, None)): iwm *= sparse.diags(mu) + A1 = pe @ curl_h(dxes) @ pm + A2 = pm @ curl_e(dxes) @ pe + A = sparse.bmat([[-iwe, A1], - [A2, +iwm]]) + [A2, iwm]]) return A