From 063b5f145936d2260c3c71cdd229e6d1da129a5b Mon Sep 17 00:00:00 2001 From: Chris Bradley Date: Wed, 29 Jul 2009 03:30:50 +0000 Subject: [PATCH] Fixes for incorrect residual and RHS for nonlinear problems with linear matrices. Fixes for new libraries with Fluid makefiles. Format tidyups. git-svn-id: https://opencmiss.svn.sourceforge.net/svnroot/opencmiss/cm/trunk@595 61e36070-c886-4827-b7c0-ed9e4dc59cc8 --- ExampleMakefile | 2 +- src/Navier_Stokes_equations_routines.f90 | 1649 +++++++++++----------- src/Stokes_equations_routines.f90 | 458 +++--- src/solver_routines.f90 | 6 +- 4 files changed, 1055 insertions(+), 1060 deletions(-) diff --git a/ExampleMakefile b/ExampleMakefile index 46fdd2fe..e885df08 100644 --- a/ExampleMakefile +++ b/ExampleMakefile @@ -5,7 +5,7 @@ # $Id: ExampleMakefile 27 2007-07-24 16:52:51Z cpb $ # #---------------------------------------------------------------------------------------------------------------------------------- -# ExampleMakefile for compiling openCMISS examples +# ExampleMakefile for compiling OpenCMISS examples # # Original by Chris Bradley adapted from the CMISS Makefile by Karl Tomlinson # Changes: diff --git a/src/Navier_Stokes_equations_routines.f90 b/src/Navier_Stokes_equations_routines.f90 index 59842fe5..d97ff174 100644 --- a/src/Navier_Stokes_equations_routines.f90 +++ b/src/Navier_Stokes_equations_routines.f90 @@ -1,7 +1,7 @@ !> \file !> $Id: Stokes_fluid_routines.f90 372 2009-04-20 !> \author Sebastian Krittian -!> \brief This module handles all Stokes fluid routines. +!> \brief This module handles all Navier-Stokes fluid routines. !> !> \section LICENSE !> @@ -40,7 +40,7 @@ !> the terms of any one of the MPL, the GPL or the LGPL. !> -!>This module handles all Stokes fluid routines. +!>This module handles all Navier-Stokes fluid routines. MODULE NAVIER_STOKES_EQUATIONS_ROUTINES USE BASE_ROUTINES @@ -126,7 +126,7 @@ SUBROUTINE NAVIER_STOKES_EQUATIONS_SET_SOLUTION_METHOD_SET(EQUATIONS_SET,SOLUTIO CASE DEFAULT LOCAL_ERROR="Equations set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " is not valid for a Stokes flow equation type of a fluid mechanics equations set class." + & " is not valid for a Navier-Stokes flow equation type of a fluid mechanics equations set class." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) END SELECT ELSE @@ -160,26 +160,26 @@ SUBROUTINE NAVIER_STOKES_EQUATIONS_SET_SUBTYPE_SET(EQUATIONS_SET,EQUATIONS_SET_S IF(ASSOCIATED(EQUATIONS_SET)) THEN SELECT CASE(EQUATIONS_SET_SUBTYPE) CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE) - EQUATIONS_SET%CLASS=EQUATIONS_SET_FLUID_MECHANICS_CLASS - EQUATIONS_SET%TYPE=EQUATIONS_SET_NAVIER_STOKES_TYPE - EQUATIONS_SET%SUBTYPE=EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE + EQUATIONS_SET%CLASS=EQUATIONS_SET_FLUID_MECHANICS_CLASS + EQUATIONS_SET%TYPE=EQUATIONS_SET_NAVIER_STOKES_TYPE + EQUATIONS_SET%SUBTYPE=EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE CASE(EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE) - EQUATIONS_SET%CLASS=EQUATIONS_SET_FLUID_MECHANICS_CLASS - EQUATIONS_SET%TYPE=EQUATIONS_SET_NAVIER_STOKES_TYPE - EQUATIONS_SET%SUBTYPE=EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE + EQUATIONS_SET%CLASS=EQUATIONS_SET_FLUID_MECHANICS_CLASS + EQUATIONS_SET%TYPE=EQUATIONS_SET_NAVIER_STOKES_TYPE + EQUATIONS_SET%SUBTYPE=EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE CASE(EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) - EQUATIONS_SET%CLASS=EQUATIONS_SET_FLUID_MECHANICS_CLASS - EQUATIONS_SET%TYPE=EQUATIONS_SET_NAVIER_STOKES_TYPE - EQUATIONS_SET%SUBTYPE=EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE + EQUATIONS_SET%CLASS=EQUATIONS_SET_FLUID_MECHANICS_CLASS + EQUATIONS_SET%TYPE=EQUATIONS_SET_NAVIER_STOKES_TYPE + EQUATIONS_SET%SUBTYPE=EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE CASE(EQUATIONS_SET_OPTIMISED_NAVIER_STOKES_SUBTYPE) - CALL FLAG_ERROR("Not implemented yet.",ERR,ERROR,*999) + CALL FLAG_ERROR("Not implemented yet.",ERR,ERROR,*999) CASE DEFAULT - LOCAL_ERROR="Equations set subtype "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SUBTYPE,"*",ERR,ERROR))// & - & " is not valid for a Stokes fluid type of a fluid mechanics equations set class." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + LOCAL_ERROR="Equations set subtype "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SUBTYPE,"*",ERR,ERROR))// & + & " is not valid for a Stokes fluid type of a fluid mechanics equations set class." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) END SELECT ELSE - CALL FLAG_ERROR("Equations set is not associated.",ERR,ERROR,*999) + CALL FLAG_ERROR("Equations set is not associated.",ERR,ERROR,*999) ENDIF CALL EXITS("NAVIER_STOKES_EQUATIONS_SET_SUBTYPE_SET") @@ -193,7 +193,7 @@ END SUBROUTINE NAVIER_STOKES_EQUATIONS_SET_SUBTYPE_SET !================================================================================================================================ ! - !>Sets up the standard Stokes fluid setup. + !>Sets up the Navier-Stokes fluid setup. SUBROUTINE NAVIER_STOKES_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*) !Argument variables @@ -221,514 +221,509 @@ SUBROUTINE NAVIER_STOKES_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SETUP,E NULLIFY(EQUATIONS_MATRICES) NULLIFY(GEOMETRIC_DECOMPOSITION) NULLIFY(BOUNDARY_CONDITIONS) - + IF(ASSOCIATED(EQUATIONS_SET)) THEN - SELECT CASE(EQUATIONS_SET%SUBTYPE) + SELECT CASE(EQUATIONS_SET%SUBTYPE) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! NAVIER_STOKES SUBTYPES !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! NAVIER_STOKES SUBTYPES !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,& - & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE, & - & EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE) + CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,& + & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE, & + & EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE) SELECT CASE(EQUATIONS_SET_SETUP%SETUP_TYPE) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! SOLUTION METHOD !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! SOLUTION METHOD !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE(EQUATIONS_SET_SETUP_INITIAL_TYPE) - SELECT CASE(EQUATIONS_SET%SUBTYPE) - CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& - & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) - - SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE) - CASE(EQUATIONS_SET_SETUP_START_ACTION) - CALL NAVIER_STOKES_EQUATIONS_SET_SOLUTION_METHOD_SET(EQUATIONS_SET,& - &EQUATIONS_SET_FEM_SOLUTION_METHOD,ERR,ERROR,*999) - EQUATIONS_SET%SOLUTION_METHOD=EQUATIONS_SET_FEM_SOLUTION_METHOD - CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) - !!TODO: Check valid setup - CASE DEFAULT - LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%ACTION_TYPE,& - &"*",ERR,ERROR))// " for a setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%& - &SETUP_TYPE,"*",ERR,ERROR))// " is invalid for a standard Stokes fluid." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT - - CASE DEFAULT - LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " for a setup sub type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " is invalid for a Stokes equation." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& + & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) + + SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE) + CASE(EQUATIONS_SET_SETUP_START_ACTION) + CALL NAVIER_STOKES_EQUATIONS_SET_SOLUTION_METHOD_SET(EQUATIONS_SET,& + &EQUATIONS_SET_FEM_SOLUTION_METHOD,ERR,ERROR,*999) + EQUATIONS_SET%SOLUTION_METHOD=EQUATIONS_SET_FEM_SOLUTION_METHOD + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) +!!TODO: Check valid setup + CASE DEFAULT + LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%ACTION_TYPE,& + &"*",ERR,ERROR))// " for a setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%& + &SETUP_TYPE,"*",ERR,ERROR))// " is invalid for a Navier-Stokes fluid." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + + CASE DEFAULT + LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " for a setup sub type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " is invalid for a Stokes equation." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! GEOMETRY FIELD !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + CASE(EQUATIONS_SET_SETUP_GEOMETRY_TYPE) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! GEOMETRY FIELD !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& + & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) - CASE(EQUATIONS_SET_SETUP_GEOMETRY_TYPE) + !Do nothing??? - SELECT CASE(EQUATIONS_SET%SUBTYPE) - CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& - & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) + CASE DEFAULT + LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes equation." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT - !Do nothing??? +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! DEPENDENT FIELD !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - CASE DEFAULT - LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " is invalid for a Stokes equation." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT + CASE(EQUATIONS_SET_SETUP_DEPENDENT_TYPE) + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& + & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) + + SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE) + !Set start action + CASE(EQUATIONS_SET_SETUP_START_ACTION) + IF(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN + !Create the auto created dependent field + !start field creation with name 'DEPENDENT_FIELD' + CALL FIELD_CREATE_START(EQUATIONS_SET_SETUP%FIELD_USER_NUMBER,EQUATIONS_SET%REGION,& + &EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,ERR,ERROR,*999) + !start creation of a new field + CALL FIELD_TYPE_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,FIELD_GENERAL_TYPE,ERR,ERROR,*999) + !define new created field to be dependent + CALL FIELD_DEPENDENT_TYPE_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,& + &FIELD_DEPENDENT_TYPE,ERR,ERROR,*999) + !look for decomposition rule already defined + CALL FIELD_MESH_DECOMPOSITION_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,GEOMETRIC_DECOMPOSITION,& + & ERR,ERROR,*999) + !apply decomposition rule found on new created field + CALL FIELD_MESH_DECOMPOSITION_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,& + &GEOMETRIC_DECOMPOSITION,ERR,ERROR,*999) + !point new field to geometric field + CALL FIELD_GEOMETRIC_FIELD_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,EQUATIONS_SET%GEOMETRY% & + & GEOMETRIC_FIELD,ERR,ERROR,*999) + !set number of variables to 2 (1 for U and one for DELUDELN) + DEPENDENT_FIELD_NUMBER_OF_VARIABLES=2 + CALL FIELD_NUMBER_OF_VARIABLES_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD, & + & DEPENDENT_FIELD_NUMBER_OF_VARIABLES,ERR,ERROR,*999) + CALL FIELD_VARIABLE_TYPES_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,(/FIELD_U_VARIABLE_TYPE, & + & FIELD_DELUDELN_VARIABLE_TYPE/),ERR,ERROR,*999) + CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_VECTOR_DIMENSION_TYPE,ERR,ERROR,*999) + CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,FIELD_DELUDELN_VARIABLE_TYPE, & + & FIELD_VECTOR_DIMENSION_TYPE,ERR,ERROR,*999) + CALL FIELD_DATA_TYPE_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_DP_TYPE,ERR,ERROR,*999) + CALL FIELD_DATA_TYPE_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,FIELD_DELUDELN_VARIABLE_TYPE, & + & FIELD_DP_TYPE,ERR,ERROR,*999) + CALL FIELD_NUMBER_OF_COMPONENTS_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE, & + & NUMBER_OF_DIMENSIONS,ERR,ERROR,*999) + + !calculate number of components with one component for each dimension and one for pressure + DEPENDENT_FIELD_NUMBER_OF_COMPONENTS=NUMBER_OF_DIMENSIONS+1 + CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE,& + DEPENDENT_FIELD_NUMBER_OF_COMPONENTS,ERR,ERROR,*999) + CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%DEPENDENT%DEPENDENT_FIELD,& + & FIELD_DELUDELN_VARIABLE_TYPE,DEPENDENT_FIELD_NUMBER_OF_COMPONENTS,ERR,ERROR,*999) + + !set first i components to velocity field (2) and the (i+1)th component to the pressure field (3) + DO I=1,DEPENDENT_FIELD_NUMBER_OF_COMPONENTS + IF(IEQUATIONS_SET%MATERIALS - IF(ASSOCIATED(EQUATIONS_MATERIALS)) THEN - IF(EQUATIONS_MATERIALS%MATERIALS_FIELD_AUTO_CREATED) THEN - !Create the auto created materials field - !start field creation with name 'MATERIAL_FIELD' - CALL FIELD_CREATE_START(EQUATIONS_SET_SETUP%FIELD_USER_NUMBER,EQUATIONS_SET%REGION,EQUATIONS_SET%MATERIALS% & - & MATERIALS_FIELD,ERR,ERROR,*999) - CALL FIELD_TYPE_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_MATERIAL_TYPE,ERR,ERROR,*999) - CALL FIELD_DEPENDENT_TYPE_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_INDEPENDENT_TYPE, & - & ERR,ERROR,*999) - CALL FIELD_MESH_DECOMPOSITION_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,GEOMETRIC_DECOMPOSITION, & - & ERR,ERROR,*999) - !apply decomposition rule found on new created field - CALL FIELD_MESH_DECOMPOSITION_SET_AND_LOCK(EQUATIONS_SET%MATERIALS%MATERIALS_FIELD,GEOMETRIC_DECOMPOSITION, & - & ERR,ERROR,*999) - !point new field to geometric field - CALL FIELD_GEOMETRIC_FIELD_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,EQUATIONS_SET%GEOMETRY% & - & GEOMETRIC_FIELD,ERR,ERROR,*999) - CALL FIELD_NUMBER_OF_VARIABLES_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,MATERIAL_FIELD_NUMBER_OF_VARIABLES,ERR & - &,ERROR,*999) - CALL FIELD_VARIABLE_TYPES_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,(/FIELD_U_VARIABLE_TYPE/), & - & ERR,ERROR,*999) - CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & - & FIELD_VECTOR_DIMENSION_TYPE,ERR,ERROR,*999) - CALL FIELD_DATA_TYPE_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & - & FIELD_DP_TYPE,ERR,ERROR,*999) - CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & - & MATERIAL_FIELD_NUMBER_OF_COMPONENTS,ERR,ERROR,*999) - CALL FIELD_COMPONENT_MESH_COMPONENT_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE, & - & 1,GEOMETRIC_COMPONENT_NUMBER,ERR,ERROR,*999) - CALL FIELD_COMPONENT_MESH_COMPONENT_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & - & 1,GEOMETRIC_COMPONENT_NUMBER,ERR,ERROR,*999) - CALL FIELD_COMPONENT_INTERPOLATION_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & - & 1,FIELD_CONSTANT_INTERPOLATION,ERR,ERROR,*999) - !Default the field scaling to that of the geometric field - CALL FIELD_SCALING_TYPE_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,GEOMETRIC_SCALING_TYPE,ERR,ERROR,*999) - CALL FIELD_SCALING_TYPE_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,GEOMETRIC_SCALING_TYPE,ERR,ERROR,*999) - ELSE - !Check the user specified field - CALL FIELD_TYPE_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_MATERIAL_TYPE,ERR,ERROR,*999) - CALL FIELD_DEPENDENT_TYPE_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_INDEPENDENT_TYPE,ERR,ERROR,*999) - CALL FIELD_NUMBER_OF_VARIABLES_CHECK(EQUATIONS_SET_SETUP%FIELD,1,ERR,ERROR,*999) - CALL FIELD_VARIABLE_TYPES_CHECK(EQUATIONS_SET_SETUP%FIELD,(/FIELD_U_VARIABLE_TYPE/),ERR,ERROR,*999) - CALL FIELD_DIMENSION_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,FIELD_VECTOR_DIMENSION_TYPE, & - & ERR,ERROR,*999) - CALL FIELD_DATA_TYPE_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,FIELD_DP_TYPE,ERR,ERROR,*999) - CALL FIELD_NUMBER_OF_COMPONENTS_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE, & - & NUMBER_OF_DIMENSIONS,ERR,ERROR,*999) - CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,1,ERR,ERROR,*999) - ENDIF - ELSE - CALL FLAG_ERROR("Equations set materials is not associated.",ERR,ERROR,*999) - END IF - - !Specify start action - CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) - EQUATIONS_MATERIALS=>EQUATIONS_SET%MATERIALS - IF(ASSOCIATED(EQUATIONS_MATERIALS)) THEN - IF(EQUATIONS_MATERIALS%MATERIALS_FIELD_AUTO_CREATED) THEN - !Finish creating the materials field - CALL FIELD_CREATE_FINISH(EQUATIONS_MATERIALS%MATERIALS_FIELD,ERR,ERROR,*999) - !Set the default values for the materials field - !First set the mu values to 0.001 - !MATERIAL_FIELD_NUMBER_OF_COMPONENTS - ! viscosity=1 - CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & - & FIELD_VALUES_SET_TYPE,1,1.0_DP,ERR,ERROR,*999) - ! density=2 - CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & - & FIELD_VALUES_SET_TYPE,2,1.0_DP,ERR,ERROR,*999) - ENDIF - ELSE - CALL FLAG_ERROR("Equations set materials is not associated.",ERR,ERROR,*999) - ENDIF - - CASE DEFAULT - LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & - & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & - & " is invalid for Stokes equation." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT - - CASE DEFAULT - LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " for a setup sub type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " is invalid for a Stokes equation." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& + & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) + + !variable X with has Y components, here Y represents viscosity only + MATERIAL_FIELD_NUMBER_OF_VARIABLES=1!X + MATERIAL_FIELD_NUMBER_OF_COMPONENTS=2!Y + + SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE) + !Specify start action + + CASE(EQUATIONS_SET_SETUP_START_ACTION) + EQUATIONS_MATERIALS=>EQUATIONS_SET%MATERIALS + IF(ASSOCIATED(EQUATIONS_MATERIALS)) THEN + IF(EQUATIONS_MATERIALS%MATERIALS_FIELD_AUTO_CREATED) THEN + !Create the auto created materials field + !start field creation with name 'MATERIAL_FIELD' + CALL FIELD_CREATE_START(EQUATIONS_SET_SETUP%FIELD_USER_NUMBER,EQUATIONS_SET%REGION,EQUATIONS_SET%MATERIALS% & + & MATERIALS_FIELD,ERR,ERROR,*999) + CALL FIELD_TYPE_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_MATERIAL_TYPE,ERR,ERROR,*999) + CALL FIELD_DEPENDENT_TYPE_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_INDEPENDENT_TYPE, & + & ERR,ERROR,*999) + CALL FIELD_MESH_DECOMPOSITION_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,GEOMETRIC_DECOMPOSITION, & + & ERR,ERROR,*999) + !apply decomposition rule found on new created field + CALL FIELD_MESH_DECOMPOSITION_SET_AND_LOCK(EQUATIONS_SET%MATERIALS%MATERIALS_FIELD,GEOMETRIC_DECOMPOSITION, & + & ERR,ERROR,*999) + !point new field to geometric field + CALL FIELD_GEOMETRIC_FIELD_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,EQUATIONS_SET%GEOMETRY% & + & GEOMETRIC_FIELD,ERR,ERROR,*999) + CALL FIELD_NUMBER_OF_VARIABLES_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,MATERIAL_FIELD_NUMBER_OF_VARIABLES,ERR & + &,ERROR,*999) + CALL FIELD_VARIABLE_TYPES_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,(/FIELD_U_VARIABLE_TYPE/), & + & ERR,ERROR,*999) + CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_VECTOR_DIMENSION_TYPE,ERR,ERROR,*999) + CALL FIELD_DATA_TYPE_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_DP_TYPE,ERR,ERROR,*999) + CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & MATERIAL_FIELD_NUMBER_OF_COMPONENTS,ERR,ERROR,*999) + CALL FIELD_COMPONENT_MESH_COMPONENT_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE, & + & 1,GEOMETRIC_COMPONENT_NUMBER,ERR,ERROR,*999) + CALL FIELD_COMPONENT_MESH_COMPONENT_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & 1,GEOMETRIC_COMPONENT_NUMBER,ERR,ERROR,*999) + CALL FIELD_COMPONENT_INTERPOLATION_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & 1,FIELD_CONSTANT_INTERPOLATION,ERR,ERROR,*999) + !Default the field scaling to that of the geometric field + CALL FIELD_SCALING_TYPE_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,GEOMETRIC_SCALING_TYPE,ERR,ERROR,*999) + CALL FIELD_SCALING_TYPE_SET(EQUATIONS_MATERIALS%MATERIALS_FIELD,GEOMETRIC_SCALING_TYPE,ERR,ERROR,*999) + ELSE + !Check the user specified field + CALL FIELD_TYPE_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_MATERIAL_TYPE,ERR,ERROR,*999) + CALL FIELD_DEPENDENT_TYPE_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_INDEPENDENT_TYPE,ERR,ERROR,*999) + CALL FIELD_NUMBER_OF_VARIABLES_CHECK(EQUATIONS_SET_SETUP%FIELD,1,ERR,ERROR,*999) + CALL FIELD_VARIABLE_TYPES_CHECK(EQUATIONS_SET_SETUP%FIELD,(/FIELD_U_VARIABLE_TYPE/),ERR,ERROR,*999) + CALL FIELD_DIMENSION_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,FIELD_VECTOR_DIMENSION_TYPE, & + & ERR,ERROR,*999) + CALL FIELD_DATA_TYPE_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,FIELD_DP_TYPE,ERR,ERROR,*999) + CALL FIELD_NUMBER_OF_COMPONENTS_GET(EQUATIONS_SET%GEOMETRY%GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE, & + & NUMBER_OF_DIMENSIONS,ERR,ERROR,*999) + CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET_SETUP%FIELD,FIELD_U_VARIABLE_TYPE,1,ERR,ERROR,*999) + ENDIF + ELSE + CALL FLAG_ERROR("Equations set materials is not associated.",ERR,ERROR,*999) + END IF + + !Specify start action + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) + EQUATIONS_MATERIALS=>EQUATIONS_SET%MATERIALS + IF(ASSOCIATED(EQUATIONS_MATERIALS)) THEN + IF(EQUATIONS_MATERIALS%MATERIALS_FIELD_AUTO_CREATED) THEN + !Finish creating the materials field + CALL FIELD_CREATE_FINISH(EQUATIONS_MATERIALS%MATERIALS_FIELD,ERR,ERROR,*999) + !Set the default values for the materials field + !First set the mu values to 0.001 + !MATERIAL_FIELD_NUMBER_OF_COMPONENTS + ! viscosity=1 + CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_VALUES_SET_TYPE,1,1.0_DP,ERR,ERROR,*999) + ! density=2 + CALL FIELD_COMPONENT_VALUES_INITIALISE(EQUATIONS_MATERIALS%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_VALUES_SET_TYPE,2,1.0_DP,ERR,ERROR,*999) + ENDIF + ELSE + CALL FLAG_ERROR("Equations set materials is not associated.",ERR,ERROR,*999) + ENDIF + CASE DEFAULT + LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & + & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & + & " is invalid for Stokes equation." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + + CASE DEFAULT + LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " for a setup sub type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " is invalid for a Stokes equation." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! SOURCE TYPE !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - CASE(EQUATIONS_SET_SETUP_SOURCE_TYPE) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! SOURCE TYPE !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - SELECT CASE(EQUATIONS_SET%SUBTYPE) - CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& - & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) - - !TO DO: INCLUDE GRAVITY AS SOURCE TYPE - SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE) - CASE(EQUATIONS_SET_SETUP_START_ACTION) - !Do nothing - CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) - !Do nothing - !? Maybe set finished flag???? - CASE DEFAULT - LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & - & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & - & " is invalid for a standard Stokes fluid." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT - - CASE DEFAULT - LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " for a setup sub type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " is invalid for a Stokes equation." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT + CASE(EQUATIONS_SET_SETUP_SOURCE_TYPE) + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& + & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) + + !TO DO: INCLUDE GRAVITY AS SOURCE TYPE + SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE) + CASE(EQUATIONS_SET_SETUP_START_ACTION) + !Do nothing + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) + !Do nothing + !? Maybe set finished flag???? + CASE DEFAULT + LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & + & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes fluid." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + + CASE DEFAULT + LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " for a setup sub type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes equation." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! EQUATIONS_TYPE !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! EQUATIONS_TYPE !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE(EQUATIONS_SET_SETUP_EQUATIONS_TYPE) - SELECT CASE(EQUATIONS_SET%SUBTYPE) - CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE) - - SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE) - CASE(EQUATIONS_SET_SETUP_START_ACTION) - EQUATIONS_MATERIALS=>EQUATIONS_SET%MATERIALS - IF(ASSOCIATED(EQUATIONS_MATERIALS)) THEN - IF(EQUATIONS_MATERIALS%MATERIALS_FINISHED) THEN - CALL EQUATIONS_CREATE_START(EQUATIONS_SET,EQUATIONS,ERR,ERROR,*999) - CALL EQUATIONS_LINEARITY_TYPE_SET(EQUATIONS,EQUATIONS_NONLINEAR,ERR,ERROR,*999) - CALL EQUATIONS_TIME_DEPENDENCE_TYPE_SET(EQUATIONS,EQUATIONS_STATIC,ERR,ERROR,*999) - ELSE - CALL FLAG_ERROR("Equations set materials has not been finished.",ERR,ERROR,*999) - ENDIF - ELSE - CALL FLAG_ERROR("Equations materials is not associated.",ERR,ERROR,*999) - ENDIF - CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) - SELECT CASE(EQUATIONS_SET%SOLUTION_METHOD) - CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD) - !Finish the equations creation - CALL EQUATIONS_SET_EQUATIONS_GET(EQUATIONS_SET,EQUATIONS,ERR,ERROR,*999) - CALL EQUATIONS_CREATE_FINISH(EQUATIONS,ERR,ERROR,*999) - !Create the equations mapping. - CALL EQUATIONS_MAPPING_CREATE_START(EQUATIONS,EQUATIONS_MAPPING,ERR,ERROR,*999) - CALL EQUATIONS_MAPPING_LINEAR_MATRICES_NUMBER_SET(EQUATIONS_MAPPING,1,ERR,ERROR,*999) - CALL EQUATIONS_MAPPING_LINEAR_MATRICES_VARIABLE_TYPES_SET(EQUATIONS_MAPPING,(/FIELD_U_VARIABLE_TYPE/), & - & ERR,ERROR,*999) - CALL EQUATIONS_MAPPING_RHS_VARIABLE_TYPE_SET(EQUATIONS_MAPPING,FIELD_DELUDELN_VARIABLE_TYPE,ERR,ERROR,*999) - CALL EQUATIONS_MAPPING_CREATE_FINISH(EQUATIONS_MAPPING,ERR,ERROR,*999) - !Create the equations matrices - CALL EQUATIONS_MATRICES_CREATE_START(EQUATIONS,EQUATIONS_MATRICES,ERR,ERROR,*999) - - SELECT CASE(EQUATIONS%SPARSITY_TYPE) - CASE(EQUATIONS_MATRICES_FULL_MATRICES) - CALL EQUATIONS_MATRICES_LINEAR_STORAGE_TYPE_SET(EQUATIONS_MATRICES,(/MATRIX_BLOCK_STORAGE_TYPE/), & - & ERR,ERROR,*999) - CALL EQUATIONS_MATRICES_NONLINEAR_STORAGE_TYPE_SET(EQUATIONS_MATRICES,MATRIX_BLOCK_STORAGE_TYPE, & - & ERR,ERROR,*999) - CASE(EQUATIONS_MATRICES_SPARSE_MATRICES) - CALL EQUATIONS_MATRICES_LINEAR_STORAGE_TYPE_SET(EQUATIONS_MATRICES,& - &(/MATRIX_COMPRESSED_ROW_STORAGE_TYPE/),ERR,ERROR,*999) - CALL EQUATIONS_MATRICES_NONLINEAR_STORAGE_TYPE_SET(EQUATIONS_MATRICES,& - &MATRIX_COMPRESSED_ROW_STORAGE_TYPE,ERR,ERROR,*999) - CALL EQUATIONS_MATRICES_LINEAR_STRUCTURE_TYPE_SET(EQUATIONS_MATRICES,(/EQUATIONS_MATRIX_FEM_STRUCTURE/), & - & ERR,ERROR,*999) - CALL EQUATIONS_MATRICES_NONLINEAR_STRUCTURE_TYPE_SET(EQUATIONS_MATRICES,EQUATIONS_MATRIX_FEM_STRUCTURE, & - & ERR,ERROR,*999) - CASE DEFAULT - LOCAL_ERROR="The equations matrices sparsity type of "// & - & TRIM(NUMBER_TO_VSTRING(EQUATIONS%SPARSITY_TYPE,"*",ERR,ERROR))//" is invalid." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT - - CALL EQUATIONS_MATRICES_CREATE_FINISH(EQUATIONS_MATRICES,ERR,ERROR,*999) - CASE DEFAULT - LOCAL_ERROR="The solution method of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SOLUTION_METHOD,& - &"*",ERR,ERROR))//" is invalid." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT - CASE DEFAULT - LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & - & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & - & " is invalid for a steady Laplace equation." + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE) + + SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE) + CASE(EQUATIONS_SET_SETUP_START_ACTION) + EQUATIONS_MATERIALS=>EQUATIONS_SET%MATERIALS + IF(ASSOCIATED(EQUATIONS_MATERIALS)) THEN + IF(EQUATIONS_MATERIALS%MATERIALS_FINISHED) THEN + CALL EQUATIONS_CREATE_START(EQUATIONS_SET,EQUATIONS,ERR,ERROR,*999) + CALL EQUATIONS_LINEARITY_TYPE_SET(EQUATIONS,EQUATIONS_NONLINEAR,ERR,ERROR,*999) + CALL EQUATIONS_TIME_DEPENDENCE_TYPE_SET(EQUATIONS,EQUATIONS_STATIC,ERR,ERROR,*999) + ELSE + CALL FLAG_ERROR("Equations set materials has not been finished.",ERR,ERROR,*999) + ENDIF + ELSE + CALL FLAG_ERROR("Equations materials is not associated.",ERR,ERROR,*999) + ENDIF + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) + SELECT CASE(EQUATIONS_SET%SOLUTION_METHOD) + CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD) + !Finish the equations creation + CALL EQUATIONS_SET_EQUATIONS_GET(EQUATIONS_SET,EQUATIONS,ERR,ERROR,*999) + CALL EQUATIONS_CREATE_FINISH(EQUATIONS,ERR,ERROR,*999) + !Create the equations mapping. + CALL EQUATIONS_MAPPING_CREATE_START(EQUATIONS,EQUATIONS_MAPPING,ERR,ERROR,*999) + CALL EQUATIONS_MAPPING_LINEAR_MATRICES_NUMBER_SET(EQUATIONS_MAPPING,1,ERR,ERROR,*999) + CALL EQUATIONS_MAPPING_LINEAR_MATRICES_VARIABLE_TYPES_SET(EQUATIONS_MAPPING,(/FIELD_U_VARIABLE_TYPE/), & + & ERR,ERROR,*999) + CALL EQUATIONS_MAPPING_RHS_VARIABLE_TYPE_SET(EQUATIONS_MAPPING,FIELD_DELUDELN_VARIABLE_TYPE,ERR,ERROR,*999) + CALL EQUATIONS_MAPPING_CREATE_FINISH(EQUATIONS_MAPPING,ERR,ERROR,*999) + !Create the equations matrices + CALL EQUATIONS_MATRICES_CREATE_START(EQUATIONS,EQUATIONS_MATRICES,ERR,ERROR,*999) + + SELECT CASE(EQUATIONS%SPARSITY_TYPE) + CASE(EQUATIONS_MATRICES_FULL_MATRICES) + CALL EQUATIONS_MATRICES_LINEAR_STORAGE_TYPE_SET(EQUATIONS_MATRICES,(/MATRIX_BLOCK_STORAGE_TYPE/), & + & ERR,ERROR,*999) + CALL EQUATIONS_MATRICES_NONLINEAR_STORAGE_TYPE_SET(EQUATIONS_MATRICES,MATRIX_BLOCK_STORAGE_TYPE, & + & ERR,ERROR,*999) + CASE(EQUATIONS_MATRICES_SPARSE_MATRICES) + CALL EQUATIONS_MATRICES_LINEAR_STORAGE_TYPE_SET(EQUATIONS_MATRICES,& + &(/MATRIX_COMPRESSED_ROW_STORAGE_TYPE/),ERR,ERROR,*999) + CALL EQUATIONS_MATRICES_NONLINEAR_STORAGE_TYPE_SET(EQUATIONS_MATRICES,& + &MATRIX_COMPRESSED_ROW_STORAGE_TYPE,ERR,ERROR,*999) + CALL EQUATIONS_MATRICES_LINEAR_STRUCTURE_TYPE_SET(EQUATIONS_MATRICES,(/EQUATIONS_MATRIX_FEM_STRUCTURE/), & + & ERR,ERROR,*999) + CALL EQUATIONS_MATRICES_NONLINEAR_STRUCTURE_TYPE_SET(EQUATIONS_MATRICES,EQUATIONS_MATRIX_FEM_STRUCTURE, & + & ERR,ERROR,*999) + CASE DEFAULT + LOCAL_ERROR="The equations matrices sparsity type of "// & + & TRIM(NUMBER_TO_VSTRING(EQUATIONS%SPARSITY_TYPE,"*",ERR,ERROR))//" is invalid." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT + END SELECT - CASE(EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) + CALL EQUATIONS_MATRICES_CREATE_FINISH(EQUATIONS_MATRICES,ERR,ERROR,*999) + CASE DEFAULT + LOCAL_ERROR="The solution method of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SOLUTION_METHOD,& + &"*",ERR,ERROR))//" is invalid." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + CASE DEFAULT + LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & + & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & + & " is invalid for a steady Laplace equation." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT - LOCAL_ERROR="Transient Navier-Stokes is not implemented yet!" - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + CASE(EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) - CASE DEFAULT - LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " for a setup sub type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " is invalid for a Stokes equation." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT + LOCAL_ERROR="Transient Navier-Stokes is not implemented yet!" + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! BOUNDARY CONDITIONS TYPE !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + CASE DEFAULT + LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " for a setup sub type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes equation." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! BOUNDARY CONDITIONS TYPE !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE(EQUATIONS_SET_SETUP_BOUNDARY_CONDITIONS_TYPE) - SELECT CASE(EQUATIONS_SET%SUBTYPE) - CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& - & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) - - SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE) - CASE(EQUATIONS_SET_SETUP_START_ACTION) - CALL EQUATIONS_SET_EQUATIONS_GET(EQUATIONS_SET,EQUATIONS,ERR,ERROR,*999) - IF(ASSOCIATED(EQUATIONS)) THEN - IF(EQUATIONS%EQUATIONS_FINISHED) THEN - CALL BOUNDARY_CONDITIONS_CREATE_START(EQUATIONS_SET,BOUNDARY_CONDITIONS,ERR,ERROR,*999) - ELSE - CALL FLAG_ERROR("Equations set equations has not been finished.",ERR,ERROR,*999) - ENDIF - ELSE - CALL FLAG_ERROR("Equations set equations is not associated.",ERR,ERROR,*999) - ENDIF - CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) - CALL EQUATIONS_SET_BOUNDARY_CONDITIONS_GET(EQUATIONS_SET,BOUNDARY_CONDITIONS,ERR,ERROR,*999) - CALL BOUNDARY_CONDITIONS_CREATE_FINISH(BOUNDARY_CONDITIONS,ERR,ERROR,*999) - CASE DEFAULT - LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & - & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & - & " is invalid for a standard Stokes fluid." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT - - CASE DEFAULT - LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " for a setup sub type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " is invalid for a Stokes equation." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& + & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) + + SELECT CASE(EQUATIONS_SET_SETUP%ACTION_TYPE) + CASE(EQUATIONS_SET_SETUP_START_ACTION) + CALL EQUATIONS_SET_EQUATIONS_GET(EQUATIONS_SET,EQUATIONS,ERR,ERROR,*999) + IF(ASSOCIATED(EQUATIONS)) THEN + IF(EQUATIONS%EQUATIONS_FINISHED) THEN + CALL BOUNDARY_CONDITIONS_CREATE_START(EQUATIONS_SET,BOUNDARY_CONDITIONS,ERR,ERROR,*999) + ELSE + CALL FLAG_ERROR("Equations set equations has not been finished.",ERR,ERROR,*999) + ENDIF + ELSE + CALL FLAG_ERROR("Equations set equations is not associated.",ERR,ERROR,*999) + ENDIF + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) + CALL EQUATIONS_SET_BOUNDARY_CONDITIONS_GET(EQUATIONS_SET,BOUNDARY_CONDITIONS,ERR,ERROR,*999) + CALL BOUNDARY_CONDITIONS_CREATE_FINISH(BOUNDARY_CONDITIONS,ERR,ERROR,*999) + CASE DEFAULT + LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & + & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes fluid." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + + CASE DEFAULT + LOCAL_ERROR="The equation set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " for a setup sub type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes equation." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! CASE DEFAULT !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! CASE DEFAULT !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - CASE DEFAULT LOCAL_ERROR="The setup type of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & - & " is invalid for a standard Stokes fluid." + & " is invalid for a Navier-Stokes fluid." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) END SELECT - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! THE END !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! THE END !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - CASE DEFAULT - LOCAL_ERROR="The equations set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " does not equal a standard Stokes fluid subtype." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT + CASE DEFAULT + LOCAL_ERROR="The equations set subtype of "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " does not equal a Navier-Stokes fluid subtype." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT ELSE CALL FLAG_ERROR("Equations set is not associated.",ERR,ERROR,*999) ENDIF - + CALL EXITS("NAVIER_STOKES_EQUATIONS_SET_SETUP") RETURN 999 CALL ERRORS("NAVIER_STOKES_EQUATIONS_SET_SETUP",ERR,ERROR) @@ -748,7 +743,7 @@ END SUBROUTINE NAVIER_STOKES_EQUATIONS_SET_SETUP !================================================================================================================================ ! - !>Sets/changes the problem subtype for a Stokes fluid type . + !>Sets/changes the problem subtype for a Navier-Stokes fluid type . SUBROUTINE NAVIER_STOKES_PROBLEM_SUBTYPE_SET(PROBLEM,PROBLEM_SUBTYPE,ERR,ERROR,*) !Argument variables @@ -764,26 +759,26 @@ SUBROUTINE NAVIER_STOKES_PROBLEM_SUBTYPE_SET(PROBLEM,PROBLEM_SUBTYPE,ERR,ERROR,* IF(ASSOCIATED(PROBLEM)) THEN SELECT CASE(PROBLEM_SUBTYPE) CASE(PROBLEM_STATIC_NAVIER_STOKES_SUBTYPE) - PROBLEM%CLASS=PROBLEM_FLUID_MECHANICS_CLASS - PROBLEM%TYPE=PROBLEM_NAVIER_STOKES_TYPE - PROBLEM%SUBTYPE=PROBLEM_STATIC_NAVIER_STOKES_SUBTYPE + PROBLEM%CLASS=PROBLEM_FLUID_MECHANICS_CLASS + PROBLEM%TYPE=PROBLEM_NAVIER_STOKES_TYPE + PROBLEM%SUBTYPE=PROBLEM_STATIC_NAVIER_STOKES_SUBTYPE CASE(PROBLEM_LAPLACE_NAVIER_STOKES_SUBTYPE) - PROBLEM%CLASS=PROBLEM_FLUID_MECHANICS_CLASS - PROBLEM%TYPE=PROBLEM_NAVIER_STOKES_TYPE - PROBLEM%SUBTYPE=PROBLEM_LAPLACE_NAVIER_STOKES_SUBTYPE + PROBLEM%CLASS=PROBLEM_FLUID_MECHANICS_CLASS + PROBLEM%TYPE=PROBLEM_NAVIER_STOKES_TYPE + PROBLEM%SUBTYPE=PROBLEM_LAPLACE_NAVIER_STOKES_SUBTYPE CASE(PROBLEM_TRANSIENT_NAVIER_STOKES_SUBTYPE) - PROBLEM%CLASS=PROBLEM_FLUID_MECHANICS_CLASS - PROBLEM%TYPE=PROBLEM_NAVIER_STOKES_TYPE - PROBLEM%SUBTYPE=PROBLEM_TRANSIENT_NAVIER_STOKES_SUBTYPE + PROBLEM%CLASS=PROBLEM_FLUID_MECHANICS_CLASS + PROBLEM%TYPE=PROBLEM_NAVIER_STOKES_TYPE + PROBLEM%SUBTYPE=PROBLEM_TRANSIENT_NAVIER_STOKES_SUBTYPE CASE(PROBLEM_OPTIMISED_NAVIER_STOKES_SUBTYPE) - CALL FLAG_ERROR("Not implemented yet.",ERR,ERROR,*999) + CALL FLAG_ERROR("Not implemented yet.",ERR,ERROR,*999) CASE DEFAULT - LOCAL_ERROR="Problem subtype "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SUBTYPE,"*",ERR,ERROR))// & - & " is not valid for a Stokes fluid type of a fluid mechanics problem class." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + LOCAL_ERROR="Problem subtype "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SUBTYPE,"*",ERR,ERROR))// & + & " is not valid for a Navier-Stokes fluid type of a fluid mechanics problem class." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) END SELECT ELSE - CALL FLAG_ERROR("Problem is not associated.",ERR,ERROR,*999) + CALL FLAG_ERROR("Problem is not associated.",ERR,ERROR,*999) ENDIF CALL EXITS("NAVIER_STOKES_PROBLEM_SUBTYPE_SET") @@ -796,12 +791,12 @@ END SUBROUTINE NAVIER_STOKES_PROBLEM_SUBTYPE_SET ! OK !================================================================================================================================ ! - - !>Sets up the Stokes problem. + + !>Sets up the Navier-Stokes problem. SUBROUTINE NAVIER_STOKES_PROBLEM_SETUP(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*) !Argument variables - TYPE(PROBLEM_TYPE), POINTER :: PROBLEM !PROBLEM%CONTROL_LOOP - CALL CONTROL_LOOP_GET(CONTROL_LOOP_ROOT,CONTROL_LOOP_NODE,CONTROL_LOOP,ERR,ERROR,*999) - CALL CONTROL_LOOP_CREATE_FINISH(CONTROL_LOOP,ERR,ERROR,*999) - CASE DEFAULT - LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & - & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & - & " is invalid for a standard Stokes fluid." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT - CASE(PROBLEM_SETUP_SOLVERS_TYPE) - !Get the control loop - CONTROL_LOOP_ROOT=>PROBLEM%CONTROL_LOOP - CALL CONTROL_LOOP_GET(CONTROL_LOOP_ROOT,CONTROL_LOOP_NODE,CONTROL_LOOP,ERR,ERROR,*999) - SELECT CASE(PROBLEM_SETUP%ACTION_TYPE) - CASE(PROBLEM_SETUP_START_ACTION) - !Start the solvers creation - CALL SOLVERS_CREATE_START(CONTROL_LOOP,SOLVERS,ERR,ERROR,*999) - CALL SOLVERS_NUMBER_SET(SOLVERS,1,ERR,ERROR,*999) - !Set the solver to be a nonlinear solver - CALL SOLVERS_SOLVER_GET(SOLVERS,1,SOLVER,ERR,ERROR,*999) - CALL SOLVER_TYPE_SET(SOLVER,SOLVER_NONLINEAR_TYPE,ERR,ERROR,*999) - !Set solver defaults - CALL SOLVER_LIBRARY_TYPE_SET(SOLVER,SOLVER_PETSC_LIBRARY,ERR,ERROR,*999) - CASE(PROBLEM_SETUP_FINISH_ACTION) - !Get the solvers - CALL CONTROL_LOOP_SOLVERS_GET(CONTROL_LOOP,SOLVERS,ERR,ERROR,*999) - !Finish the solvers creation - CALL SOLVERS_CREATE_FINISH(SOLVERS,ERR,ERROR,*999) - CASE DEFAULT - LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & - & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & - & " is invalid for a standard Stokes fluid." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT - CASE(PROBLEM_SETUP_SOLVER_EQUATIONS_TYPE) - SELECT CASE(PROBLEM_SETUP%ACTION_TYPE) - CASE(PROBLEM_SETUP_START_ACTION) - !Get the control loop - CONTROL_LOOP_ROOT=>PROBLEM%CONTROL_LOOP - CALL CONTROL_LOOP_GET(CONTROL_LOOP_ROOT,CONTROL_LOOP_NODE,CONTROL_LOOP,ERR,ERROR,*999) - !Get the solver - CALL CONTROL_LOOP_SOLVERS_GET(CONTROL_LOOP,SOLVERS,ERR,ERROR,*999) - CALL SOLVERS_SOLVER_GET(SOLVERS,1,SOLVER,ERR,ERROR,*999) - !Create the solver equations - CALL SOLVER_EQUATIONS_CREATE_START(SOLVER,SOLVER_EQUATIONS,ERR,ERROR,*999) - CALL SOLVER_EQUATIONS_LINEARITY_TYPE_SET(SOLVER_EQUATIONS,SOLVER_EQUATIONS_NONLINEAR,ERR,ERROR,*999) - CALL SOLVER_EQUATIONS_TIME_DEPENDENCE_TYPE_SET(SOLVER_EQUATIONS,SOLVER_EQUATIONS_STATIC,ERR,ERROR,*999) - CALL SOLVER_EQUATIONS_SPARSITY_TYPE_SET(SOLVER_EQUATIONS,SOLVER_SPARSE_MATRICES,ERR,ERROR,*999) - CASE(PROBLEM_SETUP_FINISH_ACTION) - !Get the control loop - CONTROL_LOOP_ROOT=>PROBLEM%CONTROL_LOOP - CALL CONTROL_LOOP_GET(CONTROL_LOOP_ROOT,CONTROL_LOOP_NODE,CONTROL_LOOP,ERR,ERROR,*999) - !Get the solver equations - CALL CONTROL_LOOP_SOLVERS_GET(CONTROL_LOOP,SOLVERS,ERR,ERROR,*999) - CALL SOLVERS_SOLVER_GET(SOLVERS,1,SOLVER,ERR,ERROR,*999) - CALL SOLVER_SOLVER_EQUATIONS_GET(SOLVER,SOLVER_EQUATIONS,ERR,ERROR,*999) - !Finish the solver equations creation - CALL SOLVER_EQUATIONS_CREATE_FINISH(SOLVER_EQUATIONS,ERR,ERROR,*999) - CASE DEFAULT - LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & - & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & - & " is invalid for a standard Stokes fluid." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT - CASE DEFAULT - LOCAL_ERROR="The setup type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & - & " is invalid for a standard Stokes fluid." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT + SELECT CASE(PROBLEM_SETUP%SETUP_TYPE) + CASE(PROBLEM_SETUP_INITIAL_TYPE) + SELECT CASE(PROBLEM_SETUP%ACTION_TYPE) + CASE(PROBLEM_SETUP_START_ACTION) + !Do nothing???? + CASE(PROBLEM_SETUP_FINISH_ACTION) + !Do nothing??? + CASE DEFAULT + LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & + & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes fluid." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + CASE(PROBLEM_SETUP_CONTROL_TYPE) + SELECT CASE(PROBLEM_SETUP%ACTION_TYPE) + CASE(PROBLEM_SETUP_START_ACTION) + !Set up a simple control loop + CALL CONTROL_LOOP_CREATE_START(PROBLEM,CONTROL_LOOP,ERR,ERROR,*999) + CASE(PROBLEM_SETUP_FINISH_ACTION) + !Finish the control loops + CONTROL_LOOP_ROOT=>PROBLEM%CONTROL_LOOP + CALL CONTROL_LOOP_GET(CONTROL_LOOP_ROOT,CONTROL_LOOP_NODE,CONTROL_LOOP,ERR,ERROR,*999) + CALL CONTROL_LOOP_CREATE_FINISH(CONTROL_LOOP,ERR,ERROR,*999) + CASE DEFAULT + LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & + & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes fluid." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + CASE(PROBLEM_SETUP_SOLVERS_TYPE) + !Get the control loop + CONTROL_LOOP_ROOT=>PROBLEM%CONTROL_LOOP + CALL CONTROL_LOOP_GET(CONTROL_LOOP_ROOT,CONTROL_LOOP_NODE,CONTROL_LOOP,ERR,ERROR,*999) + SELECT CASE(PROBLEM_SETUP%ACTION_TYPE) + CASE(PROBLEM_SETUP_START_ACTION) + !Start the solvers creation + CALL SOLVERS_CREATE_START(CONTROL_LOOP,SOLVERS,ERR,ERROR,*999) + CALL SOLVERS_NUMBER_SET(SOLVERS,1,ERR,ERROR,*999) + !Set the solver to be a nonlinear solver + CALL SOLVERS_SOLVER_GET(SOLVERS,1,SOLVER,ERR,ERROR,*999) + CALL SOLVER_TYPE_SET(SOLVER,SOLVER_NONLINEAR_TYPE,ERR,ERROR,*999) + !Set solver defaults + CALL SOLVER_LIBRARY_TYPE_SET(SOLVER,SOLVER_PETSC_LIBRARY,ERR,ERROR,*999) + CASE(PROBLEM_SETUP_FINISH_ACTION) + !Get the solvers + CALL CONTROL_LOOP_SOLVERS_GET(CONTROL_LOOP,SOLVERS,ERR,ERROR,*999) + !Finish the solvers creation + CALL SOLVERS_CREATE_FINISH(SOLVERS,ERR,ERROR,*999) + CASE DEFAULT + LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & + & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes fluid." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + CASE(PROBLEM_SETUP_SOLVER_EQUATIONS_TYPE) + SELECT CASE(PROBLEM_SETUP%ACTION_TYPE) + CASE(PROBLEM_SETUP_START_ACTION) + !Get the control loop + CONTROL_LOOP_ROOT=>PROBLEM%CONTROL_LOOP + CALL CONTROL_LOOP_GET(CONTROL_LOOP_ROOT,CONTROL_LOOP_NODE,CONTROL_LOOP,ERR,ERROR,*999) + !Get the solver + CALL CONTROL_LOOP_SOLVERS_GET(CONTROL_LOOP,SOLVERS,ERR,ERROR,*999) + CALL SOLVERS_SOLVER_GET(SOLVERS,1,SOLVER,ERR,ERROR,*999) + !Create the solver equations + CALL SOLVER_EQUATIONS_CREATE_START(SOLVER,SOLVER_EQUATIONS,ERR,ERROR,*999) + CALL SOLVER_EQUATIONS_LINEARITY_TYPE_SET(SOLVER_EQUATIONS,SOLVER_EQUATIONS_NONLINEAR,ERR,ERROR,*999) + CALL SOLVER_EQUATIONS_TIME_DEPENDENCE_TYPE_SET(SOLVER_EQUATIONS,SOLVER_EQUATIONS_STATIC,ERR,ERROR,*999) + CALL SOLVER_EQUATIONS_SPARSITY_TYPE_SET(SOLVER_EQUATIONS,SOLVER_SPARSE_MATRICES,ERR,ERROR,*999) + CASE(PROBLEM_SETUP_FINISH_ACTION) + !Get the control loop + CONTROL_LOOP_ROOT=>PROBLEM%CONTROL_LOOP + CALL CONTROL_LOOP_GET(CONTROL_LOOP_ROOT,CONTROL_LOOP_NODE,CONTROL_LOOP,ERR,ERROR,*999) + !Get the solver equations + CALL CONTROL_LOOP_SOLVERS_GET(CONTROL_LOOP,SOLVERS,ERR,ERROR,*999) + CALL SOLVERS_SOLVER_GET(SOLVERS,1,SOLVER,ERR,ERROR,*999) + CALL SOLVER_SOLVER_EQUATIONS_GET(SOLVER,SOLVER_EQUATIONS,ERR,ERROR,*999) + !Finish the solver equations creation + CALL SOLVER_EQUATIONS_CREATE_FINISH(SOLVER_EQUATIONS,ERR,ERROR,*999) + CASE DEFAULT + LOCAL_ERROR="The action type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%ACTION_TYPE,"*",ERR,ERROR))// & + & " for a setup type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes fluid." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + CASE DEFAULT + LOCAL_ERROR="The setup type of "//TRIM(NUMBER_TO_VSTRING(PROBLEM_SETUP%SETUP_TYPE,"*",ERR,ERROR))// & + & " is invalid for a Navier-Stokes fluid." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! PROBLEM_TRANSIENT_NAVIER_STOKES_SUBTYPE !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! PROBLEM_TRANSIENT_NAVIER_STOKES_SUBTYPE !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE(PROBLEM_TRANSIENT_NAVIER_STOKES_SUBTYPE) - LOCAL_ERROR="Transient Navier-Stokes is not implemented yet!" - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + LOCAL_ERROR="Transient Navier-Stokes is not implemented yet!" + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! DEFAULT !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! DEFAULT !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE DEFAULT LOCAL_ERROR="Problem subtype "//TRIM(NUMBER_TO_VSTRING(PROBLEM%SUBTYPE,"*",ERR,ERROR))// & - & " is not valid for a Stokes fluid type of a fluid mechanics problem class." + & " is not valid for a Navier-Stokes equation type of a fluid mechanics problem class." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) END SELECT ELSE @@ -998,7 +993,7 @@ SUBROUTINE NAVIER_STOKES_FINITE_ELEMENT_RESIDUAL_EVALUATE(EQUATIONS_SET,ELEMENT_ DOUBLE PRECISION:: U_VALUE(3) DOUBLE PRECISION:: U_DERIV(3,3) - + CALL ENTERS("NAVIER_STOKES_FINITE_ELEMENT_RESIDUAL_EVALUATE",ERR,ERROR,*999) out=0 @@ -1015,310 +1010,308 @@ SUBROUTINE NAVIER_STOKES_FINITE_ELEMENT_RESIDUAL_EVALUATE(EQUATIONS_SET,ELEMENT_ EQUATIONS=>EQUATIONS_SET%EQUATIONS IF(ASSOCIATED(EQUATIONS)) THEN - SELECT CASE(EQUATIONS_SET%SUBTYPE) - CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& - & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE,& + & EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! SET POINTERS !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - DEPENDENT_FIELD=>EQUATIONS%INTERPOLATION%DEPENDENT_FIELD - GEOMETRIC_FIELD=>EQUATIONS%INTERPOLATION%GEOMETRIC_FIELD - MATERIALS_FIELD=>EQUATIONS%INTERPOLATION%MATERIALS_FIELD - EQUATIONS_MATRICES=>EQUATIONS%EQUATIONS_MATRICES - GEOMETRIC_BASIS=>GEOMETRIC_FIELD%DECOMPOSITION%DOMAIN(GEOMETRIC_FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% & - & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS - DEPENDENT_BASIS=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(DEPENDENT_FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% & - & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS - QUADRATURE_SCHEME=>DEPENDENT_BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR - RHS_VECTOR=>EQUATIONS_MATRICES%RHS_VECTOR - EQUATIONS_MAPPING=>EQUATIONS%EQUATIONS_MAPPING - - SELECT CASE(EQUATIONS_SET%SUBTYPE) - CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE) - LINEAR_MATRICES=>EQUATIONS_MATRICES%LINEAR_MATRICES - NONLINEAR_MATRICES=>EQUATIONS_MATRICES%NONLINEAR_MATRICES - STIFFNESS_MATRIX=>LINEAR_MATRICES%MATRICES(1)%PTR - LINEAR_MAPPING=>EQUATIONS_MAPPING%LINEAR_MAPPING - NONLINEAR_MAPPING=>EQUATIONS_MAPPING%NONLINEAR_MAPPING -! FIELD_VARIABLE=>LINEAR_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE - FIELD_VARIABLE=>NONLINEAR_MAPPING%RESIDUAL_VARIABLE -! SOURCE_VECTOR=>EQUATIONS_MATRICES%SOURCE_VECTOR - - STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX=0.0_DP - NONLINEAR_MATRICES%ELEMENT_RESIDUAL%VECTOR=0.0_DP - - CASE(EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) - - LOCAL_ERROR="Transient Navier-Stokes is not implemented yet!" - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + DEPENDENT_FIELD=>EQUATIONS%INTERPOLATION%DEPENDENT_FIELD + GEOMETRIC_FIELD=>EQUATIONS%INTERPOLATION%GEOMETRIC_FIELD + MATERIALS_FIELD=>EQUATIONS%INTERPOLATION%MATERIALS_FIELD + EQUATIONS_MATRICES=>EQUATIONS%EQUATIONS_MATRICES + GEOMETRIC_BASIS=>GEOMETRIC_FIELD%DECOMPOSITION%DOMAIN(GEOMETRIC_FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% & + & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS + DEPENDENT_BASIS=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(DEPENDENT_FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% & + & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS + QUADRATURE_SCHEME=>DEPENDENT_BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR + RHS_VECTOR=>EQUATIONS_MATRICES%RHS_VECTOR + EQUATIONS_MAPPING=>EQUATIONS%EQUATIONS_MAPPING + + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_NAVIER_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE) + LINEAR_MATRICES=>EQUATIONS_MATRICES%LINEAR_MATRICES + NONLINEAR_MATRICES=>EQUATIONS_MATRICES%NONLINEAR_MATRICES + STIFFNESS_MATRIX=>LINEAR_MATRICES%MATRICES(1)%PTR + LINEAR_MAPPING=>EQUATIONS_MAPPING%LINEAR_MAPPING + NONLINEAR_MAPPING=>EQUATIONS_MAPPING%NONLINEAR_MAPPING + ! FIELD_VARIABLE=>LINEAR_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE + FIELD_VARIABLE=>NONLINEAR_MAPPING%RESIDUAL_VARIABLE + ! SOURCE_VECTOR=>EQUATIONS_MATRICES%SOURCE_VECTOR + + STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX=0.0_DP + NONLINEAR_MATRICES%ELEMENT_RESIDUAL%VECTOR=0.0_DP + + CASE(EQUATIONS_SET_TRANSIENT_NAVIER_STOKES_SUBTYPE) + + LOCAL_ERROR="Transient Navier-Stokes is not implemented yet!" + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - CASE DEFAULT - LOCAL_ERROR="Equations set subtype "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " is not valid for a Stokes fluid type of a fluid mechanics equations set class." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT + CASE DEFAULT + LOCAL_ERROR="Equations set subtype "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " is not valid for a Stokes fluid type of a fluid mechanics equations set class." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,EQUATIONS%INTERPOLATION% & - & DEPENDENT_INTERP_PARAMETERS,ERR,ERROR,*999) + & DEPENDENT_INTERP_PARAMETERS,ERR,ERROR,*999) CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,EQUATIONS%INTERPOLATION% & - & GEOMETRIC_INTERP_PARAMETERS,ERR,ERROR,*999) + & GEOMETRIC_INTERP_PARAMETERS,ERR,ERROR,*999) CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,EQUATIONS%INTERPOLATION% & - & MATERIALS_INTERP_PARAMETERS,ERR,ERROR,*999) + & MATERIALS_INTERP_PARAMETERS,ERR,ERROR,*999) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! LOOP OVER GAUSS POINTS !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - DO ng=1,QUADRATURE_SCHEME%NUMBER_OF_GAUSS - CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,EQUATIONS%INTERPOLATION% & + + DO ng=1,QUADRATURE_SCHEME%NUMBER_OF_GAUSS + CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,EQUATIONS%INTERPOLATION% & & DEPENDENT_INTERP_POINT,ERR,ERROR,*999) - CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,EQUATIONS%INTERPOLATION% & + CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,EQUATIONS%INTERPOLATION% & & GEOMETRIC_INTERP_POINT,ERR,ERROR,*999) - CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(GEOMETRIC_BASIS%NUMBER_OF_XI,EQUATIONS%INTERPOLATION% & + CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(GEOMETRIC_BASIS%NUMBER_OF_XI,EQUATIONS%INTERPOLATION% & & GEOMETRIC_INTERP_POINT_METRICS,ERR,ERROR,*999) - CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,EQUATIONS%INTERPOLATION% & + CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,EQUATIONS%INTERPOLATION% & & MATERIALS_INTERP_POINT,ERR,ERROR,*999) - !Define MU_PARAM, viscosity=1 - MU_PARAM=EQUATIONS%INTERPOLATION%MATERIALS_INTERP_POINT%VALUES(1,NO_PART_DERIV) - !Define RHO_PARAM, density=2 - RHO_PARAM=EQUATIONS%INTERPOLATION%MATERIALS_INTERP_POINT%VALUES(2,NO_PART_DERIV) + !Define MU_PARAM, viscosity=1 + MU_PARAM=EQUATIONS%INTERPOLATION%MATERIALS_INTERP_POINT%VALUES(1,NO_PART_DERIV) + !Define RHO_PARAM, density=2 + RHO_PARAM=EQUATIONS%INTERPOLATION%MATERIALS_INTERP_POINT%VALUES(2,NO_PART_DERIV) - IF(STIFFNESS_MATRIX%FIRST_ASSEMBLY) THEN - IF(STIFFNESS_MATRIX%UPDATE_MATRIX) THEN + IF(STIFFNESS_MATRIX%FIRST_ASSEMBLY) THEN + IF(STIFFNESS_MATRIX%UPDATE_MATRIX) THEN - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! CALCULATE PARTIAL MATRICES !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! CALCULATE PARTIAL MATRICES !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - mhs=0 - DO mh=1,(FIELD_VARIABLE%NUMBER_OF_COMPONENTS-1) + mhs=0 + DO mh=1,(FIELD_VARIABLE%NUMBER_OF_COMPONENTS-1) MESH_COMPONENT1=FIELD_VARIABLE%COMPONENTS(mh)%MESH_COMPONENT_NUMBER DEPENDENT_BASIS1=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT1)%PTR% & - & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS + & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS QUADRATURE_SCHEME1=>DEPENDENT_BASIS1%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR JGW=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN*QUADRATURE_SCHEME1%GAUSS_WEIGHTS(ng) DO ms=1,DEPENDENT_BASIS1%NUMBER_OF_ELEMENT_PARAMETERS - mhs=mhs+1 - nhs=0 - - !Loop over element columns - DO nh=1,(FIELD_VARIABLE%NUMBER_OF_COMPONENTS) - MESH_COMPONENT2=FIELD_VARIABLE%COMPONENTS(nh)%MESH_COMPONENT_NUMBER - DEPENDENT_BASIS2=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT2)%PTR% & - & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS - QUADRATURE_SCHEME2=>DEPENDENT_BASIS2%QUADRATURE%QUADRATURE_SCHEME_MAP& - &(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR - ! JGW=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN*QUADRATURE_SCHEME2%& - ! &GAUSS_WEIGHTS(ng) + mhs=mhs+1 + nhs=0 + + !Loop over element columns + DO nh=1,(FIELD_VARIABLE%NUMBER_OF_COMPONENTS) + MESH_COMPONENT2=FIELD_VARIABLE%COMPONENTS(nh)%MESH_COMPONENT_NUMBER + DEPENDENT_BASIS2=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT2)%PTR% & + & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS + QUADRATURE_SCHEME2=>DEPENDENT_BASIS2%QUADRATURE%QUADRATURE_SCHEME_MAP& + &(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR + ! JGW=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN*QUADRATURE_SCHEME2%& + ! &GAUSS_WEIGHTS(ng) + + DO ns=1,DEPENDENT_BASIS2%NUMBER_OF_ELEMENT_PARAMETERS + nhs=nhs+1 - DO ns=1,DEPENDENT_BASIS2%NUMBER_OF_ELEMENT_PARAMETERS - nhs=nhs+1 - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! PRECONDITIONS !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! PRECONDITIONS !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI - DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI - DXI_DX(mi,ni)=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%DXI_DX(mi,ni) - END DO + DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI + DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI + DXI_DX(mi,ni)=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%DXI_DX(mi,ni) + END DO - DPHIMS_DXI(ni)=QUADRATURE_SCHEME1%GAUSS_BASIS_FNS(ms,PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(ni),ng) - DPHINS_DXI(ni)=QUADRATURE_SCHEME2%GAUSS_BASIS_FNS(ns,PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(ni),ng) - END DO !ni + DPHIMS_DXI(ni)=QUADRATURE_SCHEME1%GAUSS_BASIS_FNS(ms,PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(ni),ng) + DPHINS_DXI(ni)=QUADRATURE_SCHEME2%GAUSS_BASIS_FNS(ns,PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(ni),ng) + END DO !ni - PHIMS=QUADRATURE_SCHEME1%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng) - PHINS=QUADRATURE_SCHEME2%GAUSS_BASIS_FNS(ns,NO_PART_DERIV,ng) + PHIMS=QUADRATURE_SCHEME1%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng) + PHINS=QUADRATURE_SCHEME2%GAUSS_BASIS_FNS(ns,NO_PART_DERIV,ng) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! A LAPLACE ONLY !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! A LAPLACE ONLY !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !LAPLACE TYPE - IF(nh==mh) THEN - SUM=0.0_DP - !Calculate SUM - DO x=1,DEPENDENT_BASIS1%NUMBER_OF_XI - DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI - DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI - SUM=SUM+MU_PARAM*DPHINS_DXI(ni)*DXI_DX(ni,x)*DPHIMS_DXI(mi)*DXI_DX(mi,x) - ENDDO !ni - ENDDO !mi - ENDDO !x - !Calculate MATRIX - AL_MATRIX(mhs,nhs)=AL_MATRIX(mhs,nhs)+SUM*JGW - END IF + !LAPLACE TYPE + IF(nh==mh) THEN + SUM=0.0_DP + !Calculate SUM + DO x=1,DEPENDENT_BASIS1%NUMBER_OF_XI + DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI + DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI + SUM=SUM+MU_PARAM*DPHINS_DXI(ni)*DXI_DX(ni,x)*DPHIMS_DXI(mi)*DXI_DX(mi,x) + ENDDO !ni + ENDDO !mi + ENDDO !x + !Calculate MATRIX + AL_MATRIX(mhs,nhs)=AL_MATRIX(mhs,nhs)+SUM*JGW + END IF - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! A STANDARD !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! A STANDARD !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !GRADIENT TRANSPOSE TYPE - IF(EQUATIONS_SET%SUBTYPE/=EQUATIONS_SET_LAPLACE_NAVIER_STOKES_SUBTYPE) THEN - IF(nhDEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT1)%PTR% & + MESH_COMPONENT1=FIELD_VARIABLE%COMPONENTS(mh)%MESH_COMPONENT_NUMBER + DEPENDENT_BASIS1=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT1)%PTR% & & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS - QUADRATURE_SCHEME1=>DEPENDENT_BASIS1%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR - JGW=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN*QUADRATURE_SCHEME1%GAUSS_WEIGHTS(ng) + QUADRATURE_SCHEME1=>DEPENDENT_BASIS1%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR + JGW=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN*QUADRATURE_SCHEME1%GAUSS_WEIGHTS(ng) - DO ms=1,DEPENDENT_BASIS1%NUMBER_OF_ELEMENT_PARAMETERS - mhs=mhs+1 + DO ms=1,DEPENDENT_BASIS1%NUMBER_OF_ELEMENT_PARAMETERS + mhs=mhs+1 - PHIMS=QUADRATURE_SCHEME1%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng) + PHIMS=QUADRATURE_SCHEME1%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng) - !note mh value derivative - SUM=0.0_DP - SUM=RHO_PARAM*(U_VALUE(1)*U_DERIV(mh,1)+U_VALUE(2)*U_DERIV(mh,2)+U_VALUE(3)*U_DERIV(mh,3))*PHIMS - NL_VECTOR(mhs)=NL_VECTOR(mhs)+SUM*JGW + !note mh value derivative + SUM=0.0_DP + SUM=RHO_PARAM*(U_VALUE(1)*U_DERIV(mh,1)+U_VALUE(2)*U_DERIV(mh,2)+U_VALUE(3)*U_DERIV(mh,3))*PHIMS + NL_VECTOR(mhs)=NL_VECTOR(mhs)+SUM*JGW - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! GO TO MATRIX ASSEMBLY !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ENDDO !ms - ENDDO !mh - END IF +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! GO TO MATRIX ASSEMBLY !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ENDDO !ms + ENDDO !mh + END IF - ENDDO !ng - + ENDDO !ng !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! ASSEMBLE STIFFNESS AND VECTORS !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + mhs_min=mhs + mhs_max=nhs + nhs_min=mhs + nhs_max=nhs - mhs_min=mhs - mhs_max=nhs - nhs_min=mhs - nhs_max=nhs - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! STIFFNESS MATRIX !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - IF(STIFFNESS_MATRIX%FIRST_ASSEMBLY) THEN - IF(STIFFNESS_MATRIX%UPDATE_MATRIX) THEN - - STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(1:mhs_min,1:nhs_min)=(AL_MATRIX(1:mhs_min,1:nhs_min)+AG_MATRIX(1:mhs_min,1:nhs_min)) - STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(1:mhs_min,nhs_min+1:nhs_max)=(BT_MATRIX(1:mhs_min,nhs_min+1:nhs_max)) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! STIFFNESS MATRIX !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IF(STIFFNESS_MATRIX%FIRST_ASSEMBLY) THEN + IF(STIFFNESS_MATRIX%UPDATE_MATRIX) THEN + STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(1:mhs_min,1:nhs_min)=(AL_MATRIX(1:mhs_min,1:nhs_min)+ & + & AG_MATRIX(1:mhs_min,1:nhs_min)) + STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(1:mhs_min,nhs_min+1:nhs_max)=(BT_MATRIX(1:mhs_min,nhs_min+1:nhs_max)) - DO mhs=mhs_min+1,mhs_max - DO nhs=1,nhs_min - !Transpose pressure type entries for mass equation - STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)=STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(nhs,mhs) - END DO - END DO + DO mhs=mhs_min+1,mhs_max + DO nhs=1,nhs_min + !Transpose pressure type entries for mass equation + STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)=STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(nhs,mhs) + END DO + END DO + ENDIF ENDIF - ENDIF - ! BUT WHY??? - STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX=STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX + ! BUT WHY??? + STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX=STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! RIGHT HAND SIDE VECTOR !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - IF(RHS_VECTOR%FIRST_ASSEMBLY) THEN - IF(RHS_VECTOR%UPDATE_VECTOR) THEN - RHS_VECTOR%ELEMENT_VECTOR%VECTOR(1:mhs_max)=RH_VECTOR(1:mhs_max) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! RIGHT HAND SIDE VECTOR !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IF(RHS_VECTOR%FIRST_ASSEMBLY) THEN + IF(RHS_VECTOR%UPDATE_VECTOR) THEN + RHS_VECTOR%ELEMENT_VECTOR%VECTOR(1:mhs_max)=RH_VECTOR(1:mhs_max) + ENDIF ENDIF - ENDIF - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! NONLINEAR VECTOR !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - IF(NONLINEAR_MATRICES%UPDATE_RESIDUAL) THEN - NONLINEAR_MATRICES%ELEMENT_RESIDUAL%VECTOR(1:mhs_max)=NL_VECTOR(1:mhs_max) - END IF +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! NONLINEAR VECTOR !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IF(NONLINEAR_MATRICES%UPDATE_RESIDUAL) THEN + NONLINEAR_MATRICES%ELEMENT_RESIDUAL%VECTOR(1:mhs_max)=NL_VECTOR(1:mhs_max) + END IF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/src/Stokes_equations_routines.f90 b/src/Stokes_equations_routines.f90 index b0f1d9f9..60960d60 100755 --- a/src/Stokes_equations_routines.f90 +++ b/src/Stokes_equations_routines.f90 @@ -1100,11 +1100,11 @@ END SUBROUTINE STOKES_EQUATION_PROBLEM_SETUP ! !================================================================================================================================ ! - - + + !>Calculates the element stiffness matrices and RHS for a Stokes fluid finite element equations set. SUBROUTINE STOKES_EQUATION_FINITE_ELEMENT_CALCULATE(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*) - + !Argument variables TYPE(EQUATIONS_SET_TYPE), POINTER :: EQUATIONS_SET !EQUATIONS_SET%EQUATIONS IF(ASSOCIATED(EQUATIONS)) THEN - SELECT CASE(EQUATIONS_SET%SUBTYPE) - CASE(EQUATIONS_SET_STATIC_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_STOKES_SUBTYPE,& - & EQUATIONS_SET_TRANSIENT_STOKES_SUBTYPE) + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_STOKES_SUBTYPE,& + & EQUATIONS_SET_TRANSIENT_STOKES_SUBTYPE) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! SET POINTERS !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - DEPENDENT_FIELD=>EQUATIONS%INTERPOLATION%DEPENDENT_FIELD - GEOMETRIC_FIELD=>EQUATIONS%INTERPOLATION%GEOMETRIC_FIELD - MATERIALS_FIELD=>EQUATIONS%INTERPOLATION%MATERIALS_FIELD - EQUATIONS_MATRICES=>EQUATIONS%EQUATIONS_MATRICES - GEOMETRIC_BASIS=>GEOMETRIC_FIELD%DECOMPOSITION%DOMAIN(GEOMETRIC_FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% & - & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS - DEPENDENT_BASIS=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(DEPENDENT_FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% & - & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS - QUADRATURE_SCHEME=>DEPENDENT_BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR - RHS_VECTOR=>EQUATIONS_MATRICES%RHS_VECTOR - EQUATIONS_MAPPING=>EQUATIONS%EQUATIONS_MAPPING - SELECT CASE(EQUATIONS_SET%SUBTYPE) - CASE(EQUATIONS_SET_STATIC_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_STOKES_SUBTYPE) - LINEAR_MATRICES=>EQUATIONS_MATRICES%LINEAR_MATRICES - STIFFNESS_MATRIX=>LINEAR_MATRICES%MATRICES(1)%PTR - LINEAR_MAPPING=>EQUATIONS_MAPPING%LINEAR_MAPPING - FIELD_VARIABLE=>LINEAR_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE - - STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX=0.0_DP - CASE(EQUATIONS_SET_TRANSIENT_STOKES_SUBTYPE) - DYNAMIC_MATRICES=>EQUATIONS_MATRICES%DYNAMIC_MATRICES - STIFFNESS_MATRIX=>DYNAMIC_MATRICES%MATRICES(1)%PTR - DAMPING_MATRIX=>DYNAMIC_MATRICES%MATRICES(2)%PTR - DYNAMIC_MAPPING=>EQUATIONS_MAPPING%DYNAMIC_MAPPING - FIELD_VARIABLE=>DYNAMIC_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE - - STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX=0.0_DP - DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX=0.0_DP - CASE DEFAULT - LOCAL_ERROR="Equations set subtype "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & - & " is not valid for a Stokes fluid type of a fluid mechanics equations set class." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - END SELECT - CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,EQUATIONS%INTERPOLATION% & - & GEOMETRIC_INTERP_PARAMETERS,ERR,ERROR,*999) - CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,EQUATIONS%INTERPOLATION% & - & MATERIALS_INTERP_PARAMETERS,ERR,ERROR,*999) + DEPENDENT_FIELD=>EQUATIONS%INTERPOLATION%DEPENDENT_FIELD + GEOMETRIC_FIELD=>EQUATIONS%INTERPOLATION%GEOMETRIC_FIELD + MATERIALS_FIELD=>EQUATIONS%INTERPOLATION%MATERIALS_FIELD + EQUATIONS_MATRICES=>EQUATIONS%EQUATIONS_MATRICES + GEOMETRIC_BASIS=>GEOMETRIC_FIELD%DECOMPOSITION%DOMAIN(GEOMETRIC_FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% & + & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS + DEPENDENT_BASIS=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(DEPENDENT_FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% & + & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS + QUADRATURE_SCHEME=>DEPENDENT_BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR + RHS_VECTOR=>EQUATIONS_MATRICES%RHS_VECTOR + EQUATIONS_MAPPING=>EQUATIONS%EQUATIONS_MAPPING + SELECT CASE(EQUATIONS_SET%SUBTYPE) + CASE(EQUATIONS_SET_STATIC_STOKES_SUBTYPE,EQUATIONS_SET_LAPLACE_STOKES_SUBTYPE) + LINEAR_MATRICES=>EQUATIONS_MATRICES%LINEAR_MATRICES + STIFFNESS_MATRIX=>LINEAR_MATRICES%MATRICES(1)%PTR + LINEAR_MAPPING=>EQUATIONS_MAPPING%LINEAR_MAPPING + FIELD_VARIABLE=>LINEAR_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE + + STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX=0.0_DP + CASE(EQUATIONS_SET_TRANSIENT_STOKES_SUBTYPE) + DYNAMIC_MATRICES=>EQUATIONS_MATRICES%DYNAMIC_MATRICES + STIFFNESS_MATRIX=>DYNAMIC_MATRICES%MATRICES(1)%PTR + DAMPING_MATRIX=>DYNAMIC_MATRICES%MATRICES(2)%PTR + DYNAMIC_MAPPING=>EQUATIONS_MAPPING%DYNAMIC_MAPPING + FIELD_VARIABLE=>DYNAMIC_MAPPING%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE + + STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX=0.0_DP + DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX=0.0_DP + CASE DEFAULT + LOCAL_ERROR="Equations set subtype "//TRIM(NUMBER_TO_VSTRING(EQUATIONS_SET%SUBTYPE,"*",ERR,ERROR))// & + & " is not valid for a Stokes fluid type of a fluid mechanics equations set class." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + END SELECT + CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,EQUATIONS%INTERPOLATION% & + & GEOMETRIC_INTERP_PARAMETERS,ERR,ERROR,*999) + CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,ELEMENT_NUMBER,EQUATIONS%INTERPOLATION% & + & MATERIALS_INTERP_PARAMETERS,ERR,ERROR,*999) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! LOOP OVER GAUSS POINTS !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - DO ng=1,QUADRATURE_SCHEME%NUMBER_OF_GAUSS - CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,EQUATIONS%INTERPOLATION% & + DO ng=1,QUADRATURE_SCHEME%NUMBER_OF_GAUSS + CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,EQUATIONS%INTERPOLATION% & & GEOMETRIC_INTERP_POINT,ERR,ERROR,*999) - CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(GEOMETRIC_BASIS%NUMBER_OF_XI,EQUATIONS%INTERPOLATION% & + CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(GEOMETRIC_BASIS%NUMBER_OF_XI,EQUATIONS%INTERPOLATION% & & GEOMETRIC_INTERP_POINT_METRICS,ERR,ERROR,*999) - CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,EQUATIONS%INTERPOLATION% & + CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,EQUATIONS%INTERPOLATION% & & MATERIALS_INTERP_POINT,ERR,ERROR,*999) - !Define MU_PARAM, viscosity=1 - MU_PARAM=EQUATIONS%INTERPOLATION%MATERIALS_INTERP_POINT%VALUES(1,NO_PART_DERIV) + !Define MU_PARAM, viscosity=1 + MU_PARAM=EQUATIONS%INTERPOLATION%MATERIALS_INTERP_POINT%VALUES(1,NO_PART_DERIV) - !Define RHO_PARAM, density=2 - RHO_PARAM=EQUATIONS%INTERPOLATION%MATERIALS_INTERP_POINT%VALUES(2,NO_PART_DERIV) + !Define RHO_PARAM, density=2 + RHO_PARAM=EQUATIONS%INTERPOLATION%MATERIALS_INTERP_POINT%VALUES(2,NO_PART_DERIV) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! CALCULATE PARTIAL MATRICES !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! CALCULATE PARTIAL MATRICES !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - IF(EQUATIONS_SET%SUBTYPE==EQUATIONS_SET_STATIC_STOKES_SUBTYPE.OR. & + IF(EQUATIONS_SET%SUBTYPE==EQUATIONS_SET_STATIC_STOKES_SUBTYPE.OR. & & EQUATIONS_SET%SUBTYPE==EQUATIONS_SET_LAPLACE_STOKES_SUBTYPE.OR. & & EQUATIONS_SET%SUBTYPE==EQUATIONS_SET_TRANSIENT_STOKES_SUBTYPE) THEN !Loop over field components mhs=0 DO mh=1,(FIELD_VARIABLE%NUMBER_OF_COMPONENTS-1) - MESH_COMPONENT1=FIELD_VARIABLE%COMPONENTS(mh)%MESH_COMPONENT_NUMBER - DEPENDENT_BASIS1=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT1)%PTR% & + MESH_COMPONENT1=FIELD_VARIABLE%COMPONENTS(mh)%MESH_COMPONENT_NUMBER + DEPENDENT_BASIS1=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT1)%PTR% & & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS - QUADRATURE_SCHEME1=>DEPENDENT_BASIS1%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR - JGW=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN*QUADRATURE_SCHEME1%GAUSS_WEIGHTS(ng) - - DO ms=1,DEPENDENT_BASIS1%NUMBER_OF_ELEMENT_PARAMETERS - mhs=mhs+1 - nhs=0 - IF(STIFFNESS_MATRIX%UPDATE_MATRIX.OR.DAMPING_MATRIX%UPDATE_MATRIX) THEN - !Loop over element columns - DO nh=1,(FIELD_VARIABLE%NUMBER_OF_COMPONENTS) - MESH_COMPONENT2=FIELD_VARIABLE%COMPONENTS(nh)%MESH_COMPONENT_NUMBER - DEPENDENT_BASIS2=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT2)%PTR% & - & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS - QUADRATURE_SCHEME2=>DEPENDENT_BASIS2%QUADRATURE%QUADRATURE_SCHEME_MAP& - &(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR - ! JGW=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN*QUADRATURE_SCHEME2%& - ! &GAUSS_WEIGHTS(ng) - - DO ns=1,DEPENDENT_BASIS2%NUMBER_OF_ELEMENT_PARAMETERS - nhs=nhs+1 - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! PRECONDITIONS !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - - DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI - DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI - DXI_DX(mi,ni)=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%DXI_DX(mi,ni) - END DO - - DPHIMS_DXI(ni)=QUADRATURE_SCHEME1%GAUSS_BASIS_FNS(ms,PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(ni),ng) - DPHINS_DXI(ni)=QUADRATURE_SCHEME2%GAUSS_BASIS_FNS(ns,PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(ni),ng) - END DO !ni - - PHIMS=QUADRATURE_SCHEME1%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng) - PHINS=QUADRATURE_SCHEME2%GAUSS_BASIS_FNS(ns,NO_PART_DERIV,ng) - - - -! DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI -! DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI -! SUM=SUM-MU_PARAM*DPHIMSS_DXI(mi)*DPHINSS_DXI(ni)*EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%GU(mi,ni) -! ENDDO !ni -! ENDDO !mi - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! A LAPLACE ONLY !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - IF(STIFFNESS_MATRIX%UPDATE_MATRIX) THEN - !LAPLACE TYPE - IF(nh==mh) THEN - SUM=0.0_DP - !Calculate SUM - DO x=1,DEPENDENT_BASIS1%NUMBER_OF_XI - DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI - DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI - SUM=SUM+MU_PARAM*DPHINS_DXI(ni)*DXI_DX(ni,x)*DPHIMS_DXI(mi)*DXI_DX(mi,x) - ENDDO !ni - ENDDO !mi - ENDDO !x - !Calculate MATRIX - AL_MATRIX(mhs,nhs)=AL_MATRIX(mhs,nhs)+SUM*JGW - END IF - END IF - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! A STANDARD !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - - IF(STIFFNESS_MATRIX%UPDATE_MATRIX) THEN - !GRADIENT TRANSPOSE TYPE - IF(EQUATIONS_SET%SUBTYPE/=EQUATIONS_SET_LAPLACE_STOKES_SUBTYPE) THEN - IF(nhDEPENDENT_BASIS1%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR + JGW=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN*QUADRATURE_SCHEME1%GAUSS_WEIGHTS(ng) + + DO ms=1,DEPENDENT_BASIS1%NUMBER_OF_ELEMENT_PARAMETERS + mhs=mhs+1 + nhs=0 + IF(STIFFNESS_MATRIX%UPDATE_MATRIX.OR.DAMPING_MATRIX%UPDATE_MATRIX) THEN + !Loop over element columns + DO nh=1,(FIELD_VARIABLE%NUMBER_OF_COMPONENTS) + MESH_COMPONENT2=FIELD_VARIABLE%COMPONENTS(nh)%MESH_COMPONENT_NUMBER + DEPENDENT_BASIS2=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT2)%PTR% & + & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS + QUADRATURE_SCHEME2=>DEPENDENT_BASIS2%QUADRATURE%QUADRATURE_SCHEME_MAP& + &(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR + ! JGW=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%JACOBIAN*QUADRATURE_SCHEME2%& + ! &GAUSS_WEIGHTS(ng) + + DO ns=1,DEPENDENT_BASIS2%NUMBER_OF_ELEMENT_PARAMETERS + nhs=nhs+1 + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! PRECONDITIONS !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI + DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI + DXI_DX(mi,ni)=EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%DXI_DX(mi,ni) + END DO + + DPHIMS_DXI(ni)=QUADRATURE_SCHEME1%GAUSS_BASIS_FNS(ms,PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(ni),ng) + DPHINS_DXI(ni)=QUADRATURE_SCHEME2%GAUSS_BASIS_FNS(ns,PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(ni),ng) + END DO !ni + + PHIMS=QUADRATURE_SCHEME1%GAUSS_BASIS_FNS(ms,NO_PART_DERIV,ng) + PHINS=QUADRATURE_SCHEME2%GAUSS_BASIS_FNS(ns,NO_PART_DERIV,ng) + + - ENDDO !ng + ! DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI + ! DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI + ! SUM=SUM-MU_PARAM*DPHIMSS_DXI(mi)*DPHINSS_DXI(ni)*EQUATIONS%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS%GU(mi,ni) + ! ENDDO !ni + ! ENDDO !mi + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! A LAPLACE ONLY !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + IF(STIFFNESS_MATRIX%UPDATE_MATRIX) THEN + !LAPLACE TYPE + IF(nh==mh) THEN + SUM=0.0_DP + !Calculate SUM + DO x=1,DEPENDENT_BASIS1%NUMBER_OF_XI + DO mi=1,DEPENDENT_BASIS1%NUMBER_OF_XI + DO ni=1,DEPENDENT_BASIS2%NUMBER_OF_XI + SUM=SUM+MU_PARAM*DPHINS_DXI(ni)*DXI_DX(ni,x)*DPHIMS_DXI(mi)*DXI_DX(mi,x) + ENDDO !ni + ENDDO !mi + ENDDO !x + !Calculate MATRIX + AL_MATRIX(mhs,nhs)=AL_MATRIX(mhs,nhs)+SUM*JGW + END IF + END IF + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! A STANDARD !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + IF(STIFFNESS_MATRIX%UPDATE_MATRIX) THEN + !GRADIENT TRANSPOSE TYPE + IF(EQUATIONS_SET%SUBTYPE/=EQUATIONS_SET_LAPLACE_STOKES_SUBTYPE) THEN + IF(nhDEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT1)%PTR% & - & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS + & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS DO ms=1,DEPENDENT_BASIS1%NUMBER_OF_ELEMENT_PARAMETERS mhs=mhs+1 nhs=0 IF(STIFFNESS_MATRIX%UPDATE_MATRIX.OR.DAMPING_MATRIX%UPDATE_MATRIX) THEN !Loop over element columns DO nh=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS - MESH_COMPONENT2=FIELD_VARIABLE%COMPONENTS(nh)%MESH_COMPONENT_NUMBER - DEPENDENT_BASIS2=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT2)%PTR% & - & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS + MESH_COMPONENT2=FIELD_VARIABLE%COMPONENTS(nh)%MESH_COMPONENT_NUMBER + DEPENDENT_BASIS2=>DEPENDENT_FIELD%DECOMPOSITION%DOMAIN(MESH_COMPONENT2)%PTR% & + & TOPOLOGY%ELEMENTS%ELEMENTS(ELEMENT_NUMBER)%BASIS DO ns=1,DEPENDENT_BASIS2%NUMBER_OF_ELEMENT_PARAMETERS nhs=nhs+1 IF(STIFFNESS_MATRIX%UPDATE_MATRIX)THEN - STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)=STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)* & - & EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS%SCALE_FACTORS(ms,mh)* & - & EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS%SCALE_FACTORS(ns,nh) + STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)=STIFFNESS_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)* & + & EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS%SCALE_FACTORS(ms,mh)* & + & EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS%SCALE_FACTORS(ns,nh) END IF IF(DAMPING_MATRIX%UPDATE_MATRIX)THEN - DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)=DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)* & - & EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS%SCALE_FACTORS(ms,mh)* & - & EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS%SCALE_FACTORS(ns,nh) + DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)=DAMPING_MATRIX%ELEMENT_MATRIX%MATRIX(mhs,nhs)* & + & EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS%SCALE_FACTORS(ms,mh)* & + & EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS%SCALE_FACTORS(ns,nh) END IF ENDDO !ns ENDDO !nh diff --git a/src/solver_routines.f90 b/src/solver_routines.f90 index d8fcf2ef..1104a5dd 100644 --- a/src/solver_routines.f90 +++ b/src/solver_routines.f90 @@ -1277,7 +1277,8 @@ SUBROUTINE SOLVER_DAE_CRANK_NICHOLSON_INITIALISE(DAE_SOLVER,ERR,ERROR,*) IF(ASSOCIATED(DAE_SOLVER)) THEN IF(ASSOCIATED(DAE_SOLVER%CRANK_NICHOLSON_SOLVER)) THEN - CALL FLAG_ERROR("Crank-Nicholson solver is already associated for this differential-algebraic equation solver.",ERR,ERROR,*998) + CALL FLAG_ERROR("Crank-Nicholson solver is already associated for this differential-algebraic equation solver.", & + & ERR,ERROR,*998) ELSE !Allocate the Runge-Kutta solver ALLOCATE(DAE_SOLVER%CRANK_NICHOLSON_SOLVER,STAT=ERR) @@ -6501,6 +6502,7 @@ SUBROUTINE SOLVER_MATRICES_STATIC_ASSEMBLE(SOLVER,SELECTION_TYPE,ERR,ERROR,*) RHS_VECTOR=>EQUATIONS_MATRICES%RHS_VECTOR IF(ASSOCIATED(RHS_VECTOR)) THEN LINEAR_MAPPING=>EQUATIONS_MAPPING%LINEAR_MAPPING + NONLINEAR_MAPPING=>EQUATIONS_MAPPING%NONLINEAR_MAPPING IF(ASSOCIATED(LINEAR_MAPPING)) THEN LINEAR_MATRICES=>EQUATIONS_MATRICES%LINEAR_MATRICES IF(ASSOCIATED(LINEAR_MATRICES)) THEN @@ -6570,7 +6572,7 @@ SUBROUTINE SOLVER_MATRICES_STATIC_ASSEMBLE(SOLVER,SELECTION_TYPE,ERR,ERROR,*) & ERR,ERROR,*999) ENDDO !solver_row_idx !Set Dirichlet boundary conditions - IF(ASSOCIATED(LINEAR_MAPPING)) THEN + IF(ASSOCIATED(LINEAR_MAPPING).AND..NOT.ASSOCIATED(NONLINEAR_MAPPING)) THEN !Loop over the dependent variables associated with this equations set row DO variable_idx=1,LINEAR_MAPPING%NUMBER_OF_LINEAR_MATRIX_VARIABLES variable_type=LINEAR_MAPPING%LINEAR_MATRIX_VARIABLE_TYPES(variable_idx)