fixed writing out entries in advective bc
[OpenFOAM-1.6-ext.git] / src / engine / engineTopoChangerMesh / thoboisSliding / addPistonFacesPointZonesThoboisSliding.H
blobaf773c8a18973009f930d0d1a0ceb5f0c6994c90
1     // Add the piston zone
2     if (piston().patchID().active())
3     {
5         // Piston position
6         
7         Info << "Adding face zone for piston layer addition/removal" << endl;
8         
9         label pistonPatchID = piston().patchID().index();
10         
11         scalarField zPistonValves(valves().size() + 1,max(boundary()[pistonPatchID].patch().localPoints()).z());
12         
13         scalarField zPistonValvesOffset = zPistonValves + offSet();
14         
15         labelListList zonePiston(valves().size()+1);
16         boolListList flipZonePiston(valves().size() + 1);
17         labelList nZoneFacesPiston(valves().size() + 1, 0);
18         boolList foundOneFace(valves().size() + 1, false);        
19         scalarList zHigher(valves().size() + 1, GREAT);
20         scalarList zLower(valves().size() + 1, GREAT);
21         scalarList dh(valves().size() + 1, GREAT);
22         scalarList dl(valves().size() + 1, GREAT);
23                 
24         forAll(zonePiston, zoneI)
25         {
26             zonePiston[zoneI].setSize(faceCentres().size());
27             flipZonePiston[zoneI].setSize(faceCentres().size());
28             flipZonePiston[zoneI] = false;
29         }
30         
31         forAll (faceCentres(), faceI)
32         {
33             bool foundInValve = false;
34         
35             scalar zc = faceCentres()[faceI].z();        
36             
37             vector n = faceAreas()[faceI]/mag(faceAreas()[faceI]);
38             scalar dd = n & vector(0,0,1);
39                         
40             forAll(valves(), valveI)
41             {
43                 if(inValve(faceCentres()[faceI], valveI))
44                 {
45                     foundInValve = true;
46                     
47                     scalar zPistV = zPistonValvesOffset[valveI];
48                     
49                     if (mag(dd) > 0.1)
50                     {
51                         if (zPistV - zc > 0 && zPistV - zc < dl[valveI])
52                         {
53                             zLower[valveI] = zc;
54                             dl[valveI] = zPistV - zc;
55                         }
56             
57                         if (zc - zPistV > 0 && zc - zPistV < dh[valveI])
58                         {
59                             zHigher[valveI] = zc;
60                             dh[valveI] = zc - zHigher[valveI];
61                         }
62             
63                         if
64                         (
65                             zc > zPistV - delta()
66                             && zc < zPistV + delta()
67                         )   
68                         {
70                             foundOneFace[valveI] = true;
72                             if ((faceAreas()[faceI] & vector(0,0,1)) < 0)
73                             {
74                                 flipZonePiston[valveI][nZoneFacesPiston[valveI]] = true;
75                             }
77                             zonePiston[valveI][nZoneFacesPiston[valveI]] = faceI;
78                             nZoneFacesPiston[valveI]++;
80                         }
81                     }
82                 }
83             }
85             if(!foundInValve)
86             {
87                 label valveSize = valves().size();
88             
89                 scalar zPistV = zPistonValvesOffset[valveSize];
90                     
91                 if (mag(dd) > 0.1)
92                 {
93                     if (zPistV - zc > 0 && zPistV - zc < dl[valveSize])
94                     {
95                         zLower[valveSize] = zc;
96                         dl[valveSize] = zPistV - zc;
97                     }
98             
99                     if (zc - zPistV > 0 && zc - zPistV < dh[valveSize])
100                     {
101                         zHigher[valveSize] = zc;
102                         dh[valveSize] = zc - zHigher[valveSize];
103                     }
104             
105                     if
106                     (
107                         zc > zPistV - delta()
108                         && zc < zPistV + delta()
109                     )   
110                     {
112                         foundOneFace[valveSize] = true;
114                         if ((faceAreas()[faceI] & vector(0,0,1)) < 0)
115                         {
116                             flipZonePiston[valveSize][nZoneFacesPiston[valveSize]] = true;
117                         }
119                         zonePiston[valveSize][nZoneFacesPiston[valveSize]] = faceI;
120                         nZoneFacesPiston[valveSize]++;
122                     }
123                 }
124             }
125         }
126         
128         // if no cut was found use the layer above
130         forAll(valves(), valveI)
131         {
132         
133             if (!foundOneFace[valveI])
134             {
135                 zPistonValvesOffset[valveI] = zHigher[valveI];
136             }
137         }
138         
139         if(!foundOneFace[valves().size()])
140         {
141             zPistonValvesOffset[valves().size()] = zHigher[valves().size()];
142         }
144         forAll (faceCentres(), faceI)
145         {
147             bool foundInValve = false;
149             scalar zc = faceCentres()[faceI].z();
151             vector n = faceAreas()[faceI]/mag(faceAreas()[faceI]);
152             scalar dd = n & vector(0,0,1);
153                                     
154             forAll(valves(), valveI)
155             {
157                 if(!foundOneFace[valveI])
158                 {
159                     
160                     scalar zPistV = zPistonValvesOffset[valveI];    
161                     
162                         
163                     if(inValve(faceCentres()[faceI], valveI))
164                     {
165                         foundInValve = true;
166                         
167                         if (mag(dd) > 0.1)
168                         {   
170                             if
171                             (
172                                 zc > zPistV - delta()
173                                 && zc < zPistV + delta()
174                             )
175                             {
176                                 
177                                 if ((faceAreas()[faceI] & vector(0,0,1)) < 0)
178                                 {
179                                     flipZonePiston[valveI][nZoneFacesPiston[valveI]] = true;
180                                 }
181             
182                                 zonePiston[valveI][nZoneFacesPiston[valveI]] = faceI;
183                                 nZoneFacesPiston[valveI]++;
184                             }
185                         }
186                     }
188                 }
189             }
190             
191             if(!foundInValve && !foundOneFace[valves().size()])
192             {
193                 label valveSize = valves().size();
194                 scalar zPistV = zPistonValvesOffset[valveSize];    
195                 
196                 if (mag(dd) > 0.1)
197                 {   
199                     if
200                     (
201                         zc > zPistV - delta()
202                         && zc < zPistV + delta()
203                     )
204                     {
205                         if ((faceAreas()[faceI] & vector(0,0,1)) < 0)
206                         {
207                             flipZonePiston[valveSize][nZoneFacesPiston[valveSize]] = true;
208                         }
209             
210                         zonePiston[valveSize][nZoneFacesPiston[valveSize]] = faceI;
211                         nZoneFacesPiston[valveSize]++;
212                     }
213                 }
214             }
215         }
217         forAll(valves(), valveI)
218         {
219             zonePiston[valveI].setSize(nZoneFacesPiston[valveI]);
220             flipZonePiston[valveI].setSize(nZoneFacesPiston[valveI]);
222             fz.append
223             (
224                 new faceZone
225                 (
226                     "pistonLayerFacesV" + Foam::name(valveI + 1),
227                     zonePiston[valveI],
228                     flipZonePiston[valveI],
229                     nFaceZones,
230                     faceZones()
231                 )    
232             );
233         
234             nFaceZones++;
235         }
236         
237         {
238             label valveSize = valves().size();
239             zonePiston[valveSize].setSize(nZoneFacesPiston[valveSize]);
240             flipZonePiston[valveSize].setSize(nZoneFacesPiston[valveSize]);
242             fz.append
243             (
244                 new faceZone
245                 (
246                     "pistonLayerFaces",
247                     zonePiston[valveSize],
248                     flipZonePiston[valveSize],
249                     nFaceZones,
250                     faceZones()
251                 )    
252             );
253         
254             nFaceZones++;
255         
256         }
257         
258         forAll(valves(), valveI)        
259         {
260         
261             labelList movingCells(nCells());
262                             
263             label nMovingCells = 0;
264                          
265             forAll(cellCentres(),cellI)
266             {
267                 const vector& v = cellCentres()[cellI];
268                     
269                 if
270                 (
271                     inValve(v, valveI)
272                     &&
273                     v.z() < zPistonValvesOffset[valveI] 
274                 )
275                 {
276                     movingCells[nMovingCells] = cellI;
277                     nMovingCells++;
278                 }
279         
280             }
281             
282             movingCells.setSize(nMovingCells);
283             Info << "Number of cells in the moving region poppet valve: " << nMovingCells << endl;
284                                             
285             cz.append
286             (
287                 new cellZone
288                 (
289                     "movingCellsPistonV"+ Foam::name(valveI + 1),
290                     movingCells,
291                     nCellZones,
292                     cellZones()
293                 )
294             );
296             nCellZones++;
298             // valve top points (move all with valve displacement)
299             DynamicList<label> valvePistonPoints(nPoints() / 10);
300             List<bool> valvePistonPoint(nPoints(), false);
301             bool foundOne = false;
302             
303             const cellList& c = cells();
304             const faceList& f = allFaces();
306             forAll (movingCells, cellI)
307             {
308                 const cell& curCell = c[movingCells[cellI]];
310                 forAll (curCell, faceI)
311                 {
312                     // Mark all the points as moving
313                     const face& curFace = f[curCell[faceI]];
315                     forAll (curFace, pointI)
316                     {
317                         foundOne = true;  
318                         valvePistonPoint[curFace[pointI]] = true;            
319                     }
320         
321                 }
322         
323             }
324             
325             forAll(valvePistonPoint, pointI)
326             {
327                 if(valvePistonPoint[pointI])
328                 {
329                     valvePistonPoints.append(pointI);
330                 }
331             }
332                                        
333             pz.append
334             (
335                 new pointZone
336                 (
337                     "valvePistonPointsV"+ Foam::name(valveI + 1),
338                     valvePistonPoints.shrink(),
339                     nPointZones,
340                     pointZones()
341                 )
342             );
343             
344             nPointZones++;
346         
347         }
348         
349         
350         {
351             label valveSize = valves().size();
352             labelList movingCells(nCells());
353                             
354             label nMovingCells = 0;
355                          
356             forAll(cellCentres(),cellI)
357             {
358                 const vector& v = cellCentres()[cellI];
359                 
360                 bool fallInValve = false;
361                 forAll(valves(), valveI)
362                 {
363                     if(inValve(v, valveI))
364                     {
365                         fallInValve = true;
366                     }
367                 }
368                     
369                 if
370                 (
371                     !fallInValve
372                     &&
373                     v.z() < zPistonValvesOffset[valveSize] 
374                 )
375                 {
376                     movingCells[nMovingCells] = cellI;
377                     nMovingCells++;
378                 }
379         
380             }
381             
382             movingCells.setSize(nMovingCells);
383             Info << "Number of cells in the moving region poppet valve: " << nMovingCells << endl;
384                                             
385             cz.append
386             (
387                 new cellZone
388                 (
389                     "movingCellsPiston",
390                     movingCells,
391                     nCellZones,
392                     cellZones()
393                 )
394             );
396             nCellZones++;
398             // valve top points (move all with valve displacement)
399             DynamicList<label> pistonPoints(nPoints() / 10);
400             List<bool> pistonPoint(nPoints(), false);
401             bool foundOne = false;
402             
403             const cellList& c = cells();
404             const faceList& f = allFaces();
406             forAll (movingCells, cellI)
407             {
408                 const cell& curCell = c[movingCells[cellI]];
410                 forAll (curCell, faceI)
411                 {
412                     // Mark all the points as moving
413                     const face& curFace = f[curCell[faceI]];
415                     forAll (curFace, pointI)
416                     {                            
417                         foundOne = true;  
418                         pistonPoint[curFace[pointI]] = true;            
419                     }
420         
421                 }
422         
423             }
424             
425             forAll(pistonPoint, pointI)
426             {
427                 if(pistonPoint[pointI])
428                 {
429                     pistonPoints.append(pointI);
430                 }
431             }
432                                         
433             pz.append
434             (
435                 new pointZone
436                 (
437                     "pistonPoints",
438                     pistonPoints.shrink(),
439                     nPointZones,
440                     pointZones()
441                 )
442             );
443             
444             nPointZones++;
445         }
446         
447     }