UoN / SACAMOS

Blame view

SRC/WRITE_FH2_IPFILE/set_segments_from_grid.F90 11.5 KB
Edit Raw Normal View History
ca8ab9e9   Chris Smartt   Update to proximi... 
1
2
3
4
5
6
7
8
9
10
11
12
13
 

    nxmin=conductor_data(conductor)%nxmin
    nxmax=conductor_data(conductor)%nxmax
    nymin=conductor_data(conductor)%nymin
    nymax=conductor_data(conductor)%nymax
    
    layer_number=0
      
    if (conductor_data(conductor)%mesh_to_layer_type.EQ.1) then
! each marked cell in the grid becomes a layer and hence a segment
    
      do ix=nxmin,nxmax
        do iy=nymin,nymax
197629dd   Chris Smartt   Update proximity_... 
14
15
 
          
          if ( conductor_data(conductor)%grid(ix,iy).EQ.1 ) then
ca8ab9e9   Chris Smartt   Update to proximi... 
16
197629dd   Chris Smartt   Update proximity_... 
17
 
            layer_number=layer_number+1
ca8ab9e9   Chris Smartt   Update to proximi... 
18
197629dd   Chris Smartt   Update proximity_... 
19
20
 
            apx=conductor_data(conductor)%xc+real(ix)*conductor_data(conductor)%dl
            apy=conductor_data(conductor)%yc+real(iy)*conductor_data(conductor)%dl
ca8ab9e9   Chris Smartt   Update to proximi... 
21
22
23
24
25
26
27
28
29
30
 
    
            conductor_data(conductor)%x(layer_number)=apx
            conductor_data(conductor)%y(layer_number)=apy
            conductor_data(conductor)%w(layer_number)=conductor_data(conductor)%dl
            conductor_data(conductor)%h(layer_number)=conductor_data(conductor)%dl
              
            conductor_data(conductor)%d(layer_number)=0.0
            conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
            conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
197629dd   Chris Smartt   Update proximity_... 
31
 
          end if
ca8ab9e9   Chris Smartt   Update to proximi... 
32
33
34
35
36
37
 
        end do
      end do

    else if (conductor_data(conductor)%mesh_to_layer_type.EQ.2) then
    
! each row of cells in the grid becomes a segment
197629dd   Chris Smartt   Update proximity_... 
38
 
      do iy=nymin,nymax
ca8ab9e9   Chris Smartt   Update to proximi... 
39
197629dd   Chris Smartt   Update proximity_... 
40
41
42
 
        in_conductor=.FALSE.
        cx1=0
        cx2=0
ca8ab9e9   Chris Smartt   Update to proximi... 
43
197629dd   Chris Smartt   Update proximity_... 
44
 
        do ix=nxmin,nxmax
ca8ab9e9   Chris Smartt   Update to proximi... 
45
197629dd   Chris Smartt   Update proximity_... 
46
 
          if ( (.NOT.in_conductor).AND.(conductor_data(conductor)%grid(ix,iy).EQ.1) ) then
ca8ab9e9   Chris Smartt   Update to proximi... 
47
48
 
    
! we have moved from air to conductor so save this x cell number
197629dd   Chris Smartt   Update proximity_... 
49
50
 
             cx1=ix
             in_conductor=.TRUE.
ca8ab9e9   Chris Smartt   Update to proximi... 
51
197629dd   Chris Smartt   Update proximity_... 
52
53
 
           else if ( (in_conductor).AND.(conductor_data(conductor)%grid(ix,iy).EQ.0) ) then
! we have moved from conductor to air so create a layer for this region
ca8ab9e9   Chris Smartt   Update to proximi... 
54
197629dd   Chris Smartt   Update proximity_... 
55
56
 
             cx2=ix-1
             in_conductor=.FALSE.
ca8ab9e9   Chris Smartt   Update to proximi... 
57
197629dd   Chris Smartt   Update proximity_... 
58
59
60
61
 
             layer_number=layer_number+1
         
             cw=real(cx2-cx1+1)*conductor_data(conductor)%dl
             ch=conductor_data(conductor)%dl
ca8ab9e9   Chris Smartt   Update to proximi... 
62
197629dd   Chris Smartt   Update proximity_... 
63
64
 
             apx=conductor_data(conductor)%xc+real(cx2+cx1)/2.0
             apy=conductor_data(conductor)%yc+real(iy)*conductor_data(conductor)%dl
ca8ab9e9   Chris Smartt   Update to proximi... 
65
197629dd   Chris Smartt   Update proximity_... 
66
67
68
69
 
             conductor_data(conductor)%x(layer_number)=apx
             conductor_data(conductor)%y(layer_number)=apy
             conductor_data(conductor)%w(layer_number)=cw
             conductor_data(conductor)%h(layer_number)=ch
ca8ab9e9   Chris Smartt   Update to proximi... 
70
197629dd   Chris Smartt   Update proximity_... 
71
72
73
 
             conductor_data(conductor)%d(layer_number)=0.0
             conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
             conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
ca8ab9e9   Chris Smartt   Update to proximi... 
74
197629dd   Chris Smartt   Update proximity_... 
75
76
77
78
 
             cx1=0
             cx2=0
             
          end if
ca8ab9e9   Chris Smartt   Update to proximi... 
79
197629dd   Chris Smartt   Update proximity_... 
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
 
        end do ! next col

      end do ! next row   

    else if (conductor_data(conductor)%mesh_to_layer_type.EQ.3) then

! specific process for cylinders...    
! First find the central square region
      found_square=.FALSE.
      do ix=0,nxmax-2
        iy=ix
                  
        if ( (conductor_data(conductor)%grid(ix,iy).EQ.1).AND.  &
             (conductor_data(conductor)%grid(ix+1,iy+1).EQ.0)) then                                            
           
          layer_number=layer_number+1
          cx1=ix
          
          cw=2d0*(real(ix)+0.5d0)*conductor_data(conductor)%dl
          ch=2d0*(real(iy)+0.5d0)*conductor_data(conductor)%dl
       
          apx=conductor_data(conductor)%xc     ! centre of circle
          apy=conductor_data(conductor)%yc
    
          conductor_data(conductor)%x(layer_number)=apx
          conductor_data(conductor)%y(layer_number)=apy
          conductor_data(conductor)%w(layer_number)=cw
          conductor_data(conductor)%h(layer_number)=ch
            
          conductor_data(conductor)%d(layer_number)=0.0
          
          if (conductor_data(conductor)%auto_grid_density) then
            CALL calc_nmin(conductor_data(conductor)%skin_depth,conductor_data(conductor)%w(layer_number), &
                                          conductor_data(conductor)%rw,nw_auto)
            CALL calc_nmin(conductor_data(conductor)%skin_depth,conductor_data(conductor)%h(layer_number), &
                                          conductor_data(conductor)%rh,nh_auto)
            conductor_data(conductor)%anwinc(layer_number)=nw_auto
            conductor_data(conductor)%anhinc(layer_number)=nh_auto         
          else
            conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
            conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
          end if
          
          found_square=.TRUE.
          EXIT                   ! leave the loop
    
        end if   ! found the corner
ca8ab9e9   Chris Smartt   Update to proximi... 
127
197629dd   Chris Smartt   Update proximity_... 
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
 
      end do ! next row
      
      if (.NOT.found_square) then
        write(*,*)'ERROR: not found central square...'
        STOP 1
      end if
      
! now search out in the x directiion adding layers
      do ix=cx1+1,nxmax
        
        do iy=0,cx1  ! search in y to find the extent of this layer
                  
          if ( (conductor_data(conductor)%grid(ix,iy).EQ.1).AND.  &
               (conductor_data(conductor)%grid(ix,iy+1).EQ.0)) then
             
! we have moved from metal to air so add four segments
ca8ab9e9   Chris Smartt   Update to proximi... 
144
197629dd   Chris Smartt   Update proximity_... 
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
 
! add layer on xmax side           
            layer_number=layer_number+1
         
            cw=conductor_data(conductor)%dl
            ch=2d0*real(iy+0.5d0)*conductor_data(conductor)%dl       
            apx=conductor_data(conductor)%xc+real(ix)*conductor_data(conductor)%dl
            apy=conductor_data(conductor)%yc
    
            conductor_data(conductor)%x(layer_number)=apx
            conductor_data(conductor)%y(layer_number)=apy
            conductor_data(conductor)%w(layer_number)=cw
            conductor_data(conductor)%h(layer_number)=ch
            
            conductor_data(conductor)%d(layer_number)=0.0
          
            if (conductor_data(conductor)%auto_grid_density) then
              CALL calc_nmin(conductor_data(conductor)%skin_depth,conductor_data(conductor)%w(layer_number), &
                                          conductor_data(conductor)%rw,nw_auto)
              CALL calc_nmin(conductor_data(conductor)%skin_depth,conductor_data(conductor)%h(layer_number), &
                                          conductor_data(conductor)%rh,nh_auto)
              conductor_data(conductor)%anwinc(layer_number)=nw_auto
              conductor_data(conductor)%anhinc(layer_number)=nh_auto         
            else
              conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
              conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
            end if

! add layer on xmin side           
            layer_number=layer_number+1
         
            cw=conductor_data(conductor)%dl
            ch=2d0*real(iy+0.5d0)*conductor_data(conductor)%dl       
            apx=conductor_data(conductor)%xc-real(ix)*conductor_data(conductor)%dl
            apy=conductor_data(conductor)%yc
    
            conductor_data(conductor)%x(layer_number)=apx
            conductor_data(conductor)%y(layer_number)=apy
            conductor_data(conductor)%w(layer_number)=cw
            conductor_data(conductor)%h(layer_number)=ch
            
            conductor_data(conductor)%d(layer_number)=0.0
            if (conductor_data(conductor)%auto_grid_density) then
              CALL calc_nmin(conductor_data(conductor)%skin_depth,conductor_data(conductor)%w(layer_number), &
                                          conductor_data(conductor)%rw,nw_auto)
              CALL calc_nmin(conductor_data(conductor)%skin_depth,conductor_data(conductor)%h(layer_number), &
                                          conductor_data(conductor)%rh,nh_auto)
              conductor_data(conductor)%anwinc(layer_number)=nw_auto
              conductor_data(conductor)%anhinc(layer_number)=nh_auto         
            else
              conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
              conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
            end if

! add layer on ymax side           
            layer_number=layer_number+1
         
            cw=2d0*real(iy+0.5d0)*conductor_data(conductor)%dl       
            ch=conductor_data(conductor)%dl
            apx=conductor_data(conductor)%xc
            apy=conductor_data(conductor)%yc+real(ix)*conductor_data(conductor)%dl
    
            conductor_data(conductor)%x(layer_number)=apx
            conductor_data(conductor)%y(layer_number)=apy
            conductor_data(conductor)%w(layer_number)=cw
            conductor_data(conductor)%h(layer_number)=ch
            
            conductor_data(conductor)%d(layer_number)=0.0
            if (conductor_data(conductor)%auto_grid_density) then
              CALL calc_nmin(conductor_data(conductor)%skin_depth,conductor_data(conductor)%w(layer_number), &
                                          conductor_data(conductor)%rw,nw_auto)
              CALL calc_nmin(conductor_data(conductor)%skin_depth,conductor_data(conductor)%h(layer_number), &
                                          conductor_data(conductor)%rh,nh_auto)
              conductor_data(conductor)%anwinc(layer_number)=nw_auto
              conductor_data(conductor)%anhinc(layer_number)=nh_auto         
            else
              conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
              conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
            end if

! add layer on ymin side           
            layer_number=layer_number+1
         
            cw=2d0*real(iy+0.5d0)*conductor_data(conductor)%dl       
            ch=conductor_data(conductor)%dl
            apx=conductor_data(conductor)%xc
            apy=conductor_data(conductor)%yc-real(ix)*conductor_data(conductor)%dl
    
            conductor_data(conductor)%x(layer_number)=apx
            conductor_data(conductor)%y(layer_number)=apy
            conductor_data(conductor)%w(layer_number)=cw
            conductor_data(conductor)%h(layer_number)=ch
            
            conductor_data(conductor)%d(layer_number)=0.0
            if (conductor_data(conductor)%auto_grid_density) then
              CALL calc_nmin(conductor_data(conductor)%skin_depth,conductor_data(conductor)%w(layer_number), &
                                          conductor_data(conductor)%rw,nw_auto)
              CALL calc_nmin(conductor_data(conductor)%skin_depth,conductor_data(conductor)%h(layer_number), &
                                          conductor_data(conductor)%rh,nh_auto)
              conductor_data(conductor)%anwinc(layer_number)=nw_auto
              conductor_data(conductor)%anhinc(layer_number)=nh_auto         
            else
              conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
              conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
            end if
    
            EXIT ! this layer has been set so go to the next x layer
    
          end if   ! found the corner
    
        end do ! next y
    
      end do ! next row
ca8ab9e9   Chris Smartt   Update to proximi... 
257
258
259
260
261
262
263
264
265
 

    else
      write(*,*)'Unknown mesh_to_layer_type:',conductor_data(conductor)%mesh_to_layer_type, &
                'in conductor',conductor
      STOP 1
    end if
    
! the number of layers may have changed if we have combined grid cells      
    conductor_data(conductor)%tot_n_layers=layer_number