1 /*
2 * Copyright (C) 2010 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 "engine.h"
19 #include "engine_cylinder.h"
20 #include "Common/processfields.h"
21 #include "operator_cylinder.h"
22 #include "extensions/operator_extension.h"
23 #include "extensions/operator_ext_cylinder.h"
24 #include "tools/useful.h"
25
New(unsigned int numThreads)26 Operator_Cylinder* Operator_Cylinder::New(unsigned int numThreads)
27 {
28 cout << "Create cylindrical FDTD operator" << endl;
29 Operator_Cylinder* op = new Operator_Cylinder();
30 op->setNumThreads(numThreads);
31 op->Init();
32 return op;
33 }
34
Operator_Cylinder()35 Operator_Cylinder::Operator_Cylinder() : Operator_Multithread()
36 {
37 m_MeshType = CYLINDRICAL;
38 m_Cyl_Ext = NULL;
39 }
40
~Operator_Cylinder()41 Operator_Cylinder::~Operator_Cylinder()
42 {
43 }
44
45
CreateEngine()46 Engine* Operator_Cylinder::CreateEngine()
47 {
48 //! create a special cylindrical-engine
49 m_Engine = Engine_Cylinder::New(this, m_numThreads);
50 return m_Engine;
51 }
52
Init()53 void Operator_Cylinder::Init()
54 {
55 CC_closedAlpha = false;
56 CC_R0_included = false;
57 Operator_Multithread::Init();
58 }
59
GetRawDiscDelta(int ny,const int pos) const60 double Operator_Cylinder::GetRawDiscDelta(int ny, const int pos) const
61 {
62 if (CC_closedAlpha && ny==1 && pos==-1)
63 return (discLines[1][numLines[1]-2] - discLines[1][numLines[1]-3]);
64 if (CC_closedAlpha && ny==1 && (pos==(int)numLines[ny]-1))
65 return (discLines[1][2] - discLines[1][1]);
66
67 return Operator_Multithread::GetRawDiscDelta(ny,pos);
68 }
69
GetMaterial(int ny,const double * coords,int MatType,vector<CSPrimitives * > vPrims,bool markAsUsed) const70 double Operator_Cylinder::GetMaterial(int ny, const double* coords, int MatType, vector<CSPrimitives*> vPrims, bool markAsUsed) const
71 {
72 double l_coords[] = {coords[0],coords[1],coords[2]};
73 if (CC_closedAlpha && (coords[1]>GetDiscLine(1,0,false)+2*PI))
74 l_coords[1]-=2*PI;
75 if (CC_closedAlpha && (coords[1]<GetDiscLine(1,0,false)))
76 l_coords[1] += 2*PI;
77 return Operator_Multithread::GetMaterial(ny,l_coords,MatType,vPrims,markAsUsed);
78 }
79
CalcECOperator(DebugFlags debugFlags)80 int Operator_Cylinder::CalcECOperator( DebugFlags debugFlags )
81 {
82 // debugs only work with the native vector dumps
83 bool natDump = g_settings.NativeFieldDumps();
84 g_settings.SetNativeFieldDumps(true);
85 int rc = Operator_Multithread::CalcECOperator(debugFlags);
86 // reset original settings
87 g_settings.SetNativeFieldDumps(natDump);
88 return rc;
89 }
90
CalcTimestep()91 double Operator_Cylinder::CalcTimestep()
92 {
93 if (discLines[0][0]==0.0)
94 // use conservative timestep for a mesh including the r==0 singularity
95 m_TimeStepVar = 1;
96
97 return Operator::CalcTimestep();
98 }
99
GetNumberOfLines(int ny,bool full) const100 inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny, bool full) const
101 {
102 if (full)
103 return Operator_Multithread::GetNumberOfLines(ny, full);
104
105 //this is necessary for a correct field processing... cylindrical engine has to reset this by adding +1
106 if (CC_closedAlpha && ny==1)
107 return Operator_Multithread::GetNumberOfLines(ny, true)-2;
108
109 return Operator_Multithread::GetNumberOfLines(ny, full);
110 }
111
GetDirName(int ny) const112 string Operator_Cylinder::GetDirName(int ny) const
113 {
114 if (ny==0) return "rho";
115 if (ny==1) return "alpha";
116 if (ny==2) return "z";
117 return "";
118 }
119
GetYeeCoords(int ny,unsigned int pos[3],double * coords,bool dualMesh) const120 bool Operator_Cylinder::GetYeeCoords(int ny, unsigned int pos[3], double* coords, bool dualMesh) const
121 {
122 bool ret = Operator_Multithread::GetYeeCoords(ny,pos,coords,dualMesh);
123
124 if (CC_closedAlpha && (coords[1]>=GetDiscLine(1,0,false)+2*PI))
125 coords[1]-=2*PI;
126 if (CC_closedAlpha && (coords[1]<GetDiscLine(1,0,false)))
127 coords[1]+=2*PI;
128
129 return ret;
130 }
131
GetNodeWidth(int ny,const unsigned int pos[3],bool dualMesh) const132 double Operator_Cylinder::GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh) const
133 {
134 if ((ny<0) || (ny>2)) return 0.0;
135 if (pos[ny]>=numLines[ny]) return 0.0;
136 double width = Operator_Multithread::GetEdgeLength(ny,pos,!dualMesh);
137 if (ny==1)
138 width *= GetDiscLine(0,pos[0],dualMesh);
139 return width;
140 }
141
GetNodeWidth(int ny,const int pos[3],bool dualMesh) const142 double Operator_Cylinder::GetNodeWidth(int ny, const int pos[3], bool dualMesh) const
143 {
144 if ( (pos[0]<0) || (pos[1]<0 && CC_closedAlpha==false) || (pos[2]<0) )
145 return 0.0;
146
147 unsigned int uiPos[]={(unsigned int)pos[0],(unsigned int)pos[1],(unsigned int)pos[2]};
148 if (pos[1]<0 && CC_closedAlpha==true)
149 uiPos[1]+=numLines[1]-2;
150
151 return GetNodeWidth(ny, uiPos, dualMesh);
152 }
153
GetNodeArea(int ny,const unsigned int pos[3],bool dualMesh) const154 double Operator_Cylinder::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) const
155 {
156 if (pos[ny]>=numLines[ny]) return 0.0;
157 if (pos[0]>=numLines[0]) return 0.0;
158 if (ny==2)
159 {
160 double da = Operator_Multithread::GetEdgeLength(1,pos,dualMesh)/gridDelta;
161 double r1,r2;
162
163 if (dualMesh)
164 {
165 r1 = GetDiscLine(0,pos[0],false)*gridDelta;
166 r2 = r1 + GetEdgeLength(0,pos,false);
167 }
168 else
169 {
170 r2 = GetDiscLine(0,pos[0],!dualMesh)*gridDelta;
171 r1 = r2 - GetEdgeLength(0,pos,true);
172 }
173
174 if (r1<=0)
175 return da/2 * pow(r2,2);
176 else
177 return da/2* (pow(r2,2) - pow(r1,2));
178 }
179
180 return Operator_Multithread::GetNodeArea(ny,pos,dualMesh);
181 }
182
GetNodeArea(int ny,const int pos[3],bool dualMesh) const183 double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
184 {
185 if ( (pos[0]<0) || (pos[1]<0 && CC_closedAlpha==false) || (pos[2]<0) )
186 return 0.0;
187
188 unsigned int uiPos[]={(unsigned int)pos[0],(unsigned int)pos[1],(unsigned int)pos[2]};
189 if (pos[1]<0 && CC_closedAlpha==true)
190 uiPos[1]+=numLines[1]-2;
191
192 return GetNodeArea(ny, uiPos, dualMesh);
193 }
194
GetEdgeLength(int ny,const unsigned int pos[3],bool dualMesh) const195 double Operator_Cylinder::GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh) const
196 {
197 double length = Operator_Multithread::GetEdgeLength(ny,pos,dualMesh);
198 if (ny!=1)
199 return length;
200 return length * GetDiscLine(0,pos[0],dualMesh);
201 }
202
GetCellVolume(const unsigned int pos[3],bool dualMesh) const203 double Operator_Cylinder::GetCellVolume(const unsigned int pos[3], bool dualMesh) const
204 {
205 return GetEdgeArea(2,pos,dualMesh)*GetEdgeLength(2,pos,dualMesh);
206 }
207
GetEdgeArea(int ny,const unsigned int pos[3],bool dualMesh) const208 double Operator_Cylinder::GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh) const
209 {
210 if (ny!=0)
211 return GetNodeArea(ny,pos,dualMesh);
212
213 return GetEdgeLength(1,pos,!dualMesh) * GetEdgeLength(2,pos,!dualMesh);
214 }
215
FitToAlphaRange(double a_coord,bool fullMesh) const216 double Operator_Cylinder::FitToAlphaRange(double a_coord, bool fullMesh) const
217 {
218 double min = GetDiscLine(1,0);
219 double max = GetDiscLine(1,GetNumberOfLines(1, fullMesh)-1);
220 if ((a_coord>=min) && (a_coord<=max))
221 return a_coord;
222 while (a_coord<min)
223 {
224 a_coord+=2*PI;
225 if (a_coord>max)
226 return a_coord-2*PI;
227 if (a_coord>min)
228 return a_coord;
229 }
230 while (a_coord>max)
231 {
232 a_coord-=2*PI;
233 if (a_coord<min)
234 return a_coord+2*PI;
235 if (a_coord<max)
236 return a_coord;
237 }
238 // this cannot happen
239 return a_coord;
240 }
241
MapAlphaIndex2Range(int pos) const242 int Operator_Cylinder::MapAlphaIndex2Range(int pos) const
243 {
244 if (!CC_closedAlpha)
245 return pos;
246 if (pos<0)
247 return (int)numLines[1]+pos-2;
248 else if (pos>=(int)numLines[1]-2)
249 return pos-(int)numLines[1]+2;
250 else
251 return pos;
252 }
253
GetCellCenterMaterialAvgCoord(const int pos[3],double coord[3]) const254 bool Operator_Cylinder::GetCellCenterMaterialAvgCoord(const int pos[3], double coord[3]) const
255 {
256 if (!CC_closedAlpha || ((pos[1]>=0) && (pos[1]<(int)numLines[1]-2)))
257 {
258 return Operator_Multithread::GetCellCenterMaterialAvgCoord(pos, coord);
259 }
260 if ((pos[0]<0) || (pos[2]<0))
261 return false;
262 int l_pos[3] = {pos[0], 0, pos[2]};
263 l_pos[1] = MapAlphaIndex2Range(pos[1]);
264 return Operator_Multithread::GetCellCenterMaterialAvgCoord(l_pos, coord);
265 }
266
SnapToMeshLine(int ny,double coord,bool & inside,bool dualMesh,bool fullMesh) const267 unsigned int Operator_Cylinder::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh, bool fullMesh) const
268 {
269 if (ny==1)
270 coord=FitToAlphaRange(coord);
271 return Operator_Multithread::SnapToMeshLine(ny, coord, inside, dualMesh, fullMesh);
272 }
273
SnapBox2Mesh(const double * start,const double * stop,unsigned int * uiStart,unsigned int * uiStop,bool dualMesh,bool fullMesh,int SnapMethod,bool * bStartIn,bool * bStopIn) const274 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
275 {
276 double a_min = GetDiscLine(1,0);
277 double a_max = GetDiscLine(1,GetNumberOfLines(1,fullMesh)-1);
278
279 double a_size = stop[1] - start[1];
280 double a_center = FitToAlphaRange(0.5*(stop[1]+start[1]));
281 double a_start = a_center-a_size/2;
282 double a_stop = a_start + a_size;
283 if (a_stop>a_max)
284 a_stop=a_max;
285 if (a_stop<a_min)
286 a_stop=a_min;
287 if (a_start>a_max)
288 a_start=a_max;
289 if (a_start<a_min)
290 a_start=a_min;
291
292 double l_start[3] = {start[0], a_start, start[2]};
293 double l_stop[3] = {stop[0] , a_stop , stop[2] };
294 return Operator_Multithread::SnapBox2Mesh(l_start, l_stop, uiStart, uiStop, dualMesh, fullMesh, SnapMethod, bStartIn, bStopIn);
295 }
296
SnapLine2Mesh(const double * start,const double * stop,unsigned int * uiStart,unsigned int * uiStop,bool dualMesh,bool fullMesh) const297 int Operator_Cylinder::SnapLine2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, bool fullMesh) const
298 {
299 int ret = Operator_Multithread::SnapLine2Mesh(start, stop, uiStart, uiStop, dualMesh, fullMesh);
300
301 if ((stop[1]>start[1]) && (uiStop[1]<uiStart[1]) && (uiStop[1]==0))
302 uiStop[1] = GetNumberOfLines(1, fullMesh)-1-(int)CC_closedAlpha;
303 if ((stop[1]<start[1]) && (uiStop[1]>uiStart[1]) && (uiStop[1]==GetNumberOfLines(1, fullMesh)-1-(int)CC_closedAlpha))
304 uiStop[1] = 0;
305
306 return ret;
307 }
308
309
FindPath(double start[],double stop[])310 Grid_Path Operator_Cylinder::FindPath(double start[], double stop[])
311 {
312 double l_start[3];
313 double l_stop[3];
314
315 for (int n=0;n<3;++n)
316 {
317 l_start[n] = start[n];
318 l_stop[n] = stop[n];
319 }
320
321 while (fabs(l_stop[1]-l_start[1])>PI)
322 {
323 if (l_stop[1]>l_start[1])
324 l_stop[1]-=2*PI;
325 else
326 l_stop[1]+=2*PI;
327 }
328
329 double help=0;
330 if (l_start[1]>l_stop[1])
331 {
332 for (int n=0;n<3;++n)
333 {
334 help = l_start[n];
335 l_start[n] = l_stop[n];
336 l_stop[n] = help;
337 }
338 }
339
340 double a_start = FitToAlphaRange(l_start[1]);
341 double a_stop = FitToAlphaRange(l_stop[1]);
342
343 if (a_stop >= a_start)
344 {
345 l_start[1] = a_start;
346 l_stop[1] = a_stop;
347 return Operator_Multithread::FindPath(l_start, l_stop);
348 }
349
350 // if a-stop fitted to disc range is now smaller than a-start, it must step over the a-bounds...
351
352 Grid_Path path;
353 for (int n=0;n<3;++n)
354 {
355 if ((l_start[n]<GetDiscLine(n,0)) && (l_stop[n]<GetDiscLine(n,0)))
356 return path; //lower bound violation
357 if ((l_start[n]>GetDiscLine(n,GetNumberOfLines(n,true)-1)) && (l_stop[n]>GetDiscLine(n,GetNumberOfLines(n,true)-1)))
358 return path; //upper bound violation
359 }
360
361 if (g_settings.GetVerboseLevel()>2)
362 cerr << __func__ << ": A path was leaving the alpha-direction mesh..." << endl;
363
364 // this section comes into play, if the line moves over the angulare mesh-end/start
365 // we try to have one part of the path on both "ends" of the mesh and stitch them together
366
367 Grid_Path path1;
368 Grid_Path path2;
369
370 // calculate the intersection of the line with the a-max boundary
371 double p0[3],p1[3],p2[3];
372 for (int n=0;n<3;++n)
373 {
374 p0[n] = GetDiscLine(n,0);
375 p1[n] = p0[n];
376 p2[n] = p0[n];
377 }
378 p0[1] = GetDiscLine(1,GetNumberOfLines(1,true)-1-(int)CC_closedAlpha);
379 p1[1] = p0[1];
380 p2[1] = p0[1];
381 p1[0] = discLines[0][numLines[0]-1];
382 p2[2] = discLines[2][numLines[2]-1];
383
384 TransformCoordSystem(p0,p0,m_MeshType,CARTESIAN);
385 TransformCoordSystem(p1,p1,m_MeshType,CARTESIAN);
386 TransformCoordSystem(p2,p2,m_MeshType,CARTESIAN);
387
388 double c_start[3],c_stop[3];
389 TransformCoordSystem(l_start,c_start,m_MeshType,CARTESIAN);
390 TransformCoordSystem(l_stop,c_stop,m_MeshType,CARTESIAN);
391 double intersect[3];
392 double dist;
393 int ret = LinePlaneIntersection(p0,p1,p2,c_start,c_stop,intersect,dist);
394 if (ret<0)
395 {
396 cerr << __func__ << ": Error, unable to calculate intersection, this should not happen!" << endl;
397 return path; // return empty path;
398 }
399
400 if (ret==0)
401 {
402 TransformCoordSystem(intersect,intersect,CARTESIAN,m_MeshType);
403 intersect[1] = GetDiscLine(1,GetNumberOfLines(1,true)-1-(int)CC_closedAlpha);
404 l_start[1] = FitToAlphaRange(l_start[1]);
405 path1 = Operator::FindPath(l_start, intersect);
406 if (g_settings.GetVerboseLevel()>2)
407 cerr << __func__ << ": Intersection top: " << intersect[0] << "," << intersect[1] << "," << intersect[2] << endl;
408 } //otherwise the path was not intersecting the upper a-bound...
409
410 if (CC_closedAlpha==false)
411 {
412 for (int n=0;n<3;++n)
413 {
414 p0[n] = GetDiscLine(n,0);
415 p1[n] = p0[n];
416 p2[n] = p0[n];
417 }
418 p1[0] = discLines[0][numLines[0]-1];
419 p2[2] = discLines[2][numLines[2]-1];
420
421 TransformCoordSystem(p0,p0,m_MeshType,CARTESIAN);
422 TransformCoordSystem(p1,p1,m_MeshType,CARTESIAN);
423 TransformCoordSystem(p2,p2,m_MeshType,CARTESIAN);
424
425 TransformCoordSystem(l_start,c_start,m_MeshType,CARTESIAN);
426 TransformCoordSystem(l_stop,c_stop,m_MeshType,CARTESIAN);
427
428 ret = LinePlaneIntersection(p0,p1,p2,c_start,c_stop,intersect,dist);
429 TransformCoordSystem(intersect,intersect,CARTESIAN,m_MeshType);
430 }
431
432 if (ret==0)
433 {
434 intersect[1] = GetDiscLine(1,0);
435 l_stop[1] = FitToAlphaRange(l_stop[1]);
436 path2 = Operator::FindPath(intersect, l_stop);
437 if (g_settings.GetVerboseLevel()>2)
438 cerr << __func__ << ": Intersection bottom: " << intersect[0] << "," << intersect[1] << "," << intersect[2] << endl;
439 }
440
441 //combine path
442 for (size_t t=0; t<path1.dir.size(); ++t)
443 {
444 path.posPath[0].push_back(path1.posPath[0].at(t));
445 path.posPath[1].push_back(path1.posPath[1].at(t));
446 path.posPath[2].push_back(path1.posPath[2].at(t));
447 path.dir.push_back(path1.dir.at(t));
448 }
449 for (size_t t=0; t<path2.dir.size(); ++t)
450 {
451 path.posPath[0].push_back(path2.posPath[0].at(t));
452 path.posPath[1].push_back(path2.posPath[1].at(t));
453 path.posPath[2].push_back(path2.posPath[2].at(t));
454 path.dir.push_back(path2.dir.at(t));
455 }
456
457 if (CC_closedAlpha==true)
458 for (size_t t=0; t<path.dir.size(); ++t)
459 {
460 if ( ((path.dir.at(t)==0) || (path.dir.at(t)==2)) && (path.posPath[1].at(t)==0))
461 path.posPath[1].at(t) = numLines[1]-2;
462 }
463
464 return path;
465 }
466
SetupCSXGrid(CSRectGrid * grid)467 bool Operator_Cylinder::SetupCSXGrid(CSRectGrid* grid)
468 {
469 unsigned int alphaNum;
470 double* alphaLines = NULL;
471 alphaLines = grid->GetLines(1,alphaLines,alphaNum,true);
472
473 double minmaxA = fabs(alphaLines[alphaNum-1]-alphaLines[0]);
474 if (fabs(minmaxA-2*PI) < OPERATOR_CYLINDER_CLOSED_ALPHA_THRESHOLD)
475 {
476 if (g_settings.GetVerboseLevel()>0)
477 cout << "Operator_Cylinder::SetupCSXGrid: Alpha is a full 2*PI => closed Cylinder..." << endl;
478 CC_closedAlpha = true;
479 grid->SetLine(1,alphaNum-1,2*PI+alphaLines[0]);
480 grid->AddDiscLine(1,2*PI+alphaLines[1]);
481 }
482 else if (minmaxA>2*PI)
483 {
484 cerr << "Operator_Cylinder::SetupCSXGrid: Alpha Max-Min must not be larger than 2*PI!!!" << endl;
485 Reset();
486 return false;
487 }
488 else
489 {
490 CC_closedAlpha=false;
491 }
492
493 CC_R0_included = false;
494 if (grid->GetLine(0,0)<0)
495 {
496 cerr << "Operator_Cylinder::SetupCSXGrid: r<0 not allowed in Cylinder Coordinates!!!" << endl;
497 Reset();
498 return false;
499 }
500 else if (grid->GetLine(0,0)==0.0)
501 {
502 if (g_settings.GetVerboseLevel()>0)
503 cout << "Operator_Cylinder::SetupCSXGrid: r=0 included..." << endl;
504 CC_R0_included = CC_closedAlpha; //needed for correct ec-calculation, deactivate if closed cylinder is false... --> E_r = 0 anyways
505 }
506
507 #ifdef MPI_SUPPORT
508 // Setup an MPI split in alpha direction for a closed cylinder
509 CC_MPI_Alpha = false;
510 if ((m_NeighborUp[1]>=0) || (m_NeighborDown[1]>=0)) //check for MPI split in alpha direction
511 {
512 double minmaxA = 2*PI;// fabs(m_OrigDiscLines[1][m_OrigNumLines[1]-1]-m_OrigDiscLines[1][0]);
513 if (fabs(minmaxA-2*PI) < OPERATOR_CYLINDER_CLOSED_ALPHA_THRESHOLD) //check for closed alpha MPI split
514 {
515 CC_MPI_Alpha = true;
516 if (m_OrigDiscLines[0][0]==0)
517 {
518 cerr << "Operator_Cylinder::SetupCSXGrid: Error: MPI split in alpha direction for closed cylinder including r==0 is currently not supported! Exit!" << endl;
519 exit(-2);
520 }
521
522 if (m_NeighborUp[1]<0) //check if this process is at the alpha-end
523 {
524 grid->SetLine(1,alphaNum-1,2*PI+m_OrigDiscLines[1][0]);
525 grid->AddDiscLine(1,2*PI+m_OrigDiscLines[1][1]);
526
527 SetNeighborUp(1,m_ProcTable[m_ProcTablePos[0]][0][m_ProcTablePos[2]]);
528 }
529
530 if (m_NeighborDown[1]<0) //check if this process is at the alpha-start
531 {
532 SetNeighborDown(1,m_ProcTable[m_ProcTablePos[0]][m_SplitNumber[1]-1][m_ProcTablePos[2]]);
533 }
534
535 //Note: the process table will not reflect this up/down neighbors necessary for a closed cylinder
536 }
537 }
538 #endif
539
540 if (Operator_Multithread::SetupCSXGrid(grid)==false)
541 return false;
542
543 if (CC_closedAlpha || CC_R0_included)
544 {
545 m_Cyl_Ext = new Operator_Ext_Cylinder(this);
546 this->AddExtension(m_Cyl_Ext);
547 }
548
549 return true;
550 }
551
ApplyMagneticBC(bool * dirs)552 void Operator_Cylinder::ApplyMagneticBC(bool* dirs)
553 {
554 if (dirs==NULL) return;
555 if (CC_closedAlpha)
556 {
557 dirs[2]=0;
558 dirs[3]=0; //no PMC in alpha directions...
559 }
560 if (CC_R0_included)
561 {
562 dirs[0]=0; //no PMC in r_min directions...
563 }
564 Operator_Multithread::ApplyMagneticBC(dirs);
565 }
566
AddExtension(Operator_Extension * op_ext)567 void Operator_Cylinder::AddExtension(Operator_Extension* op_ext)
568 {
569 if (op_ext->IsCylinderCoordsSave(CC_closedAlpha, CC_R0_included))
570 Operator_Multithread::AddExtension(op_ext);
571 else
572 {
573 cerr << "Operator_Cylinder::AddExtension: Warning: Operator extension \"" << op_ext->GetExtensionName() << "\" is not compatible with cylinder-coords!! skipping...!" << endl;
574 delete op_ext;
575 }
576 }
577