diff --git a/PRIVATE b/PRIVATE index fabe69749..0ff1d7540 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit fabe69749425e8a7aceb3b7c2758b40d97d8b809 +Subproject commit 0ff1d7540fede9611d46d2885bebbbe60dcbbfb0 diff --git a/src/IO.f90 b/src/IO.f90 index cd7c09c75..717493006 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -432,7 +432,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'Nconstituents mismatch between homogenization and material' !-------------------------------------------------------------------------------------------------- -! material error messages and related messages in mesh +! material error messages and related messages in geometry case (150) msg = 'index out of bounds' case (153) @@ -499,6 +499,11 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (710) msg = 'Closing quotation mark missing in string' +!------------------------------------------------------------------------------------------------- +! errors related to the mesh solver + case (821) + msg = 'order not supported' + !------------------------------------------------------------------------------------------------- ! errors related to the grid solver case (831) diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 64bcc3896..a3856ccaa 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -19,6 +19,7 @@ module FEM_utilities use IO use discretization_mesh use homogenization + use FEM_quadrature implicit none private @@ -29,8 +30,8 @@ module FEM_utilities !-------------------------------------------------------------------------------------------------- ! field labels information - character(len=*), parameter, public :: & - FIELD_MECH_label = 'mechanical' + character(len=*), parameter, public :: & + FIELD_MECH_label = 'mechanical' enum, bind(c); enumerator :: & FIELD_UNDEFINED_ID, & @@ -86,7 +87,9 @@ subroutine FEM_utilities_init class(tNode), pointer :: & num_mesh, & debug_mesh ! pointer to mesh debug options - integer :: structOrder !< order of displacement shape functions + integer :: & + p_s, & !< order of shape functions + p_i !< integration order (quadrature rule) character(len=*), parameter :: & PETSCDEBUG = ' -snes_view -snes_monitor ' PetscErrorCode :: ierr @@ -96,7 +99,14 @@ subroutine FEM_utilities_init print'(/,a)', ' <<<+- FEM_utilities init -+>>>' num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) - structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2) + + p_s = num_mesh%get_asInt('p_s',defaultVal = 2) + p_i = num_mesh%get_asInt('p_i',defaultVal = p_s) + + if (p_s < 1_pInt .or. p_s > size(FEM_nQuadrature,2)) & + call IO_error(821,ext_msg='shape function order (p_s) out of bounds') + if (p_i < max(1_pInt,p_s-1_pInt) .or. p_i > p_s) & + call IO_error(821,ext_msg='integration order (p_i) out of bounds') debug_mesh => config_debug%get('mesh',defaultVal=emptyList) debugPETSc = debug_mesh%contains('PETSc') @@ -119,7 +129,7 @@ subroutine FEM_utilities_init CHKERRQ(ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asString('PETSc_options',defaultVal=''),ierr) CHKERRQ(ierr) - write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder + write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', p_s call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr) CHKERRQ(ierr) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index e5f989484..f1d38760d 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -85,7 +85,7 @@ subroutine discretization_mesh_init(restart) materialAt class(tNode), pointer :: & num_mesh - integer :: integrationOrder !< order of quadrature rule required + integer :: p_i !< integration order (quadrature rule) type(tvec) :: coords_node0 print'(/,a)', ' <<<+- discretization_mesh init -+>>>' @@ -93,7 +93,7 @@ subroutine discretization_mesh_init(restart) !-------------------------------------------------------------------------------- ! read numerics parameter num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) - integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2) + p_i = num_mesh%get_asInt('p_i',defaultVal = 2) !--------------------------------------------------------------------------------- ! read debug parameters @@ -150,9 +150,9 @@ subroutine discretization_mesh_init(restart) call VecGetArrayF90(coords_node0, mesh_node0_temp,ierr) CHKERRQ(ierr) - mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder) + mesh_maxNips = FEM_nQuadrature(dimPlex,p_i) - call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) + call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,p_i)%p) call mesh_FEM_build_ipVolumes(dimPlex) allocate(materialAt(mesh_NcpElems)) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index d6d314a42..58429992c 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -41,7 +41,7 @@ module mesh_mechanical_FEM type, private :: tNumerics integer :: & - integrationOrder, & !< order of quadrature rule required + p_i, & !< integration order itmax logical :: & BBarStabilisation @@ -118,7 +118,7 @@ subroutine FEM_mechanical_init(fieldBC) !----------------------------------------------------------------------------- ! read numerical parametes and do sanity checks num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) - num%integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2) + num%p_i = num_mesh%get_asInt('p_i',defaultVal = 2) num%itmax = num_mesh%get_asInt('itmax',defaultVal=250) num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.) num%eps_struct_atol = num_mesh%get_asFloat('eps_struct_atol', defaultVal = 1.0e-10_pReal) @@ -135,9 +135,9 @@ subroutine FEM_mechanical_init(fieldBC) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech discretization - qPoints = FEM_quadrature_points( dimPlex,num%integrationOrder)%p - qWeights = FEM_quadrature_weights(dimPlex,num%integrationOrder)%p - nQuadrature = FEM_nQuadrature( dimPlex,num%integrationOrder) + qPoints = FEM_quadrature_points( dimPlex,num%p_i)%p + qWeights = FEM_quadrature_weights(dimPlex,num%p_i)%p + nQuadrature = FEM_nQuadrature( dimPlex,num%p_i) qPointsP => qPoints qWeightsP => qWeights call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr) @@ -146,7 +146,7 @@ subroutine FEM_mechanical_init(fieldBC) call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr) CHKERRQ(ierr) call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, & - num%integrationOrder,mechFE,ierr); CHKERRQ(ierr) + num%p_i,mechFE,ierr); CHKERRQ(ierr) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) nBasis = nBasis/nc