Demo entry 2842032

blubb

   

Submitted by anonymous on Oct 10, 2015 at 15:39
Language: Python. Code size: 9.5 kB.

# -*- coding: utf-8 -*-
"""
Created on Thu Oct  8 15:50:16 2015

@author: user0
"""

""" 
Das Suite Modul übernimmt die Erstellung des Simulators sowie dessen Init-
ialisierung. Desweiteren stellt sie alle notwendigen Hilfsfunktionen zur 
Verfügung.
"""

"""
IMPORTS
"""
import numpy as np
import scipy.sparse as sp
import sys as sys
import timeit
import time
import matplotlib.pyplot as plt
from matplotlib import cm
from multiprocessing import Pool
from DataStructures import field_type
from Core import sim_4_freq
from Core import init_Core


class Simulator(object):

    def __init__(self,SIG,EPS,MU,hx,hy,fvec,SigY=None,SigZ=None,EpsY=None,EpsZ=None,\
    MuY=None,MuZ=None,which=None,neigs=20,nr_proc=4):
        '''SIG = Leitfaehigkeitsmatrix, EPS= Permit.matrix, MU=Permeab.matrix
        hxv/hyv = Abstaende der Simulationspunkte
        fvec = Frequenzvektor
        which = Definiert das Feld, dass simuliert werden soll ('H' oder 'E')
        neigs = Anzahl der zu berechnenden Eigenwerte
        Sig*,Eps*,Mu* dienen zur Implementierung von anisotropen Material. Default: wenn nicht
        uebergeben, erhalten sie die Werte von SIG,EPS bzw. MU
       '''
        tic_init = timeit.default_timer()
       #UEBERPRUEFUNG DES INPUTS
        SigY,SigZ,EpsY,EpsZ,\
        MuY,MuZ = input_Check(SIG,EPS,MU,SigY,SigZ,EpsY,EpsZ,MuY,MuZ,hx,hy,fvec,which,neigs)
        self.fvec = fvec
        self.field = which
        self.neigs = neigs
        self.nx = len(hx)
        print self.nx
        self.ny = len(hy)
        print self.ny
        self.n = self.nx * self.ny
        self.sigx,self.sigy,self.sigz = grid_2xTo1x(SIG,SigY,SigZ,'E')
        self.epsrx,self.epsry,self.epsrz = grid_2xTo1x(EPS,EpsY,EpsZ,'E')
        self.murx,self.mury,self.murz = grid_2xTo1x(MU,MuY,MuZ,'H')
        self.hx = hx[1:]
        self.hy = hy[1:]
        self.DEx,self.DEy = Dr_2dim(self.hx,self.hy)
        self.DHx,self.DHy = Dl_2dim(self.hx,self.hy)
        self.Xplot,self.Yplot = plot_mesh(hx,hy)
        self.number_of_processes = nr_proc
        self.simlist = []
        
        print 'Initialisierung abgeschlossen! Dauer {}s:'.format(timeit.default_timer()-tic_init)
        
        
    def simulate(self,freqs=None):
        plt.close('all')
        tic_sim = timeit.default_timer()
        init_Core(self)
        pool = Pool(processes=self.number_of_processes)
        if freqs == None:
            liste = pool.map(sim_4_freq,self.fvec)
        else:
            liste = pool.map(sim_4_freq,freqs)
        pool.close()
        pool.join()
        if freqs == None:
            self.simlist = liste
            print 'Simulation von {} Frequenzen in {:.4f} s abgeschlossen'.format(len(self.fvec),timeit.default_timer()-tic_sim)

        else:
            self.simlist.append(liste)
            self.fvec.append(freqs)
            print 'Simulation von {} Frequenzen in {:.4f} s abgeschlossen'.format(len(freqs),timeit.default_timer()-tic_sim)
     

    def getSimData(self,k=None):
            '''getSimData(k=None) gibt eine Liste von SimPoint-Strukturen zurueck. die Liste ist nach den Frequenzen des
            fvec Vektors geordnet. Bspw entspricht das Simulationsergebnis fuer fvec[3] gleich dem 4 Element der Liste (Beachte 0 basierte Indices)
            wird k uebergeben, werden die Simulationsdaten des k+1'ten Frequenzpunkts zurueckgegeben (Beachte auch hier 0 basierte Indices)'''
            if len(self.simlist)==0:
                sys.stderr.write('Die Simulation wurde noch nicht durchgefuehrt! Simulator.simulate() ausfuehren\n')
                return -1
            if k == None:
                return self.simlist
            else:
                if 0<k<len(self.simlist)-1:
                    return self.simlist[k]
                else:
                    sys.stderr.write('Index out of bounds\n')
                    return -1
    def plot_gamma(self,return_vectors=False):
        '''Plotten der Eigenwerte (alle Moden) über den simulierten Frequenzen.
            Bei Uebergabe des return_vector Flags (True) wird der geplottete Vector (Frequenz und Eigenwerte) zurückgegeben
            und kein eigener Plot angezeigt'''
        gam = []
        for k in self.simlist:
            gam.append(k.gamma)

        values = []
        for l in np.arange(0,len(self.fvec)):
            values = np.r_[values,gam[l][:]]
        freq_dim = np.kron(self.fvec,np.ones(self.neigs))
        if return_vectors:
            return values, freq_dim 

        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.plot(freq_dim,np.abs(np.real(values)) , ' bo',linewidth=1.1)
        ax.plot(freq_dim,np.abs(np.imag(values)) , ' rD',linewidth=1.0)
        fig.show()

    def plot_Fieldmag(self,k,field=None):
        '''urface-Plot der Feldkomponenten (Fx,Fy) des gewählten Frequenzpunktes aus simlist[k] (beachte: 0 indexing)
        field gibt an ob 'H' oder 'E' Feld geplottet werden soll. Default ist das Plotten beider Felder (-> 4 Plots)
        '''
        fig = plt.figure()
        simK =self.simlist[k]
        Emag = np.sqrt(abs(simK.Ex)**2 +abs(simK.Ey)**2)
        Emag = np.reshape(Emag,(self.ny,self.nx))
        Hmag = np.sqrt(abs(simK.Hx)**2 +abs(simK.Hy)**2)
        Hmag = np.reshape(Hmag,(self.ny,self.nx))
        colormap = cm.rainbow
        if field == 'E':
            ax = fig.add_subplot(111)
            plt.contourf(self.Xplot,self.Yplot,Emag.T,cmap=colormap)
            plt.colorbar()
            ax.set_title('Emag')
        elif field == 'H':
            ax = fig.add_subplot(111)
            plt.contourf(self.Xplot,self.Yplot,Hmag.T,cmap=colormap)
            plt.colorbar()
            ax.set_title('Hmag')

        else:
            ax = fig.add_subplot(211)
            plt.contourf(self.Xplot.T,self.Yplot.T,Emag,cmap=colormap)
            plt.colorbar()
            ax.set_title('Emag')
            ax = fig.add_subplot(212)
            plt.contourf(self.Xplot.T,self.Yplot.T,Hmag,cmap=colormap)
            plt.colorbar()
            ax.set_title('Hmag')
        fig.show()
       
def input_Check(SIG,EPS,MU,SigY,SigZ,EpsY,EpsZ,MuY,MuZ,hxv,hyv,fvec,which,neigs):
        if SigY == None:
            SigY = SIG
        if SigZ == None:
            SigZ = SIG
        if EpsY == None:
            EpsY = EPS
        if EpsZ == None:
            EpsZ = EPS
        if MuY == None:
            MuY = MU
        if MuZ == None:
            MuZ = MU

        xsig,ysig = np.shape(SIG)
        xsigY,ysigY = np.shape(SigY)
        xsigZ,ysigZ = np.shape(SigZ)
        xeps,yeps = np.shape(EPS)
        xepsY,yepsY = np.shape(EpsY)
        xepsZ,yepsZ = np.shape(EpsZ)
        xmu, ymu = np.shape(MU)
        xmuY, ymuY = np.shape(MuY)
        xmuZ, ymuZ = np.shape(MuZ)
        if not (xsig == xeps == xmu == xsigY == xepsY == xmuY == xsigZ == xepsZ == xmuZ):
            sys.stderr.write('X dimension mismatch!\n')
            time.sleep(3)
            sys.exit()
        if not (ysig == yeps == ymu == ysigY == yepsY == ymuY == ysigZ == yepsZ == ymuZ):
            sys.stderr.write('Y dimension mismatch!\n')
            time.sleep(3)
            sys.exit()
        if len(fvec) == 0:
            sys.stderr.write('Given empty frequency vector!\n')
            time.sleep(3)
            sys.exit()
        if which not in field_type:
            sys.stderr.write('Unknown fieldtype: which={E,H,None}!\n')
            time.sleep(3)
            sys.exit()
        return SigY,SigZ,EpsY,EpsZ,MuY,MuZ
def grid_2xTo1x(matx,maty,matz,K):
    '''Greift aus den Strukturmatrizen die Feldzuordnung heraus'''
    if K == 'E':
        matx = matx[1::2,::2].flatten('F')
        maty = maty[::2,1::2].flatten('F')
        matz = matz[::2,::2].flatten('F')
        return matx,maty,matz
    elif K == 'H':
        matx = matx[::2,1::2].flatten('F')
        maty = maty[1::2,::2].flatten('F')
        matz = matz[1::2,1::2].flatten('F')
        return matx,maty,matz
    else:
        sys.stderr.write('ERROR: K must either be \'E\' or \'H\' indicating assignment of material to E or H field\n\n')
        return 0
        
"""
ABLEITUNGSOPERATOREN
"""
#Ableitungsoperatoren
def Dr_1dim(h):
    n = np.size(h)+1
    h0 = np.r_[h,h[-1]]
    Dr = sp.diags([-h0**-1,h**-1],[0,1],shape=(n,n),format="csc")
    return Dr

def Dl_1dim(h):
    n = np.size(h)+1
    h0 = np.r_[h[0],h]
    Dl = sp.diags([h0**-1,-h**-1],[0,-1],shape=(n,n),format="csc")
    return Dl

def Dr_2dim(hx,hy):
    nx = np.size(hx)+1
    ny = np.size(hy)+1
    Ix = sp.eye(nx,nx,format="csc")
    Iy = sp.eye(ny,ny,format="csc")
    Drx_1d = Dr_1dim(hx)
    Dry_1d = Dr_1dim(hy)
    Drx_2d = sp.kron(Iy,Drx_1d,format="csc")
    Dry_2d = sp.kron(Dry_1d,Ix,format="csc")
    return Drx_2d,Dry_2d

def Dl_2dim(hx,hy):
    nx = np.size(hx)+1
    ny = np.size(hy)+1
    Ix = sp.eye(nx,nx,format="csc")
    Iy = sp.eye(ny,ny,format="csc")
    Dlx_1d = Dl_1dim(hx)
    Dly_1d = Dl_1dim(hy)
    Dlx_2d = sp.kron(Iy,Dlx_1d,format="csc")
    Dly_2d = sp.kron(Dly_1d,Ix,format="csc")
    return Dlx_2d,Dly_2d

def plot_mesh(hxv,hyv):
    x = np.cumsum(hxv)
    y = np.cumsum(hyv)
    X,Y = np.meshgrid(x,y,indexing='ij')
    return X,Y

def cmpComplex(a,b):
    '''Hilfsfunktion zum sortieren von Listen/Arrays mit komplexen Werten.
    a und b sind zwei Zahlen, deren Betrag (euklidsche Vektorlaenge) verglichen wird.
    Wenn gilt a=b: 0, a<b:-1 a>b:+1 zurueck gegeben'''

    if abs(a) == abs(b):
        return 0
    if abs(a) < abs(b):
        return -1
    if abs(a) > abs(b):
        return 1

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).