From 3fc2a41af9e5d2628aafe2130af56b2ae682e2e2 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Sat, 28 Dec 2013 21:02:49 +0100 Subject: [PATCH] operator: change how to average material to allow for overloaded cylindrical handling Signed-off-by: Thorsten Liebig --- FDTD/operator.cpp | 160 ++++++++++++++++++++++++++++++++++++++++++---- FDTD/operator.h | 15 ++++- openems.cpp | 5 +- 3 files changed, 164 insertions(+), 16 deletions(-) diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index f83ca74..f34c7a5 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -95,7 +95,7 @@ void Operator::Init() m_Exc = 0; m_TimeStepFactor = 1; - m_MatCellShiftFaktor = 0.25; + SetMaterialAvgMethod(QuarterCell); } void Operator::Delete() @@ -487,6 +487,18 @@ Grid_Path Operator::FindPath(double start[], double stop[]) return path; } +void Operator::SetMaterialAvgMethod(MatAverageMethods method) +{ + switch (method) + { + default: + case QuarterCell: + return SetQuarterCellMaterialAvg(); + case CentralCell: + return SetCellConstantMaterial(); + } +} + double Operator::GetNumberCells() const { if (numLines) @@ -1265,7 +1277,113 @@ double Operator::GetMaterial(int ny, const double* coords, int MatType, bool mar } } -bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) const +bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat) 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)*A_n; + EffMat[1] += GetMaterial(n, coord, 1)*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)*A_n; + EffMat[1] += GetMaterial(n, coord, 1)*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)*A_n; + EffMat[1] += GetMaterial(n, coord, 1)*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)*A_n; + EffMat[1] += GetMaterial(n, coord, 1)*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); + sigma = GetMaterial(n, coord, 3); + 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); + sigma = GetMaterial(n, coord, 3); + 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) const { int n=ny; double coord[3]; @@ -1289,8 +1407,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c //******************************* epsilon,kappa averaging *****************************// //shift up-right shiftCoord[n] = coord[n]+delta*0.5; - shiftCoord[nP] = coord[nP]+deltaP*m_MatCellShiftFaktor; - shiftCoord[nPP] = coord[nPP]+deltaPP*m_MatCellShiftFaktor; + 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)*A_n; EffMat[1] = GetMaterial(n, shiftCoord, 1)*A_n; @@ -1298,8 +1416,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c //shift up-left shiftCoord[n] = coord[n]+delta*0.5; - shiftCoord[nP] = coord[nP]-deltaP_M*m_MatCellShiftFaktor; - shiftCoord[nPP] = coord[nPP]+deltaPP*m_MatCellShiftFaktor; + 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); @@ -1309,8 +1427,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c //shift down-right shiftCoord[n] = coord[n]+delta*0.5; - shiftCoord[nP] = coord[nP]+deltaP*m_MatCellShiftFaktor; - shiftCoord[nPP] = coord[nPP]-deltaPP_M*m_MatCellShiftFaktor; + 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); @@ -1320,8 +1438,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c //shift down-left shiftCoord[n] = coord[n]+delta*0.5; - shiftCoord[nP] = coord[nP]-deltaP_M*m_MatCellShiftFaktor; - shiftCoord[nPP] = coord[nPP]-deltaPP_M*m_MatCellShiftFaktor; + 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)*A_n; @@ -1338,7 +1456,7 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c double length=0; //shift down - shiftCoord[n] = coord[n]-delta_M*m_MatCellShiftFaktor; + 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]; @@ -1352,7 +1470,7 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c length=delta_ny; //shift up - shiftCoord[n] = coord[n]+delta*m_MatCellShiftFaktor; + shiftCoord[n] = coord[n]+delta*0.25; shiftCoord[nP] = coord[nP]+deltaP*0.5; shiftCoord[nPP] = coord[nPP]+deltaPP*0.5; ++loc_pos[n]; @@ -1371,13 +1489,29 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c 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; + 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) const +{ + switch (m_MatAverageMethod) + { + case QuarterCell: + return AverageMatQuarterCell(ny, pos, EffMat); + case CentralCell: + return AverageMatCellCenter(ny, pos, EffMat); + default: + cerr << "Operator:: " << __func__ << ": Error, unknown material averaging method... exit" << endl; + exit(1); + } + return false; +} + bool Operator::Calc_LumpedElements() { vector props = CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT); diff --git a/FDTD/operator.h b/FDTD/operator.h index f0cf824..8dd259b 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -42,6 +42,8 @@ class Operator : public Operator_Base public: enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4}; + enum MatAverageMethods {QuarterCell=0, CentralCell=1}; + //! Create a new operator static Operator* New(); virtual ~Operator(); @@ -84,8 +86,14 @@ public: //! Choose a time step method (0=auto, 1=CFL, 3=Rennings) void SetTimeStepMethod(int var) {m_TimeStepVar=var;} + //! Set the material averaging method /sa MatAverageMethods + void SetMaterialAvgMethod(MatAverageMethods method); + + //! Set material averaging method to the advanced quarter cell material interpolation (default) + void SetQuarterCellMaterialAvg() {m_MatAverageMethod=QuarterCell;} + //! Set operator to assume a constant material inside a cell (material probing in the cell center) - void SetCellConstantMaterial() {m_MatCellShiftFaktor=0.5;} + void SetCellConstantMaterial() {m_MatAverageMethod=CentralCell;} virtual double GetNumberCells() const; @@ -219,11 +227,14 @@ protected: //! Get the material at a given coordinate, direction and type from CSX (internal use only) virtual double GetMaterial(int ny, const double coords[3], int MatType, bool markAsUsed=true) const; - double m_MatCellShiftFaktor; + MatAverageMethods m_MatAverageMethod; //! Calculate the effective/averaged material properties at the given position and direction ny. virtual bool Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) const; + virtual bool AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat) const; + virtual bool AverageMatQuarterCell(int ny, const unsigned int* pos, double* EffMat) const; + //! Calc operator at certain \a pos virtual void Calc_ECOperatorPos(int n, unsigned int* pos); diff --git a/openems.cpp b/openems.cpp index 3af6a57..921b31d 100644 --- a/openems.cpp +++ b/openems.cpp @@ -675,9 +675,12 @@ int openEMS::SetupFDTD(const char* file) if (SetupOperator(FDTD_Opts)==false) return 2; + // default material averaging is quarter cell averaging + FDTD_Op->SetQuarterCellMaterialAvg(); + + // check for cell constant material averaging if (FDTD_Opts->QueryIntAttribute("CellConstantMaterial",&ihelp)==TIXML_SUCCESS) m_CellConstantMaterial=(ihelp==1); - if (m_CellConstantMaterial) { FDTD_Op->SetCellConstantMaterial();