diff --git a/PRIVATE b/PRIVATE index cea3296f7..23a19f513 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit cea3296f78f5b48ab949364a80b738d3a7cf2194 +Subproject commit 23a19f513e760e0b6d93ebb3bbf03ca3171c71d6 diff --git a/examples/config/numerics.yaml b/examples/config/numerics.yaml index 487b6cf88..cdc487589 100644 --- a/examples/config/numerics.yaml +++ b/examples/config/numerics.yaml @@ -1,5 +1,58 @@ -# Available numerical parameters -# Case sensitive keys +# Default values of all available numerical parameters +# Please note that keys are case sensitive + +solver: + grid: + N_staggered_iter_max: 10 # max number of field-level staggered iterations + N_cutback_max: 3 # maximum cutback level (0: 1, 1: 0.5, 2: 0.25, etc) + + damage: + N_iter_max: 100 # maximum iteration number + eps_abs_phi: 1.0e-2 # absolute tolerance for damage evolution + eps_rel_phi: 1.0e-6 # relative tolerance for damage evolution + phi_min: 1.0e-6 # residual integrity + + thermal: + N_iter_max: 100 # maximum iteration number + eps_abs_T: 1.0e-2 # absolute tolerance for thermal equilibrium + eps_rel_T: 1.0e-6 # relative tolerance for thermal equilibrium + + mechanical: + N_iter_min: 1 # minimum iteration number + N_iter_max: 100 # maximum iteration number + eps_abs_div(P): 1.0e-4 # absolute tolerance for fulfillment of stress equilibrium + eps_rel_div(P): 5.0e-4 # relative tolerance for fulfillment of stress equilibrium + eps_abs_P: 1.0e3 # absolute tolerance for fulfillment of stress BC + eps_rel_P: 1.0e-3 # relative tolerance for fulfillment of stress BC + update_gamma: false # update Gamma-operator with current dPdF (not possible if FFT: memory_efficient == true) + + FFT: + memory_efficient: true # precalculate Gamma-operator (81 doubles per point) + divergence_correction: size+grid # use size-independent divergence criterion {none, size, size+grid} + derivative: continuous # approximation used for derivatives in Fourier space {continuous, central_difference, FWBW_difference} + FFTW_plan_mode: FFTW_MEASURE # planning-rigor flags, see manual at https://www.fftw.org/fftw3_doc/Planner-Flags.html + FFTW_timelimit: -1.0 # time limit of plan creation for FFTW, see manual on www.fftw.org. (-1.0: disable time limit) + PETSc_options: -snes_type ngmres -snes_ngmres_anderson # PETSc solver options + alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0 (1.0: AL scheme, 2.0: accelerated scheme) + beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0 (1.0: AL scheme, 2.0: accelerated scheme) + eps_abs_curl(F): 1.0e-10 # absolute tolerance for fulfillment of strain compatibility + eps_rel_curl(F): 5.0e-4 # relative tolerance for fulfillment of strain compatibility + + mesh: + N_cutback_max: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc) + N_staggered_iter_max: 10 # max number of field level staggered iterations + p_s: 2 # order of displacement shape functions + p_i: 2 # order of quadrature rule required + bbarstabilization: false + + mechanical: + N_iter_max: 250 # Maximum iteration number + eps_abs_div(P): 1.0e-10 # absolute tolerance for mechanical equilibrium + eps_rel_div(P): 1.0e-4 # relative tolerance for mechanical equilibrium + + Marc: + unit_length: 1.0 # physical length of one computational length unit + homogenization: mechanical: @@ -11,7 +64,7 @@ homogenization: Delta_a: 1.0e-7 # perturbation for computing penalty tangent relevant_mismatch: 1.0e-5 # minimum threshold of mismatch viscosity_exponent: 1.0e+0 # power (sensitivity rate) of numerical viscosity in RGC scheme - viscosity_modulus: 0.0e+0 # stress modulus of RGC numerical viscosity (zero = without numerical viscosity) + viscosity_modulus: 0.0e+0 # stress modulus of RGC numerical viscosity (0: without numerical viscosity) # suggestion: larger than the aTol_RGC but still far below the expected flow stress of material dot_a_ref: 1.0e-3 # reference rate of relaxation (about the same magnitude as straining rate, possibly a bit higher) dot_a_max: 1.0e+0 # threshold of maximum relaxation vector increment (if exceed this then cutback) @@ -19,77 +72,29 @@ homogenization: Delta_V_modulus: 1.0e+12 Delta_V_exponent: 5.0 -solver: - grid: - N_staggered_iter_max: 10 # max number of field level staggered iterations - N_cutback_max: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc) - - damage: - N_iter_max: 100 # maximum iteration number - eps_abs_phi: 1.0e-2 # absolute tolerance for damage evolution - eps_rel_phi: 1.0e-6 # relative tolerance for damage evolution - thermal: - N_iter_max: 100 # maximum iteration number - eps_abs_T: 1.0e-2 # absolute tolerance for thermal equilibrium - eps_rel_T: 1.0e-6 # relative tolerance for thermal equilibrium - - mechanical: - eps_abs_div(P): 1.0e-4 # absolute tolerance for fulfillment of stress equilibrium - eps_rel_div(P): 5.0e-4 # relative tolerance for fulfillment of stress equilibrium - eps_abs_P: 1.0e3 # absolute tolerance for fulfillment of stress BC - eps_rel_P: 1.0e-3 # relative tolerance for fulfillment of stress BC - N_iter_min: 1 # minimum iteration number - N_iter_max: 100 # maximum iteration number - update_gamma: false # Update Gamma-operator with current dPdF (not possible if memory_efficient=1) - - FFT: - memory_efficient: True # Precalculate Gamma-operator (81 double per point) - divergence_correction: size+grid # Use size-independent divergence criterion - derivative: continuous # Approximation used for derivatives in Fourier space - FFTW_plan_mode: FFTW_MEASURE # planing-rigor flag, see manual on www.fftw.org - FFTW_timelimit: -1.0 # timelimit of plan creation for FFTW, see manual on www.fftw.org. -1.0: disable timelimit - PETSc_options: -snes_type ngmres -snes_ngmres_anderson # PETSc solver options - alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme - beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme - eps_abs_curl(F): 1.0e-10 # absolute tolerance for fulfillment of strain compatibility - eps_rel_curl(F): 5.0e-4 # relative tolerance for fulfillment of strain compatibility - -mesh: - maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc) - maxStaggeredIter: 10 # max number of field level staggered iterations - structorder: 2 # order of displacement shape functions (when mesh is defined) - bbarstabilisation: false - integrationorder: 2 # order of quadrature rule required (when mesh is defined) - itmax: 250 # Maximum iteration number - itmin: 2 # Minimum iteration number - eps_struct_atol: 1.0e-10 # absolute tolerance for mechanical equilibrium - eps_struct_rtol: 1.0e-4 # relative tolerance for mechanical equilibrium phase: mechanical: - r_cutback_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in phase state calculation - r_cutback: 0.25 # factor to decrease size of step when cutback introduced in phase state calculation (value between 0 and 1) - r_increase: 1.5 # factor to increase size of next step when previous step converged in phase state calculation - eps_rel_state: 1.0e-6 # relative tolerance in phase state loop (abs tol provided by constitutive law) - N_iter_state_max: 10 # state loop limit + r_cutback_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in phase state calculation + r_cutback: 0.25 # factor to decrease size of step when cutback introduced in phase state calculation (value between 0 and 1) + r_increase: 1.5 # factor to increase size of next step when previous step converged in phase state calculation + eps_rel_state: 1.0e-6 # relative tolerance in phase state loop (abs tol provided by constitutive law) + N_iter_state_max: 10 # state loop limit plastic: - r_linesearch_Lp: 0.5 # factor to decrease the step due to non-convergence in Lp calculation - eps_rel_Lp: 1.0e-6 # relative tolerance in Lp residuum - eps_abs_Lp: 1.0e-8 # absolute tolerance in Lp residuum - N_iter_Lp_max: 40 # stress loop limit for Lp - f_update_jacobi_Lp: 1 # frequency of Jacobian update of residuum in Lp - integrator_state: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) - eigen: - r_linesearch_Li: 0.5 # factor to decrease the step due to non-convergence in Li calculation - eps_rel_Li: 1.0e-6 # relative tolerance in Li residuum - eps_abs_Li: 1.0e-8 # absolute tolerance in Li residuum - N_iter_Li_max: 40 # stress loop limit for Li - f_update_jacobi_Li: 1 # frequency of Jacobian update of residuum in Li + r_linesearch_Lp: 0.5 # factor to decrease the step if Lp calculation fails to converge + eps_rel_Lp: 1.0e-6 # relative tolerance in Lp residuum + eps_abs_Lp: 1.0e-8 # absolute tolerance in Lp residuum + N_iter_Lp_max: 40 # stress loop limit for Lp + f_update_jacobi_Lp: 1 # frequency of Jacobian update of residuum in Lp + integrator_state: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) -commercialFEM: - unitlength: 1 # physical length of one computational length unit + eigen: + r_linesearch_Li: 0.5 # factor to decrease the step if Li calculation fails to converge + eps_rel_Li: 1.0e-6 # relative tolerance in Li residuum + eps_abs_Li: 1.0e-8 # absolute tolerance in Li residuum + N_iter_Li_max: 40 # stress loop limit for Li + f_update_jacobi_Li: 1 # frequency of updating the Jacobian of residuum in Li generic: - random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed. - phi_min: 1.0e-6 # non-zero residual damage. + random_seed: 0 # fixed seeding for pseudo-random number generator (0: use random seed) diff --git a/src/Marc/discretization_Marc.f90 b/src/Marc/discretization_Marc.f90 index 3d9778e0e..5e42180de 100644 --- a/src/Marc/discretization_Marc.f90 +++ b/src/Marc/discretization_Marc.f90 @@ -69,14 +69,16 @@ subroutine discretization_Marc_init unscaledNormals type(tDict), pointer :: & + num_solver, & num_commercialFEM print'(/,a)', ' <<<+- discretization_Marc init -+>>>'; flush(6) - num_commercialFEM => config_numerics%get_dict('commercialFEM',defaultVal = emptyDict) - mesh_unitlength = num_commercialFEM%get_asReal('unitlength',defaultVal=1.0_pREAL) ! set physical extent of a length unit in mesh - if (mesh_unitlength <= 0.0_pREAL) call IO_error(301,'unitlength') + num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict) + num_commercialFEM => num_solver%get_dict('Marc',defaultVal=emptyDict) + mesh_unitlength = num_commercialFEM%get_asReal('unit_length',defaultVal=1.0_pREAL) ! set physical extent of a length unit in mesh + if (mesh_unitlength <= 0.0_pREAL) call IO_error(301,'unit_length') call inputRead(elem,node0_elem,connectivity_elem,materialAt) nElems = size(connectivity_elem,2) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 48023cc54..648f9ebbb 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -102,7 +102,7 @@ subroutine grid_damage_spectral_init(num_grid) ! read numerical parameters and do sanity checks num_grid_damage => num_grid%get_dict('damage',defaultVal=emptyDict) - num%itmax = num_grid_damage%get_asInt('N_iter_max', defaultVal=100) + num%itmax = num_grid_damage%get_asInt ('N_iter_max', defaultVal=100) num%eps_damage_atol = num_grid_damage%get_asReal('eps_abs_phi',defaultVal=1.0e-2_pREAL) num%eps_damage_rtol = num_grid_damage%get_asReal('eps_rel_phi',defaultVal=1.0e-6_pREAL) num%phi_min = num_grid_damage%get_asReal('phi_min', defaultVal=1.0e-6_pREAL) @@ -110,7 +110,7 @@ subroutine grid_damage_spectral_init(num_grid) extmsg = '' if (num%eps_damage_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_phi' if (num%eps_damage_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_phi' - if (num%phi_min < 0.0_pREAL) extmsg = trim(extmsg)//' phi_min' + if (num%phi_min <= 0.0_pREAL) extmsg = trim(extmsg)//' phi_min' if (num%itmax < 1) extmsg = trim(extmsg)//' N_iter_max' if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg)) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 4bfd791f5..370dc3535 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -66,6 +66,7 @@ program DAMASK_mesh stagIter, & component type(tDict), pointer :: & + num_solver, & num_mesh character(len=pSTRLEN), dimension(:), allocatable :: fileContent character(len=pSTRLEN) :: & @@ -90,12 +91,13 @@ program DAMASK_mesh !--------------------------------------------------------------------- ! reading field information from numerics file and do sanity checks - num_mesh => config_numerics%get_dict('mesh', defaultVal=emptyDict) - stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10) - maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3) + num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict) + num_mesh => num_solver%get_dict('mesh',defaultVal=emptyDict) + stagItMax = num_mesh%get_asInt('N_staggered_iter_max',defaultVal=10) + maxCutBack = num_mesh%get_asInt('N_cutback_max',defaultVal=3) - if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') - if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') + if (stagItMax < 0) call IO_error(301,ext_msg='N_staggered_iter_max') + if (maxCutBack < 0) call IO_error(301,ext_msg='N_cutback_max') ! reading basic information from load case file and allocate data structure containing load cases call DMGetDimension(geomMesh,dimPlex,err_PETSc) !< dimension of mesh (2D or 3D) @@ -229,8 +231,8 @@ program DAMASK_mesh !-------------------------------------------------------------------------------------------------- ! doing initialization depending on active solvers - call FEM_Utilities_init() - call FEM_mechanical_init(loadCases(1)%fieldBC(1)) + call FEM_Utilities_init(num_mesh) + call FEM_mechanical_init(loadCases(1)%fieldBC(1),num_mesh) call config_numerics_deallocate() if (worldrank == 0) then diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index d30b223f0..4480e412d 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -16,6 +16,7 @@ module FEM_utilities use prec use config use math + use misc use IO use discretization_mesh use homogenization @@ -90,11 +91,16 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief Allocate all neccessary fields. !-------------------------------------------------------------------------------------------------- -subroutine FEM_utilities_init +subroutine FEM_utilities_init(num_mesh) - character(len=pSTRLEN) :: petsc_optionsOrder - type(tDict), pointer :: & + type(tDict), pointer, intent(in) :: & num_mesh + + type(tDict), pointer :: & + num_mech + character(len=pSTRLEN) :: petsc_optionsOrder + character(len=:), allocatable :: & + petsc_options integer :: & p_s, & !< order of shape functions p_i !< integration order (quadrature rule) @@ -103,7 +109,7 @@ subroutine FEM_utilities_init print'(/,1x,a)', '<<<+- FEM_utilities init -+>>>' - num_mesh => config_numerics%get_dict('mesh',defaultVal=emptyDict) + num_mech => num_mesh%get_dict('mechanical', defaultVal=emptyDict) p_s = num_mesh%get_asInt('p_s',defaultVal = 2) p_i = num_mesh%get_asInt('p_i',defaultVal = p_s) @@ -117,20 +123,20 @@ subroutine FEM_utilities_init call PetscOptionsClear(PETSC_NULL_OPTIONS,err_PETSc) CHKERRQ(err_PETSc) CHKERRQ(err_PETSc) - call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type newtonls & - &-mechanical_snes_linesearch_type cp -mechanical_snes_ksp_ew & - &-mechanical_snes_ksp_ew_rtol0 0.01 -mechanical_snes_ksp_ew_rtolmax 0.01 & - &-mechanical_ksp_type fgmres -mechanical_ksp_max_it 25', err_PETSc) - CHKERRQ(err_PETSc) - call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asStr('PETSc_options',defaultVal=''),err_PETSc) - CHKERRQ(err_PETSc) + + petsc_options = misc_prefixOptions('-snes_type newtonls & + &-snes_linesearch_type cp -snes_ksp_ew & + &-snes_ksp_ew_rtol0 0.01 -snes_ksp_ew_rtolmax 0.01 & + &-ksp_type fgmres -ksp_max_it 25 ' // & + num_mech%get_asStr('PETSc_options',defaultVal=''), 'mechanical_') + write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', p_s - call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),err_PETSc) + petsc_options = petsc_options // ' ' // petsc_optionsOrder + call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,err_PETSc) CHKERRQ(err_PETSc) wgt = real(mesh_maxNips*mesh_NcpElemsGlobal,pREAL)**(-1) - end subroutine FEM_utilities_init diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 3f5811a5c..be9be3b19 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -89,6 +89,7 @@ subroutine discretization_mesh_init(restart) PetscInt, dimension(:), allocatable :: & materialAt type(tDict), pointer :: & + num_solver, & num_mesh integer :: p_i, dim !< integration order (quadrature rule) type(tvec) :: coords_node0 @@ -99,8 +100,9 @@ subroutine discretization_mesh_init(restart) !-------------------------------------------------------------------------------- ! read numerics parameter - num_mesh => config_numerics%get_dict('mesh',defaultVal=emptyDict) - p_i = num_mesh%get_asInt('p_i',defaultVal = 2) + num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict) + num_mesh => num_solver%get_dict('mesh',defaultVal=emptyDict) + p_i = num_mesh%get_asInt('p_i',defaultVal=2) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>16) call DMPlexCreateFromFile(PETSC_COMM_WORLD,CLI_geomFile,'n/a',PETSC_TRUE,globalMesh,err_PETSc) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index d08906439..c566990e2 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -47,7 +47,7 @@ module mesh_mechanical_FEM p_i, & !< integration order (quadrature rule) itmax logical :: & - BBarStabilisation + BBarStabilization real(pREAL) :: & eps_struct_atol, & !< absolute tolerance for mechanical equilibrium eps_struct_rtol !< relative tolerance for mechanical equilibrium @@ -94,9 +94,10 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data !-------------------------------------------------------------------------------------------------- -subroutine FEM_mechanical_init(fieldBC) +subroutine FEM_mechanical_init(fieldBC,num_mesh) - type(tFieldBC), intent(in) :: fieldBC + type(tFieldBC), intent(in) :: fieldBC + type(tDict), pointer, intent(in) :: num_mesh DM :: mechanical_mesh PetscFE :: mechFE @@ -126,23 +127,24 @@ subroutine FEM_mechanical_init(fieldBC) character(len=*), parameter :: prefix = 'mechFE_' PetscErrorCode :: err_PETSc real(pREAL), dimension(3,3) :: devNull - type(tDict), pointer :: & - num_mesh + type(tDict), pointer :: num_mech print'(/,1x,a)', '<<<+- FEM_mech init -+>>>'; flush(IO_STDOUT) !----------------------------------------------------------------------------- ! read numerical parametes and do sanity checks - num_mesh => config_numerics%get_dict('mesh',defaultVal=emptyDict) - num%p_i = int(num_mesh%get_asInt('p_i',defaultVal = 2),pPETSCINT) - num%itmax = int(num_mesh%get_asInt('itmax',defaultVal=250),pPETSCINT) - num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.) - num%eps_struct_atol = num_mesh%get_asReal('eps_struct_atol', defaultVal = 1.0e-10_pREAL) - num%eps_struct_rtol = num_mesh%get_asReal('eps_struct_rtol', defaultVal = 1.0e-4_pREAL) + num_mech => num_mesh%get_dict('mechanical', defaultVal=emptyDict) - if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') - if (num%eps_struct_rtol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_struct_rtol') - if (num%eps_struct_atol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_struct_atol') + num%p_i = int(num_mesh%get_asInt('p_i',defaultVal=2),pPETSCINT) + num%BBarStabilization = num_mesh%get_asBool('bbarstabilization',defaultVal=.false.) + + num%itmax = int(num_mech%get_asInt('N_iter_max',defaultVal=250),pPETSCINT) + num%eps_struct_atol = num_mech%get_asReal('eps_abs_div(P)', defaultVal=1.0e-10_pREAL) + num%eps_struct_rtol = num_mech%get_asReal('eps_rel_div(P)', defaultVal=1.0e-4_pREAL) + + if (num%itmax <= 1) call IO_error(301,ext_msg='N_iter_max') + if (num%eps_struct_rtol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_rel_div(P)') + if (num%eps_struct_atol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_abs_div(P)') !-------------------------------------------------------------------------------------------------- ! Setup FEM mech mesh @@ -437,7 +439,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc end do homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) end do - if (num%BBarStabilisation) then + if (num%BBarStabilization) then detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature,pREAL)) do qPt = 0, nQuadrature-1 m = cell*nQuadrature + qPt+1 @@ -588,7 +590,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT) - if (num%BBarStabilisation) then + if (num%BBarStabilization) then F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex]) FInv = math_inv33(F) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0_pREAL/real(dimPlex,pREAL)) @@ -604,7 +606,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P K_eA = K_eA + matmul(transpose(BMat),MatA) end if end do - if (num%BBarStabilisation) then + if (num%BBarStabilization) then FInv = math_inv33(FAvg) K_e = K_eA*math_det33(FAvg/real(nQuadrature,pREAL))**(1.0_pREAL/real(dimPlex,pREAL)) + & (matmul(matmul(transpose(BMatAvg), &