2010-03-11 15:47:40 +00:00
/*
* 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/>.
*/
2010-03-01 18:35:28 +00:00
# include <fstream>
2010-03-01 08:19:39 +00:00
# include "operator.h"
2010-05-08 16:12:44 +00:00
# include "engine.h"
2010-04-25 19:59:05 +00:00
# include "operator_extension.h"
2010-12-06 09:44:25 +00:00
# include "Common/processfields.h"
2010-03-01 13:56:27 +00:00
# include "tools/array_ops.h"
2010-04-07 14:31:23 +00:00
# include "fparser.hh"
2010-03-01 08:19:39 +00:00
2010-04-09 13:51:37 +00:00
Operator * Operator : : New ( )
{
2010-05-19 19:25:15 +00:00
cout < < " Create FDTD operator " < < endl ;
2010-04-09 13:51:37 +00:00
Operator * op = new Operator ( ) ;
op - > Init ( ) ;
return op ;
}
2010-12-06 08:59:42 +00:00
Operator : : Operator ( ) : Operator_Base ( )
2010-03-01 08:19:39 +00:00
{
2010-05-03 16:33:14 +00:00
Exc = 0 ;
2010-09-08 14:07:28 +00:00
m_InvaildTimestep = false ;
2010-03-01 08:19:39 +00:00
}
Operator : : ~ Operator ( )
{
2010-04-27 21:06:42 +00:00
for ( size_t n = 0 ; n < m_Op_exts . size ( ) ; + + n )
delete m_Op_exts . at ( n ) ;
m_Op_exts . clear ( ) ;
2010-04-09 13:51:37 +00:00
Reset ( ) ;
2010-03-01 08:19:39 +00:00
}
2010-05-09 18:27:17 +00:00
Engine * Operator : : CreateEngine ( ) const
2010-05-08 16:12:44 +00:00
{
Engine * eng = Engine : : New ( this ) ;
return eng ;
}
2010-03-01 08:19:39 +00:00
void Operator : : Init ( )
{
2010-12-06 08:59:42 +00:00
Operator_Base : : Init ( ) ;
2010-03-01 08:19:39 +00:00
CSX = NULL ;
2010-03-01 13:56:27 +00:00
vv = NULL ;
vi = NULL ;
iv = NULL ;
ii = NULL ;
2010-03-05 13:20:25 +00:00
MainOp = NULL ;
DualOp = NULL ;
for ( int n = 0 ; n < 3 ; + + n )
{
EC_C [ n ] = NULL ;
EC_G [ n ] = NULL ;
EC_L [ n ] = NULL ;
EC_R [ n ] = NULL ;
}
2010-04-28 20:01:03 +00:00
2010-05-03 16:33:14 +00:00
Exc = 0 ;
2010-03-01 08:19:39 +00:00
}
void Operator : : Reset ( )
{
2010-03-01 13:56:27 +00:00
Delete_N_3DArray ( vv , numLines ) ;
Delete_N_3DArray ( vi , numLines ) ;
Delete_N_3DArray ( iv , numLines ) ;
Delete_N_3DArray ( ii , numLines ) ;
2010-03-05 13:20:25 +00:00
delete MainOp ;
delete DualOp ;
for ( int n = 0 ; n < 3 ; + + n )
{
delete [ ] EC_C [ n ] ;
delete [ ] EC_G [ n ] ;
delete [ ] EC_L [ n ] ;
delete [ ] EC_R [ n ] ;
}
2010-05-03 16:33:14 +00:00
delete Exc ;
2010-03-01 08:19:39 +00:00
2010-12-06 08:59:42 +00:00
Operator_Base : : Reset ( ) ;
2010-04-23 06:17:42 +00:00
}
2010-04-19 14:08:32 +00:00
double Operator : : GetMeshDelta ( int n , const unsigned int * pos , bool dualMesh ) const
2010-04-13 16:40:43 +00:00
{
if ( ( n < 0 ) | | ( n > 2 ) ) return 0.0 ;
2010-09-17 12:51:07 +00:00
if ( pos [ n ] > = numLines [ n ] ) return 0.0 ;
2010-04-13 16:40:43 +00:00
if ( dualMesh = = false )
return fabs ( MainOp - > GetIndexDelta ( n , pos [ n ] ) ) * gridDelta ;
else
return fabs ( MainOp - > GetIndexWidth ( n , pos [ n ] ) ) * gridDelta ;
}
double Operator : : GetDiscLine ( int n , unsigned int pos , bool dualMesh ) const
{
if ( ( n < 0 ) | | ( n > 2 ) ) return 0.0 ;
2010-09-17 12:51:07 +00:00
if ( pos > = numLines [ n ] ) return 0.0 ;
2010-04-13 16:40:43 +00:00
if ( dualMesh = = false )
return discLines [ n ] [ pos ] ;
else
2010-07-29 16:30:50 +00:00
return ( discLines [ n ] [ pos ] + 0.5 * fabs ( MainOp - > GetIndexDelta ( n , pos ) ) ) ;
}
2010-09-17 12:51:07 +00:00
double Operator : : GetNodeArea ( int ny , const unsigned int pos [ 3 ] , bool dualMesh ) const
2010-07-29 16:30:50 +00:00
{
int nyP = ( ny + 1 ) % 3 ;
int nyPP = ( ny + 2 ) % 3 ;
2010-09-02 13:35:13 +00:00
return GetNodeWidth ( nyP , pos , dualMesh ) * GetNodeWidth ( nyPP , pos , dualMesh ) ;
2010-04-13 16:40:43 +00:00
}
2010-04-09 13:58:15 +00:00
bool Operator : : SnapToMesh ( double * dcoord , unsigned int * uicoord , bool lower , bool * inside )
2010-03-01 18:35:28 +00:00
{
bool ok = true ;
2010-04-11 21:42:54 +00:00
unsigned int numLines [ 3 ] ;
2010-03-01 18:35:28 +00:00
for ( int n = 0 ; n < 3 ; + + n )
{
2010-04-11 21:42:54 +00:00
numLines [ n ] = GetNumberOfLines ( n ) ;
2010-04-09 13:58:15 +00:00
if ( inside ) //set defaults
inside [ n ] = true ;
2010-03-02 13:54:50 +00:00
uicoord [ n ] = 0 ;
2010-04-11 21:42:54 +00:00
if ( dcoord [ n ] < discLines [ n ] [ 0 ] )
{
ok = false ; uicoord [ n ] = 0 ;
if ( inside ) inside [ n ] = false ;
}
else if ( dcoord [ n ] = = discLines [ n ] [ 0 ] )
uicoord [ n ] = 0 ;
else if ( dcoord [ n ] > discLines [ n ] [ numLines [ n ] - 1 ] )
{
ok = false ;
uicoord [ n ] = numLines [ n ] - 1 ;
if ( lower ) uicoord [ n ] = numLines [ n ] - 2 ;
if ( inside ) inside [ n ] = false ;
}
else if ( dcoord [ n ] = = discLines [ n ] [ numLines [ n ] - 1 ] )
{
uicoord [ n ] = numLines [ n ] - 1 ;
if ( lower ) uicoord [ n ] = numLines [ n ] - 2 ;
}
2010-03-02 13:54:50 +00:00
else
2010-03-27 22:05:45 +00:00
for ( unsigned int i = 1 ; i < numLines [ n ] ; + + i )
2010-03-01 18:35:28 +00:00
{
2010-03-02 14:37:00 +00:00
if ( dcoord [ n ] < discLines [ n ] [ i ] )
2010-03-02 13:54:50 +00:00
{
2010-03-02 14:37:00 +00:00
if ( fabs ( dcoord [ n ] - discLines [ n ] [ i ] ) < ( fabs ( dcoord [ n ] - discLines [ n ] [ i - 1 ] ) ) )
2010-03-02 13:54:50 +00:00
uicoord [ n ] = i ;
else
2010-03-02 14:37:00 +00:00
uicoord [ n ] = i - 1 ;
if ( lower ) uicoord [ n ] = i - 1 ;
2010-03-02 13:54:50 +00:00
i = numLines [ n ] ;
}
2010-03-01 18:35:28 +00:00
}
}
2010-03-02 13:54:50 +00:00
// cerr << "Operator::SnapToMesh Wish: " << dcoord[0] << " " << dcoord[1] << " " << dcoord[2] << endl;
// cerr << "Operator::SnapToMesh Found: " << discLines[0][uicoord[0]] << " " << discLines[1][uicoord[1]] << " " << discLines[2][uicoord[2]] << endl;
// cerr << "Operator::SnapToMesh Index: " << uicoord[0] << " " << uicoord[1] << " " << uicoord[2] << endl;
2010-03-01 18:35:28 +00:00
return ok ;
}
2010-03-22 07:19:17 +00:00
struct Operator : : Grid_Path Operator : : FindPath ( double start [ ] , double stop [ ] )
{
struct Grid_Path path ;
2010-04-13 16:51:44 +00:00
// double dV[] = {stop[0]-start[0],stop[1]-start[1],stop[2]-start[2]};
2010-03-22 07:19:17 +00:00
2010-04-13 16:51:44 +00:00
unsigned int uiStart [ 3 ] , uiStop [ 3 ] , currPos [ 3 ] ;
2010-03-22 07:19:17 +00:00
SnapToMesh ( start , uiStart ) ;
SnapToMesh ( stop , uiStop ) ;
currPos [ 0 ] = uiStart [ 0 ] ;
currPos [ 1 ] = uiStart [ 1 ] ;
currPos [ 2 ] = uiStart [ 2 ] ;
double meshStart [ ] = { discLines [ 0 ] [ uiStart [ 0 ] ] , discLines [ 1 ] [ uiStart [ 1 ] ] , discLines [ 2 ] [ uiStart [ 2 ] ] } ;
double meshStop [ ] = { discLines [ 0 ] [ uiStop [ 0 ] ] , discLines [ 1 ] [ uiStop [ 1 ] ] , discLines [ 2 ] [ uiStop [ 2 ] ] } ;
2010-04-13 16:51:44 +00:00
bool UpDir = false ;
2010-03-26 18:51:09 +00:00
double foot = 0 , dist = 0 , minFoot = 0 , minDist = 0 ;
int minDir = 0 ;
2010-03-22 07:19:17 +00:00
unsigned int minPos [ 3 ] ;
double startFoot , stopFoot , currFoot ;
Point_Line_Distance ( meshStart , start , stop , startFoot , dist ) ;
Point_Line_Distance ( meshStop , start , stop , stopFoot , dist ) ;
currFoot = startFoot ;
2010-03-26 18:33:44 +00:00
minFoot = startFoot ;
2010-03-22 07:19:17 +00:00
double P [ 3 ] ;
while ( minFoot < stopFoot )
{
minDist = 1e300 ;
for ( int n = 0 ; n < 3 ; + + n ) //check all 6 surrounding points
{
P [ 0 ] = discLines [ 0 ] [ currPos [ 0 ] ] ;
P [ 1 ] = discLines [ 1 ] [ currPos [ 1 ] ] ;
P [ 2 ] = discLines [ 2 ] [ currPos [ 2 ] ] ;
2010-03-26 18:51:09 +00:00
if ( ( ( int ) currPos [ n ] - 1 ) > = 0 )
2010-03-22 07:19:17 +00:00
{
P [ n ] = discLines [ n ] [ currPos [ n ] - 1 ] ;
Point_Line_Distance ( P , start , stop , foot , dist ) ;
if ( ( foot > currFoot ) & & ( dist < minDist ) )
{
minFoot = foot ;
minDist = dist ;
minDir = n ;
2010-03-26 18:33:44 +00:00
UpDir = false ;
2010-03-22 07:19:17 +00:00
}
}
if ( ( currPos [ n ] + 1 ) < numLines [ n ] )
{
P [ n ] = discLines [ n ] [ currPos [ n ] + 1 ] ;
Point_Line_Distance ( P , start , stop , foot , dist ) ;
if ( ( foot > currFoot ) & & ( dist < minDist ) )
{
minFoot = foot ;
minDist = dist ;
minDir = n ;
2010-03-26 18:33:44 +00:00
UpDir = true ;
2010-03-22 07:19:17 +00:00
}
}
}
2010-03-26 18:33:44 +00:00
minPos [ 0 ] = currPos [ 0 ] ;
minPos [ 1 ] = currPos [ 1 ] ;
minPos [ 2 ] = currPos [ 2 ] ;
if ( UpDir )
{
currPos [ minDir ] + = 1 ;
}
else
{
currPos [ minDir ] + = - 1 ;
minPos [ minDir ] - = 1 ;
}
2010-03-22 07:19:17 +00:00
path . posPath [ 0 ] . push_back ( minPos [ 0 ] ) ;
path . posPath [ 1 ] . push_back ( minPos [ 1 ] ) ;
path . posPath [ 2 ] . push_back ( minPos [ 2 ] ) ;
2010-03-26 18:33:44 +00:00
currFoot = minFoot ;
path . dir . push_back ( minDir ) ;
2010-03-22 07:19:17 +00:00
}
2010-10-12 10:49:44 +00:00
//close missing edges, if currPos is not equal to uiStopPos
for ( int n = 0 ; n < 3 ; + + n )
{
if ( currPos [ n ] > uiStop [ n ] )
{
- - currPos [ n ] ;
path . posPath [ 0 ] . push_back ( currPos [ 0 ] ) ;
path . posPath [ 1 ] . push_back ( currPos [ 1 ] ) ;
path . posPath [ 2 ] . push_back ( currPos [ 2 ] ) ;
path . dir . push_back ( n ) ;
}
else if ( currPos [ n ] < uiStop [ n ] )
{
path . posPath [ 0 ] . push_back ( currPos [ 0 ] ) ;
path . posPath [ 1 ] . push_back ( currPos [ 1 ] ) ;
path . posPath [ 2 ] . push_back ( currPos [ 2 ] ) ;
path . dir . push_back ( n ) ;
}
}
2010-03-22 07:19:17 +00:00
return path ;
}
2010-03-31 14:35:43 +00:00
double Operator : : GetNumberCells ( ) const
2010-03-01 13:56:27 +00:00
{
if ( numLines )
2010-03-29 20:11:24 +00:00
return ( numLines [ 0 ] ) * ( numLines [ 1 ] ) * ( numLines [ 2 ] ) ; //it's more like number of nodes???
2010-03-01 13:56:27 +00:00
return 0 ;
}
2010-03-31 14:35:43 +00:00
void Operator : : ShowStat ( ) const
2010-03-01 13:56:27 +00:00
{
unsigned int OpSize = 12 * numLines [ 0 ] * numLines [ 1 ] * numLines [ 2 ] * sizeof ( FDTD_FLOAT ) ;
unsigned int FieldSize = 6 * numLines [ 0 ] * numLines [ 1 ] * numLines [ 2 ] * sizeof ( FDTD_FLOAT ) ;
double MBdiff = 1024 * 1024 ;
2010-03-29 20:11:24 +00:00
cout < < " ------- Stat: FDTD Operator ------- " < < endl ;
2010-05-10 07:14:29 +00:00
cout < < " Dimensions \t \t : " < < numLines [ 0 ] < < " x " < < numLines [ 1 ] < < " x " < < numLines [ 2 ] < < " = " < < numLines [ 0 ] * numLines [ 1 ] * numLines [ 2 ] < < " Cells ( " < < numLines [ 0 ] * numLines [ 1 ] * numLines [ 2 ] / 1e6 < < " MCells) " < < endl ;
2010-05-21 08:14:09 +00:00
cout < < " Size of Operator \t : " < < OpSize < < " Byte ( " < < ( double ) OpSize / MBdiff < < " MiB) " < < endl ;
cout < < " Size of Field-Data \t : " < < FieldSize < < " Byte ( " < < ( double ) FieldSize / MBdiff < < " MiB) " < < endl ;
2010-03-29 20:11:24 +00:00
cout < < " ----------------------------------- " < < endl ;
2010-08-16 11:28:19 +00:00
cout < < " Voltage excitations \t : " < < Exc - > Volt_Count < < " \t ( " < < Exc - > Volt_Count_Dir [ 0 ] < < " , " < < Exc - > Volt_Count_Dir [ 1 ] < < " , " < < Exc - > Volt_Count_Dir [ 2 ] < < " ) " < < endl ;
2010-08-16 11:22:40 +00:00
cout < < " Current excitations \t : " < < Exc - > Curr_Count < < " \t ( " < < Exc - > Curr_Count_Dir [ 0 ] < < " , " < < Exc - > Curr_Count_Dir [ 1 ] < < " , " < < Exc - > Curr_Count_Dir [ 2 ] < < " ) " < < endl ;
2010-05-19 09:42:56 +00:00
cout < < " ----------------------------------- " < < endl ;
2010-05-10 07:14:29 +00:00
cout < < " Number of PEC edges \t : " < < m_Nr_PEC [ 0 ] + m_Nr_PEC [ 1 ] + m_Nr_PEC [ 2 ] < < endl ;
cout < < " in " < < GetDirName ( 0 ) < < " direction \t \t : " < < m_Nr_PEC [ 0 ] < < endl ;
cout < < " in " < < GetDirName ( 1 ) < < " direction \t \t : " < < m_Nr_PEC [ 1 ] < < endl ;
cout < < " in " < < GetDirName ( 2 ) < < " direction \t \t : " < < m_Nr_PEC [ 2 ] < < endl ;
cout < < " ----------------------------------- " < < endl ;
2010-09-02 20:12:03 +00:00
cout < < " Timestep (s) \t \t : " < < dT ;
if ( opt_dT )
cout < < " \t ( " < < opt_dT < < " ) " ;
cout < < endl ;
2010-07-13 13:37:56 +00:00
cout < < " Timestep method name \t : " < < m_Used_TS_Name < < endl ;
2010-05-10 07:14:29 +00:00
cout < < " Nyquist criteria (TS) \t : " < < Exc - > GetNyquistNum ( ) < < endl ;
cout < < " Nyquist criteria (s) \t : " < < Exc - > GetNyquistNum ( ) * dT < < endl ;
cout < < " Excitation Length (TS) \t : " < < Exc - > Length < < endl ;
cout < < " Excitation Length (s) \t : " < < Exc - > Length * dT < < endl ;
2010-03-29 20:11:24 +00:00
cout < < " ----------------------------------- " < < endl ;
2010-03-01 13:56:27 +00:00
}
2010-07-11 21:45:41 +00:00
void Operator : : ShowExtStat ( ) const
{
2010-07-13 13:37:56 +00:00
if ( m_Op_exts . size ( ) = = 0 ) return ;
2010-07-11 21:45:41 +00:00
cout < < " ----------------------------------- " < < endl ;
for ( size_t n = 0 ; n < m_Op_exts . size ( ) ; + + n )
m_Op_exts . at ( n ) - > ShowStat ( cout ) ;
cout < < " ----------------------------------- " < < endl ;
}
2010-03-09 20:34:23 +00:00
2010-03-01 18:35:28 +00:00
void Operator : : DumpOperator2File ( string filename )
{
2010-07-15 10:58:48 +00:00
# ifdef OUTPUT_IN_DRAWINGUNITS
double discLines_scaling = 1 ;
# else
double discLines_scaling = GetGridDelta ( ) ;
# endif
2010-03-01 18:35:28 +00:00
ofstream file ( filename . c_str ( ) , ios_base : : out ) ;
2010-10-27 09:17:58 +00:00
if ( ! file . is_open ( ) )
2010-03-01 18:35:28 +00:00
{
2010-10-27 09:17:58 +00:00
cerr < < " Operator::DumpOperator2File(): Can't open file: " < < filename < < endl ;
2010-03-01 18:35:28 +00:00
return ;
}
2010-10-27 09:17:58 +00:00
cout < < " Operator: Dumping FDTD operator information to vtk file: " < < filename < < " ... " < < flush ;
2010-06-07 21:08:38 +00:00
2010-06-06 18:00:24 +00:00
FDTD_FLOAT * * * * exc = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
2010-10-27 09:17:58 +00:00
if ( Exc )
{
2010-08-16 11:28:19 +00:00
for ( unsigned int n = 0 ; n < Exc - > Volt_Count ; + + n )
exc [ Exc - > Volt_dir [ n ] ] [ Exc - > Volt_index [ 0 ] [ n ] ] [ Exc - > Volt_index [ 1 ] [ n ] ] [ Exc - > Volt_index [ 2 ] [ n ] ] = Exc - > Volt_amp [ n ] ;
2010-03-27 22:05:45 +00:00
}
2010-10-27 09:17:58 +00:00
FDTD_FLOAT * * * * vv_temp = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
FDTD_FLOAT * * * * vi_temp = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
FDTD_FLOAT * * * * iv_temp = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
FDTD_FLOAT * * * * ii_temp = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
unsigned int pos [ 3 ] , n ;
for ( n = 0 ; n < 3 ; n + + )
for ( pos [ 0 ] = 0 ; pos [ 0 ] < numLines [ 0 ] ; pos [ 0 ] + + )
for ( pos [ 1 ] = 0 ; pos [ 1 ] < numLines [ 1 ] ; pos [ 1 ] + + )
for ( pos [ 2 ] = 0 ; pos [ 2 ] < numLines [ 2 ] ; pos [ 2 ] + + )
{
vv_temp [ n ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetVV ( n , pos ) ;
vi_temp [ n ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetVI ( n , pos ) ;
iv_temp [ n ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetIV ( n , pos ) ;
ii_temp [ n ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetII ( n , pos ) ;
}
2010-03-27 22:05:45 +00:00
string names [ ] = { " vv " , " vi " , " iv " , " ii " , " exc " } ;
2010-10-27 09:17:58 +00:00
FDTD_FLOAT * * * * array [ ] = { vv_temp , vi_temp , iv_temp , ii_temp , exc } ;
2010-03-27 22:05:45 +00:00
2010-07-15 10:58:48 +00:00
ProcessFields : : DumpMultiVectorArray2VTK ( file , names , array , 5 , discLines , numLines , 6 , " Operator dump " , ( ProcessFields : : MeshType ) m_MeshType , discLines_scaling ) ;
2010-03-01 18:35:28 +00:00
2010-10-27 09:17:58 +00:00
Delete_N_3DArray ( ii_temp , numLines ) ;
Delete_N_3DArray ( iv_temp , numLines ) ;
Delete_N_3DArray ( vi_temp , numLines ) ;
Delete_N_3DArray ( vv_temp , numLines ) ;
2010-03-27 22:05:45 +00:00
Delete_N_3DArray ( exc , numLines ) ;
2010-03-01 18:35:28 +00:00
file . close ( ) ;
2010-06-04 12:08:42 +00:00
cout < < " done! " < < endl ;
2010-03-01 18:35:28 +00:00
}
2010-03-05 13:20:25 +00:00
2010-06-02 14:37:21 +00:00
//! \brief dump PEC (perfect electric conductor) information (into VTK-file)
//! visualization via paraview
//! visualize only one component (x, y or z)
void Operator : : DumpPEC2File ( string filename )
{
ofstream file ( filename . c_str ( ) ) ;
2010-10-27 09:17:58 +00:00
if ( ! file . is_open ( ) )
{
cerr < < " Operator::DumpPEC2File(): Can't open file: " < < filename < < endl ;
2010-06-02 14:37:21 +00:00
return ;
}
2010-10-27 09:17:58 +00:00
cout < < " Operator: Dumping PEC information to vtk file: " < < filename < < " ... " < < flush ;
2010-06-04 12:08:42 +00:00
2010-06-06 18:00:24 +00:00
FDTD_FLOAT * * * * pec = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
2010-06-02 14:37:21 +00:00
unsigned int pos [ 3 ] ;
2010-09-17 08:50:06 +00:00
# ifdef OUTPUT_IN_DRAWINGUNITS
2010-09-17 10:08:49 +00:00
double scaling = 1.0 / GetGridDelta ( ) ;
2010-09-17 08:50:06 +00:00
# else
2010-09-17 10:08:49 +00:00
double scaling = 1 ;
2010-09-17 08:50:06 +00:00
# endif
2010-10-20 07:25:50 +00:00
for ( pos [ 0 ] = 0 ; pos [ 0 ] < numLines [ 0 ] - 1 ; pos [ 0 ] + + ) {
for ( pos [ 1 ] = 0 ; pos [ 1 ] < numLines [ 1 ] - 1 ; pos [ 1 ] + + ) {
for ( pos [ 2 ] = 0 ; pos [ 2 ] < numLines [ 2 ] - 1 ; pos [ 2 ] + + ) {
if ( ( pos [ 1 ] ! = 0 ) & & ( pos [ 2 ] ! = 0 ) )
{
// PEC surrounds the computational area; do not output this
2010-10-27 09:17:58 +00:00
if ( ( GetVV ( 0 , pos ) = = 0 ) & & ( GetVI ( 0 , pos ) = = 0 ) )
2010-10-20 07:25:50 +00:00
pec [ 0 ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetEdgeLength ( 0 , pos ) * scaling ; // PEC-x found
}
if ( ( pos [ 0 ] ! = 0 ) & & ( pos [ 2 ] ! = 0 ) )
{
// PEC surrounds the computational area; do not output this
2010-10-27 09:17:58 +00:00
if ( ( GetVV ( 1 , pos ) = = 0 ) & & ( GetVI ( 1 , pos ) = = 0 ) )
2010-10-20 07:25:50 +00:00
pec [ 1 ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetEdgeLength ( 1 , pos ) * scaling ; // PEC-y found
}
if ( ( pos [ 0 ] ! = 0 ) & & ( pos [ 1 ] ! = 0 ) )
{
// PEC surrounds the computational area; do not output this
2010-10-27 09:17:58 +00:00
if ( ( GetVV ( 2 , pos ) = = 0 ) & & ( GetVI ( 2 , pos ) = = 0 ) )
2010-10-20 07:25:50 +00:00
pec [ 2 ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetEdgeLength ( 2 , pos ) * scaling ; // PEC-z found
}
}
}
}
// evaluate boundary conditions
for ( int n = 0 ; n < 3 ; n + + )
{
int nP = ( n + 1 ) % 3 ;
int nPP = ( n + 2 ) % 3 ;
for ( pos [ nP ] = 0 ; pos [ nP ] < numLines [ nP ] ; pos [ nP ] + + )
{
for ( pos [ nPP ] = 0 ; pos [ nPP ] < numLines [ nPP ] ; pos [ nPP ] + + )
{
pos [ n ] = 0 ;
if ( ( pos [ nP ] ! = numLines [ nP ] - 1 ) & & ( m_BC [ 2 * n ] = = 0 ) )
pec [ nP ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetEdgeLength ( nP , pos ) * scaling ;
if ( ( pos [ nPP ] ! = numLines [ nPP ] - 1 ) & & ( m_BC [ 2 * n ] = = 0 ) )
pec [ nPP ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetEdgeLength ( nPP , pos ) * scaling ;
pos [ n ] = numLines [ n ] - 1 ;
if ( ( pos [ nP ] ! = numLines [ nP ] - 1 ) & & ( m_BC [ 2 * n + 1 ] = = 0 ) )
pec [ nP ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetEdgeLength ( nP , pos ) * scaling ;
if ( ( pos [ nPP ] ! = numLines [ nPP ] - 1 ) & & ( m_BC [ 2 * n + 1 ] = = 0 ) )
pec [ nPP ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = GetEdgeLength ( nPP , pos ) * scaling ;
2010-06-02 14:37:21 +00:00
}
}
}
2010-09-17 10:08:49 +00:00
# ifdef OUTPUT_IN_DRAWINGUNITS
scaling = 1 ;
# else
scaling = GetGridDelta ( ) ;
# endif
2010-09-17 08:50:06 +00:00
ProcessFields : : DumpVectorArray2VTK ( file , " PEC " , pec , discLines , numLines , 6 , " PEC dump " , ( ProcessFields : : MeshType ) m_MeshType , scaling ) ;
2010-06-02 14:37:21 +00:00
file . close ( ) ;
2010-06-04 12:08:42 +00:00
cout < < " done! " < < endl ;
2010-06-02 14:37:21 +00:00
}
2010-03-12 19:39:04 +00:00
void Operator : : DumpMaterial2File ( string filename )
{
2010-07-15 10:58:48 +00:00
# ifdef OUTPUT_IN_DRAWINGUNITS
double discLines_scaling = 1 ;
# else
double discLines_scaling = GetGridDelta ( ) ;
# endif
2010-06-04 12:08:42 +00:00
ofstream file ( filename . c_str ( ) , ios_base : : out ) ;
2010-10-27 09:17:58 +00:00
if ( ! file . is_open ( ) )
2010-06-04 12:08:42 +00:00
{
2010-10-27 09:17:58 +00:00
cerr < < " Operator::DumpMaterial2File(): Can't open file: " < < filename < < endl ;
2010-06-04 12:08:42 +00:00
return ;
}
2010-10-27 09:17:58 +00:00
cout < < " Operator: Dumping material information to vtk file: " < < filename < < " ... " < < flush ;
2010-06-04 12:08:42 +00:00
2010-10-27 09:17:58 +00:00
FDTD_FLOAT * * * * epsilon = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
FDTD_FLOAT * * * * mue = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
FDTD_FLOAT * * * * kappa = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
FDTD_FLOAT * * * * sigma = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
2010-03-12 19:39:04 +00:00
2010-10-27 09:17:58 +00:00
unsigned int pos [ 3 ] ;
2010-03-12 19:39:04 +00:00
for ( pos [ 0 ] = 0 ; pos [ 0 ] < numLines [ 0 ] ; + + pos [ 0 ] )
{
for ( pos [ 1 ] = 0 ; pos [ 1 ] < numLines [ 1 ] ; + + pos [ 1 ] )
{
for ( pos [ 2 ] = 0 ; pos [ 2 ] < numLines [ 2 ] ; + + pos [ 2 ] )
{
for ( int n = 0 ; n < 3 ; + + n )
{
2010-10-27 09:17:58 +00:00
double inMat [ 4 ] ;
2010-03-12 19:39:04 +00:00
Calc_EffMatPos ( n , pos , inMat ) ;
2010-10-27 09:17:58 +00:00
epsilon [ n ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = inMat [ 0 ] / __EPS0__ ;
mue [ n ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = inMat [ 2 ] / __MUE0__ ;
kappa [ n ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = inMat [ 1 ] ;
sigma [ n ] [ pos [ 0 ] ] [ pos [ 1 ] ] [ pos [ 2 ] ] = inMat [ 3 ] ;
2010-03-12 19:39:04 +00:00
}
}
}
}
string names [ ] = { " epsilon " , " mue " , " kappa " , " sigma " } ;
2010-09-02 13:35:57 +00:00
FDTD_FLOAT * * * * array [ ] = { epsilon , mue , kappa , sigma } ;
2010-10-27 09:17:58 +00:00
2010-09-02 13:35:57 +00:00
ProcessFields : : DumpMultiVectorArray2VTK ( file , names , array , 4 , discLines , numLines , 6 , " Material dump " , ( ProcessFields : : MeshType ) m_MeshType , discLines_scaling ) ;
2010-10-27 09:17:58 +00:00
2010-09-02 13:35:57 +00:00
Delete_N_3DArray ( epsilon , numLines ) ;
Delete_N_3DArray ( mue , numLines ) ;
Delete_N_3DArray ( kappa , numLines ) ;
Delete_N_3DArray ( sigma , numLines ) ;
2010-10-27 09:17:58 +00:00
2010-03-12 19:39:04 +00:00
file . close ( ) ;
2010-06-04 12:08:42 +00:00
cout < < " done! " < < endl ;
2010-03-12 19:39:04 +00:00
}
2010-03-11 09:48:47 +00:00
bool Operator : : SetGeometryCSX ( ContinuousStructure * geo )
2010-03-05 13:20:25 +00:00
{
2010-03-11 09:48:47 +00:00
if ( geo = = NULL ) return false ;
2010-03-05 13:20:25 +00:00
CSX = geo ;
CSRectGrid * grid = CSX - > GetGrid ( ) ;
for ( int n = 0 ; n < 3 ; + + n )
{
discLines [ n ] = grid - > GetLines ( n , discLines [ n ] , numLines [ n ] , true ) ;
2010-04-13 16:40:43 +00:00
if ( n = = 1 )
2010-03-11 09:48:47 +00:00
if ( numLines [ n ] < 3 ) { cerr < < " CartOperator::SetGeometryCSX: you need at least 3 disc-lines in every direction (3D!)!!! " < < endl ; Reset ( ) ; return false ; }
2010-03-05 13:20:25 +00:00
}
MainOp = new AdrOp ( numLines [ 0 ] , numLines [ 1 ] , numLines [ 2 ] ) ;
MainOp - > SetGrid ( discLines [ 0 ] , discLines [ 1 ] , discLines [ 2 ] ) ;
2010-03-11 09:48:47 +00:00
if ( grid - > GetDeltaUnit ( ) < = 0 ) { cerr < < " CartOperator::SetGeometryCSX: grid delta unit must not be <=0 !!! " < < endl ; Reset ( ) ; return false ; }
2010-03-05 13:20:25 +00:00
else gridDelta = grid - > GetDeltaUnit ( ) ;
MainOp - > SetGridDelta ( 1 ) ;
MainOp - > AddCellAdrOp ( ) ;
2010-03-11 09:48:47 +00:00
return true ;
2010-03-05 13:20:25 +00:00
}
void Operator : : InitOperator ( )
{
Delete_N_3DArray ( vv , numLines ) ;
Delete_N_3DArray ( vi , numLines ) ;
Delete_N_3DArray ( iv , numLines ) ;
Delete_N_3DArray ( ii , numLines ) ;
2010-06-06 18:00:24 +00:00
vv = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
vi = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
iv = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
ii = Create_N_3DArray < FDTD_FLOAT > ( numLines ) ;
2010-03-05 13:20:25 +00:00
}
2010-05-03 16:33:14 +00:00
void Operator : : InitExcitation ( )
{
delete Exc ;
Exc = new Excitation ( dT ) ;
}
2010-04-21 13:43:39 +00:00
void Operator : : Calc_ECOperatorPos ( int n , unsigned int * pos )
2010-04-09 13:51:37 +00:00
{
unsigned int i = MainOp - > SetPos ( pos [ 0 ] , pos [ 1 ] , pos [ 2 ] ) ;
2010-05-03 20:36:04 +00:00
if ( EC_C [ n ] [ i ] > 0 )
{
2010-10-20 05:26:16 +00:00
SetVV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , ( 1 - dT * EC_G [ n ] [ i ] / 2 / EC_C [ n ] [ i ] ) / ( 1 + dT * EC_G [ n ] [ i ] / 2 / EC_C [ n ] [ i ] ) ) ;
SetVI ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , ( dT / EC_C [ n ] [ i ] ) / ( 1 + dT * EC_G [ n ] [ i ] / 2 / EC_C [ n ] [ i ] ) ) ;
2010-05-03 20:36:04 +00:00
}
else
{
2010-10-20 05:26:16 +00:00
SetVV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetVI ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
2010-05-03 20:36:04 +00:00
}
if ( EC_L [ n ] [ i ] > 0 )
{
2010-10-20 05:26:16 +00:00
SetII ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , ( 1 - dT * EC_R [ n ] [ i ] / 2 / EC_L [ n ] [ i ] ) / ( 1 + dT * EC_R [ n ] [ i ] / 2 / EC_L [ n ] [ i ] ) ) ;
SetIV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , ( dT / EC_L [ n ] [ i ] ) / ( 1 + dT * EC_R [ n ] [ i ] / 2 / EC_L [ n ] [ i ] ) ) ;
2010-05-03 20:36:04 +00:00
}
else
{
2010-10-20 05:26:16 +00:00
SetII ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetIV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
2010-05-03 20:36:04 +00:00
}
2010-04-09 13:51:37 +00:00
}
2010-10-27 09:17:58 +00:00
int Operator : : CalcECOperator ( DebugFlags debugFlags )
2010-03-05 13:20:25 +00:00
{
2010-06-05 09:47:21 +00:00
Init_EC ( ) ;
2010-03-05 13:20:25 +00:00
if ( Calc_EC ( ) = = 0 )
return - 1 ;
2010-09-08 14:07:28 +00:00
m_InvaildTimestep = false ;
2010-09-02 20:12:03 +00:00
opt_dT = 0 ;
2010-08-19 08:12:09 +00:00
if ( dT > 0 )
{
double save_dT = dT ;
CalcTimestep ( ) ;
2010-09-02 20:12:03 +00:00
opt_dT = dT ;
2010-08-19 08:12:09 +00:00
if ( dT < save_dT )
2010-09-08 14:07:28 +00:00
{
2010-08-19 08:12:09 +00:00
cerr < < " Operator::CalcECOperator: Warning, forced timestep: " < < save_dT < < " s is larger than calculated timestep: " < < dT < < " s! It is not recommended using this timestep!! " < < endl ;
2010-09-08 14:07:28 +00:00
m_InvaildTimestep = true ;
}
2010-08-19 08:12:09 +00:00
dT = save_dT ;
}
else
CalcTimestep ( ) ;
2010-03-05 13:20:25 +00:00
InitOperator ( ) ;
unsigned int pos [ 3 ] ;
for ( int n = 0 ; n < 3 ; + + n )
{
for ( pos [ 0 ] = 0 ; pos [ 0 ] < numLines [ 0 ] ; + + pos [ 0 ] )
{
for ( pos [ 1 ] = 0 ; pos [ 1 ] < numLines [ 1 ] ; + + pos [ 1 ] )
{
for ( pos [ 2 ] = 0 ; pos [ 2 ] < numLines [ 2 ] ; + + pos [ 2 ] )
{
2010-04-09 13:51:37 +00:00
Calc_ECOperatorPos ( n , pos ) ;
2010-03-05 13:20:25 +00:00
}
}
}
}
2010-07-16 15:25:32 +00:00
//Apply PEC to all boundary's
2010-03-05 13:20:25 +00:00
bool PEC [ 6 ] = { 1 , 1 , 1 , 1 , 1 , 1 } ;
2010-10-04 09:52:39 +00:00
//make an exception for BC == -1
2010-07-16 15:25:32 +00:00
for ( int n = 0 ; n < 6 ; + + n )
2010-10-04 09:52:39 +00:00
if ( ( m_BC [ n ] = = - 1 ) )
2010-09-08 05:36:32 +00:00
PEC [ n ] = false ;
2010-03-05 13:20:25 +00:00
ApplyElectricBC ( PEC ) ;
CalcPEC ( ) ;
2010-04-28 20:01:03 +00:00
bool PMC [ 6 ] ;
for ( int n = 0 ; n < 6 ; + + n )
PMC [ n ] = m_BC [ n ] = = 1 ;
ApplyMagneticBC ( PMC ) ;
2010-10-04 09:52:39 +00:00
InitExcitation ( ) ;
if ( CalcFieldExcitation ( ) = = false ) return - 1 ;
//all information available for extension... create now...
for ( size_t n = 0 ; n < m_Op_exts . size ( ) ; + + n )
m_Op_exts . at ( n ) - > BuildExtension ( ) ;
2010-10-27 09:17:58 +00:00
if ( debugFlags & debugMaterial )
DumpMaterial2File ( " material_dump.vtk " ) ;
if ( debugFlags & debugOperator )
DumpOperator2File ( " operator_dump.vtk " ) ;
if ( debugFlags & debugPEC )
DumpPEC2File ( " PEC_dump.vtk " ) ;
2010-10-04 09:52:39 +00:00
//cleanup
for ( int n = 0 ; n < 3 ; + + n )
{
delete [ ] EC_C [ n ] ; EC_C [ n ] = NULL ;
delete [ ] EC_G [ n ] ; EC_G [ n ] = NULL ;
delete [ ] EC_L [ n ] ; EC_L [ n ] = NULL ;
delete [ ] EC_R [ n ] ; EC_R [ n ] = NULL ;
}
2010-03-05 13:20:25 +00:00
return 0 ;
}
void Operator : : ApplyElectricBC ( bool * dirs )
{
2010-11-03 14:04:33 +00:00
if ( ! dirs )
return ;
2010-03-05 13:20:25 +00:00
unsigned int pos [ 3 ] ;
for ( int n = 0 ; n < 3 ; + + n )
{
int nP = ( n + 1 ) % 3 ;
int nPP = ( n + 2 ) % 3 ;
for ( pos [ nP ] = 0 ; pos [ nP ] < numLines [ nP ] ; + + pos [ nP ] )
{
for ( pos [ nPP ] = 0 ; pos [ nPP ] < numLines [ nPP ] ; + + pos [ nPP ] )
{
2010-11-03 14:04:33 +00:00
if ( dirs [ 2 * n ] )
{
// set to PEC
pos [ n ] = 0 ;
SetVV ( nP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetVI ( nP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetVV ( nPP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetVI ( nPP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
}
if ( dirs [ 2 * n + 1 ] )
{
// set to PEC
pos [ n ] = numLines [ n ] - 1 ;
SetVV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ; // these are outside the FDTD-domain as defined by the main disc
SetVI ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ; // these are outside the FDTD-domain as defined by the main disc
SetVV ( nP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetVI ( nP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetVV ( nPP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetVI ( nPP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
}
2010-03-05 13:20:25 +00:00
}
}
}
}
void Operator : : ApplyMagneticBC ( bool * dirs )
{
2010-11-03 14:04:33 +00:00
if ( ! dirs )
return ;
2010-03-05 13:20:25 +00:00
unsigned int pos [ 3 ] ;
for ( int n = 0 ; n < 3 ; + + n )
{
int nP = ( n + 1 ) % 3 ;
int nPP = ( n + 2 ) % 3 ;
for ( pos [ nP ] = 0 ; pos [ nP ] < numLines [ nP ] ; + + pos [ nP ] )
{
for ( pos [ nPP ] = 0 ; pos [ nPP ] < numLines [ nPP ] ; + + pos [ nPP ] )
{
2010-11-03 14:04:33 +00:00
if ( dirs [ 2 * n ] )
{
// set to PMC
pos [ n ] = 0 ;
SetII ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetIV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetII ( nP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetIV ( nP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetII ( nPP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetIV ( nPP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
}
if ( dirs [ 2 * n + 1 ] )
{
// set to PMC
pos [ n ] = numLines [ n ] - 2 ;
SetII ( nP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetIV ( nP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetII ( nPP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetIV ( nPP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
}
2010-10-04 09:52:39 +00:00
2010-10-20 05:26:16 +00:00
// the last current lines are outside the FDTD domain and cannot be iterated by the FDTD engine
2010-11-03 14:04:33 +00:00
pos [ n ] = numLines [ n ] - 1 ;
2010-10-20 05:26:16 +00:00
SetII ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetIV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetII ( nP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetIV ( nP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetII ( nPP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetIV ( nPP , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
2010-03-05 13:20:25 +00:00
}
}
}
}
2010-09-02 13:35:13 +00:00
bool Operator : : Calc_ECPos ( int ny , const unsigned int * pos , double * EC ) const
{
double EffMat [ 4 ] ;
Calc_EffMatPos ( ny , pos , EffMat ) ;
double delta = GetEdgeLength ( ny , pos ) ;
double area = GetEdgeArea ( ny , pos ) ;
// if (isnan(EffMat[0]))
// {
// cerr << ny << " " << pos[0] << " " << pos[1] << " " << pos[2] << " : " << EffMat[0] << endl;
// }
if ( delta )
{
EC [ 0 ] = EffMat [ 0 ] * area / delta ;
EC [ 1 ] = EffMat [ 1 ] * area / delta ;
}
else
{
EC [ 0 ] = 0 ;
EC [ 1 ] = 0 ;
}
delta = GetEdgeLength ( ny , pos , true ) ;
area = GetEdgeArea ( ny , pos , true ) ;
if ( delta )
{
EC [ 2 ] = EffMat [ 2 ] * area / delta ;
EC [ 3 ] = EffMat [ 3 ] * area / delta ;
}
else
{
EC [ 2 ] = 0 ;
EC [ 3 ] = 0 ;
}
2010-03-05 13:20:25 +00:00
2010-09-02 13:35:13 +00:00
return true ;
}
bool Operator : : Calc_EffMatPos ( int ny , const unsigned int * pos , double * EffMat ) const
2010-03-05 13:20:25 +00:00
{
2010-09-02 13:35:13 +00:00
int n = ny ;
2010-03-05 13:20:25 +00:00
double coord [ 3 ] ;
double shiftCoord [ 3 ] ;
int nP = ( n + 1 ) % 3 ;
int nPP = ( n + 2 ) % 3 ;
coord [ 0 ] = discLines [ 0 ] [ pos [ 0 ] ] ;
coord [ 1 ] = discLines [ 1 ] [ pos [ 1 ] ] ;
coord [ 2 ] = discLines [ 2 ] [ pos [ 2 ] ] ;
double delta = MainOp - > GetIndexDelta ( n , pos [ n ] ) ;
double deltaP = MainOp - > GetIndexDelta ( nP , pos [ nP ] ) ;
double deltaPP = MainOp - > GetIndexDelta ( nPP , pos [ nPP ] ) ;
double delta_M = MainOp - > GetIndexDelta ( n , pos [ n ] - 1 ) ;
double deltaP_M = MainOp - > GetIndexDelta ( nP , pos [ nP ] - 1 ) ;
double deltaPP_M = MainOp - > GetIndexDelta ( nPP , pos [ nPP ] - 1 ) ;
2010-09-02 13:35:13 +00:00
int loc_pos [ 3 ] = { pos [ 0 ] , pos [ 1 ] , pos [ 2 ] } ;
double A_n ;
double area = 0 ;
2010-03-05 13:20:25 +00:00
//******************************* epsilon,kappa averaging *****************************//
//shift up-right
shiftCoord [ n ] = coord [ n ] + delta * 0.5 ;
shiftCoord [ nP ] = coord [ nP ] + deltaP * 0.25 ;
shiftCoord [ nPP ] = coord [ nPP ] + deltaPP * 0.25 ;
2010-09-17 12:51:07 +00:00
A_n = GetNodeArea ( ny , ( unsigned int * ) loc_pos , true ) ;
2010-09-02 13:35:13 +00:00
// {
// cerr << ny << " " << pos[0] << " " << pos[1] << " " << pos[2] << ": " << A_n << endl;
// exit(0);
// }
2010-10-19 14:54:35 +00:00
CSProperties * prop = CSX - > GetPropertyByCoordPriority ( shiftCoord , CSProperties : : MATERIAL , true ) ;
2010-03-05 13:20:25 +00:00
if ( prop )
{
CSPropMaterial * mat = prop - > ToMaterial ( ) ;
2010-09-02 13:35:13 +00:00
EffMat [ 0 ] = mat - > GetEpsilonWeighted ( n , shiftCoord ) * A_n ;
EffMat [ 1 ] = mat - > GetKappaWeighted ( n , shiftCoord ) * A_n ;
2010-03-05 13:20:25 +00:00
}
else
{
2010-09-02 13:35:13 +00:00
EffMat [ 0 ] = 1 * A_n ;
EffMat [ 1 ] = 0 ;
2010-03-05 13:20:25 +00:00
}
2010-09-02 13:35:13 +00:00
area + = A_n ;
2010-03-05 13:20:25 +00:00
//shift up-left
shiftCoord [ n ] = coord [ n ] + delta * 0.5 ;
shiftCoord [ nP ] = coord [ nP ] - deltaP_M * 0.25 ;
shiftCoord [ nPP ] = coord [ nPP ] + deltaPP * 0.25 ;
2010-09-02 13:35:13 +00:00
- - loc_pos [ nP ] ;
2010-09-17 12:51:07 +00:00
A_n = GetNodeArea ( ny , ( unsigned int * ) loc_pos , true ) ;
2010-09-02 13:35:13 +00:00
// cerr << A_n << endl;
2010-10-19 14:54:35 +00:00
prop = CSX - > GetPropertyByCoordPriority ( shiftCoord , CSProperties : : MATERIAL , true ) ;
2010-03-05 13:20:25 +00:00
if ( prop )
{
CSPropMaterial * mat = prop - > ToMaterial ( ) ;
2010-09-02 13:35:13 +00:00
EffMat [ 0 ] + = mat - > GetEpsilonWeighted ( n , shiftCoord ) * A_n ;
EffMat [ 1 ] + = mat - > GetKappaWeighted ( n , shiftCoord ) * A_n ;
2010-03-05 13:20:25 +00:00
}
else
{
2010-09-02 13:35:13 +00:00
EffMat [ 0 ] + = 1 * A_n ;
EffMat [ 1 ] + = 0 ;
2010-03-05 13:20:25 +00:00
}
2010-09-02 13:35:13 +00:00
area + = A_n ;
2010-03-05 13:20:25 +00:00
//shift down-right
shiftCoord [ n ] = coord [ n ] + delta * 0.5 ;
shiftCoord [ nP ] = coord [ nP ] + deltaP * 0.25 ;
shiftCoord [ nPP ] = coord [ nPP ] - deltaPP_M * 0.25 ;
2010-09-02 13:35:13 +00:00
+ + loc_pos [ nP ] ;
- - loc_pos [ nPP ] ;
2010-09-17 12:51:07 +00:00
A_n = GetNodeArea ( ny , ( unsigned int * ) loc_pos , true ) ;
2010-10-19 14:54:35 +00:00
prop = CSX - > GetPropertyByCoordPriority ( shiftCoord , CSProperties : : MATERIAL , true ) ;
2010-03-05 13:20:25 +00:00
if ( prop )
{
CSPropMaterial * mat = prop - > ToMaterial ( ) ;
2010-09-02 13:35:13 +00:00
EffMat [ 0 ] + = mat - > GetEpsilonWeighted ( n , shiftCoord ) * A_n ;
EffMat [ 1 ] + = mat - > GetKappaWeighted ( n , shiftCoord ) * A_n ;
2010-03-05 13:20:25 +00:00
}
else
{
2010-09-02 13:35:13 +00:00
EffMat [ 0 ] + = 1 * A_n ;
EffMat [ 1 ] + = 0 ;
2010-03-05 13:20:25 +00:00
}
2010-09-02 13:35:13 +00:00
area + = A_n ;
2010-03-05 13:20:25 +00:00
//shift down-left
shiftCoord [ n ] = coord [ n ] + delta * 0.5 ;
shiftCoord [ nP ] = coord [ nP ] - deltaP_M * 0.25 ;
shiftCoord [ nPP ] = coord [ nPP ] - deltaPP_M * 0.25 ;
2010-09-02 13:35:13 +00:00
- - loc_pos [ nP ] ;
2010-09-17 12:51:07 +00:00
A_n = GetNodeArea ( ny , ( unsigned int * ) loc_pos , true ) ;
2010-10-19 14:54:35 +00:00
prop = CSX - > GetPropertyByCoordPriority ( shiftCoord , CSProperties : : MATERIAL , true ) ;
2010-03-05 13:20:25 +00:00
if ( prop )
{
CSPropMaterial * mat = prop - > ToMaterial ( ) ;
2010-09-02 13:35:13 +00:00
EffMat [ 0 ] + = mat - > GetEpsilonWeighted ( n , shiftCoord ) * A_n ;
EffMat [ 1 ] + = mat - > GetKappaWeighted ( n , shiftCoord ) * A_n ;
2010-03-05 13:20:25 +00:00
}
else
{
2010-09-02 13:35:13 +00:00
EffMat [ 0 ] + = 1 * A_n ;
EffMat [ 1 ] + = 0 ;
2010-03-05 13:20:25 +00:00
}
2010-09-02 13:35:13 +00:00
area + = A_n ;
2010-03-05 13:20:25 +00:00
2010-09-02 13:35:13 +00:00
EffMat [ 0 ] * = __EPS0__ / area ;
EffMat [ 1 ] / = area ;
2010-03-05 13:20:25 +00:00
//******************************* mu,sigma averaging *****************************//
2010-09-02 13:35:13 +00:00
loc_pos [ 0 ] = pos [ 0 ] ; loc_pos [ 1 ] = pos [ 1 ] ; loc_pos [ 2 ] = pos [ 2 ] ;
double length = 0 ;
2010-03-05 13:20:25 +00:00
//shift down
shiftCoord [ n ] = coord [ n ] - delta_M * 0.25 ;
shiftCoord [ nP ] = coord [ nP ] + deltaP * 0.5 ;
shiftCoord [ nPP ] = coord [ nPP ] + deltaPP * 0.5 ;
2010-09-02 13:35:13 +00:00
- - loc_pos [ n ] ;
2010-09-17 12:51:07 +00:00
double delta_ny = GetNodeWidth ( n , ( unsigned int * ) loc_pos , true ) ;
2010-10-19 14:54:35 +00:00
prop = CSX - > GetPropertyByCoordPriority ( shiftCoord , CSProperties : : MATERIAL , true ) ;
2010-03-05 13:20:25 +00:00
if ( prop )
{
CSPropMaterial * mat = prop - > ToMaterial ( ) ;
2010-09-02 13:35:13 +00:00
EffMat [ 2 ] = delta_ny / mat - > GetMueWeighted ( n , shiftCoord ) ;
if ( mat - > GetSigmaWeighted ( n , shiftCoord ) )
EffMat [ 3 ] = delta_ny / mat - > GetSigmaWeighted ( n , shiftCoord ) ;
2010-03-05 13:20:25 +00:00
else
2010-09-02 13:35:13 +00:00
EffMat [ 3 ] = 0 ;
2010-03-05 13:20:25 +00:00
}
else
{
2010-09-02 13:35:13 +00:00
EffMat [ 2 ] = delta_ny ;
EffMat [ 3 ] = 0 ;
2010-03-05 13:20:25 +00:00
}
2010-09-02 13:35:13 +00:00
length = delta_ny ;
2010-03-05 13:20:25 +00:00
//shift up
shiftCoord [ n ] = coord [ n ] + delta * 0.25 ;
shiftCoord [ nP ] = coord [ nP ] + deltaP * 0.5 ;
shiftCoord [ nPP ] = coord [ nPP ] + deltaPP * 0.5 ;
2010-09-02 13:35:13 +00:00
+ + loc_pos [ n ] ;
2010-09-17 12:51:07 +00:00
delta_ny = GetNodeWidth ( n , ( unsigned int * ) loc_pos , true ) ;
2010-10-19 14:54:35 +00:00
prop = CSX - > GetPropertyByCoordPriority ( shiftCoord , CSProperties : : MATERIAL , true ) ;
2010-03-05 13:20:25 +00:00
if ( prop )
{
CSPropMaterial * mat = prop - > ToMaterial ( ) ;
2010-09-02 13:35:13 +00:00
EffMat [ 2 ] + = delta_ny / mat - > GetMueWeighted ( n , shiftCoord ) ;
2010-03-05 13:20:25 +00:00
if ( mat - > GetSigmaWeighted ( n , shiftCoord ) )
2010-09-02 13:35:13 +00:00
EffMat [ 3 ] + = delta_ny / mat - > GetSigmaWeighted ( n , shiftCoord ) ;
2010-03-05 13:20:25 +00:00
else
2010-09-02 13:35:13 +00:00
EffMat [ 3 ] = 0 ;
2010-03-05 13:20:25 +00:00
}
else
{
2010-09-02 13:35:13 +00:00
EffMat [ 2 ] + = 1 * delta_ny ;
EffMat [ 3 ] = 0 ;
2010-03-05 13:20:25 +00:00
}
2010-09-02 13:35:13 +00:00
length + = delta_ny ;
2010-03-05 13:20:25 +00:00
2010-09-02 13:35:13 +00:00
EffMat [ 2 ] = length * __MUE0__ / EffMat [ 2 ] ;
if ( EffMat [ 3 ] ) EffMat [ 3 ] = length / EffMat [ 3 ] ;
2010-03-05 13:20:25 +00:00
2010-09-02 13:35:13 +00:00
for ( int n = 0 ; n < 4 ; + + n )
if ( isnan ( EffMat [ n ] ) | | isinf ( EffMat [ n ] ) )
{
cerr < < " Operator::Calc_EffMatPos: An effective material parameter is not a valid result, this should NOT have happend... exit... " < < endl ;
exit ( 0 ) ;
}
2010-03-05 13:20:25 +00:00
return true ;
}
2010-06-05 09:47:21 +00:00
void Operator : : Init_EC ( )
2010-03-05 13:20:25 +00:00
{
for ( int n = 0 ; n < 3 ; + + n )
{
//init x-cell-array
delete [ ] EC_C [ n ] ;
delete [ ] EC_G [ n ] ;
delete [ ] EC_L [ n ] ;
delete [ ] EC_R [ n ] ;
EC_C [ n ] = new double [ MainOp - > GetSize ( ) ] ;
EC_G [ n ] = new double [ MainOp - > GetSize ( ) ] ;
EC_L [ n ] = new double [ MainOp - > GetSize ( ) ] ;
EC_R [ n ] = new double [ MainOp - > GetSize ( ) ] ;
for ( unsigned int i = 0 ; i < MainOp - > GetSize ( ) ; i + + ) //init all
{
EC_C [ n ] [ i ] = 0 ;
EC_G [ n ] [ i ] = 0 ;
EC_L [ n ] [ i ] = 0 ;
EC_R [ n ] [ i ] = 0 ;
}
2010-06-05 09:47:21 +00:00
}
}
bool Operator : : Calc_EC ( )
{
if ( CSX = = NULL ) { cerr < < " CartOperator::Calc_EC: CSX not given or invalid!!! " < < endl ; return false ; }
unsigned int ipos ;
unsigned int pos [ 3 ] ;
double inEC [ 4 ] ;
for ( int n = 0 ; n < 3 ; + + n )
{
2010-03-05 13:20:25 +00:00
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 ] )
{
Calc_ECPos ( n , pos , inEC ) ;
ipos = MainOp - > SetPos ( pos [ 0 ] , pos [ 1 ] , pos [ 2 ] ) ;
EC_C [ n ] [ ipos ] = inEC [ 0 ] ;
EC_G [ n ] [ ipos ] = inEC [ 1 ] ;
EC_L [ n ] [ ipos ] = inEC [ 2 ] ;
EC_R [ n ] [ ipos ] = inEC [ 3 ] ;
}
}
}
}
return true ;
}
2010-06-22 10:49:51 +00:00
2010-03-05 13:20:25 +00:00
double Operator : : CalcTimestep ( )
{
2010-07-13 13:37:56 +00:00
# if 1 //use the new timestep-calc (1) or the old one (0)
return CalcTimestep_Var3 ( ) ; //the biggest one for cartesian meshes
2010-06-22 10:49:51 +00:00
# else
return CalcTimestep_Var1 ( ) ;
# endif
}
////Berechnung nach Andreas Rennings Dissertation 2008, Seite 66, Formel 4.52
double Operator : : CalcTimestep_Var1 ( )
{
2010-07-13 13:37:56 +00:00
m_Used_TS_Name = string ( " Rennings_1 " ) ;
// cout << "Operator::CalcTimestep(): Using timestep algorithm by Andreas Rennings, Dissertation @ University Duisburg-Essen, 2008, pp. 66, eq. 4.52" << endl;
2010-03-05 13:20:25 +00:00
dT = 1e200 ;
double newT ;
unsigned int pos [ 3 ] ;
unsigned int ipos ;
unsigned int ipos_PM ;
unsigned int ipos_PPM ;
MainOp - > SetReflection2Cell ( ) ;
for ( int n = 0 ; n < 3 ; + + n )
{
int nP = ( n + 1 ) % 3 ;
int nPP = ( n + 2 ) % 3 ;
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 ] )
{
ipos = MainOp - > SetPos ( pos [ 0 ] , pos [ 1 ] , pos [ 2 ] ) ;
ipos_PM = MainOp - > Shift ( nP , - 1 ) ;
MainOp - > ResetShift ( ) ;
ipos_PPM = MainOp - > Shift ( nPP , - 1 ) ;
MainOp - > ResetShift ( ) ;
newT = 2 / sqrt ( ( 4 / EC_L [ nP ] [ ipos ] + 4 / EC_L [ nP ] [ ipos_PPM ] + 4 / EC_L [ nPP ] [ ipos ] + 4 / EC_L [ nPP ] [ ipos_PM ] ) / EC_C [ n ] [ ipos ] ) ;
2010-04-13 16:40:43 +00:00
if ( ( newT < dT ) & & ( newT > 0.0 ) ) dT = newT ;
2010-03-05 13:20:25 +00:00
}
}
}
}
2010-04-13 16:40:43 +00:00
if ( dT = = 0 )
{
cerr < < " Operator::CalcTimestep: Timestep is zero... this is not supposed to happen!!! exit! " < < endl ;
exit ( 3 ) ;
}
2010-03-30 06:19:41 +00:00
// cerr << "Operator Timestep: " << dT << endl;
2010-03-05 13:20:25 +00:00
return 0 ;
}
2010-06-18 10:37:37 +00:00
double min ( double * val , unsigned int count )
{
if ( count = = 0 )
return 0.0 ;
double min = val [ 0 ] ;
for ( unsigned int n = 1 ; n < count ; + + n )
if ( val [ n ] < min )
min = val [ n ] ;
return min ;
}
//Berechnung nach Andreas Rennings Dissertation 2008, Seite 76 ff, Formel 4.77 ff
2010-06-22 10:49:51 +00:00
double Operator : : CalcTimestep_Var3 ( )
2010-06-18 10:37:37 +00:00
{
dT = 1e200 ;
2010-07-13 13:37:56 +00:00
m_Used_TS_Name = string ( " Rennings_2 " ) ;
// cout << "Operator::CalcTimestep(): Using timestep algorithm by Andreas Rennings, Dissertation @ University Duisburg-Essen, 2008, pp. 76, eq. 4.77 ff." << endl;
2010-06-18 10:37:37 +00:00
double newT ;
unsigned int pos [ 3 ] ;
unsigned int ipos ;
double w_total = 0 ;
double wqp = 0 , wt1 = 0 , wt2 = 0 ;
double wt_4 [ 4 ] = { 0 , 0 , 0 , 0 } ;
MainOp - > SetReflection2Cell ( ) ;
for ( int n = 0 ; n < 3 ; + + n )
{
int nP = ( n + 1 ) % 3 ;
int nPP = ( n + 2 ) % 3 ;
2010-06-22 10:49:51 +00:00
for ( pos [ 2 ] = 0 ; pos [ 2 ] < numLines [ 2 ] ; + + pos [ 2 ] )
2010-06-18 10:37:37 +00:00
{
2010-06-22 10:49:51 +00:00
for ( pos [ 1 ] = 0 ; pos [ 1 ] < numLines [ 1 ] ; + + pos [ 1 ] )
2010-06-18 10:37:37 +00:00
{
2010-06-22 10:49:51 +00:00
for ( pos [ 0 ] = 0 ; pos [ 0 ] < numLines [ 0 ] ; + + pos [ 0 ] )
2010-06-18 10:37:37 +00:00
{
MainOp - > ResetShift ( ) ;
ipos = MainOp - > SetPos ( pos [ 0 ] , pos [ 1 ] , pos [ 2 ] ) ;
wqp = 1 / ( EC_L [ nPP ] [ ipos ] * EC_C [ n ] [ MainOp - > GetShiftedPos ( nP , 1 ) ] ) + 1 / ( EC_L [ nPP ] [ ipos ] * EC_C [ n ] [ ipos ] ) ;
wqp + = 1 / ( EC_L [ nP ] [ ipos ] * EC_C [ n ] [ MainOp - > GetShiftedPos ( nPP , 1 ) ] ) + 1 / ( EC_L [ nP ] [ ipos ] * EC_C [ n ] [ ipos ] ) ;
ipos = MainOp - > Shift ( nP , - 1 ) ;
wqp + = 1 / ( EC_L [ nPP ] [ ipos ] * EC_C [ n ] [ MainOp - > GetShiftedPos ( nP , 1 ) ] ) + 1 / ( EC_L [ nPP ] [ ipos ] * EC_C [ n ] [ ipos ] ) ;
ipos = MainOp - > Shift ( nPP , - 1 ) ;
wqp + = 1 / ( EC_L [ nP ] [ ipos ] * EC_C [ n ] [ MainOp - > GetShiftedPos ( nPP , 1 ) ] ) + 1 / ( EC_L [ nP ] [ ipos ] * EC_C [ n ] [ ipos ] ) ;
MainOp - > ResetShift ( ) ;
ipos = MainOp - > SetPos ( pos [ 0 ] , pos [ 1 ] , pos [ 2 ] ) ;
wt_4 [ 0 ] = 1 / ( EC_L [ nPP ] [ ipos ] * EC_C [ nP ] [ ipos ] ) ;
wt_4 [ 1 ] = 1 / ( EC_L [ nPP ] [ MainOp - > GetShiftedPos ( nP , - 1 ) ] * EC_C [ nP ] [ ipos ] ) ;
wt_4 [ 2 ] = 1 / ( EC_L [ nP ] [ ipos ] * EC_C [ nPP ] [ ipos ] ) ;
wt_4 [ 3 ] = 1 / ( EC_L [ nP ] [ MainOp - > GetShiftedPos ( nPP , - 1 ) ] * EC_C [ nPP ] [ ipos ] ) ;
wt1 = wt_4 [ 0 ] + wt_4 [ 1 ] + wt_4 [ 2 ] + wt_4 [ 3 ] - 2 * min ( wt_4 , 4 ) ;
MainOp - > ResetShift ( ) ;
ipos = MainOp - > SetPos ( pos [ 0 ] , pos [ 1 ] , pos [ 2 ] ) ;
wt_4 [ 0 ] = 1 / ( EC_L [ nPP ] [ ipos ] * EC_C [ nP ] [ MainOp - > GetShiftedPos ( n , 1 ) ] ) ;
wt_4 [ 1 ] = 1 / ( EC_L [ nPP ] [ MainOp - > GetShiftedPos ( nP , - 1 ) ] * EC_C [ nP ] [ MainOp - > GetShiftedPos ( n , 1 ) ] ) ;
wt_4 [ 2 ] = 1 / ( EC_L [ nP ] [ ipos ] * EC_C [ nPP ] [ MainOp - > GetShiftedPos ( n , 1 ) ] ) ;
wt_4 [ 3 ] = 1 / ( EC_L [ nP ] [ MainOp - > GetShiftedPos ( nPP , - 1 ) ] * EC_C [ nPP ] [ MainOp - > GetShiftedPos ( n , 1 ) ] ) ;
wt2 = wt_4 [ 0 ] + wt_4 [ 1 ] + wt_4 [ 2 ] + wt_4 [ 3 ] - 2 * min ( wt_4 , 4 ) ;
w_total = wqp + wt1 + wt2 ;
newT = 2 / sqrt ( w_total ) ;
if ( ( newT < dT ) & & ( newT > 0.0 ) )
dT = newT ;
}
}
}
}
if ( dT = = 0 )
{
cerr < < " Operator::CalcTimestep: Timestep is zero... this is not supposed to happen!!! exit! " < < endl ;
exit ( 3 ) ;
}
// cerr << "Operator Timestep: " << dT << endl;
return 0 ;
}
2010-04-28 20:42:55 +00:00
bool Operator : : CalcFieldExcitation ( )
2010-03-05 13:20:25 +00:00
{
2010-05-03 16:33:14 +00:00
if ( dT = = 0 )
return false ;
if ( Exc = = 0 )
return false ;
2010-03-29 08:12:38 +00:00
unsigned int pos [ 3 ] ;
2010-03-16 20:39:37 +00:00
double delta [ 3 ] ;
double amp = 0 ;
2010-03-05 13:20:25 +00:00
2010-04-28 22:28:07 +00:00
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 ] ;
2010-05-20 13:17:30 +00:00
vector < CSProperties * > vec_prop = CSX - > GetPropertyByType ( CSProperties : : ELECTRODE ) ;
if ( vec_prop . size ( ) = = 0 )
{
cerr < < " Operator::CalcFieldExcitation: Warning, no excitation properties found " < < endl ;
return false ;
}
CSPropElectrode * elec = NULL ;
CSProperties * prop = NULL ;
int priority = 0 ;
2010-03-29 08:12:38 +00:00
for ( pos [ 2 ] = 0 ; pos [ 2 ] < numLines [ 2 ] ; + + pos [ 2 ] )
2010-03-05 13:20:25 +00:00
{
2010-03-16 20:39:37 +00:00
delta [ 2 ] = fabs ( MainOp - > GetIndexDelta ( 2 , pos [ 2 ] ) ) ;
2010-03-29 08:12:38 +00:00
for ( pos [ 1 ] = 0 ; pos [ 1 ] < numLines [ 1 ] ; + + pos [ 1 ] )
2010-03-05 13:20:25 +00:00
{
2010-03-16 20:39:37 +00:00
delta [ 1 ] = fabs ( MainOp - > GetIndexDelta ( 1 , pos [ 1 ] ) ) ;
2010-03-29 08:12:38 +00:00
for ( pos [ 0 ] = 0 ; pos [ 0 ] < numLines [ 0 ] ; + + pos [ 0 ] )
2010-03-05 13:20:25 +00:00
{
2010-03-16 20:39:37 +00:00
delta [ 0 ] = fabs ( MainOp - > GetIndexDelta ( 0 , pos [ 0 ] ) ) ;
2010-04-28 22:28:07 +00:00
//electric field excite
2010-03-16 20:39:37 +00:00
for ( int n = 0 ; n < 3 ; + + n )
2010-03-05 13:20:25 +00:00
{
2010-04-28 20:42:55 +00:00
volt_coord [ 0 ] = discLines [ 0 ] [ pos [ 0 ] ] ;
volt_coord [ 1 ] = discLines [ 1 ] [ pos [ 1 ] ] ;
volt_coord [ 2 ] = discLines [ 2 ] [ pos [ 2 ] ] ;
volt_coord [ n ] + = delta [ n ] * 0.5 ;
2010-05-20 13:17:30 +00:00
for ( size_t p = 0 ; p < vec_prop . size ( ) ; + + p )
2010-03-16 20:39:37 +00:00
{
2010-05-20 13:17:30 +00:00
prop = vec_prop . at ( p ) ;
elec = prop - > ToElectrode ( ) ;
2010-10-19 14:54:35 +00:00
if ( elec = = NULL )
continue ;
if ( prop - > CheckCoordInPrimitive ( volt_coord , priority , true ) )
2010-03-05 13:20:25 +00:00
{
2010-04-28 20:42:55 +00:00
if ( ( elec - > GetActiveDir ( n ) ) & & ( ( elec - > GetExcitType ( ) = = 0 ) | | ( elec - > GetExcitType ( ) = = 1 ) ) ) //&& (pos[n]<numLines[n]-1))
2010-03-05 13:20:25 +00:00
{
2010-04-28 20:42:55 +00:00
amp = elec - > GetWeightedExcitation ( n , volt_coord ) * GetMeshDelta ( n , pos ) ; // delta[n]*gridDelta;
2010-03-16 20:39:37 +00:00
if ( amp ! = 0 )
{
2010-04-28 22:28:07 +00:00
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 ] ) ;
2010-03-16 20:39:37 +00:00
}
if ( elec - > GetExcitType ( ) = = 1 ) //hard excite
2010-03-11 14:03:24 +00:00
{
2010-10-20 05:26:16 +00:00
SetVV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetVI ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
2010-03-11 14:03:24 +00:00
}
2010-03-05 13:20:25 +00:00
}
}
2010-03-16 20:39:37 +00:00
}
2010-03-05 13:20:25 +00:00
}
2010-04-28 22:28:07 +00:00
//magnetic field excite
for ( int n = 0 ; n < 3 ; + + n )
{
2010-08-03 15:12:33 +00:00
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
2010-04-28 22:28:07 +00:00
int nP = ( n + 1 ) % 3 ;
int nPP = ( n + 2 ) % 3 ;
curr_coord [ 0 ] = discLines [ 0 ] [ pos [ 0 ] ] ;
curr_coord [ 1 ] = discLines [ 1 ] [ pos [ 1 ] ] ;
curr_coord [ 2 ] = discLines [ 2 ] [ pos [ 2 ] ] ;
curr_coord [ nP ] + = delta [ nP ] * 0.5 ;
curr_coord [ nPP ] + = delta [ nPP ] * 0.5 ;
2010-05-20 13:17:30 +00:00
for ( size_t p = 0 ; p < vec_prop . size ( ) ; + + p )
2010-04-28 22:28:07 +00:00
{
2010-05-20 13:17:30 +00:00
prop = vec_prop . at ( p ) ;
elec = prop - > ToElectrode ( ) ;
2010-10-19 14:54:35 +00:00
if ( elec = = NULL )
continue ;
if ( prop - > CheckCoordInPrimitive ( curr_coord , priority , true ) )
2010-04-28 22:28:07 +00:00
{
2010-08-03 15:12:33 +00:00
if ( ( elec - > GetActiveDir ( n ) ) & & ( ( elec - > GetExcitType ( ) = = 2 ) | | ( elec - > GetExcitType ( ) = = 3 ) ) )
2010-04-28 22:28:07 +00:00
{
amp = elec - > GetWeightedExcitation ( n , curr_coord ) * GetMeshDelta ( 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
{
2010-10-20 05:26:16 +00:00
SetII ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetIV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
2010-04-28 22:28:07 +00:00
}
}
}
}
}
2010-03-05 13:20:25 +00:00
}
}
}
2010-03-16 20:39:37 +00:00
2010-03-27 22:05:45 +00:00
//special treatment for primitives of type curve (treated as wires) see also Calc_PEC
double p1 [ 3 ] ;
double p2 [ 3 ] ;
double deltaN = 0.0 ;
struct Grid_Path path ;
for ( size_t p = 0 ; p < vec_prop . size ( ) ; + + p )
{
prop = vec_prop . at ( p ) ;
elec = prop - > ToElectrode ( ) ;
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 ) ;
curv - > GetPoint ( i , p2 ) ;
path = FindPath ( p1 , p2 ) ;
2010-05-29 15:40:18 +00:00
if ( path . dir . size ( ) > 0 )
prim - > SetPrimitiveUsed ( true ) ;
2010-03-27 22:05:45 +00:00
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 ) ;
MainOp - > SetPos ( pos [ 0 ] , pos [ 1 ] , pos [ 2 ] ) ;
deltaN = fabs ( MainOp - > GetIndexDelta ( n , pos [ n ] ) ) ;
2010-04-28 20:42:55 +00:00
volt_coord [ 0 ] = discLines [ 0 ] [ pos [ 0 ] ] ;
volt_coord [ 1 ] = discLines [ 1 ] [ pos [ 1 ] ] ;
volt_coord [ 2 ] = discLines [ 2 ] [ pos [ 2 ] ] ;
volt_coord [ n ] + = 0.5 * deltaN ;
2010-03-27 22:05:45 +00:00
// cerr << n << " " << coord[0] << " " << coord[1] << " " << coord[2] << endl;
if ( elec ! = NULL )
{
2010-04-28 20:42:55 +00:00
if ( ( elec - > GetActiveDir ( n ) ) & & ( pos [ n ] < numLines [ n ] - 1 ) & & ( ( elec - > GetExcitType ( ) = = 0 ) | | ( elec - > GetExcitType ( ) = = 1 ) ) )
2010-03-27 22:05:45 +00:00
{
2010-04-28 20:42:55 +00:00
amp = elec - > GetWeightedExcitation ( n , volt_coord ) * deltaN * gridDelta ;
2010-03-27 22:05:45 +00:00
if ( amp ! = 0 )
{
2010-04-28 22:28:07 +00:00
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 ] ) ;
2010-03-27 22:05:45 +00:00
}
if ( elec - > GetExcitType ( ) = = 1 ) //hard excite
{
2010-10-20 05:26:16 +00:00
SetVV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetVI ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
2010-03-27 22:05:45 +00:00
}
}
}
}
}
}
}
}
2010-04-28 22:28:07 +00:00
// set voltage excitations
2010-05-03 16:33:14 +00:00
Exc - > setupVoltageExcitation ( volt_vIndex , volt_vExcit , volt_vDelay , volt_vDir ) ;
2010-04-28 22:28:07 +00:00
// set current excitations
2010-05-03 16:33:14 +00:00
Exc - > setupCurrentExcitation ( curr_vIndex , curr_vExcit , curr_vDelay , curr_vDir ) ;
2010-04-28 22:28:07 +00:00
2010-03-05 13:20:25 +00:00
return true ;
}
2010-06-05 22:53:05 +00:00
bool Operator : : CalcPEC ( )
{
m_Nr_PEC [ 0 ] = 0 ; m_Nr_PEC [ 1 ] = 0 ; m_Nr_PEC [ 2 ] = 0 ;
2010-05-03 16:33:14 +00:00
2010-06-05 23:47:32 +00:00
CalcPEC_Range ( 0 , numLines [ 0 ] - 1 , m_Nr_PEC ) ;
2010-05-03 16:33:14 +00:00
2010-06-05 22:53:05 +00:00
CalcPEC_Curves ( ) ;
return true ;
}
2010-06-05 23:47:32 +00:00
void Operator : : CalcPEC_Range ( unsigned int startX , unsigned int stopX , unsigned int * counter )
2010-03-05 13:20:25 +00:00
{
double coord [ 3 ] ;
double delta ;
2010-06-05 22:53:05 +00:00
unsigned int pos [ 3 ] ;
for ( pos [ 0 ] = startX ; pos [ 0 ] < = stopX ; + + pos [ 0 ] )
2010-03-05 13:20:25 +00:00
{
2010-06-05 22:53:05 +00:00
for ( pos [ 1 ] = 0 ; pos [ 1 ] < numLines [ 1 ] ; + + pos [ 1 ] )
2010-03-05 13:20:25 +00:00
{
2010-06-05 22:53:05 +00:00
for ( pos [ 2 ] = 0 ; pos [ 2 ] < numLines [ 2 ] ; + + pos [ 2 ] )
2010-03-05 13:20:25 +00:00
{
2010-06-05 22:53:05 +00:00
for ( int n = 0 ; n < 3 ; + + n )
2010-03-05 13:20:25 +00:00
{
coord [ 0 ] = discLines [ 0 ] [ pos [ 0 ] ] ;
2010-03-12 07:38:07 +00:00
coord [ 1 ] = discLines [ 1 ] [ pos [ 1 ] ] ;
coord [ 2 ] = discLines [ 2 ] [ pos [ 2 ] ] ;
2010-03-05 13:20:25 +00:00
delta = MainOp - > GetIndexDelta ( n , pos [ n ] ) ;
coord [ n ] = discLines [ n ] [ pos [ n ] ] + delta * 0.5 ;
2010-10-19 14:54:35 +00:00
CSProperties * prop = CSX - > GetPropertyByCoordPriority ( coord , ( CSProperties : : PropertyType ) ( CSProperties : : MATERIAL | CSProperties : : METAL ) , true ) ;
2010-03-05 13:20:25 +00:00
if ( prop )
{
if ( prop - > GetType ( ) = = CSProperties : : METAL ) //set to PEC
{
2010-10-20 05:26:16 +00:00
SetVV ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
SetVI ( n , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] , 0 ) ;
2010-06-05 23:47:32 +00:00
+ + counter [ n ] ;
2010-06-05 22:53:05 +00:00
// cerr << "CartOperator::CalcPEC: PEC found at " << pos[0] << " ; " << pos[1] << " ; " << pos[2] << endl;
2010-03-05 13:20:25 +00:00
}
}
}
}
}
}
2010-06-05 22:53:05 +00:00
}
2010-03-22 07:19:17 +00:00
2010-06-05 22:53:05 +00:00
void Operator : : CalcPEC_Curves ( )
{
2010-03-22 07:19:17 +00:00
//special treatment for primitives of type curve (treated as wires)
double p1 [ 3 ] ;
double p2 [ 3 ] ;
struct Grid_Path path ;
vector < CSProperties * > vec_prop = CSX - > GetPropertyByType ( CSProperties : : METAL ) ;
for ( size_t p = 0 ; p < vec_prop . size ( ) ; + + p )
{
CSProperties * prop = vec_prop . at ( p ) ;
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 ) ;
curv - > GetPoint ( i , p2 ) ;
path = FindPath ( p1 , p2 ) ;
2010-05-29 15:40:18 +00:00
if ( path . dir . size ( ) > 0 )
prim - > SetPrimitiveUsed ( true ) ;
2010-03-22 07:19:17 +00:00
for ( size_t t = 0 ; t < path . dir . size ( ) ; + + t )
{
2010-03-26 07:23:37 +00:00
// cerr << path.dir.at(t) << " " << path.posPath[0].at(t) << " " << path.posPath[1].at(t) << " " << path.posPath[2].at(t) << endl;
2010-10-20 05:26:16 +00:00
SetVV ( path . dir . at ( t ) , path . posPath [ 0 ] . at ( t ) , path . posPath [ 1 ] . at ( t ) , path . posPath [ 2 ] . at ( t ) , 0 ) ;
SetVI ( path . dir . at ( t ) , path . posPath [ 0 ] . at ( t ) , path . posPath [ 1 ] . at ( t ) , path . posPath [ 2 ] . at ( t ) , 0 ) ;
2010-05-10 07:14:29 +00:00
+ + m_Nr_PEC [ path . dir . at ( t ) ] ;
2010-05-29 15:40:18 +00:00
}
2010-03-26 07:23:37 +00:00
// cerr << "found path size: " << path.dir.size() << endl;
2010-03-22 07:19:17 +00:00
}
}
}
}
2010-03-05 13:20:25 +00:00
}
2010-04-25 19:59:05 +00:00
void Operator : : AddExtension ( Operator_Extension * op_ext )
{
m_Op_exts . push_back ( op_ext ) ;
}