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..afc06072fe 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,25 +263,33 @@ 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 = DGIFY
+ 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
end if
+ if (model%continuity .lt. 0) lprojection = DGIFY
if (present(number_of_partitions)) then
nparts = number_of_partitions
@@ -483,32 +498,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 +604,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 +701,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)
+
+
+
+
+
+
+