1 /*
2 * Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 #include "operator_ext_conductingsheet.h"
19 #include "tools/array_ops.h"
20 #include "tools/constants.h"
21 #include "cond_sheet_parameter.h"
22
23 #include "CSPropConductingSheet.h"
24
Operator_Ext_ConductingSheet(Operator * op,double f_max)25 Operator_Ext_ConductingSheet::Operator_Ext_ConductingSheet(Operator* op, double f_max) : Operator_Ext_LorentzMaterial(op)
26 {
27 m_f_max = f_max;
28 }
29
Operator_Ext_ConductingSheet(Operator * op,Operator_Ext_ConductingSheet * op_ext)30 Operator_Ext_ConductingSheet::Operator_Ext_ConductingSheet(Operator* op, Operator_Ext_ConductingSheet* op_ext) : Operator_Ext_LorentzMaterial(op, op_ext)
31 {
32 m_f_max = op_ext->m_f_max;
33 }
34
Clone(Operator * op)35 Operator_Extension* Operator_Ext_ConductingSheet::Clone(Operator* op)
36 {
37 if (dynamic_cast<Operator_Ext_ConductingSheet*>(this)==NULL)
38 return NULL;
39 return new Operator_Ext_ConductingSheet(op, this);
40 }
41
BuildExtension()42 bool Operator_Ext_ConductingSheet::BuildExtension()
43 {
44 double dT = m_Op->GetTimestep();
45 unsigned int pos[] = {0,0,0};
46 double coord[3];
47 unsigned int numLines[3] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)};
48
49 m_Order = 0;
50 vector<unsigned int> v_pos[3];
51 int ****tanDir = Create_N_3DArray<int>(numLines);
52 float ****Conductivity = Create_N_3DArray<float>(numLines);
53 float ****Thickness = Create_N_3DArray<float>(numLines);
54
55 CSPrimitives* cs_sheet = NULL;
56 double box[6];
57 int nP, nPP;
58 bool b_pos_on;
59 bool disable_pos;
60 for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
61 {
62 for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
63 {
64 vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL));
65 for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
66 {
67 b_pos_on = false;
68 disable_pos = false;
69 // disable conducting sheet model inside the boundary conditions, especially inside a pml
70 for (int m=0;m<3;++m)
71 if ((pos[m]<=(unsigned int)m_Op->GetBCSize(2*m)) || (pos[m]>=(numLines[m]-m_Op->GetBCSize(2*m+1)-1)))
72 disable_pos = true;
73
74 for (int n=0; n<3; ++n)
75 {
76 nP = (n+1)%3;
77 nPP = (n+2)%3;
78
79 tanDir[n][pos[0]][pos[1]][pos[2]] = -1; //deactivate by default
80 Conductivity[n][pos[0]][pos[1]][pos[2]] = 0; //deactivate by default
81 Thickness[n][pos[0]][pos[1]][pos[2]] = 0; //deactivate by default
82
83 if (m_Op->GetYeeCoords(n,pos,coord,false)==false)
84 continue;
85
86 // Ez at r==0 not supported --> set to PEC
87 if (m_CC_R0_included && (n==2) && (pos[0]==0))
88 disable_pos = true;
89
90 // CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), false, &cs_sheet);
91 CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, false, &cs_sheet);
92 CSPropConductingSheet* cs_prop = dynamic_cast<CSPropConductingSheet*>(prop);
93 if (cs_prop)
94 {
95 if (cs_sheet==NULL)
96 return false; //sanity check, this should never happen
97 if (cs_sheet->GetDimension()!=2)
98 {
99 cerr << "Operator_Ext_ConductingSheet::BuildExtension: A conducting sheet primitive (ID: " << cs_sheet->GetID() << ") with dimension: " << cs_sheet->GetDimension() << " found, fallback to PEC!" << endl;
100 m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 );
101 m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 );
102 ++m_Op->m_Nr_PEC[n];
103 continue;
104 }
105 cs_sheet->SetPrimitiveUsed(true);
106
107 if (disable_pos)
108 {
109 m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 );
110 m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 );
111 ++m_Op->m_Nr_PEC[n];
112 continue;
113 }
114
115 Conductivity[n][pos[0]][pos[1]][pos[2]] = cs_prop->GetConductivity();
116 Thickness[n][pos[0]][pos[1]][pos[2]] = cs_prop->GetThickness();
117
118 if ((Conductivity[n][pos[0]][pos[1]][pos[2]]<=0) || (Thickness[n][pos[0]][pos[1]][pos[2]]<=0))
119 {
120 cerr << "Operator_Ext_ConductingSheet::BuildExtension: Warning: Zero conductivity or thickness detected... fallback to PEC!" << endl;
121 m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 );
122 m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 );
123 ++m_Op->m_Nr_PEC[n];
124 continue;
125 }
126
127 cs_sheet->GetBoundBox(box);
128 if (box[2*nP]!=box[2*nP+1])
129 tanDir[n][pos[0]][pos[1]][pos[2]] = nP;
130 if (box[2*nPP]!=box[2*nPP+1])
131 tanDir[n][pos[0]][pos[1]][pos[2]] = nPP;
132 b_pos_on = true;
133 }
134 }
135 if (b_pos_on)
136 {
137 for (int n=0; n<3; ++n)
138 v_pos[n].push_back(pos[n]);
139 }
140 }
141 }
142 }
143
144 size_t numCS = v_pos[0].size();
145 if (numCS==0)
146 return false;
147
148 m_LM_Count.push_back(numCS);
149 m_LM_Count.push_back(numCS);
150
151 m_Order = 2;
152 m_volt_ADE_On = new bool[m_Order];
153 m_volt_ADE_On[0] = m_volt_ADE_On[1]=true;
154 m_curr_ADE_On = new bool[m_Order];
155 m_curr_ADE_On[0] = m_curr_ADE_On[1]=false;
156
157 m_volt_Lor_ADE_On = new bool[m_Order];
158 m_volt_Lor_ADE_On[0] = m_volt_Lor_ADE_On[1]=false;
159 m_curr_Lor_ADE_On = new bool[m_Order];
160 m_curr_Lor_ADE_On[0] = m_curr_Lor_ADE_On[1]=false;
161
162 m_LM_pos = new unsigned int**[m_Order];
163 m_LM_pos[0] = new unsigned int*[3];
164 m_LM_pos[1] = new unsigned int*[3];
165
166 v_int_ADE = new FDTD_FLOAT**[m_Order];
167 v_ext_ADE = new FDTD_FLOAT**[m_Order];
168
169 v_int_ADE[0] = new FDTD_FLOAT*[3];
170 v_ext_ADE[0] = new FDTD_FLOAT*[3];
171 v_int_ADE[1] = new FDTD_FLOAT*[3];
172 v_ext_ADE[1] = new FDTD_FLOAT*[3];
173
174 for (int n=0; n<3; ++n)
175 {
176 m_LM_pos[0][n] = new unsigned int[numCS];
177 m_LM_pos[1][n] = new unsigned int[numCS];
178 for (unsigned int i=0; i<numCS; ++i)
179 {
180 m_LM_pos[0][n][i] = v_pos[n].at(i);
181 m_LM_pos[1][n][i] = v_pos[n].at(i);
182 }
183 v_int_ADE[0][n] = new FDTD_FLOAT[numCS];
184 v_int_ADE[1][n] = new FDTD_FLOAT[numCS];
185 v_ext_ADE[0][n] = new FDTD_FLOAT[numCS];
186 v_ext_ADE[1][n] = new FDTD_FLOAT[numCS];
187 }
188
189 unsigned int index;
190 float w_stop = m_f_max*2*PI;
191 float Omega_max=0;
192 float G,L1,L2,R1,R2,Lmin;
193 float G0, w0;
194 float wtl; //width to length factor
195 float factor=1;
196 int t_dir=0; //tangential sheet direction
197 unsigned int tpos[] = {0,0,0};
198 unsigned int optParaPos;
199 for (unsigned int i=0;i<numCS;++i)
200 {
201 pos[0]=m_LM_pos[0][0][i];pos[1]=m_LM_pos[0][1][i];pos[2]=m_LM_pos[0][2][i];
202 tpos[0]=pos[0];tpos[1]=pos[1];tpos[2]=pos[2];
203 index = m_Op->MainOp->SetPos(pos[0],pos[1],pos[2]);
204 for (int n=0;n<3;++n)
205 {
206 tpos[0]=pos[0];tpos[1]=pos[1];tpos[2]=pos[2];
207 t_dir = tanDir[n][pos[0]][pos[1]][pos[2]];
208 G0 = Conductivity[n][pos[0]][pos[1]][pos[2]]*Thickness[n][pos[0]][pos[1]][pos[2]];
209 w0 = 8.0/ G0 / Thickness[n][pos[0]][pos[1]][pos[2]]/__MUE0__;
210 Omega_max = w_stop/w0;
211 for (optParaPos=0;optParaPos<numOptPara;++optParaPos)
212 if (omega_stop[optParaPos]>Omega_max)
213 break;
214 if (optParaPos>=numOptPara)
215 {
216 cerr << "Operator_Ext_ConductingSheet::BuildExtension(): Error, conductor thickness, conductivity or max. simulation frequency of interest is too high! Check parameter!" << endl;
217 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;
218 optParaPos = numOptPara-1;
219 }
220 v_int_ADE[0][n][i]=0;
221 v_ext_ADE[0][n][i]=0;
222 v_int_ADE[1][n][i]=0;
223 v_ext_ADE[1][n][i]=0;
224 if (t_dir>=0)
225 {
226 wtl = m_Op->GetEdgeLength(n,pos)/m_Op->GetNodeWidth(t_dir,pos);
227 factor = 1;
228 if (tanDir[t_dir][tpos[0]][tpos[1]][tpos[2]]<0)
229 factor = 2;
230 --tpos[t_dir];
231 if (tanDir[t_dir][tpos[0]][tpos[1]][tpos[2]]<0)
232 factor = 2;
233
234 L1 = l1[optParaPos]/G0/w0*factor;
235 L2 = l2[optParaPos]/G0/w0*factor;
236 R1 = r1[optParaPos]/G0*factor;
237 R2 = r2[optParaPos]/G0*factor;
238 G = G0*g[optParaPos]/factor;
239
240 L1*=wtl;
241 L2*=wtl;
242 R1*=wtl;
243 R2*=wtl;
244 G/=wtl;
245
246 Lmin = L1;
247 if (L2<L1)
248 Lmin = L2;
249 m_Op->EC_G[n][index]= G;
250 m_Op->EC_C[n][index]= dT*dT/4.0*(16.0/Lmin + 1/L1 + 1/L2);
251 m_Op->Calc_ECOperatorPos(n,pos);
252
253 v_int_ADE[0][n][i]=(2.0*L1-dT*R1)/(2.0*L1+dT*R1);
254 v_ext_ADE[0][n][i]=dT/(L1+dT*R1/2.0)*m_Op->GetVI(n,pos[0],pos[1],pos[2]);
255 v_int_ADE[1][n][i]=(2.0*L2-dT*R2)/(2.0*L2+dT*R2);
256 v_ext_ADE[1][n][i]=dT/(L2+dT*R2/2.0)*m_Op->GetVI(n,pos[0],pos[1],pos[2]);
257 }
258 }
259 }
260 return true;
261 }
262