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
141
142
143
144
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
257
|
! Copyright (C) 2006 Imperial College London and others.
!
! Please see the AUTHORS file in the main source directory for a full list
! of copyright holders.
!
! Prof. C Pain
! Applied Modelling and Computation Group
! Department of Earth Science and Engineering
! Imperial College London
!
! amcgsoftware@imperial.ac.uk
!
! This library is free software; you can redistribute it and/or
! modify it under the terms of the GNU Lesser General Public
! License as published by the Free Software Foundation; either
! version 2.1 of the License, or (at your option) any later version.
!
! This library is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
! Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with this library; if not, write to the Free Software
! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
! USA
#include "fdebug.h"
module adaptivity_1d
use fields
use fldebug
use hadapt_metric_based_extrude
use node_locking
use tictoc
implicit none
private
public :: adapt_mesh_1d, adaptivity_1d_check_options
contains
subroutine adapt_mesh_1d(old_positions, metric, new_positions, node_ownership, force_preserve_regions)
type(vector_field), intent(in) :: old_positions
type(tensor_field), intent(inout) :: metric
type(vector_field), intent(out) :: new_positions
integer, dimension(:), pointer, optional :: node_ownership
logical, optional, intent(in) :: force_preserve_regions
integer :: i
integer, dimension(:), allocatable :: locked_nodes
logical :: preserve_regions
type(element_type), pointer :: shape
type(scalar_field) :: sizing
ewrite(1, *) "In adapt_mesh_1d"
assert(old_positions%dim == 1)
assert(old_positions%mesh == metric%mesh)
preserve_regions = have_option("/mesh_adaptivity/hr_adaptivity/preserve_mesh_regions") .or. &
& present_and_true(force_preserve_regions)
if(halo_count(old_positions) > 0) then
FLExit("1D adaptivity does not work in parallel")
end if
call get_locked_nodes(old_positions, locked_nodes)
if(size(locked_nodes) > 0) then
FLExit("Node locking is not supported by 1D adaptivity")
end if
deallocate(locked_nodes)
call allocate(sizing, metric%mesh, name = "SizingFunction")
do i = 1, node_count(sizing)
call set(sizing, i, edge_length_from_eigenvalue(node_val(metric, 1, 1, i)))
end do
assert(ele_count(old_positions) > 0)
shape => ele_shape(old_positions, 1)
! since we're using the shape of the old mesh, this shape and its quadrature will survive the adapt
! and should not show up in print_tagged_refences()
shape%refcount%tagged=.false.
shape%quadrature%refcount%tagged=.false.
! TODO: Sort old_positions into descending coordinate order here
assert(descending_coordinate_ordered(old_positions))
call tic(TICTOC_ID_SERIAL_ADAPT)
call adapt_1d(old_positions, sizing, shape, new_positions, preserve_regions=preserve_regions)
call toc(TICTOC_ID_SERIAL_ADAPT)
call deallocate(sizing)
assert(descending_coordinate_ordered(new_positions))
! adapt_1d doesn't build a complete mesh. Build the rest of the mesh.
assert(ele_count(new_positions) == node_count(new_positions) - 1)
do i = 1, ele_count(new_positions)
call set_ele_nodes(new_positions%mesh, i, (/i, i + 1/))
end do
! note that adapt_1d has already allocated and inserted the region_ids
! if they were meant to be preserved
! HOWEVER adapt_1d does assume that the elements as well as the nodes
! are ordered so if this assumption changes here then
! new_positions%mesh%region_ids will have to be reordered too!
assert(surface_element_count(old_positions) == 2)
call add_faces(new_positions%mesh, sndgln = (/1, node_count(new_positions)/), boundary_ids = old_positions%mesh%faces%boundary_ids)
new_positions%name = old_positions%name
new_positions%mesh%name = old_positions%mesh%name
new_positions%option_path = old_positions%option_path
new_positions%mesh%option_path = old_positions%mesh%option_path
if(present(node_ownership)) call generate_1d_node_ownership(old_positions, new_positions, node_ownership)
ewrite(1, *) "Exiting adapt_mesh_1d"
end subroutine adapt_mesh_1d
subroutine generate_1d_node_ownership(old_positions, new_positions, node_ownership)
!!< Generate the 1d node ownership list. Assumes the old and new positions
!!< are descending coordinate ordered.
type(vector_field), intent(in) :: old_positions
type(vector_field), intent(in) :: new_positions
integer, dimension(:), pointer :: node_ownership
integer :: ele, node
assert(descending_coordinate_ordered(old_positions))
assert(descending_coordinate_ordered(new_positions))
assert(node_count(old_positions) > 1)
assert(node_count(new_positions) > 0)
assert(.not. associated(node_ownership))
allocate(node_ownership(node_count(new_positions)))
#ifdef DDEBUG
node_ownership = -1
#endif
! The first node is owned by the first element
node_ownership(1) = 1
ele = 1
new_pos_loop: do node = 2, node_count(new_positions) - 1
do while(node_val(new_positions, 1, node) < node_val(old_positions, 1, ele + 1))
ele = ele + 1
if(ele >= node_count(old_positions)) exit new_pos_loop
end do
! Intermediate nodes are owned by the first element that has a left
! coordinate less than this node coordinate
node_ownership(node) = ele
end do new_pos_loop
! The last node is owned by the last element
node_ownership(node_count(new_positions)) = ele_count(old_positions)
#ifdef DDEBUG
! All nodes are owned by an element
assert(all(node_ownership > 0))
call verify_node_ownership(old_positions, new_positions, node_ownership)
#endif
end subroutine generate_1d_node_ownership
function descending_coordinate_ordered(positions)
!!< Return whether the supplied 1D mesh is in descending coordinate order
type(vector_field), intent(in) :: positions
logical :: descending_coordinate_ordered
integer :: i
assert(positions%dim == 1)
descending_coordinate_ordered = .true.
do i = 2, node_count(positions)
if(node_val(positions, 1, i) > node_val(positions, 1, i - 1)) then
descending_coordinate_ordered = .false.
return
end if
end do
end function descending_coordinate_ordered
subroutine verify_node_ownership(old_positions, new_positions, node_ownership)
!!< Check that the supplied node ownership list is valid. Assumes simplex
!!< elements.
type(vector_field), intent(in) :: old_positions
type(vector_field), intent(in) :: new_positions
integer, dimension(node_count(new_positions)), intent(in) :: node_ownership
integer :: i
real, parameter :: tol = 1000.0 * epsilon(0.0)
do i = 1, node_count(new_positions)
call verify_node_ownership_node(node_ownership(i), i, old_positions, new_positions)
end do
contains
subroutine verify_node_ownership_node(ele, node, old_positions, new_positions)
integer, intent(in) :: ele
integer, intent(in) :: node
type(vector_field), intent(in) :: old_positions
type(vector_field), intent(in) :: new_positions
real, dimension(ele_loc(old_positions, ele)) :: l_coords
l_coords = local_coords(old_positions, ele, node_val(new_positions, node))
if(any(l_coords < -tol)) then
ewrite(-1, "(a,i0,a)") "For node ", node, " in the new positions"
ewrite(-1, *) "Claimed owner in the old positions: ", ele
ewrite(-1, *) "Local coordinates in claimed owner: ", l_coords
ewrite(-1, *) "Test tolerance: ", tol
FLAbort("Invalid node ownership")
end if
end subroutine verify_node_ownership_node
end subroutine verify_node_ownership
subroutine adaptivity_1d_check_options
!!< Checks 1D adaptivity related options
character(len = *), parameter :: base_path = "/mesh_adaptivity/hr_adaptivity"
integer :: dim, stat
if(.not. have_option(base_path)) then
! Nothing to check
return
end if
call get_option("/geometry/dimension", dim, stat)
if(stat /= SPUD_NO_ERROR) then
! This isn't the place to complain about this error
return
else if(have_option(base_path // "/adaptivity_library_adaptivity_1d") .or. dim == 1) then
if(dim /= 1) then
FLExit("1D adaptivity can only be used in 1D")
else if(isparallel()) then
FLExit("1D adaptivity can only be used in serial")
end if
end if
end subroutine adaptivity_1d_check_options
end module adaptivity_1d
|