diff --git a/FDTD/processmodematch.cpp b/FDTD/processmodematch.cpp new file mode 100644 index 0000000..546fe40 --- /dev/null +++ b/FDTD/processmodematch.cpp @@ -0,0 +1,262 @@ +/* +* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + +#include "processmodematch.h" +#include "CSFunctionParser.h" +#include "tools/array_ops.h" + +ProcessModeMatch::ProcessModeMatch(Operator* op, Engine* eng) : ProcessIntegral(op, eng) +{ + for (int n=0;n<2;++n) + { + m_ModeParser[n] = new CSFunctionParser(); + m_ModeDist[n] = NULL; + } + m_dualMesh = false; +} + +ProcessModeMatch::~ProcessModeMatch() +{ + for (int n=0;n<2;++n) + { + delete m_ModeParser[n]; + m_ModeParser[n] = NULL; + } + Reset(); +} + +void ProcessModeMatch::InitProcess() +{ + if (!Enabled) return; + + if (m_ModeFieldType==1) + { + m_TimeShift = Op->GetTimestep()/2.0; + } + + int Dump_Dim=0; + for (int n=0;n<3;++n) + { + if (start[n]>stop[n]) + { + unsigned int help=start[n]; + start[n]=stop[n]; + stop[n]=help; + } + + if (stop[n]>start[n]) + ++Dump_Dim; + + if (stop[n] == start[n]) + m_ny = n; + } + + int nP = (m_ny+1)%3; + int nPP = (m_ny+2)%3; + m_numLines[0] = stop[nP] - start[nP] + 1; + m_numLines[1] = stop[nPP] - start[nPP] + 1; + + if (Dump_Dim!=2) + { + cerr << "ProcessModeMatch::InitProcess(): Warning Mode Matching Integration Box \"" << m_filename << "\" is not a surface (found dimension: " << Dump_Dim << ")" << endl; + SetEnable(false); + Reset(); + return; + } + + for (int n=0;n<2;++n) + { + int ny = (m_ny+n+1)%3; + int res = m_ModeParser[n]->Parse(m_ModeFunction[ny], "x,y,z,rho,a,r,t"); + if (res >= 0) + { + cerr << "ProcessModeMatch::InitProcess(): Warning, an error occured parsing the mode matching function (see below) ..." << endl; + cerr << m_ModeFunction[ny] << "\n" << string(res, ' ') << "^\n" << m_ModeParser[n]->ErrorMsg() << "\n"; + SetEnable(false); + Reset(); + } + } + + for (int n=0;n<2;++n) + { + m_ModeDist[n] = Create2DArray(m_numLines); + } + + unsigned int pos[3] = {0,0,0}; + double discLine[3] = {0,0,0}; + double gridDelta = Op->GetGridDelta(); + double var[7]; + pos[m_ny] = start[m_ny]; + discLine[m_ny] = Op->GetDiscLine(m_ny,pos[m_ny],m_dualMesh); + double norm = 0; + double area = 0; + for (unsigned int posP = 0;posPGetDiscLine(nP,pos[nP],m_dualMesh); + for (unsigned int posPP = 0;posPPGetDiscLine(nPP,pos[nPP],m_dualMesh); + + var[0] = discLine[0] * gridDelta; // x + var[1] = discLine[1] * gridDelta; // y + var[2] = discLine[2] * gridDelta; // z + var[3] = sqrt(discLine[0]*discLine[0] + discLine[1]*discLine[1]) * gridDelta; // rho = sqrt(x^2 + y^2) + var[4] = atan2(discLine[1], discLine[0]); // a = atan(y,x) + var[5] = sqrt(pow(discLine[0],2)+pow(discLine[1],2)+pow(discLine[2],2)) * gridDelta; // r + var[6] = asin(1)-atan(var[2]/var[3]); //theta (t) + + if (m_Mesh_Type == CYLINDRICAL_MESH) + { + var[3] = discLine[0] * gridDelta; // rho + var[4] = discLine[1]; // a + var[0] = discLine[0] * cos(discLine[1]) * gridDelta; // x = r*cos(a) + var[1] = discLine[0] * sin(discLine[1]) * gridDelta; // y = r*sin(a) + var[5] = sqrt(pow(discLine[0],2)+pow(discLine[2],2)) * gridDelta; // r + var[6] = asin(1)-atan(var[2]/var[3]); //theta (t) + } + area = Op->GetNodeArea(m_ny,pos,m_dualMesh); + for (int n=0;n<2;++n) + { + m_ModeDist[n][posP][posPP] = m_ModeParser[n]->Eval(var); //calc mode template + if ((isnan(m_ModeDist[n][posP][posPP])) || (isinf(m_ModeDist[n][posP][posPP]))) + m_ModeDist[n][posP][posPP] = 0.0; + norm += pow(m_ModeDist[n][posP][posPP],2) * area; + } +// cerr << discLine[0] << " " << discLine[1] << " : " << m_ModeDist[0][posP][posPP] << " , " << m_ModeDist[1][posP][posPP] << endl; + } + } + + norm = sqrt(norm); +// cerr << norm << endl; + // normalize template function... + for (unsigned int posP = 0;posP(m_ModeDist[n],m_numLines); + m_ModeDist[n] = NULL; + } +} + + +void ProcessModeMatch::SetModeFunction(int ny, string function) +{ + if ((ny<0) || (ny>2)) return; + m_ModeFunction[ny] = function; +} + +void ProcessModeMatch::SetFieldType(int type) +{ + m_ModeFieldType = type; + if ((type<0) || (type>1)) + cerr << "ProcessModeMatch::SetFieldType: Warning, unknown field type..." << endl; +} + +double ProcessModeMatch::GetField(int ny, unsigned int pos[3]) +{ + if (m_ModeFieldType==0) + return GetEField(ny,pos); + if (m_ModeFieldType==1) + return GetHField(ny,pos); + return 0; +} + +double ProcessModeMatch::GetEField(int ny, unsigned int pos[3]) +{ + if ((pos[ny]==0) || (pos[ny]==Op->GetNumberOfLines(ny)-1)) + return 0.0; + unsigned int DownPos[] = {pos[0],pos[1],pos[2]}; + --DownPos[ny]; + double delta = Op->GetMeshDelta(ny,pos); + double deltaDown = Op->GetMeshDelta(ny,DownPos); + double deltaRel = delta / (delta+deltaDown); + if (delta*deltaDown) + { + return (double)Eng->GetVolt(ny,pos)*(1.0-deltaRel)/delta + (double)Eng->GetVolt(ny,DownPos)/deltaDown*deltaRel; + } + return 0.0; +} + +double ProcessModeMatch::GetHField(int ny, unsigned int pos[3]) +{ + if ((pos[ny]==0) || (pos[ny]>=Op->GetNumberOfLines(ny)-1)) + return 0.0; + + unsigned int EngPos[] = {pos[0],pos[1],pos[2]}; + + int nyP = (ny+1)%3; + if (pos[nyP] >= Op->GetNumberOfLines(nyP)-1) + return 0.0; + int nyPP = (ny+2)%3; + if (pos[nyPP] >= Op->GetNumberOfLines(nyPP)-1) + return 0.0; + + double hfield = Eng->GetCurr(ny,EngPos) / Op->GetMeshDelta(ny,EngPos,true); + EngPos[nyP]++; + hfield += Eng->GetCurr(ny,EngPos) / Op->GetMeshDelta(ny,EngPos,true); + EngPos[nyPP]++; + hfield += Eng->GetCurr(ny,EngPos) / Op->GetMeshDelta(ny,EngPos,true); + EngPos[nyP]--; + hfield += Eng->GetCurr(ny,EngPos) / Op->GetMeshDelta(ny,EngPos,true); + return hfield/4.0; +} + + +double ProcessModeMatch::CalcIntegral() +{ + double value = 0; + double area = 0; + + int nP = (m_ny+1)%3; + int nPP = (m_ny+2)%3; + + unsigned int pos[3] = {0,0,0}; + pos[m_ny] = start[m_ny]; + + for (unsigned int posP = 0;posPGetNodeArea(m_ny,pos,m_dualMesh); + + for (int n=0;n<2;++n) + { + value += GetField((m_ny+n+1)%3,pos) * m_ModeDist[n][posP][posPP] * area; + } + } + } + + return value; +} diff --git a/FDTD/processmodematch.h b/FDTD/processmodematch.h new file mode 100644 index 0000000..d35e7be --- /dev/null +++ b/FDTD/processmodematch.h @@ -0,0 +1,55 @@ +/* +* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + +#ifndef PROCESSMODEMATCH_H +#define PROCESSMODEMATCH_H + +#include "processintegral.h" + +class CSFunctionParser; + +class ProcessModeMatch : public ProcessIntegral +{ +public: + ProcessModeMatch(Operator* op, Engine* eng); + virtual ~ProcessModeMatch(); + + virtual void InitProcess(); + virtual void Reset(); + + void SetFieldType(int type); + void SetModeFunction(int ny, string function); + virtual double CalcIntegral(); + +protected: + //normal direction of the mode plane + int m_ny; + + int m_ModeFieldType; + + double GetField(int ny, unsigned int pos[3]); + double GetEField(int ny, unsigned int pos[3]); + double GetHField(int ny, unsigned int pos[3]); + + string m_ModeFunction[3]; + CSFunctionParser* m_ModeParser[2]; + + unsigned int m_numLines[2]; + double** m_ModeDist[2]; +}; + +#endif // PROCESSMODEMATCH_H diff --git a/openEMS.pro b/openEMS.pro index cf45d0c..a219e72 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -51,9 +51,9 @@ SOURCES += main.cpp \ FDTD/operator.cpp \ tools/array_ops.cpp \ FDTD/processvoltage.cpp \ - FDTD/process_efield.cpp \ - FDTD/process_hfield.cpp \ - FDTD/processing.cpp \ + FDTD/process_efield.cpp \ + FDTD/process_hfield.cpp \ + FDTD/processing.cpp \ FDTD/processintegral.cpp \ FDTD/processfields.cpp \ FDTD/processfields_td.cpp \ @@ -80,7 +80,8 @@ SOURCES += main.cpp \ FDTD/engine_ext_dispersive.cpp \ FDTD/engine_ext_lorentzmaterial.cpp \ FDTD/operator_ext_pml_sf.cpp \ - FDTD/engine_ext_pml_sf.cpp + FDTD/engine_ext_pml_sf.cpp \ + FDTD/processmodematch.cpp HEADERS += tools/ErrorMsg.h \ tools/AdrOp.h \ tools/constants.h \ @@ -88,9 +89,9 @@ HEADERS += tools/ErrorMsg.h \ FDTD/operator.h \ tools/array_ops.h \ FDTD/processvoltage.h \ - FDTD/process_efield.h \ - FDTD/process_hfield.h \ - FDTD/processing.h \ + FDTD/process_efield.h \ + FDTD/process_hfield.h \ + FDTD/processing.h \ FDTD/processintegral.h \ FDTD/processfields.h \ FDTD/processfields_td.h \ @@ -117,7 +118,8 @@ HEADERS += tools/ErrorMsg.h \ FDTD/engine_ext_dispersive.h \ FDTD/engine_ext_lorentzmaterial.h \ FDTD/operator_ext_pml_sf.h \ - FDTD/engine_ext_pml_sf.h + FDTD/engine_ext_pml_sf.h \ + FDTD/processmodematch.h QMAKE_CXXFLAGS_RELEASE = -O3 \ -g \ -march=native diff --git a/openems.cpp b/openems.cpp index 50caa4e..01d5fe2 100644 --- a/openems.cpp +++ b/openems.cpp @@ -27,6 +27,7 @@ #include "FDTD/processcurrent.h" #include "FDTD/process_efield.h" #include "FDTD/process_hfield.h" +#include "FDTD/processmodematch.h" #include "FDTD/processfields_td.h" #include #include @@ -388,6 +389,15 @@ int openEMS::SetupFDTD(const char* file) proc = new ProcessEField(FDTD_Op,FDTD_Eng); else if (pb->GetProbeType()==3) proc = new ProcessHField(FDTD_Op,FDTD_Eng); + else if ((pb->GetProbeType()==10) || (pb->GetProbeType()==11)) + { + ProcessModeMatch* pmm = new ProcessModeMatch(FDTD_Op,FDTD_Eng); + pmm->SetFieldType(pb->GetProbeType()-10); + pmm->SetModeFunction(0,pb->GetAttributeValue("ModeFunctionX")); + pmm->SetModeFunction(1,pb->GetAttributeValue("ModeFunctionY")); + pmm->SetModeFunction(2,pb->GetAttributeValue("ModeFunctionZ")); + proc = pmm; + } else { cerr << "openEMS::SetupFDTD: Warning: Probe type " << pb->GetProbeType() << " of property '" << pb->GetName() << "' is unknown..." << endl;