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
|
! 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"
subroutine pseudo_mesh(filename, filename_len, target_elements)
use adapt_state_module
use conformity_measurement
use fields
use fldebug
use global_parameters, only : current_debug_level
use interpolation_module
use limit_metric_module
use read_triangle
use parallel_tools
use reference_counting
use write_triangle
implicit none
integer, intent(in) :: filename_len
character(len = filename_len), intent(in) :: filename
integer, intent(in) :: target_elements
integer :: i
integer, parameter :: adapt_iterations = 5
real, dimension(:, :), allocatable :: unit_matrix
type(tensor_field) :: metric, target_metric
type(vector_field) :: input_positions, output_positions, stage_2_input_positions
ewrite(1, *) "In pseudo_mesh"
if(isparallel()) then
ewrite(-1,*) "You are running in parallel"
FLExit("pseudo_mesh does not work in parallel")
end if
ewrite(2, *) "Target input elements? ", target_elements /= 0
! Step 1: Input
input_positions = read_triangle_files(filename, quad_degree = 4)
! Step 2: Adapt. We do this in two stages - the first to pull us a long way
! from the target mesh, and the second to pull us back towards it. This
! ensures that adaptivity actually does something, rather than just giving us
! a copy of the input mesh.
! Adapt stage 1: Coarsen as much as possible everywhere
ewrite(2, *) "Adapt stage 1: Coarsening"
call allocate(metric, input_positions%mesh, "Metric", field_type = FIELD_TYPE_CONSTANT)
allocate(unit_matrix(metric%dim, metric%dim))
unit_matrix = 0.0
do i = 1, metric%dim
unit_matrix(i, i) = 1.0
end do
call set(metric, unit_matrix)
deallocate(unit_matrix)
call limit_metric(input_positions, metric, target_nodes = 1)
call adapt_mesh(input_positions, metric, output_positions)
call deallocate(metric)
call write_triangle_files("pseudo_mesh_stage_1", output_positions)
ewrite(2, *) "Node count = ", node_count(output_positions)
ewrite(2, *) "Element count = ", ele_count(output_positions)
! Adapt stage 2: Now adapt to the actual target metric
ewrite(2, *) "Adapt stage 2: Targetting input mesh metric"
call compute_mesh_metric(input_positions, target_metric)
! Subcycle, as each time we're targetting an interpolated metric
do i = 1, adapt_iterations
ewrite(2, "(a,i0,a,i0)") "Iteration ", i, " of ", adapt_iterations
call allocate(metric, output_positions%mesh, "InterpolatedTargetMetric")
call linear_interpolation(target_metric, input_positions, metric, output_positions)
if(target_elements /= 0) then
! Target elements, to ensure convergence of repeated pseudo-meshing
call limit_metric_elements(output_positions, metric, target_eles = ele_count(input_positions))
else
! If we're not targetting, limit the metric anyway (but loosly) as it's
! possible that there are very small or large edge lengths at the
! boundaries of the input mesh
call limit_metric(output_positions, metric, min_nodes = max(int(0.5 * node_count(input_positions)), 1), max_nodes = 2 * node_count(input_positions))
end if
! Copy the output from the last adapt, for use as input for the next adapt
call allocate(stage_2_input_positions, output_positions%dim, output_positions%mesh, output_positions%name)
call set(stage_2_input_positions, output_positions)
call deallocate(output_positions)
call adapt_mesh(stage_2_input_positions, metric, output_positions)
call deallocate(metric)
call deallocate(stage_2_input_positions)
call write_triangle_files("pseudo_mesh_stage_2_" // int2str(i), output_positions)
ewrite(2, *) "Node count = ", node_count(output_positions)
ewrite(2, *) "Element count = ", ele_count(output_positions)
end do
call deallocate(target_metric)
call deallocate(input_positions)
! Step 3: Output
call write_triangle_files("pseudo_mesh", output_positions)
call deallocate(output_positions)
call print_references(0)
ewrite(1, *) "Exiting pseudo_mesh"
end subroutine pseudo_mesh
|