1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
|
#include "confdefs.h"
#include "fdebug.h"
subroutine unifiedmesh(filename1, filename1_len, &
filename2, filename2_len, &
output, output_len)
use mpi_interfaces
use fldebug
use read_triangle
use write_triangle
use fields
use linked_lists
use intersection_finder_module
use transform_elements
use elements
use supermesh_construction
use vtk_interfaces
use unify_meshes_module
implicit none
integer, intent(in) :: filename1_len, filename2_len, output_len
character(len=filename1_len), intent(in) :: filename1
character(len=filename2_len), intent(in) :: filename2
character(len=output_len), intent(in) :: output
type(vector_field) :: positionsA, positionsB
type(ilist), dimension(:), allocatable :: map_BA
real, dimension(:), allocatable :: tri_detwei
integer :: ele_A, ele_B
type(inode), pointer :: llnode
type(vector_field) :: intersection
type(element_type) :: supermesh_shape
type(quadrature_type) :: supermesh_quad
integer :: dim
type(mesh_type) :: accum_mesh
type(vector_field) :: accum_positions, accum_positions_tmp
call set_global_debug_level(0)
positionsA = read_triangle_files(trim(filename1), quad_degree=1)
positionsB = read_triangle_files(trim(filename2), quad_degree=1)
dim = positionsA%dim
allocate(map_BA(ele_count(positionsB)))
supermesh_quad = make_quadrature(vertices=dim+1, dim=dim, degree=5)
supermesh_shape = make_element_shape(vertices=dim+1, dim=dim, degree=1, quad=supermesh_quad)
allocate(tri_detwei(supermesh_shape%ngi))
call allocate(accum_mesh, 0, 0, supermesh_shape, "AccumulatedMesh")
call allocate(accum_positions, dim, accum_mesh, "AccumulatedPositions")
map_BA = intersection_finder(positionsB, positionsA)
call intersector_set_dimension(dim)
! inputs: positionsB, map_BA, positionsA, supermesh_shape
! output: the supermesh!
call recursive_supermesh(positionsA, positionsB, map_BA, supermesh_shape, 1, ele_count(positionsB), accum_positions)
! do ele_B=1,ele_count(positionsB)
! llnode => map_BA(ele_B)%firstnode
! do while(associated(llnode))
! ele_A = llnode%value
! intersection = intersect_elements(positionsA, ele_A, ele_val(positionsB, ele_B), supermesh_shape)
! call unify_meshes(accum_positions, intersection, accum_positions_tmp)
! call deallocate(accum_positions)
! accum_positions = accum_positions_tmp
! llnode => llnode%next
! call deallocate(intersection)
! end do
! end do
! call allocate(accum_positions_tmp, ! need a way to make this continuous!
! call remap_field(accum_positions, accum_positions_tmp)
! call write_triangle_files(trim(output), accum_positions_tmp)
! call deallocate(accum_positions_tmp)
call write_triangle_files(trim(output), accum_positions)
call vtk_write_fields(trim(output), 0, accum_positions, accum_positions%mesh)
contains
recursive subroutine recursive_supermesh(positionsA, positionsB, map_BA, supermesh_shape, start_ele, end_ele, supermesh)
type(vector_field), intent(in) :: positionsA, positionsB
type(ilist), dimension(:), intent(in) :: map_BA
type(element_type), intent(in) :: supermesh_shape
integer, intent(in) :: start_ele, end_ele
type(vector_field), intent(out) :: supermesh
integer, parameter :: blocksize = 4
integer :: i
type(vector_field) :: supermesh_tmp
type(mesh_type) :: supermesh_mesh
type(vector_field) :: supermesh_accum
type(inode), pointer :: llnode
integer :: ele_A, ele_B
integer :: new_start, new_end, step
call allocate(supermesh_mesh, 0, 0, supermesh_shape, "AccumulatedMesh")
call allocate(supermesh, positionsA%dim, supermesh_mesh, "AccumulatedPositions")
call deallocate(supermesh_mesh)
if ((end_ele - start_ele) <= blocksize) then
do ele_B=start_ele,end_ele
llnode => map_BA(ele_B)%firstnode
do while(associated(llnode))
ele_A = llnode%value
supermesh_tmp = intersect_elements(positionsA, ele_A, ele_val(positionsB, ele_B), supermesh_shape)
call unify_meshes_quadratic(supermesh, supermesh_tmp, supermesh_accum)
call deallocate(supermesh)
call deallocate(supermesh_tmp)
supermesh = supermesh_accum
llnode => llnode%next
end do
end do
else
do i=1,blocksize
step = (end_ele - start_ele)/blocksize
new_start = start_ele + (i-1)*step
if (i /= blocksize) then
new_end = start_ele + (i)*step - 1
else
new_end = end_ele ! step might not divide exactly
end if
write(0,*) "Calling recursive_supermesh with element range ", new_start, new_end
call recursive_supermesh(positionsA, positionsB, map_BA, supermesh_shape, new_start, new_end, supermesh_tmp)
call unify_meshes_quadratic(supermesh, supermesh_tmp, supermesh_accum)
call deallocate(supermesh)
call deallocate(supermesh_tmp)
supermesh = supermesh_accum
end do
endif
end subroutine recursive_supermesh
end subroutine unifiedmesh
|