Merge lp:~fluidity-core/fluidity/gmsh-on-sphere into lp:fluidity
- gmsh-on-sphere
- Merge into dev-trunk
Status: | Merged | ||||
---|---|---|---|---|---|
Merged at revision: | 3987 | ||||
Proposed branch: | lp:~fluidity-core/fluidity/gmsh-on-sphere | ||||
Merge into: | lp:fluidity | ||||
Diff against target: |
2135 lines (+909/-760) 17 files modified
assemble/Adapt_State.F90 (+4/-5) assemble/Adapt_State_Prescribed.F90 (+2/-2) assemble/tests/test_empty_populate_state.F90 (+1/-0) femtools/Checkpoint.F90 (+12/-12) femtools/Mesh_Files.F90 (+20/-215) femtools/Read_GMSH.F90 (+97/-290) femtools/Read_Triangle.F90 (+1/-163) preprocessor/Populate_State.F90 (+86/-61) tests/spherical_patch/Makefile (+11/-0) tests/spherical_patch/spherical_patch.xml (+47/-0) tests/spherical_patch/spherical_segment.flml (+502/-0) tests/spherical_patch/src/spherical_segment.geo (+113/-0) tools/Checkmesh.F90 (+2/-2) tools/Fladapt.F90 (+4/-4) tools/periodise.F90 (+3/-2) tools/project_to_continuous.F90 (+2/-2) tools/test_laplacian.F90 (+2/-2) |
||||
To merge this branch: | bzr merge lp:~fluidity-core/fluidity/gmsh-on-sphere | ||||
Related bugs: |
|
Reviewer | Review Type | Date Requested | Status |
---|---|---|---|
Jon Hill | Approve | ||
Review via email: mp+99395@code.launchpad.net |
Commit message
Description of the change
This merge fixes the reading of gmsh input meshes in fluidity, where the input mesh is 2d spherical, to be extruded inside fluidity to a 3d mesh. Because this mesh is 2d topologically but has 3d coordinates, this was not supported/broken. The functionality is now tested (first short "on_sphere" test) together with gmsh writing, parallel, checkpointing and flredecomp for this type of mesh.
Some more details of changes in Read_GMSH.F90:
- prevent reading of the gmsh file twice. This was because read_gmsh_simple() called identify_
- only one version of read_gmsh_files to remain is the _simple version (as this was the only one in use). It
first reads in all the data, then establishes the necessary dimensions and then allocates the appropriate femtools
objects.
- the dimension of the read gmsh mesh is now established from the topology of the mesh, basically looking at the
highest dimension elements present, i.e. if there's tets or hexes it's 3d, otherwise 2d (1d gmsh reading is and remains untested). Previously this was done by looking at the z coordinate and seeing if it was small using some hard-coded definition of small=1e-9 (units?!).
- the dimension of the coordinate field that is returned (X%dim) is in most cases chosen to be the same as
the topological dimension. Only with the /geometry/
horizontal mesh (to be extruded inside fluidity) can be read. This was broken before. The other codes that use
embedded meshes (shallow water) use a specific option that only works for triangle, so that functionality should
remain unchanged.
- remove identify_
- because of this remove identify_
- in all subroutines of the mesh_files module make format a required argument. Before, if the argument was not present a "guess" would be made. In some tools where no option tree is present (e.g. test_laplacian) this would always default to "triangle". These therefore now directly call read_triangle_
- identify_
to more processes), where inactive processed needed to know the dimensions of the mesh. This is now communicated
via a MPI_Bcast()
- 3955. By Stephan Kramer
-
Fix for silly unittest.
Jon Hill (jon-hill) wrote : | # |
Lawrence Mitchell (wence) wrote : | # |
The fladapt issue is a real one, but it is generally applicable to all branches, rather than just this one I think. The error (with --enable-debugging) is something like:
Fortran runtime error: Actual string length is shorter than the declared one for dummy argument 'input' (-2123626932/7)
Here's a simplified test of the problematic code:
$ cat > foo.cpp <<\EOF
extern "C" {
#define F77_FUNC(name,NAME) name ## _
#define foo F77_FUNC(foo, FOO)
void foo(const char *in, const int *len);
}
#include <iostream>
#include <stdlib.h>
using namespace std;
int main(int argc, char **argv)
{
string in = "fladapt";
int len = in.size();
foo(in.c_str(), &len);
return 0;
}
EOF
$ cat > foo.f90 <<\EOF
subroutine foo(input, length)
integer, intent(in) :: length
character(len = length), intent(in) :: input
write(*,*)input
end subroutine foo
EOF
$ g++ -fbounds-check -O0 -g -c -o foo_main.o foo.cpp
$ gfortran -fbounds-check -O0 -g -c -o foo.o foo.f90
$ g++ -o foo -O0 -g foo.o foo_main.o -lgfortran
$ gdb ./foo
(gdb) r
Starting program: /home/lmitche4/
At line 1 of file foo.f90
Fortran runtime error: Actual string length is shorter than the declared one for dummy argument 'input' (-8200/7)
Program exited with code 02.
(gdb)
I think TRTD is either specified fixed length strings, or use iso_c_binding appropriately as is done in petsc_readnsolve. I think this is necessary for the majority of C++ main tools in tools/
Jon Hill (jon-hill) wrote : | # |
I fixed a bunch of these a while back. Looks like there's a bunch more to fix...
Will add a bug report and set-up a branch
Stephan Kramer (s-kramer) wrote : | # |
The unittest was already fixed in 3955. Just hasn't run again after that.
Jon Hill (jon-hill) wrote : | # |
Just kicked the buildbot for this to check for greenness and will then approve.
Stephan Kramer (s-kramer) wrote : | # |
Fladapt was fixed by Jon Hill in the trunk and merged in above. Will merge in again from trunk and restart the buildbot queue.
- 3957. By Stephan Kramer
-
Merging in from trunk.
Stephan Kramer (s-kramer) wrote : | # |
Ah, the fladapt failure actually didn't have anything to do with fortran/c binding issues, but was simply caused by my own stupidity. Fix is on the way.
- 3958. By Stephan Kramer
-
Not a good idea to nullify a pointer just before trying to access its pointee.
Jon Hill (jon-hill) wrote : | # |
Kicking buildbot again to test fix
Preview Diff
1 | === modified file 'assemble/Adapt_State.F90' |
2 | --- assemble/Adapt_State.F90 2012-04-16 12:54:50 +0000 |
3 | +++ assemble/Adapt_State.F90 2012-05-16 11:01:22 +0000 |
4 | @@ -73,13 +73,12 @@ |
5 | use diagnostic_variables |
6 | use diagnostic_fields_wrapper_new, only : calculate_diagnostic_variables_new => calculate_diagnostic_variables |
7 | use pickers |
8 | + use write_gmsh |
9 | #ifdef HAVE_ZOLTAN |
10 | use zoltan_integration |
11 | use mpi_interfaces |
12 | #endif |
13 | |
14 | - use mesh_files |
15 | - |
16 | implicit none |
17 | |
18 | private |
19 | @@ -923,10 +922,10 @@ |
20 | output_positions => extract_vector_field(states(1), "Coordinate") |
21 | |
22 | if(isparallel()) then |
23 | - call write_mesh_files(parallel_filename("first_timestep_adapted_mesh"), output_positions) |
24 | + call write_gmsh_file(parallel_filename("first_timestep_adapted_mesh"), output_positions) |
25 | call write_halos("first_timestep_adapted_mesh", output_positions%mesh) |
26 | else |
27 | - call write_mesh_files("first_timestep_adapted_mesh", output_positions) |
28 | + call write_gmsh_file("first_timestep_adapted_mesh", output_positions) |
29 | end if |
30 | |
31 | end if |
32 | @@ -1613,7 +1612,7 @@ |
33 | file_name = adapt_state_debug_file_name("adapted_mesh", mesh_dump_no) |
34 | call find_mesh_to_adapt(states(1), mesh) |
35 | positions = get_coordinate_field(states(1), mesh) |
36 | - call write_mesh_files(file_name, positions) |
37 | + call write_gmsh_file(file_name, positions) |
38 | if(isparallel()) then |
39 | file_name = adapt_state_debug_file_name("adapted_mesh", mesh_dump_no, add_parallel = .false.) ! parallel extension is added by write_halos |
40 | call write_halos(file_name, positions%mesh) |
41 | |
42 | === modified file 'assemble/Adapt_State_Prescribed.F90' |
43 | --- assemble/Adapt_State_Prescribed.F90 2011-07-22 15:34:14 +0000 |
44 | +++ assemble/Adapt_State_Prescribed.F90 2012-05-16 11:01:22 +0000 |
45 | @@ -109,8 +109,8 @@ |
46 | if(have_option(base_path // "/mesh/from_file")) then |
47 | call get_option(base_path // "/mesh/from_file/format/name", format) |
48 | call get_option("/geometry/quadrature/degree", quad_degree) |
49 | - |
50 | - new_positions = read_mesh_files(mesh_name, quad_degree = quad_degree) |
51 | + |
52 | + new_positions = read_mesh_files(mesh_name, quad_degree = quad_degree, format = format) |
53 | else |
54 | ewrite(2, *) "Extracting new mesh from state" |
55 | |
56 | |
57 | === modified file 'assemble/tests/test_empty_populate_state.F90' |
58 | --- assemble/tests/test_empty_populate_state.F90 2009-12-16 14:11:47 +0000 |
59 | +++ assemble/tests/test_empty_populate_state.F90 2012-05-16 11:01:22 +0000 |
60 | @@ -11,6 +11,7 @@ |
61 | logical :: fail |
62 | |
63 | is_active_process = .false. |
64 | + no_active_processes = 0 |
65 | call load_options("data/empty-mesh.flml") |
66 | call populate_state(states) |
67 | |
68 | |
69 | === modified file 'femtools/Checkpoint.F90' |
70 | --- femtools/Checkpoint.F90 2011-07-19 15:52:43 +0000 |
71 | +++ femtools/Checkpoint.F90 2012-05-16 11:01:22 +0000 |
72 | @@ -1,4 +1,4 @@ |
73 | -! Copyright (C) 2006 Imperial College London and others. |
74 | +! Copyrigh (C) 2006 Imperial College London and others. |
75 | ! |
76 | ! Please see the AUTHORS file in the main source directory for a full list |
77 | ! of copyright holders. |
78 | @@ -164,9 +164,9 @@ |
79 | & total_dete_lag |
80 | type(vector_field), pointer :: vfield |
81 | |
82 | - integer(KIND=MPI_OFFSET_KIND) :: location_to_write, offset, offset_to_read |
83 | + integer(KIND=MPI_OFFSET_KIND) :: location_to_write, offset |
84 | integer, ALLOCATABLE, DIMENSION(:) :: status |
85 | - integer :: nints, count, realsize, dimen, total_num_det, number_total_columns |
86 | + integer :: nints, realsize, dimen, total_num_det, number_total_columns |
87 | real, dimension(:), allocatable :: buffer |
88 | |
89 | ewrite(1, *) "Checkpointing detectors" |
90 | @@ -432,14 +432,12 @@ |
91 | logical, optional, intent(in) :: keep_initial_data |
92 | |
93 | type(vector_field), pointer:: position |
94 | - character(len = FIELD_NAME_LEN) :: mesh_name |
95 | + character(len = FIELD_NAME_LEN) :: mesh_name, mesh_format |
96 | character(len = OPTION_PATH_LEN) :: mesh_path, mesh_filename |
97 | integer :: i, n_meshes, stat1, stat2 |
98 | - type(mesh_type), pointer :: mesh |
99 | + type(mesh_type), pointer :: mesh, external_mesh |
100 | logical :: from_file, extruded |
101 | |
102 | - character(len = OPTION_PATH_LEN) :: currentMeshFormat |
103 | - |
104 | assert(len_trim(prefix) > 0) |
105 | |
106 | n_meshes = option_count("/geometry/mesh") |
107 | @@ -467,12 +465,14 @@ |
108 | ! Update the options tree (required for options tree checkpointing) |
109 | if (from_file) then |
110 | call set_option_attribute(trim(mesh_path) // "/from_file/file_name", trim(mesh_filename)) |
111 | + call get_option(trim(mesh_path) // "/from_file/format/name", mesh_format) |
112 | else if (extruded) then |
113 | |
114 | - ! Lunge forth with a wild stab at what the current mesh format is |
115 | - call guess_external_mesh_format(currentMeshFormat) |
116 | + ! the mesh format is determined from the external mesh |
117 | + external_mesh => get_external_mesh(state) |
118 | + call get_option(trim(external_mesh%option_path) // "/from_file/format/name", mesh_format) |
119 | |
120 | - call set_option_attribute(trim(mesh_path) // "/from_mesh/extrude/checkpoint_from_file/format/name", trim(currentMeshFormat), stat=stat1) |
121 | + call set_option_attribute(trim(mesh_path) // "/from_mesh/extrude/checkpoint_from_file/format/name", trim(mesh_format), stat=stat1) |
122 | call set_option_attribute(trim(mesh_path) // "/from_mesh/extrude/checkpoint_from_file/file_name", trim(mesh_filename), stat=stat2) |
123 | if ((stat1/=SPUD_NO_ERROR .and. stat1/=SPUD_NEW_KEY_WARNING) .or. & |
124 | & (stat2/=SPUD_NO_ERROR .and. stat2/=SPUD_NEW_KEY_WARNING)) then |
125 | @@ -487,13 +487,13 @@ |
126 | else |
127 | position => extract_vector_field(state(1), trim(mesh%name)//"Coordinate") |
128 | end if |
129 | - call write_mesh_files(parallel_filename(mesh_filename), position) |
130 | + call write_mesh_files(parallel_filename(mesh_filename), mesh_format, position) |
131 | ! Write out the halos |
132 | ewrite(2, *) "Checkpointing halos" |
133 | call write_halos(mesh_filename, mesh) |
134 | else |
135 | ! Write out the mesh |
136 | - call write_mesh_files(mesh_filename, state(1), mesh) |
137 | + call write_mesh_files(mesh_filename, mesh_format, state(1), mesh) |
138 | end if |
139 | end if |
140 | end do |
141 | |
142 | === modified file 'femtools/Mesh_Files.F90' |
143 | --- femtools/Mesh_Files.F90 2011-03-08 13:26:40 +0000 |
144 | +++ femtools/Mesh_Files.F90 2012-05-16 11:01:22 +0000 |
145 | @@ -62,8 +62,7 @@ |
146 | private |
147 | |
148 | interface read_mesh_files |
149 | - module procedure read_mesh_files_to_field, & |
150 | - read_mesh_files_to_state, read_mesh_simple |
151 | + module procedure read_mesh_simple |
152 | end interface |
153 | |
154 | interface write_mesh_files |
155 | @@ -71,8 +70,7 @@ |
156 | write_positions_to_file |
157 | end interface |
158 | |
159 | - public :: read_mesh_files, identify_mesh_file, write_mesh_files |
160 | - public :: guess_external_mesh_format |
161 | + public :: read_mesh_files, write_mesh_files |
162 | |
163 | |
164 | contains |
165 | @@ -82,143 +80,18 @@ |
166 | ! Read routines first |
167 | ! -------------------------------------------------------------------------- |
168 | |
169 | - subroutine identify_mesh_file(filename, dim, loc, nodes, elements, & |
170 | - node_attributes, selements, selement_boundaries, format) |
171 | - ! Discover the dimension and size of the mesh. Filename is |
172 | - ! the base name of the mesh file. |
173 | - ! In parallel, filename must *include* the process number. |
174 | - |
175 | - character(len=*), intent(in) :: filename |
176 | - character(len=*), optional, intent(in) :: format |
177 | - !! Dimension of mesh elements. |
178 | - integer, intent(out), optional :: dim |
179 | - !! Number of vertices of elements. |
180 | - integer, intent(out), optional :: loc |
181 | - !! Node and element counts. |
182 | - integer, intent(out), optional :: nodes, elements |
183 | - integer, intent(out), optional :: node_attributes |
184 | - ! Surface element meta data |
185 | - integer, optional, intent(out) :: selements |
186 | - integer, optional, intent(out) :: selement_boundaries |
187 | - |
188 | - character(len=OPTION_PATH_LEN) :: meshFormat |
189 | - |
190 | - if( .not. present(format) ) then |
191 | - call guess_external_mesh_format(meshFormat) |
192 | - else |
193 | - meshFormat = format |
194 | - end if |
195 | - |
196 | - ! Call appropriate subroutine depending upon format |
197 | - select case( trim(meshFormat) ) |
198 | - case("triangle") |
199 | - call identify_triangle_file(filename, dim, loc, nodes,& |
200 | - elements, node_attributes, selements, selement_boundaries) |
201 | - |
202 | - case("gmsh") |
203 | - call identify_gmsh_file(filename, dim, loc, nodes,& |
204 | - elements, node_attributes, selements, selement_boundaries) |
205 | - |
206 | - ! Additional mesh format subroutines go here |
207 | - |
208 | - case default |
209 | - FLExit("Identifying mesh type "//format//" not supported within Fluidity") |
210 | - end select |
211 | - |
212 | - end subroutine identify_mesh_file |
213 | - |
214 | - |
215 | - ! -------------------------------------------------------------------------- |
216 | - |
217 | - |
218 | - function read_mesh_files_to_field(filename, format, shape) result (field) |
219 | - character(len=*), intent(in) :: filename |
220 | - type(element_type), intent(in), target :: shape |
221 | - type(vector_field) :: field |
222 | - character(len=OPTION_PATH_LEN) :: meshFormat |
223 | - |
224 | - character(len=*), optional, intent(in) :: format |
225 | - |
226 | - if( .not. present(format) ) then |
227 | - call guess_external_mesh_format(meshFormat) |
228 | - else |
229 | - meshFormat = format |
230 | - end if |
231 | - |
232 | - select case( trim(meshFormat) ) |
233 | - case("triangle") |
234 | - field = read_triangle_files(filename, shape) |
235 | - |
236 | - case("gmsh") |
237 | - field = read_gmsh_file(filename, shape) |
238 | - |
239 | - ! Additional mesh format subroutines go here |
240 | - |
241 | - case default |
242 | - FLExit("Reading mesh type "//format//" not supported within Fluidity") |
243 | - end select |
244 | - end function read_mesh_files_to_field |
245 | - |
246 | - |
247 | - ! -------------------------------------------------------------------------- |
248 | - |
249 | - function read_mesh_files_to_state(filename, shape, shape_type, & |
250 | - n_states, format) & |
251 | - result (result_state) |
252 | - |
253 | - ! Filename is the base name of the mesh file without .ele, .msh, etc. |
254 | - ! In parallel the filename must *not* include the process number. |
255 | - |
256 | - character(len=*), intent(in) :: filename |
257 | - type(element_type), intent(in), target :: shape |
258 | - logical , intent(in):: shape_type |
259 | - integer, intent(in), optional :: n_states |
260 | - |
261 | - type(state_type) :: result_state |
262 | - |
263 | - character(len=*), optional, intent(in) :: format |
264 | - character(len=option_path_len) :: meshFormat |
265 | - |
266 | - if( .not. present(format) ) then |
267 | - call guess_external_mesh_format(meshFormat) |
268 | - else |
269 | - meshFormat = format |
270 | - end if |
271 | - |
272 | - |
273 | - select case( trim(meshFormat) ) |
274 | - case("triangle") |
275 | - result_state = read_triangle_files(filename, shape,shape_type,n_states) |
276 | - |
277 | - case("gmsh") |
278 | - result_state = read_gmsh_file(filename, shape,shape_type,n_states) |
279 | - |
280 | - ! Additional mesh format subroutines go here |
281 | - |
282 | - case default |
283 | - FLExit("Reading mesh type "//format//" not supported within Fluidity") |
284 | - end select |
285 | - |
286 | - end function read_mesh_files_to_state |
287 | - |
288 | - |
289 | - ! -------------------------------------------------------------------------- |
290 | - |
291 | - function read_mesh_simple(filename, quad_degree, & |
292 | - quad_ngi, no_faces, & |
293 | - quad_family, mdim, format ) & |
294 | + function read_mesh_simple(filename, format, quad_degree, & |
295 | + quad_ngi, quad_family, mdim) & |
296 | result (field) |
297 | |
298 | ! A simpler mechanism for reading a mesh file into a field. |
299 | ! In parallel the filename must *not* include the process number. |
300 | |
301 | - character(len=*), intent(in) :: filename |
302 | + character(len=*), intent(in) :: filename, format |
303 | ! The degree of the quadrature. |
304 | integer, intent(in), optional, target :: quad_degree |
305 | ! The degree of the quadrature. |
306 | integer, intent(in), optional, target :: quad_ngi |
307 | - ! Whether to add_faces on the resulting mesh. |
308 | - logical, intent(in), optional :: no_faces |
309 | ! What quadrature family to use |
310 | integer, intent(in), optional :: quad_family |
311 | ! Dimension of mesh |
312 | @@ -226,24 +99,17 @@ |
313 | |
314 | type(vector_field) :: field |
315 | |
316 | - character(len=*), optional, intent(in) :: format |
317 | - character(len=option_path_len) :: meshFormat |
318 | - |
319 | - if( .not. present(format) ) then |
320 | - call guess_external_mesh_format(meshFormat) |
321 | - else |
322 | - meshFormat = format |
323 | - end if |
324 | - |
325 | - |
326 | - select case( trim(meshFormat) ) |
327 | + select case( trim(format) ) |
328 | case("triangle") |
329 | - field = read_triangle_files(filename, quad_degree, quad_ngi, & |
330 | - no_faces, quad_family, mdim) |
331 | + field = read_triangle_files(filename, quad_degree=quad_degree, quad_ngi=quad_ngi, & |
332 | + quad_family=quad_family, mdim=mdim) |
333 | |
334 | case("gmsh") |
335 | - field = read_gmsh_file(filename, quad_degree, quad_ngi, & |
336 | - no_faces, quad_family) |
337 | + if (present(mdim)) then |
338 | + FLExit("Cannot specify dimension for gmsh format") |
339 | + end if |
340 | + field = read_gmsh_file(filename, quad_degree=quad_degree, quad_ngi=quad_ngi, & |
341 | + quad_family=quad_family) |
342 | |
343 | ! Additional mesh format subroutines go here |
344 | |
345 | @@ -260,22 +126,15 @@ |
346 | ! -------------------------------------------------------------------------- |
347 | |
348 | |
349 | - subroutine write_mesh_to_file(filename, state, mesh, format ) |
350 | + subroutine write_mesh_to_file(filename, format, state, mesh) |
351 | ! Write out the supplied mesh to the specified filename as mesh files. |
352 | |
353 | character(len = *), intent(in) :: filename |
354 | - character(len = *), intent(in), optional ::format |
355 | - character(len = option_path_len) :: meshFormat |
356 | + character(len = *), intent(in) :: format |
357 | type(state_type), intent(in) :: state |
358 | type(mesh_type), intent(in) :: mesh |
359 | |
360 | - if( .not. present(format) ) then |
361 | - call guess_external_mesh_format(meshFormat) |
362 | - else |
363 | - meshFormat = format |
364 | - end if |
365 | - |
366 | - select case(meshFormat) |
367 | + select case(format) |
368 | case("triangle") |
369 | call write_triangle_files(filename, state, mesh) |
370 | |
371 | @@ -293,24 +152,13 @@ |
372 | |
373 | ! -------------------------------------------------------------------------- |
374 | |
375 | - subroutine write_positions_to_file(filename, positions, format) |
376 | + subroutine write_positions_to_file(filename, format, positions) |
377 | !!< Write out the mesh given by the position field in mesh files |
378 | !!< In parallel, empty trailing processes are not written. |
379 | - character(len=*), intent(in):: filename |
380 | + character(len=*), intent(in):: filename, format |
381 | type(vector_field), intent(in):: positions |
382 | - !! "triangle" or "gmsh" |
383 | - character(len=*), intent(in), optional :: format |
384 | - |
385 | - character(len=option_path_len) :: meshFormat |
386 | - |
387 | - if( .not. present(format) ) then |
388 | - call guess_external_mesh_format(meshFormat) |
389 | - else |
390 | - meshFormat = format |
391 | - end if |
392 | - |
393 | - |
394 | - select case( trim(meshFormat) ) |
395 | + |
396 | + select case( trim(format) ) |
397 | case("triangle") |
398 | call write_triangle_files( trim(filename), positions) |
399 | |
400 | @@ -325,48 +173,5 @@ |
401 | |
402 | end subroutine write_positions_to_file |
403 | |
404 | - |
405 | - |
406 | - |
407 | - ! -------------------------------------------------------------------------- |
408 | - ! Subroutine which finds external mesh format and puts it in 'meshFormat' |
409 | - ! Follows similar logic to Field_Options::get_external_mesh() |
410 | - |
411 | - |
412 | - subroutine guess_external_mesh_format(meshFormat) |
413 | - character(len=*) :: meshFormat |
414 | - character(len=OPTION_PATH_LEN) :: meshPath, formatPath |
415 | - integer :: numMeshes, meshFound, i |
416 | - |
417 | - numMeshes = option_count("/geometry/mesh") |
418 | - |
419 | - ! Search for external meshes, and exit loop when found one |
420 | - meshFound=0 |
421 | - do i = 0, numMeshes-1 |
422 | - meshPath = "/geometry/mesh["//int2str(i)//"]/from_file" |
423 | - |
424 | - if(have_option(trim(meshPath))) then |
425 | - meshFound=1 |
426 | - exit |
427 | - end if |
428 | - end do |
429 | - |
430 | - ! We've found a mesh, now get its format |
431 | - if (meshFound==1) then |
432 | - formatPath = trim(meshPath) // "/format/name" |
433 | - |
434 | - if(have_option(trim(formatPath))) then |
435 | - call get_option( trim(formatPath), meshFormat ) |
436 | - else |
437 | - ! If mesh exists but no format found, default to triangle |
438 | - meshFormat = "triangle" |
439 | - end if |
440 | - else |
441 | - ! If we can't find any external meshes... |
442 | - meshFormat = "triangle" |
443 | - end if |
444 | - |
445 | - end subroutine guess_external_mesh_format |
446 | - |
447 | end module mesh_files |
448 | |
449 | |
450 | === modified file 'femtools/Read_GMSH.F90' |
451 | --- femtools/Read_GMSH.F90 2011-10-21 10:05:25 +0000 |
452 | +++ femtools/Read_GMSH.F90 2012-05-16 11:01:22 +0000 |
453 | @@ -44,160 +44,46 @@ |
454 | private |
455 | |
456 | interface read_gmsh_file |
457 | - module procedure read_gmsh_file_to_field, read_gmsh_simple, & |
458 | - read_gmsh_file_to_state |
459 | + module procedure read_gmsh_simple |
460 | end interface |
461 | |
462 | - public :: read_gmsh_file, identify_gmsh_file |
463 | + public :: read_gmsh_file |
464 | + |
465 | + integer, parameter:: GMSH_LINE=1, GMSH_TRIANGLE=2, GMSH_QUAD=3, GMSH_TET=4, GMSH_HEX=5, GMSH_NODE=15 |
466 | |
467 | contains |
468 | |
469 | ! ----------------------------------------------------------------- |
470 | - ! GMSH version of triangle equivalent. |
471 | - |
472 | - subroutine identify_gmsh_file(filename, numDimenOut, locOut, & |
473 | - numNodesOut, numElementsOut, & |
474 | - nodeAttributesOut, selementsOut, boundaryFlagOut) |
475 | - ! Discover the dimension and size of the GMSH mesh. |
476 | - ! Filename is the base name of the file without .msh. |
477 | - ! In parallel, filename must *include* the process number. |
478 | - |
479 | - character(len=*), intent(in) :: filename |
480 | - character(len=option_path_len) :: lfilename |
481 | - |
482 | - !! Number of vertices of elements. |
483 | - integer, intent(out), optional :: numDimenOut, locOut, numNodesOut |
484 | - integer, intent(out), optional :: numElementsOut, nodeAttributesOut |
485 | - integer, intent(out), optional :: selementsOut, boundaryFlagOut |
486 | - |
487 | - |
488 | - logical :: fileExists |
489 | - |
490 | - integer :: fd, gmshFormat |
491 | - type(GMSHnode), pointer :: nodes(:) |
492 | - type(GMSHelement), pointer :: elements(:), faces(:) |
493 | - integer :: numElements, boundaryFlag, numNodes, numDimen |
494 | - integer :: loc, effDimen, nodeAttributes |
495 | - integer :: i |
496 | - |
497 | - logical :: haveBounds, haveInternalBounds |
498 | - |
499 | - |
500 | - lfilename = trim(filename) // ".msh" |
501 | - |
502 | - ! Read node file header |
503 | - inquire(file = trim(lfilename), exist = fileExists) |
504 | - if(.not. fileExists) then |
505 | - FLExit("gmsh file " // trim(lfilename) // " not found") |
506 | - end if |
507 | - |
508 | - ewrite(2, *) "Opening " // trim(lfilename) // " for reading." |
509 | - fd = free_unit() |
510 | - open(unit = fd, file = trim(lfilename), & |
511 | - err=42, action="read", access="stream", form="formatted" ) |
512 | - |
513 | - call read_header(fd, lfilename, gmshFormat) |
514 | - |
515 | - ! Read in the nodes |
516 | - call read_nodes_coords( fd, lfilename, gmshFormat, nodes ) |
517 | - |
518 | - numDimen = mesh_dimensions( nodes ) |
519 | - numNodes = size(nodes) |
520 | - |
521 | - ! Read in elements |
522 | - call read_faces_and_elements( fd, lfilename, gmshFormat, & |
523 | - elements, faces ) |
524 | - |
525 | - |
526 | - ! Try reading in node column ID data (if there is any) |
527 | - call read_node_column_IDs( fd, lfilename, gmshFormat, nodes) |
528 | - |
529 | - close( fd ) |
530 | - |
531 | - ! We're assuming all elements have the same number of vertices/nodes |
532 | - loc = size(elements(1)%nodeIDs) |
533 | - |
534 | - do i=1, numNodes |
535 | - if( nodes(i)%columnID>0 ) nodeAttributes=1 |
536 | - end do |
537 | - |
538 | - ! NOTE: 'boundaries' variable |
539 | - ! (see Read_Triangle.F90) is a flag which indicates whether |
540 | - ! faces have boundaries (physical IDs). Values: |
541 | - ! =0 : no boundaries |
542 | - ! =1 : boundaries, normal |
543 | - ! =2 : boundaries, periodic mesh (internal boundary) |
544 | - |
545 | - haveBounds=.false. |
546 | - haveInternalBounds=.false. |
547 | - |
548 | - do i=1, size(faces) |
549 | - if(faces(i)%numTags > 0) then |
550 | - ! We have boundaries, |
551 | - haveBounds=.true. |
552 | - |
553 | - ! We have internal boundaries (at the join of a period mesh) |
554 | - if(faces(i)%numTags >= 4) then |
555 | - if (faces(i)%tags(4) > 0) then |
556 | - haveInternalBounds=.true. |
557 | - end if |
558 | - end if |
559 | - end if |
560 | - end do |
561 | - |
562 | - if(haveInternalBounds) then |
563 | - boundaryFlag=2 |
564 | - elseif(haveBounds) then |
565 | - boundaryFlag=1 |
566 | - else |
567 | - boundaryFlag=0 |
568 | - end if |
569 | - |
570 | - if( numDimen.eq.2 .and. have_option("/geometry/spherical_earth/") ) then |
571 | - effDimen = numDimen+1 |
572 | - else |
573 | - effDimen = numDimen |
574 | - end if |
575 | - |
576 | - |
577 | - ! Return optional variables requested |
578 | - |
579 | - if(present(nodeAttributesOut)) nodeAttributesOut=nodeAttributes |
580 | - if(present(numDimenOut)) numDimenOut=effDimen |
581 | - if(present(numElementsOut)) numElementsOut=numElements |
582 | - if(present(locOut)) locOut=loc |
583 | - if(present(boundaryFlagOut)) boundaryFlagOut=boundaryFlag |
584 | - |
585 | - deallocate(nodes) |
586 | - |
587 | - return |
588 | - |
589 | -42 FLExit("Unable to open "//trim(lfilename)) |
590 | - |
591 | - end subroutine identify_gmsh_file |
592 | - |
593 | - |
594 | - |
595 | - ! ----------------------------------------------------------------- |
596 | ! The main function for reading GMSH files |
597 | |
598 | - function read_gmsh_file_to_field(filename, shape) result (field) |
599 | - ! Filename is the base name of the GMSH file without .msh |
600 | - ! In parallel the filename must *not* include the process number. |
601 | + function read_gmsh_simple( filename, quad_degree, & |
602 | + quad_ngi, quad_family ) & |
603 | + result (field) |
604 | + !!< Read a GMSH file into a coordinate field. |
605 | + !!< In parallel the filename must *not* include the process number. |
606 | |
607 | character(len=*), intent(in) :: filename |
608 | - type(element_type), intent(in), target :: shape |
609 | - type(vector_field) :: field |
610 | + !! The degree of the quadrature. |
611 | + integer, intent(in), optional, target :: quad_degree |
612 | + !! The degree of the quadrature. |
613 | + integer, intent(in), optional, target :: quad_ngi |
614 | + !! What quadrature family to use |
615 | + integer, intent(in), optional :: quad_family |
616 | + !! result: a coordinate field |
617 | + type(vector_field) :: field |
618 | + |
619 | + type(quadrature_type):: quad |
620 | + type(element_type):: shape |
621 | + type(mesh_type):: mesh |
622 | |
623 | integer :: fd |
624 | integer, pointer, dimension(:) :: sndglno, boundaryIDs, faceOwner |
625 | |
626 | character(len = parallel_filename_len(filename)) :: lfilename |
627 | integer :: loc, sloc |
628 | - type(mesh_type) :: mesh |
629 | integer :: numNodes, numElements, numFaces |
630 | logical :: haveBounds, haveElementOwners, haveRegionIDs |
631 | - integer :: numDimen, effDimen |
632 | + integer :: dim, coordinate_dim |
633 | integer :: gmshFormat |
634 | integer :: n, d, e, f, nodeID |
635 | |
636 | @@ -225,11 +111,9 @@ |
637 | ! Read in the nodes |
638 | call read_nodes_coords( fd, lfilename, gmshFormat, nodes ) |
639 | |
640 | - numDimen = mesh_dimensions( nodes ) |
641 | - |
642 | ! Read in elements |
643 | call read_faces_and_elements( fd, lfilename, gmshFormat, & |
644 | - elements, faces ) |
645 | + elements, faces, dim) |
646 | |
647 | call read_node_column_IDs( fd, lfilename, gmshFormat, nodes ) |
648 | |
649 | @@ -285,39 +169,51 @@ |
650 | |
651 | end if |
652 | |
653 | - if( numDimen.eq.2 .and. have_option("/geometry/spherical_earth/") ) then |
654 | - effDimen = numDimen+1 |
655 | - else |
656 | - effDimen = numDimen |
657 | - end if |
658 | - |
659 | - |
660 | + if(have_option("/geometry/spherical_earth/") ) then |
661 | + ! on the sphere the input mesh may be 2d (extrusion), or 3d but |
662 | + ! Coordinate is always 3-dimensional |
663 | + coordinate_dim = 3 |
664 | + else |
665 | + coordinate_dim = dim |
666 | + end if |
667 | + |
668 | + loc = size( elements(1)%nodeIDs ) |
669 | + if (numFaces>0) then |
670 | + sloc = size( faces(1)%nodeIDs ) |
671 | + else |
672 | + sloc = 0 |
673 | + end if |
674 | + |
675 | + ! Now construct within Fluidity data structures |
676 | + |
677 | + if (present(quad_degree)) then |
678 | + quad = make_quadrature(loc, dim, degree=quad_degree, family=quad_family) |
679 | + else if (present(quad_ngi)) then |
680 | + quad = make_quadrature(loc, dim, ngi=quad_ngi, family=quad_family) |
681 | + else |
682 | + FLAbort("Need to specify either quadrature degree or ngi") |
683 | + end if |
684 | + shape=make_element_shape(loc, dim, 1, quad) |
685 | call allocate(mesh, numNodes, numElements, shape, name="CoordinateMesh") |
686 | - call allocate( field, effDimen, mesh, name="Coordinate") |
687 | - call deallocate( mesh ) |
688 | - |
689 | - ! Now construct within Fluidity data structures |
690 | + call allocate( field, coordinate_dim, mesh, name="Coordinate") |
691 | + |
692 | + ! deallocate our references of mesh, shape and quadrature: |
693 | + call deallocate(mesh) |
694 | + call deallocate(shape) |
695 | + call deallocate(quad) |
696 | + |
697 | |
698 | if (haveRegionIDs) then |
699 | allocate( field%mesh%region_ids(numElements) ) |
700 | end if |
701 | - if(nodes(1)%columnID.ge.0) allocate(field%mesh%columns(1:numNodes)) |
702 | - |
703 | - loc = size( elements(1)%nodeIDs ) |
704 | - if (numFaces>0) then |
705 | - sloc = size( faces(1)%nodeIDs ) |
706 | - else |
707 | - sloc = 0 |
708 | - end if |
709 | - |
710 | - assert(loc==shape%loc) |
711 | + if(nodes(1)%columnID>=0) allocate(field%mesh%columns(1:numNodes)) |
712 | |
713 | ! Loop round nodes copying across coords and column IDs to field mesh, |
714 | ! if they exist |
715 | do n=1, numNodes |
716 | |
717 | nodeID = nodes(n)%nodeID |
718 | - forall (d = 1:effDimen) |
719 | + forall (d = 1:field%dim) |
720 | field%val(d,nodeID) = nodes(n)%x(d) |
721 | end forall |
722 | |
723 | @@ -380,86 +276,8 @@ |
724 | |
725 | 43 FLExit("Unable to open "//trim(lfilename)) |
726 | |
727 | - end function read_gmsh_file_to_field |
728 | - |
729 | - |
730 | - |
731 | - |
732 | - |
733 | - |
734 | - ! ----------------------------------------------------------------- |
735 | - ! Simplified interface to reading GMSH files |
736 | - function read_gmsh_simple( filename, quad_degree, & |
737 | - quad_ngi, no_faces, quad_family ) & |
738 | - result (field) |
739 | - !!< A simpler mechanism for reading a GMSH file into a field. |
740 | - !!< In parallel the filename must *not* include the process number. |
741 | - |
742 | - character(len=*), intent(in) :: filename |
743 | - !! The degree of the quadrature. |
744 | - integer, intent(in), optional, target :: quad_degree |
745 | - !! The degree of the quadrature. |
746 | - integer, intent(in), optional, target :: quad_ngi |
747 | - !! Whether to add_faces on the resulting mesh. |
748 | - logical, intent(in), optional :: no_faces |
749 | - !! What quadrature family to use |
750 | - integer, intent(in), optional :: quad_family |
751 | - |
752 | - type(vector_field) :: field |
753 | - type(quadrature_type) :: quad |
754 | - type(element_type) :: shape |
755 | - |
756 | - integer :: dim, loc |
757 | - |
758 | - |
759 | - if(isparallel()) then |
760 | - call identify_gmsh_file(parallel_filename(filename), dim, loc) |
761 | - else |
762 | - call identify_gmsh_file(filename, dim, loc) |
763 | - end if |
764 | - |
765 | - |
766 | - if (present(quad_degree)) then |
767 | - quad = make_quadrature(loc, dim, degree=quad_degree, family=quad_family) |
768 | - else if (present(quad_ngi)) then |
769 | - quad = make_quadrature(loc, dim, ngi=quad_ngi, family=quad_family) |
770 | - else |
771 | - FLAbort("Need to specify either quadrature degree or ngi") |
772 | - end if |
773 | - |
774 | - |
775 | - shape=make_element_shape(loc, dim, 1, quad) |
776 | - |
777 | - field=read_gmsh_file(filename, shape) |
778 | - |
779 | - ! deallocate our references of shape and quadrature: |
780 | - call deallocate_element(shape) |
781 | - call deallocate(quad) |
782 | - |
783 | end function read_gmsh_simple |
784 | |
785 | - |
786 | - |
787 | - |
788 | - ! ----------------------------------------------------------------- |
789 | - ! Read GMSH file to state object. |
790 | - |
791 | - function read_gmsh_file_to_state(filename, shape,shape_type,n_states) & |
792 | - result (result_state) |
793 | - ! Filename is the base name of the GMSH file without .node or .ele. |
794 | - ! In parallel the filename must *not* include the process number. |
795 | - |
796 | - character(len=*), intent(in) :: filename |
797 | - type(element_type), intent(in), target :: shape |
798 | - logical , intent(in):: shape_type |
799 | - integer, intent(in), optional :: n_states |
800 | - type(state_type) :: result_state |
801 | - |
802 | - FLAbort("read_gmsh_file_to_state() not implemented yet") |
803 | - |
804 | - end function read_gmsh_file_to_state |
805 | - |
806 | - |
807 | ! ----------------------------------------------------------------- |
808 | ! Read through the head to decide whether binary or ASCII, and decide |
809 | ! whether this looks like a GMSH mesh file or not. |
810 | @@ -666,37 +484,17 @@ |
811 | |
812 | |
813 | ! ----------------------------------------------------------------- |
814 | - ! Guesstimate mesh dimensions. If any node z-coords are non-zero, then |
815 | - ! it's a 3D mesh; otherwise 2D. Not pretty, but there seems no other |
816 | - ! way to do this. |
817 | - |
818 | - integer function mesh_dimensions( nodes ) |
819 | - type(GMSHnode), pointer :: nodes(:) |
820 | - real, pointer :: absZ(:) |
821 | - |
822 | - mesh_dimensions = 3 |
823 | - |
824 | - allocate ( absZ(size(nodes)) ) |
825 | - absZ(:) = nodes(:)%x(3) |
826 | - |
827 | - if( all( absZ > -verySmall .and. absZ < verySmall) ) mesh_dimensions=2 |
828 | - |
829 | - deallocate(absZ) |
830 | - end function mesh_dimensions |
831 | - |
832 | - |
833 | - |
834 | - |
835 | - ! ----------------------------------------------------------------- |
836 | - ! Read in element header data |
837 | + ! Read in element header data and establish topological dimension |
838 | |
839 | subroutine read_faces_and_elements( fd, filename, gmshFormat, & |
840 | - elements, faces ) |
841 | - |
842 | - integer :: fd, gmshFormat |
843 | - character(len=*) :: filename |
844 | - |
845 | - type(GMSHelement), pointer :: allElements(:), elements(:), faces(:) |
846 | + elements, faces, dim) |
847 | + |
848 | + integer, intent(in) :: fd, gmshFormat |
849 | + character(len=*), intent(in) :: filename |
850 | + type(GMSHelement), pointer :: elements(:), faces(:) |
851 | + integer, intent(out) :: dim |
852 | + |
853 | + type(GMSHelement), pointer :: allElements(:) |
854 | |
855 | integer :: numAllElements |
856 | character(len=longStringLen) :: charBuf |
857 | @@ -708,14 +506,14 @@ |
858 | |
859 | |
860 | read(fd,*) charBuf |
861 | - if( trim(charBuf) .ne. "$Elements" ) then |
862 | + if( trim(charBuf)/="$Elements" ) then |
863 | FLExit("Error: cannot find '$Elements' in GMSH mesh file") |
864 | end if |
865 | |
866 | read(fd,*) numAllElements |
867 | |
868 | ! Sanity check. |
869 | - if(numAllElements .lt. 1) then |
870 | + if(numAllElements<1) then |
871 | FLExit("Error: number of elements in GMSH file < 1") |
872 | end if |
873 | |
874 | @@ -799,23 +597,21 @@ |
875 | allElements(e)%type ) |
876 | |
877 | select case ( allElements(e)%type ) |
878 | - ! A line |
879 | - case (1) |
880 | + case (GMSH_LINE) |
881 | numEdges = numEdges+1 |
882 | - ! A triangle |
883 | - case (2) |
884 | + case (GMSH_TRIANGLE) |
885 | numTriangles = numTriangles+1 |
886 | - ! A quad |
887 | - case (3) |
888 | + case (GMSH_QUAD) |
889 | numQuads = numQuads+1 |
890 | - case (4) |
891 | + case (GMSH_TET) |
892 | numTets = numTets+1 |
893 | - case (5) |
894 | + case (GMSH_HEX) |
895 | numHexes = numHexes+1 |
896 | - case (15) |
897 | + case (GMSH_NODE) |
898 | ! Do nothing |
899 | case default |
900 | - FLExit("read_faces_and_elements(): unsupported element type") |
901 | + ewrite(0,*) "element id,type: ", allElements(e)%elementID, allElements(e)%type |
902 | + FLExit("Unsupported element type in gmsh .msh file") |
903 | end select |
904 | |
905 | end do |
906 | @@ -833,29 +629,40 @@ |
907 | ! meshes are verboten: |
908 | ! tet/hex, tet/quad, triangle/hex and triangle/quad |
909 | |
910 | - if (numTets .gt. 0) then |
911 | + if (numTets>0) then |
912 | numElements = numTets |
913 | - elementType = 4 |
914 | + elementType = GMSH_TET |
915 | numFaces = numTriangles |
916 | - faceType = 2 |
917 | + faceType = GMSH_TRIANGLE |
918 | + dim = 3 |
919 | + if (numQuads>0 .or. numHexes>0) then |
920 | + FLExit("Cannot combine hexes or quads with tetrahedrals in one gmsh .msh file") |
921 | + end if |
922 | |
923 | - elseif (numTriangles .gt. 0) then |
924 | + elseif (numTriangles>0) then |
925 | numElements = numTriangles |
926 | - elementType = 2 |
927 | + elementType = GMSH_TRIANGLE |
928 | numFaces = numEdges |
929 | - faceType = 1 |
930 | + faceType = GMSH_LINE |
931 | + dim = 2 |
932 | + if (numQuads>0 .or. numHexes>0) then |
933 | + FLExit("Cannot combine hexes or quads with triangles in one gmsh .msh file") |
934 | + end if |
935 | |
936 | elseif (numHexes .gt. 0) then |
937 | numElements = numHexes |
938 | - elementType = 5 |
939 | + elementType = GMSH_HEX |
940 | numFaces = numQuads |
941 | - faceType = 3 |
942 | + faceType = GMSH_QUAD |
943 | + dim = 3 |
944 | |
945 | elseif (numQuads .gt. 0) then |
946 | numElements = numQuads |
947 | - elementType = 3 |
948 | + elementType = GMSH_QUAD |
949 | numFaces = numEdges |
950 | - faceType = 1 |
951 | + faceType = GMSH_LINE |
952 | + dim = 2 |
953 | + |
954 | else |
955 | FLExit("Unsupported mixture of face/element types") |
956 | end if |
957 | |
958 | === modified file 'femtools/Read_Triangle.F90' |
959 | --- femtools/Read_Triangle.F90 2011-03-08 13:26:40 +0000 |
960 | +++ femtools/Read_Triangle.F90 2012-05-16 11:01:22 +0000 |
961 | @@ -42,8 +42,7 @@ |
962 | private |
963 | |
964 | interface read_triangle_files |
965 | - module procedure read_triangle_files_to_field, & |
966 | - read_triangle_files_to_state, read_triangle_simple |
967 | + module procedure read_triangle_files_to_field, read_triangle_simple |
968 | end interface |
969 | |
970 | public :: read_triangle_files, identify_triangle_file, read_elemental_mappings, read_triangle_serial |
971 | @@ -414,167 +413,6 @@ |
972 | |
973 | end function read_triangle_files_to_field |
974 | |
975 | - function read_triangle_files_to_state(filename, shape,shape_type,n_states) result (result_state) |
976 | - !!< Filename is the base name of the triangle file without .node or .ele. |
977 | - !!< In parallel the filename must *not* include the process number. |
978 | - |
979 | - character(len=*), intent(in) :: filename |
980 | - type(element_type), intent(in), target :: shape |
981 | - logical , intent(in):: shape_type |
982 | - integer, intent(in), optional :: n_states |
983 | - type(state_type) :: result_state |
984 | - type(vector_field) :: field |
985 | - type(scalar_field), allocatable, dimension(:) :: attribs |
986 | - type(scalar_field) :: bndry_mark |
987 | - |
988 | - character(len = parallel_filename_len(filename)) :: lfilename |
989 | - integer :: node_unit, ele_unit |
990 | - real, allocatable, dimension(:) :: read_buffer |
991 | - character(len=6) :: write_buffer |
992 | - integer :: i, j, nodes, dim, xdim, node_attributes, boundaries,& |
993 | - & ele_attributes, loc, elements |
994 | - integer, allocatable, dimension(:) :: node_order |
995 | - |
996 | - type(mesh_type), dimension(:), allocatable, save :: meshes |
997 | - integer, save :: counter=0 |
998 | - |
999 | - assert(shape_type) |
1000 | - |
1001 | - ! If running in parallel, add the process number |
1002 | - if(isparallel()) then |
1003 | - lfilename = parallel_filename(filename) |
1004 | - else |
1005 | - lfilename = trim(filename) |
1006 | - end if |
1007 | - |
1008 | - node_unit=free_unit() |
1009 | - |
1010 | - counter=counter+1 |
1011 | - if (counter==1) then |
1012 | - if (present(n_states)) then |
1013 | - allocate(meshes(n_states)) |
1014 | - else |
1015 | - allocate(meshes(1)) |
1016 | - end if |
1017 | - end if |
1018 | - assert(counter .le. size(meshes)) |
1019 | - |
1020 | - call nullify(result_state) |
1021 | - call deallocate(result_state) |
1022 | - |
1023 | - ewrite(2, *) "Opening "//trim(lfilename)//".node for reading." |
1024 | - ! Open node file |
1025 | - open(unit=node_unit, file=trim(lfilename)//".node", err=42, action="read") |
1026 | - |
1027 | - ! Read node file header. |
1028 | - read (node_unit, *) nodes, xdim, node_attributes, boundaries |
1029 | - |
1030 | - ele_unit=free_unit() |
1031 | - |
1032 | - ewrite(2, *) "Opening "//trim(lfilename)//".ele for reading." |
1033 | - ! Open element file |
1034 | - open(unit=ele_unit, file=trim(lfilename)//".ele", err=43, action="read") |
1035 | - |
1036 | - ! Read element file header. |
1037 | - read (ele_unit, *) elements, loc, ele_attributes |
1038 | - |
1039 | - assert(loc==shape%loc) |
1040 | - allocate(node_order(loc)) |
1041 | - select case(loc) |
1042 | - case(3) |
1043 | - node_order = (/1,2,3/) |
1044 | - case(6) |
1045 | - node_order = (/1,6,2,5,4,3/) |
1046 | - case default |
1047 | - do j=1,loc |
1048 | - node_order(:)=j |
1049 | - enddo |
1050 | - end select |
1051 | - |
1052 | - call allocate(meshes(counter), nodes, elements, shape, name="CoordinateMesh") |
1053 | - call allocate(field, xdim, meshes(counter), name=filename//' Coordinate') |
1054 | - allocate(attribs(node_attributes)) |
1055 | - do j=1,node_attributes |
1056 | - write(write_buffer,'(i0)') j |
1057 | - call allocate(attribs(j),meshes(counter),& |
1058 | - name='Attribute '//trim(write_buffer)) |
1059 | - end do |
1060 | - call allocate(bndry_mark,meshes(counter),& |
1061 | - name='Boundary Marker') |
1062 | - |
1063 | - allocate(read_buffer(xdim+node_attributes+boundaries+1)) |
1064 | - |
1065 | - do i=1,nodes |
1066 | - read(node_unit,*) read_buffer |
1067 | - |
1068 | - forall (j=1:xdim) |
1069 | - field%val(j,i)=read_buffer(j+1) |
1070 | - end forall |
1071 | - |
1072 | - forall (j=1:node_attributes) |
1073 | - attribs(j)%val(i)=read_buffer(j+1+xdim) |
1074 | - end forall |
1075 | - |
1076 | - if (boundaries>0) then |
1077 | - bndry_mark%val(i)=read_buffer(1+1+xdim+node_attributes) |
1078 | - end if |
1079 | - |
1080 | - end do |
1081 | - |
1082 | - |
1083 | - deallocate(read_buffer) |
1084 | - allocate(read_buffer(loc+ele_attributes+1)) |
1085 | - |
1086 | - if(ele_attributes==1) then ! this assumes the element attribute is a region id |
1087 | - allocate(field%mesh%region_ids(1:elements)) |
1088 | - do j=1,node_attributes |
1089 | - allocate(attribs(j)%mesh%region_ids(1:elements)) |
1090 | - end do |
1091 | - allocate(bndry_mark%mesh%region_ids(1:elements)) |
1092 | - allocate(meshes(counter)%region_ids(1:elements)) |
1093 | - end if |
1094 | - |
1095 | - do i=1,elements |
1096 | - read(ele_unit,*) read_buffer |
1097 | - |
1098 | - |
1099 | - field%mesh%ndglno((i-1)*loc+1:i*loc)=read_buffer(node_order+1) |
1100 | - forall (j=1:node_attributes) |
1101 | - attribs(j)%mesh%ndglno((i-1)*loc+1:i*loc)=read_buffer(node_order+1) |
1102 | - end forall |
1103 | - bndry_mark%mesh%ndglno((i-1)*loc+1:i*loc)=read_buffer(node_order+1) |
1104 | - meshes(counter)%ndglno((i-1)*loc+1:i*loc)=read_buffer(node_order+1) |
1105 | - |
1106 | - if(ele_attributes==1) then |
1107 | - field%mesh%region_ids(i)=read_buffer(loc+2) |
1108 | - forall (j=1:node_attributes) |
1109 | - attribs(j)%mesh%region_ids(i)=read_buffer(loc+2) |
1110 | - end forall |
1111 | - bndry_mark%mesh%region_ids(i)=read_buffer(loc+2) |
1112 | - meshes(counter)%region_ids(i)=read_buffer(loc+2) |
1113 | - end if |
1114 | - end do |
1115 | - |
1116 | - close(node_unit) |
1117 | - close(ele_unit) |
1118 | - |
1119 | - |
1120 | - call insert(result_state,field,"Coordinate") |
1121 | - do j=1,node_attributes |
1122 | - write(write_buffer,'(i0)') j |
1123 | - call insert(result_state,attribs(j),"Attributes "//trim(write_buffer)) |
1124 | - end do |
1125 | - if(boundaries>0) call insert(result_state,bndry_mark,"Boundary Marker") |
1126 | - call insert(result_state,meshes(counter),filename) |
1127 | - |
1128 | - return |
1129 | - |
1130 | -42 FLExit("Unable to open "//trim(lfilename)//".node") |
1131 | - |
1132 | -43 FLExit("Unable to open "//trim(lfilename)//".ele") |
1133 | - |
1134 | - end function read_triangle_files_to_state |
1135 | - |
1136 | function read_triangle_simple(filename, quad_degree, quad_ngi, no_faces, quad_family, mdim) result (field) |
1137 | !!< A simpler mechanism for reading a triangle file into a field. |
1138 | !!< In parallel the filename must *not* include the process number. |
1139 | |
1140 | === modified file 'preprocessor/Populate_State.F90' |
1141 | --- preprocessor/Populate_State.F90 2012-01-30 15:36:17 +0000 |
1142 | +++ preprocessor/Populate_State.F90 2012-05-16 11:01:22 +0000 |
1143 | @@ -174,6 +174,7 @@ |
1144 | character(len=OPTION_PATH_LEN) :: mesh_path, mesh_file_name,& |
1145 | & mesh_file_format, from_file_path |
1146 | integer, dimension(:), pointer :: coplanar_ids |
1147 | + integer, dimension(3) :: mesh_dims |
1148 | integer :: i, j, nmeshes, nstates, quad_degree, stat |
1149 | type(element_type), pointer :: shape |
1150 | type(quadrature_type), pointer :: quad |
1151 | @@ -219,22 +220,96 @@ |
1152 | call get_option("/geometry/quadrature/degree", quad_degree) |
1153 | quad_family = get_quad_family() |
1154 | |
1155 | + |
1156 | + if (is_active_process) then |
1157 | + select case (mesh_file_format) |
1158 | + case ("triangle", "gmsh") |
1159 | + ! Get mesh dimension if present |
1160 | + call get_option(trim(mesh_path)//"/from_file/dimension", mdim, stat) |
1161 | + ! Read mesh |
1162 | + if(stat==0) then |
1163 | + position=read_mesh_files(trim(mesh_file_name), & |
1164 | + quad_degree=quad_degree, & |
1165 | + quad_family=quad_family, mdim=mdim, & |
1166 | + format=mesh_file_format) |
1167 | + else |
1168 | + position=read_mesh_files(trim(mesh_file_name), & |
1169 | + quad_degree=quad_degree, & |
1170 | + quad_family=quad_family, & |
1171 | + format=mesh_file_format) |
1172 | + end if |
1173 | + mesh=position%mesh |
1174 | + case ("vtu") |
1175 | + position_ptr => vtk_cache_read_positions_field(mesh_file_name) |
1176 | + ! No hybrid mesh support here |
1177 | + assert(ele_count(position_ptr) > 0) |
1178 | + dim = position_ptr%dim |
1179 | + loc = ele_loc(position_ptr, 1) |
1180 | + |
1181 | + ! Generate a copy, and swap the quadrature degree |
1182 | + ! Note: Even if positions_ptr has the correct quadrature degree, it |
1183 | + ! won't have any faces and hence a copy is still required (as |
1184 | + ! add_faces is a construction routine only) |
1185 | + allocate(quad) |
1186 | + allocate(shape) |
1187 | + quad = make_quadrature(loc, dim, degree = quad_degree, family=quad_family) |
1188 | + shape = make_element_shape(loc, dim, 1, quad) |
1189 | + call allocate(mesh, nodes = node_count(position_ptr), elements = ele_count(position_ptr), shape = shape, name = position_ptr%mesh%name) |
1190 | + do j = 1, ele_count(mesh) |
1191 | + call set_ele_nodes(mesh, j, ele_nodes(position_ptr%mesh, j)) |
1192 | + end do |
1193 | + call add_faces(mesh) |
1194 | + call allocate(position, dim, mesh, position_ptr%name) |
1195 | + call set(position, position_ptr) |
1196 | + call deallocate(mesh) |
1197 | + call deallocate(shape) |
1198 | + call deallocate(quad) |
1199 | + deallocate(quad) |
1200 | + deallocate(shape) |
1201 | + |
1202 | + mesh = position%mesh |
1203 | + case default |
1204 | + ewrite(-1,*) trim(mesh_file_format), " is not a valid format for a mesh file" |
1205 | + FLAbort("Invalid format for mesh file") |
1206 | + end select |
1207 | + end if |
1208 | + |
1209 | + if (no_active_processes /= getnprocs()) then |
1210 | + ! not all processes are active, they need to be told the mesh dimensions |
1211 | + |
1212 | + ! receive the mesh dimension from rank 0 |
1213 | + if (getrank()==0) then |
1214 | + if (is_active_process) then |
1215 | + ! normally rank 0 should always be active, so it knows the dimensions |
1216 | + mesh_dims(1)=mesh_dim(mesh) |
1217 | + mesh_dims(2)=ele_loc(mesh,1) |
1218 | + if (associated(mesh%columns)) then |
1219 | + mesh_dims(3)=1 |
1220 | + else |
1221 | + mesh_dims(3)=0 |
1222 | + end if |
1223 | + else |
1224 | + ! this is a special case for a unit test with 1 inactive process |
1225 | + call get_option('/geometry/dimension', mesh_dims(1)) |
1226 | + mesh_dims(2)=mesh_dims(1)+1 |
1227 | + mesh_dims(3)=0 |
1228 | + end if |
1229 | + end if |
1230 | + call MPI_bcast(mesh_dims, 3, getpinteger(), 0, MPI_COMM_FEMTOOLS, stat) |
1231 | + end if |
1232 | + |
1233 | + |
1234 | if (.not. is_active_process) then |
1235 | ! is_active_process records whether we have data on disk or not |
1236 | ! see the comment in Global_Parameters. In this block, |
1237 | ! we want to allocate an empty mesh and positions. |
1238 | |
1239 | + dim=mesh_dims(1) |
1240 | + loc=mesh_dims(2) |
1241 | + column_ids=mesh_dims(3) |
1242 | + |
1243 | allocate(quad) |
1244 | allocate(shape) |
1245 | - if (no_active_processes == 1) then |
1246 | - call identify_mesh_file(trim(mesh_file_name), dim, loc, & |
1247 | - node_attributes=column_ids, & |
1248 | - format=mesh_file_format ) |
1249 | - else |
1250 | - call identify_mesh_file(trim(mesh_file_name) // "_0", dim, loc, & |
1251 | - node_attributes=column_ids, & |
1252 | - format=mesh_file_format ) |
1253 | - end if |
1254 | quad = make_quadrature(loc, dim, degree=quad_degree, family=quad_family) |
1255 | shape=make_element_shape(loc, dim, 1, quad) |
1256 | call allocate(mesh, nodes=0, elements=0, shape=shape, name="EmptyMesh") |
1257 | @@ -252,58 +327,8 @@ |
1258 | |
1259 | deallocate(quad) |
1260 | deallocate(shape) |
1261 | - |
1262 | - else if(trim(mesh_file_format)=="triangle" & |
1263 | - .or. trim(mesh_file_format)=="gmsh") then |
1264 | - ! Get mesh dimension if present |
1265 | - call get_option(trim(mesh_path)//"/from_file/dimension", mdim, stat) |
1266 | - ! Read mesh |
1267 | - if(stat==0) then |
1268 | - position=read_mesh_files(trim(mesh_file_name), & |
1269 | - quad_degree=quad_degree, & |
1270 | - quad_family=quad_family, mdim=mdim, & |
1271 | - format=mesh_file_format) |
1272 | - else |
1273 | - position=read_mesh_files(trim(mesh_file_name), & |
1274 | - quad_degree=quad_degree, & |
1275 | - quad_family=quad_family, & |
1276 | - format=mesh_file_format) |
1277 | - end if |
1278 | - mesh=position%mesh |
1279 | - else if(trim(mesh_file_format) == "vtu") then |
1280 | - position_ptr => vtk_cache_read_positions_field(mesh_file_name) |
1281 | - ! No hybrid mesh support here |
1282 | - assert(ele_count(position_ptr) > 0) |
1283 | - dim = position_ptr%dim |
1284 | - loc = ele_loc(position_ptr, 1) |
1285 | - |
1286 | - ! Generate a copy, and swap the quadrature degree |
1287 | - ! Note: Even if positions_ptr has the correct quadrature degree, it |
1288 | - ! won't have any faces and hence a copy is still required (as |
1289 | - ! add_faces is a construction routine only) |
1290 | - allocate(quad) |
1291 | - allocate(shape) |
1292 | - quad = make_quadrature(loc, dim, degree = quad_degree, family=quad_family) |
1293 | - shape = make_element_shape(loc, dim, 1, quad) |
1294 | - call allocate(mesh, nodes = node_count(position_ptr), elements = ele_count(position_ptr), shape = shape, name = position_ptr%mesh%name) |
1295 | - do j = 1, ele_count(mesh) |
1296 | - call set_ele_nodes(mesh, j, ele_nodes(position_ptr%mesh, j)) |
1297 | - end do |
1298 | - call add_faces(mesh) |
1299 | - call allocate(position, dim, mesh, position_ptr%name) |
1300 | - call set(position, position_ptr) |
1301 | - call deallocate(mesh) |
1302 | - call deallocate(shape) |
1303 | - call deallocate(quad) |
1304 | - deallocate(quad) |
1305 | - deallocate(shape) |
1306 | - |
1307 | - mesh = position%mesh |
1308 | - else |
1309 | - ewrite(-1,*) trim(mesh_file_format), " is not a valid format for a mesh file" |
1310 | - FLAbort("Invalid format for mesh file") |
1311 | - end if |
1312 | - |
1313 | + end if |
1314 | + |
1315 | ! if there is a derived mesh which specifies periodic bcs |
1316 | ! to be *removed*, we assume the external mesh is periodic |
1317 | mesh%periodic = option_count("/geometry/mesh/from_mesh/& |
1318 | |
1319 | === added directory 'tests/spherical_patch' |
1320 | === added file 'tests/spherical_patch/Makefile' |
1321 | --- tests/spherical_patch/Makefile 1970-01-01 00:00:00 +0000 |
1322 | +++ tests/spherical_patch/Makefile 2012-05-16 11:01:22 +0000 |
1323 | @@ -0,0 +1,11 @@ |
1324 | +GMSH_GEO_FILE=src/spherical_segment.geo |
1325 | +NPROCS=3 |
1326 | +input: clean |
1327 | + gmsh -2 -bin $(GMSH_GEO_FILE) -o patch.msh |
1328 | + ../../bin/fldecomp -n $(NPROCS) -s -m gmsh patch |
1329 | + |
1330 | +clean: |
1331 | + rm -rf *.msh *.halo *.vtu *.stat fluidity.* |
1332 | + rm -rf spherical_segment_* *.pvtu |
1333 | + rm -rf flredecomp.* matrixdump* |
1334 | + |
1335 | |
1336 | === added file 'tests/spherical_patch/spherical_patch.xml' |
1337 | --- tests/spherical_patch/spherical_patch.xml 1970-01-01 00:00:00 +0000 |
1338 | +++ tests/spherical_patch/spherical_patch.xml 2012-05-16 11:01:22 +0000 |
1339 | @@ -0,0 +1,47 @@ |
1340 | +<?xml version='1.0' encoding='utf-8'?> |
1341 | +<testproblem> |
1342 | + <name>spherical patch<comment>Tests gmsh reading and writing on the sphere with 2d horizontal spherical input mesh, checkpointing and flredecomp.</comment></name> |
1343 | + <owner userid="skramer"/> |
1344 | + <problem_definition length="short" nprocs="1"> |
1345 | + <command_line>mpirun -np 3 fluidity -v2 -l spherical_segment.flml && |
1346 | +mpirun -np 3 ../../bin/flredecomp -v -l -i 3 -o 1 spherical_segment_2_checkpoint serial && |
1347 | +fluidity -v2 -l serial.flml</command_line> |
1348 | + </problem_definition> |
1349 | + <variables> |
1350 | + <variable name="volume_before" language="python">from fluidity_tools import stat_parser |
1351 | +stat=stat_parser('spherical_segment.stat') |
1352 | +# this only works because we don't move the free surface |
1353 | +# we actually compute H \int_top \eta, i.e. the volume difference with |
1354 | +# equilibrium state multiplied by depth H |
1355 | +volume_before=stat['Fields']['FreeSurface']['integral']</variable> |
1356 | + <variable name="volume_after" language="python">from fluidity_tools import stat_parser |
1357 | +stat=stat_parser('spherical_segment_checkpoint.stat') |
1358 | +# this only works because we don't move the free surface |
1359 | +# we actually compute H \int_top \eta, i.e. the volume difference with |
1360 | +# equilibrium state multiplied by depth H |
1361 | +volume_after=stat['Fields']['FreeSurface']['integral']</variable> |
1362 | + <variable name="solvers_converged" language="python">import os |
1363 | +files = os.listdir("./") |
1364 | +solvers_converged = not "matrixdump" in files and not "matrixdump.info" in files</variable> |
1365 | + <variable name="max_velocity_before" language="python">from fluidity_tools import stat_parser |
1366 | +stat=stat_parser('spherical_segment.stat') |
1367 | +max_velocity_before=stat['Fields']['Velocity%magnitude']['max']</variable> |
1368 | + <variable name="max_velocity_after" language="python">from fluidity_tools import stat_parser |
1369 | +stat=stat_parser('spherical_segment_checkpoint.stat') |
1370 | +max_velocity_after=stat['Fields']['Velocity%magnitude']['max']</variable> |
1371 | + </variables> |
1372 | + <pass_tests> |
1373 | + <test name="volume_conservation_run1" language="python">from fluidity_tools import compare_variable |
1374 | +compare_variable(volume_before[0], volume_before[-1], 1e-5)</test> |
1375 | + <test name="volume_conservation_run12" language="python">from fluidity_tools import compare_variable |
1376 | +# compares volume at the end of first run with that of end first timestep second run |
1377 | +compare_variable(volume_before[-1], volume_after[0], 1e-5)</test> |
1378 | + <test name="volume_conservation_run2" language="python">from fluidity_tools import compare_variable |
1379 | +compare_variable(volume_after[0], volume_after[-1], 1e-5)</test> |
1380 | + <test name="solvers_converged" language="python">assert solvers_converged</test> |
1381 | + <test name="max_velocity_run1" language="python"># make sure we're actually doing something and not blowing up |
1382 | +assert all(max_velocity_before>0.01) and all(max_velocity_before<0.02)</test> |
1383 | + <test name="max_velocity_run2" language="python"># make sure we're actually doing something and not blowing up |
1384 | +assert all(max_velocity_after>0.01) and all(max_velocity_after<0.02)</test> |
1385 | + </pass_tests> |
1386 | +</testproblem> |
1387 | |
1388 | === added file 'tests/spherical_patch/spherical_segment.flml' |
1389 | --- tests/spherical_patch/spherical_segment.flml 1970-01-01 00:00:00 +0000 |
1390 | +++ tests/spherical_patch/spherical_segment.flml 2012-05-16 11:01:22 +0000 |
1391 | @@ -0,0 +1,502 @@ |
1392 | +<?xml version='1.0' encoding='utf-8'?> |
1393 | +<fluidity_options> |
1394 | + <simulation_name> |
1395 | + <string_value lines="1">spherical_segment</string_value> |
1396 | + </simulation_name> |
1397 | + <problem_type> |
1398 | + <string_value lines="1">oceans</string_value> |
1399 | + </problem_type> |
1400 | + <geometry> |
1401 | + <dimension> |
1402 | + <integer_value rank="0">3</integer_value> |
1403 | + </dimension> |
1404 | + <mesh name="CoordinateMesh"> |
1405 | + <from_mesh> |
1406 | + <mesh name="BaseMesh"/> |
1407 | + <mesh_shape> |
1408 | + <polynomial_degree> |
1409 | + <integer_value rank="0">2</integer_value> |
1410 | + </polynomial_degree> |
1411 | + </mesh_shape> |
1412 | + <stat> |
1413 | + <include_in_stat/> |
1414 | + </stat> |
1415 | + </from_mesh> |
1416 | + </mesh> |
1417 | + <mesh name="VelocityMesh"> |
1418 | + <from_mesh> |
1419 | + <mesh name="BaseMesh"/> |
1420 | + <mesh_continuity> |
1421 | + <string_value>discontinuous</string_value> |
1422 | + </mesh_continuity> |
1423 | + <stat> |
1424 | + <exclude_from_stat/> |
1425 | + </stat> |
1426 | + </from_mesh> |
1427 | + </mesh> |
1428 | + <mesh name="PressureMesh"> |
1429 | + <from_mesh> |
1430 | + <mesh name="BaseMesh"/> |
1431 | + <mesh_shape> |
1432 | + <polynomial_degree> |
1433 | + <integer_value rank="0">2</integer_value> |
1434 | + </polynomial_degree> |
1435 | + </mesh_shape> |
1436 | + <stat> |
1437 | + <exclude_from_stat/> |
1438 | + </stat> |
1439 | + </from_mesh> |
1440 | + </mesh> |
1441 | + <mesh name="InputMesh"> |
1442 | + <from_file file_name="patch"> |
1443 | + <format name="gmsh"/> |
1444 | + <stat> |
1445 | + <exclude_from_stat/> |
1446 | + </stat> |
1447 | + </from_file> |
1448 | + </mesh> |
1449 | + <mesh name="BaseMesh"> |
1450 | + <from_mesh> |
1451 | + <mesh name="InputMesh"/> |
1452 | + <extrude> |
1453 | + <regions name="WholeMesh"> |
1454 | + <bottom_depth> |
1455 | + <constant> |
1456 | + <real_value rank="0">90</real_value> |
1457 | + </constant> |
1458 | + </bottom_depth> |
1459 | + <sizing_function> |
1460 | + <constant> |
1461 | + <real_value rank="0">30</real_value> |
1462 | + </constant> |
1463 | + </sizing_function> |
1464 | + <top_surface_id> |
1465 | + <integer_value rank="0">1</integer_value> |
1466 | + </top_surface_id> |
1467 | + <bottom_surface_id> |
1468 | + <integer_value rank="0">2</integer_value> |
1469 | + </bottom_surface_id> |
1470 | + </regions> |
1471 | + </extrude> |
1472 | + <stat> |
1473 | + <include_in_stat/> |
1474 | + </stat> |
1475 | + </from_mesh> |
1476 | + </mesh> |
1477 | + <mesh name="P0Mesh"> |
1478 | + <from_mesh> |
1479 | + <mesh name="BaseMesh"/> |
1480 | + <mesh_shape> |
1481 | + <polynomial_degree> |
1482 | + <integer_value rank="0">0</integer_value> |
1483 | + </polynomial_degree> |
1484 | + </mesh_shape> |
1485 | + <mesh_continuity> |
1486 | + <string_value>discontinuous</string_value> |
1487 | + </mesh_continuity> |
1488 | + <stat> |
1489 | + <exclude_from_stat/> |
1490 | + </stat> |
1491 | + </from_mesh> |
1492 | + </mesh> |
1493 | + <quadrature> |
1494 | + <degree> |
1495 | + <integer_value rank="0">4</integer_value> |
1496 | + </degree> |
1497 | + <surface_degree> |
1498 | + <integer_value rank="0">3</integer_value> |
1499 | + </surface_degree> |
1500 | + </quadrature> |
1501 | + <spherical_earth> |
1502 | + <superparametric_mapping/> |
1503 | + </spherical_earth> |
1504 | + <ocean_boundaries> |
1505 | + <top_surface_ids> |
1506 | + <integer_value shape="1" rank="1">1</integer_value> |
1507 | + </top_surface_ids> |
1508 | + <bottom_surface_ids> |
1509 | + <integer_value shape="1" rank="1">2</integer_value> |
1510 | + </bottom_surface_ids> |
1511 | + <scalar_field name="DistanceToTop" rank="0"> |
1512 | + <diagnostic> |
1513 | + <algorithm name="Internal" material_phase_support="multiple"/> |
1514 | + <mesh name="CoordinateMesh"/> |
1515 | + <output/> |
1516 | + <stat/> |
1517 | + <convergence> |
1518 | + <include_in_convergence/> |
1519 | + </convergence> |
1520 | + <detectors> |
1521 | + <include_in_detectors/> |
1522 | + </detectors> |
1523 | + <steady_state> |
1524 | + <include_in_steady_state/> |
1525 | + </steady_state> |
1526 | + </diagnostic> |
1527 | + </scalar_field> |
1528 | + <scalar_field name="DistanceToBottom" rank="0"> |
1529 | + <diagnostic> |
1530 | + <algorithm name="Internal" material_phase_support="multiple"/> |
1531 | + <mesh name="CoordinateMesh"/> |
1532 | + <output/> |
1533 | + <stat/> |
1534 | + <convergence> |
1535 | + <include_in_convergence/> |
1536 | + </convergence> |
1537 | + <detectors> |
1538 | + <include_in_detectors/> |
1539 | + </detectors> |
1540 | + <steady_state> |
1541 | + <include_in_steady_state/> |
1542 | + </steady_state> |
1543 | + </diagnostic> |
1544 | + </scalar_field> |
1545 | + </ocean_boundaries> |
1546 | + </geometry> |
1547 | + <io> |
1548 | + <dump_format> |
1549 | + <string_value>vtk</string_value> |
1550 | + </dump_format> |
1551 | + <dump_period> |
1552 | + <constant> |
1553 | + <real_value rank="0">0.0</real_value> |
1554 | + <comment>43200: dump results every 12 hours (0.5 days) of simulated time. Given a finish-time of 2592000 s (1 month), 60 dumps are expected. Given a time-step size of 300s, a dump every 144 time steps is expected.</comment> |
1555 | + </constant> |
1556 | + </dump_period> |
1557 | + <output_mesh name="VelocityMesh"/> |
1558 | + <checkpointing> |
1559 | + <checkpoint_period_in_dumps> |
1560 | + <integer_value rank="0">10</integer_value> |
1561 | + </checkpoint_period_in_dumps> |
1562 | + <checkpoint_at_end/> |
1563 | + </checkpointing> |
1564 | + <stat/> |
1565 | + </io> |
1566 | + <timestepping> |
1567 | + <current_time> |
1568 | + <real_value rank="0">0.0</real_value> |
1569 | + <time_units date="seconds since 1987-01-05 00:00.0"/> |
1570 | + </current_time> |
1571 | + <timestep> |
1572 | + <real_value rank="0">3600.0</real_value> |
1573 | + </timestep> |
1574 | + <finish_time> |
1575 | + <real_value rank="0">2592000</real_value> |
1576 | + <comment>(2592000)1 month, in seconds. Given a timestep size of 300 s, the simulation should perform 8640 time-steps.</comment> |
1577 | + </finish_time> |
1578 | + <final_timestep> |
1579 | + <integer_value rank="0">2</integer_value> |
1580 | + </final_timestep> |
1581 | + <nonlinear_iterations> |
1582 | + <integer_value rank="0">2</integer_value> |
1583 | + </nonlinear_iterations> |
1584 | + </timestepping> |
1585 | + <physical_parameters> |
1586 | + <gravity> |
1587 | + <magnitude> |
1588 | + <real_value rank="0">9.81</real_value> |
1589 | + </magnitude> |
1590 | + <vector_field name="GravityDirection" rank="1"> |
1591 | + <prescribed> |
1592 | + <mesh name="CoordinateMesh"/> |
1593 | + <value name="WholeMesh"> |
1594 | + <python> |
1595 | + <string_value lines="20" type="code" language="python">def val(X, t): |
1596 | + |
1597 | + a = X[0] |
1598 | + b = X[1] |
1599 | + c = X[2] |
1600 | + |
1601 | + x_component = -(a/((a**2+b**2+c**2)**0.5)) |
1602 | + y_component = -(b/((a**2+b**2+c**2)**0.5)) |
1603 | + z_component = -(c/((a**2+b**2+c**2)**0.5)) |
1604 | + |
1605 | + return [x_component, y_component, z_component]</string_value> |
1606 | + </python> |
1607 | + </value> |
1608 | + <output/> |
1609 | + <stat> |
1610 | + <include_in_stat/> |
1611 | + </stat> |
1612 | + <detectors> |
1613 | + <exclude_from_detectors/> |
1614 | + </detectors> |
1615 | + </prescribed> |
1616 | + </vector_field> |
1617 | + </gravity> |
1618 | + <coriolis> |
1619 | + <on_sphere> |
1620 | + <omega> |
1621 | + <real_value rank="0">7.27220522e-5</real_value> |
1622 | + </omega> |
1623 | + </on_sphere> |
1624 | + </coriolis> |
1625 | + </physical_parameters> |
1626 | + <material_phase name="Fields"> |
1627 | + <equation_of_state> |
1628 | + <fluids> |
1629 | + <linear> |
1630 | + <reference_density> |
1631 | + <real_value rank="0">1.0</real_value> |
1632 | + </reference_density> |
1633 | + <subtract_out_hydrostatic_level/> |
1634 | + </linear> |
1635 | + </fluids> |
1636 | + </equation_of_state> |
1637 | + <scalar_field name="Pressure" rank="0"> |
1638 | + <prognostic> |
1639 | + <mesh name="PressureMesh"/> |
1640 | + <spatial_discretisation> |
1641 | + <continuous_galerkin> |
1642 | + <remove_stabilisation_term/> |
1643 | + <integrate_continuity_by_parts/> |
1644 | + </continuous_galerkin> |
1645 | + </spatial_discretisation> |
1646 | + <scheme> |
1647 | + <poisson_pressure_solution> |
1648 | + <string_value lines="1">only first timestep</string_value> |
1649 | + </poisson_pressure_solution> |
1650 | + <use_projection_method/> |
1651 | + </scheme> |
1652 | + <solver> |
1653 | + <iterative_method name="cg"/> |
1654 | + <preconditioner name="mg"> |
1655 | + <vertical_lumping/> |
1656 | + </preconditioner> |
1657 | + <relative_error> |
1658 | + <real_value rank="0">1.0e-8</real_value> |
1659 | + </relative_error> |
1660 | + <max_iterations> |
1661 | + <integer_value rank="0">5000</integer_value> |
1662 | + </max_iterations> |
1663 | + <never_ignore_solver_failures/> |
1664 | + <diagnostics> |
1665 | + <monitors/> |
1666 | + </diagnostics> |
1667 | + </solver> |
1668 | + <initial_condition name="WholeMesh"> |
1669 | + <free_surface> |
1670 | + <python> |
1671 | + <string_value type="code" lines="20" language="python">def val(X,t): |
1672 | + from math import asin, atan2, sqrt, pi, exp |
1673 | + r=sqrt(sum([x**2 for x in X])) |
1674 | + lon=180/pi*atan2(X[1],X[0]) |
1675 | + lat=180/pi*asin(X[2]/r) |
1676 | + d=(lat-52.75)**2+(lon+5.875)**2 |
1677 | + return exp(-d)</string_value> |
1678 | + </python> |
1679 | + </free_surface> |
1680 | + </initial_condition> |
1681 | + <output/> |
1682 | + <stat/> |
1683 | + <convergence> |
1684 | + <include_in_convergence/> |
1685 | + </convergence> |
1686 | + <detectors> |
1687 | + <exclude_from_detectors/> |
1688 | + </detectors> |
1689 | + <steady_state> |
1690 | + <include_in_steady_state/> |
1691 | + </steady_state> |
1692 | + <consistent_interpolation/> |
1693 | + </prognostic> |
1694 | + </scalar_field> |
1695 | + <scalar_field name="Density" rank="0"> |
1696 | + <diagnostic> |
1697 | + <algorithm name="Internal" material_phase_support="multiple"/> |
1698 | + <mesh name="VelocityMesh"/> |
1699 | + <output/> |
1700 | + <stat/> |
1701 | + <convergence> |
1702 | + <include_in_convergence/> |
1703 | + </convergence> |
1704 | + <detectors> |
1705 | + <include_in_detectors/> |
1706 | + </detectors> |
1707 | + <steady_state> |
1708 | + <include_in_steady_state/> |
1709 | + </steady_state> |
1710 | + </diagnostic> |
1711 | + </scalar_field> |
1712 | + <vector_field name="Velocity" rank="1"> |
1713 | + <prognostic> |
1714 | + <mesh name="VelocityMesh"/> |
1715 | + <equation name="Boussinesq"/> |
1716 | + <spatial_discretisation> |
1717 | + <discontinuous_galerkin> |
1718 | + <viscosity_scheme> |
1719 | + <compact_discontinuous_galerkin/> |
1720 | + </viscosity_scheme> |
1721 | + <advection_scheme> |
1722 | + <upwind/> |
1723 | + <project_velocity_to_continuous> |
1724 | + <mesh name="BaseMesh"/> |
1725 | + </project_velocity_to_continuous> |
1726 | + <integrate_advection_by_parts> |
1727 | + <twice/> |
1728 | + </integrate_advection_by_parts> |
1729 | + </advection_scheme> |
1730 | + </discontinuous_galerkin> |
1731 | + <conservative_advection> |
1732 | + <real_value rank="0">0.0</real_value> |
1733 | + </conservative_advection> |
1734 | + </spatial_discretisation> |
1735 | + <temporal_discretisation> |
1736 | + <theta> |
1737 | + <real_value rank="0">1.0</real_value> |
1738 | + </theta> |
1739 | + <relaxation> |
1740 | + <real_value rank="0">1.0</real_value> |
1741 | + </relaxation> |
1742 | + <discontinuous_galerkin> |
1743 | + <maximum_courant_number_per_subcycle> |
1744 | + <real_value rank="0">0.1</real_value> |
1745 | + </maximum_courant_number_per_subcycle> |
1746 | + </discontinuous_galerkin> |
1747 | + </temporal_discretisation> |
1748 | + <solver> |
1749 | + <iterative_method name="gmres"> |
1750 | + <restart> |
1751 | + <integer_value rank="0">30</integer_value> |
1752 | + </restart> |
1753 | + </iterative_method> |
1754 | + <preconditioner name="sor"/> |
1755 | + <relative_error> |
1756 | + <real_value rank="0">1.0e-8</real_value> |
1757 | + </relative_error> |
1758 | + <max_iterations> |
1759 | + <integer_value rank="0">10000</integer_value> |
1760 | + </max_iterations> |
1761 | + <never_ignore_solver_failures/> |
1762 | + <diagnostics> |
1763 | + <monitors/> |
1764 | + </diagnostics> |
1765 | + </solver> |
1766 | + <initial_condition name="WholeMesh"> |
1767 | + <constant> |
1768 | + <real_value shape="3" dim1="dim" rank="1">0.0 0.0 0.0</real_value> |
1769 | + </constant> |
1770 | + </initial_condition> |
1771 | + <boundary_conditions name="FreeSurface"> |
1772 | + <surface_ids> |
1773 | + <integer_value shape="1" rank="1">1</integer_value> |
1774 | + </surface_ids> |
1775 | + <type name="free_surface"/> |
1776 | + </boundary_conditions> |
1777 | + <boundary_conditions name="NoNormalFlow"> |
1778 | + <surface_ids> |
1779 | + <integer_value shape="5" rank="1">2 1009 1010 1011 1013</integer_value> |
1780 | + </surface_ids> |
1781 | + <type name="no_normal_flow"/> |
1782 | + </boundary_conditions> |
1783 | + <tensor_field name="Viscosity" rank="2"> |
1784 | + <prescribed> |
1785 | + <value name="WholeMesh"> |
1786 | + <anisotropic_symmetric> |
1787 | + <python> |
1788 | + <string_value lines="20" type="code" language="python">def val(X, t): |
1789 | + a = X[0] |
1790 | + b = X[1] |
1791 | + c = X[2] |
1792 | + from math import sqrt, sin, cos, atan2, acos |
1793 | + r=sqrt(a**2+b**2+c**2) |
1794 | + #A1=2000.0 |
1795 | + A1=20000.0 |
1796 | + A2=A1 |
1797 | + A3=1.0 |
1798 | + phi=atan2(b,a) |
1799 | + theta=acos(c/r) |
1800 | + T11=A1*sin(phi)**2+A2*cos(theta)**2*cos(phi)**2+A3*sin(theta)**2*cos(phi)**2 |
1801 | + T12=-A1*sin(phi)*cos(phi)+A2*cos(theta)**2*sin(phi)*cos(phi)+A3*sin(theta)**2*sin(phi)*cos(phi) |
1802 | + T13=-A2*sin(theta)*cos(theta)*cos(phi)+A3*sin(theta)*cos(theta)*cos(phi) |
1803 | + T21=T12 |
1804 | + T22=A1*cos(phi)**2+A2*cos(theta)**2*sin(phi)**2+A3*sin(theta)**2*sin(phi)**2 |
1805 | + T23=-A2*sin(theta)*cos(theta)*sin(phi)+A3*sin(theta)*cos(theta)*sin(phi) |
1806 | + T31=T13 |
1807 | + T32=T23 |
1808 | + T33=A2*sin(theta)**2+A3*cos(theta)**2 |
1809 | + |
1810 | + return [[T11, T12, T13], |
1811 | + [T21, T22, T23], |
1812 | + [T31, T32, T33]]</string_value> |
1813 | + </python> |
1814 | + </anisotropic_symmetric> |
1815 | + </value> |
1816 | + <output/> |
1817 | + </prescribed> |
1818 | + </tensor_field> |
1819 | + <output/> |
1820 | + <stat> |
1821 | + <include_in_stat/> |
1822 | + <previous_time_step> |
1823 | + <exclude_from_stat/> |
1824 | + </previous_time_step> |
1825 | + <nonlinear_field> |
1826 | + <exclude_from_stat/> |
1827 | + </nonlinear_field> |
1828 | + </stat> |
1829 | + <convergence> |
1830 | + <include_in_convergence/> |
1831 | + </convergence> |
1832 | + <detectors> |
1833 | + <include_in_detectors/> |
1834 | + </detectors> |
1835 | + <steady_state> |
1836 | + <include_in_steady_state/> |
1837 | + </steady_state> |
1838 | + <consistent_interpolation/> |
1839 | + </prognostic> |
1840 | + </vector_field> |
1841 | + <scalar_field name="FreeSurface" rank="0"> |
1842 | + <diagnostic> |
1843 | + <algorithm name="Internal" material_phase_support="multiple"/> |
1844 | + <mesh name="PressureMesh"/> |
1845 | + <output/> |
1846 | + <stat/> |
1847 | + <convergence> |
1848 | + <include_in_convergence/> |
1849 | + </convergence> |
1850 | + <detectors> |
1851 | + <include_in_detectors/> |
1852 | + </detectors> |
1853 | + <steady_state> |
1854 | + <include_in_steady_state/> |
1855 | + </steady_state> |
1856 | + </diagnostic> |
1857 | + </scalar_field> |
1858 | + <scalar_field name="DG_CourantNumber" rank="0"> |
1859 | + <diagnostic> |
1860 | + <algorithm name="Internal" material_phase_support="multiple"/> |
1861 | + <mesh name="VelocityMesh"/> |
1862 | + <output/> |
1863 | + <stat/> |
1864 | + <convergence> |
1865 | + <include_in_convergence/> |
1866 | + </convergence> |
1867 | + <detectors> |
1868 | + <include_in_detectors/> |
1869 | + </detectors> |
1870 | + <steady_state> |
1871 | + <include_in_steady_state/> |
1872 | + </steady_state> |
1873 | + </diagnostic> |
1874 | + </scalar_field> |
1875 | + <scalar_field name="Partition" rank="0"> |
1876 | + <diagnostic> |
1877 | + <algorithm name="element_ownership" material_phase_support="single"/> |
1878 | + <mesh name="P0Mesh"/> |
1879 | + <output/> |
1880 | + <stat/> |
1881 | + <convergence> |
1882 | + <include_in_convergence/> |
1883 | + </convergence> |
1884 | + <detectors> |
1885 | + <include_in_detectors/> |
1886 | + </detectors> |
1887 | + <steady_state> |
1888 | + <include_in_steady_state/> |
1889 | + </steady_state> |
1890 | + </diagnostic> |
1891 | + </scalar_field> |
1892 | + </material_phase> |
1893 | +</fluidity_options> |
1894 | |
1895 | === added directory 'tests/spherical_patch/src' |
1896 | === added file 'tests/spherical_patch/src/spherical_segment.geo' |
1897 | --- tests/spherical_patch/src/spherical_segment.geo 1970-01-01 00:00:00 +0000 |
1898 | +++ tests/spherical_patch/src/spherical_segment.geo 2012-05-16 11:01:22 +0000 |
1899 | @@ -0,0 +1,113 @@ |
1900 | +Printf("Reading function and point definitions..."); |
1901 | + |
1902 | +Function DrawParallel |
1903 | + alpha = parallelSectionStartingY/parallelSectionStartingX; |
1904 | + parallelRadius = Sqrt(parallelSectionStartingX^2 + parallelSectionStartingY^2); |
1905 | + deltaAlpha = (parallelSectionEndingY/parallelSectionEndingX - parallelSectionStartingY/parallelSectionStartingX)/pointsOnParallel; |
1906 | + For t In {1:pointsOnParallel} |
1907 | + alpha += deltaAlpha; |
1908 | + newParallelPointX = (parallelSectionStartingX/Fabs(parallelSectionStartingX))*parallelRadius/(Sqrt(alpha^2 + 1)); |
1909 | + newParallelPointY = newParallelPointX*alpha; |
1910 | + newPointOnParallel = newp; Point(newPointOnParallel) = {newParallelPointX, newParallelPointY, 0}; |
1911 | + EndFor |
1912 | + BSpline(newParallelID) = {firstPointOnParallel, newPointOnParallel-(pointsOnParallel-1):newPointOnParallel, lastPointOnParallel}; |
1913 | +Return |
1914 | + |
1915 | +//The coordinates of various points on the Surface of the Earth are converted to radians |
1916 | +//and stereographic coordinatesbelow. |
1917 | + |
1918 | +//Point A is located 2 deg 00'W, 48 deg 45'N. |
1919 | +Point_A_longitude_rad = -(2 + (00/60))*(Pi/180); |
1920 | +Point_A_latitude_rad = (48 + (45/60))*(Pi/180); |
1921 | +Point_A_stereographic_x = -Cos(Fabs(Point_A_longitude_rad))*Cos(Fabs(Point_A_latitude_rad))/( 1 + Sin(Fabs(Point_A_latitude_rad)) ); |
1922 | +Point_A_stereographic_y = Cos(Fabs(Point_A_latitude_rad))*Sin(Fabs(Point_A_longitude_rad))/( 1 + Sin(Fabs(Point_A_latitude_rad)) ); |
1923 | + |
1924 | +//Point B is located at 9 deg 46'W, 48 deg 45'N |
1925 | +Point_B_longitude_rad = -(9 + (46/60))*(Pi/180); |
1926 | +Point_B_latitude_rad = (48 + (45/60))*(Pi/180); |
1927 | +Point_B_stereographic_x = -Cos(Fabs(Point_B_longitude_rad))*Cos(Fabs(Point_B_latitude_rad))/( 1 + Sin(Fabs(Point_B_latitude_rad)) ); |
1928 | +Point_B_stereographic_y = Cos(Fabs(Point_B_latitude_rad))*Sin(Fabs(Point_B_longitude_rad))/( 1 + Sin(Fabs(Point_B_latitude_rad)) ); |
1929 | + |
1930 | +//Point C is located at 9 deg 45'W, 56 deg 45'N |
1931 | +Point_C_longitude_rad = -(9 + (45/60))*(Pi/180); |
1932 | +Point_C_latitude_rad = (56 + (45/60))*(Pi/180); |
1933 | +Point_C_stereographic_x = -Cos(Fabs(Point_C_longitude_rad))*Cos(Fabs(Point_C_latitude_rad))/( 1 + Sin(Fabs(Point_C_latitude_rad)) ); |
1934 | +Point_C_stereographic_y = Cos(Fabs(Point_C_latitude_rad))*Sin(Fabs(Point_C_longitude_rad))/( 1 + Sin(Fabs(Point_C_latitude_rad)) ); |
1935 | + |
1936 | +//Point D is located at 2 deg 00'W, 56 deg 45'N |
1937 | +Point_D_longitude_rad = -(2 + (00/60))*(Pi/180); |
1938 | +Point_D_latitude_rad = (56 + (45/60))*(Pi/180); |
1939 | +Point_D_stereographic_x = -Cos(Fabs(Point_D_longitude_rad))*Cos(Fabs(Point_D_latitude_rad))/( 1 + Sin(Fabs(Point_D_latitude_rad)) ); |
1940 | +Point_D_stereographic_y = Cos(Fabs(Point_D_latitude_rad))*Sin(Fabs(Point_D_longitude_rad))/( 1 + Sin(Fabs(Point_D_latitude_rad)) ); |
1941 | + |
1942 | +//Create point, line, line-loop surface and field ID's |
1943 | +IP = newp; |
1944 | +IL = newl; |
1945 | +ILL = newll; |
1946 | +IS = news; |
1947 | +IFI = newf; |
1948 | + |
1949 | +//Draw a point on the centre of the earth and one on the North pole, then create a sphere using |
1950 | +// these points, note that remaining points will be drawin in the stereographic projection plane. |
1951 | +Point ( IP + 0 ) = {0, 0, 0 }; |
1952 | +Point ( IP + 1 ) = {0, 0,6.37101e+06}; |
1953 | +PolarSphere ( IS + 0 ) = {IP , IP+1}; |
1954 | + |
1955 | +Printf("Creating Meridian & Parallel segments to close the domain..."); |
1956 | +//So far, only points on the shorelies and curves approximating shorelines have been |
1957 | +//constructed. we now add points in the sea, to close the domain boundaries |
1958 | +Point(IP + 60000) = {Point_A_stereographic_x, Point_A_stereographic_y, 0}; |
1959 | +Point(IP + 60001) = {Point_B_stereographic_x, Point_B_stereographic_y, 0}; |
1960 | +Point(IP + 60002) = {Point_C_stereographic_x, Point_C_stereographic_y, 0}; |
1961 | +Point(IP + 60003) = {Point_D_stereographic_x, Point_D_stereographic_y, 0}; |
1962 | + |
1963 | +//Draw South-most parallel of the Domain. |
1964 | +pointsOnParallel = 200; //Assign parameters to variables and then call function DrawParallel, |
1965 | +parallelSectionStartingX = Point_A_stereographic_x; |
1966 | +parallelSectionStartingY = Point_A_stereographic_y; |
1967 | +firstPointOnParallel = IP + 60000; |
1968 | +parallelSectionEndingX = Point_B_stereographic_x; |
1969 | +parallelSectionEndingY = Point_B_stereographic_y; |
1970 | +lastPointOnParallel = IP + 60001; |
1971 | +newParallelID = IL + 1000; |
1972 | +Call DrawParallel; |
1973 | + |
1974 | +BSpline(IL + 1001) = {IP + 60001, IP + 60002}; |
1975 | +BSpline(IL + 1002) = {IP + 60000, IP + 60003}; |
1976 | + |
1977 | +////Draw North-most parallel of the Domain. |
1978 | +pointsOnParallel = 200; //Assign parameters to variables and then call function DrawParallel, |
1979 | +parallelSectionStartingX = Point_C_stereographic_x; |
1980 | +parallelSectionStartingY = Point_C_stereographic_y; |
1981 | +firstPointOnParallel = IP + 60002; |
1982 | +parallelSectionEndingX = Point_D_stereographic_x; |
1983 | +parallelSectionEndingY = Point_D_stereographic_y; |
1984 | +lastPointOnParallel = IP + 60003; |
1985 | +newParallelID = IL + 1003; |
1986 | +Call DrawParallel; |
1987 | + |
1988 | +Printf("Declaring surface on sphere to be meshed..."); |
1989 | +Line Loop(1006) = {1003, -1004, -1002, -1001}; |
1990 | +Plane Surface(1007) = {1006}; |
1991 | +Physical Surface(1008) = {1007}; |
1992 | + |
1993 | +Physical Line(1009) = {1001}; //Physical line corresponding to south-most parallel. |
1994 | +Physical Line(1010) = {1002}; //Physical line corresponding to meridian segment through point A. |
1995 | +Physical Line(1011) = {1003}; //Physical line corresponding to meridian segment through point B. |
1996 | +Physical Line(1013) = {1004}; //Physical line corresponding to north-most parallel. |
1997 | + |
1998 | +Printf("Assigning characteristic mesh sizes..."); |
1999 | +Field[ IFI + 1] = MathEval; |
2000 | +Field[ IFI + 1].F = "50000"; |
2001 | + |
2002 | +Background Field = IFI + 1; |
2003 | + |
2004 | +// Dont extent the elements sizes from the boundary inside the domain |
2005 | +Mesh.CharacteristicLengthExtendFromBoundary = 0; |
2006 | + |
2007 | +//Set some options for better png output |
2008 | +General.Color.Background = {255,255,255}; |
2009 | +General.Color.BackgroundGradient = {255,255,255}; |
2010 | +General.Color.Foreground = Black; |
2011 | +Mesh.Color.Lines = {0,0,0}; |
2012 | + |
2013 | |
2014 | === modified file 'tools/Checkmesh.F90' |
2015 | --- tools/Checkmesh.F90 2012-04-03 12:24:30 +0000 |
2016 | +++ tools/Checkmesh.F90 2012-05-16 11:01:22 +0000 |
2017 | @@ -15,7 +15,7 @@ |
2018 | use linked_lists |
2019 | use meshdiagnostics |
2020 | use metric_tools |
2021 | - use mesh_files |
2022 | + use read_triangle |
2023 | use supermesh_construction |
2024 | use tetrahedron_intersection_module |
2025 | use iso_c_binding |
2026 | @@ -38,7 +38,7 @@ |
2027 | rformat = real_format() |
2028 | |
2029 | print "(a)", "Reading in mesh mesh with base name " // trim(filename) |
2030 | - positions = read_mesh_files(filename, quad_degree = 4) |
2031 | + positions = read_triangle_files(filename, quad_degree = 4) |
2032 | if(isparallel()) call read_halos(filename, positions) |
2033 | print "(a)", "Read successful" |
2034 | |
2035 | |
2036 | === modified file 'tools/Fladapt.F90' |
2037 | --- tools/Fladapt.F90 2012-04-16 12:54:50 +0000 |
2038 | +++ tools/Fladapt.F90 2012-05-16 11:01:22 +0000 |
2039 | @@ -76,6 +76,7 @@ |
2040 | type(vector_field) :: new_mesh_field |
2041 | type(vector_field), pointer :: new_mesh_field_ptr, old_mesh_field |
2042 | type(tensor_field) :: metric, t_edge_lengths |
2043 | + character(len=FIELD_NAME_LEN) :: mesh_format |
2044 | |
2045 | ! now turn into proper fortran strings (is there an easier way to do this?) |
2046 | do i=1, input_basename_len |
2047 | @@ -144,11 +145,10 @@ |
2048 | else |
2049 | call adapt_mesh(old_mesh_field, metric, new_mesh_field) |
2050 | end if |
2051 | - old_mesh_field => null() |
2052 | - old_mesh => null() |
2053 | |
2054 | + call get_option(trim(old_mesh_field%mesh%option_path)//"/from_file/format/name", mesh_format) |
2055 | ! Write the output mesh |
2056 | - call write_mesh_files(output_basename, new_mesh_field) |
2057 | + call write_mesh_files(output_basename, mesh_format, new_mesh_field) |
2058 | |
2059 | ! Deallocate |
2060 | do i = 1, size(states) |
2061 | @@ -183,7 +183,7 @@ |
2062 | FLAbort("External mesh field " // trim(mesh_field_name) // " not found in the system state") |
2063 | end if |
2064 | |
2065 | - mesh_field => extract_vector_field(state, mesh_field_name) |
2066 | + mesh_field => extract_vector_field(state, mesh_field_name) |
2067 | |
2068 | end subroutine find_mesh_field_to_adapt |
2069 | |
2070 | |
2071 | === modified file 'tools/periodise.F90' |
2072 | --- tools/periodise.F90 2011-02-28 14:43:05 +0000 |
2073 | +++ tools/periodise.F90 2012-05-16 11:01:22 +0000 |
2074 | @@ -20,7 +20,7 @@ |
2075 | integer :: status |
2076 | type(state_type), dimension(:), pointer :: states |
2077 | integer :: ierr |
2078 | - character(len=FIELD_NAME_LEN) :: external_name, periodic_name |
2079 | + character(len=FIELD_NAME_LEN) :: external_name, periodic_name, mesh_format |
2080 | type(mesh_type), pointer :: periodic_mesh, external_mesh |
2081 | type(vector_field) :: periodic_positions, external_positions |
2082 | integer :: stat, i, nstates |
2083 | @@ -83,7 +83,8 @@ |
2084 | |
2085 | ! Dump out the periodic mesh to disk: |
2086 | new_external_filename = trim(external_filename) // '_periodic' |
2087 | - call write_mesh_files(new_external_filename, periodic_positions) |
2088 | + call get_option(trim(external_mesh%option_path)//"/from_file/format/name", mesh_format) |
2089 | + call write_mesh_files(new_external_filename, mesh_format, periodic_positions) |
2090 | |
2091 | ! OK! Now we need to do some setting of options. |
2092 | call manipulate_options(external_mesh, trim(external_mesh%option_path), periodic_mesh, trim(periodic_mesh%option_path), new_external_filename) |
2093 | |
2094 | === modified file 'tools/project_to_continuous.F90' |
2095 | --- tools/project_to_continuous.F90 2011-09-09 13:24:20 +0000 |
2096 | +++ tools/project_to_continuous.F90 2012-05-16 11:01:22 +0000 |
2097 | @@ -35,7 +35,7 @@ |
2098 | use state_module |
2099 | use elements |
2100 | use fields |
2101 | - use mesh_files |
2102 | + use read_triangle |
2103 | use vtk_interfaces |
2104 | use sparse_tools |
2105 | use fefields |
2106 | @@ -72,7 +72,7 @@ |
2107 | |
2108 | call vtk_read_state(vtuname, dg_state, quad_degree=6) |
2109 | |
2110 | - cg_coordinate= read_mesh_files(meshname, quad_degree=6) |
2111 | + cg_coordinate= read_triangle_files(meshname, quad_degree=6) |
2112 | cg_mesh=cg_coordinate%mesh |
2113 | |
2114 | call allocate(lumped_mass, cg_mesh, "LumpedMass") |
2115 | |
2116 | === modified file 'tools/test_laplacian.F90' |
2117 | --- tools/test_laplacian.F90 2010-07-02 08:26:03 +0000 |
2118 | +++ tools/test_laplacian.F90 2012-05-16 11:01:22 +0000 |
2119 | @@ -9,7 +9,7 @@ |
2120 | ! 1x1x1 cube with the boundary conditions 1 and 2 applied on the |
2121 | ! sides of the first coordinates direction (e.g. test_laplacian.poly |
2122 | ! - create a mesh with 'triangle -a0.01 -e test_laplacian.poly' ) |
2123 | - use mesh_files |
2124 | + use read_triangle |
2125 | use fields |
2126 | use FEtools |
2127 | use elements |
2128 | @@ -70,7 +70,7 @@ |
2129 | |
2130 | call read_command_line(filename, degree, quad_degree) |
2131 | |
2132 | - positions=read_mesh_files(filename, quad_degree=quad_degree) |
2133 | + positions=read_triangle_files(filename, quad_degree=quad_degree) |
2134 | |
2135 | call insert(state, positions, "Coordinate") |
2136 |
Code looks good at a quick eyeball.
Build queue is showing red:
unittests: empty_populate_ state'] empty_populate_ state']
Warnings: 1; tests = ['./test_
Failures: 1; tests = ['./test_
Medium tests:
Summary of test problems with failures or warnings:
fladapt.xml: F
Any ideas?