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