diff --git a/anisotropic_test.py b/anisotropic_test.py index 1ba3f99..7d17a58 100644 --- a/anisotropic_test.py +++ b/anisotropic_test.py @@ -2,17 +2,30 @@ from scattering import * -mesh_file = "mesh/homo.msh" -epsilon = [[5.47781441, 0. ], [0., 5.78054277]] -permittivity_dict = {1: epsilon, 2: epsilon, 3: np.identity(2)} -s = np.array([1, 2]) -p = np.array([-2, 1]) -k0L = np.pi +def run_anisotropic_test(mesh_file, s, p): + dim = s.shape[0] + epsilon = { + 2: [[5.47781441, 0. ], [0., 5.78054277]], + 3: [[5.4844069, 0., 0.], [0., 5.78683953, 0.], [0., 0., 7.33240001]] + } + permittivity_dict = {1: epsilon[dim], 2: epsilon[dim], 3: np.identity(dim)} + k0L = np.pi -problem = AnisotropicScattering(mesh_file, permittivity_dict, k0L) -pw = PlaneWave(s, p) + problem = AnisotropicScattering(mesh_file, permittivity_dict, k0L) + pw = PlaneWave(s, p) -E = problem.solve(pw) + E = problem.solve(pw) + save_field(E, str(dim) + "D_anisotropic_field.pvd") -#phi, FF = problem.get_far_field(E, 40) -#plot_far_field(phi, FF, "anisotropic") + phi, FF = problem.get_far_field(E, 40) + #plot_far_field(phi, FF, "anisotropic") + +def run_2d_anisotropic_test(): + run_anisotropic_test("mesh/hexa.msh", np.array([-1, -2]), np.array([2, -1])) + +def run_3d_anisotropic_test(): + run_anisotropic_test("mesh/het3D.msh", np.array([-1, -2, 0]), np.array([2, -1, 0])) + +if __name__ == "__main__": + run_2d_anisotropic_test() + run_3d_anisotropic_test() diff --git a/experiment_1.py b/experiment_1.py index 8100db2..80dba12 100644 --- a/experiment_1.py +++ b/experiment_1.py @@ -24,8 +24,8 @@ def run_command(args): "mesh/hexa.geo"]) permittivity_dict = {1: 1, 2: 11.8, 3: 1} - s = np.array([1, 2]) - p = np.array([-2, 1]) + s = np.array([-1, -2]) + p = np.array([2, -1]) k0L = np.pi print("Isotropic Scattering with permittivity {} and n {}".format(permittivity_dict, i)) diff --git a/experiment_2.py b/experiment_2.py index 68f9641..7f510df 100644 --- a/experiment_2.py +++ b/experiment_2.py @@ -13,8 +13,8 @@ HOMO_FILE = "mesh/homo.geo" HEXA_MESH = "mesh/hexa.msh" HEXA_FILE = "mesh/hexa.geo" -s = np.array([1, 2]) -p = np.array([-2, 1]) +s = np.array([-1, -2]) +p = np.array([2, -1]) pw = PlaneWave(s, p) k0L = np.pi diff --git a/isotropic_test.py b/isotropic_test.py index 127c4a0..000d4b1 100644 --- a/isotropic_test.py +++ b/isotropic_test.py @@ -1,16 +1,27 @@ import numpy as np from scattering import * -mesh_file = "mesh/hexa.msh" -permittivity_dict = {1: 1, 2: 11.8, 3: 1} -s = np.array([1, 2]) -p = np.array([-2, 1]) -k0L = np.pi +def run_isotropic_test(mesh_file, s, p): + dim = s.shape[0] + permittivity_dict = {1: 1, 2: 11.8, 3: 1} + k0L = np.pi -problem = IsotropicScattering(mesh_file, permittivity_dict, k0L) -pw = PlaneWave(s, p) + problem = IsotropicScattering(mesh_file, permittivity_dict, k0L) + pw = PlaneWave(s, p) -E = problem.solve(pw) + E = problem.solve(pw) -#phi, FF = problem.get_far_field(E, 40) -#plot_far_field(phi, FF, "isotropic") + save_field(E, str(dim) + "D_field.pvd") + phi, FF = problem.get_far_field(E, 40) + #plot_far_field(phi, FF, "isotropic") + +def run_2d_isotropic_test(): + run_isotropic_test("mesh/hexa.msh", np.array([-1, -2]), np.array([2, -1])) + + +def run_3d_isotropic_test(): + run_isotropic_test("mesh/het3D.msh", np.array([-1, -2, 0]), np.array([2, -1, 0])) + +if __name__ == "__main__": + run_2d_isotropic_test() + run_3d_isotropic_test() diff --git a/mesh/anisotropic.h5 b/mesh/anisotropic.h5 deleted file mode 100644 index 90614f5..0000000 Binary files a/mesh/anisotropic.h5 and /dev/null differ diff --git a/mesh/het3D.geo b/mesh/het3D.geo new file mode 100644 index 0000000..d3eae31 --- /dev/null +++ b/mesh/het3D.geo @@ -0,0 +1,163 @@ +n = 6; + +/********************************************************************* + * + * HEXA square unit cell in 2D inside circle + * + *********************************************************************/ + +// parameters +a = 0.4; +kor = 0; +v = a * Sqrt( 3 ) / 2; + +// Extrude Height +hg = 0.1; +alt = 1; + +// Scaling parameter +h = 1 / n; + +// Meshing parameters +lc1 = 20 * 1e-2 * h; +lc2 = 120 * 1e-2; + + +// Macro defining unit cell +Macro HEXA + + // Points in (x, y, 0) Plane + p1 = newp; Point(p1) = {h*(x0 - a/2) + t1 * h -0.5, h*(y0 - v + kor) + t2 * h + alt * h/2 -0.5, -0.5, lc1}; + p2 = newp; Point(p2) = {h*(x0 + a/2) + t1 * h -0.5, h*(y0 - v + kor) + t2 * h + alt * h/2 -0.5, -0.5, lc1}; + p3 = newp; Point(p3) = {h*(x0 + a) + t1 * h -0.5, h*(y0) + t2 * h + alt * h/2 -0.5, -0.5, lc1}; + p4 = newp; Point(p4) = {h*(x0 + a/2) + t1 * h -0.5, h*(y0 + v - kor) + t2 * h + alt * h/2 -0.5, -0.5, lc1}; + p5 = newp; Point(p5) = {h*(x0 - a/2) + t1 * h -0.5, h*(y0 + v - kor) + t2 * h + alt * h/2 -0.5, -0.5, lc1}; + p6 = newp; Point(p6) = {h*(x0 - a) + t1 * h -0.5, h*(y0) + t2 * h + alt * h/2 -0.5, -0.5, lc1}; + + + // Lines in (x, y, 0) Plane + l1 = newl; Line(l1) = {p1, p2}; + l2 = newl; Line(l2) = {p2, p3}; + l3 = newl; Line(l3) = {p3, p4}; + l4 = newl; Line(l4) = {p4, p5}; + l5 = newl; Line(l5) = {p5, p6}; + l6 = newl; Line(l6) = {p6, p1}; + + // Line loop in (x, y, 0) Plane + ll = newreg; Line Loop(ll) = {l1, l2, l3, l4, l5, l6}; + lloops[t1 * n + t2] = ll; + + // Surface in (x, y, 0) Plane + s = news; Plane Surface(s) = {ll}; + surfaces[t1 * n + t2] = s; + +Return + +// CALL MACRO +x0 = .5; y0 = 0.5; + +For t1 In {0 : n - 1} + alt = t1 % 2; + + For t2 In {0 : n - 1} + Call HEXA; + + EndFor +EndFor + + +// Surface in (x, y, 0) Plane +s = news; Plane Surface(s) = {lloops[]}; + + +Color Blue{ Surface{ surfaces[] }; } + + +// Si Material Box + +// Si Points in (x, y, 0) Plane +p10 = newp; Point(p10) = {-0.5, -0.5, -0.5, lc1}; +p20 = newp; Point(p20) = {0.5, -0.5, -0.5, lc1}; +p30 = newp; Point(p30) = {0.5, 0.5 + h / 2, -0.5, lc1}; +p40 = newp; Point(p40) = {-0.5, 0.5 + h / 2, -0.5, lc1}; + +// Si Lines in (x, y, 0) Plane +l10 = newl; Line(l10) = {p10, p20}; +l20 = newl; Line(l20) = {p20, p30}; +l30 = newl; Line(l30) = {p30, p40}; +l40 = newl; Line(l40) = {p40, p10}; + +// Si Line Loop in (x, y, 0) Plane +llSi = newreg; Line Loop(llSi) = {l10, l20, l30, l40}; + +// Surface in (x, y, 0) Plane +sSi = news; Plane Surface(sSi) = {llSi, -lloops[]}; +Color Red{ Surface{ sSi }; } + + +For t1 In {0 : n - 1} + For t2 In {0 : n - 1} + out[] = Extrude{0, 0, hg}{ Surface{ surfaces[t1 * n + t2] }; }; + phy[t1 * n + t2] = out[1]; + top[t1 * n + t2] = out[0]; + + EndFor +EndFor + + + +outSi[] = Extrude{0, 0, hg}{ Surface{ sSi }; }; + +Coherence Mesh; + +//Outer surface Loop +osl = newreg; Surface Loop(osl) = {top[], surfaces[], sSi, outSi[0], outSi[2], outSi[3], outSi[4], outSi[5]}; + + +Color Blue{ Volume{ phy[] }; } +Color Yellow{ Volume{ outSi[1] }; } + +// Volume 1 is air +Physical Volume(1) = {phy[]}; + +// Volume 2 is Si +Physical Volume(2) = {outSi[1]}; + + +// Outer Box +x = 0; +y = 0; +z = 0; +r = 4; + +p1 = newp; Point(p1) = {x, y, z, lc2} ; +p2 = newp; Point(p2) = {x+r,y, z, lc2} ; +p3 = newp; Point(p3) = {x, y+r,z, lc2} ; +p4 = newp; Point(p4) = {x, y, z+r, lc2} ; +p5 = newp; Point(p5) = {x-r,y, z, lc2} ; +p6 = newp; Point(p6) = {x, y-r,z, lc2} ; +p7 = newp; Point(p7) = {x, y, z-r, lc2} ; + +c1 = newreg; Circle(c1) = {p2,p1,p7}; c2 = newreg; Circle(c2) = {p7,p1,p5}; +c3 = newreg; Circle(c3) = {p5,p1,p4}; c4 = newreg; Circle(c4) = {p4,p1,p2}; +c5 = newreg; Circle(c5) = {p2,p1,p3}; c6 = newreg; Circle(c6) = {p3,p1,p5}; +c7 = newreg; Circle(c7) = {p5,p1,p6}; c8 = newreg; Circle(c8) = {p6,p1,p2}; +c9 = newreg; Circle(c9) = {p7,p1,p3}; c10 = newreg; Circle(c10) = {p3,p1,p4}; +c11 = newreg; Circle(c11) = {p4,p1,p6}; c12 = newreg; Circle(c12) = {p6,p1,p7}; + +l1 = newreg; Line Loop(l1) = {c5,c10,c4}; Surface(newreg) = {l1}; +l2 = newreg; Line Loop(l2) = {c9,-c5,c1}; Surface(newreg) = {l2}; +l3 = newreg; Line Loop(l3) = {c12,-c8,-c1}; Surface(newreg) = {l3}; +l4 = newreg; Line Loop(l4) = {c8,-c4,c11}; Surface(newreg) = {l4}; +l5 = newreg; Line Loop(l5) = {-c10,c6,c3}; Surface(newreg) = {l5}; +l6 = newreg; Line Loop(l6) = {-c11,-c3,c7}; Surface(newreg) = {l6}; +l7 = newreg; Line Loop(l7) = {-c2,-c7,-c12}; Surface(newreg) = {l7}; +l8 = newreg; Line Loop(l8) = {-c6,-c9,c2}; Surface(newreg) = {l8}; + +loop = newreg; Surface Loop(loop) = {l8 + 1, l5 + 1, l1 + 1, l2 + 1, l3 + 1, l7 + 1, l6 + 1, l4 + 1}; +slmb = newreg; Volume(slmb) = {loop, osl}; + +Coherence Mesh; + +// Volume 3 is outer box air +Physical Volume(3) = slmb; diff --git a/mesh/hom3D.geo b/mesh/hom3D.geo new file mode 100644 index 0000000..3803565 --- /dev/null +++ b/mesh/hom3D.geo @@ -0,0 +1,97 @@ +n = 6; + +/********************************************************************* + * + * HEXA square unit cell in 2D inside circle + * + *********************************************************************/ + +// parameters +a = 0.4; +kor = 0; +v = a * Sqrt( 3 ) / 2; + +// Extrude Height +hg = 0.1; +alt = 1; + +// Scaling parameter +h = 1 / n; + +// Meshing parameters +lc1 = 20 * 1e-2; +lc2 = 120 * 1e-2; + +// Si Material Box + +// Si Points in (x, y, 0) Plane +p10 = newp; Point(p10) = {-0.5, -0.5, -0.5, lc1}; +p20 = newp; Point(p20) = {0.5, -0.5, -0.5, lc1}; +p30 = newp; Point(p30) = {0.5, 0.5 + h / 2, -0.5, lc1}; +p40 = newp; Point(p40) = {-0.5, 0.5 + h / 2, -0.5, lc1}; + +// Si Lines in (x, y, 0) Plane +l10 = newl; Line(l10) = {p10, p20}; +l20 = newl; Line(l20) = {p20, p30}; +l30 = newl; Line(l30) = {p30, p40}; +l40 = newl; Line(l40) = {p40, p10}; + +// Si Line Loop in (x, y, 0) Plane +llSi = newreg; Line Loop(llSi) = {l10, l20, l30, l40}; + +// Surface in (x, y, 0) Plane +sSi = news; Plane Surface(sSi) = {llSi}; +Color Red{ Surface{ sSi }; } + + +outSi[] = Extrude{0, 0, hg}{ Surface{ sSi }; }; + +Coherence Mesh; + +//Outer surface Loop +osl = newreg; Surface Loop(osl) = {sSi, outSi[0], outSi[2], outSi[3], outSi[4], outSi[5]}; + + +Color Yellow{ Volume{ outSi[1] }; } + +// Volume 20 is Si +Physical Volume(20) = {outSi[1]}; + + +// Outer Box +x = 0; +y = 0; +z = 0; +r = 4; + +p1 = newp; Point(p1) = {x, y, z, lc2} ; +p2 = newp; Point(p2) = {x+r,y, z, lc2} ; +p3 = newp; Point(p3) = {x, y+r,z, lc2} ; +p4 = newp; Point(p4) = {x, y, z+r, lc2} ; +p5 = newp; Point(p5) = {x-r,y, z, lc2} ; +p6 = newp; Point(p6) = {x, y-r,z, lc2} ; +p7 = newp; Point(p7) = {x, y, z-r, lc2} ; + +c1 = newreg; Circle(c1) = {p2,p1,p7}; c2 = newreg; Circle(c2) = {p7,p1,p5}; +c3 = newreg; Circle(c3) = {p5,p1,p4}; c4 = newreg; Circle(c4) = {p4,p1,p2}; +c5 = newreg; Circle(c5) = {p2,p1,p3}; c6 = newreg; Circle(c6) = {p3,p1,p5}; +c7 = newreg; Circle(c7) = {p5,p1,p6}; c8 = newreg; Circle(c8) = {p6,p1,p2}; +c9 = newreg; Circle(c9) = {p7,p1,p3}; c10 = newreg; Circle(c10) = {p3,p1,p4}; +c11 = newreg; Circle(c11) = {p4,p1,p6}; c12 = newreg; Circle(c12) = {p6,p1,p7}; + +l1 = newreg; Line Loop(l1) = {c5,c10,c4}; Surface(newreg) = {l1}; +l2 = newreg; Line Loop(l2) = {c9,-c5,c1}; Surface(newreg) = {l2}; +l3 = newreg; Line Loop(l3) = {c12,-c8,-c1}; Surface(newreg) = {l3}; +l4 = newreg; Line Loop(l4) = {c8,-c4,c11}; Surface(newreg) = {l4}; +l5 = newreg; Line Loop(l5) = {-c10,c6,c3}; Surface(newreg) = {l5}; +l6 = newreg; Line Loop(l6) = {-c11,-c3,c7}; Surface(newreg) = {l6}; +l7 = newreg; Line Loop(l7) = {-c2,-c7,-c12}; Surface(newreg) = {l7}; +l8 = newreg; Line Loop(l8) = {-c6,-c9,c2}; Surface(newreg) = {l8}; + +loop = newreg; Surface Loop(loop) = {l8 + 1, l5 + 1, l1 + 1, l2 + 1, l3 + 1, l7 + 1, l6 + 1, l4 + 1}; +slmb = newreg; Volume(slmb) = {loop, osl}; + +Coherence Mesh; + +// Volume 30 is outer box air +Physical Volume(30) = slmb; diff --git a/mesh/isotropic.h5 b/mesh/isotropic.h5 deleted file mode 100644 index 53447e0..0000000 Binary files a/mesh/isotropic.h5 and /dev/null differ diff --git a/scattering/anisotropic_scattering.py b/scattering/anisotropic_scattering.py index 0562d6c..46a7c9f 100644 --- a/scattering/anisotropic_scattering.py +++ b/scattering/anisotropic_scattering.py @@ -1,15 +1,22 @@ -import firedrake as fd +import numpy as np +from firedrake import as_matrix, Function, TensorFunctionSpace from .scattering import Scattering class AnisotropicScattering(Scattering): - def __init__(self, mesh, permittivity_dict, k0L, **kvargs): - super().__init__(mesh, k0L, **kvargs) - self.II = fd.as_matrix(((1, 0), (0,1))) - self.permittivity = fd.Function(fd.TensorFunctionSpace(self.mesh, "DG", 0)) + """ + Derived class for AnisotropicScattering + + Computation using this object models material as an anisotropic scatterer. Permittivity of + each material is modeled using tensor. + """ + def __init__(self, mesh, permittivity_dict, k0L, **kwargs): + super().__init__(mesh, k0L, **kwargs) + # identity here is matrix + self.II = as_matrix(np.identity(self.mesh.topological_dimension())) + self.permittivity = Function(TensorFunctionSpace(self.mesh, "DG", 0)) for (subd_id, epsilon_tensor) in permittivity_dict.items(): - epsilon = fd.as_matrix(epsilon_tensor) - self.permittivity.interpolate(epsilon, + self.permittivity.interpolate(as_matrix(epsilon_tensor), self.mesh.measure_set("cell", subd_id)) diff --git a/scattering/isotropic_scattering.py b/scattering/isotropic_scattering.py index e445a99..8544962 100644 --- a/scattering/isotropic_scattering.py +++ b/scattering/isotropic_scattering.py @@ -3,11 +3,19 @@ class IsotropicScattering(Scattering): - def __init__(self, mesh, permittivity_dict, k0L, **kvargs): - super().__init__(mesh, k0L, **kvargs) + """ + Derived class for IsotropicScattering + + Computation using this object models material as in isotropic scatterer. Permittivity of + each material is modeled using one real number. + """ + def __init__(self, mesh, permittivity_dict, k0L, **kwargs): + super().__init__(mesh, k0L, **kwargs) + # identity for isotropic problem is just 1 self.II = 1 + # permittivity is DG0 function that changes between each subdomain self.permittivity = fd.Function(fd.FunctionSpace(self.mesh, "DG", 0)) - for (subdomain_id, epsilon) in permittivity_dict.items(): + for (subd_id, epsilon) in permittivity_dict.items(): self.permittivity.interpolate(fd.Constant(epsilon), - self.mesh.measure_set("cell", subdomain_id)) + self.mesh.measure_set("cell", subd_id)) diff --git a/scattering/plane_wave.py b/scattering/plane_wave.py index 34f7b34..50b8c78 100644 --- a/scattering/plane_wave.py +++ b/scattering/plane_wave.py @@ -1,13 +1,26 @@ -import firedrake as fd import numpy as np +from abc import ABC, abstractmethod + +from firedrake import VectorFunctionSpace, SpatialCoordinate, interpolate, dot, exp, as_vector + +from .scattering import ExcitationBase def orthogonal(s, p): + """Return true if np.arrays that represent vectors are orthogonal""" return np.isclose(np.dot(s, p), 0) def normalize(a): + """Normalize vector a represented as np.array""" return a / np.linalg.norm(a) -class PlaneWave(): +class PlaneWave(ExcitationBase): + """ + Plane wave excitation + + We define plane wave as: + E_i(x) = \hat{p} * exp(i k0L * \hat{s} \cdot x) + """ + def __init__(self, s, p): # If vectors aren't orthogonal, make them orthogonal by ignoring # electric field component in the direction of propagation @@ -20,14 +33,13 @@ def __init__(self, s, p): s = normalize(s) p = normalize(p) - self.s = fd.as_vector(s) - self.p = fd.as_vector(p) - - #TODO: this should be made abtract and other excitements should implement - def interpolate(self, mesh, k0L): - V = fd.VectorFunctionSpace(mesh, "CG", 1) - x = fd.SpatialCoordinate(mesh) + self.s = as_vector(s) + self.p = as_vector(p) - pw = fd.interpolate(self.p * fd.exp(1j * k0L * fd.dot(self.s, x)), V) + def interpolate(self, mesh, *args): + """Interpolate plane wave on the mesh""" + (k0L,) = args + V = VectorFunctionSpace(mesh, "CG", 1) + x = SpatialCoordinate(mesh) - return pw + return interpolate(self.p * exp(1j * k0L * dot(self.s, x)), V) diff --git a/scattering/scattering.py b/scattering/scattering.py index 154f77f..ba8175c 100644 --- a/scattering/scattering.py +++ b/scattering/scattering.py @@ -2,22 +2,69 @@ import numpy as np from abc import ABC, abstractmethod +class ExcitationBase(ABC): + """Abstract Excitation object""" + + @abstractmethod + def interpolate(self, mesh, *args): + """Interpolate Excitation on the mesh""" + pass + + class Scattering(ABC): - def __init__(self, mesh, k0L, output_dir="output"): + """ + Base class for all Scattering problems. + Scattering object is wrapper above problem description and domain where problem + is defined. + """ + + def __init__(self, mesh, k0L): + """ + Initialize Scattering object + + Parameters: + ----------- + mesh: either path to the file with extension of some supported mesh format[1] or + already created firedrake mesh. One should use firedrake mesh in order to + do a errornorm between two solutions on the same mesh. + k0L: wave number + + [1] - https://www.firedrakeproject.org/firedrake.html#firedrake.mesh.Mesh + """ if isinstance(mesh, str): self.mesh = fd.Mesh(mesh) elif isinstance(mesh, fd.firedrake.mesh.MeshGeometry): self.mesh = mesh else: raise ValueError("Mesh object unknown") - self.output_dir = output_dir self.k0L = k0L def solve(self, incident, method="lu"): - Ei = incident.interpolate(self.mesh, self.k0L) + """ + Solve Maxwell Scattering problem using incident excitation + + We are solving: + + curl(curl(E)) - k^2 epsilon_r * E = 0 + E = E_i + E_s where E_i is incident wave + + As a boundary conditions we are using first order Silver-Muller BC. + + Function will compute E_s and return E = E_i + E_s. - V = fd.FunctionSpace(self.mesh, 'N1curl', 1) + Parameters: + ----------- + incident: object that should inherit from ExcitationBase + method: method for solving Ax=b system. Currently we default to MUMPS. + """ + if isinstance(incident, ExcitationBase): + Ei = incident.interpolate(self.mesh, self.k0L) + else: + raise ValueError("incident object unknown") + + V_el = fd.FiniteElement('N1curl', self.mesh.ufl_cell(), 1, variant="integral") + V = fd.FunctionSpace(self.mesh, V_el) Es = fd.Function(V) u = fd.TrialFunction(V) v = fd.TestFunction(V) @@ -29,8 +76,13 @@ def solve(self, incident, method="lu"): ik = (1j * k) a = (fd.inner(fd.curl(u), fd.curl(v)) * fd.dx - - k**2 * fd.inner(epsilon * u, v) * fd.dx - - ik * ((fd.dot(u, n) * fd.inner(n, v)) - fd.inner(u, v)) * fd.ds) + - k**2 * fd.inner(epsilon * u, v) * fd.dx) + if self.mesh.topological_dimension() == 3: + a -= ik * fd.inner((fd.cross(n, fd.cross(n, u))), v) * fd.ds + else: + # for 2D case we have to write a \times b \times c as: + # (a \cdot c) \cdot b - (a \cdot b) \cdot c + a -= ik * ((fd.dot(n, u) * fd.inner(n, v)) - fd.inner(u, v)) * fd.ds L = - k**2 * fd.inner((self.II - epsilon) * Ei, v) * fd.dx solvers = { @@ -41,6 +93,7 @@ def solve(self, incident, method="lu"): "mat_mumps_icntl_4": "1"}, } + #TODO: LinearVariationalSolver A = fd.assemble(a) b = fd.assemble(L) @@ -50,12 +103,26 @@ def solve(self, incident, method="lu"): return fd.project(Es + Ei, V) - def get_far_field(self, E, FF_n): - phi = np.linspace(0, 2 * np.pi, num = FF_n, endpoint = False) - FF = np.zeros(FF_n) + def get_far_field(self, E, ff_num_samples, theta=np.pi/2): + """Compute Far field for the electric field - cos_values = np.cos(phi) - sin_values = np.sin(phi) + Compute far field components using Mischenkos: Electromagnetic scattering by Particles + and Particle Groups. Cambridge University Press. + + Parameters: + ----------- + E: electromagnetic field for which we are computing far field. + ff_num_samples: number of far field components. In other word, how much values will + we sample along the circle. Points on the circle are equidistant. + theta: angle with respect to polar axis. pi/2 generates xOy plane. + + Returns: + -------- + phi: np.array of ff_num_samples angles where we computed far field. + ff_data: np.array of ff_num_samples which represent far field at each angle phi. + """ + phi = np.linspace(0, 2 * np.pi, num = ff_num_samples, endpoint = False) + ff_data = np.zeros(ff_num_samples) V3 = fd.VectorFunctionSpace(self.mesh, "CG", 1) V = fd.FunctionSpace(self.mesh, "CG", 1) @@ -65,23 +132,31 @@ def get_far_field(self, E, FF_n): epsilon = self.permittivity tdim = self.mesh.topological_dimension() - r = fd.Constant([1, 0]) + if tdim == 2: + r = fd.Constant([1, 0]) + else: + r = fd.Constant([1, 0, 0]) - for n in range(FF_n): - rx = cos_values[n] - ry = sin_values[n] + for n in range(ff_num_samples): + if tdim == 2: + r_vals = [np.cos(phi[n]), + np.sin(phi[n])] + else: + r_vals = [np.sin(theta) * np.cos(phi[n]), + np.sin(theta) * np.sin(phi[n]), + np.cos(theta)] - r_vals = np.array([rx, ry]) r.assign(r_vals) A = np.identity(tdim) - np.outer(r_vals, r_vals) - f = fd.exp(1j * self.k0L * fd.dot(r, x)) + f = fd.exp(-1j * self.k0L * fd.dot(r, x)) #TODO: don't integrate if integral is zero ffi = fd.assemble(fd.inner((epsilon - self.II) * E * f, v) * fd.dx) # here we have [real_0 + imag_0, ..., real_d + imag_d] - ff_components = np.sum(ffi.dat.data_ro, axis=0) + ff_components = np.sum(ffi.vector().gather().reshape(-1, tdim), axis=0) - FF[n] = np.linalg.norm(k**2 /(4 * np.pi) * A.dot(ff_components)) + # ff_components = [sum_along_0, sum_along_1, ..., sum_along_d] + ff_data[n] = np.linalg.norm(k**2 /(4 * np.pi) * A.dot(ff_components)) - return phi, FF + return phi, ff_data