zcy
/
openEMS
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
openEMS/FDTD/operator.cpp

2132 lines
56 KiB
C++
Raw Blame History

/*
* 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 <fstream>
#include <algorithm>
#include "operator.h"
#include "engine.h"
#include "extensions/operator_extension.h"
#include "extensions/operator_ext_excitation.h"
#include "Common/processfields.h"
#include "tools/array_ops.h"
#include "tools/vtk_file_writer.h"
#include "fparser.hh"
#include "extensions/operator_ext_excitation.h"
#include "vtkPolyData.h"
#include "vtkCellArray.h"
#include "vtkPoints.h"
#include "vtkXMLPolyDataWriter.h"
#include "CSPrimBox.h"
#include "CSPrimCurve.h"
#include "CSPropMaterial.h"
#include "CSPropLumpedElement.h"
Operator* Operator::New()
{
cout << "Create FDTD operator" << endl;
Operator* op = new Operator();
op->Init();
return op;
}
Operator::Operator() : Operator_Base()
{
m_Exc = 0;
m_InvaildTimestep = false;
m_TimeStepVar = 3;
}
Operator::~Operator()
{
for (size_t n=0; n<m_Op_exts.size(); ++n)
delete m_Op_exts.at(n);
m_Op_exts.clear();
Delete();
}
Engine* Operator::CreateEngine()
{
m_Engine = Engine::New(this);
return m_Engine;
}
void Operator::Init()
{
CSX = NULL;
m_Engine = NULL;
Operator_Base::Init();
vv=NULL;
vi=NULL;
iv=NULL;
ii=NULL;
m_epsR=NULL;
m_kappa=NULL;
m_mueR=NULL;
m_sigma=NULL;
MainOp=NULL;
for (int n=0; n<3; ++n)
{
EC_C[n]=NULL;
EC_G[n]=NULL;
EC_L[n]=NULL;
EC_R[n]=NULL;
}
m_Exc = 0;
m_TimeStepFactor = 1;
SetMaterialAvgMethod(QuarterCell);
}
void Operator::Delete()
{
CSX = NULL;
Delete_N_3DArray(vv,numLines);
Delete_N_3DArray(vi,numLines);
Delete_N_3DArray(iv,numLines);
Delete_N_3DArray(ii,numLines);
vv=vi=iv=ii=0;
delete MainOp; MainOp=0;
for (int n=0; n<3; ++n)
{
delete[] EC_C[n];EC_C[n]=0;
delete[] EC_G[n];EC_G[n]=0;
delete[] EC_L[n];EC_L[n]=0;
delete[] EC_R[n];EC_R[n]=0;
}
Delete_N_3DArray(m_epsR,numLines);
m_epsR=0;
Delete_N_3DArray(m_kappa,numLines);
m_kappa=0;
Delete_N_3DArray(m_mueR,numLines);
m_mueR=0;
Delete_N_3DArray(m_sigma,numLines);
m_sigma=0;
}
void Operator::Reset()
{
Delete();
Operator_Base::Reset();
}
double Operator::GetDiscLine(int n, unsigned int pos, bool dualMesh) const
{
if ((n<0) || (n>2)) return 0.0;
if (pos>=numLines[n]) return 0.0;
if (dualMesh==false)
return discLines[n][pos];
// return dual mesh node
if (pos<numLines[n]-1)
return 0.5*(discLines[n][pos] + discLines[n][pos+1]);
// dual node for the last line (outside the field domain)
return discLines[n][pos] + 0.5*(discLines[n][pos] - discLines[n][pos-1]);
}
double Operator::GetDiscDelta(int n, unsigned int pos, bool dualMesh) const
{
if ((n<0) || (n>2)) return 0.0;
if (pos>=numLines[n]) return 0.0;
double delta=0;
if (dualMesh==false)
{
if (pos<numLines[n]-1)
delta = GetDiscLine(n,pos+1,false) - GetDiscLine(n,pos,false);
else
delta = GetDiscLine(n,pos,false) - GetDiscLine(n,pos-1,false);
return delta;
}
else
{
if (pos>0)
delta = GetDiscLine(n,pos,true) - GetDiscLine(n,pos-1,true);
else
delta = GetDiscLine(n,1,false) - GetDiscLine(n,0,false);
return delta;
}
}
bool Operator::GetYeeCoords(int ny, unsigned int pos[3], double* coords, bool dualMesh) const
{
for (int n=0;n<3;++n)
coords[n]=GetDiscLine(n,pos[n],dualMesh);
coords[ny]=GetDiscLine(ny,pos[ny],!dualMesh);
//check if position is inside the FDTD domain
if (dualMesh==false) //main grid
{
if (pos[ny]>=numLines[ny]-1)
return false;
}
else //dual grid
{
int nP = (ny+1)%3;
int nPP = (ny+2)%3;
if ((pos[nP]>=numLines[nP]-1) || (pos[nPP]>=numLines[nPP]-1))
return false;
}
return true;
}
bool Operator::GetNodeCoords(const unsigned int pos[3], double* coords, bool dualMesh, CoordinateSystem c_system) const
{
for (int n=0;n<3;++n)
coords[n]=GetDiscLine(n,pos[n],dualMesh);
TransformCoordSystem(coords,coords,m_MeshType,c_system);
return true;
}
double Operator::GetEdgeLength(int n, const unsigned int* pos, bool dualMesh) const
{
return GetDiscDelta(n,pos[n],dualMesh)*gridDelta;
}
double Operator::GetCellVolume(const unsigned int pos[3], bool dualMesh) const
{
double vol=1;
for (int n=0;n<3;++n)
vol*=GetEdgeLength(n,pos,dualMesh);
return vol;
}
double Operator::GetNodeWidth(int ny, const int pos[3], bool dualMesh) const
{
if ( (pos[0]<0) || (pos[1]<0) || (pos[2]<0) )
return 0.0;
//call the unsigned int version of GetNodeWidth
unsigned int uiPos[]={(unsigned int)pos[0],(unsigned int)pos[1],(unsigned int)pos[2]};
return GetNodeWidth(ny, uiPos, dualMesh);
}
double Operator::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) const
{
int nyP = (ny+1)%3;
int nyPP = (ny+2)%3;
return GetNodeWidth(nyP,pos,dualMesh) * GetNodeWidth(nyPP,pos,dualMesh);
}
double Operator::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
{
if ( (pos[0]<0) || (pos[1]<0) || (pos[2]<0) )
return 0.0;
//call the unsigned int version of GetNodeArea
unsigned int uiPos[]={(unsigned int)pos[0],(unsigned int)pos[1],(unsigned int)pos[2]};
return GetNodeArea(ny, uiPos, dualMesh);
}
unsigned int Operator::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh, bool fullMesh) const
{
inside = false;
if ((ny<0) || (ny>2))
return 0;
if (coord<GetDiscLine(ny,0))
return 0;
unsigned int numLines = GetNumberOfLines(ny, fullMesh);
if (coord>GetDiscLine(ny,numLines-1))
return numLines-1;
inside=true;
if (dualMesh==false)
{
for (unsigned int n=0;n<numLines;++n)
{
if (coord<=GetDiscLine(ny,n,true))
return n;
}
}
else
{
for (unsigned int n=1;n<numLines;++n)
{
if (coord<=GetDiscLine(ny,n,false))
return n-1;
}
}
//should not happen
return 0;
}
bool Operator::SnapToMesh(const double* dcoord, unsigned int* uicoord, bool dualMesh, bool fullMesh, bool* inside) const
{
bool meshInside=false;
bool ok=true;
for (int n=0; n<3; ++n)
{
uicoord[n] = SnapToMeshLine(n,dcoord[n],meshInside,dualMesh,fullMesh);
ok &= meshInside;
if (inside)
inside[n]=meshInside;
}
return ok;
}
int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, bool fullMesh, int SnapMethod, bool* bStartIn, bool* bStopIn) const
{
double l_start[3], l_stop[3];
for (int n=0;n<3;++n)
{
l_start[n] = fmin(start[n],stop[n]);
l_stop[n] = fmax(start[n], stop[n]);
double min = GetDiscLine(n,0);
double max = GetDiscLine(n,GetNumberOfLines(n, fullMesh)-1);
if ( ((l_start[n]<min) && (l_stop[n]<min)) || ((l_start[n]>max) && (l_stop[n]>max)) )
{
return -2;
}
}
SnapToMesh(l_start, uiStart, dualMesh, fullMesh, bStartIn);
SnapToMesh(l_stop, uiStop, dualMesh, fullMesh, bStopIn);
int iDim = 0;
if (SnapMethod==0)
{
for (int n=0;n<3;++n)
if (uiStop[n]>uiStart[n])
++iDim;
return iDim;
}
else if (SnapMethod==1)
{
for (int n=0;n<3;++n)
{
if (uiStop[n]>uiStart[n])
{
if ((GetDiscLine( n, uiStart[n], dualMesh ) > l_start[n]) && (uiStart[n]>0))
--uiStart[n];
if ((GetDiscLine( n, uiStop[n], dualMesh ) < l_stop[n]) && (uiStop[n]<GetNumberOfLines(n, fullMesh)-1))
++uiStop[n];
}
if (uiStop[n]>uiStart[n])
++iDim;
}
return iDim;
}
else if (SnapMethod==2)
{
for (int n=0;n<3;++n)
{
if (uiStop[n]>uiStart[n])
{
if ((GetDiscLine( n, uiStart[n], dualMesh ) < l_start[n]) && (uiStart[n]<GetNumberOfLines(n, fullMesh)-1))
++uiStart[n];
if ((GetDiscLine( n, uiStop[n], dualMesh ) > l_stop[n]) && (uiStop[n]>0))
--uiStop[n];
}
if (uiStop[n]>uiStart[n])
++iDim;
}
return iDim;
}
else
cerr << "Operator::SnapBox2Mesh: Unknown snapping method!" << endl;
return -1;
}
int Operator::SnapLine2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, bool fullMesh) const
{
bool bStartIn[3];
bool bStopIn[3];
SnapToMesh(start, uiStart, dualMesh, fullMesh, bStartIn);
SnapToMesh(stop, uiStop, dualMesh, fullMesh, bStopIn);
for (int n=0;n<3;++n)
{
if ((start[n]<GetDiscLine(n,0)) && (stop[n]<GetDiscLine(n,0)))
return -1; //lower bound violation
if ((start[n]>GetDiscLine(n,GetNumberOfLines(n,true)-1)) && (stop[n]>GetDiscLine(n,GetNumberOfLines(n,true)-1)))
return -1; //upper bound violation
}
int ret = 0;
if (!(bStartIn[0] && bStartIn[1] && bStartIn[2]))
ret = ret + 1;
if (!(bStopIn[0] && bStopIn[1] && bStopIn[2]))
ret = ret + 2;
if (ret==0)
return ret;
//fixme, do we need to do something about start or stop being outside the field domain?
//maybe caclulate the intersection point and snap to that?
//it seems to work like this as well...
return ret;
}
Grid_Path Operator::FindPath(double start[], double stop[])
{
Grid_Path path;
unsigned int uiStart[3],uiStop[3],currPos[3];
int ret = SnapLine2Mesh(start, stop, uiStart, uiStop, false, true);
if (ret<0)
return path;
currPos[0]=uiStart[0];
currPos[1]=uiStart[1];
currPos[2]=uiStart[2];
double meshStart[3] = {discLines[0][uiStart[0]], discLines[1][uiStart[1]], discLines[2][uiStart[2]]};
double meshStop[3] = {discLines[0][uiStop[0]], discLines[1][uiStop[1]], discLines[2][uiStop[2]]};
bool UpDir = false;
double foot=0,dist=0,minFoot=0,minDist=0;
int minDir=0;
unsigned int minPos[3];
double startFoot,stopFoot,currFoot;
Point_Line_Distance(meshStart,start,stop,startFoot,dist, m_MeshType);
Point_Line_Distance(meshStop,start,stop,stopFoot,dist, m_MeshType);
currFoot=startFoot;
minFoot=startFoot;
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]];
if (((int)currPos[n]-1)>=0)
{
P[n] = discLines[n][currPos[n]-1];
Point_Line_Distance(P,start,stop,foot,dist, m_MeshType);
if ((foot>currFoot) && (dist<minDist))
{
minFoot=foot;
minDist=dist;
minDir = n;
UpDir = false;
}
}
if ((currPos[n]+1)<numLines[n])
{
P[n] = discLines[n][currPos[n]+1];
Point_Line_Distance(P,start,stop,foot,dist, m_MeshType);
if ((foot>currFoot) && (dist<minDist))
{
minFoot=foot;
minDist=dist;
minDir = n;
UpDir = true;
}
}
}
minPos[0]=currPos[0];
minPos[1]=currPos[1];
minPos[2]=currPos[2];
if (UpDir)
{
currPos[minDir]+=1;
}
else
{
currPos[minDir]+=-1;
minPos[minDir]-=1;
}
//check validity of current postion
for (int n=0;n<3;++n)
if (currPos[n]>=numLines[n])
{
cerr << __func__ << ": Error, path went out of simulation domain, skipping path!" << endl;
Grid_Path empty;
return empty;
}
path.posPath[0].push_back(minPos[0]);
path.posPath[1].push_back(minPos[1]);
path.posPath[2].push_back(minPos[2]);
currFoot=minFoot;
path.dir.push_back(minDir);
}
//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);
}
}
return path;
}
void Operator::SetMaterialAvgMethod(MatAverageMethods method)
{
switch (method)
{
default:
case QuarterCell:
return SetQuarterCellMaterialAvg();
case CentralCell:
return SetCellConstantMaterial();
}
}
double Operator::GetNumberCells() const
{
if (numLines)
return (numLines[0])*(numLines[1])*(numLines[2]); //it's more like number of nodes???
return 0;
}
void Operator::ShowStat() const
{
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;
cout << "------- Stat: FDTD Operator -------" << endl;
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;
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;
cout << "-----------------------------------" << endl;
cout << "Background materials (epsR/mueR/kappa/sigma): " << GetBackgroundEpsR() << "/" << GetBackgroundMueR() << "/" << GetBackgroundKappa() << "/" << GetBackgroundSigma() << endl;
cout << "-----------------------------------" << endl;
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;
cout << "Timestep (s)\t\t: " << dT ;
if (opt_dT)
cout <<"\t(" << opt_dT << ")";
cout << endl;
cout << "Timestep method name\t: " << m_Used_TS_Name << endl;
cout << "Nyquist criteria (TS)\t: " << m_Exc->GetNyquistNum() << endl;
cout << "Nyquist criteria (s)\t: " << m_Exc->GetNyquistNum()*dT << endl;
cout << "-----------------------------------" << endl;
}
void Operator::ShowExtStat() const
{
if (m_Op_exts.size()==0) return;
cout << "-----------------------------------" << endl;
for (size_t n=0; n<m_Op_exts.size(); ++n)
m_Op_exts.at(n)->ShowStat(cout);
cout << "-----------------------------------" << endl;
}
void Operator::DumpOperator2File(string filename)
{
#ifdef OUTPUT_IN_DRAWINGUNITS
double discLines_scaling = 1;
#else
double discLines_scaling = GetGridDelta();
#endif
cout << "Operator: Dumping FDTD operator information to vtk file: " << filename << " ..." << flush;
VTK_File_Writer* vtk_Writer = new VTK_File_Writer(filename.c_str(), m_MeshType);
vtk_Writer->SetMeshLines(discLines,numLines,discLines_scaling);
vtk_Writer->SetHeader("openEMS - Operator dump");
vtk_Writer->SetNativeDump(true);
//find excitation extension
Operator_Ext_Excitation* Op_Ext_Exc=GetExcitationExtension();
if (Op_Ext_Exc)
{
FDTD_FLOAT**** exc = NULL;
if (Op_Ext_Exc->Volt_Count>0)
{
exc = Create_N_3DArray<FDTD_FLOAT>(numLines);
for (unsigned int n=0; n< Op_Ext_Exc->Volt_Count; ++n)
exc[ Op_Ext_Exc->Volt_dir[n]][ Op_Ext_Exc->Volt_index[0][n]][ Op_Ext_Exc->Volt_index[1][n]][ Op_Ext_Exc->Volt_index[2][n]] = Op_Ext_Exc->Volt_amp[n];
vtk_Writer->AddVectorField("exc_volt",exc);
Delete_N_3DArray(exc,numLines);
}
if ( Op_Ext_Exc->Curr_Count>0)
{
exc = Create_N_3DArray<FDTD_FLOAT>(numLines);
for (unsigned int n=0; n< Op_Ext_Exc->Curr_Count; ++n)
exc[ Op_Ext_Exc->Curr_dir[n]][ Op_Ext_Exc->Curr_index[0][n]][ Op_Ext_Exc->Curr_index[1][n]][ Op_Ext_Exc->Curr_index[2][n]] = Op_Ext_Exc->Curr_amp[n];
vtk_Writer->AddVectorField("exc_curr",exc);
Delete_N_3DArray(exc,numLines);
}
}
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);
}
vtk_Writer->AddVectorField("vv",vv_temp);
Delete_N_3DArray(vv_temp,numLines);
vtk_Writer->AddVectorField("vi",vi_temp);
Delete_N_3DArray(vi_temp,numLines);
vtk_Writer->AddVectorField("iv",iv_temp);
Delete_N_3DArray(iv_temp,numLines);
vtk_Writer->AddVectorField("ii",ii_temp);
Delete_N_3DArray(ii_temp,numLines);
if (vtk_Writer->Write()==false)
cerr << "Operator::DumpOperator2File: Error: Can't write file... skipping!" << endl;
delete vtk_Writer;
}
//! \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 , unsigned int *range)
{
cout << "Operator: Dumping PEC information to vtk file: " << filename << " ..." << flush;
#ifdef OUTPUT_IN_DRAWINGUNITS
double scaling = 1.0;
#else
double scaling = GetGridDelta();;
#endif
unsigned int start[3] = {0, 0, 0};
unsigned int stop[3] = {numLines[0]-1,numLines[1]-1,numLines[2]-1};
if (range!=NULL)
for (int n=0;n<3;++n)
{
start[n] = range[2*n];
stop[n] = range[2*n+1];
}
vtkPolyData* polydata = vtkPolyData::New();
vtkCellArray *poly = vtkCellArray::New();
vtkPoints *points = vtkPoints::New();
int* pointIdx[2];
pointIdx[0] = new int[numLines[0]*numLines[1]];
pointIdx[1] = new int[numLines[0]*numLines[1]];
// init point idx
for (unsigned int n=0;n<numLines[0]*numLines[1];++n)
{
pointIdx[0][n]=-1;
pointIdx[1][n]=-1;
}
int nP,nPP;
double coord[3];
unsigned int pos[3],rpos[3];
unsigned int mesh_idx=0;
for (pos[2]=start[2];pos[2]<stop[2];++pos[2])
{ // each xy-plane
for (unsigned int n=0;n<numLines[0]*numLines[1];++n)
{
pointIdx[0][n]=pointIdx[1][n];
pointIdx[1][n]=-1;
}
for (pos[0]=start[0];pos[0]<stop[0];++pos[0])
for (pos[1]=start[1];pos[1]<stop[1];++pos[1])
{
for (int n=0;n<3;++n)
{
nP = (n+1)%3;
nPP = (n+2)%3;
if ((GetVV(n,pos) == 0) && (GetVI(n,pos) == 0) && (pos[nP]>0) && (pos[nPP]>0))
{
rpos[0]=pos[0];
rpos[1]=pos[1];
rpos[2]=pos[2];
poly->InsertNextCell(2);
mesh_idx = rpos[0] + rpos[1]*numLines[0];
if (pointIdx[0][mesh_idx]<0)
{
for (int m=0;m<3;++m)
coord[m] = discLines[m][rpos[m]];
TransformCoordSystem(coord, coord, m_MeshType, CARTESIAN);
for (int m=0;m<3;++m)
coord[m] *= scaling;
pointIdx[0][mesh_idx] = (int)points->InsertNextPoint(coord);
}
poly->InsertCellPoint(pointIdx[0][mesh_idx]);
++rpos[n];
mesh_idx = rpos[0] + rpos[1]*numLines[0];
if (pointIdx[n==2][mesh_idx]<0)
{
for (int m=0;m<3;++m)
coord[m] = discLines[m][rpos[m]];
TransformCoordSystem(coord, coord, m_MeshType, CARTESIAN);
for (int m=0;m<3;++m)
coord[m] *= scaling;
pointIdx[n==2][mesh_idx] = (int)points->InsertNextPoint(coord);
}
poly->InsertCellPoint(pointIdx[n==2][mesh_idx]);
}
}
}
}
delete[] pointIdx[0];
delete[] pointIdx[1];
polydata->SetPoints(points);
points->Delete();
polydata->SetLines(poly);
poly->Delete();
vtkXMLPolyDataWriter* writer = vtkXMLPolyDataWriter::New();
filename += ".vtp";
writer->SetFileName(filename.c_str());
#if VTK_MAJOR_VERSION==6
writer->SetInputData(polydata);
#else
writer->SetInput(polydata);
#endif
writer->Write();
writer->Delete();
polydata->Delete();
cout << " done." << endl;
}
void Operator::DumpMaterial2File(string filename)
{
#ifdef OUTPUT_IN_DRAWINGUNITS
double discLines_scaling = 1;
#else
double discLines_scaling = GetGridDelta();
#endif
cout << "Operator: Dumping material information to vtk file: " << filename << " ..." << flush;
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);
unsigned int pos[3];
for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{
vector<CSPrimitives*> vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL);
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{
for (int n=0; n<3; ++n)
{
double inMat[4];
Calc_EffMatPos(n, pos, inMat, vPrims);
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];
}
}
}
}
VTK_File_Writer* vtk_Writer = new VTK_File_Writer(filename.c_str(), m_MeshType);
vtk_Writer->SetMeshLines(discLines,numLines,discLines_scaling);
vtk_Writer->SetHeader("openEMS - material dump");
vtk_Writer->SetNativeDump(true);
vtk_Writer->AddVectorField("epsilon",epsilon);
Delete_N_3DArray(epsilon,numLines);
vtk_Writer->AddVectorField("mue",mue);
Delete_N_3DArray(mue,numLines);
vtk_Writer->AddVectorField("kappa",kappa);
Delete_N_3DArray(kappa,numLines);
vtk_Writer->AddVectorField("sigma",sigma);
Delete_N_3DArray(sigma,numLines);
if (vtk_Writer->Write()==false)
cerr << "Operator::DumpMaterial2File: Error: Can't write file... skipping!" << endl;
delete vtk_Writer;
}
bool Operator::SetupCSXGrid(CSRectGrid* grid)
{
for (int n=0; n<3; ++n)
{
discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true);
if (numLines[n]<3)
{
cerr << "CartOperator::SetupCSXGrid: you need at least 3 disc-lines in every direction (3D!)!!!" << endl;
Reset();
return false;
}
}
MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]);
MainOp->SetGrid(discLines[0],discLines[1],discLines[2]);
if (grid->GetDeltaUnit()<=0)
{
cerr << "CartOperator::SetupCSXGrid: grid delta unit must not be <=0 !!!" << endl;
Reset();
return false;
}
else gridDelta=grid->GetDeltaUnit();
MainOp->SetGridDelta(1);
MainOp->AddCellAdrOp();
//delete the grid clone...
delete grid;
return true;
}
bool Operator::SetGeometryCSX(ContinuousStructure* geo)
{
if (geo==NULL) return false;
CSX = geo;
CSBackgroundMaterial* bg_mat=CSX->GetBackgroundMaterial();
SetBackgroundEpsR(bg_mat->GetEpsilon());
SetBackgroundMueR(bg_mat->GetMue());
SetBackgroundKappa(bg_mat->GetKappa());
SetBackgroundSigma(bg_mat->GetSigma());
SetBackgroundDensity(0);
CSRectGrid* grid=CSX->GetGrid();
return SetupCSXGrid(CSRectGrid::Clone(grid));
}
void Operator::InitOperator()
{
Delete_N_3DArray(vv,numLines);
Delete_N_3DArray(vi,numLines);
Delete_N_3DArray(iv,numLines);
Delete_N_3DArray(ii,numLines);
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);
}
void Operator::InitDataStorage()
{
if (m_StoreMaterial[0])
{
if (g_settings.GetVerboseLevel()>0)
cerr << "Operator::InitDataStorage(): Storing epsR material data..." << endl;
Delete_N_3DArray(m_epsR,numLines);
m_epsR = Create_N_3DArray<float>(numLines);
}
if (m_StoreMaterial[1])
{
if (g_settings.GetVerboseLevel()>0)
cerr << "Operator::InitDataStorage(): Storing kappa material data..." << endl;
Delete_N_3DArray(m_kappa,numLines);
m_kappa = Create_N_3DArray<float>(numLines);
}
if (m_StoreMaterial[2])
{
if (g_settings.GetVerboseLevel()>0)
cerr << "Operator::InitDataStorage(): Storing muR material data..." << endl;
Delete_N_3DArray(m_mueR,numLines);
m_mueR = Create_N_3DArray<float>(numLines);
}
if (m_StoreMaterial[3])
{
if (g_settings.GetVerboseLevel()>0)
cerr << "Operator::InitDataStorage(): Storing sigma material data..." << endl;
Delete_N_3DArray(m_sigma,numLines);
m_sigma = Create_N_3DArray<float>(numLines);
}
}
void Operator::CleanupMaterialStorage()
{
if (!m_StoreMaterial[0] && m_epsR)
{
if (g_settings.GetVerboseLevel()>0)
cerr << "Operator::CleanupMaterialStorage(): Delete epsR material data..." << endl;
Delete_N_3DArray(m_epsR,numLines);
m_epsR = NULL;
}
if (!m_StoreMaterial[1] && m_kappa)
{
if (g_settings.GetVerboseLevel()>0)
cerr << "Operator::CleanupMaterialStorage(): Delete kappa material data..." << endl;
Delete_N_3DArray(m_kappa,numLines);
m_kappa = NULL;
}
if (!m_StoreMaterial[2] && m_mueR)
{
if (g_settings.GetVerboseLevel()>0)
cerr << "Operator::CleanupMaterialStorage(): Delete mueR material data..." << endl;
Delete_N_3DArray(m_mueR,numLines);
m_mueR = NULL;
}
if (!m_StoreMaterial[3] && m_sigma)
{
if (g_settings.GetVerboseLevel()>0)
cerr << "Operator::CleanupMaterialStorage(): Delete sigma material data..." << endl;
Delete_N_3DArray(m_sigma,numLines);
m_sigma = NULL;
}
}
double Operator::GetDiscMaterial(int type, int n, const unsigned int pos[3]) const
{
switch (type)
{
case 0:
if (m_epsR==0)
return 0;
return m_epsR[n][pos[0]][pos[1]][pos[2]];
case 1:
if (m_kappa==0)
return 0;
return m_kappa[n][pos[0]][pos[1]][pos[2]];
case 2:
if (m_mueR==0)
return 0;
return m_mueR[n][pos[0]][pos[1]][pos[2]];
case 3:
if (m_sigma==0)
return 0;
return m_sigma[n][pos[0]][pos[1]][pos[2]];
}
return 0;
}
void Operator::SetExcitationSignal(Excitation* exc)
{
m_Exc=exc;
}
void Operator::Calc_ECOperatorPos(int n, unsigned int* pos)
{
unsigned int i = MainOp->SetPos(pos[0],pos[1],pos[2]);
double C = EC_C[n][i];
double G = EC_G[n][i];
if (C>0)
{
SetVV(n,pos[0],pos[1],pos[2], (1.0-dT*G/2.0/C)/(1.0+dT*G/2.0/C) );
SetVI(n,pos[0],pos[1],pos[2], (dT/C)/(1.0+dT*G/2.0/C) );
}
else
{
SetVV(n,pos[0],pos[1],pos[2], 0 );
SetVI(n,pos[0],pos[1],pos[2], 0 );
}
double L = EC_L[n][i];
double R = EC_R[n][i];
if (L>0)
{
SetII(n,pos[0],pos[1],pos[2], (1.0-dT*R/2.0/L)/(1.0+dT*R/2.0/L) );
SetIV(n,pos[0],pos[1],pos[2], (dT/L)/(1.0+dT*R/2.0/L) );
}
else
{
SetII(n,pos[0],pos[1],pos[2], 0 );
SetIV(n,pos[0],pos[1],pos[2], 0 );
}
}
int Operator::CalcECOperator( DebugFlags debugFlags )
{
Init_EC();
InitDataStorage();
if (Calc_EC()==0)
return -1;
m_InvaildTimestep = false;
opt_dT = 0;
if (dT>0)
{
double save_dT = dT;
CalcTimestep();
opt_dT = dT;
if (dT<save_dT)
{
cerr << "Operator::CalcECOperator: Warning, forced timestep: " << save_dT << "s is larger than calculated timestep: " << dT << "s! It is not recommended using this timestep!! " << endl;
m_InvaildTimestep = true;
}
dT = save_dT;
}
else
CalcTimestep();
dT*=m_TimeStepFactor;
if (m_Exc->GetSignalPeriod()>0)
{
unsigned int TS = ceil(m_Exc->GetSignalPeriod()/dT);
double new_dT = m_Exc->GetSignalPeriod()/TS;
cout << "Operartor::CalcECOperator: Decreasing timestep by " << round((dT-new_dT)/dT*1000)/10.0 << "% to " << new_dT << " (" << dT << ") to match periodic signal" << endl;
dT = new_dT;
}
m_Exc->Reset(dT);
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])
{
Calc_ECOperatorPos(n,pos);
}
}
}
}
//Apply PEC to all boundary's
bool PEC[6]={1,1,1,1,1,1};
//make an exception for BC == -1
for (int n=0; n<6; ++n)
if ((m_BC[n]==-1))
PEC[n] = false;
ApplyElectricBC(PEC);
CalcPEC();
Calc_LumpedElements();
bool PMC[6];
for (int n=0; n<6; ++n)
PMC[n] = m_BC[n]==1;
ApplyMagneticBC(PMC);
//all information available for extension... create now...
for (size_t n=0; n<m_Op_exts.size(); ++n)
m_Op_exts.at(n)->BuildExtension();
//remove inactive extensions
vector<Operator_Extension*>::iterator it = m_Op_exts.begin();
while (it!=m_Op_exts.end())
{
if ( (*it)->IsActive() == false)
{
DeleteExtension((*it));
it = m_Op_exts.begin(); //restart search for inactive extension
}
else
++it;
}
if (debugFlags & debugMaterial)
DumpMaterial2File( "material_dump" );
if (debugFlags & debugOperator)
DumpOperator2File( "operator_dump" );
if (debugFlags & debugPEC)
DumpPEC2File( "PEC_dump" );
//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;
}
return 0;
}
void Operator::ApplyElectricBC(bool* dirs)
{
if (!dirs)
return;
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])
{
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 );
}
}
}
}
}
void Operator::ApplyMagneticBC(bool* dirs)
{
if (!dirs)
return;
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])
{
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 );
}
// the last current lines are outside the FDTD domain and cannot be iterated by the FDTD engine
pos[n] = numLines[n]-1;
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 );
}
}
}
}
bool Operator::Calc_ECPos(int ny, const unsigned int* pos, double* EC, vector<CSPrimitives*> vPrims) const
{
double EffMat[4];
Calc_EffMatPos(ny,pos,EffMat, vPrims);
if (m_epsR)
m_epsR[ny][pos[0]][pos[1]][pos[2]] = EffMat[0];
if (m_kappa)
m_kappa[ny][pos[0]][pos[1]][pos[2]] = EffMat[1];
if (m_mueR)
m_mueR[ny][pos[0]][pos[1]][pos[2]] = EffMat[2];
if (m_sigma)
m_sigma[ny][pos[0]][pos[1]][pos[2]] = EffMat[3];
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;
}
return true;
}
double Operator::GetRawDiscDelta(int ny, const int pos) const
{
//numLines[ny] is expected to be larger then 1 !
if (pos<0)
return (discLines[ny][0] - discLines[ny][1]);
if (pos>=(int)numLines[ny]-1)
return (discLines[ny][numLines[ny]-2] - discLines[ny][numLines[ny]-1]);
return (discLines[ny][pos+1] - discLines[ny][pos]);
}
bool Operator::GetCellCenterMaterialAvgCoord(const int pos[], double coord[3]) const
{
unsigned int ui_pos[3];
for (int n=0;n<3;++n)
{
if ((pos[n]<0) || (pos[n]>=(int)numLines[n]))
return false;
ui_pos[n] = pos[n];
}
GetNodeCoords(ui_pos, coord, true);
return true;
}
double Operator::GetMaterial(int ny, const double* coords, int MatType, vector<CSPrimitives*> vPrims, bool markAsUsed) const
{
CSProperties* prop = CSX->GetPropertyByCoordPriority(coords,vPrims,markAsUsed);
// CSProperties* old_prop = CSX->GetPropertyByCoordPriority(coords,CSProperties::MATERIAL,markAsUsed);
// if (old_prop!=prop)
// {
// cerr << "ERROR: Unequal properties!" << endl;
// exit(-1);
// }
CSPropMaterial* mat = dynamic_cast<CSPropMaterial*>(prop);
if (mat)
{
switch (MatType)
{
case 0:
return mat->GetEpsilonWeighted(ny,coords);
case 1:
return mat->GetKappaWeighted(ny,coords);
case 2:
return mat->GetMueWeighted(ny,coords);
case 3:
return mat->GetSigmaWeighted(ny,coords);
case 4:
return mat->GetDensityWeighted(coords);
default:
cerr << "Operator::GetMaterial: Error: unknown material type" << endl;
return 0;
}
}
switch (MatType)
{
case 0:
return GetBackgroundEpsR();
case 1:
return GetBackgroundKappa();
case 2:
return GetBackgroundMueR();
case 3:
return GetBackgroundSigma();
case 4:
return GetBackgroundDensity();
default:
cerr << "Operator::GetMaterial: Error: unknown material type" << endl;
return 0;
}
}
bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives *> vPrims) const
{
int n=ny;
double coord[3];
int nP = (n+1)%3;
int nPP = (n+2)%3;
int loc_pos[3] = {(int)pos[0],(int)pos[1],(int)pos[2]};
double A_n;
double area = 0;
EffMat[0] = 0;
EffMat[1] = 0;
EffMat[2] = 0;
EffMat[3] = 0;
//******************************* epsilon,kappa averaging *****************************//
//shift up-right
if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{
A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
area+=A_n;
}
//shift up-left
--loc_pos[nP];
if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{
A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
area+=A_n;
}
//shift down-right
++loc_pos[nP];
--loc_pos[nPP];
if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{
A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
area+=A_n;
}
//shift down-left
--loc_pos[nP];
if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{
A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
area+=A_n;
}
EffMat[0]*=__EPS0__/area;
EffMat[1]/=area;
//******************************* mu,sigma averaging *****************************//
loc_pos[0]=pos[0];
loc_pos[1]=pos[1];
loc_pos[2]=pos[2];
double length=0;
double delta_ny,sigma;
//shift down
--loc_pos[n];
if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{
delta_ny = GetNodeWidth(n,loc_pos,true);
EffMat[2] += delta_ny / GetMaterial(n, coord, 2, vPrims);
sigma = GetMaterial(n, coord, 3, vPrims);
if (sigma)
EffMat[3] += delta_ny / sigma;
else
EffMat[3] = 0;
length+=delta_ny;
}
//shift up
++loc_pos[n];
if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{
delta_ny = GetNodeWidth(n,loc_pos,true);
EffMat[2] += delta_ny / GetMaterial(n, coord, 2, vPrims);
sigma = GetMaterial(n, coord, 3, vPrims);
if (sigma)
EffMat[3] += delta_ny / sigma;
else
EffMat[3] = 0;
length+=delta_ny;
}
EffMat[2] = length * __MUE0__ / EffMat[2];
if (EffMat[3]) EffMat[3]=length / EffMat[3];
for (int n=0; n<4; ++n)
if (isnan(EffMat[n]) || isinf(EffMat[n]))
{
cerr << "Operator::" << __func__ << ": Error, an effective material parameter is not a valid result, this should NOT have happend... exit..." << endl;
cerr << ny << "@" << n << " : " << pos[0] << "," << pos[1] << "," << pos[2] << endl;
exit(0);
}
return true;
}
bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives*> vPrims) const
{
int n=ny;
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=GetRawDiscDelta(n,pos[n]);
double deltaP=GetRawDiscDelta(nP,pos[nP]);
double deltaPP=GetRawDiscDelta(nPP,pos[nPP]);
double delta_M=GetRawDiscDelta(n,pos[n]-1);
double deltaP_M=GetRawDiscDelta(nP,pos[nP]-1);
double deltaPP_M=GetRawDiscDelta(nPP,pos[nPP]-1);
int loc_pos[3] = {(int)pos[0],(int)pos[1],(int)pos[2]};
double A_n;
double area = 0;
//******************************* 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;
A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] = GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
EffMat[1] = GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
area+=A_n;
//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;
--loc_pos[nP];
A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
area+=A_n;
//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;
++loc_pos[nP];
--loc_pos[nPP];
A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
area+=A_n;
//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;
--loc_pos[nP];
A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
area+=A_n;
EffMat[0]*=__EPS0__/area;
EffMat[1]/=area;
//******************************* mu,sigma averaging *****************************//
loc_pos[0]=pos[0];
loc_pos[1]=pos[1];
loc_pos[2]=pos[2];
double length=0;
//shift down
shiftCoord[n] = coord[n]-delta_M*0.25;
shiftCoord[nP] = coord[nP]+deltaP*0.5;
shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
--loc_pos[n];
double delta_ny = GetNodeWidth(n,loc_pos,true);
EffMat[2] = delta_ny / GetMaterial(n, shiftCoord, 2, vPrims);
double sigma = GetMaterial(n, shiftCoord, 3, vPrims);
if (sigma)
EffMat[3] = delta_ny / sigma;
else
EffMat[3] = 0;
length=delta_ny;
//shift up
shiftCoord[n] = coord[n]+delta*0.25;
shiftCoord[nP] = coord[nP]+deltaP*0.5;
shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
++loc_pos[n];
delta_ny = GetNodeWidth(n,loc_pos,true);
EffMat[2] += delta_ny / GetMaterial(n, shiftCoord, 2, vPrims);
sigma = GetMaterial(n, shiftCoord, 3, vPrims);
if (sigma)
EffMat[3] += delta_ny / sigma;
else
EffMat[3] = 0;
length+=delta_ny;
EffMat[2] = length * __MUE0__ / EffMat[2];
if (EffMat[3]) EffMat[3]=length / EffMat[3];
for (int n=0; n<4; ++n)
if (isnan(EffMat[n]) || isinf(EffMat[n]))
{
cerr << "Operator::" << __func__ << ": Error, An effective material parameter is not a valid result, this should NOT have happend... exit..." << endl;
cerr << ny << "@" << n << " : " << pos[0] << "," << pos[1] << "," << pos[2] << endl;
exit(0);
}
return true;
}
bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives *> vPrims) const
{
switch (m_MatAverageMethod)
{
case QuarterCell:
return AverageMatQuarterCell(ny, pos, EffMat, vPrims);
case CentralCell:
return AverageMatCellCenter(ny, pos, EffMat, vPrims);
default:
cerr << "Operator:: " << __func__ << ": Error, unknown material averaging method... exit" << endl;
exit(1);
}
return false;
}
bool Operator::Calc_LumpedElements()
{
vector<CSProperties*> props = CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT);
for (size_t i=0;i<props.size();++i)
{
CSPropLumpedElement* PLE = dynamic_cast<CSPropLumpedElement*>(props.at(i));
if (PLE==NULL)
return false; //sanity check: this should never happen!
vector<CSPrimitives*> prims = PLE->GetAllPrimitives();
for (size_t bn=0;bn<prims.size();++bn)
{
CSPrimBox* box = dynamic_cast<CSPrimBox*>(prims.at(bn));
if (box)
{ //calculate lumped element parameter
double C = PLE->GetCapacity();
if (C<=0)
C = NAN;
double R = PLE->GetResistance();
if (R<0)
R = NAN;
if ((isnan(R)) && (isnan(C)))
{
cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element R or C not specified! skipping. "
<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
continue;
}
int ny = PLE->GetDirection();
if ((ny<0) || (ny>2))
{
cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element direction is invalid! skipping. "
<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
continue;
}
int nyP = (ny+1)%3;
int nyPP = (ny+2)%3;
unsigned int uiStart[3];
unsigned int uiStop[3];
// snap to the native coordinate system
int Snap_Dimension = Operator::SnapBox2Mesh(box->GetStartCoord()->GetCoords(m_MeshType), box->GetStopCoord()->GetCoords(m_MeshType), uiStart, uiStop, false, true);
if (Snap_Dimension<=0)
{
if (Snap_Dimension>=-1)
cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element snapping failed! Dimension is: " << Snap_Dimension << " skipping. "
<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
// Snap_Dimension == -2 means outside the simulation domain --> no special warning, but box probably marked as unused!
continue;
}
if (uiStart[ny]==uiStop[ny])
{
cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element with zero (snapped) length is invalid! skipping. "
<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
continue;
}
//calculate geometric property for this lumped element
unsigned int pos[3];
double unitGC=0;
int ipos=0;
for (pos[ny]=uiStart[ny];pos[ny]<uiStop[ny];++pos[ny])
{
double unitGC_Plane=0;
for (pos[nyP]=uiStart[nyP];pos[nyP]<=uiStop[nyP];++pos[nyP])
{
for (pos[nyPP]=uiStart[nyPP];pos[nyPP]<=uiStop[nyPP];++pos[nyPP])
{
// capacity/conductivity in parallel: add values
unitGC_Plane += GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
}
}
//capacity/conductivity in series: add reciprocal values
unitGC += 1/unitGC_Plane;
}
unitGC = 1/unitGC;
bool caps = PLE->GetCaps();
double kappa = 0;
double epsilon = 0;
if (R>0)
kappa = 1 / R / unitGC;
if (C>0)
{
epsilon = C / unitGC;
if (epsilon< __EPS0__)
{
cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element capacity is too small for its size! skipping. "
<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
C = 0;
}
}
for (pos[ny]=uiStart[ny];pos[ny]<uiStop[ny];++pos[ny])
{
for (pos[nyP]=uiStart[nyP];pos[nyP]<=uiStop[nyP];++pos[nyP])
{
for (pos[nyPP]=uiStart[nyPP];pos[nyPP]<=uiStop[nyPP];++pos[nyPP])
{
ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
if (C>0)
EC_C[ny][ipos] = epsilon * GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
if (R>0)
EC_G[ny][ipos] = kappa * GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
if (R==0) //make lumped element a PEC if resistance is zero
{
SetVV(ny,pos[0],pos[1],pos[2], 0 );
SetVI(ny,pos[0],pos[1],pos[2], 0 );
}
else //recalculate operator inside the lumped element
Calc_ECOperatorPos(ny,pos);
}
}
}
// setup metal caps
if (caps)
{
for (pos[nyP]=uiStart[nyP];pos[nyP]<=uiStop[nyP];++pos[nyP])
{
for (pos[nyPP]=uiStart[nyPP];pos[nyPP]<=uiStop[nyPP];++pos[nyPP])
{
pos[ny]=uiStart[ny];
if (pos[nyP]<uiStop[nyP])
{
SetVV(nyP,pos[0],pos[1],pos[2], 0 );
SetVI(nyP,pos[0],pos[1],pos[2], 0 );
++m_Nr_PEC[nyP];
}
if (pos[nyPP]<uiStop[nyPP])
{
SetVV(nyPP,pos[0],pos[1],pos[2], 0 );
SetVI(nyPP,pos[0],pos[1],pos[2], 0 );
++m_Nr_PEC[nyPP];
}
pos[ny]=uiStop[ny];
if (pos[nyP]<uiStop[nyP])
{
SetVV(nyP,pos[0],pos[1],pos[2], 0 );
SetVI(nyP,pos[0],pos[1],pos[2], 0 );
++m_Nr_PEC[nyP];
}
if (pos[nyPP]<uiStop[nyPP])
{
SetVV(nyPP,pos[0],pos[1],pos[2], 0 );
SetVI(nyPP,pos[0],pos[1],pos[2], 0 );
++m_Nr_PEC[nyPP];
}
}
}
}
box->SetPrimitiveUsed(true);
}
else
cerr << "Operator::Calc_LumpedElements(): Warning: Primitves other than boxes are not supported for lumped elements! skipping "
<< prims.at(bn)->GetTypeName() << " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
}
}
return true;
}
void Operator::Init_EC()
{
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 FDTD_FLOAT[MainOp->GetSize()];
EC_G[n] = new FDTD_FLOAT[MainOp->GetSize()];
EC_L[n] = new FDTD_FLOAT[MainOp->GetSize()];
EC_R[n] = new FDTD_FLOAT[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;
}
}
}
bool Operator::Calc_EC()
{
if (CSX==NULL)
{
cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl;
return false;
}
MainOp->SetPos(0,0,0);
Calc_EC_Range(0,numLines[0]-1);
return true;
}
vector<CSPrimitives*> Operator::GetPrimitivesBoundBox(int posX, int posY, int posZ, CSProperties::PropertyType type) const
{
double boundBox[6];
int BBpos[3] = {posX, posY, posZ};
for (int n=0;n<3;++n)
{
if (BBpos[n]<0)
{
boundBox[2*n] = this->GetDiscLine(n,0);
boundBox[2*n+1] = this->GetDiscLine(n,numLines[n]-1);
}
else
{
boundBox[2*n] = this->GetDiscLine(n, max(0, BBpos[n]-1));
boundBox[2*n+1] = this->GetDiscLine(n, min(int(numLines[n])-1, BBpos[n]+1));
}
}
vector<CSPrimitives*> vPrim = this->CSX->GetPrimitivesByBoundBox(boundBox, true, type);
return vPrim;
}
void Operator::Calc_EC_Range(unsigned int xStart, unsigned int xStop)
{
// vector<CSPrimitives*> vPrims = this->CSX->GetAllPrimitives(true, CSProperties::MATERIAL);
unsigned int ipos;
unsigned int pos[3];
double inEC[4];
for (pos[0]=xStart; pos[0]<=xStop; ++pos[0])
{
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{
vector<CSPrimitives*> vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL);
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{
ipos = MainOp->GetPos(pos[0],pos[1],pos[2]);
for (int n=0; n<3; ++n)
{
Calc_ECPos(n,pos,inEC,vPrims);
EC_C[n][ipos]=inEC[0];
EC_G[n][ipos]=inEC[1];
EC_L[n][ipos]=inEC[2];
EC_R[n][ipos]=inEC[3];
}
}
}
}
}
void Operator::SetTimestepFactor(double factor)
{
if ((factor<=0) || (factor>1))
{
cerr << "Operator::SetTimestepFactor: Warning, invalid timestep factor, skipping!" << endl;
return;
}
cout << "Operator::SetTimestepFactor: Setting timestep factor to " << factor << endl;
m_TimeStepFactor=factor;
}
double Operator::CalcTimestep()
{
if (m_TimeStepVar==3)
return CalcTimestep_Var3(); //the biggest one for cartesian meshes
//variant 1 is default
return CalcTimestep_Var1();
}
////Berechnung nach Andreas Rennings Dissertation 2008, Seite 66, Formel 4.52
double Operator::CalcTimestep_Var1()
{
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;
dT=1e200;
double newT;
unsigned int pos[3];
unsigned int smallest_pos[3] = {0, 0, 0};
unsigned int smallest_n = 0;
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] );
if ((newT<dT) && (newT>0.0))
{
dT=newT;
smallest_pos[0]=pos[0];smallest_pos[1]=pos[1];smallest_pos[2]=pos[2];
smallest_n = n;
}
}
}
}
}
if (dT==0)
{
cerr << "Operator::CalcTimestep: Timestep is zero... this is not supposed to happen!!! exit!" << endl;
exit(3);
}
if (g_settings.GetVerboseLevel()>1)
{
cout << "Operator::CalcTimestep_Var1: Smallest timestep (" << dT << "s) found at position: " << smallest_n << " : " << smallest_pos[0] << ";" << smallest_pos[1] << ";" << smallest_pos[2] << endl;
}
return 0;
}
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
double Operator::CalcTimestep_Var3()
{
dT=1e200;
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;
double newT;
unsigned int pos[3];
unsigned int smallest_pos[3] = {0, 0, 0};
unsigned int smallest_n = 0;
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;
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])
{
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;
smallest_pos[0]=pos[0];smallest_pos[1]=pos[1];smallest_pos[2]=pos[2];
smallest_n = n;
}
}
}
}
}
if (dT==0)
{
cerr << "Operator::CalcTimestep: Timestep is zero... this is not supposed to happen!!! exit!" << endl;
exit(3);
}
if (g_settings.GetVerboseLevel()>1)
{
cout << "Operator::CalcTimestep_Var3: Smallest timestep (" << dT << "s) found at position: " << smallest_n << " : " << smallest_pos[0] << ";" << smallest_pos[1] << ";" << smallest_pos[2] << endl;
}
return 0;
}
bool Operator::CalcPEC()
{
m_Nr_PEC[0]=0;
m_Nr_PEC[1]=0;
m_Nr_PEC[2]=0;
CalcPEC_Range(0,numLines[0]-1,m_Nr_PEC);
CalcPEC_Curves();
return true;
}
void Operator::CalcPEC_Range(unsigned int startX, unsigned int stopX, unsigned int* counter)
{
double coord[3];
unsigned int pos[3];
for (pos[0]=startX; pos[0]<=stopX; ++pos[0])
{
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{
vector<CSPrimitives*> vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL));
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{
for (int n=0; n<3; ++n)
{
GetYeeCoords(n,pos,coord,false);
CSProperties* prop = CSX->GetPropertyByCoordPriority(coord, vPrims, true);
// CSProperties* old_prop = CSX->GetPropertyByCoordPriority(coord, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL), true);
// if (old_prop!=prop)
// {
// cerr << "CalcPEC_Range: " << old_prop << " vs " << prop << endl;
// exit(-1);
// }
if (prop)
{
if (prop->GetType()==CSProperties::METAL) //set to PEC
{
SetVV(n,pos[0],pos[1],pos[2], 0 );
SetVI(n,pos[0],pos[1],pos[2], 0 );
++counter[n];
}
}
}
}
}
}
}
void Operator::CalcPEC_Curves()
{
//special treatment for primitives of type curve (treated as wires)
double p1[3];
double p2[3];
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,m_MeshType);
curv->GetPoint(i,p2,m_MeshType);
path = FindPath(p1,p2);
if (path.dir.size()>0)
prim->SetPrimitiveUsed(true);
for (size_t t=0; t<path.dir.size(); ++t)
{
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 );
++m_Nr_PEC[path.dir.at(t)];
}
}
}
}
}
}
Operator_Ext_Excitation* Operator::GetExcitationExtension() const
{
//search for excitation extension
Operator_Ext_Excitation* Op_Ext_Exc=0;
for (size_t n=0; n<m_Op_exts.size(); ++n)
{
Op_Ext_Exc = dynamic_cast<Operator_Ext_Excitation*>(m_Op_exts.at(n));
if (Op_Ext_Exc)
break;
}
return Op_Ext_Exc;
}
void Operator::AddExtension(Operator_Extension* op_ext)
{
m_Op_exts.push_back(op_ext);
}
void Operator::DeleteExtension(Operator_Extension* op_ext)
{
for (size_t n=0;n<m_Op_exts.size();++n)
{
if (m_Op_exts.at(n)==op_ext)
{
m_Op_exts.erase(m_Op_exts.begin()+n);
return;
}
}
}
double Operator::CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const
{
double average_mesh_disc[3];
double c0 = __C0__/sqrt(GetBackgroundEpsR()*GetBackgroundMueR());
//calculate average mesh deltas
for (int n=0;n<3;++n)
{
average_mesh_disc[n] = fabs(GetDiscLine(n,start[n])-GetDiscLine(n,stop[n]))*GetGridDelta() / (fabs(stop[n]-start[n]));
}
// if propagation is in a Cartesian direction, return analytic solution
for (int n=0;n<3;++n)
{
int nP = (n+1)%3;
int nPP = (n+2)%3;
if ((fabs(propDir[n])==1) && (propDir[nP]==0) && (propDir[nPP]==0))
{
double kx = asin(average_mesh_disc[0]/c0/dT*sin(2*PI*freq*dT/2))*2/average_mesh_disc[0];
return 2*PI*freq/kx;
}
}
// else, do an newton iterative estimation
double k0=2*PI*freq/c0;
double k=k0;
double RHS = pow(sin(2*PI*freq*dT/2)/c0/dT,2);
double fk=1,fdk=0;
double old_phv=0;
double phv=c0;
double err_est = 1e-6;
int it_count=0;
while (fabs(old_phv-phv)>err_est)
{
++it_count;
old_phv=phv;
fk=0;
fdk=0;
for (int n=0;n<3;++n)
{
fk+= pow(sin(propDir[n]*k*average_mesh_disc[n]/2)/average_mesh_disc[n],2);
fdk+= propDir[n]*sin(propDir[n]*k*average_mesh_disc[n]/2)*cos(propDir[n]*k*average_mesh_disc[n]/2)/average_mesh_disc[n];
}
fk -= RHS;
k-=fk/fdk;
// do not allow a speed greater than c0 due to a numerical inaccuracy
if (k<k0)
k=k0;
phv=2*PI*freq/k;
//abort if iteration count is getting to high!
if (it_count>99)
{
cerr << "Operator::CalcNumericPhaseVelocity: Error, newton iteration estimation can't find a solution!!" << endl;
break;
}
}
if (g_settings.GetVerboseLevel()>1)
cerr << "Operator::CalcNumericPhaseVelocity: Newton iteration estimated solution: " << phv/__C0__ << "*c0 in " << it_count << " iterations." << endl;
return phv;
}
Reference in New Issue View Git Blame Copy Permalink