234 lines
7.2 KiB
C++
234 lines
7.2 KiB
C++
/*
|
|
* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
|
|
*
|
|
* This program is free software: you can redistribute it and/or modify
|
|
* it under the terms of the GNU General Public License as published by
|
|
* the Free Software Foundation, either version 3 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU General Public License
|
|
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#include "operator_ext_conductingsheet.h"
|
|
#include "tools/array_ops.h"
|
|
#include "tools/constants.h"
|
|
#include "cond_sheet_parameter.h"
|
|
|
|
Operator_Ext_ConductingSheet::Operator_Ext_ConductingSheet(Operator* op, double f_max) : Operator_Ext_LorentzMaterial(op)
|
|
{
|
|
m_f_max = f_max;
|
|
}
|
|
|
|
bool Operator_Ext_ConductingSheet::BuildExtension()
|
|
{
|
|
double dT = m_Op->GetTimestep();
|
|
unsigned int pos[] = {0,0,0};
|
|
double coord[3];
|
|
unsigned int numLines[3] = {m_Op->GetNumberOfLines(0),m_Op->GetNumberOfLines(1),m_Op->GetNumberOfLines(2)};
|
|
|
|
m_Order = 0;
|
|
vector<unsigned int> v_pos[3];
|
|
int ****tanDir = Create_N_3DArray<int>(numLines);
|
|
float ****Conductivity = Create_N_3DArray<float>(numLines);
|
|
float ****Thickness = Create_N_3DArray<float>(numLines);
|
|
|
|
CSPrimitives* cs_sheet = NULL;
|
|
double box[6];
|
|
int nP, nPP;
|
|
bool b_pos_on;
|
|
bool disable_pos;
|
|
for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
|
|
{
|
|
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
|
|
{
|
|
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
|
|
{
|
|
b_pos_on = false;
|
|
disable_pos = false;
|
|
// disable conducting sheet model inside the boundary conditions, especially inside a pml
|
|
for (int m=0;m<3;++m)
|
|
if ((pos[m]<=(m_Op->GetBCSize(2*m))) || (pos[m]>=(numLines[m]-m_Op->GetBCSize(2*m+1)-1)))
|
|
disable_pos = true;
|
|
|
|
for (int n=0; n<3; ++n)
|
|
{
|
|
nP = (n+1)%3;
|
|
nPP = (n+2)%3;
|
|
|
|
tanDir[n][pos[0]][pos[1]][pos[2]] = -1; //deactivate by default
|
|
Conductivity[n][pos[0]][pos[1]][pos[2]] = 0; //deactivate by default
|
|
Thickness[n][pos[0]][pos[1]][pos[2]] = 0; //deactivate by default
|
|
|
|
coord[0] = m_Op->GetDiscLine(0,pos[0]);
|
|
coord[1] = m_Op->GetDiscLine(1,pos[1]);
|
|
coord[2] = m_Op->GetDiscLine(2,pos[2]);
|
|
coord[n] = m_Op->GetDiscLine(n,pos[n],true); //pos of E_n
|
|
CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,CSProperties::CONDUCTINGSHEET, false, &cs_sheet);
|
|
if (prop)
|
|
{
|
|
CSPropConductingSheet* cs_prop = dynamic_cast<CSPropConductingSheet*>(prop);
|
|
if (cs_prop==NULL)
|
|
return false; //sanity check, this should never happen
|
|
if (cs_sheet==NULL)
|
|
return false; //sanity check, this should never happen
|
|
if (cs_sheet->GetDimension()!=2)
|
|
continue;
|
|
cs_sheet->SetPrimitiveUsed(true);
|
|
|
|
if (disable_pos)
|
|
{
|
|
m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 );
|
|
m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 );
|
|
continue;
|
|
}
|
|
|
|
Conductivity[n][pos[0]][pos[1]][pos[2]] = cs_prop->GetConductivity();
|
|
Thickness[n][pos[0]][pos[1]][pos[2]] = cs_prop->GetThickness();
|
|
|
|
if ((Conductivity[n][pos[0]][pos[1]][pos[2]]==0) || (Thickness[n][pos[0]][pos[1]][pos[2]]==0))
|
|
{
|
|
cerr << "Operator_Ext_ConductingSheet::BuildExtension: Warning: Zero conductivity or thickness detected... skipping node!" << endl;
|
|
continue;
|
|
}
|
|
|
|
cs_sheet->GetBoundBox(box);
|
|
if (box[2*nP]!=box[2*nP+1])
|
|
tanDir[n][pos[0]][pos[1]][pos[2]] = nP;
|
|
if (box[2*nPP]!=box[2*nPP+1])
|
|
tanDir[n][pos[0]][pos[1]][pos[2]] = nPP;
|
|
b_pos_on = true;
|
|
}
|
|
}
|
|
if (b_pos_on)
|
|
{
|
|
for (int n=0; n<3; ++n)
|
|
v_pos[n].push_back(pos[n]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
size_t numCS = v_pos[0].size();
|
|
if (numCS==0)
|
|
return false;
|
|
|
|
m_LM_Count.push_back(numCS);
|
|
m_LM_Count.push_back(numCS);
|
|
|
|
m_Order = 2;
|
|
m_volt_ADE_On = new bool[m_Order];
|
|
m_volt_ADE_On[0] = m_volt_ADE_On[1]=true;
|
|
m_curr_ADE_On = new bool[m_Order];
|
|
m_curr_ADE_On[0] = m_curr_ADE_On[1]=false;
|
|
|
|
m_LM_pos = new unsigned int**[m_Order];
|
|
m_LM_pos[0] = new unsigned int*[3];
|
|
m_LM_pos[1] = new unsigned int*[3];
|
|
|
|
v_int_ADE = new FDTD_FLOAT**[m_Order];
|
|
v_ext_ADE = new FDTD_FLOAT**[m_Order];
|
|
i_int_ADE = NULL;
|
|
i_ext_ADE = NULL;
|
|
|
|
v_int_ADE[0] = new FDTD_FLOAT*[3];
|
|
v_ext_ADE[0] = new FDTD_FLOAT*[3];
|
|
v_int_ADE[1] = new FDTD_FLOAT*[3];
|
|
v_ext_ADE[1] = new FDTD_FLOAT*[3];
|
|
|
|
for (int n=0; n<3; ++n)
|
|
{
|
|
m_LM_pos[0][n] = new unsigned int[numCS];
|
|
m_LM_pos[1][n] = new unsigned int[numCS];
|
|
for (unsigned int i=0; i<numCS; ++i)
|
|
{
|
|
m_LM_pos[0][n][i] = v_pos[n].at(i);
|
|
m_LM_pos[1][n][i] = v_pos[n].at(i);
|
|
}
|
|
v_int_ADE[0][n] = new FDTD_FLOAT[numCS];
|
|
v_int_ADE[1][n] = new FDTD_FLOAT[numCS];
|
|
v_ext_ADE[0][n] = new FDTD_FLOAT[numCS];
|
|
v_ext_ADE[1][n] = new FDTD_FLOAT[numCS];
|
|
}
|
|
|
|
unsigned int index;
|
|
float w_stop = m_f_max*2*PI;
|
|
float Omega_max=0;
|
|
float G,L1,L2,R1,R2,Lmin;
|
|
float G0, w0;
|
|
float wtl; //width to length factor
|
|
float factor=1;
|
|
int t_dir=0; //tangential sheet direction
|
|
unsigned int tpos[] = {0,0,0};
|
|
unsigned int optParaPos;
|
|
for (unsigned int i=0;i<numCS;++i)
|
|
{
|
|
pos[0]=m_LM_pos[0][0][i];pos[1]=m_LM_pos[0][1][i];pos[2]=m_LM_pos[0][2][i];
|
|
tpos[0]=pos[0];tpos[1]=pos[1];tpos[2]=pos[2];
|
|
index = m_Op->MainOp->SetPos(pos[0],pos[1],pos[2]);
|
|
for (int n=0;n<3;++n)
|
|
{
|
|
t_dir = tanDir[n][pos[0]][pos[1]][pos[2]];
|
|
G0 = Conductivity[n][pos[0]][pos[1]][pos[2]]*Thickness[n][pos[0]][pos[1]][pos[2]];
|
|
w0 = 8.0/ G0 / Thickness[n][pos[0]][pos[1]][pos[2]]/__MUE0__;
|
|
Omega_max = w_stop/w0;
|
|
for (optParaPos=0;optParaPos<numOptPara;++optParaPos)
|
|
if (omega_stop[optParaPos]>Omega_max)
|
|
break;
|
|
if (optParaPos>=numOptPara)
|
|
{
|
|
cerr << "Operator_Ext_ConductingSheet::BuildExtension(): Error, conductor thickness, conductivity or max. simulation frequency of interest is too high! Check parameter!" << endl;
|
|
cerr << " --> max f: " << m_f_max << "Hz, Conductivity: " << Conductivity[n][pos[0]][pos[1]][pos[2]] << "S/m, Thickness " << Thickness[n][pos[0]][pos[1]][pos[2]]*1e6 << "um" << endl;
|
|
optParaPos = numOptPara-1;
|
|
}
|
|
v_int_ADE[0][n][i]=0;
|
|
v_ext_ADE[0][n][i]=0;
|
|
v_int_ADE[1][n][i]=0;
|
|
v_ext_ADE[1][n][i]=0;
|
|
if (t_dir>=0)
|
|
{
|
|
wtl = m_Op->GetEdgeLength(n,pos)/m_Op->GetNodeWidth(t_dir,pos);
|
|
factor = 1;
|
|
if (tanDir[t_dir][tpos[0]][tpos[1]][tpos[2]]<0)
|
|
factor = 2;
|
|
--tpos[t_dir];
|
|
if (tanDir[t_dir][tpos[0]][tpos[1]][tpos[2]]<0)
|
|
factor = 2;
|
|
|
|
if (n==1)
|
|
factor=1;
|
|
|
|
L1 = l1[optParaPos]/G0/w0*factor;
|
|
L2 = l2[optParaPos]/G0/w0*factor;
|
|
R1 = r1[optParaPos]/G0*factor;
|
|
R2 = r2[optParaPos]/G0*factor;
|
|
G = G0*g[optParaPos]/factor;
|
|
|
|
L1*=wtl;
|
|
L2*=wtl;
|
|
R1*=wtl;
|
|
R2*=wtl;
|
|
G/=wtl;
|
|
|
|
Lmin = L1;
|
|
if (L2<L1)
|
|
Lmin = L2;
|
|
m_Op->EC_G[n][index]= G;
|
|
m_Op->EC_C[n][index]= dT*dT/4.0*(16.0/Lmin + 1/L1 + 1/L2);
|
|
m_Op->Calc_ECOperatorPos(n,pos);
|
|
|
|
v_int_ADE[0][n][i]=(2.0*L1-dT*R1)/(2.0*L1+dT*R1);
|
|
v_ext_ADE[0][n][i]=dT/(L1+dT*R1/2.0)*m_Op->GetVI(n,pos[0],pos[1],pos[2]);
|
|
v_int_ADE[1][n][i]=(2.0*L2-dT*R2)/(2.0*L2+dT*R2);
|
|
v_ext_ADE[1][n][i]=dT/(L2+dT*R2/2.0)*m_Op->GetVI(n,pos[0],pos[1],pos[2]);
|
|
}
|
|
}
|
|
}
|
|
return true;
|
|
}
|