Código ejecutor Electromagnetic Scattering in Microwires

In [ ]:

    

#Preambulo
import numpy as np
import bempp.api
omega = 2.*np.pi*10.e9
e0 = 8.854*1e-12*1e-18
mu0 = 4.*np.pi*1e-7*1e6
mue = (10.)*mu0
ee = (10.-10.j)*e0
mui = (-2.9214+0.5895j)*mu0
ei = (-25589.156-131017.27j)*e0
k = omega*np.sqrt(e0*mu0)
lam = 2*np.pi/k
nm = np.sqrt((ee*mue)/(e0*mu0))
nc = np.sqrt((ei*mui)/(e0*mu0))
alfa_m = mue/mu0
alfa_c = mui/mue
antena = np.array([[1e4],[0.],[0.]])
print "Numero de onda exterior:", k
print "Indice de refraccion matriz:", nm
print "Indice de refraccion conductor:", nc
print "Numero de onda interior matriz:", nm*k
print "Numero de onda interior conductor:", nm*nc*k
print "Indice de transmision matriz:", alfa_m
print "Indice de transmision conductor:", alfa_c
print "Longitud de onda:", lam, "micras"

#Importando mallas
matriz = bempp.api.import_grid('/home/milan/matriz_12x12x300_E16772.msh')

#Funciones de dirichlet y neumann
def dirichlet_fun(x, n, domain_index, result):
    result[0] = 1.*np.exp(1j*k*x[0])
def neumann_fun(x, n, domain_index, result):
    result[0] = 1.*1j*k*n[0]*np.exp(1j*k*x[0])
 
#Operadores multitrazo
Ai_m = bempp.api.operators.boundary.helmholtz.multitrace_operator(matriz, nm*k)
Ae_m = bempp.api.operators.boundary.helmholtz.multitrace_operator(matriz, k)

#Transmision en Multitrazo
Ae_m[0,1] = Ae_m[0,1]*(1./alfa_m)
Ae_m[1,1] = Ae_m[1,1]*(1./alfa_m)

#Acople interior y exterior
op_m = (Ai_m + Ae_m)

#Espacios
dirichlet_space_m = Ai_m[0,0].domain
neumann_space_m = Ai_m[0,1].domain

#Operadores identidad
ident_m = bempp.api.operators.boundary.sparse.identity(neumann_space_m, neumann_space_m, neumann_space_m)

#Operadores diagonales
op_m[1,1] = op_m[1,1] + 0.5 * ident_m * ((alfa_m -1)/alfa_m)
 
#Operadores entre mallas
 
#Matriz de operadores
blocked = bempp.api.BlockedOperator(2,2)
 
#Diagonal
blocked[0,0] = op_m[0,0]
blocked[0,1] = op_m[0,1]
blocked[1,0] = op_m[1,0]
blocked[1,1] = op_m[1,1]

#Contribucion hilos-matriz

#Condiciones de borde
dirichlet_grid_fun_m = bempp.api.GridFunction(dirichlet_space_m, fun=dirichlet_fun)
neumann_grid_fun_m = bempp.api.GridFunction(neumann_space_m, fun=neumann_fun)

#Discretizacion lado izquierdo
blocked_discretizado = blocked.strong_form()

#Discretizacion lado derecho
rhs = np.concatenate([dirichlet_grid_fun_m.coefficients, neumann_grid_fun_m.coefficients,])

#Sistema de ecuaciones
import inspect
from scipy.sparse.linalg import gmres
array_it = np.array([])
array_frame = np.array([])
it_count = 0
def iteration_counter(x):
    global array_it
    global array_frame
    global it_count
    it_count += 1
    frame = inspect.currentframe().f_back
    array_it = np.append(array_it, it_count)
    array_frame = np.append(array_frame, frame.f_locals["resid"])
    print it_count, frame.f_locals["resid"]
print("Shape of matrix: {0}".format(blocked_discretizado.shape))
x,info = gmres(blocked_discretizado, rhs, tol=1e-5, callback = iteration_counter, maxiter = 50000)
print("El sistema fue resuelto en {0} iteraciones".format(it_count))
np.savetxt("Solucion.out", x, delimiter=",")

#Campo interior
interior_field_dirichlet_m = bempp.api.GridFunction(dirichlet_space_m, coefficients=x[:dirichlet_space_m.global_dof_count])
interior_field_neumann_m = bempp.api.GridFunction(neumann_space_m,coefficients=x[dirichlet_space_m.global_dof_count:dirichlet_space_m.globa    l_dof_count + neumann_space_m.global_dof_count])

#Campo exterior
exterior_field_dirichlet_m = interior_field_dirichlet_m
exterior_field_neumann_m = interior_field_neumann_m*(1./alfa_m)

#Calculo campo en antena
slp_pot_ext_m = bempp.api.operators.potential.helmholtz.single_layer(dirichlet_space_m, antena, k)
dlp_pot_ext_m = bempp.api.operators.potential.helmholtz.double_layer(dirichlet_space_m, antena, k)
Campo_en_antena = (dlp_pot_ext_m * exterior_field_dirichlet_m - slp_pot_ext_m * exterior_field_neumann_m).ravel() + np.exp(1j*k*antena[0])
print "Valor del campo en receptor:", Campo_en_antena

import matplotlib
matplotlib.use("Agg")
from matplotlib import pyplot
from matplotlib import rcParams
rcParams["font.family"] = "serif"
rcParams["font.size"] = 20
pyplot.figure(figsize = (15,10))
pyplot.title("Convergence")
pyplot.plot(array_it, array_frame, lw=2)
pyplot.xlabel("iteration")
pyplot.ylabel("residual")
pyplot.grid()
pyplot.savefig("Convergence.pdf")