From dceddcb7b53c1b33268f382ee07e2cf158ea7150 Mon Sep 17 00:00:00 2001 From: James Percival Date: Thu, 9 Mar 2017 20:17:43 +0000 Subject: [PATCH 1/2] Add option to project DG fields when outputting. --- assemble/Makefile.dependencies | 14 +- femtools/FEFields.F90 | 440 ++++++++++++++++++++--- femtools/Fields_Manipulation.F90 | 18 + femtools/Halos_Diagnostics.F90 | 28 +- femtools/Makefile.dependencies | 18 +- femtools/State_Fields.F90 | 99 ----- femtools/VTK_interfaces.F90 | 196 +++++----- femtools/Write_State.F90 | 17 +- main/Makefile.dependencies | 28 +- population_balance/Makefile.dependencies | 18 +- schemas/fluidity_options.rnc | 18 + schemas/fluidity_options.rng | 21 ++ 12 files changed, 655 insertions(+), 260 deletions(-) diff --git a/assemble/Makefile.dependencies b/assemble/Makefile.dependencies index 5809574f9e..a0bcd390a8 100644 --- a/assemble/Makefile.dependencies +++ b/assemble/Makefile.dependencies @@ -23,9 +23,10 @@ Adapt_State.o ../include/adapt_state_module.mod: Adapt_State.F90 \ ../include/detector_parallel.mod ../include/diagnostic_fields_wrapper.mod \ ../include/diagnostic_fields_wrapper_new.mod \ ../include/diagnostic_variables.mod \ - ../include/discrete_properties_module.mod ../include/edge_length_module.mod \ - ../include/elements.mod ../include/eventcounter.mod ../include/fdebug.h \ - ../include/fefields.mod ../include/field_options.mod ../include/fields.mod \ + ../include/discrete_properties_module.mod ../include/dqmom.mod \ + ../include/edge_length_module.mod ../include/elements.mod \ + ../include/eventcounter.mod ../include/fdebug.h ../include/fefields.mod \ + ../include/field_options.mod ../include/fields.mod \ ../include/fields_halos.mod ../include/fldebug.mod ../include/futils.mod \ ../include/global_parameters.mod ../include/hadapt_extrude.mod \ ../include/hadapt_metric_based_extrude.mod ../include/halos.mod \ @@ -214,9 +215,10 @@ Diagnostic_Fields_Matrices.o ../include/diagnostic_fields_matrices.mod: \ Diagnostic_fields_wrapper.o ../include/diagnostic_fields_wrapper.mod: \ Diagnostic_fields_wrapper.F90 ../include/diagnostic_fields.mod \ - ../include/diagnostic_fields_matrices.mod ../include/equation_of_state.mod \ - ../include/fdebug.h ../include/fetools.mod ../include/field_derivatives.mod \ - ../include/field_options.mod ../include/fields.mod ../include/fldebug.mod \ + ../include/diagnostic_fields_matrices.mod ../include/dqmom.mod \ + ../include/equation_of_state.mod ../include/fdebug.h ../include/fetools.mod \ + ../include/field_derivatives.mod ../include/field_options.mod \ + ../include/fields.mod ../include/fldebug.mod \ ../include/free_surface_module.mod ../include/futils.mod \ ../include/geostrophic_pressure.mod ../include/global_parameters.mod \ ../include/momentum_diagnostic_fields.mod \ diff --git a/femtools/FEFields.F90 b/femtools/FEFields.F90 index e6a7fe1e0e..b0e0a7aa6d 100644 --- a/femtools/FEFields.F90 +++ b/femtools/FEFields.F90 @@ -3,18 +3,23 @@ module fefields !!< Module containing general tools for discretising Finite Element problems. use fldebug + use global_parameters, only: FIELD_NAME_LEN use data_structures use element_numbering use elements, only: element_type use parallel_tools use sparse_tools use transform_elements, only: transform_to_physical, element_volume - use fetools, only: shape_shape, shape_rhs, shape_vector_rhs + use fetools, only: shape_shape, shape_rhs, shape_vector_rhs, shape_tensor_rhs use fields use state_module use field_options, only: get_coordinate_field use halos use sparse_matrices_fields + use sparsity_patterns + use sparsity_patterns_meshes + use eventcounter + use solvers, only: petsc_solve implicit none interface add_source_to_rhs @@ -22,14 +27,22 @@ module fefields end interface add_source_to_rhs interface project_field - module procedure project_scalar_field, project_vector_field + module procedure project_scalar_field, project_vector_field, project_tensor_field end interface + + interface get_lumped_mass + module procedure get_lumped_mass_single_state, get_lumped_mass_multiple_states + end interface get_lumped_mass + + interface get_mass_matrix + module procedure get_mass_matrix_single_state, get_mass_matrix_multiple_states + end interface get_mass_matrix private public :: compute_lumped_mass, compute_mass, compute_projection_matrix, add_source_to_rhs, & compute_lumped_mass_on_submesh, compute_cv_mass, project_field - public :: create_subdomain_mesh - + public :: create_subdomain_mesh, get_lumped_mass, get_mass_matrix + contains subroutine compute_cv_mass(positions, cv_mass, porosity) @@ -441,17 +454,24 @@ end subroutine projection_matrix_element end function compute_projection_matrix - subroutine project_scalar_field(from_field, to_field, X) + subroutine project_scalar_field(from_field, to_field, X, option_path, state) !!< Project from_field onto to_field. If to_field is discontinuous then !!< this will be calculated using the full mass matrix locally. - !!< Otherwise, the mass will be lumped. + !!< Otherwise, the mass will be lumped, unless an option path for a + !!< solver node is provided, in which case a linear solve occurs + type(scalar_field), intent(in) :: from_field type(scalar_field), intent(inout) :: to_field type(vector_field), intent(in) :: X + character(len=*), optional :: option_path + type(state_type), optional, intent(inout) :: state type(scalar_field) :: masslump + type(scalar_field), pointer :: pmass integer :: ele + type(csr_sparsity) :: sparsity type(csr_matrix) :: P + type(csr_matrix), pointer :: mass type(mesh_type) :: dg_mesh, cg_mesh @@ -472,26 +492,62 @@ subroutine project_scalar_field(from_field, to_field, X) ! DG to CG case cg_mesh=to_field%mesh - - call allocate(masslump, cg_mesh, "LumpedMass") - - call compute_lumped_mass(X, masslump) - ! Invert lumped mass. - masslump%val=1./masslump%val - dg_mesh=from_field%mesh - + P=compute_projection_matrix(cg_mesh, dg_mesh , X) + + if (present(option_path)) then + !! do a full Galerkin projection + + call zero(to_field) + if (present(state)) then + mass => get_mass_matrix_single_state(state,cg_mesh) + else + allocate(mass) + sparsity = make_sparsity(cg_mesh, cg_mesh, "MassSparsity") + call allocate(mass, sparsity, name="MassMatrix") + + call compute_mass(X, to_field%mesh, mass) + end if + + ! Perform projection. + call mult(to_field, P, from_field) + + call petsc_solve(to_field, mass, to_field, option_path=option_path) + + if (.not. present(state)) then + call deallocate(mass) + deallocate(mass) + call deallocate(sparsity) + end if + + else + + !! use lumped mass + + call allocate(masslump, cg_mesh, "LumpedMass") + + if (present(state)) then + pmass=>get_lumped_mass_single_state(state, cg_mesh) + call invert(pmass, masslump) + else + call compute_lumped_mass(X, masslump) + ! Invert lumped mass. + masslump%val=1./masslump%val + end if - call zero(to_field) - ! Perform projection. - call mult(to_field, P, from_field) - ! Apply inverted lumped mass to projected quantity. - call scale(to_field, masslump) + call zero(to_field) + ! Perform projection. + call mult(to_field, P, from_field) + ! Apply inverted lumped mass to projected quantity. + call scale(to_field, masslump) + + call deallocate(masslump) - call deallocate(masslump) - call deallocate(P) + end if + call deallocate(P) + end if end if @@ -524,17 +580,22 @@ end subroutine dg_projection_ele end subroutine project_scalar_field - subroutine project_vector_field(from_field, to_field, X) + subroutine project_vector_field(from_field, to_field, X, option_path, state) !!< Project from_field onto to_field. If to_field is discontinuous then !!< this will be calculated using the full mass matrix locally. !!< Otherwise, the mass will be lumped. type(vector_field), intent(in) :: from_field type(vector_field), intent(inout) :: to_field type(vector_field), intent(in) :: X + character(len=*), optional :: option_path + type(state_type), optional, intent(inout) :: state type(scalar_field) :: masslump, dg_scalar, cg_scalar + type(scalar_field), pointer :: pmass integer :: ele + type(csr_sparsity) :: sparsity type(csr_matrix) :: P + type(csr_matrix), pointer :: mass type(mesh_type) :: dg_mesh, cg_mesh integer :: j @@ -556,34 +617,70 @@ subroutine project_vector_field(from_field, to_field, X) ! CG case cg_mesh=to_field%mesh + dg_mesh=from_field%mesh + P=compute_projection_matrix(cg_mesh, dg_mesh , X) - call allocate(masslump, cg_mesh, "LumpedMass") - - call compute_lumped_mass(X, masslump) - ! Invert lumped mass. - masslump%val=1./masslump%val + if (present(option_path)) then + !! do a full Galerkin projection + + call zero(to_field) + if (present(state)) then + mass => get_mass_matrix_single_state(state,cg_mesh) + else + allocate(mass) + sparsity = make_sparsity(cg_mesh, cg_mesh, "MassSparsity") + call allocate(mass, sparsity, name="MassMatrix") + + call compute_mass(X, to_field%mesh, mass) + end if + + ! Perform projection. + do j=1,to_field%dim + cg_scalar=extract_scalar_field_from_vector_field(to_field, j) + dg_scalar=extract_scalar_field_from_vector_field(from_field, j) + call mult(cg_scalar, P, dg_scalar) + call petsc_solve(cg_scalar, mass, cg_scalar, option_path=option_path) + call set(to_field, j, cg_scalar) + end do + + if (.not. present(state)) then + call deallocate(mass) + deallocate(mass) + call deallocate(sparsity) + end if + + else + + !! use lumped mass - dg_mesh=from_field%mesh + call allocate(masslump, cg_mesh, "LumpedMass") - P=compute_projection_matrix(cg_mesh, dg_mesh, X) + if (present(state)) then + pmass=>get_lumped_mass_single_state(state, cg_mesh) + call invert(pmass, masslump) + else + call compute_lumped_mass(X, masslump) + ! Invert lumped mass. + masslump%val=1./masslump%val + end if - call zero(to_field) - ! Perform projection. - do j=1,to_field%dim - cg_scalar=extract_scalar_field_from_vector_field(to_field, j) - dg_scalar=extract_scalar_field_from_vector_field(from_field, j) - call mult(cg_scalar, P, dg_scalar) - call set(to_field, j, cg_scalar) - end do - - ! Apply inverted lumped mass to projected quantity. - call scale(to_field, masslump) - - call deallocate(masslump) - call deallocate(P) + call zero(to_field) + ! Perform projection. + do j=1,to_field%dim + cg_scalar=extract_scalar_field_from_vector_field(to_field, j) + dg_scalar=extract_scalar_field_from_vector_field(from_field, j) + call mult(cg_scalar, P, dg_scalar) + call set(to_field, j, cg_scalar) + end do + + ! Apply inverted lumped mass to projected quantity. + call scale(to_field, masslump) + + call deallocate(masslump) + end if + call deallocate(P) end if - end if contains @@ -639,6 +736,170 @@ subroutine cg_projection_ele(ele, from_field, to_field, masslump, X) end subroutine cg_projection_ele end subroutine project_vector_field + + subroutine project_tensor_field(from_field, to_field, X, option_path, state) + !!< Project from_field onto to_field. If to_field is discontinuous then + !!< this will be calculated using the full mass matrix locally. + !!< Otherwise, the mass will be lumped. + type(tensor_field), intent(in) :: from_field + type(tensor_field), intent(inout) :: to_field + type(vector_field), intent(in) :: X + character(len=*), optional :: option_path + type(state_type), optional, intent(inout) :: state + + type(scalar_field) :: masslump, dg_scalar, cg_scalar + type(scalar_field), pointer :: pmass + integer :: ele + type(csr_sparsity) :: sparsity + type(csr_matrix) :: P + type(csr_matrix), pointer :: mass + type(mesh_type) :: dg_mesh, cg_mesh + integer :: j, k + + + if(from_field%mesh==to_field%mesh) then + + call set(to_field, from_field) + + else + + if (to_field%mesh%continuity<0) then + ! DG case + + do ele=1,element_count(to_field) + call dg_projection_ele(ele, from_field, to_field, X) + end do + + else + ! CG case + + cg_mesh=to_field%mesh + dg_mesh=from_field%mesh + P=compute_projection_matrix(cg_mesh, dg_mesh , X) + + if (present(option_path)) then + !! do a full Galerkin projection + + call zero(to_field) + if (present(state)) then + mass => get_mass_matrix_single_state(state,cg_mesh) + else + allocate(mass) + sparsity = make_sparsity(cg_mesh, cg_mesh, "MassSparsity") + call allocate(mass, sparsity, name="MassMatrix") + + call compute_mass(X, to_field%mesh, mass) + end if + + ! Perform projection. + do k=1,to_field%dim(2) + do j=1,to_field%dim(1) + cg_scalar=extract_scalar_field_from_tensor_field(to_field, j, k) + dg_scalar=extract_scalar_field_from_tensor_field(from_field, j, k) + call mult(cg_scalar, P, dg_scalar) + call petsc_solve(cg_scalar, mass, cg_scalar, option_path=option_path) + call set(to_field, j, k, cg_scalar) + end do + end do + + if (.not. present(state)) then + call deallocate(mass) + deallocate(mass) + call deallocate(sparsity) + end if + + else + + !! use lumped mass + call allocate(masslump, cg_mesh, "LumpedMass") + + if (present(state)) then + pmass=>get_lumped_mass_single_state(state, cg_mesh) + call invert(pmass, masslump) + else + call compute_lumped_mass(X, masslump) + ! Invert lumped mass. + masslump%val=1./masslump%val + end if + + call zero(to_field) + ! Perform projection. + do k=1,to_field%dim(2) + do j=1,to_field%dim(1) + cg_scalar=extract_scalar_field_from_tensor_field(to_field, j, k) + dg_scalar=extract_scalar_field_from_tensor_field(from_field, j, k) + call mult(cg_scalar, P, dg_scalar) + call set(to_field, j, k, cg_scalar) + end do + end do + + ! Apply inverted lumped mass to projected quantity. + call scale(to_field, masslump) + + call deallocate(masslump) + + end if + + call deallocate(P) + + end if + + end if + + contains + + subroutine dg_projection_ele(ele, from_field, to_field, X) + integer :: ele + type(tensor_field), intent(in) :: from_field + type(tensor_field), intent(inout) :: to_field + type(vector_field), intent(in) :: X + + real, dimension(ele_loc(to_field,ele), ele_loc(to_field,ele)) :: mass + real, dimension(ele_ngi(to_field,ele)) :: detwei + type(element_type), pointer :: to_shape + integer :: dim1, dim2 + + call transform_to_physical(X, ele, detwei) + + to_shape=>ele_shape(to_field, ele) + + mass=shape_shape(to_shape, to_shape, detwei) + + call invert(mass) + + do dim2 = 1, to_field%dim(2) + do dim1 = 1, to_field%dim(1) + call set(to_field, dim1, dim2, ele_nodes(to_field, ele), & + matmul(mass, shape_rhs(to_shape, & + ele_val_at_quad(from_field, dim1, dim2, ele)* detwei))) + end do + end do + + end subroutine dg_projection_ele + + subroutine cg_projection_ele(ele, from_field, to_field, masslump, X) + integer :: ele + type(tensor_field), intent(in) :: from_field + type(tensor_field), intent(inout) :: to_field + type(scalar_field), intent(inout) :: masslump + type(vector_field), intent(in) :: X + + real, dimension(ele_ngi(to_field,ele)) :: detwei + type(element_type), pointer :: to_shape + + to_shape=>ele_shape(to_field, ele) + + call transform_to_physical(X, ele, detwei) + + call addto(masslump, ele_nodes(to_field, ele), & + shape_rhs(to_shape, detwei)) + + call addto(to_field, ele_nodes(to_field, ele), & + shape_tensor_rhs(to_shape, ele_val_at_quad(from_field, ele), detwei)) + + end subroutine cg_projection_ele + + end subroutine project_tensor_field subroutine add_source_to_rhs_scalar(rhs, source, positions) !!< Add in a source field to the rhs of a FE equation, @@ -893,4 +1154,95 @@ subroutine generate_subdomain_halos(external_mesh,subdomain_mesh,node_list,inver end subroutine generate_subdomain_halos + function get_lumped_mass_single_state(state, mesh) result(lumped_mass) + !!< extracts the lumped mass from states or creates it if it doesn't find it + type(scalar_field), pointer :: lumped_mass + type(state_type), intent(inout) :: state + type(mesh_type), intent(inout) :: mesh + + type(state_type), dimension(1) :: states + + states = (/state/) + lumped_mass => get_lumped_mass(states, mesh) + state = states(1) + + end function get_lumped_mass_single_state + + function get_lumped_mass_multiple_states(states, mesh) result(lumped_mass) + !!< extracts the lumped mass from states or creates it if it doesn't find it + type(scalar_field), pointer :: lumped_mass + type(state_type), dimension(:), intent(inout) :: states + type(mesh_type), intent(inout) :: mesh + + integer :: stat + character(len=FIELD_NAME_LEN) :: name + type(scalar_field) :: temp_lumped_mass + type(vector_field), pointer :: positions + integer, save :: last_mesh_movement = -1 + + name = trim(mesh%name)//"LumpedMass" + + lumped_mass => extract_scalar_field(states, trim(name), stat) + + if((stat/=0).or.(eventcount(EVENT_MESH_MOVEMENT)/=last_mesh_movement)) then + + positions => extract_vector_field(states(1), "Coordinate") + call allocate(temp_lumped_mass, mesh, name=trim(name)) + call compute_lumped_mass(positions, temp_lumped_mass) + call insert(states, temp_lumped_mass, trim(name)) + call deallocate(temp_lumped_mass) + + lumped_mass => extract_scalar_field(states, trim(name)) + last_mesh_movement = eventcount(EVENT_MESH_MOVEMENT) + end if + + end function get_lumped_mass_multiple_states + + function get_mass_matrix_single_state(state, mesh) result(mass) + !!< extracts the mass from states or creates it if it doesn't find it + type(csr_matrix), pointer :: mass + type(state_type), intent(inout) :: state + type(mesh_type), intent(inout) :: mesh + + type(state_type), dimension(1) :: states + + states = (/state/) + mass => get_mass_matrix(states, mesh) + state = states(1) + + end function get_mass_matrix_single_state + + function get_mass_matrix_multiple_states(states, mesh) result(mass) + !!< extracts the mass from states or creates it if it doesn't find it + type(csr_matrix), pointer :: mass + type(state_type), dimension(:), intent(inout) :: states + type(mesh_type), intent(inout) :: mesh + + integer :: stat + character(len=FIELD_NAME_LEN) :: name + type(csr_matrix) :: temp_mass + type(csr_sparsity), pointer :: temp_mass_sparsity + type(vector_field), pointer :: positions + + integer, save :: last_mesh_movement = -1 + + name = trim(mesh%name)//"MassMatrix" + + mass => extract_csr_matrix(states, trim(name), stat) + + if((stat/=0).or.(eventcount(EVENT_MESH_MOVEMENT)/=last_mesh_movement)) then + positions => extract_vector_field(states(1), "Coordinate") + + temp_mass_sparsity => get_csr_sparsity_firstorder(states, mesh, mesh) + call allocate(temp_mass, temp_mass_sparsity, name=trim(name)) + call compute_mass(positions, mesh, temp_mass) + call insert(states, temp_mass, trim(name)) + call deallocate(temp_mass) + + mass => extract_csr_matrix(states, trim(name)) + last_mesh_movement = eventcount(EVENT_MESH_MOVEMENT) + end if + + end function get_mass_matrix_multiple_states + end module fefields diff --git a/femtools/Fields_Manipulation.F90 b/femtools/Fields_Manipulation.F90 index 5ee93b0926..c3a430e688 100644 --- a/femtools/Fields_Manipulation.F90 +++ b/femtools/Fields_Manipulation.F90 @@ -97,6 +97,7 @@ module fields_manipulation & set_vector_field_nodes, & & set_vector_field_nodes_dim, & & set_tensor_field_nodes, & + & set_tensor_field_nodes_dim, & & set_scalar_field_field, & & set_scalar_field_from_vector_field, & & set_vector_field_field, & @@ -1369,6 +1370,23 @@ subroutine set_tensor_field_nodes(field, node_numbers, val) field%val(:, :, node_numbers) = val end subroutine set_tensor_field_nodes + + subroutine set_tensor_field_nodes_dim(field, dim1, dim2, node_numbers, val) + !!< Set the tensor field at the specified nodes + !!< Does not work for constant fields + type(tensor_field), intent(inout) :: field + integer, intent(in) :: dim1, dim2 + integer, dimension(:), intent(in) :: node_numbers + real, intent(in), dimension(:) :: val + + assert(field%field_type==FIELD_TYPE_NORMAL) + assert(field%field_type==FIELD_TYPE_NORMAL) + assert(dim1>=1 .and. dim1<=field%dim(1)) + assert(dim2>=1 .and. dim2<=field%dim(2)) + + field%val(dim1, dim2, node_numbers) = val + + end subroutine set_tensor_field_nodes_dim subroutine set_tensor_field(field, val) !!< Sets tensor with constant value diff --git a/femtools/Halos_Diagnostics.F90 b/femtools/Halos_Diagnostics.F90 index c59a1391e4..35a2bda4fb 100644 --- a/femtools/Halos_Diagnostics.F90 +++ b/femtools/Halos_Diagnostics.F90 @@ -35,12 +35,38 @@ module halos_diagnostics use halos_numbering use fields_allocates use fields_manipulation - use vtk_interfaces +! use vtk_interfaces implicit none private public write_universal_numbering + interface + subroutine vtk_write_fields(filename, index, position, model, sfields,& + & vfields, tfields, write_region_ids, write_columns, number_of_partitions,& + projection, stat) + + use fields + + implicit none + + character(len=*), intent(in) :: filename ! Base filename with no + ! trailing _number.vtu + integer, intent(in), optional :: index ! Index number of dump for filename. + type(vector_field), intent(in) :: position + type(mesh_type), intent(in) :: model + type(scalar_field), dimension(:), intent(in), optional :: sfields + type(vector_field), dimension(:), intent(in), optional :: vfields + type(tensor_field), dimension(:), intent(in), optional :: tfields + logical, intent(in), optional :: write_region_ids + logical, intent(in), optional :: write_columns + !!< If present, only write for processes 1:number_of_partitions (assumes the other partitions are empty) + integer, optional, intent(in):: number_of_partitions + integer, optional, intent(in) :: projection + integer, intent(out), optional :: stat + end subroutine vtk_write_fields +end interface + contains subroutine write_universal_numbering(halo, mesh, position, name) diff --git a/femtools/Makefile.dependencies b/femtools/Makefile.dependencies index 8a2b7de90e..2af3e618a0 100644 --- a/femtools/Makefile.dependencies +++ b/femtools/Makefile.dependencies @@ -609,11 +609,10 @@ Halos_Derivation.o ../include/halos_derivation.mod: Halos_Derivation.F90 \ @true Halos_Diagnostics.o ../include/halos_diagnostics.mod: Halos_Diagnostics.F90 \ - ../include/fdebug.h ../include/fields_allocates.mod \ + ../include/fdebug.h ../include/fields.mod ../include/fields_allocates.mod \ ../include/fields_base.mod ../include/fields_data_types.mod \ ../include/fields_manipulation.mod ../include/halo_data_types.mod \ - ../include/halos_numbering.mod ../include/shape_functions.mod \ - ../include/vtk_interfaces.mod + ../include/halos_numbering.mod ../include/shape_functions.mod ../include/halos_numbering.mod: Halos_Numbering.o @true @@ -1260,12 +1259,13 @@ Unittest_tools.o ../include/unittest_tools.mod: Unittest_tools.F90 \ VTK_interfaces.o ../include/vtk_interfaces.mod: VTK_interfaces.F90 \ ../include/data_structures.mod ../include/elements.mod ../include/fdebug.h \ - ../include/fetools.mod ../include/fields.mod ../include/fldebug.mod \ - ../include/futils.mod ../include/global_numbering.mod \ - ../include/global_parameters.mod ../include/mpi_interfaces.mod \ - ../include/parallel_fields.mod ../include/parallel_tools.mod \ - ../include/quadrature.mod ../include/sparse_tools.mod \ - ../include/state_module.mod ../include/vtkfortran.mod + ../include/fefields.mod ../include/fetools.mod ../include/fields.mod \ + ../include/fldebug.mod ../include/futils.mod \ + ../include/global_numbering.mod ../include/global_parameters.mod \ + ../include/mpi_interfaces.mod ../include/parallel_fields.mod \ + ../include/parallel_tools.mod ../include/quadrature.mod \ + ../include/sparse_tools.mod ../include/state_module.mod \ + ../include/vtkfortran.mod ../include/vector_tools.mod: Vector_Tools.o @true diff --git a/femtools/State_Fields.F90 b/femtools/State_Fields.F90 index afd1c34f79..037f5ed303 100644 --- a/femtools/State_Fields.F90 +++ b/femtools/State_Fields.F90 @@ -41,14 +41,6 @@ module state_fields_module interface get_cv_mass module procedure get_cv_mass_single_state, get_cv_mass_multiple_states end interface get_cv_mass - - interface get_lumped_mass - module procedure get_lumped_mass_single_state, get_lumped_mass_multiple_states - end interface get_lumped_mass - - interface get_mass_matrix - module procedure get_mass_matrix_single_state, get_mass_matrix_multiple_states - end interface get_mass_matrix interface get_dg_inverse_mass module procedure get_dg_inverse_mass_single_state, get_dg_inverse_mass_multiple_states @@ -107,97 +99,6 @@ function get_cv_mass_multiple_states(states, mesh) result(cv_mass) end function get_cv_mass_multiple_states - function get_lumped_mass_single_state(state, mesh) result(lumped_mass) - !!< extracts the lumped mass from states or creates it if it doesn't find it - type(scalar_field), pointer :: lumped_mass - type(state_type), intent(inout) :: state - type(mesh_type), intent(inout) :: mesh - - type(state_type), dimension(1) :: states - - states = (/state/) - lumped_mass => get_lumped_mass(states, mesh) - state = states(1) - - end function get_lumped_mass_single_state - - function get_lumped_mass_multiple_states(states, mesh) result(lumped_mass) - !!< extracts the lumped mass from states or creates it if it doesn't find it - type(scalar_field), pointer :: lumped_mass - type(state_type), dimension(:), intent(inout) :: states - type(mesh_type), intent(inout) :: mesh - - integer :: stat - character(len=FIELD_NAME_LEN) :: name - type(scalar_field) :: temp_lumped_mass - type(vector_field), pointer :: positions - integer, save :: last_mesh_movement = -1 - - name = trim(mesh%name)//"LumpedMass" - - lumped_mass => extract_scalar_field(states, trim(name), stat) - - if((stat/=0).or.(eventcount(EVENT_MESH_MOVEMENT)/=last_mesh_movement)) then - - positions => extract_vector_field(states(1), "Coordinate") - call allocate(temp_lumped_mass, mesh, name=trim(name)) - call compute_lumped_mass(positions, temp_lumped_mass) - call insert(states, temp_lumped_mass, trim(name)) - call deallocate(temp_lumped_mass) - - lumped_mass => extract_scalar_field(states, trim(name)) - last_mesh_movement = eventcount(EVENT_MESH_MOVEMENT) - end if - - end function get_lumped_mass_multiple_states - - function get_mass_matrix_single_state(state, mesh) result(mass) - !!< extracts the mass from states or creates it if it doesn't find it - type(csr_matrix), pointer :: mass - type(state_type), intent(inout) :: state - type(mesh_type), intent(inout) :: mesh - - type(state_type), dimension(1) :: states - - states = (/state/) - mass => get_mass_matrix(states, mesh) - state = states(1) - - end function get_mass_matrix_single_state - - function get_mass_matrix_multiple_states(states, mesh) result(mass) - !!< extracts the mass from states or creates it if it doesn't find it - type(csr_matrix), pointer :: mass - type(state_type), dimension(:), intent(inout) :: states - type(mesh_type), intent(inout) :: mesh - - integer :: stat - character(len=FIELD_NAME_LEN) :: name - type(csr_matrix) :: temp_mass - type(csr_sparsity), pointer :: temp_mass_sparsity - type(vector_field), pointer :: positions - - integer, save :: last_mesh_movement = -1 - - name = trim(mesh%name)//"MassMatrix" - - mass => extract_csr_matrix(states, trim(name), stat) - - if((stat/=0).or.(eventcount(EVENT_MESH_MOVEMENT)/=last_mesh_movement)) then - positions => extract_vector_field(states(1), "Coordinate") - - temp_mass_sparsity => get_csr_sparsity_firstorder(states, mesh, mesh) - call allocate(temp_mass, temp_mass_sparsity, name=trim(name)) - call compute_mass(positions, mesh, temp_mass) - call insert(states, temp_mass, trim(name)) - call deallocate(temp_mass) - - mass => extract_csr_matrix(states, trim(name)) - last_mesh_movement = eventcount(EVENT_MESH_MOVEMENT) - end if - - end function get_mass_matrix_multiple_states - function get_dg_inverse_mass_single_state(state, mesh) result(inverse_mass) !!< extracts the dg inverse mass from state or creates it if it doesn't find it type(csr_matrix), pointer :: inverse_mass diff --git a/femtools/VTK_interfaces.F90 b/femtools/VTK_interfaces.F90 index 17bdedce04..eb874f3c55 100644 --- a/femtools/VTK_interfaces.F90 +++ b/femtools/VTK_interfaces.F90 @@ -46,6 +46,7 @@ module vtk_interfaces use parallel_fields use fields use state_module + use fefields use vtkfortran implicit none @@ -55,6 +56,9 @@ module vtk_interfaces public :: vtk_write_state, vtk_write_fields, vtk_read_state, & vtk_write_surface_mesh, vtk_write_internal_face_mesh, & vtk_get_sizes, vtk_read_file + + integer, parameter, public :: DGIFY=0, GAUSSIAN_PROJECTION=1, MASS_LUMPED_PROJECTION=2 + interface subroutine vtk_read_file(& @@ -219,7 +223,8 @@ subroutine vtk_write_state(filename, index, model, state, write_region_ids, writ end subroutine vtk_write_state subroutine vtk_write_fields(filename, index, position, model, sfields,& - & vfields, tfields, write_region_ids, write_columns, number_of_partitions, stat) + & vfields, tfields, write_region_ids, write_columns, number_of_partitions,& + projection, state, stat) !!< Write the state variables out to a vtu file. Two different elements !!< are supported along with fields corresponding to each of them. !!< @@ -240,6 +245,8 @@ subroutine vtk_write_fields(filename, index, position, model, sfields,& logical, intent(in), optional :: write_columns !!< If present, only write for processes 1:number_of_partitions (assumes the other partitions are empty) integer, optional, intent(in):: number_of_partitions + integer, optional, intent(in) :: projection + type(state_type), optional, intent(inout) :: state integer, intent(out), optional :: stat integer :: NNod, sz_enlist, i, dim, j, k, nparts @@ -256,21 +263,28 @@ subroutine vtk_write_fields(filename, index, position, model, sfields,& integer, allocatable, dimension(:)::ghost_levels real, allocatable, dimension(:,:) :: tempval integer :: lstat + integer :: lprojection + if (present(projection)) then + lprojection = projection + else + lprojection = 0 + end if + if (present(stat)) stat = 0 dgify_fields = .false. - if (present(sfields)) then + if (present(sfields) .and. lprojection == DGIFY) then do i=1,size(sfields) if ( (sfields(i)%mesh%continuity .lt. 0 .and. sfields(i)%mesh%shape%degree /= 0) ) dgify_fields = .true. end do end if - if (present(vfields)) then + if (present(vfields) .and. lprojection == DGIFY) then do i=1,size(vfields) if ( (vfields(i)%mesh%continuity .lt. 0 .and. vfields(i)%mesh%shape%degree /= 0) ) dgify_fields = .true. end do end if - if (present(tfields)) then + if (present(tfields) .and. lprojection == DGIFY ) then do i=1,size(tfields) if ( (tfields(i)%mesh%continuity .lt. 0 .and. tfields(i)%mesh%shape%degree /= 0) ) dgify_fields = .true. end do @@ -483,32 +497,40 @@ subroutine vtk_write_fields(filename, index, position, model, sfields,& if (sfields(i)%mesh%shape%degree /= 0) then - call remap_field(from_field=sfields(i), to_field=l_model, stat=lstat) - ! if this is being called from something other than the main output routines - ! then these tests can be disabled by passing in the optional stat argument - ! to vtk_write_fields - if(lstat==REMAP_ERR_DISCONTINUOUS_CONTINUOUS) then - if(present(stat)) then - stat = lstat - else - FLAbort("Just remapped from a discontinuous to a continuous field!") - end if - else if(lstat==REMAP_ERR_UNPERIODIC_PERIODIC) then - if(present(stat)) then - stat = lstat - else - FLAbort("Just remapped from an unperiodic to a periodic continuous field!") - end if - else if ((lstat/=0).and. & - (lstat/=REMAP_ERR_BUBBLE_LAGRANGE).and. & - (lstat/=REMAP_ERR_HIGHER_LOWER_CONTINUOUS)) then - if(present(stat)) then - stat = lstat - else - FLAbort("Unknown error when remapping field.") - end if + if (lprojection == DGIFY) then + call remap_field(from_field=sfields(i), to_field=l_model, stat=lstat) + ! if this is being called from something other than the main output routines + ! then these tests can be disabled by passing in the optional stat argument + ! to vtk_write_fields + if(lstat==REMAP_ERR_DISCONTINUOUS_CONTINUOUS) then + if(present(stat)) then + stat = lstat + else + FLAbort("Just remapped from a discontinuous to a continuous field!") + end if + else if(lstat==REMAP_ERR_UNPERIODIC_PERIODIC) then + if(present(stat)) then + stat = lstat + else + FLAbort("Just remapped from an unperiodic to a periodic continuous field!") + end if + else if ((lstat/=0).and. & + (lstat/=REMAP_ERR_BUBBLE_LAGRANGE).and. & + (lstat/=REMAP_ERR_HIGHER_LOWER_CONTINUOUS)) then + if(present(stat)) then + stat = lstat + else + FLAbort("Unknown error when remapping field.") + end if + end if + ! we've just allowed remapping from a higher order to a lower order continuous field + else if (lprojection == GAUSSIAN_PROJECTION) then + call project_field(from_field=sfields(i), to_field=l_model, X=position,& + state=state, option_path="/io/project/full_projection") + else if (lprojection == MASS_LUMPED_PROJECTION) then + call project_field(from_field=sfields(i), to_field=l_model, & + X=position, state=state) end if - ! we've just allowed remapping from a higher order to a lower order continuous field call vtkwritesn(l_model%val, trim(sfields(i)%name)) @@ -581,33 +603,40 @@ subroutine vtk_write_fields(filename, index, position, model, sfields,& if(mesh_dim(vfields(i))/=mesh_dim(v_model(vfields(i)%dim))) cycle if(vfields(i)%mesh%shape%degree /= 0) then - - call remap_field(from_field=vfields(i), to_field=v_model(vfields(i)%dim), stat=lstat) - ! if this is being called from something other than the main output routines - ! then these tests can be disabled by passing in the optional stat argument - ! to vtk_write_fields - if(lstat==REMAP_ERR_DISCONTINUOUS_CONTINUOUS) then - if(present(stat)) then - stat = lstat - else - FLAbort("Just remapped from a discontinuous to a continuous field!") - end if - else if(lstat==REMAP_ERR_UNPERIODIC_PERIODIC) then - if(present(stat)) then - stat = lstat - else - FLAbort("Just remapped from an unperiodic to a periodic continuous field!") - end if - else if ((lstat/=0).and. & + if (lprojection == DGIFY) then + call remap_field(from_field=vfields(i), to_field=v_model(vfields(i)%dim), stat=lstat) + ! if this is being called from something other than the main output routines + ! then these tests can be disabled by passing in the optional stat argument + ! to vtk_write_fields + if(lstat==REMAP_ERR_DISCONTINUOUS_CONTINUOUS) then + if(present(stat)) then + stat = lstat + else + FLAbort("Just remapped from a discontinuous to a continuous field!") + end if + else if(lstat==REMAP_ERR_UNPERIODIC_PERIODIC) then + if(present(stat)) then + stat = lstat + else + FLAbort("Just remapped from an unperiodic to a periodic continuous field!") + end if + else if ((lstat/=0).and. & (lstat/=REMAP_ERR_BUBBLE_LAGRANGE).and. & (lstat/=REMAP_ERR_HIGHER_LOWER_CONTINUOUS)) then - if(present(stat)) then - stat = lstat - else - FLAbort("Unknown error when remapping field.") - end if + if(present(stat)) then + stat = lstat + else + FLAbort("Unknown error when remapping field.") + end if + end if + ! we've just allowed remapping from a higher order to a lower order continuous field + else if (lprojection == GAUSSIAN_PROJECTION) then + call project_field(from_field=vfields(i), to_field=v_model(vfields(i)%dim), X=position,& + state=state, option_path="/io/project/full_projection") + else if (lprojection == MASS_LUMPED_PROJECTION) then + call project_field(from_field=vfields(i), to_field=v_model(vfields(i)%dim), & + X=position, state=state) end if - ! we've just allowed remapping from a higher order to a lower order continuous field do k=1, vfields(i)%dim v_field_buffer(:,k)=v_model(vfields(i)%dim)%val(k,:) @@ -671,33 +700,40 @@ subroutine vtk_write_fields(filename, index, position, model, sfields,& if(tfields(i)%mesh%shape%degree /= 0) then - call remap_field(from_field=tfields(i), to_field=t_model, stat=lstat) - ! if this is being called from something other than the main output routines - ! then these tests can be disabled by passing in the optional stat argument - ! to vtk_write_fields - if(lstat==REMAP_ERR_DISCONTINUOUS_CONTINUOUS) then - if(present(stat)) then - stat = lstat - else - FLAbort("Just remapped from a discontinuous to a continuous field!") - end if - else if(lstat==REMAP_ERR_UNPERIODIC_PERIODIC) then - if(present(stat)) then - stat = lstat - else - FLAbort("Just remapped from an unperiodic to a periodic continuous field!") - end if - else if ((lstat/=0).and. & - (lstat/=REMAP_ERR_BUBBLE_LAGRANGE).and. & - (lstat/=REMAP_ERR_HIGHER_LOWER_CONTINUOUS)) then - if(present(stat)) then - stat = lstat - else - FLAbort("Unknown error when remapping field.") - end if - end if - ! we've just allowed remapping from a higher order to a lower order continuous field - + if (lprojection == DGIFY) then + call remap_field(from_field=tfields(i), to_field=t_model, stat=lstat) + ! if this is being called from something other than the main output routines + ! then these tests can be disabled by passing in the optional stat argument + ! to vtk_write_fields + if(lstat==REMAP_ERR_DISCONTINUOUS_CONTINUOUS) then + if(present(stat)) then + stat = lstat + else + FLAbort("Just remapped from a discontinuous to a continuous field!") + end if + else if(lstat==REMAP_ERR_UNPERIODIC_PERIODIC) then + if(present(stat)) then + stat = lstat + else + FLAbort("Just remapped from an unperiodic to a periodic continuous field!") + end if + else if ((lstat/=0).and. & + (lstat/=REMAP_ERR_BUBBLE_LAGRANGE).and. & + (lstat/=REMAP_ERR_HIGHER_LOWER_CONTINUOUS)) then + if(present(stat)) then + stat = lstat + else + FLAbort("Unknown error when remapping field.") + end if + end if + ! we've just allowed remapping from a higher order to a lower order continuous field + else if (lprojection == GAUSSIAN_PROJECTION) then + call project_field(from_field=tfields(i), to_field=t_model, X=position,& + state=state, option_path="/io/project/full_projection") + else if (lprojection == MASS_LUMPED_PROJECTION) then + call project_field(from_field=tfields(i), to_field=t_model, & + X=position, state=state) + end if allocate(tensor_values(node_count(t_model), 3, 3)) tensor_values=0.0 do j=1,dim diff --git a/femtools/Write_State.F90 b/femtools/Write_State.F90 index bae52d1938..fb2f7e6220 100644 --- a/femtools/Write_State.F90 +++ b/femtools/Write_State.F90 @@ -195,7 +195,7 @@ subroutine write_state(dump_no, state) type(state_type), dimension(:), intent(inout) :: state character(len = OPTION_PATH_LEN) :: dump_filename, dump_format - integer :: max_dump_no, stat + integer :: max_dump_no, stat, projection ewrite(1, *) "In write_state" call profiler_tic("I/O") @@ -205,11 +205,18 @@ subroutine write_state(dump_no, state) dump_no = modulo(dump_no, max_dump_no) + projection = DGIFY + if (have_option("/io/project/mass_lumped")) then + projection = MASS_LUMPED_PROJECTION + else if (have_option("/io/project/full_projection")) then + projection = GAUSSIAN_PROJECTION + end if + call get_option("/io/dump_format", dump_format) select case(trim(dump_format)) case("vtk") ewrite(2, *) "Writing output " // int2str(dump_no) // " to vtu" - call vtk_write_state_new_options(dump_filename, dump_no, state) + call vtk_write_state_new_options(dump_filename, dump_no, state, projection=projection) case default FLAbort("Unrecognised dump file format.") end select @@ -222,7 +229,7 @@ subroutine write_state(dump_no, state) end subroutine write_state - subroutine vtk_write_state_new_options(filename, index, state, write_region_ids) + subroutine vtk_write_state_new_options(filename, index, state, write_region_ids, projection) !!< Write the state variables out to a vtu file according to options !!< set in the options tree. Only fields present in the option tree !!< will be written, except for those disabled in the same options tree. @@ -234,6 +241,7 @@ subroutine vtk_write_state_new_options(filename, index, state, write_region_ids) integer, intent(in), optional :: index !! Index number of dump for filename. type(state_type), dimension(:), intent(inout) :: state logical, intent(in), optional :: write_region_ids + integer, intent(in), optional :: projection type(vector_field), pointer :: model_coordinate type(mesh_type), pointer :: model_mesh @@ -358,7 +366,8 @@ subroutine vtk_write_state_new_options(filename, index, state, write_region_ids) sfields=lsfields, & vfields=lvfields, & tfields=ltfields, & - write_region_ids=write_region_ids) + write_region_ids=write_region_ids, & + projection=projection, state=state(1)) ewrite(1, *) "Exiting vtk_write_state_new_options" diff --git a/main/Makefile.dependencies b/main/Makefile.dependencies index 4172d0e0cd..23fc0651a2 100644 --- a/main/Makefile.dependencies +++ b/main/Makefile.dependencies @@ -15,20 +15,20 @@ Fluids.o ../include/fluids_module.mod: Fluids.F90 \ ../include/diagnostic_children.mod ../include/diagnostic_fields_new.mod \ ../include/diagnostic_fields_wrapper.mod \ ../include/diagnostic_variables.mod \ - ../include/discrete_properties_module.mod ../include/elements.mod \ - ../include/equation_of_state.mod ../include/eventcounter.mod \ - ../include/fdebug.h ../include/field_equations_cv.mod \ - ../include/field_priority_lists.mod ../include/fields.mod \ - ../include/fldebug.mod ../include/foam_flow_module.mod \ - ../include/free_surface_module.mod ../include/futils.mod \ - ../include/global_parameters.mod ../include/gls.mod ../include/goals.mod \ - ../include/halos.mod ../include/iceshelf_meltrate_surf_normal.mod \ - ../include/implicit_solids.mod ../include/k_epsilon.mod \ - ../include/memory_diagnostics.mod ../include/meshdiagnostics.mod \ - ../include/meshmovement.mod ../include/momentum_diagnostic_fields.mod \ - ../include/momentum_equation.mod ../include/multimaterial_module.mod \ - ../include/multiphase_module.mod ../include/parallel_tools.mod \ - ../include/populate_state_module.mod \ + ../include/discrete_properties_module.mod ../include/dqmom.mod \ + ../include/elements.mod ../include/equation_of_state.mod \ + ../include/eventcounter.mod ../include/fdebug.h \ + ../include/field_equations_cv.mod ../include/field_priority_lists.mod \ + ../include/fields.mod ../include/fldebug.mod \ + ../include/foam_flow_module.mod ../include/free_surface_module.mod \ + ../include/futils.mod ../include/global_parameters.mod ../include/gls.mod \ + ../include/goals.mod ../include/halos.mod \ + ../include/iceshelf_meltrate_surf_normal.mod ../include/implicit_solids.mod \ + ../include/k_epsilon.mod ../include/memory_diagnostics.mod \ + ../include/meshdiagnostics.mod ../include/meshmovement.mod \ + ../include/momentum_diagnostic_fields.mod ../include/momentum_equation.mod \ + ../include/multimaterial_module.mod ../include/multiphase_module.mod \ + ../include/parallel_tools.mod ../include/populate_state_module.mod \ ../include/populate_sub_state_module.mod ../include/qmesh_module.mod \ ../include/reduced_model_runtime.mod ../include/reference_counting.mod \ ../include/reserve_state_module.mod \ diff --git a/population_balance/Makefile.dependencies b/population_balance/Makefile.dependencies index a4f50cde3b..770c52b14c 100644 --- a/population_balance/Makefile.dependencies +++ b/population_balance/Makefile.dependencies @@ -2,9 +2,21 @@ ../include/dqmom.mod: DQMOM.o @true -DQMOM.o ../include/dqmom.mod: DQMOM.F90 ../include/fdebug.h \ - ../include/fields.mod ../include/global_parameters.mod \ - ../include/initialise_fields_module.mod ../include/spud.mod \ +DQMOM.o ../include/dqmom.mod: DQMOM.F90 ../include/elements.mod \ + ../include/fdebug.h ../include/fetools.mod ../include/fields.mod \ + ../include/fldebug.mod ../include/futils.mod \ + ../include/global_parameters.mod ../include/initialise_fields_module.mod \ + ../include/solvers.mod ../include/sparse_tools.mod \ ../include/state_fields_module.mod ../include/state_module.mod \ ../include/vector_tools.mod +../include/dqmom.mod: DQMOM_node_source.o + @true + +DQMOM_node_source.o ../include/dqmom.mod: DQMOM_node_source.F90 \ + ../include/fdebug.h ../include/fields.mod ../include/fldebug.mod \ + ../include/futils.mod ../include/global_parameters.mod \ + ../include/initialise_fields_module.mod ../include/solvers.mod \ + ../include/sparse_tools.mod ../include/state_fields_module.mod \ + ../include/state_module.mod ../include/vector_tools.mod + diff --git a/schemas/fluidity_options.rnc b/schemas/fluidity_options.rnc index 893f13ea61..790b333b6b 100644 --- a/schemas/fluidity_options.rnc +++ b/schemas/fluidity_options.rnc @@ -135,6 +135,24 @@ start = attribute name { xsd:string } } ), + ## If a continuous model output mesh is selected, and choosing to + ## output discontinuous fields, the default behaviour is to + ## output on a discontinuous mesh of the same degree. + ## + ## Turning this option on causes a projection to be used instead. + ## Select either a mass lumped method, or a full Galerkin projection + ## (The full projection requires a linear solve) + element project { + ( + element mass_lumped { + empty + }| + element full_projection { + element solver { + linear_solver_options_sym + } + }) + }?, ## Options for convergence analysis. element convergence { ## Whether to enable the creation of a convergence diff --git a/schemas/fluidity_options.rng b/schemas/fluidity_options.rng index 3b7487a698..27c7b3574f 100644 --- a/schemas/fluidity_options.rng +++ b/schemas/fluidity_options.rng @@ -164,6 +164,27 @@ interpolated for VTK output. + + + If a continuous model output mesh is selected, and choosing to +output discontinuous fields, the default behaviour is to +output on a discontinuous mesh of the same degree. + +Turning this option on causes a projection to be used instead. +Select either a mass lumped method, or a full Galerkin projection +(The full projection requires a linear solve) + + + + + + + + + + + + Options for convergence analysis. From cbd03429b86e3bc4b91bcb8cee5495f9a31c674f Mon Sep 17 00:00:00 2001 From: James Percival Date: Mon, 20 Mar 2017 15:19:53 +0000 Subject: [PATCH 2/2] Add a test and make output on a DG mesh safe. --- femtools/VTK_interfaces.F90 | 3 +- tests/io_projection/io_projection.xml | 36 +++++ tests/io_projection/io_projection_full.flml | 152 ++++++++++++++++++ tests/io_projection/io_projection_lumped.flml | 137 ++++++++++++++++ tests/io_projection/src/square.geo | 16 ++ 5 files changed, 343 insertions(+), 1 deletion(-) create mode 100644 tests/io_projection/io_projection.xml create mode 100644 tests/io_projection/io_projection_full.flml create mode 100644 tests/io_projection/io_projection_lumped.flml create mode 100644 tests/io_projection/src/square.geo diff --git a/femtools/VTK_interfaces.F90 b/femtools/VTK_interfaces.F90 index eb874f3c55..afc06072fe 100644 --- a/femtools/VTK_interfaces.F90 +++ b/femtools/VTK_interfaces.F90 @@ -268,7 +268,7 @@ subroutine vtk_write_fields(filename, index, position, model, sfields,& if (present(projection)) then lprojection = projection else - lprojection = 0 + lprojection = DGIFY end if if (present(stat)) stat = 0 @@ -289,6 +289,7 @@ subroutine vtk_write_fields(filename, index, position, model, sfields,& if ( (tfields(i)%mesh%continuity .lt. 0 .and. tfields(i)%mesh%shape%degree /= 0) ) dgify_fields = .true. end do end if + if (model%continuity .lt. 0) lprojection = DGIFY if (present(number_of_partitions)) then nparts = number_of_partitions diff --git a/tests/io_projection/io_projection.xml b/tests/io_projection/io_projection.xml new file mode 100644 index 0000000000..702c731920 --- /dev/null +++ b/tests/io_projection/io_projection.xml @@ -0,0 +1,36 @@ + + + io_projection + + + fluidity io_projection_lumped.flml +fluidity io_projection_full.flml + + + from fluidity_tools import stat_parser +s=stat_parser('io_projection_lumped.stat') +sdata=s['Fluid']['Scalar'] +vdata=s['Fluid']['Velocity%magnitude'] +tdata=s['Fluid']['Tensor%magnitude'] +lumped_stats=[[sdata['min'],vdata['min'],tdata['min']], + [sdata['max'],vdata['max'],tdata['max']], + sdata['integral']] + from fluidity_tools import stat_parser +s=stat_parser('io_projection_full.stat') +sdata=s['Fluid']['Scalar'] +vdata=s['Fluid']['Velocity%magnitude'] +tdata=s['Fluid']['Tensor%magnitude'] +full_stats=[[sdata['min'],vdata['min'],tdata['min']], + [sdata['max'],vdata['max'],tdata['max']], + sdata['integral']] + + + assert(all([abs(x[0]-x[1])<1.0e-8 for x in zip(lumped_stats[0],[0, 1, 0])]) and + all([abs(x[0]-x[1])<1.0e-8 for x in zip(full_stats[0],[0, 1, 0])])) + from numpy import sqrt +assert(all([abs(x[0]-x[1])<1.0e-8 for x in zip(lumped_stats[1],[1, 1, sqrt(2)])]) and + all([abs(x[0]-x[1])<1.0e-8 for x in zip(full_stats[1],[1, 1, sqrt(2)])])) + from numpy import sqrt +assert(abs(lumped_stats[2]-0.5)<1.0e-8 and abs(full_stats[2]-0.5)<1.0e-8) + + diff --git a/tests/io_projection/io_projection_full.flml b/tests/io_projection/io_projection_full.flml new file mode 100644 index 0000000000..20fa186cd9 --- /dev/null +++ b/tests/io_projection/io_projection_full.flml @@ -0,0 +1,152 @@ + + + + io_projection_full + + + fluids + + + + 2 + + + + + + + + + + + + + + discontinuous + + + + + + + + + 5 + + + + + + vtk + + + + 0 + + + + + + + + + + + 1.0e-7 + + + 1000 + + + + + + + + + + + + + 0.0 + + + 1.0 + + + 1.0 + + + + + + + + + def val(X,t): + if 'n' in persistent: + n=persistent['n'] + else: + n=0 + + persistent['n'] = n+1 + return (n%2,(n+1)%2) + + + + + + + + + + + + + + + + + def val(X,t): + if 'n' in persistent: + n=persistent['n'] + else: + n=0 + + persistent['n'] = n+1 + return n%2 + + + + + + + + + + + + + + + def val(X,t): + if 'n' in persistent: + n=persistent['n'] + else: + n=0 + + persistent['n'] = n+1 + return [n%2,n%2] + + + + + + + + + + + + + diff --git a/tests/io_projection/io_projection_lumped.flml b/tests/io_projection/io_projection_lumped.flml new file mode 100644 index 0000000000..82e5429256 --- /dev/null +++ b/tests/io_projection/io_projection_lumped.flml @@ -0,0 +1,137 @@ + + + + io_projection_lumped + + + fluids + + + + 2 + + + + + + + + + + + + + + discontinuous + + + + + + + + + 5 + + + + + + vtk + + + + 0 + + + + + + + + + + + + 0.0 + + + 1.0 + + + 1.0 + + + + + + + + + def val(X,t): + if 'n' in persistent: + n=persistent['n'] + else: + n=0 + + persistent['n'] = n+1 + return (n%2,(n+1)%2) + + + + + + + + + + + + + + + + + def val(X,t): + if 'n' in persistent: + n=persistent['n'] + else: + n=0 + + persistent['n'] = n+1 + return n%2 + + + + + + + + + + + + + + + def val(X,t): + if 'n' in persistent: + n=persistent['n'] + else: + n=0 + + persistent['n'] = n+1 + return [n%2,n%2] + + + + + + + + + + + + + diff --git a/tests/io_projection/src/square.geo b/tests/io_projection/src/square.geo new file mode 100644 index 0000000000..4282042c19 --- /dev/null +++ b/tests/io_projection/src/square.geo @@ -0,0 +1,16 @@ +Point(1) = {0,0,0}; +Point(2) = {1,0,0}; +Point(3) = {1,1,0}; +Point(4) = {0,1,0}; + +Line(1) = {1,2}; +Line(2) = {2,3}; +Line(3) = {3,4}; +Line(4) = {4,1}; + +Line Loop(1) = {1,2,3,4}; + +Plane Surface(1) = 1; + +Transfinite Line {1,2,3,4} = 8; +Transfinite Surface{1}; \ No newline at end of file