symbolic notation in numerics.yaml

- p_i: integration order
- p_s: shape function order

ensure working combination (p_s = p_i: full integration, p_s = p_i+1:
reduced integration)

PRIVATE

@@ -1 +1 @@
-Subproject commit fabe69749425e8a7aceb3b7c2758b40d97d8b809
+Subproject commit 0ff1d7540fede9611d46d2885bebbbe60dcbbfb0

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)

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)
 

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))

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