caiyu/openEMS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Operator: replacement for GetOriginalNumLines + use full or simple mesh for snapping

Browse Source 
Note: these are a lot and dangerous changes --> require a lot of testing

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>
This commit is contained in:
Thorsten Liebig 2013-02-06 16:33:12 +01:00
parent 0cb33c3753
commit 183ea9f776
18 changed files with 100 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
Common/operator_base.h
View File

@ -41,8 +41,8 @@ public:
//! Get the number of timesteps satisfying the nyquist condition (may depend on the excitation) //! Get the number of timesteps satisfying the nyquist condition (may depend on the excitation)
virtual unsigned int GetNumberOfNyquistTimesteps() const =0; virtual unsigned int GetNumberOfNyquistTimesteps() const =0;
//! Returns the number of lines as needed for post-processing etc. (for the engine, use GetOriginalNumLines()) //! Returns the number of lines as needed for post-processing etc.
virtual unsigned int GetNumberOfLines(int ny) const =0; virtual unsigned int GetNumberOfLines(int ny, bool full=false) const =0;
//! Get the name for the given direction: 0 -> x, 1 -> y, 2 -> z //! Get the name for the given direction: 0 -> x, 1 -> y, 2 -> z
virtual std::string GetDirName(int ny) const; virtual std::string GetDirName(int ny) const;
@ -76,7 +76,7 @@ public:
virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const =0; virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const =0;
//! Snap the given coodinates to mesh indices, return box dimension //! Snap the given coodinates to mesh indices, return box dimension
virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool* inside=NULL) const =0; virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool fullMesh=false, bool* inside=NULL) const =0;
//! Snap a given box to the operator mesh, uiStart will be always <= uiStop //! Snap a given box to the operator mesh, uiStart will be always <= uiStop
/*! /*!
@ -88,7 +88,7 @@ public:
\param[in] SnapMethod Snapping method, 0=snap to closest line, 1/(2)=snap such that given box is inside (outside) the snapped lines \param[in] SnapMethod Snapping method, 0=snap to closest line, 1/(2)=snap such that given box is inside (outside) the snapped lines
\return returns the box dimension or -1 if box is not inside the simulation domain \return returns the box dimension or -1 if box is not inside the simulation domain
*/ */
virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const =0; virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, bool fullMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const =0;
//! Set the boundary conditions //! Set the boundary conditions
virtual void SetBoundaryCondition(int* BCs) {for (int n=0; n<6; ++n) m_BC[n]=BCs[n];} virtual void SetBoundaryCondition(int* BCs) {for (int n=0; n<6; ++n) m_BC[n]=BCs[n];}

2
Common/processing.cpp
View File

@ -170,7 +170,7 @@ void Processing::AddFrequency(vector<double> *freqs)
void Processing::DefineStartStopCoord(double* dstart, double* dstop) void Processing::DefineStartStopCoord(double* dstart, double* dstop)
{ {
m_Dimension = Op->SnapBox2Mesh(dstart,dstop,start,stop,m_dualMesh,m_SnapMethod, m_start_inside, m_stop_inside); m_Dimension = Op->SnapBox2Mesh(dstart,dstop,start,stop,m_dualMesh,false,m_SnapMethod, m_start_inside, m_stop_inside);
if (m_Dimension<0) if (m_Dimension<0)
{ {
cerr << "Processing::DefineStartStopCoord: Warning in " << m_Name << " (" << GetProcessingName() << ") : Box is outside the field domain!! Disabling" << endl; cerr << "Processing::DefineStartStopCoord: Warning in " << m_Name << " (" << GetProcessingName() << ") : Box is outside the field domain!! Disabling" << endl;

2
FDTD/engine.cpp
View File

@ -36,7 +36,7 @@ Engine::Engine(const Operator* op)
numTS = 0; numTS = 0;
Op = op; Op = op;
for (int n=0; n<3; ++n) for (int n=0; n<3; ++n)
numLines[n] = Op->GetOriginalNumLines(n); numLines[n] = Op->GetNumberOfLines(n, true);
volt=NULL; volt=NULL;
curr=NULL; curr=NULL;
} }

6
FDTD/engine_interface_fdtd.cpp
View File

@ -81,7 +81,7 @@ double* Engine_Interface_FDTD::GetRawInterpolatedField(const unsigned int* pos,
{ {
nP = (n+1)%3; nP = (n+1)%3;
nPP = (n+2)%3; nPP = (n+2)%3;
if ((pos[0]==m_Op->GetOriginalNumLines(0)-1) || (pos[1]==m_Op->GetOriginalNumLines(1)-1) || (pos[2]==m_Op->GetOriginalNumLines(2)-1)) if ((pos[0]==m_Op->GetNumberOfLines(0,true)-1) || (pos[1]==m_Op->GetNumberOfLines(1,true)-1) || (pos[2]==m_Op->GetNumberOfLines(2,true)-1))
{ {
out[n] = 0; //electric field outside the field domain is always zero out[n] = 0; //electric field outside the field domain is always zero
continue; continue;
@ -119,7 +119,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
{ {
nP = (n+1)%3; nP = (n+1)%3;
nPP = (n+2)%3; nPP = (n+2)%3;
if ((pos[0]==m_Op->GetOriginalNumLines(0)-1) || (pos[1]==m_Op->GetOriginalNumLines(1)-1) || (pos[2]==m_Op->GetOriginalNumLines(2)-1) || (pos[nP]==0) || (pos[nPP]==0)) if ((pos[0]==m_Op->GetNumberOfLines(0,true)-1) || (pos[1]==m_Op->GetNumberOfLines(1,true)-1) || (pos[2]==m_Op->GetNumberOfLines(2,true)-1) || (pos[nP]==0) || (pos[nPP]==0))
{ {
out[n] = 0; out[n] = 0;
continue; continue;
@ -140,7 +140,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
{ {
delta = m_Op->GetEdgeLength(n,iPos,true); delta = m_Op->GetEdgeLength(n,iPos,true);
out[n] = m_Eng->GetCurr(n,iPos); out[n] = m_Eng->GetCurr(n,iPos);
if ((pos[n]>=m_Op->GetOriginalNumLines(n)-1)) if ((pos[n]>=m_Op->GetNumberOfLines(n,true)-1))
{ {
out[n] = 0; //magnetic field on the outer boundaries is always zero out[n] = 0; //magnetic field on the outer boundaries is always zero
continue; continue;

2
FDTD/extensions/engine_ext_cylinder.cpp
View File

@ -29,7 +29,7 @@ Engine_Ext_Cylinder::Engine_Ext_Cylinder(Operator_Ext_Cylinder* op_ext) : Engine
CC_R0_included = op_ext->CC_R0_included; CC_R0_included = op_ext->CC_R0_included;
for (int n=0; n<3; ++n) for (int n=0; n<3; ++n)
numLines[n] = op_ext->m_Op->GetOriginalNumLines(n); numLines[n] = op_ext->m_Op->GetNumberOfLines(n,true);
//this cylindrical extension should be executed first? //this cylindrical extension should be executed first?
m_Priority = ENG_EXT_PRIO_CYLINDER; m_Priority = ENG_EXT_PRIO_CYLINDER;

2
FDTD/extensions/operator_ext_conductingsheet.cpp
View File

@ -44,7 +44,7 @@ bool Operator_Ext_ConductingSheet::BuildExtension()
double dT = m_Op->GetTimestep(); double dT = m_Op->GetTimestep();
unsigned int pos[] = {0,0,0}; unsigned int pos[] = {0,0,0};
double coord[3]; double coord[3];
unsigned int numLines[3] = {m_Op->GetOriginalNumLines(0),m_Op->GetOriginalNumLines(1),m_Op->GetOriginalNumLines(2)}; unsigned int numLines[3] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)};
m_Order = 0; m_Order = 0;
vector<unsigned int> v_pos[3]; vector<unsigned int> v_pos[3];

10
FDTD/extensions/operator_ext_cylinder.cpp
View File

@ -50,19 +50,19 @@ bool Operator_Ext_Cylinder::BuildExtension()
if (CC_R0_included==false) if (CC_R0_included==false)
return true; return true;
vv_R0 = new FDTD_FLOAT[m_Op->GetOriginalNumLines(2)]; vv_R0 = new FDTD_FLOAT[m_Op->GetNumberOfLines(2,true)];
vi_R0 = new FDTD_FLOAT[m_Op->GetOriginalNumLines(2)]; vi_R0 = new FDTD_FLOAT[m_Op->GetNumberOfLines(2,true)];
unsigned int pos[3]; unsigned int pos[3];
double coord[3]; double coord[3];
double inEC[4]; double inEC[4];
double dT = m_Op->GetTimestep(); double dT = m_Op->GetTimestep();
pos[0]=0; pos[0]=0;
for (pos[2]=0; pos[2]<m_Op->GetOriginalNumLines(2); ++pos[2]) for (pos[2]=0; pos[2]<m_Op->GetNumberOfLines(2,true); ++pos[2])
{ {
double C=0; double C=0;
double G=0; double G=0;
for (pos[1]=0; pos[1]<m_Op->GetOriginalNumLines(1)-2; ++pos[1]) for (pos[1]=0; pos[1]<m_Op->GetNumberOfLines(1,true)-2; ++pos[1])
{ {
m_Op_Cyl->Calc_ECPos(2,pos,inEC); m_Op_Cyl->Calc_ECPos(2,pos,inEC);
C+=inEC[0]; C+=inEC[0];
@ -72,7 +72,7 @@ bool Operator_Ext_Cylinder::BuildExtension()
vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C); vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C);
vi_R0[pos[2]] = (dT/C)/(1+dT*G/2/C); vi_R0[pos[2]] = (dT/C)/(1+dT*G/2/C);
for (unsigned int i=0; i<m_Op->GetOriginalNumLines(1); ++i) for (unsigned int i=0; i<m_Op->GetNumberOfLines(1,true); ++i)
{ {
m_Op->EC_C[2][m_Op->MainOp->SetPos(0,i,pos[2])] = C; m_Op->EC_C[2][m_Op->MainOp->SetPos(0,i,pos[2])] = C;
m_Op->EC_G[2][m_Op->MainOp->SetPos(0,i,pos[2])] = G; m_Op->EC_G[2][m_Op->MainOp->SetPos(0,i,pos[2])] = G;

6
FDTD/extensions/operator_ext_excitation.cpp
View File

@ -137,7 +137,7 @@ bool Operator_Ext_Excitation::BuildExtension()
CSProperties* prop=NULL; CSProperties* prop=NULL;
int priority=0; int priority=0;
unsigned int numLines[] = {m_Op->GetOriginalNumLines(0),m_Op->GetOriginalNumLines(1),m_Op->GetOriginalNumLines(2)}; unsigned int numLines[] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)};
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
@ -243,8 +243,8 @@ bool Operator_Ext_Excitation::BuildExtension()
{ {
for (size_t i=1; i<curv->GetNumberOfPoints(); ++i) for (size_t i=1; i<curv->GetNumberOfPoints(); ++i)
{ {
curv->GetPoint(i-1,p1); curv->GetPoint(i-1,p1,m_Op->m_MeshType);
curv->GetPoint(i,p2); curv->GetPoint(i,p2,m_Op->m_MeshType);
path = m_Op->FindPath(p1,p2); path = m_Op->FindPath(p1,p2);
if (path.dir.size()>0) if (path.dir.size()>0)
prim->SetPrimitiveUsed(true); prim->SetPrimitiveUsed(true);

2
FDTD/extensions/operator_ext_lorentzmaterial.cpp
View File

@ -88,7 +88,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
double dT = m_Op->GetTimestep(); double dT = m_Op->GetTimestep();
unsigned int pos[] = {0,0,0}; unsigned int pos[] = {0,0,0};
double coord[3]; double coord[3];
unsigned int numLines[3] = {m_Op->GetOriginalNumLines(0),m_Op->GetOriginalNumLines(1),m_Op->GetOriginalNumLines(2)}; unsigned int numLines[3] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)};
CSPropLorentzMaterial* mat = NULL; CSPropLorentzMaterial* mat = NULL;
double w_plasma,t_relax; double w_plasma,t_relax;

8
FDTD/extensions/operator_ext_mur_abc.cpp
View File

@ -104,12 +104,12 @@ void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny)
} }
else else
{ {
m_LineNr = m_Op->GetOriginalNumLines(m_ny)-1; m_LineNr = m_Op->GetNumberOfLines(m_ny,true)-1;
m_LineNr_Shift = m_Op->GetOriginalNumLines(m_ny) - 2; m_LineNr_Shift = m_Op->GetNumberOfLines(m_ny,true) - 2;
} }
m_numLines[0] = m_Op->GetOriginalNumLines(m_nyP); m_numLines[0] = m_Op->GetNumberOfLines(m_nyP,true);
m_numLines[1] = m_Op->GetOriginalNumLines(m_nyPP); m_numLines[1] = m_Op->GetNumberOfLines(m_nyPP,true);
m_Mur_Coeff_nyP = Create2DArray<FDTD_FLOAT>(m_numLines); m_Mur_Coeff_nyP = Create2DArray<FDTD_FLOAT>(m_numLines);
m_Mur_Coeff_nyPP = Create2DArray<FDTD_FLOAT>(m_numLines); m_Mur_Coeff_nyPP = Create2DArray<FDTD_FLOAT>(m_numLines);

2
FDTD/extensions/operator_ext_tfsf.cpp
View File

@ -176,7 +176,7 @@ bool Operator_Ext_TFST::BuildExtension()
m_ActiveDir[n][0]=false; m_ActiveDir[n][0]=false;
else else
m_ActiveDir[n][0]=true; m_ActiveDir[n][0]=true;
if (m_Stop[n]==m_Op->GetOriginalNumLines(n)-1) if (m_Stop[n]==m_Op->GetNumberOfLines(n,true)-1)
m_ActiveDir[n][1]=false; m_ActiveDir[n][1]=false;
else else
m_ActiveDir[n][1]=true; m_ActiveDir[n][1]=true;

32
FDTD/extensions/operator_ext_upml.cpp
View File

@ -78,7 +78,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
//check if mesh is large enough to support the pml //check if mesh is large enough to support the pml
for (int n=0; n<3; ++n) for (int n=0; n<3; ++n)
if ( (size[2*n]*(BC[2*n]==3)+size[2*n+1]*(BC[2*n+1]==3)) >= op->GetOriginalNumLines(n) ) if ( (size[2*n]*(BC[2*n]==3)+size[2*n+1]*(BC[2*n+1]==3)) >= op->GetNumberOfLines(n,true) )
{ {
cerr << "Operator_Ext_UPML::Create_UPML: Warning: Not enough lines in direction: " << n << ", resetting to PEC" << endl; cerr << "Operator_Ext_UPML::Create_UPML: Warning: Not enough lines in direction: " << n << ", resetting to PEC" << endl;
BC[2*n]=0; BC[2*n]=0;
@ -131,7 +131,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
Operator_Ext_UPML* op_ext_upml=NULL; Operator_Ext_UPML* op_ext_upml=NULL;
unsigned int start[3]={0 ,0 ,0}; unsigned int start[3]={0 ,0 ,0};
unsigned int stop[3] ={op->GetOriginalNumLines(0)-1,op->GetOriginalNumLines(1)-1,op->GetOriginalNumLines(2)-1}; unsigned int stop[3] ={op->GetNumberOfLines(0,true)-1,op->GetNumberOfLines(1,true)-1,op->GetNumberOfLines(2,true)-1};
//create a pml in x-direction over the full width of yz-space //create a pml in x-direction over the full width of yz-space
if (BC[0]==3) if (BC[0]==3)
@ -148,8 +148,8 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
{ {
op_ext_upml = new Operator_Ext_UPML(op); op_ext_upml = new Operator_Ext_UPML(op);
op_ext_upml->SetGradingFunction(gradFunc); op_ext_upml->SetGradingFunction(gradFunc);
start[0]=op->GetOriginalNumLines(0)-1-size[1]; start[0]=op->GetNumberOfLines(0,true)-1-size[1];
stop[0] =op->GetOriginalNumLines(0)-1; stop[0] =op->GetNumberOfLines(0,true)-1;
op_ext_upml->SetBoundaryCondition(BC, size); op_ext_upml->SetBoundaryCondition(BC, size);
op_ext_upml->SetRange(start,stop); op_ext_upml->SetRange(start,stop);
op->AddExtension(op_ext_upml); op->AddExtension(op_ext_upml);
@ -157,7 +157,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
//create a pml in y-direction over the xz-space (if a pml in x-direction already exists, skip that corner regions) //create a pml in y-direction over the xz-space (if a pml in x-direction already exists, skip that corner regions)
start[0]=(size[0]+1)*(BC[0]==3); start[0]=(size[0]+1)*(BC[0]==3);
stop[0] =op->GetOriginalNumLines(0)-1-(size[0]+1)*(BC[1]==3); stop[0] =op->GetNumberOfLines(0,true)-1-(size[0]+1)*(BC[1]==3);
if (BC[2]==3) if (BC[2]==3)
{ {
@ -173,8 +173,8 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
{ {
op_ext_upml = new Operator_Ext_UPML(op); op_ext_upml = new Operator_Ext_UPML(op);
op_ext_upml->SetGradingFunction(gradFunc); op_ext_upml->SetGradingFunction(gradFunc);
start[1]=op->GetOriginalNumLines(1)-1-size[3]; start[1]=op->GetNumberOfLines(1,true)-1-size[3];
stop[1] =op->GetOriginalNumLines(1)-1; stop[1] =op->GetNumberOfLines(1,true)-1;
op_ext_upml->SetBoundaryCondition(BC, size); op_ext_upml->SetBoundaryCondition(BC, size);
op_ext_upml->SetRange(start,stop); op_ext_upml->SetRange(start,stop);
op->AddExtension(op_ext_upml); op->AddExtension(op_ext_upml);
@ -182,7 +182,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
//create a pml in z-direction over the xy-space (if a pml in x- and/or y-direction already exists, skip that corner/edge regions) //create a pml in z-direction over the xy-space (if a pml in x- and/or y-direction already exists, skip that corner/edge regions)
start[1]=(size[2]+1)*(BC[2]==3); start[1]=(size[2]+1)*(BC[2]==3);
stop[1] =op->GetOriginalNumLines(1)-1-(size[3]+1)*(BC[3]==3); stop[1] =op->GetNumberOfLines(1,true)-1-(size[3]+1)*(BC[3]==3);
//exclude x-lines that does not belong to the base multi-grid operator; //exclude x-lines that does not belong to the base multi-grid operator;
Operator_CylinderMultiGrid* op_cyl_MG = dynamic_cast<Operator_CylinderMultiGrid*>(op); Operator_CylinderMultiGrid* op_cyl_MG = dynamic_cast<Operator_CylinderMultiGrid*>(op);
@ -203,8 +203,8 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
{ {
op_ext_upml = new Operator_Ext_UPML(op); op_ext_upml = new Operator_Ext_UPML(op);
op_ext_upml->SetGradingFunction(gradFunc); op_ext_upml->SetGradingFunction(gradFunc);
start[2]=op->GetOriginalNumLines(2)-1-size[5]; start[2]=op->GetNumberOfLines(2,true)-1-size[5];
stop[2] =op->GetOriginalNumLines(2)-1; stop[2] =op->GetNumberOfLines(2,true)-1;
op_ext_upml->SetBoundaryCondition(BC, size); op_ext_upml->SetBoundaryCondition(BC, size);
op_ext_upml->SetRange(start,stop); op_ext_upml->SetRange(start,stop);
op->AddExtension(op_ext_upml); op->AddExtension(op_ext_upml);
@ -220,7 +220,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
for (int n=0; n<2; ++n) for (int n=0; n<2; ++n)
{ {
start[n]=0; start[n]=0;
stop[n]=op_child->GetOriginalNumLines(n)-1; stop[n]=op_child->GetNumberOfLines(n,true)-1;
} }
if (op_cyl_MG) if (op_cyl_MG)
@ -240,8 +240,8 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
{ {
op_ext_upml = new Operator_Ext_UPML(op_child); op_ext_upml = new Operator_Ext_UPML(op_child);
op_ext_upml->SetGradingFunction(gradFunc); op_ext_upml->SetGradingFunction(gradFunc);
start[2]=op->GetOriginalNumLines(2)-1-size[5]; start[2]=op->GetNumberOfLines(2,true)-1-size[5];
stop[2] =op->GetOriginalNumLines(2)-1; stop[2] =op->GetNumberOfLines(2,true)-1;
op_ext_upml->SetBoundaryCondition(BC, size); op_ext_upml->SetBoundaryCondition(BC, size);
op_ext_upml->SetRange(start,stop); op_ext_upml->SetRange(start,stop);
op_child->AddExtension(op_ext_upml); op_child->AddExtension(op_ext_upml);
@ -320,10 +320,10 @@ void Operator_Ext_UPML::CalcGradingKappa(int ny, unsigned int pos[3], double Zm,
else else
kappa_i[n] = 0; kappa_i[n] = 0;
} }
else if ((pos[n] >= m_Op->GetOriginalNumLines(n) -1 -m_Size[2*n+1]) && (m_BC[2*n+1]==3)) //upper pml in n-dir else if ((pos[n] >= m_Op->GetNumberOfLines(n,true) -1 -m_Size[2*n+1]) && (m_BC[2*n+1]==3)) //upper pml in n-dir
{ {
width = (m_Op->GetDiscLine(n,m_Op->GetOriginalNumLines(n)-1) - m_Op->GetDiscLine(n,m_Op->GetOriginalNumLines(n)-m_Size[2*n+1]-1))*m_Op->GetGridDelta(); width = (m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-1) - m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-m_Size[2*n+1]-1))*m_Op->GetGridDelta();
depth = width - (m_Op->GetDiscLine(n,m_Op->GetOriginalNumLines(n)-1) - m_Op->GetDiscLine(n,pos[n]))*m_Op->GetGridDelta(); depth = width - (m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-1) - m_Op->GetDiscLine(n,pos[n]))*m_Op->GetGridDelta();
if ((m_Op_Cyl) && (n==1)) if ((m_Op_Cyl) && (n==1))
{ {

58
FDTD/operator.cpp
View File

@ -227,14 +227,14 @@ double Operator::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
return GetNodeArea(ny, uiPos, dualMesh); return GetNodeArea(ny, uiPos, dualMesh);
} }
unsigned int Operator::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh) const unsigned int Operator::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh, bool fullMesh) const
{ {
inside = false; inside = false;
if ((ny<0) || (ny>2)) if ((ny<0) || (ny>2))
return 0; return 0;
if (coord<GetDiscLine(ny,0)) if (coord<GetDiscLine(ny,0))
return 0; return 0;
unsigned int numLines = GetNumberOfLines(ny); unsigned int numLines = GetNumberOfLines(ny, fullMesh);
if (coord>GetDiscLine(ny,numLines-1)) if (coord>GetDiscLine(ny,numLines-1))
return numLines-1; return numLines-1;
inside=true; inside=true;
@ -258,24 +258,21 @@ unsigned int Operator::SnapToMeshLine(int ny, double coord, bool &inside, bool d
return 0; return 0;
} }
bool Operator::SnapToMesh(const double* dcoord, unsigned int* uicoord, bool dualMesh, bool* inside) const bool Operator::SnapToMesh(const double* dcoord, unsigned int* uicoord, bool dualMesh, bool fullMesh, bool* inside) const
{ {
bool meshInside=false; bool meshInside=false;
bool ok=true; bool ok=true;
for (int n=0; n<3; ++n) for (int n=0; n<3; ++n)
{ {
uicoord[n] = SnapToMeshLine(n,dcoord[n],meshInside,dualMesh); uicoord[n] = SnapToMeshLine(n,dcoord[n],meshInside,dualMesh,fullMesh);
ok &= meshInside; ok &= meshInside;
if (inside) if (inside)
inside[n]=meshInside; inside[n]=meshInside;
} }
// 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;
return ok; return ok;
} }
int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, int SnapMethod, bool* bStartIn, bool* bStopIn) const 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]; double l_start[3], l_stop[3];
for (int n=0;n<3;++n) for (int n=0;n<3;++n)
@ -283,15 +280,15 @@ int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int
l_start[n] = fmin(start[n],stop[n]); l_start[n] = fmin(start[n],stop[n]);
l_stop[n] = fmax(start[n], stop[n]); l_stop[n] = fmax(start[n], stop[n]);
double min = GetDiscLine(n,0); double min = GetDiscLine(n,0);
double max = GetDiscLine(n,GetNumberOfLines(n)-1); 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)) ) if ( ((l_start[n]<min) && (l_stop[n]<min)) || ((l_start[n]>max) && (l_stop[n]>max)) )
{ {
return -2; return -2;
} }
} }
SnapToMesh(l_start, uiStart, dualMesh, bStartIn); SnapToMesh(l_start, uiStart, dualMesh, fullMesh, bStartIn);
SnapToMesh(l_stop, uiStop, dualMesh, bStopIn); SnapToMesh(l_stop, uiStop, dualMesh, fullMesh, bStopIn);
int iDim = 0; int iDim = 0;
if (SnapMethod==0) if (SnapMethod==0)
@ -309,7 +306,7 @@ int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int
{ {
if ((GetDiscLine( n, uiStart[n], dualMesh ) > l_start[n]) && (uiStart[n]>0)) if ((GetDiscLine( n, uiStart[n], dualMesh ) > l_start[n]) && (uiStart[n]>0))
--uiStart[n]; --uiStart[n];
if ((GetDiscLine( n, uiStop[n], dualMesh ) < l_stop[n]) && (uiStop[n]<GetNumberOfLines(n)-1)) if ((GetDiscLine( n, uiStop[n], dualMesh ) < l_stop[n]) && (uiStop[n]<GetNumberOfLines(n, fullMesh)-1))
++uiStop[n]; ++uiStop[n];
} }
if (uiStop[n]>uiStart[n]) if (uiStop[n]>uiStart[n])
@ -323,7 +320,7 @@ int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int
{ {
if (uiStop[n]>uiStart[n]) if (uiStop[n]>uiStart[n])
{ {
if ((GetDiscLine( n, uiStart[n], dualMesh ) < l_start[n]) && (uiStart[n]<GetNumberOfLines(n)-1)) if ((GetDiscLine( n, uiStart[n], dualMesh ) < l_start[n]) && (uiStart[n]<GetNumberOfLines(n, fullMesh)-1))
++uiStart[n]; ++uiStart[n];
if ((GetDiscLine( n, uiStop[n], dualMesh ) > l_stop[n]) && (uiStop[n]>0)) if ((GetDiscLine( n, uiStop[n], dualMesh ) > l_stop[n]) && (uiStop[n]>0))
--uiStop[n]; --uiStop[n];
@ -341,24 +338,23 @@ int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int
struct Operator::Grid_Path Operator::FindPath(double start[], double stop[]) struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
{ {
struct Grid_Path path; struct Grid_Path path;
// double dV[] = {stop[0]-start[0],stop[1]-start[1],stop[2]-start[2]};
unsigned int uiStart[3],uiStop[3],currPos[3]; unsigned int uiStart[3],uiStop[3],currPos[3];
SnapToMesh(start,uiStart);
SnapToMesh(stop,uiStop); SnapToMesh(start,uiStart,false,true);
SnapToMesh(stop,uiStop,false,true);
currPos[0]=uiStart[0]; currPos[0]=uiStart[0];
currPos[1]=uiStart[1]; currPos[1]=uiStart[1];
currPos[2]=uiStart[2]; currPos[2]=uiStart[2];
double meshStart[] = {discLines[0][uiStart[0]], discLines[1][uiStart[1]], discLines[2][uiStart[2]]}; double meshStart[3] = {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]]}; double meshStop[3] = {discLines[0][uiStop[0]], discLines[1][uiStop[1]], discLines[2][uiStop[2]]};
bool UpDir = false; bool UpDir = false;
double foot=0,dist=0,minFoot=0,minDist=0; double foot=0,dist=0,minFoot=0,minDist=0;
int minDir=0; int minDir=0;
unsigned int minPos[3]; unsigned int minPos[3];
double startFoot,stopFoot,currFoot; double startFoot,stopFoot,currFoot;
Point_Line_Distance(meshStart,start,stop,startFoot,dist); Point_Line_Distance(meshStart,start,stop,startFoot,dist, m_MeshType);
Point_Line_Distance(meshStop,start,stop,stopFoot,dist); Point_Line_Distance(meshStop,start,stop,stopFoot,dist, m_MeshType);
currFoot=startFoot; currFoot=startFoot;
minFoot=startFoot; minFoot=startFoot;
double P[3]; double P[3];
@ -374,7 +370,7 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
if (((int)currPos[n]-1)>=0) if (((int)currPos[n]-1)>=0)
{ {
P[n] = discLines[n][currPos[n]-1]; P[n] = discLines[n][currPos[n]-1];
Point_Line_Distance(P,start,stop,foot,dist); Point_Line_Distance(P,start,stop,foot,dist, m_MeshType);
if ((foot>currFoot) && (dist<minDist)) if ((foot>currFoot) && (dist<minDist))
{ {
minFoot=foot; minFoot=foot;
@ -386,7 +382,7 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
if ((currPos[n]+1)<numLines[n]) if ((currPos[n]+1)<numLines[n])
{ {
P[n] = discLines[n][currPos[n]+1]; P[n] = discLines[n][currPos[n]+1];
Point_Line_Distance(P,start,stop,foot,dist); Point_Line_Distance(P,start,stop,foot,dist, m_MeshType);
if ((foot>currFoot) && (dist<minDist)) if ((foot>currFoot) && (dist<minDist))
{ {
minFoot=foot; minFoot=foot;
@ -408,6 +404,14 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
currPos[minDir]+=-1; currPos[minDir]+=-1;
minPos[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[0].push_back(minPos[0]);
path.posPath[1].push_back(minPos[1]); path.posPath[1].push_back(minPos[1]);
path.posPath[2].push_back(minPos[2]); path.posPath[2].push_back(minPos[2]);
@ -1324,7 +1328,7 @@ bool Operator::Calc_LumpedElements()
unsigned int uiStart[3]; unsigned int uiStart[3];
unsigned int uiStop[3]; unsigned int uiStop[3];
// snap to the native coordinate system // snap to the native coordinate system
int Snap_Dimension = Operator::SnapBox2Mesh(box->GetStartCoord()->GetCoords(m_MeshType), box->GetStopCoord()->GetCoords(m_MeshType), uiStart, uiStop); 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<=0)
{ {
if (Snap_Dimension>=-1) if (Snap_Dimension>=-1)
@ -1717,8 +1721,8 @@ void Operator::CalcPEC_Curves()
{ {
for (size_t i=1; i<curv->GetNumberOfPoints(); ++i) for (size_t i=1; i<curv->GetNumberOfPoints(); ++i)
{ {
curv->GetPoint(i-1,p1); curv->GetPoint(i-1,p1,m_MeshType);
curv->GetPoint(i,p2); curv->GetPoint(i,p2,m_MeshType);
path = FindPath(p1,p2); path = FindPath(p1,p2);
if (path.dir.size()>0) if (path.dir.size()>0)
prim->SetPrimitiveUsed(true); prim->SetPrimitiveUsed(true);

13
FDTD/operator.h
View File

@ -85,10 +85,7 @@ public:
virtual unsigned int GetNumberOfNyquistTimesteps() const {return m_Exc->GetNyquistNum();} virtual unsigned int GetNumberOfNyquistTimesteps() const {return m_Exc->GetNyquistNum();}
virtual unsigned int GetNumberOfLines(int ny) const {return numLines[ny];} virtual unsigned int GetNumberOfLines(int ny, bool full=false) const {UNUSED(full);return numLines[ny];}
//! Returns the number of lines as needed for the engine etc. (for post-processing etc, use GetNumLines())
virtual unsigned int GetOriginalNumLines(int ny) const {return numLines[ny];}
virtual void ShowStat() const; virtual void ShowStat() const;
virtual void ShowExtStat() const; virtual void ShowExtStat() const;
@ -127,13 +124,13 @@ public:
*/ */
virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetNodeArea(ny,pos,dualMesh);} virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetNodeArea(ny,pos,dualMesh);}
virtual unsigned int SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh=false) const; virtual unsigned int SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh=false, bool fullMesh=false) const;
//! Snap the given coodinates to mesh indices //! Snap the given coodinates to mesh indices
virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool* inside=NULL) const; virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool fullMesh=false, bool* inside=NULL) const;
//! Snap a given box to the FDTD mesh //! Snap a given box to the FDTD mesh
virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const; virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, bool fullMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const;
virtual void AddExtension(Operator_Extension* op_ext); virtual void AddExtension(Operator_Extension* op_ext);
virtual void DeleteExtension(Operator_Extension* op_ext); virtual void DeleteExtension(Operator_Extension* op_ext);
@ -168,7 +165,7 @@ protected:
vector<unsigned int> posPath[3]; vector<unsigned int> posPath[3];
vector<unsigned short> dir; vector<unsigned short> dir;
}; };
struct Grid_Path FindPath(double start[], double stop[]); virtual struct Grid_Path FindPath(double start[], double stop[]);
// debug // debug
virtual void DumpOperator2File(string filename); virtual void DumpOperator2File(string filename);

23
FDTD/operator_cylinder.cpp
View File

@ -88,13 +88,16 @@ int Operator_Cylinder::CalcECOperator( DebugFlags debugFlags )
return rc; return rc;
} }
inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny, bool full) const
{ {
if (full)
return Operator_Multithread::GetNumberOfLines(ny, full);
//this is necessary for a correct field processing... cylindrical engine has to reset this by adding +1 //this is necessary for a correct field processing... cylindrical engine has to reset this by adding +1
if (CC_closedAlpha && ny==1) if (CC_closedAlpha && ny==1)
return Operator_Multithread::GetNumberOfLines(ny)-2; return Operator_Multithread::GetNumberOfLines(ny, true)-2;
return Operator_Multithread::GetNumberOfLines(ny); return Operator_Multithread::GetNumberOfLines(ny, full);
} }
string Operator_Cylinder::GetDirName(int ny) const string Operator_Cylinder::GetDirName(int ny) const
@ -201,10 +204,10 @@ double Operator_Cylinder::GetEdgeArea(int ny, const unsigned int pos[3], bool du
return GetEdgeLength(1,pos,!dualMesh) * GetEdgeLength(2,pos,!dualMesh); return GetEdgeLength(1,pos,!dualMesh) * GetEdgeLength(2,pos,!dualMesh);
} }
double Operator_Cylinder::FitToAlphaRange(double a_coord) const double Operator_Cylinder::FitToAlphaRange(double a_coord, bool fullMesh) const
{ {
double min = GetDiscLine(1,0); double min = GetDiscLine(1,0);
double max = GetDiscLine(1,GetOriginalNumLines(1)-1); double max = GetDiscLine(1,GetNumberOfLines(1, fullMesh)-1);
if ((a_coord>=min) && (a_coord<=max)) if ((a_coord>=min) && (a_coord<=max))
return a_coord; return a_coord;
while (a_coord<min) while (a_coord<min)
@ -227,17 +230,17 @@ double Operator_Cylinder::FitToAlphaRange(double a_coord) const
return a_coord; return a_coord;
} }
unsigned int Operator_Cylinder::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh) const unsigned int Operator_Cylinder::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh, bool fullMesh) const
{ {
if (ny==1) if (ny==1)
coord=FitToAlphaRange(coord); coord=FitToAlphaRange(coord);
return Operator_Multithread::SnapToMeshLine(ny, coord, inside, dualMesh); return Operator_Multithread::SnapToMeshLine(ny, coord, inside, dualMesh, fullMesh);
} }
int Operator_Cylinder::SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, int SnapMethod, bool* bStartIn, bool* bStopIn) const int Operator_Cylinder::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 a_min = GetDiscLine(1,0); double a_min = GetDiscLine(1,0);
double a_max = GetDiscLine(1,GetNumberOfLines(1)-1); double a_max = GetDiscLine(1,GetNumberOfLines(1,fullMesh)-1);
double a_size = stop[1] - start[1]; double a_size = stop[1] - start[1];
double a_center = FitToAlphaRange(0.5*(stop[1]+start[1])); double a_center = FitToAlphaRange(0.5*(stop[1]+start[1]));
@ -254,7 +257,7 @@ int Operator_Cylinder::SnapBox2Mesh(const double* start, const double* stop, uns
double l_start[3] = {start[0], a_start, start[2]}; double l_start[3] = {start[0], a_start, start[2]};
double l_stop[3] = {stop[0] , a_stop , stop[2] }; double l_stop[3] = {stop[0] , a_stop , stop[2] };
return Operator_Multithread::SnapBox2Mesh(l_start, l_stop, uiStart, uiStop, dualMesh, SnapMethod, bStartIn, bStopIn); return Operator_Multithread::SnapBox2Mesh(l_start, l_stop, uiStart, uiStop, dualMesh, fullMesh, SnapMethod, bStartIn, bStopIn);
} }
bool Operator_Cylinder::SetupCSXGrid(CSRectGrid* grid) bool Operator_Cylinder::SetupCSXGrid(CSRectGrid* grid)

8
FDTD/operator_cylinder.h
View File

@ -40,7 +40,7 @@ public:
virtual void ApplyMagneticBC(bool* dirs); virtual void ApplyMagneticBC(bool* dirs);
virtual unsigned int GetNumberOfLines(int ny) const; virtual unsigned int GetNumberOfLines(int ny, bool full=false) const;
//! Get the name for the given direction: 0 -> rho, 1 -> alpha, 2 -> z //! Get the name for the given direction: 0 -> rho, 1 -> alpha, 2 -> z
virtual string GetDirName(int ny) const; virtual string GetDirName(int ny) const;
@ -71,12 +71,12 @@ public:
*/ */
virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const; virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const;
virtual double FitToAlphaRange(double a_coord) const; virtual double FitToAlphaRange(double a_coord, bool fullMesh=false) const;
virtual unsigned int SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh=false) const; virtual unsigned int SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh=false, bool fullMesh=false) const;
//! Snap a given box to the FDTD mesh //! Snap a given box to the FDTD mesh
virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const; virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, bool fullMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const;
bool GetClosedAlpha() const {return CC_closedAlpha;} bool GetClosedAlpha() const {return CC_closedAlpha;}
bool GetR0Included() const {return CC_R0_included;} bool GetR0Included() const {return CC_R0_included;}

7
FDTD/operator_mpi.cpp
View File

@ -156,10 +156,13 @@ void Operator_MPI::SetOriginalMesh(CSRectGrid* orig_Mesh)
} }
} }
unsigned int Operator_MPI::GetNumberOfLines(int ny) const unsigned int Operator_MPI::GetNumberOfLines(int ny, bool fullMesh) const
{ {
if (fullMesh)
return Operator_SSE_Compressed::GetNumberOfLines(ny,fullMesh);
if ((!m_MPI_Enabled) || (m_NeighborUp[ny]<0)) if ((!m_MPI_Enabled) || (m_NeighborUp[ny]<0))
return Operator_SSE_Compressed::GetNumberOfLines(ny); return Operator_SSE_Compressed::GetNumberOfLines(ny,fullMesh);
return Operator_SSE_Compressed::GetNumberOfLines(ny)-1; return Operator_SSE_Compressed::GetNumberOfLines(ny)-1;
} }

2
FDTD/operator_mpi.h
View File

@ -50,7 +50,7 @@ public:
virtual void SetSplitPos(int ny, unsigned int pos) {m_SplitPos[ny]=pos;} virtual void SetSplitPos(int ny, unsigned int pos) {m_SplitPos[ny]=pos;}
virtual void SetOriginalMesh(CSRectGrid* orig_Mesh); virtual void SetOriginalMesh(CSRectGrid* orig_Mesh);
virtual unsigned int GetNumberOfLines(int ny) const; virtual unsigned int GetNumberOfLines(int ny, bool fullMesh=false) const;
virtual void AddExtension(Operator_Extension* op_ext); virtual void AddExtension(Operator_Extension* op_ext);