diff options
Diffstat (limited to 'openEMS/FDTD/extensions')
43 files changed, 6064 insertions, 0 deletions
diff --git a/openEMS/FDTD/extensions/CMakeLists.txt b/openEMS/FDTD/extensions/CMakeLists.txt new file mode 100644 index 0000000..3566944 --- /dev/null +++ b/openEMS/FDTD/extensions/CMakeLists.txt @@ -0,0 +1,32 @@ + +#INCLUDE_DIRECTORIES( ${openEMS_SOURCE_DIR} ) +#INCLUDE_DIRECTORIES( ${CSXCAD_SOURCE_DIR}/src ) + +set(SOURCES + ${SOURCES} + ${CMAKE_CURRENT_SOURCE_DIR}/engine_extension.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_dispersive.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_lorentzmaterial.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_conductingsheet.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_dispersive.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_lorentzmaterial.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_cylindermultigrid.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_upml.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_upml.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_extension.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_mur_abc.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_mur_abc.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_cylinder.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_cylinder.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_excitation.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_excitation.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_tfsf.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_tfsf.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_steadystate.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_steadystate.cpp + PARENT_SCOPE +) + +# FDTD/extensions lib +#add_library( extensions STATIC ${SOURCES} ) + diff --git a/openEMS/FDTD/extensions/OptimizeCondSheetParameter.m b/openEMS/FDTD/extensions/OptimizeCondSheetParameter.m new file mode 100644 index 0000000..12ac8ec --- /dev/null +++ b/openEMS/FDTD/extensions/OptimizeCondSheetParameter.m @@ -0,0 +1,107 @@ +function OptimizeCondSheetParameter +% function OptimizeCondSheetParameter +% +% internal openEMS function to create the "cond_sheet_parameter.h" header +% file containing optimized parameter values for the conducting sheet model +% +% (c) 2012: Thorsten Liebig, thorsten.liebig@gmx.de + + +close all +clear +clc + + +X = [1 1 1 1 1]; +lb = zeros(size(X)); +ub = ones(size(X))*1000; + + +Omega = linspace(0,20,51); +Omega = Omega(2:end); + +Omega_fine = linspace(0,20,15001); +Omega_fine = Omega_fine(2:end); + +options = optimset('fzero'); +% options = optimset(options,'Display','iter') +options = optimset(options,'MaxFunEvals',5000); +options = optimset(options,'MaxIter',5000); + +omg_stop_str = []; +omg_critical_20 = []; +omg_critical_5 = []; +g_str = []; +r1_str= []; +r2_str= []; +l1_str= []; +l2_str= []; + +numPara = 30; +factor = 1.5; + +for p = 1:numPara + + X = lsqnonlin(@(X)CalcYDiff(Omega,X),X,lb,ub,options); + + omg_stop_str = [omg_stop_str num2str(Omega(end),10) ',']; + g_str = [g_str num2str(X(1),10) ',']; + r1_str = [r1_str num2str(X(2),10) ',']; + l1_str = [l1_str num2str(X(3),10) ',']; + r2_str = [r2_str num2str(X(4),10) ',']; + l2_str = [l2_str num2str(X(5),10) ',']; + + Ys = tanh((1+1j)*sqrt(Omega_fine))/(1+1j)./sqrt(Omega_fine); + + err = (CalcYDiff(Omega_fine,X))./real(1./Ys); + + fc = Omega_fine(find(abs(err)>0.2,1,'last')); + if (isempty(fc)) + fc = 0; + end + omg_critical_20 = [omg_critical_20 num2str(fc,10) ',']; + + fc = Omega_fine(find(abs(err)>0.05,1,'last')); + if (isempty(fc)) + fc = 0; + end + omg_critical_5 = [omg_critical_5 num2str(fc,10) ',']; + +% disp(['max error: ' num2str(max(abs(err)*100)) ]) +% figure +% plot(Omega_fine,err*100,'g--'); + + Omega_fine = Omega_fine*factor; + Omega = Omega*factor; +end + +%% write to file +fid = fopen('cond_sheet_parameter.h','w'); + +fprintf(fid,'// This is a list of conducting sheet model parameter for different ranges of omega = w/w0\n'); +fprintf(fid,'// This file was created automatically using Matlab: OptimizeCondSheetParameter.m \n'); +fprintf(fid,'// Do not change this file! \n'); +fprintf(fid,'// Creation: %s \n\n',datestr(now)); + + + +fprintf(fid,'unsigned int numOptPara=%d;\n',numPara); + +fprintf(fid,'double omega_stop[%d]={%s};\n',numPara,omg_stop_str(1:end-1)); +fprintf(fid,'double omega_critical_5[%d]={%s};\n',numPara,omg_critical_5(1:end-1)); +fprintf(fid,'double omega_critical_20[%d]={%s};\n',numPara,omg_critical_20(1:end-1)); +fprintf(fid,'double g[%d]={%s};\n',numPara,g_str(1:end-1)); +fprintf(fid,'double r1[%d]={%s};\n',numPara,r1_str(1:end-1)); +fprintf(fid,'double l1[%d]={%s};\n',numPara,l1_str(1:end-1)); +fprintf(fid,'double r2[%d]={%s};\n',numPara,r2_str(1:end-1)); +fprintf(fid,'double l2[%d]={%s};\n',numPara,l2_str(1:end-1)); + +fclose(fid); + +function Ydiff = CalcYDiff(omega, X) + +Ys = tanh((1+1j)*sqrt(omega))/(1+1j)./sqrt(omega); + +Y = X(1) + 1./(X(2)+1j*omega*X(3)) + 1./(X(4)+1j*omega*X(5)); + +Ydiff = real(1./Ys)-real(1./Y);
\ No newline at end of file diff --git a/openEMS/FDTD/extensions/cond_sheet_parameter.h b/openEMS/FDTD/extensions/cond_sheet_parameter.h new file mode 100644 index 0000000..b3fc77b --- /dev/null +++ b/openEMS/FDTD/extensions/cond_sheet_parameter.h @@ -0,0 +1,14 @@ +// This is a list of conducting sheet model parameter for different ranges of omega = w/w0 +// This file was created automatically using Matlab: OptimizeCondSheetParameter.m +// Do not change this file! +// Creation: 08-May-2012 09:49:53 + +unsigned int numOptPara=30; +double omega_stop[30]={20,30,45,67.5,101.25,151.875,227.8125,341.71875,512.578125,768.8671875,1153.300781,1729.951172,2594.926758,3892.390137,5838.585205,8757.877808,13136.81671,19705.22507,29557.8376,44336.7564,66505.1346,99757.7019,149636.5529,224454.8293,336682.2439,505023.3659,757535.0488,1136302.573,1704453.86,2556680.79}; +double omega_critical_5[30]={0,0,0,1.566,6.8445,10.4085,15.1419375,22.30284375,33.45426562,50.28391406,75.42587109,113.1388066,169.70821,254.5623149,381.8434724,572.7652086,859.1478129,1288.721719,1933.082579,2899.623869,4349.435803,6524.153704,9786.230557,14679.34583,22019.01875,33028.52813,35503.14262,119084.5097,180444.8486,704280.3349}; +double omega_critical_20[30]={0,0,0,0,1.41075,2.7135,4.100625,5.90034375,8.6113125,12.86571094,19.37545312,29.06317969,43.59476953,65.3921543,98.08823145,147.1323472,220.6985208,331.0477811,496.5716717,744.8575075,1117.286261,1675.929392,2513.894088,3770.841132,5656.261698,8484.392547,10504.48601,11363.02573,11135.76522,7499.596983}; +double g[30]={0.1370843401,0.1244327709,0.1087425925,0.09215450299,0.07639249525,0.06239206861,0.05068701937,0.04124161681,0.03364925947,0.0274825647,0.02244156564,0.01832317373,0.01496080786,0.01221545053,0.009973874666,0.008143635713,0.006649251918,0.005429092739,0.00443283702,0.003619397792,0.002955227653,0.002412935361,0.001970156167,0.001608629812,0.001313423397,0.001072422468,0.0009031383439,0.0007674973627,0.0006482579571,0.0005450699565}; +double r1[30]={1.302754546,1.2992612,1.320033314,1.41829213,1.674326257,2.138770494,2.769126112,3.498802047,4.311446869,5.266204065,6.443980095,7.893402541,9.667414736,11.84010502,14.50110421,17.76014698,21.75163911,26.64019322,32.6274156,39.96021835,48.94100192,59.94011376,73.41109983,89.90930992,110.1195785,134.8631576,149.2770459,153.5258961,148.9373215,119.3522908}; +double l1[30]={0.9542567116,0.9412051444,0.8999655165,0.8189903271,0.702028543,0.5718817624,0.4577262324,0.3689027719,0.3004106661,0.2455514861,0.2005633225,0.1637496639,0.1337010005,0.1091664663,0.0891340699,0.07277768839,0.05942276018,0.04851850861,0.03961522618,0.03234573021,0.02641021614,0.02156389712,0.0176069076,0.01437606898,0.0117376269,0.009584103516,0.008750848841,0.008547706786,0.00854127488,0.008801399365}; +double r2[30]={10.62245057,10.21952962,10.24755832,10.99147713,12.70450822,15.59830718,19.55714986,24.32580904,29.88733852,36.55680626,44.75374078,54.81575133,67.13533424,82.22360504,100.7029064,123.3353165,151.0542147,185.0027444,226.5809553,277.5035137,339.8704104,416.2534638,509.8022059,624.3729288,764.7279941,936.5525496,1000,999.9999981,999.999984,999.9999985}; +double l2[30]={0.7925524826,0.6926195458,0.5814718257,0.4784644339,0.3916667501,0.3204661897,0.2620413998,0.2140939121,0.1748272154,0.1427375388,0.1165424396,0.09515681664,0.07769522035,0.06343788358,0.05179681795,0.04229192785,0.03453121766,0.02819462457,0.0230208182,0.01879642336,0.01534721999,0.01253095815,0.01023149158,0.008353988465,0.006820957973,0.005569332359,0.004553963932,0.003696380632,0.002978504727,0.002391295771}; diff --git a/openEMS/FDTD/extensions/engine_ext_cylinder.cpp b/openEMS/FDTD/extensions/engine_ext_cylinder.cpp new file mode 100644 index 0000000..e240f1b --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_cylinder.cpp @@ -0,0 +1,101 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "FDTD/engine.h" +#include "engine_ext_cylinder.h" +#include "operator_ext_cylinder.h" +#include "FDTD/engine_sse.h" + +Engine_Ext_Cylinder::Engine_Ext_Cylinder(Operator_Ext_Cylinder* op_ext) : Engine_Extension(op_ext) +{ + cyl_Op = op_ext; + m_Eng_SSE = NULL; + + CC_closedAlpha = op_ext->CC_closedAlpha; + CC_R0_included = op_ext->CC_R0_included; + + for (int n=0; n<3; ++n) + numLines[n] = op_ext->m_Op->GetNumberOfLines(n,true); + + //this cylindrical extension should be executed first? + m_Priority = ENG_EXT_PRIO_CYLINDER; +} + +void Engine_Ext_Cylinder::SetEngine(Engine* eng) +{ + Engine_Extension::SetEngine(eng); + m_Eng_SSE = dynamic_cast<Engine_sse*>(m_Eng); +} + +void Engine_Ext_Cylinder::DoPostVoltageUpdates() +{ + if (CC_closedAlpha==false) return; + + if (CC_R0_included) + { + unsigned int pos[3]; + pos[0] = 0; + FDTD_FLOAT volt=0; + for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) + { + volt =m_Eng_SSE->Engine_sse::GetVolt(2,0,0,pos[2])*cyl_Op->vv_R0[pos[2]]; + for (pos[1]=0; pos[1]<numLines[1]-1; ++pos[1]) + volt +=cyl_Op->vi_R0[pos[2]] * m_Eng_SSE->Engine_sse::GetCurr(1,0,pos[1],pos[2]); + m_Eng_SSE->Engine_sse::SetVolt(2,0,0,pos[2], volt); + } + + for (pos[1]=0; pos[1]<numLines[1]; ++pos[1]) + { + for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) + { + m_Eng_SSE->Engine_sse::SetVolt(1,0,pos[1],pos[2], 0); //no voltage in alpha-direction at r=0 + m_Eng_SSE->Engine_sse::SetVolt(2,0,pos[1],pos[2], m_Eng_SSE->Engine_sse::GetVolt(2,0,0,pos[2]) ); + } + } + } + + //close alpha + unsigned int pos[3]; + // copy tangential voltages from last alpha-plane to first + unsigned int last_A_Line = numLines[1]-2; + for (pos[0]=0; pos[0]<numLines[0]; ++pos[0]) + { + for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) + { + m_Eng_SSE->Engine_sse::SetVolt(0,pos[0],0,pos[2], m_Eng_SSE->Engine_sse::GetVolt(0,pos[0],last_A_Line,pos[2]) ); + m_Eng_SSE->Engine_sse::SetVolt(2,pos[0],0,pos[2], m_Eng_SSE->Engine_sse::GetVolt(2,pos[0],last_A_Line,pos[2]) ); + } + } +} + +void Engine_Ext_Cylinder::DoPostCurrentUpdates() +{ + if (CC_closedAlpha==false) return; + + //close alpha + unsigned int pos[3]; + // copy tangential currents from first alpha-plane to last + for (pos[0]=0; pos[0]<numLines[0]-1; ++pos[0]) + { + unsigned int last_A_Line = numLines[1]-2; + for (pos[2]=0; pos[2]<numLines[2]-1; ++pos[2]) + { + m_Eng_SSE->Engine_sse::SetCurr(0,pos[0],last_A_Line,pos[2], m_Eng_SSE->Engine_sse::GetCurr(0,pos[0],0,pos[2]) ); + m_Eng_SSE->Engine_sse::SetCurr(2,pos[0],last_A_Line,pos[2], m_Eng_SSE->Engine_sse::GetCurr(2,pos[0],0,pos[2]) ); + } + } +} diff --git a/openEMS/FDTD/extensions/engine_ext_cylinder.h b/openEMS/FDTD/extensions/engine_ext_cylinder.h new file mode 100644 index 0000000..c996a9d --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_cylinder.h @@ -0,0 +1,49 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef ENGINE_CYLINDER_H +#define ENGINE_CYLINDER_H + +#include "FDTD/engine.h" +#include "engine_extension.h" +#include "FDTD/operator_cylinder.h" + +class Operator_Ext_Cylinder; +class Engine_sse; + +class Engine_Ext_Cylinder : public Engine_Extension +{ +public: + Engine_Ext_Cylinder(Operator_Ext_Cylinder* op_ext); + + virtual void DoPostVoltageUpdates(); + + virtual void DoPostCurrentUpdates(); + + virtual void SetEngine(Engine* eng); + +protected: + Operator_Ext_Cylinder* cyl_Op; + Engine_sse* m_Eng_SSE; + + unsigned int numLines[3]; + + bool CC_closedAlpha; + bool CC_R0_included; +}; + +#endif // ENGINE_CYLINDER_H diff --git a/openEMS/FDTD/extensions/engine_ext_cylindermultigrid.cpp b/openEMS/FDTD/extensions/engine_ext_cylindermultigrid.cpp new file mode 100644 index 0000000..eca2c9e --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_cylindermultigrid.cpp @@ -0,0 +1,165 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "FDTD/engine.h" +#include "engine_ext_cylindermultigrid.h" +#include "FDTD/engine_cylindermultigrid.h" + +Engine_Ext_CylinderMultiGrid::Engine_Ext_CylinderMultiGrid(Operator_Extension* op_ext, bool isBase) : Engine_Extension(op_ext) +{ + m_IsBase = isBase; + m_Eng_MG = NULL; + + // the multi-grid should be applies last? + m_Priority = ENG_EXT_PRIO_CYLINDERMULTIGRID; +} + +Engine_Ext_CylinderMultiGrid::~Engine_Ext_CylinderMultiGrid() +{ +} + +void Engine_Ext_CylinderMultiGrid::SetBarrier(boost::barrier* waitBase, boost::barrier* waitChild, boost::barrier* waitSync) +{ + m_WaitOnBase = waitBase; + m_WaitOnChild = waitChild; + m_WaitOnSync = waitSync; +} + +void Engine_Ext_CylinderMultiGrid::SetEngine(Engine* eng) +{ + m_Eng_MG = dynamic_cast<Engine_CylinderMultiGrid*>(eng); + if (m_Eng_MG==NULL) + { + cerr << "Engine_Ext_CylinderMultiGrid::SetEngine(): Error" << endl; + exit(0); + } +} + +void Engine_Ext_CylinderMultiGrid::DoPreVoltageUpdates() +{ + //cerr << "Engine_Ext_CylinderMultiGrid::DoPreVoltageUpdates() for " << m_IsBase << endl; + if (!m_IsBase) + { + //cerr << "child: volt wait on base " << endl; + m_WaitOnBase->wait(); //wait on base to finisch current sync and/or to finisch voltage updates, than start child voltage updates + } +} + +void Engine_Ext_CylinderMultiGrid::DoPostVoltageUpdates() +{ + +} + +void Engine_Ext_CylinderMultiGrid::Apply2Voltages() +{ + if (m_IsBase) + { + m_WaitOnBase->wait(); //base voltage updates are done, tell child to start its voltage updates + m_WaitOnChild->wait(); //wait for child to finisch its updates + SyncVoltages(); //child is finisch, run sync and go to current updates next + m_WaitOnSync->wait(); //sync is done... move on and tell child to move on... + } + else + { + m_WaitOnChild->wait(); //child is finished voltage updates, will tell base to run sync + m_WaitOnSync->wait(); //wait for base to finisch sync before going to wait for current updates + } +} + +void Engine_Ext_CylinderMultiGrid::SyncVoltages() +{ + if (m_Eng_MG==NULL) + { + cerr << "Engine_Ext_CylinderMultiGrid::SyncVoltages: Error engine is NULL" << endl; + return; + } + + unsigned int* numLines = m_Eng_MG->numLines; + + Engine_Multithread* m_InnerEng = m_Eng_MG->m_InnerEngine; + + //interpolate voltages from base engine to child engine... + unsigned int pos[3]; + pos[0] = m_Eng_MG->Op_CMG->GetSplitPos()-1; + unsigned int pos1_half = 0; + f4vector v_null; + v_null.f[0] = 0; + v_null.f[1] = 0; + v_null.f[2] = 0; + v_null.f[3] = 0; + for (pos[1]=0; pos[1]<numLines[1]-1; pos[1]+=2) + { + pos1_half = pos[1]/2; + for (pos[2]=0; pos[2]<m_Eng_MG->numVectors; ++pos[2]) + { + //r - direczion + m_InnerEng->f4_volt[0][pos[0]][pos1_half][pos[2]].v = v_null.v; + + //z - direction + m_InnerEng->f4_volt[2][pos[0]][pos1_half][pos[2]].v = m_Eng_MG->f4_volt[2][pos[0]][pos[1]][pos[2]].v; + + //alpha - direction + m_InnerEng->f4_volt[1][pos[0]][pos1_half][pos[2]].v = m_Eng_MG->f4_volt[1][pos[0]][pos[1]][pos[2]].v; + m_InnerEng->f4_volt[1][pos[0]][pos1_half][pos[2]].v += m_Eng_MG->f4_volt[1][pos[0]][pos[1]+1][pos[2]].v; + } + } +} + +void Engine_Ext_CylinderMultiGrid::DoPreCurrentUpdates() +{ + //cerr << "Engine_Ext_CylinderMultiGrid::DoPreCurrentUpdates() for " << m_IsBase << endl; + if (!m_IsBase) + { + //cerr << "child: curr wait on base " << endl; + m_WaitOnBase->wait(); //wait on base to finisch voltage sync and current updates, than start child current updates + } +} + +void Engine_Ext_CylinderMultiGrid::DoPostCurrentUpdates() +{ + +} + +void Engine_Ext_CylinderMultiGrid::Apply2Current() +{ + if (m_IsBase) + { + //cerr << "Base: curr wait on base done, wait on sync" << endl; + m_WaitOnBase->wait(); //base current updates are done, tell child to start its current updates + m_WaitOnChild->wait(); //wait for child to finisch its updates + SyncCurrents(); //child is finisch, run sync and go to voltage updates next + m_WaitOnSync->wait(); //sync is done... move on and tell child to move on... + } + else + { + m_WaitOnChild->wait(); //child is finished current updates, will tell base to run sync... + m_WaitOnSync->wait(); //wait for base to finisch sync before going to wait for next voltage updates + //cerr << "Child: curr done, wait on sync" << endl; + } +} + +void Engine_Ext_CylinderMultiGrid::SyncCurrents() +{ + if (m_Eng_MG==NULL) + { + cerr << "Engine_Ext_CylinderMultiGrid::SyncCurrents: Error engine is NULL" << endl; + return; + } + + m_Eng_MG->InterpolCurrChild2Base(m_Eng_MG->Op_CMG->GetSplitPos()-2); + return; +} diff --git a/openEMS/FDTD/extensions/engine_ext_cylindermultigrid.h b/openEMS/FDTD/extensions/engine_ext_cylindermultigrid.h new file mode 100644 index 0000000..e9226f5 --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_cylindermultigrid.h @@ -0,0 +1,58 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef ENGINE_EXT_CYLINDERMULTIGRID_H +#define ENGINE_EXT_CYLINDERMULTIGRID_H + +#include "FDTD/engine_cylindermultigrid.h" +#include "engine_extension.h" +#include "FDTD/operator_cylindermultigrid.h" + +class Operator_Ext_CylinderMultiGrid; + +class Engine_Ext_CylinderMultiGrid : public Engine_Extension +{ +public: + Engine_Ext_CylinderMultiGrid(Operator_Extension* op_ext, bool isBase); + virtual ~Engine_Ext_CylinderMultiGrid(); + + void SetBarrier(boost::barrier* waitBase, boost::barrier* waitChild, boost::barrier* waitSync); + + virtual void DoPreVoltageUpdates(); + virtual void DoPostVoltageUpdates(); + virtual void Apply2Voltages(); + + virtual void DoPreCurrentUpdates(); + virtual void DoPostCurrentUpdates(); + virtual void Apply2Current(); + + virtual void SetEngine(Engine* eng); + +protected: + void SyncVoltages(); + void SyncCurrents(); + + Engine_CylinderMultiGrid* m_Eng_MG; + + boost::barrier *m_WaitOnBase; + boost::barrier *m_WaitOnChild; + boost::barrier *m_WaitOnSync; + + bool m_IsBase; +}; + +#endif // ENGINE_EXT_CYLINDERMULTIGRID_H diff --git a/openEMS/FDTD/extensions/engine_ext_dispersive.cpp b/openEMS/FDTD/extensions/engine_ext_dispersive.cpp new file mode 100644 index 0000000..d9c431b --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_dispersive.cpp @@ -0,0 +1,162 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "engine_ext_dispersive.h" +#include "operator_ext_dispersive.h" +#include "FDTD/engine_sse.h" + +Engine_Ext_Dispersive::Engine_Ext_Dispersive(Operator_Ext_Dispersive* op_ext_disp) : Engine_Extension(op_ext_disp) +{ + m_Op_Ext_Disp = op_ext_disp; + int order = m_Op_Ext_Disp->m_Order; + curr_ADE = new FDTD_FLOAT**[order]; + volt_ADE = new FDTD_FLOAT**[order]; + for (int o=0;o<order;++o) + { + curr_ADE[o] = new FDTD_FLOAT*[3]; + volt_ADE[o] = new FDTD_FLOAT*[3]; + for (int n=0; n<3; ++n) + { + if (m_Op_Ext_Disp->m_curr_ADE_On[o]==true) + { + curr_ADE[o][n] = new FDTD_FLOAT[m_Op_Ext_Disp->m_LM_Count[o]]; + for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count[o]; ++i) + curr_ADE[o][n][i]=0.0; + } + else + curr_ADE[o][n] = NULL; + if (m_Op_Ext_Disp->m_volt_ADE_On[o]==true) + { + volt_ADE[o][n] = new FDTD_FLOAT[m_Op_Ext_Disp->m_LM_Count[o]]; + for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count[o]; ++i) + volt_ADE[o][n][i]=0.0; + } + else + volt_ADE[o][n] = NULL; + } + } +} + +Engine_Ext_Dispersive::~Engine_Ext_Dispersive() +{ + if (curr_ADE==NULL && volt_ADE==NULL) + return; + + for (int o=0;o<m_Op_Ext_Disp->m_Order;++o) + { + for (int n=0; n<3; ++n) + { + delete[] curr_ADE[o][n]; + delete[] volt_ADE[o][n]; + } + delete[] curr_ADE[o]; + delete[] volt_ADE[o]; + } + delete[] curr_ADE; + curr_ADE=NULL; + + delete[] volt_ADE; + volt_ADE=NULL; +} + +void Engine_Ext_Dispersive::Apply2Voltages() +{ + for (int o=0;o<m_Op_Ext_Disp->m_Order;++o) + { + if (m_Op_Ext_Disp->m_volt_ADE_On[o]==false) continue; + + unsigned int **pos = m_Op_Ext_Disp->m_LM_pos[o]; + + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i) + { + m_Eng->Engine::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][0][i]); + m_Eng->Engine::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][1][i]); + m_Eng->Engine::SetVolt(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][2][i]); + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i) + { + eng_sse->Engine_sse::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][0][i]); + eng_sse->Engine_sse::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][1][i]); + eng_sse->Engine_sse::SetVolt(2,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][2][i]); + } + break; + } + default: + for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i) + { + m_Eng->SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][0][i]); + m_Eng->SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][1][i]); + m_Eng->SetVolt(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][2][i]); + } + break; + } + } +} + +void Engine_Ext_Dispersive::Apply2Current() +{ + for (int o=0;o<m_Op_Ext_Disp->m_Order;++o) + { + if (m_Op_Ext_Disp->m_curr_ADE_On[o]==false) continue; + + unsigned int **pos = m_Op_Ext_Disp->m_LM_pos[o]; + + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i) + { + m_Eng->Engine::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][0][i]); + m_Eng->Engine::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][1][i]); + m_Eng->Engine::SetCurr(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][2][i]); + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i) + { + eng_sse->Engine_sse::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][0][i]); + eng_sse->Engine_sse::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][1][i]); + eng_sse->Engine_sse::SetCurr(2,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][2][i]); + } + break; + } + default: + for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i) + { + m_Eng->SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][0][i]); + m_Eng->SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][1][i]); + m_Eng->SetCurr(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][2][i]); + } + break; + } + } +} diff --git a/openEMS/FDTD/extensions/engine_ext_dispersive.h b/openEMS/FDTD/extensions/engine_ext_dispersive.h new file mode 100644 index 0000000..c1b84ea --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_dispersive.h @@ -0,0 +1,51 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef ENGINE_EXT_DISPERSIVE_H +#define ENGINE_EXT_DISPERSIVE_H + +#include "engine_extension.h" +#include "FDTD/engine.h" +#include "FDTD/operator.h" + +class Operator_Ext_Dispersive; + +class Engine_Ext_Dispersive : public Engine_Extension +{ +public: + Engine_Ext_Dispersive(Operator_Ext_Dispersive* op_ext_disp); + virtual ~Engine_Ext_Dispersive(); + + virtual void Apply2Voltages(); + virtual void Apply2Current(); + +protected: + Operator_Ext_Dispersive* m_Op_Ext_Disp; + + //! Dispersive order + int m_Order; + + //! ADE currents + // Array setup: curr_ADE[N_order][direction][mesh_pos] + FDTD_FLOAT ***curr_ADE; + + //! ADE voltages + // Array setup: volt_ADE[N_order][direction][mesh_pos] + FDTD_FLOAT ***volt_ADE; +}; + +#endif // ENGINE_EXT_DISPERSIVE_H diff --git a/openEMS/FDTD/extensions/engine_ext_excitation.cpp b/openEMS/FDTD/extensions/engine_ext_excitation.cpp new file mode 100644 index 0000000..b5bd8bd --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_excitation.cpp @@ -0,0 +1,170 @@ +/* +* Copyright (C) 2011 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 <http://www.gnu.org/licenses/>. +*/ + +#include "engine_ext_excitation.h" +#include "operator_ext_excitation.h" +#include "FDTD/engine_sse.h" + +Engine_Ext_Excitation::Engine_Ext_Excitation(Operator_Ext_Excitation* op_ext) : Engine_Extension(op_ext) +{ + m_Op_Exc = op_ext; + m_Priority = ENG_EXT_PRIO_EXCITATION; +} + +Engine_Ext_Excitation::~Engine_Ext_Excitation() +{ + +} + +void Engine_Ext_Excitation::Apply2Voltages() +{ + //soft voltage excitation here (E-field excite) + int exc_pos; + unsigned int ny; + unsigned int pos[3]; + int numTS = m_Eng->GetNumberOfTimesteps(); + unsigned int length = m_Op_Exc->m_Exc->GetLength(); + FDTD_FLOAT* exc_volt = m_Op_Exc->m_Exc->GetVoltageSignal(); + + int p = numTS+1; + if (m_Op_Exc->m_Exc->GetSignalPeriod()>0) + p = int(m_Op_Exc->m_Exc->GetSignalPeriod()/m_Op_Exc->m_Exc->GetTimestep()); + + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int n=0; n<m_Op_Exc->Volt_Count; ++n) + { + exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n]; + exc_pos *= (exc_pos>0); + exc_pos %= p; + exc_pos *= (exc_pos<=(int)length); + ny = m_Op_Exc->Volt_dir[n]; + pos[0]=m_Op_Exc->Volt_index[0][n]; + pos[1]=m_Op_Exc->Volt_index[1][n]; + pos[2]=m_Op_Exc->Volt_index[2][n]; + m_Eng->Engine::SetVolt(ny,pos, m_Eng->Engine::GetVolt(ny,pos) + m_Op_Exc->Volt_amp[n]*exc_volt[exc_pos]); + } + break; + } + case Engine::SSE: + { + for (unsigned int n=0; n<m_Op_Exc->Volt_Count; ++n) + { + Engine_sse* eng_sse = (Engine_sse*) m_Eng; + exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n]; + exc_pos *= (exc_pos>0); + exc_pos %= p; + exc_pos *= (exc_pos<=(int)length); + ny = m_Op_Exc->Volt_dir[n]; + pos[0]=m_Op_Exc->Volt_index[0][n]; + pos[1]=m_Op_Exc->Volt_index[1][n]; + pos[2]=m_Op_Exc->Volt_index[2][n]; + eng_sse->Engine_sse::SetVolt(ny,pos, eng_sse->Engine_sse::GetVolt(ny,pos) + m_Op_Exc->Volt_amp[n]*exc_volt[exc_pos]); + } + break; + } + default: + { + for (unsigned int n=0; n<m_Op_Exc->Volt_Count; ++n) + { + exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n]; + exc_pos *= (exc_pos>0); + exc_pos %= p; + exc_pos *= (exc_pos<=(int)length); + ny = m_Op_Exc->Volt_dir[n]; + pos[0]=m_Op_Exc->Volt_index[0][n]; + pos[1]=m_Op_Exc->Volt_index[1][n]; + pos[2]=m_Op_Exc->Volt_index[2][n]; + m_Eng->SetVolt(ny,pos, m_Eng->GetVolt(ny,pos) + m_Op_Exc->Volt_amp[n]*exc_volt[exc_pos]); + } + break; + } + } +} + +void Engine_Ext_Excitation::Apply2Current() +{ + //soft current excitation here (H-field excite) + + int exc_pos; + unsigned int ny; + unsigned int pos[3]; + int numTS = m_Eng->GetNumberOfTimesteps(); + unsigned int length = m_Op_Exc->m_Exc->GetLength(); + FDTD_FLOAT* exc_curr = m_Op_Exc->m_Exc->GetCurrentSignal(); + + int p = numTS+1; + if (m_Op_Exc->m_Exc->GetSignalPeriod()>0) + p = int(m_Op_Exc->m_Exc->GetSignalPeriod()/m_Op_Exc->m_Exc->GetTimestep()); + + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int n=0; n<m_Op_Exc->Curr_Count; ++n) + { + exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n]; + exc_pos *= (exc_pos>0); + exc_pos %= p; + exc_pos *= (exc_pos<=(int)length); + ny = m_Op_Exc->Curr_dir[n]; + pos[0]=m_Op_Exc->Curr_index[0][n]; + pos[1]=m_Op_Exc->Curr_index[1][n]; + pos[2]=m_Op_Exc->Curr_index[2][n]; + m_Eng->Engine::SetCurr(ny,pos, m_Eng->Engine::GetCurr(ny,pos) + m_Op_Exc->Curr_amp[n]*exc_curr[exc_pos]); + } + break; + } + case Engine::SSE: + { + for (unsigned int n=0; n<m_Op_Exc->Curr_Count; ++n) + { + Engine_sse* eng_sse = (Engine_sse*) m_Eng; + exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n]; + exc_pos *= (exc_pos>0); + exc_pos %= p; + exc_pos *= (exc_pos<=(int)length); + ny = m_Op_Exc->Curr_dir[n]; + pos[0]=m_Op_Exc->Curr_index[0][n]; + pos[1]=m_Op_Exc->Curr_index[1][n]; + pos[2]=m_Op_Exc->Curr_index[2][n]; + eng_sse->Engine_sse::SetCurr(ny,pos, eng_sse->Engine_sse::GetCurr(ny,pos) + m_Op_Exc->Curr_amp[n]*exc_curr[exc_pos]); + } + break; + } + default: + { + for (unsigned int n=0; n<m_Op_Exc->Curr_Count; ++n) + { + exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n]; + exc_pos *= (exc_pos>0); + exc_pos %= p; + exc_pos *= (exc_pos<=(int)length); + ny = m_Op_Exc->Curr_dir[n]; + pos[0]=m_Op_Exc->Curr_index[0][n]; + pos[1]=m_Op_Exc->Curr_index[1][n]; + pos[2]=m_Op_Exc->Curr_index[2][n]; + m_Eng->SetCurr(ny,pos, m_Eng->GetCurr(ny,pos) + m_Op_Exc->Curr_amp[n]*exc_curr[exc_pos]); + } + break; + } + } +} diff --git a/openEMS/FDTD/extensions/engine_ext_excitation.h b/openEMS/FDTD/extensions/engine_ext_excitation.h new file mode 100644 index 0000000..aae58a0 --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_excitation.h @@ -0,0 +1,40 @@ +/* +* Copyright (C) 2011 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef ENGINE_EXT_EXCITATION_H +#define ENGINE_EXT_EXCITATION_H + +#include "engine_extension.h" +#include "FDTD/engine.h" +#include "FDTD/operator.h" + +class Operator_Ext_Excitation; + +class Engine_Ext_Excitation : public Engine_Extension +{ +public: + Engine_Ext_Excitation(Operator_Ext_Excitation* op_ext); + virtual ~Engine_Ext_Excitation(); + + virtual void Apply2Voltages(); + virtual void Apply2Current(); + +protected: + Operator_Ext_Excitation* m_Op_Exc; +}; + +#endif // ENGINE_EXT_EXCITATION_H diff --git a/openEMS/FDTD/extensions/engine_ext_lorentzmaterial.cpp b/openEMS/FDTD/extensions/engine_ext_lorentzmaterial.cpp new file mode 100644 index 0000000..9aaa1ff --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_lorentzmaterial.cpp @@ -0,0 +1,322 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "engine_ext_lorentzmaterial.h" +#include "operator_ext_lorentzmaterial.h" +#include "FDTD/engine_sse.h" + +Engine_Ext_LorentzMaterial::Engine_Ext_LorentzMaterial(Operator_Ext_LorentzMaterial* op_ext_lorentz) : Engine_Ext_Dispersive(op_ext_lorentz) +{ + m_Op_Ext_Lor = op_ext_lorentz; + m_Order = m_Op_Ext_Lor->GetDispersionOrder(); + int order = m_Op_Ext_Lor->m_Order; + + curr_Lor_ADE = new FDTD_FLOAT**[order]; + volt_Lor_ADE = new FDTD_FLOAT**[order]; + for (int o=0;o<order;++o) + { + curr_Lor_ADE[o] = new FDTD_FLOAT*[3]; + volt_Lor_ADE[o] = new FDTD_FLOAT*[3]; + for (int n=0; n<3; ++n) + { + if (m_Op_Ext_Lor->m_curr_Lor_ADE_On[o]==true) + { + curr_Lor_ADE[o][n] = new FDTD_FLOAT[m_Op_Ext_Lor->m_LM_Count[o]]; + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count[o]; ++i) + curr_Lor_ADE[o][n][i]=0.0; + } + else + curr_Lor_ADE[o][n] = NULL; + + if (m_Op_Ext_Lor->m_volt_Lor_ADE_On[o]==true) + { + volt_Lor_ADE[o][n] = new FDTD_FLOAT[m_Op_Ext_Lor->m_LM_Count[o]]; + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count[o]; ++i) + volt_Lor_ADE[o][n][i]=0.0; + } + else + volt_Lor_ADE[o][n] = NULL; + } + } +} + +Engine_Ext_LorentzMaterial::~Engine_Ext_LorentzMaterial() +{ + if (curr_Lor_ADE==NULL && volt_Lor_ADE==NULL) + return; + + for (int o=0;o<m_Op_Ext_Lor->m_Order;++o) + { + for (int n=0; n<3; ++n) + { + delete[] curr_Lor_ADE[o][n]; + delete[] volt_Lor_ADE[o][n]; + } + delete[] curr_Lor_ADE[o]; + delete[] volt_Lor_ADE[o]; + } + delete[] curr_Lor_ADE; + curr_Lor_ADE=NULL; + + delete[] volt_Lor_ADE; + volt_Lor_ADE=NULL; +} + +void Engine_Ext_LorentzMaterial::DoPreVoltageUpdates() +{ + for (int o=0;o<m_Order;++o) + { + if (m_Op_Ext_Lor->m_volt_ADE_On[o]==false) continue; + + unsigned int **pos = m_Op_Ext_Lor->m_LM_pos[o]; + + if (m_Op_Ext_Lor->m_volt_Lor_ADE_On[o]) + { + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + volt_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][0][i]*volt_ADE[o][0][i]; + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * (m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][0][i]); + + volt_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][1][i]*volt_ADE[o][1][i]; + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * (m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][2][i]); + + volt_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][2][i]*volt_ADE[o][2][i]; + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * (m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][2][i]); + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + volt_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][0][i]*volt_ADE[o][0][i]; + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * (eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][0][i]); + + volt_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][1][i]*volt_ADE[o][1][i]; + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * (eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][1][i]); + + volt_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][2][i]*volt_ADE[o][2][i]; + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * (eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][2][i]); + } + break; + } + default: + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + volt_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][0][i]*volt_ADE[o][0][i]; + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * (m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][0][i]); + + volt_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][1][i]*volt_ADE[o][1][i]; + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * (m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][1][i]); + + volt_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][2][i]*volt_ADE[o][2][i]; + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * (m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][2][i]); + } + break; + } + } + else + { + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); + + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]); + + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); + + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]); + + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; + } + default: + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); + + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]); + + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; + } + } + } +} + +void Engine_Ext_LorentzMaterial::DoPreCurrentUpdates() +{ + for (int o=0;o<m_Order;++o) + { + if (m_Op_Ext_Lor->m_curr_ADE_On[o]==false) continue; + + unsigned int **pos = m_Op_Ext_Lor->m_LM_pos[o]; + + if (m_Op_Ext_Lor->m_curr_Lor_ADE_On[o]) + { + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + curr_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][0][i]*curr_ADE[o][0][i]; + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * (m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][0][i]); + + curr_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][1][i]*curr_ADE[o][1][i]; + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * (m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][1][i]); + + curr_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][2][i]*curr_ADE[o][2][i]; + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * (m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][2][i]); + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + curr_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][0][i]*curr_ADE[o][0][i]; + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * (eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][0][i]); + + curr_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][1][i]*curr_ADE[o][1][i]; + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * (eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][1][i]); + + curr_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][2][i]*curr_ADE[o][2][i]; + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * (eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][2][i]); + } + break; + } + default: + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + curr_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][0][i]*curr_ADE[o][0][i]; + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * (m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][0][i]); + + curr_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][1][i]*curr_ADE[o][1][i]; + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * (m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][1][i]); + + curr_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][2][i]*curr_ADE[o][2][i]; + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * (m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][2][i]); + } + break; + } + } + else + { + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); + + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]); + + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); + + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]); + + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; + } + default: + for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i) + { + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); + + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]); + + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; + } + } + } +} + diff --git a/openEMS/FDTD/extensions/engine_ext_lorentzmaterial.h b/openEMS/FDTD/extensions/engine_ext_lorentzmaterial.h new file mode 100644 index 0000000..804f49a --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_lorentzmaterial.h @@ -0,0 +1,48 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef ENGINE_EXT_LORENTZMATERIAL_H +#define ENGINE_EXT_LORENTZMATERIAL_H + +#include "engine_ext_dispersive.h" + +class Operator_Ext_LorentzMaterial; + +class Engine_Ext_LorentzMaterial : public Engine_Ext_Dispersive +{ +public: + Engine_Ext_LorentzMaterial(Operator_Ext_LorentzMaterial* op_ext_lorentz); + virtual ~Engine_Ext_LorentzMaterial(); + + virtual void DoPreVoltageUpdates(); + + virtual void DoPreCurrentUpdates(); + +protected: + Operator_Ext_LorentzMaterial* m_Op_Ext_Lor; + + //! ADE Lorentz voltages + // Array setup: volt_Lor_ADE[N_order][direction][mesh_pos] + FDTD_FLOAT ***volt_Lor_ADE; + + //! ADE Lorentz currents + // Array setup: curr_Lor_ADE[N_order][direction][mesh_pos] + FDTD_FLOAT ***curr_Lor_ADE; + +}; + +#endif // ENGINE_EXT_LORENTZMATERIAL_H diff --git a/openEMS/FDTD/extensions/engine_ext_mur_abc.cpp b/openEMS/FDTD/extensions/engine_ext_mur_abc.cpp new file mode 100644 index 0000000..d9f1b17 --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_mur_abc.cpp @@ -0,0 +1,263 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "engine_ext_mur_abc.h" +#include "operator_ext_mur_abc.h" +#include "FDTD/engine.h" +#include "FDTD/engine_sse.h" +#include "tools/array_ops.h" +#include "tools/useful.h" +#include "operator_ext_excitation.h" + +Engine_Ext_Mur_ABC::Engine_Ext_Mur_ABC(Operator_Ext_Mur_ABC* op_ext) : Engine_Extension(op_ext) +{ + m_Op_mur = op_ext; + m_numLines[0] = m_Op_mur->m_numLines[0]; + m_numLines[1] = m_Op_mur->m_numLines[1]; + m_ny = m_Op_mur->m_ny; + m_nyP = m_Op_mur->m_nyP; + m_nyPP = m_Op_mur->m_nyPP; + m_LineNr = m_Op_mur->m_LineNr; + m_LineNr_Shift = m_Op_mur->m_LineNr_Shift; + + m_Mur_Coeff_nyP = m_Op_mur->m_Mur_Coeff_nyP; + m_Mur_Coeff_nyPP = m_Op_mur->m_Mur_Coeff_nyPP; + + m_volt_nyP = Create2DArray<FDTD_FLOAT>(m_numLines); + m_volt_nyPP = Create2DArray<FDTD_FLOAT>(m_numLines); + + //find if some excitation is on this mur-abc and find the max length of this excite, so that the abc can start after the excitation is done... + int maxDelay=-1; + Operator_Ext_Excitation* Exc_ext = m_Op_mur->m_Op->GetExcitationExtension(); + for (unsigned int n=0; n<Exc_ext->GetVoltCount(); ++n) + { + if ( ((Exc_ext->Volt_dir[n]==m_nyP) || (Exc_ext->Volt_dir[n]==m_nyPP)) && (Exc_ext->Volt_index[m_ny][n]==m_LineNr) ) + { + if ((int)Exc_ext->Volt_delay[n]>maxDelay) + maxDelay = (int)Exc_ext->Volt_delay[n]; + } + } + m_start_TS = 0; + if (maxDelay>=0) + { + m_start_TS = maxDelay + m_Op_mur->m_Op->GetExcitationSignal()->GetLength() + 10; //give it some extra timesteps, for the excitation to travel at least one cell away + cerr << "Engine_Ext_Mur_ABC::Engine_Ext_Mur_ABC: Warning: Excitation inside the Mur-ABC #" << m_ny << "-" << (int)(m_LineNr>0) << " found!!!! Mur-ABC will be switched on after excitation is done at " << m_start_TS << " timesteps!!! " << endl; + } + + SetNumberOfThreads(1); +} + +Engine_Ext_Mur_ABC::~Engine_Ext_Mur_ABC() +{ + Delete2DArray(m_volt_nyP,m_numLines); + m_volt_nyP = NULL; + Delete2DArray(m_volt_nyPP,m_numLines); + m_volt_nyPP = NULL; +} + + +void Engine_Ext_Mur_ABC::SetNumberOfThreads(int nrThread) +{ + Engine_Extension::SetNumberOfThreads(nrThread); + + m_numX = AssignJobs2Threads(m_numLines[0],m_NrThreads,false); + m_start.resize(m_NrThreads,0); + m_start.at(0)=0; + for (size_t n=1; n<m_numX.size(); ++n) + m_start.at(n) = m_start.at(n-1) + m_numX.at(n-1); +} + + +void Engine_Ext_Mur_ABC::DoPreVoltageUpdates(int threadID) +{ + if (IsActive()==false) return; + if (m_Eng==NULL) return; + if (threadID>=m_NrThreads) + return; + unsigned int pos[] = {0,0,0}; + unsigned int pos_shift[] = {0,0,0}; + pos[m_ny] = m_LineNr; + pos_shift[m_ny] = m_LineNr_Shift; + + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + pos[m_nyP]=lineX+m_start.at(threadID); + pos_shift[m_nyP] = pos[m_nyP]; + for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) + { + pos_shift[m_nyPP] = pos[m_nyPP]; + m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->Engine::GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->Engine::GetVolt(m_nyP,pos); + m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->Engine::GetVolt(m_nyPP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->Engine::GetVolt(m_nyPP,pos); + } + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*) m_Eng; + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + pos[m_nyP]=lineX+m_start.at(threadID); + pos_shift[m_nyP] = pos[m_nyP]; + for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) + { + pos_shift[m_nyPP] = pos[m_nyPP]; + m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = eng_sse->Engine_sse::GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * eng_sse->Engine_sse::GetVolt(m_nyP,pos); + m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] = eng_sse->Engine_sse::GetVolt(m_nyPP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] * eng_sse->Engine_sse::GetVolt(m_nyPP,pos); + } + } + break; + } + default: + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + pos[m_nyP]=lineX+m_start.at(threadID); + pos_shift[m_nyP] = pos[m_nyP]; + for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) + { + pos_shift[m_nyPP] = pos[m_nyPP]; + m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyP,pos); + m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyPP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyPP,pos); + } + } + break; + } +} + +void Engine_Ext_Mur_ABC::DoPostVoltageUpdates(int threadID) +{ + if (IsActive()==false) return; + if (m_Eng==NULL) return; + if (threadID>=m_NrThreads) + return; + unsigned int pos[] = {0,0,0}; + unsigned int pos_shift[] = {0,0,0}; + pos[m_ny] = m_LineNr; + pos_shift[m_ny] = m_LineNr_Shift; + + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + pos[m_nyP]=lineX+m_start.at(threadID); + pos_shift[m_nyP] = pos[m_nyP]; + for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) + { + pos_shift[m_nyPP] = pos[m_nyPP]; + m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->Engine::GetVolt(m_nyP,pos_shift); + m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->Engine::GetVolt(m_nyPP,pos_shift); + } + } + break; + } + + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*) m_Eng; + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + pos[m_nyP]=lineX+m_start.at(threadID); + pos_shift[m_nyP] = pos[m_nyP]; + for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) + { + pos_shift[m_nyPP] = pos[m_nyPP]; + m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * eng_sse->Engine_sse::GetVolt(m_nyP,pos_shift); + m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] * eng_sse->Engine_sse::GetVolt(m_nyPP,pos_shift); + } + } + break; + } + + default: + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + pos[m_nyP]=lineX+m_start.at(threadID); + pos_shift[m_nyP] = pos[m_nyP]; + for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) + { + pos_shift[m_nyPP] = pos[m_nyPP]; + m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyP,pos_shift); + m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyPP,pos_shift); + } + } + break; + } +} + +void Engine_Ext_Mur_ABC::Apply2Voltages(int threadID) +{ + if (IsActive()==false) return; + if (threadID>=m_NrThreads) + return; + if (m_Eng==NULL) return; + unsigned int pos[] = {0,0,0}; + pos[m_ny] = m_LineNr; + + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + pos[m_nyP]=lineX+m_start.at(threadID); + for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) + { + m_Eng->Engine::SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]); + m_Eng->Engine::SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]); + } + } + break; + } + + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*) m_Eng; + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + pos[m_nyP]=lineX+m_start.at(threadID); + for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) + { + eng_sse->Engine_sse::SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]); + eng_sse->Engine_sse::SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]); + } + } + break; + } + + default: + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + pos[m_nyP]=lineX+m_start.at(threadID); + for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) + { + m_Eng->SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]); + m_Eng->SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]); + } + } + break; + } + +} diff --git a/openEMS/FDTD/extensions/engine_ext_mur_abc.h b/openEMS/FDTD/extensions/engine_ext_mur_abc.h new file mode 100644 index 0000000..aadb2fd --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_mur_abc.h @@ -0,0 +1,63 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef ENGINE_EXT_MUR_ABC_H +#define ENGINE_EXT_MUR_ABC_H + +#include "engine_extension.h" +#include "FDTD/engine.h" +#include "FDTD/operator.h" + +class Operator_Ext_Mur_ABC; + +class Engine_Ext_Mur_ABC : public Engine_Extension +{ +public: + Engine_Ext_Mur_ABC(Operator_Ext_Mur_ABC* op_ext); + virtual ~Engine_Ext_Mur_ABC(); + + virtual void SetNumberOfThreads(int nrThread); + + virtual void DoPreVoltageUpdates() {Engine_Ext_Mur_ABC::DoPreVoltageUpdates(0);} + virtual void DoPreVoltageUpdates(int threadID); + virtual void DoPostVoltageUpdates() {Engine_Ext_Mur_ABC::DoPostVoltageUpdates(0);} + virtual void DoPostVoltageUpdates(int threadID); + virtual void Apply2Voltages() {Engine_Ext_Mur_ABC::Apply2Voltages(0);} + virtual void Apply2Voltages(int threadID); + +protected: + Operator_Ext_Mur_ABC* m_Op_mur; + + inline bool IsActive() {if (m_Eng->GetNumberOfTimesteps()<m_start_TS) return false; return true;} + unsigned int m_start_TS; + + int m_ny; + int m_nyP,m_nyPP; + unsigned int m_LineNr; + int m_LineNr_Shift; + unsigned int m_numLines[2]; + + vector<unsigned int> m_start; + vector<unsigned int> m_numX; + + FDTD_FLOAT** m_Mur_Coeff_nyP; + FDTD_FLOAT** m_Mur_Coeff_nyPP; + FDTD_FLOAT** m_volt_nyP; //n+1 direction + FDTD_FLOAT** m_volt_nyPP; //n+2 direction +}; + +#endif // ENGINE_EXT_MUR_ABC_H diff --git a/openEMS/FDTD/extensions/engine_ext_steadystate.cpp b/openEMS/FDTD/extensions/engine_ext_steadystate.cpp new file mode 100644 index 0000000..9e5d0ce --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_steadystate.cpp @@ -0,0 +1,112 @@ +/* +* Copyright (C) 2015 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 <http://www.gnu.org/licenses/>. +*/ + +#include "engine_ext_steadystate.h" +#include "operator_ext_steadystate.h" +#include "FDTD/engine_sse.h" +#include "FDTD/engine_interface_fdtd.h" + +Engine_Ext_SteadyState::Engine_Ext_SteadyState(Operator_Ext_SteadyState* op_ext): Engine_Extension(op_ext) +{ + m_Op_SS = op_ext; + m_Priority = ENG_EXT_PRIO_STEADYSTATE; + + for (size_t n=0;n<m_Op_SS->m_E_probe_dir.size();++n) + { + double* rec = new double[m_Op_SS->m_TS_period*2]; + m_E_records.push_back(rec); + } + m_last_max_diff = 1; + last_total_energy = 0; + m_Eng_Interface = NULL; +} + +Engine_Ext_SteadyState::~Engine_Ext_SteadyState() +{ + for (size_t n=0;n<m_E_records.size();++n) + { + delete[] m_E_records.at(n); + m_E_records.at(n) = NULL; + } + m_E_records.clear(); + delete m_Eng_Interface; + m_Eng_Interface = NULL; +} + +void Engine_Ext_SteadyState::Apply2Voltages() +{ + unsigned int p = m_Op_SS->m_TS_period; + unsigned int TS = m_Eng->GetNumberOfTimesteps(); + unsigned int rel_pos = m_Eng->GetNumberOfTimesteps()%(2*p); + for (size_t n=0;n<m_E_records.size();++n) + m_E_records.at(n)[rel_pos] = m_Eng->GetVolt(m_Op_SS->m_E_probe_dir.at(n), m_Op_SS->m_E_probe_pos[0].at(n), m_Op_SS->m_E_probe_pos[1].at(n), m_Op_SS->m_E_probe_pos[2].at(n)); + if ((TS%(m_Op_SS->m_TS_period)==0) && (TS>=2*p)) + { + bool no_valid = true; + m_last_max_diff = 0; + double curr_total_energy = m_Eng_Interface->CalcFastEnergy(); + if (last_total_energy>0) + { + m_last_max_diff = abs(curr_total_energy-last_total_energy)/last_total_energy; + no_valid = false; + } + //cerr << curr_total_energy << "/" << last_total_energy << "=" << abs(curr_total_energy-last_total_energy)/last_total_energy << endl; + last_total_energy = curr_total_energy; + unsigned int old_pos = 0; + unsigned int new_pos = p; + if (rel_pos<=p) + { + new_pos = 0; + old_pos = p; + } + //cerr << TS << "/" << rel_pos << ": one period complete, new_pos" << new_pos << " old pos: " << old_pos << endl; + double *curr_pow = new double[m_E_records.size()]; + double *diff_pow = new double[m_E_records.size()]; + double max_pow = 0; + for (size_t n=0;n<m_E_records.size();++n) + { + double *buf = m_E_records.at(n); + curr_pow[n] = 0; + diff_pow[n] = 0; + for (unsigned int nt=0;nt<p;++nt) + { + curr_pow[n] += buf[nt+new_pos]*buf[nt+new_pos]; + diff_pow[n] += (buf[nt+old_pos]-buf[nt+new_pos])*(buf[nt+old_pos]-buf[nt+new_pos]); + } + max_pow = max(max_pow, curr_pow[n]); + } + for (size_t n=0;n<m_E_records.size();++n) + { + //cerr << "curr_pow: " << curr_pow[n] << " diff_pow: " << diff_pow[n] << " diff: " << diff_pow[n]/curr_pow[n] << endl; + if (curr_pow[n]>max_pow*1e-2) + { + m_last_max_diff = max(m_last_max_diff, diff_pow[n]/curr_pow[n]); + //cerr << m_last_max_diff << endl; + no_valid = false; + } + } + if ((no_valid) || (m_last_max_diff>1)) + m_last_max_diff = 1; + delete[] curr_pow; curr_pow = NULL; + //cerr << m_last_max_diff << endl; + } +} + +void Engine_Ext_SteadyState::Apply2Current() +{ + +} diff --git a/openEMS/FDTD/extensions/engine_ext_steadystate.h b/openEMS/FDTD/extensions/engine_ext_steadystate.h new file mode 100644 index 0000000..b66596e --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_steadystate.h @@ -0,0 +1,51 @@ +/* +* Copyright (C) 2015 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef ENGINE_EXT_STEADYSTATE_H +#define ENGINE_EXT_STEADYSTATE_H + +#include "engine_extension.h" +#include "FDTD/engine.h" +#include "FDTD/operator.h" + +class Operator_Ext_SteadyState; +class Engine_Interface_FDTD; + +class Engine_Ext_SteadyState : public Engine_Extension +{ +public: + Engine_Ext_SteadyState(Operator_Ext_SteadyState* op_ext); + virtual ~Engine_Ext_SteadyState(); + + virtual void Apply2Voltages(); + virtual void Apply2Current(); + + void SetEngineInterface(Engine_Interface_FDTD* eng_if) {m_Eng_Interface=eng_if;} + double GetLastDiff() {return m_last_max_diff;} + +protected: + Operator_Ext_SteadyState* m_Op_SS; + double m_last_max_diff; + vector<double*> m_E_records; + vector<double*> m_H_records; + + double last_total_energy; + Engine_Interface_FDTD* m_Eng_Interface; +}; + + +#endif // ENGINE_EXT_STEADYSTATE_H diff --git a/openEMS/FDTD/extensions/engine_ext_tfsf.cpp b/openEMS/FDTD/extensions/engine_ext_tfsf.cpp new file mode 100644 index 0000000..cd3b107 --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_tfsf.cpp @@ -0,0 +1,215 @@ +/* +* Copyright (C) 2012 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 <http://www.gnu.org/licenses/>. +*/ + +#include "engine_ext_tfsf.h" +#include "operator_ext_tfsf.h" +#include "FDTD/engine_sse.h" + +Engine_Ext_TFSF::Engine_Ext_TFSF(Operator_Ext_TFSF* op_ext) : Engine_Extension(op_ext) +{ + m_Op_TFSF = op_ext; + m_Priority = ENG_EXT_PRIO_TFSF; + + m_DelayLookup = new unsigned int[m_Op_TFSF->m_maxDelay+1]; +} + +Engine_Ext_TFSF::~Engine_Ext_TFSF() +{ + delete[] m_DelayLookup; + m_DelayLookup = NULL; +} + +void Engine_Ext_TFSF::DoPostVoltageUpdates() +{ + unsigned int numTS = m_Eng->GetNumberOfTimesteps(); + unsigned int length = m_Op_TFSF->m_Exc->GetLength(); + + int p = int(m_Op_TFSF->m_Exc->GetSignalPeriod()/m_Op_TFSF->m_Exc->GetTimestep()); + + for (unsigned int n=0;n<=m_Op_TFSF->m_maxDelay;++n) + { + if ( numTS < n ) + m_DelayLookup[n]=0; + else if ((numTS-n > length) && (p==0)) + m_DelayLookup[n]=0; + else + m_DelayLookup[n] = numTS - n; + if (p>0) + m_DelayLookup[n] = (m_DelayLookup[n] % p); + } + + //get the current signal since an H-field is added ... + FDTD_FLOAT* signal = m_Op_TFSF->m_Exc->GetCurrentSignal(); + + int nP,nPP; + unsigned int ui_pos; + unsigned int pos[3]; + for (int n=0;n<3;++n) + { + nP = (n+1)%3; + nPP = (n+2)%3; + + // lower plane + pos[nP] = m_Op_TFSF->m_Start[nP]; + ui_pos = 0; + if (m_Op_TFSF->m_ActiveDir[n][0]) + { + + for (unsigned int i=0;i<m_Op_TFSF->m_numLines[nP];++i) + { + pos[nPP] = m_Op_TFSF->m_Start[nPP]; + for (unsigned int j=0;j<m_Op_TFSF->m_numLines[nPP];++j) + { + // current updates + pos[n] = m_Op_TFSF->m_Start[n]; + + m_Eng->SetVolt(nP,pos, m_Eng->GetVolt(nP,pos) + + (1.0-m_Op_TFSF->m_VoltDelayDelta[n][0][0][ui_pos])*m_Op_TFSF->m_VoltAmp[n][0][0][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_VoltDelay[n][0][0][ui_pos]]] + + m_Op_TFSF->m_VoltDelayDelta[n][0][0][ui_pos] *m_Op_TFSF->m_VoltAmp[n][0][0][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_VoltDelay[n][0][0][ui_pos]]] ); + + m_Eng->SetVolt(nPP,pos, m_Eng->GetVolt(nPP,pos) + + (1.0-m_Op_TFSF->m_VoltDelayDelta[n][0][1][ui_pos])*m_Op_TFSF->m_VoltAmp[n][0][1][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_VoltDelay[n][0][1][ui_pos]]] + + m_Op_TFSF->m_VoltDelayDelta[n][0][1][ui_pos] *m_Op_TFSF->m_VoltAmp[n][0][1][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_VoltDelay[n][0][1][ui_pos]]] ); + + ++pos[nPP]; + ++ui_pos; + } + ++pos[nP]; + } + } + + // upper plane + pos[nP] = m_Op_TFSF->m_Start[nP]; + ui_pos = 0; + if (m_Op_TFSF->m_ActiveDir[n][1]) + { + + for (unsigned int i=0;i<m_Op_TFSF->m_numLines[nP];++i) + { + pos[nPP] = m_Op_TFSF->m_Start[nPP]; + for (unsigned int j=0;j<m_Op_TFSF->m_numLines[nPP];++j) + { + // current updates + pos[n] = m_Op_TFSF->m_Stop[n]; + + m_Eng->SetVolt(nP,pos, m_Eng->GetVolt(nP,pos) + + (1.0-m_Op_TFSF->m_VoltDelayDelta[n][1][0][ui_pos])*m_Op_TFSF->m_VoltAmp[n][1][0][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_VoltDelay[n][1][0][ui_pos]]] + + m_Op_TFSF->m_VoltDelayDelta[n][1][0][ui_pos] *m_Op_TFSF->m_VoltAmp[n][1][0][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_VoltDelay[n][1][0][ui_pos]]] ); + + m_Eng->SetVolt(nPP,pos, m_Eng->GetVolt(nPP,pos) + + (1.0-m_Op_TFSF->m_VoltDelayDelta[n][1][1][ui_pos])*m_Op_TFSF->m_VoltAmp[n][1][1][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_VoltDelay[n][1][1][ui_pos]]] + + m_Op_TFSF->m_VoltDelayDelta[n][1][1][ui_pos] *m_Op_TFSF->m_VoltAmp[n][1][1][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_VoltDelay[n][1][1][ui_pos]]] ); + + ++pos[nPP]; + ++ui_pos; + } + ++pos[nP]; + } + } + } +} + +void Engine_Ext_TFSF::DoPostCurrentUpdates() +{ + unsigned int numTS = m_Eng->GetNumberOfTimesteps(); + unsigned int length = m_Op_TFSF->m_Exc->GetLength(); + + int p = int(m_Op_TFSF->m_Exc->GetSignalPeriod()/m_Op_TFSF->m_Exc->GetTimestep()); + + for (unsigned int n=0;n<m_Op_TFSF->m_maxDelay;++n) + { + if ( numTS < n ) + m_DelayLookup[n]=0; + else if ((numTS-n > length) && (p==0)) + m_DelayLookup[n]=0; + else + m_DelayLookup[n] = numTS - n; + if (p>0) + m_DelayLookup[n] = (m_DelayLookup[n] % p); + } + + //get the current signal since an E-field is added ... + FDTD_FLOAT* signal = m_Op_TFSF->m_Exc->GetVoltageSignal(); + + int nP,nPP; + unsigned int ui_pos; + unsigned int pos[3]; + for (int n=0;n<3;++n) + { + if (!m_Op_TFSF->m_ActiveDir[n][0] && !m_Op_TFSF->m_ActiveDir[n][1]) + continue; + + nP = (n+1)%3; + nPP = (n+2)%3; + + // lower plane + pos[nP] = m_Op_TFSF->m_Start[nP]; + ui_pos = 0; + if (m_Op_TFSF->m_ActiveDir[n][0]) + { + for (unsigned int i=0;i<m_Op_TFSF->m_numLines[nP];++i) + { + pos[nPP] = m_Op_TFSF->m_Start[nPP]; + for (unsigned int j=0;j<m_Op_TFSF->m_numLines[nPP];++j) + { + // current updates + pos[n] = m_Op_TFSF->m_Start[n]-1; + + m_Eng->SetCurr(nP,pos, m_Eng->GetCurr(nP,pos) + + (1.0-m_Op_TFSF->m_CurrDelayDelta[n][0][0][ui_pos])*m_Op_TFSF->m_CurrAmp[n][0][0][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_CurrDelay[n][0][0][ui_pos]]] + + m_Op_TFSF->m_CurrDelayDelta[n][0][0][ui_pos] *m_Op_TFSF->m_CurrAmp[n][0][0][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_CurrDelay[n][0][0][ui_pos]]] ); + + m_Eng->SetCurr(nPP,pos, m_Eng->GetCurr(nPP,pos) + + (1.0-m_Op_TFSF->m_CurrDelayDelta[n][0][1][ui_pos])*m_Op_TFSF->m_CurrAmp[n][0][1][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_CurrDelay[n][0][1][ui_pos]]] + + m_Op_TFSF->m_CurrDelayDelta[n][0][1][ui_pos] *m_Op_TFSF->m_CurrAmp[n][0][1][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_CurrDelay[n][0][1][ui_pos]]] ); + + ++pos[nPP]; + ++ui_pos; + } + ++pos[nP]; + } + } + + // upper plane + pos[nP] = m_Op_TFSF->m_Start[nP]; + ui_pos = 0; + if (m_Op_TFSF->m_ActiveDir[n][1]) + { + for (unsigned int i=0;i<m_Op_TFSF->m_numLines[nP];++i) + { + pos[nPP] = m_Op_TFSF->m_Start[nPP]; + for (unsigned int j=0;j<m_Op_TFSF->m_numLines[nPP];++j) + { + // current updates + pos[n] = m_Op_TFSF->m_Stop[n]; + + m_Eng->SetCurr(nP,pos, m_Eng->GetCurr(nP,pos) + + (1.0-m_Op_TFSF->m_CurrDelayDelta[n][1][0][ui_pos])*m_Op_TFSF->m_CurrAmp[n][1][0][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_CurrDelay[n][1][0][ui_pos]]] + + m_Op_TFSF->m_CurrDelayDelta[n][1][0][ui_pos] *m_Op_TFSF->m_CurrAmp[n][1][0][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_CurrDelay[n][1][0][ui_pos]]] ); + + m_Eng->SetCurr(nPP,pos, m_Eng->GetCurr(nPP,pos) + + (1.0-m_Op_TFSF->m_CurrDelayDelta[n][1][1][ui_pos])*m_Op_TFSF->m_CurrAmp[n][1][1][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_CurrDelay[n][1][1][ui_pos]]] + + m_Op_TFSF->m_CurrDelayDelta[n][1][1][ui_pos] *m_Op_TFSF->m_CurrAmp[n][1][1][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_CurrDelay[n][1][1][ui_pos]]] ); + + ++pos[nPP]; + ++ui_pos; + } + ++pos[nP]; + } + } + } +} diff --git a/openEMS/FDTD/extensions/engine_ext_tfsf.h b/openEMS/FDTD/extensions/engine_ext_tfsf.h new file mode 100644 index 0000000..d4e662a --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_tfsf.h @@ -0,0 +1,40 @@ +/* +* Copyright (C) 2012 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef ENGINE_EXT_TFSF_H +#define ENGINE_EXT_TFSF_H + +#include "engine_extension.h" + +class Operator_Ext_TFSF; + +class Engine_Ext_TFSF : public Engine_Extension +{ +public: + Engine_Ext_TFSF(Operator_Ext_TFSF* op_ext); + virtual ~Engine_Ext_TFSF(); + + virtual void DoPostVoltageUpdates(); + virtual void DoPostCurrentUpdates(); + +protected: + Operator_Ext_TFSF* m_Op_TFSF; + + unsigned int* m_DelayLookup; +}; + +#endif // ENGINE_EXT_TFSF_H diff --git a/openEMS/FDTD/extensions/engine_ext_upml.cpp b/openEMS/FDTD/extensions/engine_ext_upml.cpp new file mode 100644 index 0000000..8cb365f --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_upml.cpp @@ -0,0 +1,493 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "engine_ext_upml.h" +#include "operator_ext_upml.h" +#include "FDTD/engine.h" +#include "FDTD/engine_sse.h" +#include "tools/array_ops.h" +#include "tools/useful.h" + +Engine_Ext_UPML::Engine_Ext_UPML(Operator_Ext_UPML* op_ext) : Engine_Extension(op_ext) +{ + m_Op_UPML = op_ext; + + //this ABC extension should be executed first! + m_Priority = ENG_EXT_PRIO_UPML; + + volt_flux = Create_N_3DArray<FDTD_FLOAT>(m_Op_UPML->m_numLines); + curr_flux = Create_N_3DArray<FDTD_FLOAT>(m_Op_UPML->m_numLines); + + SetNumberOfThreads(1); +} + +Engine_Ext_UPML::~Engine_Ext_UPML() +{ + Delete_N_3DArray<FDTD_FLOAT>(volt_flux,m_Op_UPML->m_numLines); + volt_flux=NULL; + Delete_N_3DArray<FDTD_FLOAT>(curr_flux,m_Op_UPML->m_numLines); + curr_flux=NULL; +} + +void Engine_Ext_UPML::SetNumberOfThreads(int nrThread) +{ + Engine_Extension::SetNumberOfThreads(nrThread); + + m_numX = AssignJobs2Threads(m_Op_UPML->m_numLines[0],m_NrThreads,false); + m_start.resize(m_NrThreads,0); + m_start.at(0)=0; + for (size_t n=1; n<m_numX.size(); ++n) + m_start.at(n) = m_start.at(n-1) + m_numX.at(n-1); +} + + +void Engine_Ext_UPML::DoPreVoltageUpdates(int threadID) +{ + if (m_Eng==NULL) + return; + + if (threadID>=m_NrThreads) + return; + + unsigned int pos[3]; + unsigned int loc_pos[3]; + FDTD_FLOAT f_help; + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = m_Op_UPML->vv[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->Engine::GetVolt(0,pos) + - m_Op_UPML->vvfo[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->Engine::SetVolt(0,pos, volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->vv[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->Engine::GetVolt(1,pos) + - m_Op_UPML->vvfo[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->Engine::SetVolt(1,pos, volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->vv[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->Engine::GetVolt(2,pos) + - m_Op_UPML->vvfo[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->Engine::SetVolt(2,pos, volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + } + } + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*) m_Eng; + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = m_Op_UPML->vv[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * eng_sse->Engine_sse::GetVolt(0,pos) + - m_Op_UPML->vvfo[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + eng_sse->Engine_sse::SetVolt(0,pos, volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->vv[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * eng_sse->Engine_sse::GetVolt(1,pos) + - m_Op_UPML->vvfo[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + eng_sse->Engine_sse::SetVolt(1,pos, volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->vv[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * eng_sse->Engine_sse::GetVolt(2,pos) + - m_Op_UPML->vvfo[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + eng_sse->Engine_sse::SetVolt(2,pos, volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + } + } + } + break; + } + default: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = m_Op_UPML->vv[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetVolt(0,pos) + - m_Op_UPML->vvfo[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->SetVolt(0,pos, volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->vv[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetVolt(1,pos) + - m_Op_UPML->vvfo[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->SetVolt(1,pos, volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->vv[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetVolt(2,pos) + - m_Op_UPML->vvfo[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->SetVolt(2,pos, volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + } + } + } + break; + } + } + +} + +void Engine_Ext_UPML::DoPostVoltageUpdates(int threadID) +{ + if (m_Eng==NULL) + return; + if (threadID>=m_NrThreads) + return; + + unsigned int pos[3]; + unsigned int loc_pos[3]; + FDTD_FLOAT f_help; + + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->Engine::GetVolt(0,pos); + m_Eng->Engine::SetVolt(0,pos, f_help + m_Op_UPML->vvfn[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->Engine::GetVolt(1,pos); + m_Eng->Engine::SetVolt(1,pos, f_help + m_Op_UPML->vvfn[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->Engine::GetVolt(2,pos); + m_Eng->Engine::SetVolt(2,pos, f_help + m_Op_UPML->vvfn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + } + } + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*) m_Eng; + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = eng_sse->Engine_sse::GetVolt(0,pos); + eng_sse->Engine_sse::SetVolt(0,pos, f_help + m_Op_UPML->vvfn[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = eng_sse->Engine_sse::GetVolt(1,pos); + eng_sse->Engine_sse::SetVolt(1,pos, f_help + m_Op_UPML->vvfn[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = eng_sse->Engine_sse::GetVolt(2,pos); + eng_sse->Engine_sse::SetVolt(2,pos, f_help + m_Op_UPML->vvfn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + } + } + } + break; + } + default: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetVolt(0,pos); + m_Eng->SetVolt(0,pos, f_help + m_Op_UPML->vvfn[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetVolt(1,pos); + m_Eng->SetVolt(1,pos, f_help + m_Op_UPML->vvfn[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetVolt(2,pos); + m_Eng->SetVolt(2,pos, f_help + m_Op_UPML->vvfn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + } + } + } + break; + } + } + +} + +void Engine_Ext_UPML::DoPreCurrentUpdates(int threadID) +{ + if (m_Eng==NULL) + return; + if (threadID>=m_NrThreads) + return; + + unsigned int pos[3]; + unsigned int loc_pos[3]; + FDTD_FLOAT f_help; + + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = m_Op_UPML->ii[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->Engine::GetCurr(0,pos) + - m_Op_UPML->iifo[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->Engine::SetCurr(0,pos, curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->ii[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->Engine::GetCurr(1,pos) + - m_Op_UPML->iifo[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->Engine::SetCurr(1,pos, curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->ii[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->Engine::GetCurr(2,pos) + - m_Op_UPML->iifo[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->Engine::SetCurr(2,pos, curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + } + } + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*) m_Eng; + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = m_Op_UPML->ii[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * eng_sse->Engine_sse::GetCurr(0,pos) + - m_Op_UPML->iifo[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + eng_sse->Engine_sse::SetCurr(0,pos, curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->ii[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * eng_sse->Engine_sse::GetCurr(1,pos) + - m_Op_UPML->iifo[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + eng_sse->Engine_sse::SetCurr(1,pos, curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->ii[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * eng_sse->Engine_sse::GetCurr(2,pos) + - m_Op_UPML->iifo[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + eng_sse->Engine_sse::SetCurr(2,pos, curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + } + } + } + break; + } + default: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = m_Op_UPML->ii[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetCurr(0,pos) + - m_Op_UPML->iifo[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->SetCurr(0,pos, curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->ii[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetCurr(1,pos) + - m_Op_UPML->iifo[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->SetCurr(1,pos, curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + + f_help = m_Op_UPML->ii[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetCurr(2,pos) + - m_Op_UPML->iifo[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + m_Eng->SetCurr(2,pos, curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = f_help; + } + } + } + break; + } + } +} + +void Engine_Ext_UPML::DoPostCurrentUpdates(int threadID) +{ + if (m_Eng==NULL) + return; + if (threadID>=m_NrThreads) + return; + + unsigned int pos[3]; + unsigned int loc_pos[3]; + FDTD_FLOAT f_help; + + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->Engine::GetCurr(0,pos); + m_Eng->Engine::SetCurr(0,pos, f_help + m_Op_UPML->iifn[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->Engine::GetCurr(1,pos); + m_Eng->Engine::SetCurr(1,pos, f_help + m_Op_UPML->iifn[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->Engine::GetCurr(2,pos); + m_Eng->Engine::SetCurr(2,pos, f_help + m_Op_UPML->iifn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + } + } + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*) m_Eng; + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = eng_sse->Engine_sse::GetCurr(0,pos); + eng_sse->Engine_sse::SetCurr(0,pos, f_help + m_Op_UPML->iifn[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = eng_sse->Engine_sse::GetCurr(1,pos); + eng_sse->Engine_sse::SetCurr(1,pos, f_help + m_Op_UPML->iifn[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = eng_sse->Engine_sse::GetCurr(2,pos); + eng_sse->Engine_sse::SetCurr(2,pos, f_help + m_Op_UPML->iifn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + } + } + } + break; + } + default: + { + for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX) + { + loc_pos[0]=lineX+m_start.at(threadID); + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + f_help = curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetCurr(0,pos); + m_Eng->SetCurr(0,pos, f_help + m_Op_UPML->iifn[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetCurr(1,pos); + m_Eng->SetCurr(1,pos, f_help + m_Op_UPML->iifn[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + + f_help = curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetCurr(2,pos); + m_Eng->SetCurr(2,pos, f_help + m_Op_UPML->iifn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + } + } + } + break; + } + } +} diff --git a/openEMS/FDTD/extensions/engine_ext_upml.h b/openEMS/FDTD/extensions/engine_ext_upml.h new file mode 100644 index 0000000..283c886 --- /dev/null +++ b/openEMS/FDTD/extensions/engine_ext_upml.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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef ENGINE_EXT_UPML_H +#define ENGINE_EXT_UPML_H + +#include "engine_extension.h" +#include "FDTD/engine.h" +#include "FDTD/operator.h" + +class Operator_Ext_UPML; + +class Engine_Ext_UPML : public Engine_Extension +{ +public: + Engine_Ext_UPML(Operator_Ext_UPML* op_ext); + virtual ~Engine_Ext_UPML(); + + virtual void SetNumberOfThreads(int nrThread); + + virtual void DoPreVoltageUpdates() {Engine_Ext_UPML::DoPreVoltageUpdates(0);}; + virtual void DoPreVoltageUpdates(int threadID); + virtual void DoPostVoltageUpdates() {Engine_Ext_UPML::DoPostVoltageUpdates(0);}; + virtual void DoPostVoltageUpdates(int threadID); + + virtual void DoPreCurrentUpdates() {Engine_Ext_UPML::DoPreCurrentUpdates(0);}; + virtual void DoPreCurrentUpdates(int threadID); + virtual void DoPostCurrentUpdates() {Engine_Ext_UPML::DoPostCurrentUpdates(0);}; + virtual void DoPostCurrentUpdates(int threadID); + +protected: + Operator_Ext_UPML* m_Op_UPML; + + vector<unsigned int> m_start; + vector<unsigned int> m_numX; + + FDTD_FLOAT**** volt_flux; + FDTD_FLOAT**** curr_flux; +}; + +#endif // ENGINE_EXT_UPML_H diff --git a/openEMS/FDTD/extensions/engine_extension.cpp b/openEMS/FDTD/extensions/engine_extension.cpp new file mode 100644 index 0000000..6f688a3 --- /dev/null +++ b/openEMS/FDTD/extensions/engine_extension.cpp @@ -0,0 +1,95 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "engine_extension.h" +#include "operator_extension.h" + +#include "FDTD/engine.h" + +Engine_Extension::Engine_Extension(Operator_Extension* op_ext) +{ + m_Op_ext = op_ext; + m_Eng = NULL; + m_Priority = ENG_EXT_PRIO_DEFAULT; + m_NrThreads = 1; +} + +Engine_Extension::~Engine_Extension() +{ +} + +void Engine_Extension::SetNumberOfThreads(int nrThread) +{ + if (nrThread<1) + return; + m_NrThreads=nrThread; +} + +string Engine_Extension::GetExtensionName() const +{ + if (m_Op_ext) + return m_Op_ext->GetExtensionName(); + else + return "Unknown Extension"; +} + +void Engine_Extension::DoPreVoltageUpdates(int threadID) +{ + //if this method gets called the derived extension obviously doesn't support multithrading, calling non-MT method... + if (threadID==0) + DoPreVoltageUpdates(); +} + +void Engine_Extension::DoPostVoltageUpdates(int threadID) +{ + //if this method gets called the derived extension obviously doesn't support multithrading, calling non-MT method... + if (threadID==0) + DoPostVoltageUpdates(); +} + +void Engine_Extension::Apply2Voltages(int threadID) +{ + //if this method gets called the derived extension obviously doesn't support multithrading, calling non-MT method... + if (threadID==0) + Apply2Voltages(); +} + +void Engine_Extension::DoPreCurrentUpdates(int threadID) +{ + //if this method gets called the derived extension obviously doesn't support multithrading, calling non-MT method... + if (threadID==0) + DoPreCurrentUpdates(); +} + +void Engine_Extension::DoPostCurrentUpdates(int threadID) +{ + //if this method gets called the derived extension obviously doesn't support multithrading, calling non-MT method... + if (threadID==0) + DoPostCurrentUpdates(); +} + +void Engine_Extension::Apply2Current(int threadID) +{ + //if this method gets called the derived extension obviously doesn't support multithrading, calling non-MT method... + if (threadID==0) + Apply2Current(); +} + +bool Engine_Extension::operator< (const Engine_Extension& other) +{ + return (GetPriority()<other.GetPriority()); +} diff --git a/openEMS/FDTD/extensions/engine_extension.h b/openEMS/FDTD/extensions/engine_extension.h new file mode 100644 index 0000000..12ef712 --- /dev/null +++ b/openEMS/FDTD/extensions/engine_extension.h @@ -0,0 +1,88 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef ENGINE_EXTENSION_H +#define ENGINE_EXTENSION_H + +#define ENG_EXT_PRIO_DEFAULT 0 //default engine extension priority + +// priority definitions for some important extensions +#define ENG_EXT_PRIO_STEADYSTATE +2e6 //steady state extension priority +#define ENG_EXT_PRIO_UPML +1e6 //unaxial pml extension priority +#define ENG_EXT_PRIO_CYLINDER +1e5 //cylindrial extension priority +#define ENG_EXT_PRIO_TFSF +5e4 //total-field/scattered-field extension priority +#define ENG_EXT_PRIO_EXCITATION -1000 //excitation priority +#define ENG_EXT_PRIO_CYLINDERMULTIGRID -3000 //cylindrial multi-grid extension priority + +#include <string> + +class Operator_Extension; +class Engine; + +//! Abstract base-class for all engine extensions +class Engine_Extension +{ +public: + virtual ~Engine_Extension(); + + virtual void SetNumberOfThreads(int nrThread); + + //! This methode will be called __before__ the main engine does the usual voltage updates. This methode may __not__ change the engine voltages!!! + virtual void DoPreVoltageUpdates() {} + virtual void DoPreVoltageUpdates(int threadID); + //! This methode will be called __after__ the main engine does the usual voltage updates. This methode may __not__ change the engine voltages!!! + virtual void DoPostVoltageUpdates() {} + virtual void DoPostVoltageUpdates(int threadID); + //! This methode will be called __after__ all updates to the voltages and extensions and may add/set its results to the engine voltages, but may __not__ rely on the current value of the engine voltages!!! + virtual void Apply2Voltages() {} + virtual void Apply2Voltages(int threadID); + + //! This methode will be called __before__ the main engine does the usual current updates. This methode may __not__ change the engine current!!! + virtual void DoPreCurrentUpdates() {} + virtual void DoPreCurrentUpdates(int threadID); + //! This methode will be called __after__ the main engine does the usual current updates. This methode may __not__ change the engine current!!! + virtual void DoPostCurrentUpdates() {} + virtual void DoPostCurrentUpdates(int threadID); + //! This methode will be called __after__ all updates to the current and extensions and may add/set its results to the engine current, but may __not__ rely on the current value of the engine current!!! + virtual void Apply2Current() {} + virtual void Apply2Current(int threadID); + + //! Set the Engine to this extention. This will usually done automatically by Engine::AddExtension + virtual void SetEngine(Engine* eng) {m_Eng=eng;} + + //! Get the priority for this extension + virtual int GetPriority() const {return m_Priority;} + + //! Set the priority for this extension + virtual void SetPriority(int val) {m_Priority=val;} + + virtual bool operator< (const Engine_Extension& other); + + virtual std::string GetExtensionName() const; + +protected: + Engine_Extension(Operator_Extension* op_ext); + + Operator_Extension* m_Op_ext; + Engine* m_Eng; + + int m_Priority; + + int m_NrThreads; +}; + +#endif // ENGINE_EXTENSION_H diff --git a/openEMS/FDTD/extensions/operator_ext_conductingsheet.cpp b/openEMS/FDTD/extensions/operator_ext_conductingsheet.cpp new file mode 100644 index 0000000..28c51bb --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_conductingsheet.cpp @@ -0,0 +1,261 @@ +/* +* Copyright (C) 2012 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 <http://www.gnu.org/licenses/>. +*/ + +#include "operator_ext_conductingsheet.h" +#include "tools/array_ops.h" +#include "tools/constants.h" +#include "cond_sheet_parameter.h" + +#include "CSPropConductingSheet.h" + +Operator_Ext_ConductingSheet::Operator_Ext_ConductingSheet(Operator* op, double f_max) : Operator_Ext_LorentzMaterial(op) +{ + m_f_max = f_max; +} + +Operator_Ext_ConductingSheet::Operator_Ext_ConductingSheet(Operator* op, Operator_Ext_ConductingSheet* op_ext) : Operator_Ext_LorentzMaterial(op, op_ext) +{ + m_f_max = op_ext->m_f_max; +} + +Operator_Extension* Operator_Ext_ConductingSheet::Clone(Operator* op) +{ + if (dynamic_cast<Operator_Ext_ConductingSheet*>(this)==NULL) + return NULL; + return new Operator_Ext_ConductingSheet(op, this); +} + +bool Operator_Ext_ConductingSheet::BuildExtension() +{ + double dT = m_Op->GetTimestep(); + unsigned int pos[] = {0,0,0}; + double coord[3]; + unsigned int numLines[3] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)}; + + m_Order = 0; + vector<unsigned int> v_pos[3]; + int ****tanDir = Create_N_3DArray<int>(numLines); + float ****Conductivity = Create_N_3DArray<float>(numLines); + float ****Thickness = Create_N_3DArray<float>(numLines); + + CSPrimitives* cs_sheet = NULL; + double box[6]; + int nP, nPP; + bool b_pos_on; + bool disable_pos; + for (pos[0]=0; pos[0]<numLines[0]; ++pos[0]) + { + for (pos[1]=0; pos[1]<numLines[1]; ++pos[1]) + { + vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL)); + for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) + { + b_pos_on = false; + disable_pos = false; + // disable conducting sheet model inside the boundary conditions, especially inside a pml + for (int m=0;m<3;++m) + if ((pos[m]<=(unsigned int)m_Op->GetBCSize(2*m)) || (pos[m]>=(numLines[m]-m_Op->GetBCSize(2*m+1)-1))) + disable_pos = true; + + for (int n=0; n<3; ++n) + { + nP = (n+1)%3; + nPP = (n+2)%3; + + tanDir[n][pos[0]][pos[1]][pos[2]] = -1; //deactivate by default + Conductivity[n][pos[0]][pos[1]][pos[2]] = 0; //deactivate by default + Thickness[n][pos[0]][pos[1]][pos[2]] = 0; //deactivate by default + + if (m_Op->GetYeeCoords(n,pos,coord,false)==false) + continue; + + // Ez at r==0 not supported --> set to PEC + if (m_CC_R0_included && (n==2) && (pos[0]==0)) + disable_pos = true; + +// CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), false, &cs_sheet); + CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, false, &cs_sheet); + CSPropConductingSheet* cs_prop = dynamic_cast<CSPropConductingSheet*>(prop); + if (cs_prop) + { + if (cs_sheet==NULL) + return false; //sanity check, this should never happen + if (cs_sheet->GetDimension()!=2) + { + cerr << "Operator_Ext_ConductingSheet::BuildExtension: A conducting sheet primitive (ID: " << cs_sheet->GetID() << ") with dimension: " << cs_sheet->GetDimension() << " found, fallback to PEC!" << endl; + m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 ); + m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 ); + ++m_Op->m_Nr_PEC[n]; + continue; + } + cs_sheet->SetPrimitiveUsed(true); + + if (disable_pos) + { + m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 ); + m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 ); + ++m_Op->m_Nr_PEC[n]; + continue; + } + + Conductivity[n][pos[0]][pos[1]][pos[2]] = cs_prop->GetConductivity(); + Thickness[n][pos[0]][pos[1]][pos[2]] = cs_prop->GetThickness(); + + if ((Conductivity[n][pos[0]][pos[1]][pos[2]]<=0) || (Thickness[n][pos[0]][pos[1]][pos[2]]<=0)) + { + cerr << "Operator_Ext_ConductingSheet::BuildExtension: Warning: Zero conductivity or thickness detected... fallback to PEC!" << endl; + m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 ); + m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 ); + ++m_Op->m_Nr_PEC[n]; + continue; + } + + cs_sheet->GetBoundBox(box); + if (box[2*nP]!=box[2*nP+1]) + tanDir[n][pos[0]][pos[1]][pos[2]] = nP; + if (box[2*nPP]!=box[2*nPP+1]) + tanDir[n][pos[0]][pos[1]][pos[2]] = nPP; + b_pos_on = true; + } + } + if (b_pos_on) + { + for (int n=0; n<3; ++n) + v_pos[n].push_back(pos[n]); + } + } + } + } + + size_t numCS = v_pos[0].size(); + if (numCS==0) + return false; + + m_LM_Count.push_back(numCS); + m_LM_Count.push_back(numCS); + + m_Order = 2; + m_volt_ADE_On = new bool[m_Order]; + m_volt_ADE_On[0] = m_volt_ADE_On[1]=true; + m_curr_ADE_On = new bool[m_Order]; + m_curr_ADE_On[0] = m_curr_ADE_On[1]=false; + + m_volt_Lor_ADE_On = new bool[m_Order]; + m_volt_Lor_ADE_On[0] = m_volt_Lor_ADE_On[1]=false; + m_curr_Lor_ADE_On = new bool[m_Order]; + m_curr_Lor_ADE_On[0] = m_curr_Lor_ADE_On[1]=false; + + m_LM_pos = new unsigned int**[m_Order]; + m_LM_pos[0] = new unsigned int*[3]; + m_LM_pos[1] = new unsigned int*[3]; + + v_int_ADE = new FDTD_FLOAT**[m_Order]; + v_ext_ADE = new FDTD_FLOAT**[m_Order]; + + v_int_ADE[0] = new FDTD_FLOAT*[3]; + v_ext_ADE[0] = new FDTD_FLOAT*[3]; + v_int_ADE[1] = new FDTD_FLOAT*[3]; + v_ext_ADE[1] = new FDTD_FLOAT*[3]; + + for (int n=0; n<3; ++n) + { + m_LM_pos[0][n] = new unsigned int[numCS]; + m_LM_pos[1][n] = new unsigned int[numCS]; + for (unsigned int i=0; i<numCS; ++i) + { + m_LM_pos[0][n][i] = v_pos[n].at(i); + m_LM_pos[1][n][i] = v_pos[n].at(i); + } + v_int_ADE[0][n] = new FDTD_FLOAT[numCS]; + v_int_ADE[1][n] = new FDTD_FLOAT[numCS]; + v_ext_ADE[0][n] = new FDTD_FLOAT[numCS]; + v_ext_ADE[1][n] = new FDTD_FLOAT[numCS]; + } + + unsigned int index; + float w_stop = m_f_max*2*PI; + float Omega_max=0; + float G,L1,L2,R1,R2,Lmin; + float G0, w0; + float wtl; //width to length factor + float factor=1; + int t_dir=0; //tangential sheet direction + unsigned int tpos[] = {0,0,0}; + unsigned int optParaPos; + for (unsigned int i=0;i<numCS;++i) + { + pos[0]=m_LM_pos[0][0][i];pos[1]=m_LM_pos[0][1][i];pos[2]=m_LM_pos[0][2][i]; + tpos[0]=pos[0];tpos[1]=pos[1];tpos[2]=pos[2]; + index = m_Op->MainOp->SetPos(pos[0],pos[1],pos[2]); + for (int n=0;n<3;++n) + { + tpos[0]=pos[0];tpos[1]=pos[1];tpos[2]=pos[2]; + t_dir = tanDir[n][pos[0]][pos[1]][pos[2]]; + G0 = Conductivity[n][pos[0]][pos[1]][pos[2]]*Thickness[n][pos[0]][pos[1]][pos[2]]; + w0 = 8.0/ G0 / Thickness[n][pos[0]][pos[1]][pos[2]]/__MUE0__; + Omega_max = w_stop/w0; + for (optParaPos=0;optParaPos<numOptPara;++optParaPos) + if (omega_stop[optParaPos]>Omega_max) + break; + if (optParaPos>=numOptPara) + { + cerr << "Operator_Ext_ConductingSheet::BuildExtension(): Error, conductor thickness, conductivity or max. simulation frequency of interest is too high! Check parameter!" << endl; + cerr << " --> max f: " << m_f_max << "Hz, Conductivity: " << Conductivity[n][pos[0]][pos[1]][pos[2]] << "S/m, Thickness " << Thickness[n][pos[0]][pos[1]][pos[2]]*1e6 << "um" << endl; + optParaPos = numOptPara-1; + } + v_int_ADE[0][n][i]=0; + v_ext_ADE[0][n][i]=0; + v_int_ADE[1][n][i]=0; + v_ext_ADE[1][n][i]=0; + if (t_dir>=0) + { + wtl = m_Op->GetEdgeLength(n,pos)/m_Op->GetNodeWidth(t_dir,pos); + factor = 1; + if (tanDir[t_dir][tpos[0]][tpos[1]][tpos[2]]<0) + factor = 2; + --tpos[t_dir]; + if (tanDir[t_dir][tpos[0]][tpos[1]][tpos[2]]<0) + factor = 2; + + L1 = l1[optParaPos]/G0/w0*factor; + L2 = l2[optParaPos]/G0/w0*factor; + R1 = r1[optParaPos]/G0*factor; + R2 = r2[optParaPos]/G0*factor; + G = G0*g[optParaPos]/factor; + + L1*=wtl; + L2*=wtl; + R1*=wtl; + R2*=wtl; + G/=wtl; + + Lmin = L1; + if (L2<L1) + Lmin = L2; + m_Op->EC_G[n][index]= G; + m_Op->EC_C[n][index]= dT*dT/4.0*(16.0/Lmin + 1/L1 + 1/L2); + m_Op->Calc_ECOperatorPos(n,pos); + + v_int_ADE[0][n][i]=(2.0*L1-dT*R1)/(2.0*L1+dT*R1); + v_ext_ADE[0][n][i]=dT/(L1+dT*R1/2.0)*m_Op->GetVI(n,pos[0],pos[1],pos[2]); + v_int_ADE[1][n][i]=(2.0*L2-dT*R2)/(2.0*L2+dT*R2); + v_ext_ADE[1][n][i]=dT/(L2+dT*R2/2.0)*m_Op->GetVI(n,pos[0],pos[1],pos[2]); + } + } + } + return true; +} diff --git a/openEMS/FDTD/extensions/operator_ext_conductingsheet.h b/openEMS/FDTD/extensions/operator_ext_conductingsheet.h new file mode 100644 index 0000000..f84dc4b --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_conductingsheet.h @@ -0,0 +1,51 @@ +/* +* Copyright (C) 2012 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef OPERATOR_EXT_CONDUCTINGSHEET_H +#define OPERATOR_EXT_CONDUCTINGSHEET_H + +#include "FDTD/operator.h" +#include "operator_ext_lorentzmaterial.h" + +/*! + FDTD extension for a conducting sheet model as described in: + Lauer, A.; Wolff, I.; , "A conducting sheet model for efficient wide band FDTD analysis of planar waveguides and circuits," Microwave Symposium Digest, 1999 IEEE MTT-S International , vol.4, no., pp.1589-1592 vol.4, 1999 + doi: 10.1109/MWSYM.1999.780262 + URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=780262&isnumber=16934 + */ +class Operator_Ext_ConductingSheet : public Operator_Ext_LorentzMaterial +{ +public: + Operator_Ext_ConductingSheet(Operator* op, double f_max); + + virtual Operator_Extension* Clone(Operator* op); + + virtual bool BuildExtension(); + + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return true;} + virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return true;} + virtual bool IsMPISave() const {return true;} + + virtual string GetExtensionName() const {return string("Conducting Sheet Extension");} + +protected: + //! Copy constructor + Operator_Ext_ConductingSheet(Operator* op, Operator_Ext_ConductingSheet* op_ext); + double m_f_max; +}; + +#endif // OPERATOR_EXT_CONDUCTINGSHEET_H diff --git a/openEMS/FDTD/extensions/operator_ext_cylinder.cpp b/openEMS/FDTD/extensions/operator_ext_cylinder.cpp new file mode 100644 index 0000000..d400bf8 --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_cylinder.cpp @@ -0,0 +1,114 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "operator_ext_cylinder.h" +#include "FDTD/operator_cylinder.h" +#include "engine_ext_cylinder.h" + +Operator_Ext_Cylinder::Operator_Ext_Cylinder(Operator_Cylinder* op) : Operator_Extension(op) +{ + m_Op_Cyl = op; + + CC_R0_included=m_Op_Cyl->GetR0Included(); + CC_closedAlpha=m_Op_Cyl->GetClosedAlpha(); + + vv_R0 = NULL; + vi_R0 = NULL; +} + +Operator_Ext_Cylinder::~Operator_Ext_Cylinder() +{ + delete[] vv_R0; + vv_R0=NULL; + delete[] vi_R0; + vi_R0=NULL; +} + +bool Operator_Ext_Cylinder::BuildExtension() +{ + delete[] vv_R0; + vv_R0=NULL; + delete[] vi_R0; + vi_R0=NULL; + + //if r=0 is not included -> obviously no special treatment for r=0 + //if alpha direction is not closed, PEC-BC at r=0 necessary and already set... + if (CC_R0_included==false) + return true; + + vv_R0 = new FDTD_FLOAT[m_Op->GetNumberOfLines(2,true)]; + vi_R0 = new FDTD_FLOAT[m_Op->GetNumberOfLines(2,true)]; + + unsigned int pos[3]; + double coord[3]; + double inEC[4]; + double dT = m_Op->GetTimestep(); + pos[0]=0; + vector<CSPrimitives*> vPrims_metal = m_Op->GetPrimitivesBoundBox(pos[0], -1, -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL)); + for (pos[2]=0; pos[2]<m_Op->GetNumberOfLines(2,true); ++pos[2]) + { + double C=0; + double G=0; + vector<CSPrimitives*> vPrims_mat = m_Op->GetPrimitivesBoundBox(pos[0], -1, pos[2], CSProperties::MATERIAL); + for (pos[1]=0; pos[1]<m_Op->GetNumberOfLines(1,true)-2; ++pos[1]) + { + m_Op_Cyl->Calc_ECPos(2,pos,inEC,vPrims_mat); + C+=inEC[0]; + G+=inEC[1]; + } + m_Op->SetVV(2,0,0,pos[2], 1); + vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C); + vi_R0[pos[2]] = (dT/C)/(1+dT*G/2/C); + + for (unsigned int i=0; i<m_Op->GetNumberOfLines(1,true); ++i) + { + m_Op->EC_C[2][m_Op->MainOp->SetPos(0,i,pos[2])] = C; + m_Op->EC_G[2][m_Op->MainOp->SetPos(0,i,pos[2])] = G; + } + + //search for metal on z-axis + m_Op_Cyl->GetYeeCoords(2,pos,coord,false); + CSProperties* prop = m_Op->CSX->GetPropertyByCoordPriority(coord, vPrims_metal, true); + if (prop) + { + if (prop->GetType()==CSProperties::METAL) //set to PEC + { + m_Op->SetVV(2,0,0,pos[2], 0); + vv_R0[pos[2]] = 0; + vi_R0[pos[2]] = 0; + m_Op->EC_C[2][m_Op->MainOp->SetPos(0,0,pos[2])] = 0; + m_Op->EC_G[2][m_Op->MainOp->SetPos(0,0,pos[2])] = 0; + } + } + } + return true; +} + +Engine_Extension* Operator_Ext_Cylinder::CreateEngineExtention() +{ + Engine_Ext_Cylinder* eng_ext = new Engine_Ext_Cylinder(this); + return eng_ext; +} + + +void Operator_Ext_Cylinder::ShowStat(ostream &ostr) const +{ + Operator_Extension::ShowStat(ostr); + string On_Off[2] = {"Off", "On"}; + ostr << " Zeroth Radius\t\t: " << On_Off[CC_R0_included] << endl; + ostr << " Closed Rotation\t: " << On_Off[CC_closedAlpha] << endl; +} diff --git a/openEMS/FDTD/extensions/operator_ext_cylinder.h b/openEMS/FDTD/extensions/operator_ext_cylinder.h new file mode 100644 index 0000000..c8dff2d --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_cylinder.h @@ -0,0 +1,60 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef OPERATOR_EXT_CYLINDER_H +#define OPERATOR_EXT_CYLINDER_H + +#include "operator_extension.h" +#include "FDTD/operator.h" + +class Operator_Cylinder; + +class Operator_Ext_Cylinder : public Operator_Extension +{ + friend class Engine_Ext_Cylinder; + friend class Operator_Ext_LorentzMaterial; +public: + Operator_Ext_Cylinder(Operator_Cylinder* op); + ~Operator_Ext_Cylinder(); + + virtual bool BuildExtension(); + + virtual Engine_Extension* CreateEngineExtention(); + + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return true;} + virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return true;} + + // FIXME, this extension is not save or unknown to be save to use with MPI + virtual bool IsMPISave() const {return false;} + + virtual std::string GetExtensionName() const {return std::string("Extension for the Cylinder-Coords Operator");} + + virtual void ShowStat(ostream &ostr) const; + +protected: + Operator_Cylinder* m_Op_Cyl; + + bool CC_closedAlpha; + bool CC_R0_included; + + //special EC operator for R0 + FDTD_FLOAT* vv_R0; //calc new voltage from old voltage + FDTD_FLOAT* vi_R0; //calc new voltage from old current + +}; + +#endif // OPERATOR_EXT_CYLINDER_H diff --git a/openEMS/FDTD/extensions/operator_ext_dispersive.cpp b/openEMS/FDTD/extensions/operator_ext_dispersive.cpp new file mode 100644 index 0000000..ed01f1d --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_dispersive.cpp @@ -0,0 +1,78 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "operator_ext_dispersive.h" + +#include "tools/array_ops.h" + +using namespace std; + +Operator_Ext_Dispersive::Operator_Ext_Dispersive(Operator* op) : Operator_Extension(op) +{ + m_curr_ADE_On = NULL; + m_volt_ADE_On = NULL; + + m_LM_pos=NULL; + m_curr_ADE_On=NULL; + m_volt_ADE_On=NULL; + + m_Order = 0; +} + +Operator_Ext_Dispersive::Operator_Ext_Dispersive(Operator* op, Operator_Ext_Dispersive* op_ext) : Operator_Extension(op,op_ext) +{ + m_curr_ADE_On = NULL; + m_volt_ADE_On = NULL; + + m_LM_pos=NULL; + m_curr_ADE_On=NULL; + m_volt_ADE_On=NULL; + + m_Order = 0; +} + +Operator_Ext_Dispersive::~Operator_Ext_Dispersive() +{ + delete[] m_curr_ADE_On; + delete[] m_volt_ADE_On; + m_curr_ADE_On=NULL; + m_volt_ADE_On=NULL; + + for (int n=0;n<m_Order;++n) + { + delete[] m_LM_pos[n][0]; + delete[] m_LM_pos[n][1]; + delete[] m_LM_pos[n][2]; + } + delete[] m_LM_pos; + m_LM_pos=NULL; + m_Order=0; + m_LM_Count.clear(); +} + +void Operator_Ext_Dispersive::ShowStat(ostream &ostr) const +{ + Operator_Extension::ShowStat(ostr); + string On_Off[2] = {"Off", "On"}; + ostr << " Max. Dispersion Order N = " << m_Order << endl; + for (int i=0;i<m_Order;++i) + { + ostr << " N=" << i << ":\t Active cells\t\t: " << m_LM_Count.at(i) << endl; + ostr << " N=" << i << ":\t Voltage ADE is \t: " << On_Off[m_volt_ADE_On[i]] << endl; + ostr << " N=" << i << ":\t Current ADE is \t: " << On_Off[m_curr_ADE_On[i]] << endl; + } +} diff --git a/openEMS/FDTD/extensions/operator_ext_dispersive.h b/openEMS/FDTD/extensions/operator_ext_dispersive.h new file mode 100644 index 0000000..fab354c --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_dispersive.h @@ -0,0 +1,56 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef OPERATOR_EXT_DISPERSIVE_H +#define OPERATOR_EXT_DISPERSIVE_H + +//#include "operator.h" +#include "operator_extension.h" +#include "vector" + +//! Abstract base class for all dispersive material models, based on an ADE (additional differential equation) +class Operator_Ext_Dispersive : public Operator_Extension +{ + friend class Engine_Ext_Dispersive; +public: + virtual ~Operator_Ext_Dispersive(); + + virtual int GetDispersionOrder() {return m_Order;} + + virtual std::string GetExtensionName() const {return std::string("Dispersive Material Abstract Base class");} + + virtual void ShowStat(std::ostream &ostr) const; + +protected: + Operator_Ext_Dispersive(Operator* op); + //! Copy constructor + Operator_Ext_Dispersive(Operator* op, Operator_Ext_Dispersive* op_ext); + + //! Dispersive order + int m_Order; + + //! Dispersive material count + std::vector<unsigned int> m_LM_Count; + //! Index with dispersive material + // Array setup: m_LM_pos[N_order][direction][mesh_pos] + unsigned int ***m_LM_pos; + + bool *m_curr_ADE_On; + bool *m_volt_ADE_On; +}; + +#endif // OPERATOR_EXT_DISPERSIVE_H diff --git a/openEMS/FDTD/extensions/operator_ext_excitation.cpp b/openEMS/FDTD/extensions/operator_ext_excitation.cpp new file mode 100644 index 0000000..d2085ff --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_excitation.cpp @@ -0,0 +1,372 @@ +/* +* Copyright (C) 2011 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 <http://www.gnu.org/licenses/>. +*/ + +#include "operator_ext_excitation.h" +#include "engine_ext_excitation.h" +#include "FDTD/excitation.h" +#include "ContinuousStructure.h" + +#include "CSPrimCurve.h" +#include "CSPropExcitation.h" + +Operator_Ext_Excitation::Operator_Ext_Excitation(Operator* op) : Operator_Extension(op) +{ + Init(); +} + +Operator_Ext_Excitation::~Operator_Ext_Excitation() +{ + Reset(); +} + +Operator_Extension* Operator_Ext_Excitation::Clone(Operator* op) +{ + Operator_Ext_Excitation* clone = new Operator_Ext_Excitation(op, this); + return clone; +} + +void Operator_Ext_Excitation::Init() +{ + Operator_Extension::Init(); + Volt_delay = 0; + Volt_amp = 0; + Volt_dir = 0; + Volt_Count = 0; + Curr_delay = 0; + Curr_amp = 0; + Curr_dir = 0; + Curr_Count = 0; + + for (int n=0; n<3; ++n) + { + Volt_index[n] = 0; + Curr_index[n] = 0; + Volt_Count_Dir[n] = 0; + Curr_Count_Dir[n] = 0; + } + m_Exc = 0; +} + +void Operator_Ext_Excitation::Reset() +{ + Operator_Extension::Reset(); + delete[] Volt_delay; + Volt_delay = 0; + delete[] Volt_dir; + Volt_dir = 0; + delete[] Volt_amp; + Volt_amp = 0; + delete[] Curr_delay; + Curr_delay = 0; + delete[] Curr_dir; + Curr_dir = 0; + delete[] Curr_amp; + Curr_amp = 0; + + Volt_Count = 0; + Curr_Count = 0; + + for (int n=0; n<3; ++n) + { + delete[] Volt_index[n]; + Volt_index[n] = 0; + delete[] Curr_index[n]; + Curr_index[n] = 0; + + Volt_Count_Dir[n] = 0; + Curr_Count_Dir[n] = 0; + } +} + + +Operator_Ext_Excitation::Operator_Ext_Excitation(Operator* op, Operator_Ext_Excitation* op_ext) : Operator_Extension(op, op_ext) +{ + Init(); +} + +bool Operator_Ext_Excitation::BuildExtension() +{ + m_Exc = m_Op->GetExcitationSignal(); + double dT = m_Op->GetTimestep(); + if (dT==0) + return false; + if (m_Exc==0) + return false; + + Reset(); + ContinuousStructure* CSX = m_Op->GetGeometryCSX(); + + unsigned int pos[3]; + double amp=0; + + vector<unsigned int> volt_vIndex[3]; + vector<FDTD_FLOAT> volt_vExcit; + vector<unsigned int> volt_vDelay; + vector<unsigned int> volt_vDir; + double volt_coord[3]; + + vector<unsigned int> curr_vIndex[3]; + vector<FDTD_FLOAT> curr_vExcit; + vector<unsigned int> curr_vDelay; + vector<unsigned int> curr_vDir; + double curr_coord[3]; + + vector<CSProperties*> vec_prop = CSX->GetPropertyByType(CSProperties::EXCITATION); + + if (vec_prop.size()==0) + { + cerr << "Operator::CalcFieldExcitation: Warning, no excitation properties found" << endl; + return false; + } + + CSPropExcitation* elec=NULL; + CSProperties* prop=NULL; + int priority=0; + + unsigned int numLines[] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)}; + + for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) + { + for (pos[1]=0; pos[1]<numLines[1]; ++pos[1]) + { + for (pos[0]=0; pos[0]<numLines[0]; ++pos[0]) + { + //electric field excite + for (int n=0; n<3; ++n) + { + if (m_Op->GetYeeCoords(n,pos,volt_coord,false)==false) + continue; + if (m_CC_R0_included && (n==2) && (pos[0]==0)) + volt_coord[1] = m_Op->GetDiscLine(1,0); + + if (m_CC_R0_included && (n==1) && (pos[0]==0)) + continue; + + for (size_t p=0; p<vec_prop.size(); ++p) + { + prop = vec_prop.at(p); + elec = prop->ToExcitation(); + if (elec==NULL) + continue; + if (prop->CheckCoordInPrimitive(volt_coord,priority,true)) + { + if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==0) || (elec->GetExcitType()==1) ))//&& (pos[n]<numLines[n]-1)) + { + amp = elec->GetWeightedExcitation(n,volt_coord)*m_Op->GetEdgeLength(n,pos);// delta[n]*gridDelta; + if (amp!=0) + { + volt_vExcit.push_back(amp); + volt_vDelay.push_back((unsigned int)(elec->GetDelay()/dT)); + volt_vDir.push_back(n); + volt_vIndex[0].push_back(pos[0]); + volt_vIndex[1].push_back(pos[1]); + volt_vIndex[2].push_back(pos[2]); + } + if (elec->GetExcitType()==1) //hard excite + { + m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 ); + m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 ); + } + } + } + } + } + + //magnetic field excite + for (int n=0; n<3; ++n) + { + if ((pos[0]>=numLines[0]-1) || (pos[1]>=numLines[1]-1) || (pos[2]>=numLines[2]-1)) + continue; //skip the last H-Line which is outside the FDTD-domain + if (m_Op->GetYeeCoords(n,pos,curr_coord,true)==false) + continue; + for (size_t p=0; p<vec_prop.size(); ++p) + { + prop = vec_prop.at(p); + elec = prop->ToExcitation(); + if (elec==NULL) + continue; + if (prop->CheckCoordInPrimitive(curr_coord,priority,true)) + { + if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==2) || (elec->GetExcitType()==3) )) + { + amp = elec->GetWeightedExcitation(n,curr_coord)*m_Op->GetEdgeLength(n,pos,true);// delta[n]*gridDelta; + if (amp!=0) + { + curr_vExcit.push_back(amp); + curr_vDelay.push_back((unsigned int)(elec->GetDelay()/dT)); + curr_vDir.push_back(n); + curr_vIndex[0].push_back(pos[0]); + curr_vIndex[1].push_back(pos[1]); + curr_vIndex[2].push_back(pos[2]); + } + if (elec->GetExcitType()==3) //hard excite + { + m_Op->SetII(n,pos[0],pos[1],pos[2], 0 ); + m_Op->SetIV(n,pos[0],pos[1],pos[2], 0 ); + } + } + } + } + } + + } + } + } + + //special treatment for primitives of type curve (treated as wires) see also Calc_PEC + double p1[3]; + double p2[3]; + Grid_Path path; + for (size_t p=0; p<vec_prop.size(); ++p) + { + prop = vec_prop.at(p); + elec = prop->ToExcitation(); + for (size_t n=0; n<prop->GetQtyPrimitives(); ++n) + { + CSPrimitives* prim = prop->GetPrimitive(n); + CSPrimCurve* curv = prim->ToCurve(); + if (curv) + { + for (size_t i=1; i<curv->GetNumberOfPoints(); ++i) + { + curv->GetPoint(i-1,p1,m_Op->m_MeshType); + curv->GetPoint(i,p2,m_Op->m_MeshType); + path = m_Op->FindPath(p1,p2); + if (path.dir.size()>0) + prim->SetPrimitiveUsed(true); + for (size_t t=0; t<path.dir.size(); ++t) + { + n = path.dir.at(t); + pos[0] = path.posPath[0].at(t); + pos[1] = path.posPath[1].at(t); + pos[2] = path.posPath[2].at(t); + m_Op->GetYeeCoords(n,pos,volt_coord,false); + if (elec!=NULL) + { + if ((elec->GetActiveDir(n)) && (pos[n]<numLines[n]-1) && ( (elec->GetExcitType()==0) || (elec->GetExcitType()==1) )) + { + amp = elec->GetWeightedExcitation(n,volt_coord)*m_Op->GetEdgeLength(n,pos); + if (amp!=0) + { + volt_vExcit.push_back(amp); + volt_vDelay.push_back((unsigned int)(elec->GetDelay()/dT)); + volt_vDir.push_back(n); + volt_vIndex[0].push_back(pos[0]); + volt_vIndex[1].push_back(pos[1]); + volt_vIndex[2].push_back(pos[2]); + } + if (elec->GetExcitType()==1) //hard excite + { + m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 ); + m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 ); + } + } + } + } + } + } + } + } + + // set voltage excitations + setupVoltageExcitation( volt_vIndex, volt_vExcit, volt_vDelay, volt_vDir ); + + // set current excitations + setupCurrentExcitation( curr_vIndex, curr_vExcit, curr_vDelay, curr_vDir ); + + return true; +} + +void Operator_Ext_Excitation::setupVoltageExcitation( vector<unsigned int> const volt_vIndex[3], vector<FDTD_FLOAT> const& volt_vExcit, + vector<unsigned int> const& volt_vDelay, vector<unsigned int> const& volt_vDir ) +{ + Volt_Count = volt_vIndex[0].size(); + for (int n=0; n<3; n++) + { + Volt_Count_Dir[n]=0; + delete[] Volt_index[n]; + Volt_index[n] = new unsigned int[Volt_Count]; + } + delete[] Volt_delay; + delete[] Volt_amp; + delete[] Volt_dir; + Volt_delay = new unsigned int[Volt_Count]; + Volt_amp = new FDTD_FLOAT[Volt_Count]; + Volt_dir = new unsigned short[Volt_Count]; + +// cerr << "Excitation::setupVoltageExcitation(): Number of voltage excitation points: " << Volt_Count << endl; +// if (Volt_Count==0) +// cerr << "No E-Field/voltage excitation found!" << endl; + for (int n=0; n<3; n++) + for (unsigned int i=0; i<Volt_Count; i++) + Volt_index[n][i] = volt_vIndex[n].at(i); + for (unsigned int i=0; i<Volt_Count; i++) + { + Volt_delay[i] = volt_vDelay.at(i); + Volt_amp[i] = volt_vExcit.at(i); + Volt_dir[i] = volt_vDir.at(i); + ++Volt_Count_Dir[Volt_dir[i]]; + } +} + +void Operator_Ext_Excitation::setupCurrentExcitation( vector<unsigned int> const curr_vIndex[3], vector<FDTD_FLOAT> const& curr_vExcit, + vector<unsigned int> const& curr_vDelay, vector<unsigned int> const& curr_vDir ) +{ + Curr_Count = curr_vIndex[0].size(); + for (int n=0; n<3; n++) + { + Curr_Count_Dir[n]=0; + delete[] Curr_index[n]; + Curr_index[n] = new unsigned int[Curr_Count]; + } + delete[] Curr_delay; + delete[] Curr_amp; + delete[] Curr_dir; + Curr_delay = new unsigned int[Curr_Count]; + Curr_amp = new FDTD_FLOAT[Curr_Count]; + Curr_dir = new unsigned short[Curr_Count]; + +// cerr << "Excitation::setupCurrentExcitation(): Number of current excitation points: " << Curr_Count << endl; +// if (Curr_Count==0) +// cerr << "No H-Field/current excitation found!" << endl; + for (int n=0; n<3; ++n) + for (unsigned int i=0; i<Curr_Count; i++) + Curr_index[n][i] = curr_vIndex[n].at(i); + for (unsigned int i=0; i<Curr_Count; i++) + { + Curr_delay[i] = curr_vDelay.at(i); + Curr_amp[i] = curr_vExcit.at(i); + Curr_dir[i] = curr_vDir.at(i); + ++Curr_Count_Dir[Curr_dir[i]]; + } + +} + +Engine_Extension* Operator_Ext_Excitation::CreateEngineExtention() +{ + return new Engine_Ext_Excitation(this); +} + +void Operator_Ext_Excitation::ShowStat(ostream &ostr) const +{ + Operator_Extension::ShowStat(ostr); + cout << "Voltage excitations\t: " << Volt_Count << "\t (" << Volt_Count_Dir[0] << ", " << Volt_Count_Dir[1] << ", " << Volt_Count_Dir[2] << ")" << endl; + cout << "Current excitations\t: " << Curr_Count << "\t (" << Curr_Count_Dir[0] << ", " << Curr_Count_Dir[1] << ", " << Curr_Count_Dir[2] << ")" << endl; + cout << "Excitation Length (TS)\t: " << m_Exc->GetLength() << endl; + cout << "Excitation Length (s)\t: " << m_Exc->GetLength()*m_Op->GetTimestep() << endl; +} + diff --git a/openEMS/FDTD/extensions/operator_ext_excitation.h b/openEMS/FDTD/extensions/operator_ext_excitation.h new file mode 100644 index 0000000..2abcef0 --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_excitation.h @@ -0,0 +1,84 @@ +/* +* Copyright (C) 2011 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef OPERATOR_EXT_EXCITATION_H +#define OPERATOR_EXT_EXCITATION_H + +#include "operator_extension.h" +#include "FDTD/operator.h" + +class Excitation; + +class Operator_Ext_Excitation : public Operator_Extension +{ + friend class Engine_Ext_Excitation; + friend class Engine_Ext_Mur_ABC; + friend class Operator; +public: + Operator_Ext_Excitation(Operator* op); + ~Operator_Ext_Excitation(); + + virtual Operator_Extension* Clone(Operator* op); + + virtual bool BuildExtension(); + + virtual Engine_Extension* CreateEngineExtention(); + + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return true;} + virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return true;} + virtual bool IsMPISave() const {return true;} + + virtual string GetExtensionName() const {return string("Excitation Extension");} + + virtual void ShowStat(ostream &ostr) const; + + virtual void Init(); + virtual void Reset(); + + unsigned int GetVoltCount() const {return Volt_Count;} + unsigned int GetVoltCount(int ny) const {return Volt_Count_Dir[ny];} + + unsigned int GetCurrCount() const {return Curr_Count;} + unsigned int GetCurrCount(int ny) const {return Curr_Count_Dir[ny];} + +protected: + Operator_Ext_Excitation(Operator* op, Operator_Ext_Excitation* op_ext); + + Excitation* m_Exc; + + void setupVoltageExcitation( vector<unsigned int> const volt_vIndex[3], vector<FDTD_FLOAT> const& volt_vExcit, + vector<unsigned int> const& volt_vDelay, vector<unsigned int> const& volt_vDir ); + void setupCurrentExcitation( vector<unsigned int> const curr_vIndex[3], vector<FDTD_FLOAT> const& curr_vExcit, + vector<unsigned int> const& curr_vDelay, vector<unsigned int> const& curr_vDir ); + //E-Field/voltage Excitation + unsigned int Volt_Count; + unsigned int Volt_Count_Dir[3]; + unsigned int* Volt_index[3]; + unsigned short* Volt_dir; + FDTD_FLOAT* Volt_amp; //represented as edge-voltages!! + unsigned int* Volt_delay; + + //H-Field/current Excitation + unsigned int Curr_Count; + unsigned int Curr_Count_Dir[3]; + unsigned int* Curr_index[3]; + unsigned short* Curr_dir; + FDTD_FLOAT* Curr_amp; //represented as edge-currents!! + unsigned int* Curr_delay; +}; + +#endif // OPERATOR_EXT_EXCITATION_H diff --git a/openEMS/FDTD/extensions/operator_ext_lorentzmaterial.cpp b/openEMS/FDTD/extensions/operator_ext_lorentzmaterial.cpp new file mode 100644 index 0000000..9f386fc --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_lorentzmaterial.cpp @@ -0,0 +1,453 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "operator_ext_lorentzmaterial.h" +#include "engine_ext_lorentzmaterial.h" +#include "operator_ext_cylinder.h" +#include "../operator_cylinder.h" + +#include "CSPropLorentzMaterial.h" +#include "CSPropDebyeMaterial.h" + +Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op) : Operator_Ext_Dispersive(op) +{ + v_int_ADE = NULL; + v_ext_ADE = NULL; + i_int_ADE = NULL; + i_ext_ADE = NULL; + + v_Lor_ADE = NULL; + i_Lor_ADE = NULL; + + m_curr_Lor_ADE_On = NULL; + m_curr_Lor_ADE_On = NULL; +} + +Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op, Operator_Ext_LorentzMaterial* op_ext) : Operator_Ext_Dispersive(op,op_ext) +{ + v_int_ADE = NULL; + v_ext_ADE = NULL; + i_int_ADE = NULL; + i_ext_ADE = NULL; + + v_Lor_ADE = NULL; + i_Lor_ADE = NULL; + + m_curr_Lor_ADE_On = NULL; + m_curr_Lor_ADE_On = NULL; +} + +Operator_Ext_LorentzMaterial::~Operator_Ext_LorentzMaterial() +{ + for (int i=0;i<m_Order;++i) + { + for (int n=0; n<3; ++n) + { + if (m_volt_ADE_On[i]) + { + delete[] v_int_ADE[i][n]; + delete[] v_ext_ADE[i][n]; + } + if (m_curr_ADE_On[i]) + { + delete[] i_int_ADE[i][n]; + delete[] i_ext_ADE[i][n]; + } + if (m_volt_Lor_ADE_On[i]) + delete[] v_Lor_ADE[i][n]; + if (m_curr_Lor_ADE_On[i]) + delete[] i_Lor_ADE[i][n]; + } + if (m_volt_ADE_On[i]) + { + delete[] v_int_ADE[i]; + delete[] v_ext_ADE[i]; + } + if (m_curr_ADE_On[i]) + { + delete[] i_int_ADE[i]; + delete[] i_ext_ADE[i]; + } + if (m_volt_Lor_ADE_On[i]) + delete[] v_Lor_ADE[i]; + if (m_curr_Lor_ADE_On[i]) + delete[] i_Lor_ADE[i]; + } + delete[] v_int_ADE; + delete[] v_ext_ADE; + delete[] i_int_ADE; + delete[] i_ext_ADE; + v_int_ADE = NULL; + v_ext_ADE = NULL; + i_int_ADE = NULL; + i_ext_ADE = NULL; + + delete[] v_Lor_ADE; + delete[] i_Lor_ADE; + v_Lor_ADE = NULL; + i_Lor_ADE = NULL; + + delete[] m_curr_Lor_ADE_On; + delete[] m_volt_Lor_ADE_On; + m_curr_Lor_ADE_On = NULL; + m_curr_Lor_ADE_On = NULL; +} + +Operator_Extension* Operator_Ext_LorentzMaterial::Clone(Operator* op) +{ + if (dynamic_cast<Operator_Ext_LorentzMaterial*>(this)==NULL) + return NULL; + return new Operator_Ext_LorentzMaterial(op, this); +} + +bool Operator_Ext_LorentzMaterial::BuildExtension() +{ + double dT = m_Op->GetTimestep(); + unsigned int pos[] = {0,0,0}; + double coord[3]; + unsigned int numLines[3] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)}; + CSPropLorentzMaterial* mat = NULL; + CSPropDebyeMaterial* debye_mat = NULL; + + bool warn_once = true; + + bool b_pos_on; + vector<unsigned int> v_pos[3]; + + // drude material parameter + double w_plasma,t_relax; + double L_D[3], C_D[3]; + double R_D[3], G_D[3]; + vector<double> v_int[3]; + vector<double> v_ext[3]; + vector<double> i_int[3]; + vector<double> i_ext[3]; + + //additional Dorentz material parameter + double w_Lor_Pol; + double C_L[3]; + double L_L[3]; + vector<double> v_Lor[3]; + vector<double> i_Lor[3]; + + m_Order = 0; + vector<CSProperties*> LD_props = m_Op->CSX->GetPropertyByType(CSProperties::LORENTZMATERIAL); + for (size_t n=0;n<LD_props.size();++n) + { + CSPropLorentzMaterial* LorMat = dynamic_cast<CSPropLorentzMaterial*>(LD_props.at(n)); + if (LorMat==NULL) + return false; //sanity check, this should not happen + if (LorMat->GetDispersionOrder()>m_Order) + m_Order=LorMat->GetDispersionOrder(); + } + LD_props = m_Op->CSX->GetPropertyByType(CSProperties::DEBYEMATERIAL); + for (size_t n=0;n<LD_props.size();++n) + { + CSPropDebyeMaterial* DebyeMat = dynamic_cast<CSPropDebyeMaterial*>(LD_props.at(n)); + if (DebyeMat==NULL) + return false; //sanity check, this should not happen + if (DebyeMat->GetDispersionOrder()>m_Order) + m_Order=DebyeMat->GetDispersionOrder(); + } + + m_LM_pos = new unsigned int**[m_Order]; + + m_volt_ADE_On = new bool[m_Order]; + m_curr_ADE_On = new bool[m_Order]; + m_volt_Lor_ADE_On = new bool[m_Order]; + m_curr_Lor_ADE_On = new bool[m_Order]; + + v_int_ADE = new FDTD_FLOAT**[m_Order]; + v_ext_ADE = new FDTD_FLOAT**[m_Order]; + i_int_ADE = new FDTD_FLOAT**[m_Order]; + i_ext_ADE = new FDTD_FLOAT**[m_Order]; + + v_Lor_ADE = new FDTD_FLOAT**[m_Order]; + i_Lor_ADE = new FDTD_FLOAT**[m_Order]; + + for (int order=0;order<m_Order;++order) + { + m_volt_ADE_On[order]=false; + m_curr_ADE_On[order]=false; + + m_volt_Lor_ADE_On[order]=false; + m_curr_Lor_ADE_On[order]=false; + + for (int n=0;n<3;++n) + { + v_pos[n].clear(); + + v_int[n].clear(); + v_ext[n].clear(); + i_int[n].clear(); + i_ext[n].clear(); + + v_Lor[n].clear(); + i_Lor[n].clear(); + } + + for (pos[0]=0; pos[0]<numLines[0]; ++pos[0]) + { + for (pos[1]=0; pos[1]<numLines[1]; ++pos[1]) + { + vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL)); + for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) + { + unsigned int index = m_Op->MainOp->SetPos(pos[0],pos[1],pos[2]); + //calc epsilon lorentz material + b_pos_on = false; + for (int n=0; n<3; ++n) + { + L_D[n]=0; + R_D[n]=0; + C_L[n]=0; + if (m_Op->GetYeeCoords(n,pos,coord,false)==false) + continue; + if (m_CC_R0_included && (n==2) && (pos[0]==0)) + coord[1] = m_Op->GetDiscLine(1,0); + + if (m_Op->GetVI(n,pos[0],pos[1],pos[2])==0) + continue; + +// CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), true); + CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, true); + if ((mat = prop->ToLorentzMaterial())) + { + w_plasma = mat->GetEpsPlasmaFreqWeighted(order,n,coord) * 2 * PI; + if ((w_plasma>0) && (m_Op->EC_C[n][index]>0)) + { + b_pos_on = true; + m_volt_ADE_On[order] = true; + L_D[n] = 1/(w_plasma*w_plasma*m_Op->EC_C[n][index]); + } + t_relax = mat->GetEpsRelaxTimeWeighted(order,n,coord); + if ((t_relax>0) && m_volt_ADE_On[order]) + { + R_D[n] = L_D[n]/t_relax; + } + w_Lor_Pol = mat->GetEpsLorPoleFreqWeighted(order,n,coord) * 2 * PI; + if ((w_Lor_Pol>0) && (L_D[n]>0)) + { + m_volt_Lor_ADE_On[order] = true; + C_L[n] = 1/(w_Lor_Pol*w_Lor_Pol*L_D[n]); + } + } + if ((debye_mat = prop->ToDebyeMaterial())) + { + C_L[n] = 8.85418781762e-12*debye_mat->GetEpsDeltaWeighted(order,n,coord) * m_Op->GetEdgeArea(n, pos) / m_Op->GetEdgeLength(n,pos); + t_relax = debye_mat->GetEpsRelaxTimeWeighted(order,n,coord); + if ((t_relax<2.0*dT) && warn_once) + { + warn_once = false; + cerr << "Operator_Ext_LorentzMaterial::BuildExtension(): Warning, debye relaxation time is to small, skipping..." << endl; + } + if ((C_L[n]>0) && (t_relax>0) && (t_relax>2.0*dT)) + { + R_D[n] = t_relax/C_L[n]; + b_pos_on = true; + m_volt_ADE_On[order] = true; + m_volt_Lor_ADE_On[order] = true; + } + } + } + + for (int n=0; n<3; ++n) + { + C_D[n]=0; + G_D[n]=0; + L_L[n]=0; + if (m_Op->GetYeeCoords(n,pos,coord,true)==false) + continue; + if (m_Op->GetIV(n,pos[0],pos[1],pos[2])==0) + continue; + +// CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), true); + CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, true); + if ((mat = prop->ToLorentzMaterial())) + { + w_plasma = mat->GetMuePlasmaFreqWeighted(order,n,coord) * 2 * PI; + if ((w_plasma>0) && (m_Op->EC_L[n][index]>0)) + { + b_pos_on = true; + m_curr_ADE_On[order] = true; + C_D[n] = 1/(w_plasma*w_plasma*m_Op->EC_L[n][index]); + } + t_relax = mat->GetMueRelaxTimeWeighted(order,n,coord); + if ((t_relax>0) && m_curr_ADE_On[order]) + { + G_D[n] = C_D[n]/t_relax; + } + w_Lor_Pol = mat->GetMueLorPoleFreqWeighted(order,n,coord) * 2 * PI; + if ((w_Lor_Pol>0) && (C_D[n]>0)) + { + m_curr_Lor_ADE_On[order] = true; + L_L[n] = 1/(w_Lor_Pol*w_Lor_Pol*C_D[n]); + } + } + } + + if (b_pos_on) //this position has active drude material + { + for (unsigned int n=0; n<3; ++n) + { + v_pos[n].push_back(pos[n]); + if (L_D[n]>0) + { + v_int[n].push_back((2.0*L_D[n]-dT*R_D[n])/(2.0*L_D[n]+dT*R_D[n])); + // check for r==0 in clyindrical coords and get special VI cooefficient + if (m_CC_R0_included && n==2 && pos[0]==0) + v_ext[n].push_back(dT/(L_D[n]+dT*R_D[n]/2.0)*m_Op_Cyl->m_Cyl_Ext->vi_R0[pos[2]]); + else + v_ext[n].push_back(dT/(L_D[n]+dT*R_D[n]/2.0)*m_Op->GetVI(n,pos[0],pos[1],pos[2])); + } + else if ((R_D[n]>0) && (C_L[n]>0)) + { + v_int[n].push_back((2.0*dT-R_D[n]*C_L[n])/(C_L[n]*R_D[n])); + v_ext[n].push_back(2.0/R_D[n]*m_Op->GetVI(n,pos[0],pos[1],pos[2])); + } + else + { + v_int[n].push_back(1); + v_ext[n].push_back(0); + } + if (C_D[n]>0) + { + i_int[n].push_back((2.0*C_D[n]-dT*G_D[n])/(2.0*C_D[n]+dT*G_D[n])); + i_ext[n].push_back(dT/(C_D[n]+dT*G_D[n]/2.0)*m_Op->GetIV(n,pos[0],pos[1],pos[2])); + } + else + { + i_int[n].push_back(1); + i_ext[n].push_back(0); + } + if (C_L[n]>0) + v_Lor[n].push_back(dT/C_L[n]/m_Op->GetVI(n,pos[0],pos[1],pos[2])); + else + v_Lor[n].push_back(0); + if (L_L[n]>0) + i_Lor[n].push_back(dT/L_L[n]/m_Op->GetIV(n,pos[0],pos[1],pos[2])); + else + i_Lor[n].push_back(0); + } + } + } + } + } + + //copy all vectors into the array's + m_LM_Count.push_back(v_pos[0].size()); + + m_LM_pos[order] = new unsigned int*[3]; + + if (m_volt_ADE_On[order]) + { + v_int_ADE[order] = new FDTD_FLOAT*[3]; + v_ext_ADE[order] = new FDTD_FLOAT*[3]; + } + else + { + v_int_ADE[order] = NULL; + v_ext_ADE[order] = NULL; + } + + if (m_curr_ADE_On[order]) + { + i_int_ADE[order] = new FDTD_FLOAT*[3]; + i_ext_ADE[order] = new FDTD_FLOAT*[3]; + } + else + { + i_int_ADE[order] = NULL; + i_ext_ADE[order] = NULL; + } + + if (m_volt_Lor_ADE_On[order]) + v_Lor_ADE[order] = new FDTD_FLOAT*[3]; + else + v_Lor_ADE[order] = NULL; + + if (m_curr_Lor_ADE_On[order]) + i_Lor_ADE[order] = new FDTD_FLOAT*[3]; + else + i_Lor_ADE[order] = NULL; + + for (int n=0; n<3; ++n) + { + m_LM_pos[order][n] = new unsigned int[m_LM_Count.at(order)]; + for (unsigned int i=0; i<m_LM_Count.at(order); ++i) + m_LM_pos[order][n][i] = v_pos[n].at(i); + if (m_volt_ADE_On[order]) + { + v_int_ADE[order][n] = new FDTD_FLOAT[m_LM_Count.at(order)]; + v_ext_ADE[order][n] = new FDTD_FLOAT[m_LM_Count.at(order)]; + + for (unsigned int i=0; i<m_LM_Count.at(order); ++i) + { + v_int_ADE[order][n][i] = v_int[n].at(i); + v_ext_ADE[order][n][i] = v_ext[n].at(i); + } + } + if (m_curr_ADE_On[order]) + { + i_int_ADE[order][n] = new FDTD_FLOAT[m_LM_Count.at(order)]; + i_ext_ADE[order][n] = new FDTD_FLOAT[m_LM_Count.at(order)]; + + for (unsigned int i=0; i<m_LM_Count.at(order); ++i) + { + i_int_ADE[order][n][i] = i_int[n].at(i); + i_ext_ADE[order][n][i] = i_ext[n].at(i); + } + } + + if (m_volt_Lor_ADE_On[order]) + { + v_Lor_ADE[order][n] = new FDTD_FLOAT[m_LM_Count.at(order)]; + for (unsigned int i=0; i<m_LM_Count.at(order); ++i) + v_Lor_ADE[order][n][i] = v_Lor[n].at(i); + } + if (m_curr_Lor_ADE_On[order]) + { + i_Lor_ADE[order][n] = new FDTD_FLOAT[m_LM_Count.at(order)]; + for (unsigned int i=0; i<m_LM_Count.at(order); ++i) + i_Lor_ADE[order][n][i] = i_Lor[n].at(i); + } + } + } + + return true; +} + +Engine_Extension* Operator_Ext_LorentzMaterial::CreateEngineExtention() +{ + Engine_Ext_LorentzMaterial* eng_ext_lor = new Engine_Ext_LorentzMaterial(this); + return eng_ext_lor; +} + +void Operator_Ext_LorentzMaterial::ShowStat(ostream &ostr) const +{ + Operator_Extension::ShowStat(ostr); + string On_Off[2] = {"Off", "On"}; + ostr << " Max. Dispersion Order N = " << m_Order << endl; + for (int i=0;i<m_Order;++i) + { + ostr << " N=" << i << ":\t Active cells\t\t: " << m_LM_Count.at(i) << endl; + ostr << " N=" << i << ":\t Voltage ADE is \t: " << On_Off[m_volt_ADE_On[i]] << endl; + ostr << " N=" << i << ":\t Voltage Lor-ADE is \t: " << On_Off[m_volt_Lor_ADE_On[i]] << endl; + ostr << " N=" << i << ":\t Current ADE is \t: " << On_Off[m_curr_ADE_On[i]] << endl; + ostr << " N=" << i << ":\t Current Lor-ADE is \t: " << On_Off[m_curr_Lor_ADE_On[i]] << endl; + } +} diff --git a/openEMS/FDTD/extensions/operator_ext_lorentzmaterial.h b/openEMS/FDTD/extensions/operator_ext_lorentzmaterial.h new file mode 100644 index 0000000..d9e600b --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_lorentzmaterial.h @@ -0,0 +1,62 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef OPERATOR_EXT_LORENTZMATERIAL_H +#define OPERATOR_EXT_LORENTZMATERIAL_H + +#include "FDTD/operator.h" +#include "operator_ext_dispersive.h" + +class Operator_Ext_LorentzMaterial : public Operator_Ext_Dispersive +{ + friend class Engine_Ext_LorentzMaterial; +public: + Operator_Ext_LorentzMaterial(Operator* op); + virtual ~Operator_Ext_LorentzMaterial(); + + virtual Operator_Extension* Clone(Operator* op); + + virtual bool BuildExtension(); + + virtual Engine_Extension* CreateEngineExtention(); + + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return true;} + virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return true;} + virtual bool IsMPISave() const {return true;} + + virtual string GetExtensionName() const {return string("Drude/Lorentz Dispersive Material Extension");} + + virtual void ShowStat(ostream &ostr) const; + +protected: + //! Copy constructor + Operator_Ext_LorentzMaterial(Operator* op, Operator_Ext_LorentzMaterial* op_ext); + + //ADE update coefficients, array setup: coeff[N_order][direction][mesh_pos_index] + FDTD_FLOAT ***v_int_ADE; + FDTD_FLOAT ***v_ext_ADE; + FDTD_FLOAT ***i_int_ADE; + FDTD_FLOAT ***i_ext_ADE; + + bool *m_curr_Lor_ADE_On; + bool *m_volt_Lor_ADE_On; + + FDTD_FLOAT ***v_Lor_ADE; + FDTD_FLOAT ***i_Lor_ADE; +}; + +#endif // OPERATOR_EXT_LORENTZMATERIAL_H diff --git a/openEMS/FDTD/extensions/operator_ext_mur_abc.cpp b/openEMS/FDTD/extensions/operator_ext_mur_abc.cpp new file mode 100644 index 0000000..1df8583 --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_mur_abc.cpp @@ -0,0 +1,210 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "operator_ext_mur_abc.h" +#include "engine_ext_mur_abc.h" + +#include "tools/array_ops.h" + +#include "CSPropMaterial.h" + +Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op) : Operator_Extension(op) +{ + Initialize(); +} + +Operator_Ext_Mur_ABC::~Operator_Ext_Mur_ABC() +{ + Delete2DArray(m_Mur_Coeff_nyP,m_numLines); + m_Mur_Coeff_nyP = NULL; + Delete2DArray(m_Mur_Coeff_nyPP,m_numLines); + m_Mur_Coeff_nyPP = NULL; +} + +Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op, Operator_Ext_Mur_ABC* op_ext) : Operator_Extension(op, op_ext) +{ + Initialize(); + m_v_phase = op_ext->m_v_phase; + SetDirection(op_ext->m_ny,op_ext->m_top); +} + +Operator_Extension* Operator_Ext_Mur_ABC::Clone(Operator* op) +{ + if (dynamic_cast<Operator_Ext_Mur_ABC*>(this)==NULL) + return NULL; + return new Operator_Ext_Mur_ABC(op, this); +} + +bool Operator_Ext_Mur_ABC::IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const +{ + if ((m_ny==0) && (!m_top) && (R0_included || closedAlpha)) + return false; + if ((m_ny==1) && (closedAlpha)) + return false; + return true; +} + +bool Operator_Ext_Mur_ABC::IsCylindricalMultiGridSave(bool child) const +{ + if (m_ny==2) //always allow in z-direction + return true; + if ((m_ny==0) && (m_top) && (!child)) //if top r-direction and is not a child grid allow Mur... + return true; + //in all other cases this ABC is not save to use in CylindricalMultiGrid + return false; +} + +void Operator_Ext_Mur_ABC::Initialize() +{ + m_ny = -1; + m_nyP = -1; + m_nyPP = -1; + m_LineNr = 0; + m_LineNr_Shift = 0; + + m_v_phase = 0.0; + + m_Mur_Coeff_nyP = NULL; + m_Mur_Coeff_nyPP = NULL; + + m_numLines[0]=0; + m_numLines[1]=0; +} + +void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny) +{ + if ((ny<0) || (ny>2)) + return; + + Delete2DArray(m_Mur_Coeff_nyP,m_numLines); + Delete2DArray(m_Mur_Coeff_nyPP,m_numLines); + + m_ny = ny; + m_top = top_ny; + m_nyP = (ny+1)%3; + m_nyPP = (ny+2)%3; + if (!top_ny) + { + m_LineNr = 0; + m_LineNr_Shift = 1; + } + else + { + m_LineNr = m_Op->GetNumberOfLines(m_ny,true)-1; + m_LineNr_Shift = m_Op->GetNumberOfLines(m_ny,true) - 2; + } + + m_numLines[0] = m_Op->GetNumberOfLines(m_nyP,true); + m_numLines[1] = m_Op->GetNumberOfLines(m_nyPP,true); + + m_Mur_Coeff_nyP = Create2DArray<FDTD_FLOAT>(m_numLines); + m_Mur_Coeff_nyPP = Create2DArray<FDTD_FLOAT>(m_numLines); + +} + +bool Operator_Ext_Mur_ABC::BuildExtension() +{ + if (m_ny<0) + { + cerr << "Operator_Ext_Mur_ABC::BuildExtension: Warning, Extension not initialized! Use SetDirection!! Abort build!!" << endl; + return false; + } + double dT = m_Op->GetTimestep(); + unsigned int pos[] = {0,0,0}; + pos[m_ny] = m_LineNr; + double delta = fabs(m_Op->GetEdgeLength(m_ny,pos)); + double coord[] = {0,0,0}; + coord[0] = m_Op->GetDiscLine(0,pos[0]); + coord[1] = m_Op->GetDiscLine(1,pos[1]); + coord[2] = m_Op->GetDiscLine(2,pos[2]); + + double eps,mue; + double c0t; + + if (m_LineNr==0) + coord[m_ny] = m_Op->GetDiscLine(m_ny,pos[m_ny]) + delta/2 / m_Op->GetGridDelta(); + else + coord[m_ny] = m_Op->GetDiscLine(m_ny,pos[m_ny]) - delta/2 / m_Op->GetGridDelta(); + + int posBB[3]; + posBB[m_ny] =pos[m_ny]; + posBB[m_nyPP]=-1; + + for (pos[m_nyP]=0; pos[m_nyP]<m_numLines[0]; ++pos[m_nyP]) + { + posBB[m_nyP]=pos[m_nyP]; + vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(posBB[0], posBB[1], posBB[2], CSProperties::MATERIAL); + coord[m_nyP] = m_Op->GetDiscLine(m_nyP,pos[m_nyP]); + for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) + { + coord[m_nyPP] = m_Op->GetDiscLine(m_nyPP,pos[m_nyPP]); +// CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, CSProperties::MATERIAL, false); + CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, false); + if (prop) + { + CSPropMaterial* mat = prop->ToMaterial(); + + //nP + eps = mat->GetEpsilonWeighted(m_nyP,coord); + mue = mat->GetMueWeighted(m_nyP,coord); + if (m_v_phase>0.0) + c0t = m_v_phase * dT; + else + c0t = __C0__ * dT / sqrt(eps*mue); + m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta); + + //nPP + eps = mat->GetEpsilonWeighted(m_nyPP,coord); + mue = mat->GetMueWeighted(m_nyPP,coord); + if (m_v_phase>0.0) + c0t = m_v_phase * dT; + else + c0t = __C0__ * dT / sqrt(eps*mue); + m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta); + + } + else + { + if (m_v_phase>0.0) + c0t = m_v_phase * dT; + else + c0t = __C0__ / sqrt(m_Op->GetBackgroundEpsR()*m_Op->GetBackgroundMueR()) * dT; + m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta); + m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]]; + } +// cerr << m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] << " : " << m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] << endl; + } + } +// cerr << "Operator_Ext_Mur_ABC::BuildExtension(): " << m_ny << " @ " << m_LineNr << endl; + return true; +} + +Engine_Extension* Operator_Ext_Mur_ABC::CreateEngineExtention() +{ + Engine_Ext_Mur_ABC* eng_ext = new Engine_Ext_Mur_ABC(this); + return eng_ext; +} + + +void Operator_Ext_Mur_ABC::ShowStat(ostream &ostr) const +{ + Operator_Extension::ShowStat(ostr); + string XYZ[3] = {"x","y","z"}; + ostr << " Active direction\t: " << XYZ[m_ny] << " at line: " << m_LineNr << endl; + if (m_v_phase>0.0) + ostr << " Used phase velocity\t: " << m_v_phase << " (" << m_v_phase/__C0__ << " * c_0)" <<endl; +} diff --git a/openEMS/FDTD/extensions/operator_ext_mur_abc.h b/openEMS/FDTD/extensions/operator_ext_mur_abc.h new file mode 100644 index 0000000..ddc43c1 --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_mur_abc.h @@ -0,0 +1,68 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef OPERATOR_EXT_MUR_ABC_H +#define OPERATOR_EXT_MUR_ABC_H + +#include "FDTD/operator.h" +#include "operator_extension.h" + +class Operator_Ext_Mur_ABC : public Operator_Extension +{ + friend class Engine_Ext_Mur_ABC; +public: + Operator_Ext_Mur_ABC(Operator* op); + ~Operator_Ext_Mur_ABC(); + + virtual Operator_Extension* Clone(Operator* op); + + //! Define the direction of this ABC: \a ny=0,1,2 -> x,y,z and if at bottom_ny -> e.g. x=0 or x=end + void SetDirection(int ny, bool top_ny); + + //! Set (override) the expected phase velocity of the incoming wave + void SetPhaseVelocity(double c_phase) {m_v_phase=c_phase;}; + + virtual bool BuildExtension(); + + virtual Engine_Extension* CreateEngineExtention(); + + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const; + virtual bool IsCylindricalMultiGridSave(bool child) const; + virtual bool IsMPISave() const {return true;} + + virtual string GetExtensionName() const {return string("Mur ABC Extension");} + + virtual void ShowStat(ostream &ostr) const; + +protected: + Operator_Ext_Mur_ABC(Operator* op, Operator_Ext_Mur_ABC* op_ext); + void Initialize(); + int m_ny; + int m_nyP,m_nyPP; + bool m_top; + unsigned int m_LineNr; + int m_LineNr_Shift; + + double m_v_phase; + + unsigned int m_numLines[2]; + + FDTD_FLOAT** m_Mur_Coeff_nyP; + FDTD_FLOAT** m_Mur_Coeff_nyPP; +}; + +#endif // OPERATOR_EXT_MUR_ABC_H diff --git a/openEMS/FDTD/extensions/operator_ext_steadystate.cpp b/openEMS/FDTD/extensions/operator_ext_steadystate.cpp new file mode 100644 index 0000000..b8395a8 --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_steadystate.cpp @@ -0,0 +1,104 @@ +/* +* Copyright (C) 2015 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 <http://www.gnu.org/licenses/>. +*/ + +#include "operator_ext_steadystate.h" +#include "engine_ext_steadystate.h" + +Operator_Ext_SteadyState::Operator_Ext_SteadyState(Operator* op, double period): Operator_Extension(op) +{ + this->Reset(); + m_T_period = period; +} + +Operator_Ext_SteadyState::Operator_Ext_SteadyState(Operator* op, Operator_Ext_SteadyState* op_ext): Operator_Extension(op, op_ext) +{ + this->Reset(); + m_T_period = op_ext->m_T_period; +} + +Operator_Ext_SteadyState::~Operator_Ext_SteadyState() +{ +} + +Operator_Extension* Operator_Ext_SteadyState::Clone(Operator* op) +{ + //disable cloning, only the main operator needs to have a steady state detection + UNUSED(op); + return NULL; +} + +bool Operator_Ext_SteadyState::BuildExtension() +{ + double dT = m_Op->GetTimestep(); + m_TS_period = round(m_T_period/dT); + m_T_period = m_TS_period*dT; + return true; +} + +void Operator_Ext_SteadyState::Reset() +{ + for (int n=0;n<3;++n) + { + m_E_probe_pos[n].clear(); + m_H_probe_pos[n].clear(); + } + m_E_probe_dir.clear(); + m_H_probe_dir.clear(); + m_T_period = 0; + m_TS_period = 0; +} + +bool Operator_Ext_SteadyState::Add_E_Probe(unsigned int pos[3], int dir) +{ + if ((dir<0) || (dir>2)) + return false; + for (int n=0;n<3;++n) + if (pos[n]>=m_Op->GetNumberOfLines(n)) + return false; + for (int n=0;n<3;++n) + m_E_probe_pos[n].push_back(pos[n]); + m_E_probe_dir.push_back(dir); + return true; +} + +bool Operator_Ext_SteadyState::Add_H_Probe(unsigned int pos[3], int dir) +{ + if ((dir<0) || (dir>2)) + return false; + for (int n=0;n<3;++n) + if (pos[n]>=m_Op->GetNumberOfLines(n)) + return false; + for (int n=0;n<3;++n) + m_H_probe_pos[n].push_back(pos[n]); + m_H_probe_dir.push_back(dir); + return true; +} + +Engine_Extension* Operator_Ext_SteadyState::CreateEngineExtention() +{ + m_Eng_Ext = new Engine_Ext_SteadyState(this); + return m_Eng_Ext; +} + +void Operator_Ext_SteadyState::ShowStat(ostream &ostr) const +{ + Operator_Extension::ShowStat(ostr); + cout << "Period time (s): " << m_T_period << "\t Period TS: " << m_TS_period << endl; + cout << "Number of E probes\t: " << m_E_probe_dir.size() << endl; + cout << "Number of H probes\t: " << m_H_probe_dir.size() << endl; +} + diff --git a/openEMS/FDTD/extensions/operator_ext_steadystate.h b/openEMS/FDTD/extensions/operator_ext_steadystate.h new file mode 100644 index 0000000..b735ca3 --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_steadystate.h @@ -0,0 +1,61 @@ +/* +* Copyright (C) 2015 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef OPERATOR_EXT_STEADYSTATE_H +#define OPERATOR_EXT_STEADYSTATE_H + +#include "operator_extension.h" +#include "FDTD/operator.h" + +class Engine_Ext_SteadyState; + +class Operator_Ext_SteadyState : public Operator_Extension +{ + friend class Engine_Ext_SteadyState; +public: + Operator_Ext_SteadyState(Operator* op, double period); + virtual ~Operator_Ext_SteadyState(); + + virtual Operator_Extension* Clone(Operator* op); + + virtual bool BuildExtension(); + virtual Engine_Extension* CreateEngineExtention(); + + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return true;} + virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return true;} + virtual bool IsMPISave() const {return true;} + + virtual string GetExtensionName() const {return string("Steady-State Detection Extension");} + + virtual void ShowStat(ostream &ostr) const; + + virtual void Reset(); + + bool Add_E_Probe(unsigned int pos[3], int dir); + bool Add_H_Probe(unsigned int pos[3], int dir); + +protected: + Operator_Ext_SteadyState(Operator* op, Operator_Ext_SteadyState* op_ext); + double m_T_period; + unsigned int m_TS_period; + vector<unsigned int> m_E_probe_pos[3]; + vector<unsigned int> m_E_probe_dir; + vector<unsigned int> m_H_probe_pos[3]; + vector<unsigned int> m_H_probe_dir; +}; + +#endif // OPERATOR_EXT_STEADYSTATE_H diff --git a/openEMS/FDTD/extensions/operator_ext_tfsf.cpp b/openEMS/FDTD/extensions/operator_ext_tfsf.cpp new file mode 100644 index 0000000..663dc1f --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_tfsf.cpp @@ -0,0 +1,429 @@ +/* +* Copyright (C) 2012 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 <http://www.gnu.org/licenses/>. +*/ + +#include "operator_ext_tfsf.h" +#include "engine_ext_tfsf.h" +#include <cmath> + +#include "CSPrimBox.h" +#include "CSPropExcitation.h" + +Operator_Ext_TFSF::Operator_Ext_TFSF(Operator* op) : Operator_Extension(op) +{ + Init(); +} + +Operator_Ext_TFSF::~Operator_Ext_TFSF() +{ + Reset(); +} + +void Operator_Ext_TFSF::Init() +{ + for (int n=0;n<3;++n) + for (int l=0;l<2;++l) + for (int c=0;c<2;++c) + { + m_VoltDelay[n][l][c]=NULL; + m_VoltDelayDelta[n][l][c]=NULL; + m_VoltAmp[n][l][c]=NULL; + m_CurrDelay[n][l][c]=NULL; + m_CurrDelayDelta[n][l][c]=NULL; + m_CurrAmp[n][l][c]=NULL; + } + + m_Frequency = 0.0; + m_PhVel = __C0__; + Operator_Extension::Init(); +} + +void Operator_Ext_TFSF::Reset() +{ + for (int n=0;n<3;++n) + for (int l=0;l<2;++l) + for (int c=0;c<2;++c) + { + delete[] m_VoltDelay[n][l][c]; + m_VoltDelay[n][l][c]=NULL; + delete[] m_VoltDelayDelta[n][l][c]; + m_VoltDelayDelta[n][l][c]=NULL; + delete[] m_VoltAmp[n][l][c]; + m_VoltAmp[n][l][c]=NULL; + delete[] m_CurrDelay[n][l][c]; + m_CurrDelay[n][l][c]=NULL; + delete[] m_CurrDelayDelta[n][l][c]; + m_CurrDelayDelta[n][l][c]=NULL; + delete[] m_CurrAmp[n][l][c]; + m_CurrAmp[n][l][c]=NULL; + } + Operator_Extension::Reset(); +} + +Operator_Extension* Operator_Ext_TFSF::Clone(Operator* op) +{ + UNUSED(op); + return NULL; +} + +bool Operator_Ext_TFSF::BuildExtension() +{ + m_Exc = m_Op->GetExcitationSignal(); + double dT = m_Op->GetTimestep(); + if (dT==0) + return false; + if (m_Exc==0) + return false; + + Reset(); + ContinuousStructure* CSX = m_Op->GetGeometryCSX(); + + vector<CSProperties*> vec_prop = CSX->GetPropertyByType(CSProperties::EXCITATION); + + if (vec_prop.size()==0) + { + cerr << "Operator_Ext_TFSF::BuildExtension: Warning, no excitation properties found" << endl; + SetActive(false); + return false; + } + + double ref_index = sqrt(m_Op->GetBackgroundEpsR()*m_Op->GetBackgroundMueR()); + m_PhVel = __C0__/ref_index; + CSPropExcitation* elec=NULL; + CSProperties* prop=NULL; + CSPrimitives* prim=NULL; + CSPrimBox* box=NULL; + for (size_t p=0; p<vec_prop.size(); ++p) + { + prop = vec_prop.at(p); + elec = prop->ToExcitation(); + if (elec->GetExcitType()!=10) + continue; + if (prop->GetQtyPrimitives()!=1) + { + cerr << "Operator_Ext_TFSF::BuildExtension: Warning, plane wave excitation found with more (or less) than one primitive, skipping..." << endl; + continue; + } + prim = prop->GetPrimitive(0); + if (prim->GetType()!=CSPrimitives::BOX) + { + cerr << "Operator_Ext_TFSF::BuildExtension: Warning, plane wave excitation found with false non-Box primitive, skipping..." << endl; + continue; + } + box = prim->ToBox(); + if (box==NULL) //sanity check, should not happen! + { + SetActive(false); + return false; + } + + // found a plane-wave excite with exactly one box + for (int n=0;n<3;++n) + m_PropDir[n] = elec->GetPropagationDir(n); + double dir_norm = sqrt(m_PropDir[0]*m_PropDir[0]+m_PropDir[1]*m_PropDir[1]+m_PropDir[2]*m_PropDir[2]); + if (dir_norm==0) + { + cerr << "Operator_Ext_TFSF::BuildExtension: Warning, plane wave direction is zero, skipping..." << endl; + SetActive(false); + return false; + } + + //make it a unit vector + m_PropDir[0]/=dir_norm;m_PropDir[1]/=dir_norm;m_PropDir[2]/=dir_norm; + + if (m_Op->SnapBox2Mesh(box->GetStartCoord()->GetNativeCoords(), box->GetStopCoord()->GetNativeCoords(), m_Start, m_Stop)!=3) + { + cerr << "Operator_Ext_TFSF::BuildExtension: Warning, plane wave box dimension is invalid, skipping..." << endl; + SetActive(false); + return false; + } + + m_Frequency = elec->GetFrequency(); +// if (m_Frequency<=0) +// m_Frequency = m_Op->GetExcitationSignal()->GetFrequencyOfInterest(); + if (m_Frequency<=0) + m_PhVel=__C0__/ref_index; + else + m_PhVel=m_Op->CalcNumericPhaseVelocity(m_Start,m_Stop,m_PropDir,m_Frequency); + + if ((m_PhVel<0) || (m_PhVel>__C0__/ref_index) || isnan(m_PhVel)) + { + cerr << "Operator_Ext_TFSF::BuildExtension: Warning, invalid phase velocity found, resetting to c0! " << endl; + m_PhVel = __C0__/ref_index; + } + + double origin[3]; + unsigned int ui_origin[3]; + for (int n=0;n<3;++n) + { + m_E_Amp[n] = elec->GetExcitation(n); + m_numLines[n] = m_Stop[n]-m_Start[n]+1; + m_IncLow[n] = m_PropDir[n]>=0; + + if (m_Start[n]==0) + m_ActiveDir[n][0]=false; + else + m_ActiveDir[n][0]=true; + if (m_Stop[n]==m_Op->GetNumberOfLines(n,true)-1) + m_ActiveDir[n][1]=false; + else + m_ActiveDir[n][1]=true; + + if (m_IncLow[n]) + { + ui_origin[n] = m_Start[n]-1; + } + else + { + ui_origin[n] = m_Stop[n]+1; + } + origin[n] = m_Op->GetDiscLine(n,ui_origin[n]); + } + + double dotEk = (m_E_Amp[0]*m_PropDir[0] + m_E_Amp[1]*m_PropDir[1] + m_E_Amp[2]*m_PropDir[2]); + double angle = acos( dotEk / (m_E_Amp[0]*m_E_Amp[0] + m_E_Amp[1]*m_E_Amp[1] + m_E_Amp[2]*m_E_Amp[2]) ) / PI * 180; + + if (angle==0) + { + cerr << "Operator_Ext_TFSF::BuildExtension: Warning, plane wave direction and polarization is identical, skipping..." << endl; + SetActive(false); + return false; + } + if (angle!=90) + { + cerr << "Operator_Ext_TFSF::BuildExtension: Warning, angle between propagation direction and polarization is not 90deg, resetting E-polarization to : ("; + for (int n=0;n<3;++n) + m_E_Amp[n]-=m_PropDir[n]*dotEk; + cerr << m_E_Amp[0] << "," << m_E_Amp[1] << "," << m_E_Amp[2] << ")" << endl; + } + + int nP,nPP; + for (int n=0;n<3;++n) + { + nP = (n+1)%3; + nPP = (n+2)%3; + m_H_Amp[n] = m_PropDir[nP]*m_E_Amp[nPP] - m_PropDir[nPP]*m_E_Amp[nP]; + m_H_Amp[n] /= __Z0__*sqrt(m_Op->GetBackgroundMueR()/m_Op->GetBackgroundEpsR()); + } + + double coord[3]; + double unit = m_Op->GetGridDelta(); + double delay; + double dist; + unsigned int pos[3]; + unsigned int numP, ui_pos; + m_maxDelay = 0; + for (int n=0;n<3;++n) + { + nP = (n+1)%3; + nPP = (n+2)%3; + pos[n] = 0; //unused + pos[nP] = m_Start[nP]; + + numP = m_numLines[nP]*m_numLines[nPP]; + + if (!m_ActiveDir[n][0] && !m_ActiveDir[n][1]) + continue; + + for (int l=0;l<2;++l) + for (int c=0;c<2;++c) + { + if (m_ActiveDir[n][l]) + { + m_VoltDelay[n][l][c]=new unsigned int[numP]; + m_VoltDelayDelta[n][l][c]=new FDTD_FLOAT[numP]; + m_VoltAmp[n][l][c]=new FDTD_FLOAT[numP]; + + m_CurrDelay[n][l][c]=new unsigned int[numP]; + m_CurrDelayDelta[n][l][c]=new FDTD_FLOAT[numP]; + m_CurrAmp[n][l][c]=new FDTD_FLOAT[numP]; + } + } + + ui_pos = 0; + for (unsigned int i=0;i<m_numLines[nP];++i) + { + pos[nPP] = m_Start[nPP]; + for (unsigned int j=0;j<m_numLines[nPP];++j) + { + // current updates + pos[n] = m_Start[n]; + + if (m_ActiveDir[n][0]) + { + m_Op->GetYeeCoords(nP,pos,coord,false); + dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]); + delay = dist*unit/m_PhVel/dT; + m_maxDelay = max((unsigned int)delay,m_maxDelay); + m_CurrDelay[n][0][1][ui_pos] = floor(delay); + m_CurrDelayDelta[n][0][1][ui_pos] = delay - floor(delay); + m_CurrAmp[n][0][1][ui_pos] = m_E_Amp[nP]*m_Op->GetEdgeLength(nP,pos); + + m_Op->GetYeeCoords(nPP,pos,coord,false); + dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]); + delay = dist*unit/m_PhVel/dT; + m_maxDelay = max((unsigned int)delay,m_maxDelay); + m_CurrDelay[n][0][0][ui_pos] = floor(delay); + m_CurrDelayDelta[n][0][0][ui_pos] = delay - floor(delay); + m_CurrAmp[n][0][0][ui_pos] = m_E_Amp[nPP]*m_Op->GetEdgeLength(nPP,pos); + + --pos[n]; + m_CurrAmp[n][0][0][ui_pos]*=m_Op->GetIV(nP,pos); + m_CurrAmp[n][0][1][ui_pos]*=m_Op->GetIV(nPP,pos); + } + + if (m_ActiveDir[n][1]) + { + pos[n] = m_Stop[n]; + m_Op->GetYeeCoords(nP,pos,coord,false); + dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]); + delay = dist*unit/m_PhVel/dT; + m_maxDelay = max((unsigned int)delay,m_maxDelay); + m_CurrDelay[n][1][1][ui_pos] = floor(delay); + m_CurrDelayDelta[n][1][1][ui_pos] = delay - floor(delay); + m_CurrAmp[n][1][1][ui_pos] = m_E_Amp[nP]*m_Op->GetEdgeLength(nP,pos); + + m_Op->GetYeeCoords(nPP,pos,coord,false); + dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]); + delay = dist*unit/m_PhVel/dT; + m_maxDelay = max((unsigned int)delay,m_maxDelay); + m_CurrDelay[n][1][0][ui_pos] = floor(delay); + m_CurrDelayDelta[n][1][0][ui_pos] = delay - floor(delay); + m_CurrAmp[n][1][0][ui_pos] = m_E_Amp[nPP]*m_Op->GetEdgeLength(nPP,pos); + + m_CurrAmp[n][1][0][ui_pos]*=m_Op->GetIV(nP,pos); + m_CurrAmp[n][1][1][ui_pos]*=m_Op->GetIV(nPP,pos); + } + + if (m_ActiveDir[n][0]) + m_CurrAmp[n][0][0][ui_pos]*=-1; + if (m_ActiveDir[n][1]) + m_CurrAmp[n][1][1][ui_pos]*=-1; + + if (pos[nP]==m_Stop[nP]) + { + if (m_ActiveDir[n][0]) + m_CurrAmp[n][0][1][ui_pos]=0; + if (m_ActiveDir[n][1]) + m_CurrAmp[n][1][1][ui_pos]=0; + } + if (pos[nPP]==m_Stop[nPP]) + { + if (m_ActiveDir[n][0]) + m_CurrAmp[n][0][0][ui_pos]=0; + if (m_ActiveDir[n][1]) + m_CurrAmp[n][1][0][ui_pos]=0; + } + + // voltage updates + pos[n] = m_Start[n]-1; + if (m_ActiveDir[n][0]) + { + m_Op->GetYeeCoords(nP,pos,coord,true); + dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]); + delay = dist*unit/m_PhVel/dT + 1.0; + m_maxDelay = max((unsigned int)delay,m_maxDelay); + m_VoltDelay[n][0][1][ui_pos] = floor(delay); + m_VoltDelayDelta[n][0][1][ui_pos] = delay - floor(delay); + m_VoltAmp[n][0][1][ui_pos] = m_H_Amp[nP]*m_Op->GetEdgeLength(nP,pos,true); + + m_Op->GetYeeCoords(nPP,pos,coord,true); + dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]); + delay = dist*unit/m_PhVel/dT + 1.0; + m_maxDelay = max((unsigned int)delay,m_maxDelay); + m_VoltDelay[n][0][0][ui_pos] = floor(delay); + m_VoltDelayDelta[n][0][0][ui_pos] = delay - floor(delay); + m_VoltAmp[n][0][0][ui_pos] = m_H_Amp[nPP]*m_Op->GetEdgeLength(nPP,pos,true); + + ++pos[n]; + m_VoltAmp[n][0][0][ui_pos]*=m_Op->GetVI(nP,pos); + m_VoltAmp[n][0][1][ui_pos]*=m_Op->GetVI(nPP,pos); + } + + pos[n] = m_Stop[n]; + if (m_ActiveDir[n][1]) + { + m_Op->GetYeeCoords(nP,pos,coord,true); + dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]); + delay = dist*unit/m_PhVel/dT + 1.0; + m_maxDelay = max((unsigned int)delay,m_maxDelay); + m_VoltDelay[n][1][1][ui_pos] = floor(delay); + m_VoltDelayDelta[n][1][1][ui_pos] = delay - floor(delay); + m_VoltAmp[n][1][1][ui_pos] = m_H_Amp[nP]*m_Op->GetEdgeLength(nP,pos,true); + + m_Op->GetYeeCoords(nPP,pos,coord,true); + dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]); + delay = dist*unit/m_PhVel/dT + 1.0; + m_maxDelay = max((unsigned int)delay,m_maxDelay); + m_VoltDelay[n][1][0][ui_pos] = floor(delay); + m_VoltDelayDelta[n][1][0][ui_pos] = delay - floor(delay); + m_VoltAmp[n][1][0][ui_pos] = m_H_Amp[nPP]*m_Op->GetEdgeLength(nPP,pos,true); + + m_VoltAmp[n][1][0][ui_pos]*=m_Op->GetVI(nP,pos); + m_VoltAmp[n][1][1][ui_pos]*=m_Op->GetVI(nPP,pos); + } + + if (m_ActiveDir[n][1]) + m_VoltAmp[n][1][0][ui_pos]*=-1; + if (m_ActiveDir[n][0]) + m_VoltAmp[n][0][1][ui_pos]*=-1; + + if (pos[nP]==m_Stop[nP]) + { + if (m_ActiveDir[n][0]) + m_VoltAmp[n][0][0][ui_pos]=0; + if (m_ActiveDir[n][1]) + m_VoltAmp[n][1][0][ui_pos]=0; + } + if (pos[nPP]==m_Stop[nPP]) + { + if (m_ActiveDir[n][0]) + m_VoltAmp[n][0][1][ui_pos]=0; + if (m_ActiveDir[n][1]) + m_VoltAmp[n][1][1][ui_pos]=0; + } + + ++pos[nPP]; + ++ui_pos; + } + ++pos[nP]; + } + } + ++m_maxDelay; + return true; + } + SetActive(false); + return false; +} + +Engine_Extension* Operator_Ext_TFSF::CreateEngineExtention() +{ + return new Engine_Ext_TFSF(this); +} + +void Operator_Ext_TFSF::ShowStat(ostream &ostr) const +{ + Operator_Extension::ShowStat(ostr); + cout << "Active directions\t: " << "(" << m_ActiveDir[0][0] << "/" << m_ActiveDir[0][1] << ", " << m_ActiveDir[1][0] << "/" << m_ActiveDir[1][1] << ", " << m_ActiveDir[2][0] << "/" << m_ActiveDir[2][1] << ")" << endl; + cout << "Propagation direction\t: " << "(" << m_PropDir[0] << ", " << m_PropDir[1] << ", " << m_PropDir[2] << ")" << endl; + cout << "Rel. propagation speed\t: " << m_PhVel/__C0__ << "*c0 @ " << m_Frequency << " Hz" << endl; + cout << "E-field amplitude (V/m)\t: " << "(" << m_E_Amp[0] << ", " << m_E_Amp[1] << ", " << m_E_Amp[2] << ")" << endl; + cout << "H-field amplitude (A/m)\t: " << "(" << m_H_Amp[0] << ", " << m_H_Amp[1] << ", " << m_H_Amp[2] << ")" << endl; + cout << "Box Dimensions\t\t: " << m_numLines[0] << " x " << m_numLines[1] << " x " << m_numLines[2] << endl; + cout << "Max. Delay (TS)\t\t: " << m_maxDelay << endl; + int dirs = m_ActiveDir[0][0] + m_ActiveDir[0][1] + m_ActiveDir[1][0] + m_ActiveDir[1][1] + m_ActiveDir[2][0] + m_ActiveDir[2][1] ; + cout << "Memory usage (est.)\t: ~" << m_numLines[0] * m_numLines[1] * m_numLines[2] * dirs * 4 * 4 / 1024 << " kiB" << endl; +} diff --git a/openEMS/FDTD/extensions/operator_ext_tfsf.h b/openEMS/FDTD/extensions/operator_ext_tfsf.h new file mode 100644 index 0000000..e1e714e --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_tfsf.h @@ -0,0 +1,82 @@ +/* +* Copyright (C) 2012 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef OPERATOR_EXT_TFSF_H +#define OPERATOR_EXT_TFSF_H + +#include "operator_extension.h" +#include "FDTD/operator.h" +#include "tools/constants.h" + +class Excitation; + +class Operator_Ext_TFSF : public Operator_Extension +{ + friend class Engine_Ext_TFSF; +public: + Operator_Ext_TFSF(Operator* op); + ~Operator_Ext_TFSF(); + + virtual Operator_Extension* Clone(Operator* op); + + virtual bool BuildExtension(); + + virtual Engine_Extension* CreateEngineExtention(); + + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return false;} + virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return false;} + + // FIXME, this extension is not save to use with MPI + virtual bool IsMPISave() const {return false;} + + virtual string GetExtensionName() const {return string("Total-Field/Scattered-Field Extension");} + + virtual void ShowStat(ostream &ostr) const; + + virtual void Init(); + virtual void Reset(); + +protected: + Excitation* m_Exc; + + bool m_IncLow[3]; + bool m_ActiveDir[3][2]; // m_ActiveDir[direction][low/high] + unsigned int m_Start[3]; + unsigned int m_Stop[3]; + unsigned int m_numLines[3]; + + double m_PropDir[3]; + double m_E_Amp[3]; + double m_H_Amp[3]; + + double m_Frequency; + double m_PhVel; + + unsigned int m_maxDelay; + + // array setup [direction][low/high][component][ <mesh_position> ] + unsigned int* m_VoltDelay[3][2][2]; + FDTD_FLOAT* m_VoltDelayDelta[3][2][2]; + FDTD_FLOAT* m_VoltAmp[3][2][2]; + + unsigned int* m_CurrDelay[3][2][2]; + FDTD_FLOAT* m_CurrDelayDelta[3][2][2]; + FDTD_FLOAT* m_CurrAmp[3][2][2]; + +}; + +#endif // OPERATOR_EXT_TFSF_H diff --git a/openEMS/FDTD/extensions/operator_ext_upml.cpp b/openEMS/FDTD/extensions/operator_ext_upml.cpp new file mode 100644 index 0000000..fc46d85 --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_upml.cpp @@ -0,0 +1,478 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "operator_ext_upml.h" +#include "FDTD/operator_cylindermultigrid.h" +#include "engine_ext_upml.h" +#include "tools/array_ops.h" +#include "fparser.hh" + +using namespace std; + +Operator_Ext_UPML::Operator_Ext_UPML(Operator* op) : Operator_Extension(op) +{ + m_GradingFunction = new FunctionParser(); + //default grading function + SetGradingFunction(" -log(1e-6)*log(2.5)/(2*dl*Z*(pow(2.5,W/dl)-1)) * pow(2.5, D/dl) "); + + for (int n=0; n<6; ++n) + { + m_BC[n]=0; + m_Size[n]=0; + } + for (int n=0; n<3; ++n) + { + m_StartPos[n]=0; + m_numLines[n]=0; + } + + vv = NULL; + vvfo = NULL; + vvfn = NULL; + ii = NULL; + iifo = NULL; + iifn = NULL; +} + +Operator_Ext_UPML::~Operator_Ext_UPML() +{ + delete m_GradingFunction; + m_GradingFunction = NULL; + DeleteOp(); +} + +void Operator_Ext_UPML::SetBoundaryCondition(const int* BCs, const unsigned int size[6]) +{ + for (int n=0; n<6; ++n) + { + m_BC[n]=BCs[n]; + m_Size[n]=size[n]; + } +} + +void Operator_Ext_UPML::SetRange(const unsigned int start[3], const unsigned int stop[3]) +{ + for (int n=0; n<3; ++n) + { + m_StartPos[n]=start[n]; + m_numLines[n]=stop[n]-start[n]+1; + } +} + +bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsigned int ui_size[6], string gradFunc) +{ + int BC[6]={ui_BC[0],ui_BC[1],ui_BC[2],ui_BC[3],ui_BC[4],ui_BC[5]}; + unsigned int size[6]={ui_size[0],ui_size[1],ui_size[2],ui_size[3],ui_size[4],ui_size[5]}; + + //check if mesh is large enough to support the pml + for (int n=0; n<3; ++n) + if ( (size[2*n]*(BC[2*n]==3)+size[2*n+1]*(BC[2*n+1]==3)) >= op->GetNumberOfLines(n,true) ) + { + cerr << "Operator_Ext_UPML::Create_UPML: Warning: Not enough lines in direction: " << n << ", resetting to PEC" << endl; + BC[2*n]=0; + size[2*n]=0; + BC[2*n+1]=0; + size[2*n+1]=0; + } + + //check cylindrical coord compatiblility + Operator_Cylinder* op_cyl = dynamic_cast<Operator_Cylinder*>(op); + if (op_cyl) + { + if ((BC[0]==3) && (op_cyl->GetClosedAlpha() || op_cyl->GetR0Included())) + { + BC[0]=0; + size[0]=0; + cerr << "Operator_Ext_UPML::Create_UPML: Warning: An upml in r-min direction is not possible, resetting to PEC..." << endl; + } + if ( (BC[2]==3) && (op_cyl->GetClosedAlpha()) ) + { + BC[2]=0; + size[2]=0; + cerr << "Operator_Ext_UPML::Create_UPML: Warning: An upml in alpha-min direction is not possible, resetting to PEC..." << endl; + } + if ( (BC[3]==3) && (op_cyl->GetClosedAlpha()) ) + { + BC[3]=0; + size[3]=0; + cerr << "Operator_Ext_UPML::Create_UPML: Warning: An upml in alpha-max direction is not possible, resetting to PEC..." << endl; + } + } + + //check cylindrical coord compatiblility + if (dynamic_cast<Operator_CylinderMultiGrid*>(op)) + { + if (BC[2]==3) + { + BC[2]=0; + size[2]=0; + cerr << "Operator_Ext_UPML::Create_UPML: Warning: An upml in alpha direction is not possible for a cylindrical multi-grid, resetting to PEC..." << endl; + } + if (BC[3]==3) + { + BC[3]=0; + size[3]=0; + cerr << "Operator_Ext_UPML::Create_UPML: Warning: An upml in alpha direction is not possible for a cylindrical multi-grid, resetting to PEC..." << endl; + } + } + + + Operator_Ext_UPML* op_ext_upml=NULL; + unsigned int start[3]={0 ,0 ,0}; + unsigned int stop[3] ={op->GetNumberOfLines(0,true)-1,op->GetNumberOfLines(1,true)-1,op->GetNumberOfLines(2,true)-1}; + + //create a pml in x-direction over the full width of yz-space + if (BC[0]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[0]=0; + stop[0] =size[0]; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + if (BC[1]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[0]=op->GetNumberOfLines(0,true)-1-size[1]; + stop[0] =op->GetNumberOfLines(0,true)-1; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + + //create a pml in y-direction over the xz-space (if a pml in x-direction already exists, skip that corner regions) + start[0]=(size[0]+1)*(BC[0]==3); + stop[0] =op->GetNumberOfLines(0,true)-1-(size[0]+1)*(BC[1]==3); + + if (BC[2]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[1]=0; + stop[1] =size[2]; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + if (BC[3]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[1]=op->GetNumberOfLines(1,true)-1-size[3]; + stop[1] =op->GetNumberOfLines(1,true)-1; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + + //create a pml in z-direction over the xy-space (if a pml in x- and/or y-direction already exists, skip that corner/edge regions) + start[1]=(size[2]+1)*(BC[2]==3); + stop[1] =op->GetNumberOfLines(1,true)-1-(size[3]+1)*(BC[3]==3); + + //exclude x-lines that does not belong to the base multi-grid operator; + Operator_CylinderMultiGrid* op_cyl_MG = dynamic_cast<Operator_CylinderMultiGrid*>(op); + if (op_cyl_MG) + start[0] = op_cyl_MG->GetSplitPos()-1; + + if (BC[4]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[2]=0; + stop[2] =size[4]; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + if (BC[5]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[2]=op->GetNumberOfLines(2,true)-1-size[5]; + stop[2] =op->GetNumberOfLines(2,true)-1; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + + BC[1]=0; + size[1]=0; + //create pml extensions (in z-direction only) for child operators in cylindrical multigrid operators + while (op_cyl_MG) + { + Operator_Cylinder* op_child = op_cyl_MG->GetInnerOperator(); + op_cyl_MG = dynamic_cast<Operator_CylinderMultiGrid*>(op_child); + for (int n=0; n<2; ++n) + { + start[n]=0; + stop[n]=op_child->GetNumberOfLines(n,true)-1; + } + + if (op_cyl_MG) + start[0] = op_cyl_MG->GetSplitPos()-1; + + if (BC[4]==3) + { + op_ext_upml = new Operator_Ext_UPML(op_child); + op_ext_upml->SetGradingFunction(gradFunc); + start[2]=0; + stop[2] =size[4]; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op_child->AddExtension(op_ext_upml); + } + if (BC[5]==3) + { + op_ext_upml = new Operator_Ext_UPML(op_child); + op_ext_upml->SetGradingFunction(gradFunc); + start[2]=op->GetNumberOfLines(2,true)-1-size[5]; + stop[2] =op->GetNumberOfLines(2,true)-1; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op_child->AddExtension(op_ext_upml); + } + } + + return true; +} + + +void Operator_Ext_UPML::DeleteOp() +{ + Delete_N_3DArray<FDTD_FLOAT>(vv,m_numLines); + vv = NULL; + Delete_N_3DArray<FDTD_FLOAT>(vvfo,m_numLines); + vvfo = NULL; + Delete_N_3DArray<FDTD_FLOAT>(vvfn,m_numLines); + vvfn = NULL; + Delete_N_3DArray<FDTD_FLOAT>(ii,m_numLines); + ii = NULL; + Delete_N_3DArray<FDTD_FLOAT>(iifo,m_numLines); + iifo = NULL; + Delete_N_3DArray<FDTD_FLOAT>(iifn,m_numLines); + iifn = NULL; +} + + +bool Operator_Ext_UPML::SetGradingFunction(string func) +{ + if (func.empty()) + return true; + + m_GradFunc = func; + int res = m_GradingFunction->Parse(m_GradFunc.c_str(), "D,dl,W,Z,N"); + if (res < 0) return true; + + cerr << "Operator_Ext_UPML::SetGradingFunction: Warning, an error occured parsing the pml grading function (see below) ..." << endl; + cerr << func << "\n" << string(res, ' ') << "^\n" << m_GradingFunction->ErrorMsg() << "\n"; + return false; +} + +void Operator_Ext_UPML::CalcGradingKappa(int ny, unsigned int pos[3], double Zm, double kappa_v[3], double kappa_i[3]) +{ + double depth=0; + double width=0; + for (int n=0; n<3; ++n) + { + if ((pos[n] <= m_Size[2*n]) && (m_BC[2*n]==3)) //lower pml in n-dir + { + width = (m_Op->GetDiscLine(n,m_Size[2*n]) - m_Op->GetDiscLine(n,0))*m_Op->GetGridDelta(); + depth = width - (m_Op->GetDiscLine(n,pos[n]) - m_Op->GetDiscLine(n,0))*m_Op->GetGridDelta(); + + if ((m_Op_Cyl) && (n==1)) + { + width *= m_Op_Cyl->GetDiscLine(0,pos[0]); + depth *= m_Op_Cyl->GetDiscLine(0,pos[0]); + } + + if (n==ny) + depth-=m_Op->GetEdgeLength(n,pos)/2; + double vars[5] = {depth, width/m_Size[2*n], width, Zm, (double)m_Size[2*n]}; + if (depth>0) + kappa_v[n] = m_GradingFunction->Eval(vars); + else + kappa_v[n]=0; + if (n==ny) + depth+=m_Op->GetEdgeLength(n,pos)/2; + + if (n!=ny) + depth-=m_Op->GetEdgeLength(n,pos)/2; + if (depth<0) + depth=0; + vars[0]=depth; + if (depth>0) + kappa_i[n] = m_GradingFunction->Eval(vars); + else + kappa_i[n] = 0; + } + else if ((pos[n] >= m_Op->GetNumberOfLines(n,true) -1 -m_Size[2*n+1]) && (m_BC[2*n+1]==3)) //upper pml in n-dir + { + width = (m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-1) - m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-m_Size[2*n+1]-1))*m_Op->GetGridDelta(); + depth = width - (m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-1) - m_Op->GetDiscLine(n,pos[n]))*m_Op->GetGridDelta(); + + if ((m_Op_Cyl) && (n==1)) + { + width *= m_Op_Cyl->GetDiscLine(0,pos[0]); + depth *= m_Op_Cyl->GetDiscLine(0,pos[0]); + } + + if (n==ny) + depth+=m_Op->GetEdgeLength(n,pos)/2; + double vars[5] = {depth, width/(m_Size[2*n]), width, Zm, (double)m_Size[2*n]}; + if (depth>0) + kappa_v[n] = m_GradingFunction->Eval(vars); + else + kappa_v[n]=0; + if (n==ny) + depth-=m_Op->GetEdgeLength(n,pos)/2; + + if (n!=ny) + depth+=m_Op->GetEdgeLength(n,pos)/2; + if (depth>width) + depth=0; + vars[0]=depth; + if (depth>0) + kappa_i[n] = m_GradingFunction->Eval(vars); + else + kappa_i[n]=0; + } + else + { + kappa_v[n] = 0; + kappa_i[n] = 0; + } + } +} + +bool Operator_Ext_UPML::BuildExtension() +{ + /*Calculate the upml coefficients as defined in: + Allen Taflove, computational electrodynamics - the FDTD method, third edition, chapter 7.8, pages 297-300 + - modified by Thorsten Liebig to match the equivalent circuit (EC) FDTD method + - kappa is used for conductivities (instead of sigma) + */ + if (m_Op==NULL) + return false; + + DeleteOp(); + vv = Create_N_3DArray<FDTD_FLOAT>(m_numLines); + vvfo = Create_N_3DArray<FDTD_FLOAT>(m_numLines); + vvfn = Create_N_3DArray<FDTD_FLOAT>(m_numLines); + ii = Create_N_3DArray<FDTD_FLOAT>(m_numLines); + iifo = Create_N_3DArray<FDTD_FLOAT>(m_numLines); + iifn = Create_N_3DArray<FDTD_FLOAT>(m_numLines); + + unsigned int pos[3]; + unsigned int loc_pos[3]; + int nP,nPP; + double kappa_v[3]={0,0,0}; + double kappa_i[3]={0,0,0}; + double eff_Mat[4]; + double dT = m_Op->GetTimestep(); + + for (loc_pos[0]=0; loc_pos[0]<m_numLines[0]; ++loc_pos[0]) + { + pos[0] = loc_pos[0] + m_StartPos[0]; + for (loc_pos[1]=0; loc_pos[1]<m_numLines[1]; ++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_StartPos[1]; + vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL); + for (loc_pos[2]=0; loc_pos[2]<m_numLines[2]; ++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_StartPos[2]; + for (int n=0; n<3; ++n) + { + m_Op->Calc_EffMatPos(n,pos,eff_Mat,vPrims); + CalcGradingKappa(n, pos,__Z0__ ,kappa_v ,kappa_i); + nP = (n+1)%3; + nPP = (n+2)%3; + if ((kappa_v[0]+kappa_v[1]+kappa_v[2])!=0) + { + //check if pos is on PEC + if ( (m_Op->GetVV(n,pos[0],pos[1],pos[2]) + m_Op->GetVI(n,pos[0],pos[1],pos[2])) != 0 ) + { + //modify the original operator to perform eq. (7.85) by the main engine (EC-FDTD: equation is multiplied by delta_n) + //the engine extension will replace the original voltages with the "voltage flux" (volt*eps0) prior to the voltage updates + //after the updates are done the extension will calculate the new voltages (see below) and place them back into the main field domain + m_Op->SetVV(n,pos[0],pos[1],pos[2], (2*__EPS0__ - kappa_v[nP]*dT) / (2*__EPS0__ + kappa_v[nP]*dT) ); + m_Op->SetVI(n,pos[0],pos[1],pos[2], (2*__EPS0__*dT) / (2*__EPS0__ + kappa_v[nP]*dT) * m_Op->GetEdgeLength(n,pos) / m_Op->GetEdgeArea(n,pos) ); + + + //operators needed by eq. (7.88) to calculate new voltages from old voltages and old and new "voltage fluxes" + GetVV(n,loc_pos) = (2*__EPS0__ - kappa_v[nPP]*dT) / (2*__EPS0__ + kappa_v[nPP]*dT); + GetVVFN(n,loc_pos) = (2*__EPS0__ + kappa_v[n]*dT) / (2*__EPS0__ + kappa_v[nPP]*dT)/eff_Mat[0]; + GetVVFO(n,loc_pos) = (2*__EPS0__ - kappa_v[n]*dT) / (2*__EPS0__ + kappa_v[nPP]*dT)/eff_Mat[0]; + } + } + else + { + //disable upml + GetVV(n,loc_pos) = m_Op->GetVV(n,pos[0],pos[1],pos[2]); + m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 ); + GetVVFO(n,loc_pos) = 0; + GetVVFN(n,loc_pos) = 1; + } + + if ((kappa_i[0]+kappa_i[1]+kappa_i[2])!=0) + { + //check if pos is on PMC + if ( (m_Op->GetII(n,pos[0],pos[1],pos[2]) + m_Op->GetIV(n,pos[0],pos[1],pos[2])) != 0 ) + { + //modify the original operator to perform eq. (7.89) by the main engine (EC-FDTD: equation is multiplied by delta_n) + //the engine extension will replace the original currents with the "current flux" (curr*mu0) prior to the current updates + //after the updates are done the extension will calculate the new currents (see below) and place them back into the main field domain + m_Op->SetII(n,pos[0],pos[1],pos[2], (2*__EPS0__ - kappa_i[nP]*dT) / (2*__EPS0__ + kappa_i[nP]*dT) ); + m_Op->SetIV(n,pos[0],pos[1],pos[2], (2*__EPS0__*dT) / (2*__EPS0__ + kappa_i[nP]*dT) * m_Op->GetEdgeLength(n,pos,true) / m_Op->GetEdgeArea(n,pos,true) ); + + //operators needed by eq. (7.90) to calculate new currents from old currents and old and new "current fluxes" + GetII(n,loc_pos) = (2*__EPS0__ - kappa_i[nPP]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT); + GetIIFN(n,loc_pos) = (2*__EPS0__ + kappa_i[n]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT)/eff_Mat[2]; + GetIIFO(n,loc_pos) = (2*__EPS0__ - kappa_i[n]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT)/eff_Mat[2]; + } + } + else + { + //disable upml + GetII(n,loc_pos) = m_Op->GetII(n,pos[0],pos[1],pos[2]); + m_Op->SetII(n,pos[0],pos[1],pos[2], 0 ); + GetIIFO(n,loc_pos) = 0; + GetIIFN(n,loc_pos) = 1; + } + } + } + } + } + return true; +} + +Engine_Extension* Operator_Ext_UPML::CreateEngineExtention() +{ + Engine_Ext_UPML* eng_ext = new Engine_Ext_UPML(this); + return eng_ext; +} + +void Operator_Ext_UPML::ShowStat(ostream &ostr) const +{ + Operator_Extension::ShowStat(ostr); + + ostr << " PML range\t\t: " << "[" << m_StartPos[0]<< "," << m_StartPos[1]<< "," << m_StartPos[2]<< "] to [" + << m_StartPos[0]+m_numLines[0]-1 << "," << m_StartPos[1]+m_numLines[1]-1 << "," << m_StartPos[2]+m_numLines[2]-1 << "]" << endl; + ostr << " Grading function\t: \"" << m_GradFunc << "\"" << endl; +} diff --git a/openEMS/FDTD/extensions/operator_ext_upml.h b/openEMS/FDTD/extensions/operator_ext_upml.h new file mode 100644 index 0000000..5256655 --- /dev/null +++ b/openEMS/FDTD/extensions/operator_ext_upml.h @@ -0,0 +1,106 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef OPERATOR_EXT_UPML_H +#define OPERATOR_EXT_UPML_H + +#include "FDTD/operator.h" +#include "operator_extension.h" + +class FunctionParser; + +//! Operator extension implemention an uniaxial perfectly matched layer (upml) +/* + The priority for this extension should be the highest of all extensions since this operator will use the main engine to perform vital parts in the upml implementation. + Therefore the voltages and currents as well as the operator are replaced during these update process. + This extension is propably incompatible with the most other extensions operating in the same regions. + */ +class Operator_Ext_UPML : public Operator_Extension +{ + friend class Engine_Ext_UPML; +public: + virtual ~Operator_Ext_UPML(); + + //! Returns always true, Create_UPML method will take care of creating a valid pml for the cylindrical fdtd + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const { UNUSED(closedAlpha); UNUSED(R0_included); return true;} + + //! Returns always true if base grid, Create_UPML will create proper child pml extensions. + virtual bool IsCylindricalMultiGridSave(bool child) const {if (child) return false; return true;} + + virtual bool IsMPISave() const {return true;} + + void SetBoundaryCondition(const int* BCs, const unsigned int size[6]); + + void SetRange(const unsigned int start[3], const unsigned int stop[3]); + + //! Set the grading function for the pml + /*! + Define the pml grading grading function. + Predefined variables in this grading function are: + D = depth in the pml in meter + dl = mesh delta inside the pml in meter + W = width (length) of the pml in meter + N = number of cells for the pml + Z = wave impedance at the current depth and position + example: SetGradingFunction(" -log(1e-6)*log(2.5)/(2*dl*Z*(pow(2.5,W/dl)-1)) * pow(2.5, D/dl) "); + + An empty function string will be ignored. + */ + virtual bool SetGradingFunction(string func); + + virtual bool BuildExtension(); + + virtual Engine_Extension* CreateEngineExtention(); + + virtual string GetExtensionName() const {return string("Uniaxial PML Extension");} + + virtual void ShowStat(ostream &ostr) const; + + //! Create the UPML + static bool Create_UPML(Operator* op, const int ui_BC[6], const unsigned int ui_size[6], const string gradFunc); + +protected: + Operator_Ext_UPML(Operator* op); + int m_BC[6]; + unsigned int m_Size[6]; + + unsigned int m_StartPos[3]; + unsigned int m_numLines[3]; + + string m_GradFunc; + FunctionParser* m_GradingFunction; + + void CalcGradingKappa(int ny, unsigned int pos[3], double Zm, double kappa_v[3], double kappa_i[3]); + + void DeleteOp(); + + virtual FDTD_FLOAT& GetVV(int ny, unsigned int pos[3]) {return vv[ny][pos[0]][pos[1]][pos[2]];} + virtual FDTD_FLOAT& GetVVFO(int ny, unsigned int pos[3]) {return vvfo[ny][pos[0]][pos[1]][pos[2]];} + virtual FDTD_FLOAT& GetVVFN(int ny, unsigned int pos[3]) {return vvfn[ny][pos[0]][pos[1]][pos[2]];} + virtual FDTD_FLOAT& GetII(int ny, unsigned int pos[3]) {return ii[ny][pos[0]][pos[1]][pos[2]];} + virtual FDTD_FLOAT& GetIIFO(int ny, unsigned int pos[3]) {return iifo[ny][pos[0]][pos[1]][pos[2]];} + virtual FDTD_FLOAT& GetIIFN(int ny, unsigned int pos[3]) {return iifn[ny][pos[0]][pos[1]][pos[2]];} + + FDTD_FLOAT**** vv; //calc new voltage from old voltage + FDTD_FLOAT**** vvfo; //calc new voltage from old voltage flux + FDTD_FLOAT**** vvfn; //calc new voltage from new voltage flux + FDTD_FLOAT**** ii; //calc new current from old current + FDTD_FLOAT**** iifo; //calc new current from old current flux + FDTD_FLOAT**** iifn; //calc new current from new current flux +}; + +#endif // OPERATOR_EXT_UPML_H diff --git a/openEMS/FDTD/extensions/operator_extension.cpp b/openEMS/FDTD/extensions/operator_extension.cpp new file mode 100644 index 0000000..50e7068 --- /dev/null +++ b/openEMS/FDTD/extensions/operator_extension.cpp @@ -0,0 +1,54 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#include "operator_extension.h" +#include "FDTD/operator.h" +#include "FDTD/operator_cylinder.h" + +using namespace std; + +Operator_Extension::Operator_Extension(Operator* op) +{ + m_Op = op; + m_Active = true; + + m_CC_R0_included = false; + m_Op_Cyl = dynamic_cast<Operator_Cylinder*>(op); + if (m_Op_Cyl) + m_CC_R0_included=m_Op_Cyl->GetR0Included(); + m_Eng_Ext = NULL; +} + +Operator_Extension::~Operator_Extension() +{ +} + +Operator_Extension::Operator_Extension(Operator* op, Operator_Extension* op_ext) +{ + UNUSED(op_ext); + m_Op = op; + m_Op_Cyl = dynamic_cast<Operator_Cylinder*>(op); + m_Active = op_ext->m_Active; + if (m_Op_Cyl) + m_CC_R0_included=m_Op_Cyl->GetR0Included(); + m_Eng_Ext = NULL; +} + +void Operator_Extension::ShowStat(ostream &ostr) const +{ + ostr << "--- " << GetExtensionName() << " ---" << endl; +} diff --git a/openEMS/FDTD/extensions/operator_extension.h b/openEMS/FDTD/extensions/operator_extension.h new file mode 100644 index 0000000..21e42af --- /dev/null +++ b/openEMS/FDTD/extensions/operator_extension.h @@ -0,0 +1,87 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +*/ + +#ifndef OPERATOR_EXTENSION_H +#define OPERATOR_EXTENSION_H + +#include <string> + +#include <iostream> +#include <stdio.h> +#include <stdlib.h> + +#include "tools/global.h" + +class Operator; +class Operator_Cylinder; +class Engine_Extension; + +//! Abstract base-class for all operator extensions +class Operator_Extension +{ + friend class Engine_Extension; +public: + virtual ~Operator_Extension(); + + //! Create a clone of this extension, will return NULL if this is impossible + /*! + Create a clone of this extension, will return NULL if this is impossible (e.g. derived extension has no clone method and copy-constructor)... + BuildExtension has to be called separatly! + */ + virtual Operator_Extension* Clone(Operator* op) {UNUSED(op); return NULL;} + + virtual bool BuildExtension() {return true;} + + virtual Engine_Extension* CreateEngineExtention() {return 0;} + virtual Engine_Extension* GetEngineExtention() {return m_Eng_Ext;} + + //! The cylindrical operator will check whether the extension is save to use. Default is false. Derive this method to override. + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return false;} + + //! The cylindrical multi grid operator will check whether the extension is save to use. Default is false. Derive this method to override. + virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return false;} + + //! The MPI operator (if enabled) will check whether the extension is compatible with MPI. Default is false. Derive this method to override. + virtual bool IsMPISave() const {return false;} + + virtual std::string GetExtensionName() const {return std::string("Abstract Operator Extension Base Class");} + + virtual void ShowStat(std::ostream &ostr) const; + + virtual bool IsActive() const {return m_Active;} + virtual void SetActive(bool active=true) {m_Active=active;} + + virtual void Init() {} + virtual void Reset() {} + +protected: + Operator_Extension(Operator* op); + //! Copy constructor + Operator_Extension(Operator* op, Operator_Extension* op_ext); + + bool m_Active; + + //FDTD Operator + Operator* m_Op; + Engine_Extension* m_Eng_Ext; + + //Cylindrical FDTD Operator (not NULL if a cylindrical FDTD is used) + Operator_Cylinder* m_Op_Cyl; + bool m_CC_R0_included; +}; + +#endif // OPERATOR_EXTENSION_H |