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
|
#include "fdebug.h"
program mesh_readnadapt
! A small program to read in a triangle mesh, and adapt to an (isotropic)
! edgelength function provided by the user
use read_triangle
use fields
use FEtools
use DGtools
use elements
use sparse_tools
use vtk_interfaces
use transform_elements
use solvers
use petsc_tools
use sparsity_patterns
use signal_vars
use state_module
!use diagnostic_variables
use solvers
use python_state
use global_parameters, only : current_debug_level, PYTHON_FUNC_LEN
use spud
use adapt_integration
use mesh_files
implicit none
type(vector_field), target :: positions, new_positions
integer :: degree, quad_degree
type(quadrature_type), target :: quad
type(element_type), target :: X_shape
type(scalar_field) :: h, new_h
type(tensor_field) :: metric
! Arguments for handling the command line
character(len=256) :: filename
character(len=100) :: fmt,buffer
character(len=PYTHON_FUNC_LEN), save :: func
logical :: file_exists
integer :: dim, loc, nnodes, nelements, node_attributes, unit, io1, node,&
& status
real, dimension(3,3) :: EYE
current_debug_level = 2
ewrite(1,*) 'program '
call Initialize_Petsc()
call python_init
call get_command_argument(1, value=filename, status=status)
select case(status)
case(1:)
print *, 'ERK'
stop
case(:-1)
write(0,*) "Warning: truncating filename"
end select
filename=trim(filename)
ewrite(2,*) 'Getting mesh file information'
call identify_mesh_file(trim(filename), dim, loc, nnodes, nelements, &
node_attributes)
quad_degree = 6
ewrite(2,*) 'Getting quadrature'
quad=make_quadrature(loc=loc, dimension=dim, degree=quad_degree)
ewrite(2,*) 'Getting shape functions'
! Shape functions for positions (linear)
X_shape=make_element_shape(loc=loc, dimension=dim, &
degree=1, quad=quad)
ewrite(2,*) 'reading mesh'
ewrite(2,*) 'loc = ',loc,'dim = ',dim
positions=read_mesh_files(trim(filename), X_shape)
call allocate(h,positions%mesh,'h')
ewrite(2,*) 'Setting h'
ewrite(2,*) 'setting h from python'
inquire(file=trim(filename)//".py",exist=file_exists)
if (.not.file_exists) FLAbort('Couldnt find ' // trim(filename) // '.py file')
unit=free_unit()
open(unit, file=trim(filename)//".py", action="read",&
& status="old")
read(unit, '(a)', end=43) func
! Read all the lines of the file and put in newlines between them.
do
read(unit, '(a)', end=43) buffer
func=trim(func)//achar(10)//trim(buffer)
end do
43 func=trim(func)//achar(10)
close(unit)
ewrite(2,*) func
call set_from_python_function(h,trim(func), positions, 0.0)
call allocate(metric,positions%mesh,name='metric')
call zero(metric)
EYE = 0.0
EYE(1,1) = 1.0
EYE(2,2) = 1.0
EYE(3,3) = 1.0
do node = 1, node_count(h)
call set(metric, node, 1/(h%val(node)**2)*EYE)
end do
call vtk_write_fields(trim(filename)//'_b4', &
index=0, position=positions, &
sfields=(/h/), vfields=(/positions/), & !tfields= (/metric/), &
& model=positions%mesh)
call adapt_mesh(positions, metric, new_positions)
write(0,*) "node_count(positions) == ", node_count(positions)
write(0,*) "node_count(new_positions) == ", node_count(new_positions)
call deallocate(h)
call allocate(h,new_positions%mesh,'h')
call set_from_python_function(h,trim(func), new_positions, 0.0)
call write_mesh_files(trim(filename)//'_out', &
new_positions)
call vtk_write_fields(trim(filename)//'_out', &
index=0, position=new_positions, &
sfields=(/h/), vfields=(/new_positions/),model=new_positions%mesh)
end program mesh_readnadapt
|