Line data Source code
1 : module mesh_mod
2 :
3 : use shr_kind_mod, only: r8=>shr_kind_r8
4 : use physconst, only: PI
5 : use control_mod, only: MAX_FILE_LEN
6 : use cam_abortutils, only: endrun
7 :
8 : use netcdf, only: nf90_strerror, nf90_open, nf90_close
9 : use netcdf, only: NF90_NOWRITE, nf90_NoErr
10 : use netcdf, only: nf90_inq_dimid, nf90_inquire_dimension
11 : use netcdf, only: nf90_inq_varid, nf90_get_var
12 :
13 : implicit none
14 : logical, public :: MeshUseMeshFile = .false.
15 :
16 : public :: MeshOpen ! Must be called first
17 :
18 : integer, parameter :: MXSTLN = 32
19 :
20 : ! ===============================
21 : ! Public methods for mesh_mod
22 : ! ===============================
23 :
24 : public :: MeshCubeEdgeCount ! called anytime afer MeshOpen
25 : public :: MeshCubeElemCount ! called anytime afer MeshOpen
26 : public :: MeshCubeTopology ! called afer MeshOpen
27 : public :: MeshSetCoordinates ! called after MeshCubeTopology
28 : public :: MeshPrint ! show the contents of the Mesh after it has been loaded into the module
29 : public :: MeshClose
30 : ! ===============================
31 : ! Private members
32 : ! ===============================
33 :
34 : integer,private,parameter :: nfaces = 6 ! number of faces on the cube
35 : integer,private,parameter :: nInnerElemEdge = 8 ! number of edges for an interior element
36 :
37 : character (len=MAX_FILE_LEN), private :: p_mesh_file_name
38 : integer , private :: p_ncid
39 : integer , private :: p_number_elements
40 : integer , private :: p_number_elements_per_face
41 : integer , private :: p_number_blocks
42 : integer , private :: p_number_nodes
43 : integer , private :: p_number_dimensions
44 : integer , private :: p_number_neighbor_edges
45 : real(kind=r8) , private, allocatable :: p_node_coordinates(:,:)
46 : integer , private, allocatable :: p_connectivity(:,:)
47 :
48 : ! ===============================
49 : ! Private methods
50 : ! ===============================
51 :
52 : private :: create_index_table
53 : private :: find_side_neighbors
54 : private :: find_corner_neighbors
55 : private :: get_node_coordinates
56 : private :: get_2D_sub_coordinate_indexes
57 : private :: mesh_connectivity
58 : private :: cube_face_element_centroids
59 : private :: smallest_diameter_element
60 : private :: cube_to_cube_coordinates
61 : private :: sphere_to_cube_coordinates
62 : private :: initialize_space_filling_curve
63 :
64 : private :: handle_error
65 : private :: open_mesh_file
66 : private :: close_mesh_file
67 : private :: get_number_of_elements
68 : private :: get_number_of_dimensions
69 : private :: get_number_of_elements_per_face
70 : private :: get_number_of_nodes
71 : private :: get_number_of_element_blocks
72 : private :: get_node_multiplicity
73 : private :: get_face_connectivity
74 :
75 : CONTAINS
76 :
77 : !======================================================================
78 : ! subroutine handle_error
79 : !======================================================================
80 0 : subroutine handle_error (status, file, line)
81 :
82 : integer, intent(in) :: status
83 : character (len=*), intent(in) :: file
84 : integer, intent(in) :: line
85 0 : print *, file,':', line, ': ', trim(nf90_strerror(status))
86 0 : call endrun("Terminating program due to netcdf error while obtaining mesh information, please see message in standard output.")
87 0 : end subroutine handle_error
88 :
89 : !======================================================================
90 : ! open_mesh_file()
91 : !
92 : !> Open the netcdf file containing the mesh.
93 : !! Assign the holder to the file to p_ncid so everyone else knows
94 : !! how to use it without passing the argument around.
95 : !======================================================================
96 0 : subroutine open_mesh_file()
97 : implicit none
98 : integer :: status
99 :
100 0 : status = nf90_open(p_mesh_file_name, NF90_NOWRITE, p_ncid)
101 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
102 :
103 0 : MeshUseMeshFile = .true.
104 :
105 0 : end subroutine open_mesh_file
106 :
107 : !======================================================================
108 : ! subroutine close_mesh_file()
109 : !======================================================================
110 :
111 0 : subroutine close_mesh_file()
112 : implicit none
113 : integer :: status
114 :
115 0 : status = nf90_close(p_ncid)
116 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
117 :
118 0 : end subroutine close_mesh_file
119 :
120 : !======================================================================
121 : ! function get_number_of_dimensions()
122 : !======================================================================
123 :
124 0 : function get_number_of_dimensions() result(number_dimensions)
125 : implicit none
126 : integer :: number_dimensions
127 :
128 : ! local variables
129 : integer :: status, number_of_dim_id
130 :
131 : ! Get the id of 'num_elem', if such dimension is not there panic and quit :P
132 0 : status = nf90_inq_dimid(p_ncid, "num_dim", number_of_dim_id)
133 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
134 :
135 : ! How many values for 'num_elem' are there?
136 0 : status = nf90_inquire_dimension(p_ncid, number_of_dim_id, len = number_dimensions)
137 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
138 :
139 0 : end function get_number_of_dimensions
140 :
141 : !======================================================================
142 : ! function get_number_of_elements()
143 : !======================================================================
144 :
145 0 : function get_number_of_elements() result(number_elements)
146 : implicit none
147 : integer :: number_elements
148 : ! local variables
149 : integer :: status, number_of_elements_id
150 :
151 : ! Get the id of 'num_elem', if such dimension is not there panic and quit :P
152 0 : status = nf90_inq_dimid(p_ncid, "num_elem", number_of_elements_id)
153 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
154 :
155 : ! How many values for 'num_elem' are there?
156 0 : status = nf90_inquire_dimension(p_ncid, number_of_elements_id, len = number_elements)
157 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
158 :
159 0 : end function get_number_of_elements
160 :
161 : !======================================================================
162 : ! function get_number_of_nodes()
163 : !======================================================================
164 0 : function get_number_of_nodes() result(number_nodes)
165 : implicit none
166 : integer :: number_nodes
167 : ! local variables
168 : integer :: status, number_of_nodes_id
169 :
170 : ! Get the id of 'num_nodes', if such dimension is not there panic and quit :P
171 0 : status = nf90_inq_dimid(p_ncid, "num_nodes", number_of_nodes_id)
172 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
173 :
174 : ! How many values for 'num_nodes' are there?
175 0 : status = nf90_inquire_dimension(p_ncid, number_of_nodes_id, len = number_nodes)
176 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
177 :
178 0 : end function get_number_of_nodes
179 :
180 :
181 : !======================================================================
182 : ! function get_number_of_element_blocks()
183 : !======================================================================
184 0 : function get_number_of_element_blocks() result(number_element_blocks)
185 :
186 : integer :: number_element_blocks
187 : ! local variables
188 : integer :: status, number_of_element_blocks_id
189 :
190 : ! Get the id of 'num_el_blk', if such dimension is not there panic and quit :P
191 0 : status = nf90_inq_dimid(p_ncid, "num_el_blk", number_of_element_blocks_id)
192 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
193 :
194 : ! How many values for 'num_el_blk' are there?
195 0 : status = nf90_inquire_dimension(p_ncid, number_of_element_blocks_id, len = number_element_blocks)
196 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
197 :
198 0 : if (number_element_blocks /= 1) then
199 0 : if (number_element_blocks /= 6 ) then
200 0 : call endrun('Reading cube-sphere from input file is not supported')
201 : else
202 0 : call endrun('Number of elements blocks not exactly 1 (sphere) or 6 (cube)')
203 : endif
204 : endif
205 :
206 0 : end function get_number_of_element_blocks
207 :
208 : !======================================================================
209 : ! function get_number_of_elements_per_face()
210 : !======================================================================
211 0 : function get_number_of_elements_per_face() result(number_elements_per_face)
212 :
213 : integer :: number_elements_per_face
214 :
215 : integer :: face_num ! For each of the face, we get the information
216 : character(len=MXSTLN) :: element_type ! Each face is composed of elements of certain type
217 : integer :: number_elements_in_face ! How many elements in this face
218 : integer :: num_nodes_per_elem ! How many nodes in each element
219 : integer :: number_of_attributes ! How many attributes in the face
220 :
221 : integer :: status, dimension_id
222 :
223 0 : if (p_number_blocks == 0) then
224 0 : call endrun('get_number_of_elements_per_face called before MeshOpen')
225 0 : else if (p_number_blocks == 1) then ! we are in the presence of a sphere
226 : ! First we get sure the number of nodes per element is four
227 0 : status = nf90_inq_dimid(p_ncid, "num_nod_per_el1", dimension_id)
228 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
229 0 : status = nf90_inquire_dimension(p_ncid, dimension_id, len = num_nodes_per_elem)
230 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
231 0 : if (num_nodes_per_elem /= 4) call endrun('Number of nodes per element is not four')
232 : ! now we check how many elements there are in the face
233 0 : status = nf90_inq_dimid(p_ncid, "num_el_in_blk1", dimension_id)
234 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
235 0 : status = nf90_inquire_dimension(p_ncid, dimension_id, len = number_elements_in_face)
236 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
237 0 : number_elements_per_face = number_elements_in_face
238 0 : else if (p_number_blocks == 6) then ! we are in the presence of a cube-sphere
239 0 : call endrun('Reading a mesh for a cube-sphere is not supported')
240 : else
241 0 : call endrun('Number of elements blocks not exactly 1 (sphere) or 6 (cube)')
242 : end if
243 :
244 0 : end function get_number_of_elements_per_face
245 :
246 : !======================================================================
247 : ! subroutine get_face_connectivity
248 : !======================================================================
249 0 : subroutine get_face_connectivity()
250 :
251 : integer :: var_id, status
252 :
253 0 : status = nf90_inq_varid(p_ncid, "connect1", var_id)
254 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
255 0 : status = nf90_get_var(p_ncid, var_id, p_connectivity)
256 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
257 0 : end subroutine get_face_connectivity
258 :
259 : !======================================================================
260 : ! subroutine get_node_multiplicity
261 : !======================================================================
262 0 : subroutine get_node_multiplicity(node_multiplicity)
263 : use dimensions_mod, only : max_elements_attached_to_node
264 :
265 : integer, intent(out) :: node_multiplicity(:)
266 : integer :: node_num(4)
267 :
268 : integer :: k, number_nodes
269 :
270 0 : node_multiplicity(:) = 0
271 0 : number_nodes = SIZE(node_multiplicity)
272 : ! check this external buffer was allocated correctly
273 0 : if (number_nodes /= p_number_nodes) call endrun('Number of nodes does not matches size of node multiplicity array')
274 : ! for each node, we have for four other nodes
275 :
276 0 : if (minval(p_connectivity) < 1 .or. number_nodes < maxval(p_connectivity)) then
277 0 : call endrun('get_node_multiplicity: Node number less than 1 or greater than max.')
278 : end if
279 :
280 0 : do k=1,p_number_elements_per_face
281 0 : node_num = p_connectivity(:,k)
282 0 : node_multiplicity(node_num) = node_multiplicity(node_num) + 1
283 : enddo
284 :
285 0 : if (minval(node_multiplicity) < 3 .or. max_elements_attached_to_node < maxval(node_multiplicity)) then
286 0 : print *, 'minval(node_multiplicity)', minval(node_multiplicity)
287 0 : print *, 'maxval(node_multiplicity)', maxval(node_multiplicity),&
288 0 : ' and max_elements_attached_to_node ',max_elements_attached_to_node
289 0 : call endrun('get_node_multiplicity: Number of elements attached to node less than 3 or greater than maximum.')
290 : endif
291 :
292 0 : end subroutine get_node_multiplicity
293 :
294 : !======================================================================
295 : ! subroutine get_node_coordinates ()
296 : !======================================================================
297 0 : subroutine get_node_coordinates ()
298 :
299 : integer :: var_id, status
300 :
301 0 : status = nf90_inq_varid(p_ncid, "coord", var_id)
302 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
303 0 : status = nf90_get_var(p_ncid, var_id, p_node_coordinates)
304 0 : if(status /= nf90_NoErr) call handle_error(status, __FILE__, __LINE__)
305 0 : end subroutine get_node_coordinates
306 :
307 : ! ================================================================================
308 : !
309 : ! -----------------Internal private routines that do not use netCDF IO -----------
310 : !
311 : ! ================================================================================
312 :
313 : !======================================================================
314 : ! subroutine get_2D_sub_coordinate_indexes
315 : !======================================================================
316 0 : subroutine get_2D_sub_coordinate_indexes(x, y, sgnx, sgny, face_no)
317 : implicit none
318 : integer, intent(in) :: face_no
319 : integer, intent(out) :: x,y
320 : integer, intent(out) :: sgnx, sgny
321 0 : if (face_no == 1 .or. face_no == 3) then
322 0 : x = 2
323 0 : y = 3
324 0 : else if (face_no == 2 .or. face_no == 4) then
325 0 : x = 1
326 0 : y = 3
327 : else
328 0 : x = 2
329 0 : y = 1
330 : endif
331 0 : if (face_no == 1 .or. face_no == 4 .or. face_no == 5) then
332 0 : sgnx = 1
333 0 : sgny = 1
334 0 : else if (face_no == 2 .or. face_no == 3) then
335 0 : sgnx = -1
336 0 : sgny = 1
337 : else
338 0 : sgnx = 1
339 0 : sgny = -1
340 : endif
341 0 : end subroutine get_2D_sub_coordinate_indexes
342 :
343 :
344 :
345 : !======================================================================
346 : ! subroutine mesh_connectivity(connect)
347 : !
348 : ! puts the transpose of p_connectivity into connect
349 : !======================================================================
350 :
351 0 : subroutine mesh_connectivity (connect)
352 :
353 : integer, intent(out) :: connect(p_number_elements,4)
354 :
355 : integer :: k, j
356 :
357 0 : if (0 == p_number_blocks) call endrun('mesh_connectivity called before MeshOpen')
358 0 : j=0
359 0 : do k=1, p_number_elements_per_face
360 0 : j=j+1
361 0 : connect(j,:) = p_connectivity(:,k)
362 : enddo
363 :
364 0 : if (j /= p_number_elements) call endrun('mesh_connectivity: Number of elements in side sets not equal to total elements')
365 :
366 0 : if (minval(connect) < 1 .or. maxval(connect) > p_number_nodes) then
367 0 : call endrun('mesh_connectivity: Node number out of bounds')
368 : end if
369 :
370 0 : end subroutine mesh_connectivity
371 : !======================================================================
372 : ! subroutine create_index_table()
373 : !
374 : ! this is needed to detremine side and corner neighbors
375 : !======================================================================
376 :
377 0 : subroutine create_index_table(index_table, element_nodes)
378 : use dimensions_mod, only : max_elements_attached_to_node
379 :
380 : integer, allocatable, intent(inout) :: index_table(:,:)
381 : integer , intent(in) :: element_nodes(p_number_elements, 4)
382 : integer :: cnt, cnt_index, node
383 : integer :: k, ll
384 :
385 : !Create an index table so that we can find neighbors on O(n)
386 : ! so for each node, we want to know which elements it is part of
387 0 : allocate(index_table(p_number_nodes, max_elements_attached_to_node + 1))
388 :
389 : !the last column in the index table is a count of the number of elements
390 0 : index_table = 0
391 :
392 0 : cnt_index = max_elements_attached_to_node + 1
393 :
394 0 : do k=1,p_number_elements
395 0 : do ll=1,4
396 0 : node = element_nodes(k, ll) !the node
397 0 : cnt = index_table(node, cnt_index) !how many elements for that node already in table
398 0 : cnt = cnt + 1 !increment since we are adding an element
399 0 : if (cnt > max_elements_attached_to_node) then
400 0 : call endrun('Found a node in too many elements.')
401 : endif
402 0 : index_table(node, cnt_index) = cnt
403 0 : index_table(node, cnt) = k !put the element in the indextable
404 : enddo
405 : enddo
406 :
407 0 : end subroutine create_index_table
408 :
409 : !======================================================================
410 : ! subroutine find_side_neighbors()
411 : !
412 : ! find the element neighbors to the n,s,e,w and put them in GridVertex_t
413 : ! (only 1 neighbor to the n,s,e,w)
414 : !======================================================================
415 0 : subroutine find_side_neighbors (GridVertex, normal_to_homme_ordering, element_nodes, edge_wgt, index_table)
416 : use coordinate_systems_mod, only : cartesian3D_t
417 : use gridgraph_mod, only : GridVertex_t
418 : use dimensions_mod, only : max_elements_attached_to_node
419 :
420 : integer , intent(in) :: normal_to_homme_ordering(8)
421 : integer , intent(in) :: element_nodes(p_number_elements, 4)
422 : integer , intent(in) :: edge_wgt
423 : integer , intent(in) :: index_table(:,:)
424 : type (GridVertex_t) , intent(inout) :: GridVertex(:)
425 :
426 : integer :: i_node(2), my_node(2)
427 : integer :: neighbor, direction
428 : integer :: j,k,ll,i, m
429 : integer :: i_elem, jump, end_i
430 : integer :: loc, cnt_index, a_count(2)
431 : logical :: found
432 0 : if (0 == p_number_blocks) call endrun('find_side_neighbors called before MeshOpen')
433 :
434 :
435 : !the last column in the index table is a count of the number of elements
436 0 : cnt_index = max_elements_attached_to_node + 1
437 :
438 : !use index table to find neighbors
439 0 : do k=1,p_number_elements ! for each element k
440 : !set the side weights
441 0 : GridVertex(k)%nbrs_wgt(1:4) = edge_wgt
442 0 : do ll=1,4 !loop through the four sides
443 :
444 0 : jump = normal_to_homme_ordering(ll)
445 0 : loc = GridVertex(k)%nbrs_ptr(jump)
446 :
447 0 : if (GridVertex(k)%nbrs(loc) == 0) then !if side is not set yet, then
448 : !look for side element
449 0 : found = .false.
450 0 : neighbor = 0
451 :
452 0 : my_node(1) = element_nodes(k, ll)
453 0 : a_count(1) = index_table(my_node(1), cnt_index)
454 0 : my_node(2) = element_nodes(k, mod(ll,4)+1)
455 0 : a_count(2) = index_table(my_node(2), cnt_index)
456 :
457 : !loop through the elements that are in the index table for each node
458 : !and find the element number and direction of the side neighbor
459 0 : do m = 1,2
460 0 : if (found) exit
461 0 : end_i = a_count(m)
462 0 : do i = 1, end_i
463 0 : if (found) exit
464 0 : i_elem = index_table(my_node(m),i)
465 0 : if (i_elem /= k) then !k is the element we are setting sides for
466 0 : do j=1,4 !loop through each of i_elem's four sides
467 0 : i_node(1) = element_nodes(i_elem, j)
468 0 : i_node(2) = element_nodes(i_elem, mod(j,4)+1)
469 0 : if ( (i_node(1) == my_node(2) .and. i_node(2) == my_node(1)) .or. &
470 0 : (i_node(1) == my_node(1) .and. i_node(2) == my_node(2)) ) then
471 : neighbor = i_elem
472 : direction = j
473 : found = .true.
474 : !found a match
475 : exit
476 : end if
477 : end do ! j loop
478 : end if
479 : enddo ! i loop
480 : enddo !m loop
481 :
482 0 : if (neighbor == 0) call endrun('find_side_neighbor: Neighbor not found! Every side should have a neighbor.')
483 :
484 0 : GridVertex(k)%nbrs(loc) = neighbor
485 0 : jump = normal_to_homme_ordering(direction)
486 0 : loc = GridVertex(neighbor)%nbrs_ptr(jump)
487 0 : GridVertex(neighbor)%nbrs(loc)= k
488 : endif
489 : enddo ! ll loop => 4 sides
490 : enddo ! k loop: each element
491 :
492 0 : do k=1,p_number_elements
493 0 : do ll=1,4
494 0 : if ( 0 == GridVertex(k)%nbrs(ll)) then
495 0 : call endrun('Found one side of one element witout a neighbor. Bummer!')
496 : end if
497 : end do
498 : end do
499 :
500 0 : end subroutine find_side_neighbors
501 :
502 : !======================================================================
503 : ! function smallest_diameter_element
504 : !======================================================================
505 :
506 0 : function smallest_diameter_element(element_nodes) result(min_diameter)
507 :
508 : integer ,intent(in) :: element_nodes(:,:)
509 :
510 : integer :: i, j
511 : integer :: node_numbers(4)
512 : real(kind=r8) :: coordinates (4,3)
513 : real(kind=r8) :: x(3), y(3), r(3), d, min_diameter
514 :
515 0 : if (SIZE(element_nodes,dim=1) /= p_number_elements) then
516 : call endrun('smallest_diameter_element:Element count check failed in &
517 0 : &exodus_mesh. Connectivity array length not equal to number of elements.')
518 : end if
519 0 : if ( p_number_elements_per_face /= p_number_elements) then
520 : call endrun('smallest_diameter_element: Element count check failed in &
521 0 : &exodus_mesh. Element array length not equal to sum of face.')
522 : end if
523 :
524 0 : min_diameter = 9999999.0_r8
525 0 : do i=1, p_number_elements
526 0 : node_numbers = element_nodes(i,:)
527 0 : coordinates = p_node_coordinates(node_numbers,:)
528 : ! smallest side length
529 0 : do j=1,4
530 0 : x = coordinates(j ,:)
531 0 : y = coordinates(1+MOD(j,4),:)
532 0 : r = x-y
533 0 : d = dot_product(r,r)
534 0 : if (d < min_diameter ) then
535 0 : min_diameter = d
536 : end if
537 : end do
538 : ! smallest diameter length
539 0 : do j=1,2
540 0 : x = coordinates(j ,:)
541 0 : y = coordinates(2+MOD(j,4),:)
542 0 : r = x-y
543 0 : d = dot_product(r,r)
544 0 : if (d < min_diameter ) then
545 0 : min_diameter = d
546 : end if
547 : end do
548 : enddo
549 0 : min_diameter = SQRT(min_diameter)
550 0 : end function smallest_diameter_element
551 :
552 : !======================================================================
553 : ! subroutine cube_to_cube_coordinates
554 : !======================================================================
555 :
556 0 : subroutine cube_to_cube_coordinates (cube_coor, node_coor, face_number)
557 :
558 : real(kind=r8), intent(in) :: node_coor(4,3)
559 : integer, intent(in) :: face_number
560 : real(kind=r8), intent(out) :: cube_coor(4,2)
561 :
562 : integer :: x_index, y_index, sgnx, sgny
563 0 : call get_2D_sub_coordinate_indexes(x_index, y_index, sgnx, sgny, face_number)
564 0 : cube_coor(:,1) = sgnx*node_coor(:,x_index)
565 0 : cube_coor(:,2) = sgny*node_coor(:,y_index)
566 0 : end subroutine cube_to_cube_coordinates
567 :
568 :
569 : !======================================================================
570 : ! subroutine sphere_to_cube_coordinates
571 : !======================================================================
572 :
573 0 : subroutine sphere_to_cube_coordinates (cube_coor, node_coor, face_number)
574 : use coordinate_systems_mod, only : cartesian2d_t, change_coordinates, sphere2cubedsphere
575 : implicit none
576 : real(kind=r8), intent(in) :: node_coor(4,3)
577 : integer, intent(in) :: face_number
578 : real(kind=r8), intent(out) :: cube_coor(4,2)
579 : integer :: i
580 : type(cartesian2d_t) :: cart(4)
581 :
582 0 : do i=1,4
583 0 : cart(i) = sphere2cubedsphere(change_coordinates(node_coor(i,:)), face_number)
584 : end do
585 0 : cube_coor(:,1) = cart(:)%x
586 0 : cube_coor(:,2) = cart(:)%y
587 0 : end subroutine sphere_to_cube_coordinates
588 :
589 :
590 : !======================================================================
591 : ! subroutine cube_face_element_centroids
592 : !======================================================================
593 :
594 0 : subroutine cube_face_element_centroids(centroids, face_numbers, element_nodes)
595 :
596 : integer , intent(in) :: element_nodes(:,:)
597 : integer, intent(in) :: face_numbers (p_number_elements)
598 : real(kind=r8),intent(out) :: centroids (p_number_elements,2)
599 : real(kind=r8) :: coordinates(4,3)
600 : real(kind=r8) :: cube_coor (4,2)
601 : integer :: i, node_numbers(4)
602 :
603 0 : if (0 == p_number_blocks) call endrun('cube_face_element_centroids called before MeshOpen')
604 0 : if (SIZE(element_nodes,dim=1) /= p_number_elements) then
605 : call endrun('cube_face_element_centroids:Element count check failed in &
606 0 : &exodus_mesh. Connectivity array length not equal to number of elements.')
607 : end if
608 0 : if ( p_number_elements_per_face /= p_number_elements ) then
609 : call endrun('cube_face_element_centroids: Element count check failed in &
610 0 : &exodus_mesh. Element array length not equal to sum of face.')
611 : end if
612 :
613 0 : do i=1, p_number_elements
614 0 : node_numbers = element_nodes(i,:)
615 0 : coordinates = p_node_coordinates(node_numbers,:)
616 0 : if (6 == p_number_blocks) then
617 0 : call cube_to_cube_coordinates (cube_coor, coordinates, face_numbers(i))
618 : else
619 0 : call sphere_to_cube_coordinates (cube_coor, coordinates, face_numbers(i))
620 : end if
621 0 : centroids(i,:) = SUM(cube_coor,dim=1)/4.0_r8
622 : enddo
623 0 : end subroutine cube_face_element_centroids
624 :
625 : !======================================================================
626 : ! subroutine initialize_space_filling_curve
627 : !======================================================================
628 0 : subroutine initialize_space_filling_curve(GridVertex, element_nodes)
629 : use gridgraph_mod, only : GridVertex_t
630 : use spacecurve_mod, only : GenspaceCurve
631 :
632 : type (GridVertex_t), intent(inout) :: GridVertex(:)
633 : integer , intent(in) :: element_nodes(:,:)
634 :
635 0 : integer,allocatable :: Mesh2(:,:),Mesh2_map(:,:),sfcij(:,:)
636 :
637 0 : real(kind=r8) :: centroids(p_number_elements,2)
638 0 : integer :: face_numbers(p_number_elements)
639 : real(kind=r8) :: x, y, h
640 : integer :: i, j, i2, j2, ne, ne2
641 : integer :: sfc_index, face
642 :
643 0 : if (SIZE(GridVertex) /= p_number_elements) then
644 : call endrun('initialize_space_filling_curve:Element count check failed &
645 0 : &in exodus_mesh. Vertex array length not equal to number of elements.')
646 : end if
647 0 : if (SIZE(element_nodes,dim=1) /= p_number_elements) then
648 : call endrun('initialize_space_filling_curve:Element count check failed &
649 0 : &in exodus_mesh. Connectivity array length not equal to number of elements.')
650 : end if
651 :
652 0 : face_numbers(:) = GridVertex(:)%face_number
653 0 : h = smallest_diameter_element ( element_nodes)
654 :
655 0 : call cube_face_element_centroids (centroids, face_numbers, element_nodes)
656 :
657 0 : if (h<.00001_r8) then
658 0 : call endrun('initialize_space_filling_curve: Unreasonably small element found. less than .00001')
659 : end if
660 :
661 0 : ne = CEILING(0.5_r8*PI/(h/2));
662 :
663 : ! find the smallest ne2 which is a power of 2 and ne2>ne
664 0 : ne2=2**ceiling( log(real(ne))/log(2._r8) )
665 0 : if (ne2<ne) call endrun('initialize_space_filling_curve: Fatel SFC error')
666 :
667 0 : allocate(Mesh2(ne2,ne2))
668 0 : allocate(Mesh2_map(ne2,ne2))
669 0 : allocate(sfcij(0:ne2*ne2,2))
670 :
671 : ! create a reverse index array for Mesh2
672 : ! j = Mesh2(i,j)
673 : ! (i,j) = (sfcij(j,1),sfci(j,2))
674 0 : call GenspaceCurve(Mesh2) ! SFC partition for ne2
675 0 : do j2=1,ne2
676 0 : do i2=1,ne2
677 0 : j=Mesh2(i2,j2)
678 0 : sfcij(j,1)=i2
679 0 : sfcij(j,2)=j2
680 : enddo
681 : enddo
682 :
683 :
684 0 : GridVertex(:)%SpaceCurve=-1
685 : sfc_index = 0
686 0 : do face = 1,nfaces
687 : ! associate every element on the ne x ne mesh (Mesh)
688 : ! with its closest element on the ne2 x ne2 mesh (Mesh2)
689 : ! Store this as a map from Mesh2 -> Mesh in Mesh2_map.
690 : ! elements in Mesh2 which are not mapped get assigned a value of 0
691 0 : Mesh2_map=0
692 0 : do i=1,p_number_elements
693 0 : if (face_numbers(i) == face ) then
694 0 : x = centroids(i,1)
695 0 : y = centroids(i,2)
696 : ! map this element to an (i2,j2) element
697 : ! [ -PI/4, PI/4 ] -> [ 0, ne2 ]
698 0 : i2=nint( (0.5_r8 + 2.0_r8*x/PI)*ne2 + 0.5_r8 )
699 0 : j2=nint( (0.5_r8 + 2.0_r8*y/PI)*ne2 + 0.5_r8 )
700 0 : if (face == 4 .or. face == 6 ) i2 = ne2-i2+1
701 0 : if (face == 1 .or. face == 2 .or. face == 6) j2 = ne2-j2+1
702 0 : if (i2<1 ) i2=1
703 0 : if (i2>ne2) i2=ne2
704 0 : if (j2<1 ) j2=1
705 0 : if (j2>ne2) j2=ne2
706 0 : Mesh2_map(i2,j2)=i
707 : end if
708 : end do
709 :
710 : ! generate a SFC for Mesh with the same ordering as the
711 : ! elements in Mesh2 which map to Mesh.
712 0 : do j=0,ne2*ne2-1
713 0 : i2=sfcij(j,1)
714 0 : j2=sfcij(j,2)
715 0 : i=Mesh2_map(i2,j2)
716 0 : if (i/=0) then
717 : ! (i2,j2) element maps to element
718 0 : GridVertex(i)%SpaceCurve=sfc_index
719 0 : sfc_index=sfc_index+1
720 : endif
721 : enddo
722 : enddo
723 0 : deallocate(Mesh2)
724 0 : deallocate(Mesh2_map)
725 0 : deallocate(sfcij)
726 :
727 0 : if (minval(GridVertex(:)%SpaceCurve) == -1) then
728 0 : do i=1,p_number_elements
729 0 : if (-1==GridVertex(i)%SpaceCurve) then
730 0 : write (*,*) " Error in projecting element ",i," to space filling curve."
731 0 : write (*,*) " Face:",face_numbers(i)
732 0 : write (*,*) " Centroid:",centroids(i,:)
733 : end if
734 : end do
735 0 : call endrun('initialize_space_filling_curve: Vertex not on SpaceCurve')
736 : end if
737 :
738 0 : end subroutine initialize_space_filling_curve
739 :
740 : !======================================================================
741 : ! subroutine find_corner_neighbors
742 : !======================================================================
743 :
744 0 : subroutine find_corner_neighbors (GridVertex, normal_to_homme_ordering, element_nodes, corner_wgt, index_table)
745 : use gridgraph_mod, only : GridVertex_t
746 : use dimensions_mod, only : max_elements_attached_to_node, max_corner_elem
747 : use control_mod, only: north, south, east, west, neast,seast, nwest,swest
748 :
749 : type (GridVertex_t), intent(inout) :: GridVertex(:)
750 : integer , intent(in) :: normal_to_homme_ordering(8)
751 : integer , intent(in) :: element_nodes(p_number_elements, 4)
752 : integer , intent(in) :: corner_wgt
753 : integer , intent(in) :: index_table(:,:)
754 :
755 0 : integer :: node_elements (2*max_elements_attached_to_node)
756 0 : integer :: elem_neighbor (4*max_elements_attached_to_node)
757 : integer :: nbr_cnt(4)
758 : integer :: elem_nbr_start, start
759 : integer :: i, j, k, ll, jj, kk
760 : integer :: node, loc, cnt, cnt_index
761 0 : integer :: corner_array(max_corner_elem), orig_pos(max_corner_elem)
762 0 : integer :: face_array(max_corner_elem), a_corner_elems(max_corner_elem)
763 : integer :: corner_sides(2)
764 : integer :: side_elem, corner_elem, tmp_s
765 :
766 : !the last column in the index table is a count of the number of elements
767 0 : cnt_index = max_elements_attached_to_node + 1
768 :
769 0 : do i=1, p_number_elements !loop through all elements
770 0 : node_elements(:) = 0
771 0 : elem_neighbor(:) = 0
772 0 : elem_nbr_start = 0
773 0 : nbr_cnt(:) = 0
774 :
775 0 : do j=1,4 !check each of the 4 nodes at the element corners
776 0 : node = element_nodes(i,j)
777 0 : cnt = index_table(node, cnt_index)
778 0 : if (cnt < 3 .or. max_elements_attached_to_node < cnt) then
779 0 : call endrun('find_corner_neighbors: Number of elements attached to node less than 3 or greater than maximum.')
780 : endif
781 :
782 0 : node_elements(1:cnt) = index_table(node, 1:cnt)
783 :
784 : !now node_elements contains the element neighbors to that node - so grab the
785 : ! corner neighbors - these are the ones that are not already side neighbors (or myself)
786 0 : k = 0
787 0 : do ll=1,cnt
788 0 : if ( i /= node_elements(ll) .and. & !not me
789 0 : GridVertex(i)%nbrs(1) /= node_elements(ll) .and. & !not side 1
790 0 : GridVertex(i)%nbrs(2) /= node_elements(ll) .and. & ! etc ...
791 0 : GridVertex(i)%nbrs(3) /= node_elements(ll) .and. &
792 0 : GridVertex(i)%nbrs(4) /= node_elements(ll)) then
793 0 : k = k + 1
794 0 : elem_neighbor(elem_nbr_start + k) = node_elements(ll)
795 : end if
796 : end do ! end of ll loop for multiplicity
797 :
798 : !keep track of where we are starting in elem_neighbor for each corner j
799 0 : elem_nbr_start = elem_nbr_start + k
800 0 : nbr_cnt(j) = k !how many neighbors in this corner
801 : end do ! end of j loop through 4 nodes
802 :
803 :
804 : ! now that we have done the 4 corners we can populate nbrs and nbrs_ptr
805 : ! with the corners in the proper order (clockwise) in neighbors
806 : ! also we can add the corner weight
807 :
808 0 : do j=5,8 !loop through 4 corners
809 : elem_nbr_start = 1
810 : !easiest to do the corner in ascending order - find loc
811 0 : do jj = 5,8
812 0 : ll = normal_to_homme_ordering(jj)
813 0 : if (j == ll) then
814 : loc = jj
815 : exit
816 : end if
817 0 : elem_nbr_start = elem_nbr_start + nbr_cnt(jj-4)
818 : end do
819 :
820 0 : start = GridVertex(i)%nbrs_ptr(j)
821 0 : cnt = nbr_cnt(loc - 4)
822 0 : GridVertex(i)%nbrs_ptr(j+1) = start + cnt
823 :
824 0 : if (cnt > 0) then
825 0 : GridVertex(i)%nbrs(start : start + cnt-1) = &
826 0 : elem_neighbor(elem_nbr_start : elem_nbr_start + cnt -1)
827 0 : GridVertex(i)%nbrs_face(start : start + cnt - 1) = &
828 0 : GridVertex(elem_neighbor(elem_nbr_start : elem_nbr_start + cnt -1))%face_number
829 0 : GridVertex(i)%nbrs_wgt(start : start + cnt-1) = corner_wgt
830 :
831 : end if
832 :
833 : ! within each corner neighbor, lets list the corners in clockwise order
834 0 : if (cnt > 1) then !cnt is the number of neighbors in this corner j
835 : !there can be at most max_corner element of these
836 :
837 0 : a_corner_elems = 0
838 0 : a_corner_elems(1:cnt) = elem_neighbor(elem_nbr_start : elem_nbr_start + cnt -1)
839 : !corner-sides(2) is clockwise of corner_side(1)
840 0 : corner_array= 0
841 0 : orig_pos = 0
842 0 : select case (j)
843 : case(neast)
844 0 : corner_sides(1) = north
845 0 : corner_sides(2) = east
846 : case(seast)
847 0 : corner_sides(1) = east
848 0 : corner_sides(2) = south
849 : case(swest)
850 0 : corner_sides(1) = south
851 0 : corner_sides(2) = west
852 : case(nwest)
853 0 : corner_sides(1) = west
854 0 : corner_sides(2) = north
855 : end select
856 :
857 : !so the first element to list touches corner_sides(1) element
858 0 : side_elem = GridVertex(i)%nbrs(corner_sides(1))
859 :
860 : !loop though the corner elements and see if any have a side neighbor
861 : !that = side_elem
862 0 : do k = 1,cnt !number of corner elements
863 0 : corner_elem = a_corner_elems(k)
864 0 : do kk = 1,4 !number of sides to check
865 0 : loc = GridVertex(corner_elem)%nbrs_ptr(kk)
866 0 : tmp_s = GridVertex(corner_elem)%nbrs(loc)
867 0 : if (tmp_s == side_elem) then
868 0 : corner_array(1) = corner_elem
869 0 : orig_pos(1) = k
870 0 : exit
871 : endif
872 : enddo
873 0 : if (corner_array(1)> 0) exit
874 : enddo
875 0 : if (corner_array(1)==0) then
876 0 : print *, i, cnt
877 0 : call endrun('find_corner_neighbors (1) : mistake finding corner neighbor order')
878 : endif
879 :
880 : !if cnt == 2, we are done (we know the order of neighbors)
881 0 : if (cnt ==2) then
882 0 : if (corner_array(1) == a_corner_elems(1)) then
883 0 : corner_array(2) = a_corner_elems(2)
884 0 : orig_pos(2) = 2
885 : else
886 0 : corner_array(2) = a_corner_elems(1)
887 0 : orig_pos(2) = 1
888 : end if
889 : else !cnt = 3 or 4
890 : !find which corner element borders corner_sides(2)
891 0 : side_elem = GridVertex(i)%nbrs(corner_sides(2))
892 0 : do k = 1,cnt
893 0 : corner_elem = a_corner_elems(k)
894 0 : do kk = 1,4
895 0 : loc = GridVertex(corner_elem)%nbrs_ptr(kk)
896 0 : tmp_s = GridVertex(corner_elem)%nbrs(loc)
897 0 : if (tmp_s == side_elem) then
898 0 : corner_array(4) = corner_elem
899 0 : orig_pos(4) = k
900 0 : exit
901 : endif
902 : enddo
903 0 : if (corner_array(4)> 0) exit
904 : enddo
905 0 : if (corner_array(4)==0 .or. corner_array(4) == corner_array(1)) then
906 0 : print *, i, cnt
907 0 : call endrun('find_corner_neighbors (2) : mistake finding corner neighbor order')
908 : endif
909 :
910 : !now if cnt = 3 then we are done
911 0 : if (cnt ==3) then
912 0 : corner_array(3) = corner_array(4)
913 0 : orig_pos(3) = orig_pos(4)
914 :
915 0 : do k = 1,cnt !find the "middle" element
916 0 : if (k /= orig_pos(1) .and. k /= orig_pos(3)) then
917 0 : orig_pos(2) = k
918 0 : corner_array(2) = a_corner_elems(k)
919 0 : exit
920 : endif
921 : enddo
922 : else !cnt = 4
923 : !which of the two unassigned elements borders the element in
924 : !corner_array(1) => put in corner_array(2)
925 0 : side_elem = corner_array(1)
926 :
927 0 : do k = 1,cnt
928 0 : corner_elem = a_corner_elems(k)
929 0 : if (corner_elem == corner_array(4) .or. corner_elem == corner_array(1)) then
930 : cycle
931 : else
932 0 : do kk = 1,4 !check each side
933 0 : loc = GridVertex(corner_elem)%nbrs_ptr(kk)
934 0 : tmp_s = GridVertex(corner_elem)%nbrs(loc)
935 0 : if (tmp_s == side_elem) then
936 0 : corner_array(2) = corner_elem
937 0 : orig_pos(2) = k
938 0 : exit
939 : endif
940 : enddo
941 : endif
942 0 : if (corner_array(2)> 0) exit
943 : enddo
944 : !now put the remaining one in pos 3
945 0 : do k = 1,cnt
946 0 : corner_elem = a_corner_elems(k)
947 : if (corner_elem /= corner_array(4) .and. corner_elem /= &
948 0 : corner_array(2) .and. corner_elem /= corner_array(1)) then
949 0 : corner_array(3) = corner_elem
950 0 : orig_pos(3) = k
951 0 : exit
952 : endif
953 : enddo
954 : endif ! end of cnt=4
955 : endif! end of not cnt=2
956 :
957 : !now re-set the elements in this corner
958 0 : GridVertex(i)%nbrs(start : start + cnt-1) = corner_array(1:cnt)
959 : !nbrs_wgt are the same - nothing to do
960 : !fix neighbors face
961 0 : do k = 1,cnt
962 0 : face_array(k) = GridVertex(i)%nbrs_face(start + orig_pos(k) - 1)
963 : end do
964 0 : GridVertex(i)%nbrs_face(start : start + cnt - 1) = face_array(1:cnt)
965 : endif !end of cnt > 1 loop for corners
966 :
967 : end do !j loop through each corner
968 :
969 : end do ! end of i loop through elements
970 0 : end subroutine find_corner_neighbors
971 :
972 : ! ================================================================================
973 : !
974 : ! -------------------------------Public Methods-----------------------------------
975 : !
976 : ! ================================================================================
977 :
978 : !======================================================================
979 : ! subroutine MeshOpen
980 : !======================================================================
981 :
982 0 : subroutine MeshOpen(mesh_file_name, par)
983 : use parallel_mod, only: parallel_t
984 : use cam_logfile, only: iulog
985 :
986 : character (len=*), intent(in) :: mesh_file_name
987 : type (parallel_t), intent(in) :: par
988 :
989 0 : integer, allocatable :: node_multiplicity(:)
990 : integer :: k
991 :
992 0 : p_mesh_file_name = mesh_file_name
993 0 : call open_mesh_file ()
994 :
995 0 : p_number_elements = get_number_of_elements ()
996 0 : p_number_nodes = get_number_of_nodes ()
997 0 : p_number_blocks = get_number_of_element_blocks ()
998 0 : p_number_dimensions = get_number_of_dimensions ()
999 :
1000 0 : if (p_number_dimensions /= 3) then
1001 0 : call endrun('The number of dimensions must be 3, otherwise the mesh algorithms will not work')
1002 : endif
1003 :
1004 : ! Only spheres are allowed in input files.
1005 0 : if (par%masterproc) then
1006 0 : if (p_number_blocks == 1) then
1007 0 : write(iulog,*) "Since the mesh file has only one block, it is assumed to be a sphere."
1008 : endif
1009 : end if
1010 :
1011 0 : if (p_number_blocks /= 1) then
1012 0 : call endrun('Number of elements blocks not exactly 1 (sphere)')
1013 : end if
1014 :
1015 0 : p_number_elements_per_face = get_number_of_elements_per_face()
1016 : ! Because all elements are in one face, this value must match p_number_elements
1017 0 : if ( p_number_elements /= p_number_elements_per_face) then
1018 0 : call endrun('The value of the total number of elements does not match all the elements found in face 1')
1019 : end if
1020 :
1021 0 : allocate( p_connectivity(4,p_number_elements_per_face) )
1022 0 : p_connectivity(:,:)=0
1023 : ! extract the connectivity from the netcdf file
1024 0 : call get_face_connectivity()
1025 :
1026 0 : allocate(node_multiplicity(p_number_nodes))
1027 0 : call get_node_multiplicity(node_multiplicity)
1028 :
1029 : ! tricky: For each node with multiplicity n, there are n(n-1) neighbor links
1030 : ! created. But this counts each edge twice, so: n(n-1) -n
1031 : ! Should be the same as SUM(SIZE(GridVertex(i)%nbrs(j)%n),i=1:p_number_elements,j=1:8)
1032 : ! p_number_neighbor_edges = dot_product(mult,mult) - 2*sum(mult)
1033 0 : p_number_neighbor_edges = 0
1034 0 : do k=1,p_number_nodes
1035 0 : p_number_neighbor_edges = p_number_neighbor_edges + node_multiplicity(k)*(node_multiplicity(k)-2)
1036 : end do
1037 :
1038 0 : deallocate(node_multiplicity)
1039 :
1040 : ! allocate the space for the coordinates, this is used in many functions
1041 0 : allocate(p_node_coordinates(p_number_nodes, p_number_dimensions))
1042 0 : call get_node_coordinates()
1043 :
1044 0 : if (p_number_elements_per_face /= p_number_elements) then
1045 0 : call endrun('MeshOpen: Total number of elements not equal to the number of elements on face 1!')
1046 : end if
1047 :
1048 0 : end subroutine MeshOpen
1049 :
1050 : !======================================================================
1051 : ! subroutine MeshClose
1052 : !
1053 : ! This routine acts as a destructor cleaning the memory allocated in MeshOpen
1054 : ! which acts as a constructor allocated dynamical memory for the nodes coordinates.
1055 : !======================================================================
1056 :
1057 0 : subroutine MeshClose
1058 :
1059 : ! release memory
1060 0 : deallocate(p_node_coordinates)
1061 0 : deallocate(p_connectivity)
1062 : ! let the file go
1063 0 : call close_mesh_file ()
1064 :
1065 0 : end subroutine MeshClose
1066 :
1067 :
1068 : !======================================================================
1069 : ! subroutine MeshPrint
1070 : !======================================================================
1071 :
1072 :
1073 0 : subroutine MeshPrint(par)
1074 : use parallel_mod, only: parallel_t
1075 : use cam_logfile, only: iulog
1076 :
1077 : type (parallel_t), intent(in) :: par
1078 0 : if (par%masterproc) then
1079 0 : write(iulog,*) 'This are the values for file ', trim(p_mesh_file_name)
1080 0 : write(iulog,*) 'The value for the number of dimensions (num_dim) is ', p_number_dimensions
1081 0 : write(iulog,*) 'The number of elements in the mesh file is ', p_number_elements
1082 0 : write(iulog,*) 'The number of nodes in the mesh file is ', p_number_nodes
1083 0 : write(iulog,*) 'The number of blocks in the mesh file is ', p_number_blocks
1084 0 : write(iulog,*) 'The number of elements in the face 1 (sphere) is ', p_number_elements_per_face
1085 : if ( p_number_elements == p_number_elements) then
1086 0 : write(iulog,*) 'The value of the total number of elements does match all the elements found in face 1 (the only face)'
1087 : else
1088 : write(iulog,*) 'The value of the total number of elements does not match all the elements found in face 1'
1089 : write(iulog,*) 'This message should not be appearing, there is something wrong in the code'
1090 : endif
1091 0 : write(iulog,*) 'The number of neighbor edges ', p_number_neighbor_edges
1092 : end if
1093 :
1094 0 : end subroutine MeshPrint
1095 :
1096 : !======================================================================
1097 : ! subroutine MeshCubeTopology
1098 : !======================================================================
1099 0 : subroutine MeshCubeTopology(GridEdge, GridVertex)
1100 : use dimensions_mod, only : np
1101 : use coordinate_systems_mod, only : cartesian3D_t, cube_face_number_from_cart
1102 : use gridgraph_mod, only : GridVertex_t
1103 : use gridgraph_mod, only : GridEdge_t
1104 : use cube_mod, only : CubeSetupEdgeIndex
1105 : use gridgraph_mod, only : initgridedge, num_neighbors
1106 : use control_mod, only : north, south, east, west, neast, seast, swest, nwest
1107 :
1108 : type (GridEdge_t), intent(inout), target :: GridEdge(:)
1109 : type (GridVertex_t), intent(inout), target :: GridVertex(:)
1110 :
1111 : real(kind=r8) :: coordinates(4,3)
1112 : real(kind=r8) :: centroid(3)
1113 : type (cartesian3D_t) :: face_center
1114 :
1115 : integer :: i, j, k, ll, loc
1116 0 : integer :: element_nodes(p_number_elements, 4)
1117 : integer :: EdgeWgtP,CornerWgt
1118 : integer :: normal_to_homme_ordering(8)
1119 : integer :: node_numbers(4)
1120 0 : integer, allocatable :: index_table(:,:)
1121 :
1122 0 : normal_to_homme_ordering(1) = south
1123 0 : normal_to_homme_ordering(2) = east
1124 0 : normal_to_homme_ordering(3) = north
1125 0 : normal_to_homme_ordering(4) = west
1126 0 : normal_to_homme_ordering(5) = swest
1127 0 : normal_to_homme_ordering(6) = seast
1128 0 : normal_to_homme_ordering(7) = neast
1129 0 : normal_to_homme_ordering(8) = nwest
1130 :
1131 0 : if (SIZE(GridVertex) /= p_number_elements) then
1132 : call endrun('MeshCubeTopology: Element count check failed in exodus_mesh. &
1133 0 : &Vertex array length not equal to number of elements.')
1134 : end if
1135 0 : if (p_number_elements_per_face /= p_number_elements) then
1136 : call endrun('MeshCubeTopology: Element count check failed in exodus_mesh. &
1137 0 : &Element array length not equal to sum of face.')
1138 : end if
1139 :
1140 0 : EdgeWgtP = np
1141 0 : CornerWgt = 1
1142 :
1143 :
1144 0 : call mesh_connectivity (element_nodes)
1145 :
1146 0 : do i=1, p_number_elements
1147 0 : GridVertex(i)%number = i
1148 0 : GridVertex(i)%face_number = 0
1149 0 : GridVertex(i)%processor_number = 0
1150 0 : GridVertex(i)%SpaceCurve = 0
1151 :
1152 0 : GridVertex(i)%nbrs(:) = 0
1153 0 : GridVertex(i)%nbrs_face(:) = 0
1154 0 : GridVertex(i)%nbrs_wgt(:) = 0
1155 0 : GridVertex(i)%nbrs_wgt_ghost(:) = 1
1156 :
1157 : !each elements has one side neighbor (first 4)
1158 0 : GridVertex(i)%nbrs_ptr(1) = 1
1159 0 : GridVertex(i)%nbrs_ptr(2) = 2
1160 0 : GridVertex(i)%nbrs_ptr(3) = 3
1161 0 : GridVertex(i)%nbrs_ptr(4) = 4
1162 : !don't know about corners yet
1163 0 : GridVertex(i)%nbrs_ptr(5:num_neighbors+1) = 5
1164 :
1165 : end do
1166 :
1167 : !create index table to find neighbors
1168 0 : call create_index_table(index_table, element_nodes)
1169 :
1170 : ! side neighbors
1171 0 : call find_side_neighbors(GridVertex, normal_to_homme_ordering, element_nodes, EdgeWgtP, index_table)
1172 :
1173 : ! set vertex faces
1174 0 : do i=1, p_number_elements
1175 0 : node_numbers = element_nodes(i,:)
1176 0 : coordinates = p_node_coordinates(node_numbers,:)
1177 0 : centroid = SUM(coordinates, dim=1)/4.0_r8
1178 0 : face_center%x = centroid(1)
1179 0 : face_center%y = centroid(2)
1180 0 : face_center%z = centroid(3)
1181 0 : GridVertex(i)%face_number = cube_face_number_from_cart(face_center)
1182 : end do
1183 :
1184 : ! set side neighbor faces
1185 0 : do i=1, p_number_elements
1186 0 : do j=1,4 !look at each side
1187 0 : k = normal_to_homme_ordering(j)
1188 0 : loc = GridVertex(i)%nbrs_ptr(k)
1189 0 : ll = GridVertex(i)%nbrs(loc)
1190 0 : GridVertex(i)%nbrs_face(loc) = GridVertex(ll)%face_number
1191 : end do
1192 : end do
1193 :
1194 : ! find corner neighbor and faces (weights added also)
1195 0 : call find_corner_neighbors (GridVertex, normal_to_homme_ordering, element_nodes, CornerWgt, index_table)
1196 :
1197 : !done with the index table
1198 0 : deallocate(index_table)
1199 :
1200 :
1201 0 : call initgridedge(GridEdge,GridVertex)
1202 0 : do i=1,SIZE(GridEdge)
1203 0 : call CubeSetupEdgeIndex(GridEdge(i))
1204 : enddo
1205 :
1206 0 : call initialize_space_filling_curve(GridVertex, element_nodes)
1207 0 : end subroutine MeshCubeTopology
1208 :
1209 : !======================================================================
1210 : ! subroutine MeshSetCoordinates(elem)
1211 : !======================================================================
1212 :
1213 0 : subroutine MeshSetCoordinates(elem)
1214 : use element_mod, only: element_t
1215 :
1216 : type (element_t), intent(inout) :: elem(:)
1217 :
1218 0 : integer :: connectivity(p_number_elements,4)
1219 0 : integer :: node_multiplicity(p_number_nodes)
1220 : integer :: face_no, i, k, l
1221 : integer :: number
1222 : integer :: node_num(4)
1223 : real(kind=r8) :: coordinates(4,3)
1224 : real(kind=r8) :: cube_coor (4,2)
1225 :
1226 0 : connectivity =0
1227 0 : node_multiplicity=0
1228 :
1229 0 : call mesh_connectivity (connectivity)
1230 :
1231 0 : do k=1,p_number_elements
1232 0 : node_num = connectivity(k,:)
1233 0 : node_multiplicity(node_num(:)) = node_multiplicity(node_num(:)) + 1
1234 : end do
1235 :
1236 0 : do k=1,SIZE(elem)
1237 0 : number = elem(k)%vertex%number
1238 0 : face_no = elem(k)%vertex%face_number
1239 0 : node_num = connectivity(number,:)
1240 0 : coordinates = p_node_coordinates(node_num,:)
1241 :
1242 0 : if (6 == p_number_blocks) then
1243 0 : call cube_to_cube_coordinates (cube_coor, coordinates, face_no)
1244 : else
1245 0 : call sphere_to_cube_coordinates (cube_coor, coordinates, face_no)
1246 : end if
1247 : ! elem(k)%node_numbers = node_num
1248 : ! elem(k)%node_multiplicity(:) = node_multiplicity(node_num(:))
1249 0 : elem(k)%corners(:)%x = cube_coor(:,1)
1250 0 : elem(k)%corners(:)%y = cube_coor(:,2)
1251 : end do
1252 0 : end subroutine MeshSetCoordinates
1253 :
1254 : !======================================================================
1255 : !function MeshCubeEdgeCount()
1256 : !======================================================================
1257 0 : function MeshCubeEdgeCount() result(nedge)
1258 :
1259 : integer :: nedge
1260 0 : if (0 == p_number_blocks) call endrun('MeshCubeEdgeCount called before MeshOpenMesh')
1261 0 : if (MeshUseMeshFile) then
1262 : ! should be the same as SUM(SIZE(GridVertex(i)%nbrs(j)%n),i=1:p_number_elements,j=1:nInnerElemEdge)
1263 : ! the total number of neighbors.
1264 0 : nedge = p_number_neighbor_edges
1265 : else
1266 0 : call endrun('Error in MeshCubeEdgeCount: Should not call for non-exodus mesh file.')
1267 : endif
1268 :
1269 0 : end function MeshCubeEdgeCount
1270 :
1271 0 : function MeshCubeElemCount() result(nelem)
1272 :
1273 : integer :: nelem
1274 0 : if (0 == p_number_blocks) call endrun('MeshCubeElemCount called before MeshOpenMesh')
1275 0 : if (MeshUseMeshFile) then
1276 0 : nelem = p_number_elements
1277 : else
1278 0 : call endrun('Error in MeshCubeElemCount: Should not call for non-exodus mesh file.')
1279 : end if
1280 0 : end function MeshCubeElemCount
1281 :
1282 0 : subroutine test_private_methods
1283 : implicit none
1284 0 : integer :: element_nodes(p_number_elements, 4)
1285 0 : call mesh_connectivity (element_nodes)
1286 0 : end subroutine test_private_methods
1287 :
1288 :
1289 : end module mesh_mod
|
