diff --git a/.gitlab/issue_templates/new_release.md b/.gitlab/issue_templates/new_release.md index a97731a44..25792df72 100644 --- a/.gitlab/issue_templates/new_release.md +++ b/.gitlab/issue_templates/new_release.md @@ -13,6 +13,7 @@ Detailed instructions for each step are documented in the [wiki](https://git.dam - [ ] PyPi package (wheel) - [ ] Arch Linux (PKGBUILD) - [ ] Ubuntu (Launchpad) + - [ ] Flatpak (Flathub) - Open build service - [ ] Debian - [ ] Fedora/Suse diff --git a/LICENSE b/LICENSE index b1c966ca8..28cc732e4 100644 --- a/LICENSE +++ b/LICENSE @@ -1,3 +1,4 @@ +Copyright 2024 Max-Planck-Institut für Nachhaltige Materialien GmbH Copyright 2011-2024 Max-Planck-Institut für Eisenforschung GmbH DAMASK is free software: you can redistribute it and/or modify diff --git a/PRIVATE b/PRIVATE index 890c64d34..d228c6bba 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 890c64d34dbc38baa01d350ab7ddd02fa15d9081 +Subproject commit d228c6bbafa8ed121d27cc54063c49140a5cfd12 diff --git a/README.md b/README.md index d1aedd38f..0da1b776f 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,7 @@ The site is primarily meant to provide a forum for [Discussions](https://github. ## Contact Information -Max-Planck-Institut für Eisenforschung GmbH +Max-Planck-Institut für Nachhaltige Materialien GmbH Max-Planck-Str. 1 40237 Düsseldorf Germany diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake index 079574c4b..388115594 100644 --- a/cmake/Compiler-IntelLLVM.cmake +++ b/cmake/Compiler-IntelLLVM.cmake @@ -19,7 +19,7 @@ endif () set (STANDARD_CHECK "-stand f18 -assume nostd_mod_proc_name") set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel") # Link against shared Intel libraries instead of static ones -set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel -fc=ifx") +set (LINKER_FLAGS "${LINKER_FLAGS} -fc=ifx") # enforce use of ifx for MPI wrapper #------------------------------------------------------------------------------------------------ diff --git a/env/DAMASK.sh b/env/DAMASK.sh index e26389257..d73d94285 100644 --- a/env/DAMASK.sh +++ b/env/DAMASK.sh @@ -27,8 +27,10 @@ cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd PATH=${DAMASK_ROOT}/bin:$PATH -SOLVER=$(type -p DAMASK_grid || true 2>/dev/null) -[ "x$SOLVER" == "x" ] && SOLVER=$(blink 'Not found!') +SOLVER_GRID=$(type -p DAMASK_grid || true 2>/dev/null) +[ "x$SOLVER_GRID" == "x" ] && SOLVER_GRID=$(blink 'Not found!') +SOLVER_MESH=$(type -p DAMASK_mesh || true 2>/dev/null) +[ "x$SOLVER_MESH" == "x" ] && SOLVER_MESH=$(blink 'Not found!') # currently, there is no information that unlimited stack size causes problems @@ -42,12 +44,13 @@ ulimit -s unlimited 2>/dev/null # maximum stack size (kB) if [ ! -z "$PS1" ]; then echo echo Düsseldorf Advanced Materials Simulation Kit --- DAMASK - echo Max-Planck-Institut für Eisenforschung GmbH, Düsseldorf + echo Max-Planck-Institut für Nachhaltige Materialien GmbH, Düsseldorf echo https://damask.mpie.de echo echo Using environment with ... echo "DAMASK $DAMASK_ROOT $BRANCH" - echo "Grid Solver $SOLVER" + echo "Grid Solver $SOLVER_GRID" + echo "Mesh Solver $SOLVER_MESH" if [ "x$PETSC_DIR" != "x" ]; then echo -n "PETSc location " [ -d $PETSC_DIR ] && echo $PETSC_DIR || blink $PETSC_DIR diff --git a/env/DAMASK.zsh b/env/DAMASK.zsh index 57efa7667..68666ca0c 100644 --- a/env/DAMASK.zsh +++ b/env/DAMASK.zsh @@ -17,8 +17,10 @@ cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd PATH=${DAMASK_ROOT}/bin:$PATH -SOLVER=$(which DAMASK_grid || true 2>/dev/null) -[[ "x$SOLVER" == "x" ]] && SOLVER=$(blink 'Not found!') +SOLVER_GRID=$(which DAMASK_grid || true 2>/dev/null) +[[ "x$SOLVER_GRID" == "x" ]] && SOLVER_GRID=$(blink 'Not found!') +SOLVER_MESH=$(which DAMASK_mesh || true 2>/dev/null) +[[ "x$SOLVER_MESH" == "x" ]] && SOLVER_MESH=$(blink 'Not found!') # currently, there is no information that unlimited stack size causes problems @@ -32,12 +34,13 @@ ulimit -s unlimited 2>/dev/null # maximum stack size (kB) if [ ! -z "$PS1" ]; then echo echo Düsseldorf Advanced Materials Simulation Kit --- DAMASK - echo Max-Planck-Institut für Eisenforschung GmbH, Düsseldorf + echo Max-Planck-Institut für Nachhaltige Materialien GmbH, Düsseldorf echo https://damask.mpie.de echo echo "Using environment with ..." echo "DAMASK $DAMASK_ROOT $BRANCH" - echo "Grid Solver $SOLVER" + echo "Grid Solver $SOLVER_GRID" + echo "Mesh Solver $SOLVER_MESH" if [ "x$PETSC_DIR" != "x" ]; then echo -n "PETSc location " [ -d $PETSC_DIR ] && echo $PETSC_DIR || blink $PETSC_DIR diff --git a/img/DAMASK_Logo.svg b/img/DAMASK_Logo.svg index 87e3a219b..0df1bcfea 100644 --- a/img/DAMASK_Logo.svg +++ b/img/DAMASK_Logo.svg @@ -24,7 +24,7 @@ id="metadata32">image/svg+xmlDAMASK logoPhilip EisenlohrDAMASK; Crystal Plasticity; Multi-PhysicsMax-Planck-Institut für Eisenforschung GmbHDAMASK + rdf:resource="http://creativecommons.org/licenses/by-nd/4.0/" />DAMASK; Crystal Plasticity; Multi-PhysicsMax-Planck-Institut für Nachhaltige Materialien GmbHDAMASK The Düsseldorf Advanced Material Simulation Kit Rotation: """ Return symmetry operations. + Notes + ----- + The symmetry operations defined here only consider Rotations. + More specifically, for each crystal family, an enantiomorphic + point symmetry is selected. In case that there are multiple + point groups with enantiomorphic point symmetry, the one with + the highest order is chosen: + + Overview of crystal classes and point group in Hermann-Mauguin + notation used for definition of symmetry operations. + - tricinic: 1 + - monoclinic: 2 + - orthorhombic: 222 + - tetragonal: 422 + - hexagonal: 622 + - cubic: 432 + + References ---------- U.F. Kocks et al., - Texture and Anisotropy: - Preferred Orientations in Polycrystals and their Effect on Materials Properties. + Texture and Anisotropy: Preferred Orientations in Polycrystals + and their Effect on Materials Properties. Cambridge University Press 1998. Table II + https://en.wikipedia.org/wiki/Crystal_system#Crystal_classes + """ _symmetry_operations: Dict[CrystalFamily, List] = { 'cubic': [ @@ -719,7 +739,7 @@ def symmetry_operations(self) -> Rotation: [-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ], [-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ], [-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ], - ], + ], # 432 'hexagonal': [ [ 1.0, 0.0, 0.0, 0.0 ], [-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ], @@ -733,7 +753,7 @@ def symmetry_operations(self) -> Rotation: [ 0.0, 0.0, 1.0, 0.0 ], [ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ], [ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ], - ], + ], # 622 'tetragonal': [ [ 1.0, 0.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0, 0.0 ], @@ -743,20 +763,20 @@ def symmetry_operations(self) -> Rotation: [ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ], [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], - ], + ], # 422 'orthorhombic': [ [ 1.0,0.0,0.0,0.0 ], [ 0.0,1.0,0.0,0.0 ], [ 0.0,0.0,1.0,0.0 ], [ 0.0,0.0,0.0,1.0 ], - ], + ], # 222 'monoclinic': [ [ 1.0,0.0,0.0,0.0 ], [ 0.0,0.0,1.0,0.0 ], - ], + ], # 2 'triclinic': [ [ 1.0,0.0,0.0,0.0 ], - ]} + ]} # 1 return Rotation.from_quaternion(_symmetry_operations[self.family],accept_homomorph=True) diff --git a/python/damask/_geomgrid.py b/python/damask/_geomgrid.py index 9d03da505..898f4fce4 100644 --- a/python/damask/_geomgrid.py +++ b/python/damask/_geomgrid.py @@ -476,6 +476,10 @@ def from_Laguerre_tessellation(cells: IntSequence, periodic : bool, optional Assume grid to be periodic. Defaults to True. + Notes + ----- + damask.seeds contains functionality for seed generation. + Returns ------- new : damask.GeomGrid @@ -537,6 +541,27 @@ def from_Voronoi_tessellation(cells: IntSequence, new : damask.GeomGrid Grid-based geometry from tessellation. + Notes + ----- + damask.seeds contains functionality for seed generation. + + Examples + -------- + Generate microstructure with three grains. + + >>> import numpy as np + >>> import damask + >>> seeds = np.array([[0.30828, 0.64020, 0.65237], + ... [0.62200, 0.56858, 0.32842], + ... [0.57315, 0.94534, 0.87531]])*1e-6 + >>> damask.GeomGrid.from_Voronoi_tessellation(cells=[10,10,10], + ... size=[1e-6]*3, + ... seeds=seeds) + cells: 10 × 10 × 10 + size: 1e-06 × 1e-06 × 1e-06 m³ + origin: 0.0 0.0 0.0 m + # materials: 3 + """ coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3) tree = spatial.cKDTree(seeds,boxsize=size) if periodic else \ @@ -882,7 +907,7 @@ def flip(self, Invariance of flipping a (fully) mirrored grid. - >>> g.mirror('x',True) == g.mirror('x',True).flip('x') + >>> g.mirror('x',reflect=True) == g.mirror('x',reflect=True).flip('x') True """ @@ -1183,6 +1208,14 @@ def add_primitive(self, """ Insert a primitive geometric object at a given position. + The shape of the object is defined as + + .. math:: + + |x|^{2^a} + |y|^{2^b} + |z|^{2^c} < 1 + + where dimension = (x,y,z) and exponent = (a,b,c). + Parameters ---------- dimension : sequence of int or float, len (3) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 653cf95c0..96d3cb782 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -97,15 +97,15 @@ program DAMASK_grid type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tSolutionState), allocatable, dimension(:) :: solres procedure(grid_mechanical_spectral_basic_init), pointer :: & - mechanical_init + grid_mechanical_init procedure(grid_mechanical_spectral_basic_forward), pointer :: & - mechanical_forward + grid_mechanical_forward procedure(grid_mechanical_spectral_basic_solution), pointer :: & - mechanical_solution + grid_mechanical_solution procedure(grid_mechanical_spectral_basic_updateCoords), pointer :: & - mechanical_updateCoords + grid_mechanical_updateCoords procedure(grid_mechanical_spectral_basic_restartWrite), pointer :: & - mechanical_restartWrite + grid_mechanical_restartWrite external :: & quit @@ -160,35 +160,35 @@ program DAMASK_grid nActiveFields = 1 select case (solver%get_asStr('mechanical')) case ('spectral_basic') - mechanical_init => grid_mechanical_spectral_basic_init - mechanical_forward => grid_mechanical_spectral_basic_forward - mechanical_solution => grid_mechanical_spectral_basic_solution - mechanical_updateCoords => grid_mechanical_spectral_basic_updateCoords - mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite + grid_mechanical_init => grid_mechanical_spectral_basic_init + grid_mechanical_forward => grid_mechanical_spectral_basic_forward + grid_mechanical_solution => grid_mechanical_spectral_basic_solution + grid_mechanical_updateCoords => grid_mechanical_spectral_basic_updateCoords + grid_mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite active_Gamma = .true. case ('spectral_polarization') - mechanical_init => grid_mechanical_spectral_polarization_init - mechanical_forward => grid_mechanical_spectral_polarization_forward - mechanical_solution => grid_mechanical_spectral_polarization_solution - mechanical_updateCoords => grid_mechanical_spectral_polarization_updateCoords - mechanical_restartWrite => grid_mechanical_spectral_polarization_restartWrite + grid_mechanical_init => grid_mechanical_spectral_polarization_init + grid_mechanical_forward => grid_mechanical_spectral_polarization_forward + grid_mechanical_solution => grid_mechanical_spectral_polarization_solution + grid_mechanical_updateCoords => grid_mechanical_spectral_polarization_updateCoords + grid_mechanical_restartWrite => grid_mechanical_spectral_polarization_restartWrite active_Gamma = .true. case ('spectral_Galerkin') - mechanical_init => grid_mechanical_spectral_Galerkin_init - mechanical_forward => grid_mechanical_spectral_Galerkin_forward - mechanical_solution => grid_mechanical_spectral_Galerkin_solution - mechanical_updateCoords => grid_mechanical_spectral_Galerkin_updateCoords - mechanical_restartWrite => grid_mechanical_spectral_Galerkin_restartWrite + grid_mechanical_init => grid_mechanical_spectral_Galerkin_init + grid_mechanical_forward => grid_mechanical_spectral_Galerkin_forward + grid_mechanical_solution => grid_mechanical_spectral_Galerkin_solution + grid_mechanical_updateCoords => grid_mechanical_spectral_Galerkin_updateCoords + grid_mechanical_restartWrite => grid_mechanical_spectral_Galerkin_restartWrite active_G = .true. case ('FEM') - mechanical_init => grid_mechanical_FEM_init - mechanical_forward => grid_mechanical_FEM_forward - mechanical_solution => grid_mechanical_FEM_solution - mechanical_updateCoords => grid_mechanical_FEM_updateCoords - mechanical_restartWrite => grid_mechanical_FEM_restartWrite + grid_mechanical_init => grid_mechanical_FEM_init + grid_mechanical_forward => grid_mechanical_FEM_forward + grid_mechanical_solution => grid_mechanical_FEM_solution + grid_mechanical_updateCoords => grid_mechanical_FEM_updateCoords + grid_mechanical_restartWrite => grid_mechanical_FEM_restartWrite case default call IO_error(error_ID = 891, ext_msg = trim(solver%get_asStr('mechanical'))) @@ -232,7 +232,7 @@ program DAMASK_grid end select end do - call mechanical_init(num_grid) + call grid_mechanical_init(num_grid) call config_numerics_deallocate() !-------------------------------------------------------------------------------------------------- @@ -302,7 +302,7 @@ program DAMASK_grid do field = 1, nActiveFields select case(ID(field)) case(FIELD_MECH_ID) - call mechanical_forward (& + call grid_mechanical_forward (& cutBack,guess,Delta_t,Delta_t_prev,t_remaining, & deformation_BC = loadCases(l)%deformation, & stress_BC = loadCases(l)%stress, & @@ -324,7 +324,7 @@ program DAMASK_grid do field = 1, nActiveFields select case(ID(field)) case(FIELD_MECH_ID) - solres(field) = mechanical_solution(incInfo) + solres(field) = grid_mechanical_solution(incInfo) case(FIELD_THERMAL_ID) solres(field) = grid_thermal_spectral_solution(Delta_t) case(FIELD_DAMAGE_ID) @@ -344,7 +344,7 @@ program DAMASK_grid ! check solution and either advance or retry with smaller timestep if (all(solres(:)%converged .and. solres(:)%stagConverged)) then ! converged and acceptable solution found - call mechanical_updateCoords() + call grid_mechanical_updateCoords() Delta_t_prev = Delta_t cutBack = .false. guess = .true. ! start guessing after first converged (sub)inc @@ -389,7 +389,7 @@ program DAMASK_grid do field = 1, nActiveFields select case (ID(field)) case(FIELD_MECH_ID) - call mechanical_restartWrite() + call grid_mechanical_restartWrite() case(FIELD_THERMAL_ID) call grid_thermal_spectral_restartWrite() case(FIELD_DAMAGE_ID) @@ -538,7 +538,7 @@ function parseLoadsteps(load_steps) result(loadCases) transpose(loadCases(l)%rot%asMatrix()) if (loadCases(l)%r <= 0.0_pREAL) errorID = 833 - if (loadCases(l)%t < 0.0_pREAL) errorID = 834 + if (loadCases(l)%t <= 0.0_pREAL) errorID = 834 if (loadCases(l)%N < 1) errorID = 835 if (loadCases(l)%f_out < 1) errorID = 836 if (loadCases(l)%f_restart < 1) errorID = 839 diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 8efa371e3..7978a19fe 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -237,7 +237,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution) call SNESGetDM(SNES_damage,DM_damage,err_PETSc) CHKERRQ(err_PETSc) - call DMDAVecGetArrayF90(DM_damage,phi_PETSc,phi,err_PETSc) ! returns 0-indexed phi + call DMDAVecGetArrayReadF90(DM_damage,phi_PETSc,phi,err_PETSc) ! returns 0-indexed phi CHKERRQ(err_PETSc) stagNorm = maxval(abs(phi - phi_stagInc)) @@ -250,7 +250,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution) call homogenization_set_phi(reshape(phi,[product(cells(1:2))*cells3])) - call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc) + call DMDAVecRestoreArrayReadF90(DM_damage,phi_PETSc,phi,err_PETSc) CHKERRQ(err_PETSc) if (solution%converged) & diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 37921bf09..4725dafd1 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -218,7 +218,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) call SNESGetDM(SNES_thermal,DM_thermal,err_PETSc) CHKERRQ(err_PETSc) - call DMDAVecGetArrayF90(DM_thermal,T_PETSc,T,err_PETSc) ! returns 0-indexed T + call DMDAVecGetArrayReadF90(DM_thermal,T_PETSc,T,err_PETSc) ! returns 0-indexed T CHKERRQ(err_PETSc) stagNorm = maxval(abs(T - T_stagInc)) @@ -232,7 +232,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), & reshape(T-T_lastInc,[product(cells(1:2))*cells3])/Delta_t_) - call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc) + call DMDAVecRestoreArrayReadF90(DM_thermal,T_PETSc,T,err_PETSc) CHKERRQ(err_PETSc) if (solution%converged) & diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index fa09c0c44..488ad3c88 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -503,12 +503,12 @@ function G_hat_init() result(G_hat_) complex(pREAL), dimension(3,3), parameter :: delta = cmplx(math_I3,0.0_pREAL,pREAL) - allocate(G_hat_(3,3,3,3,cells1Red,cells(3),cells2),source=cmplx(.0_pReal,.0_pReal,pREAL)) + allocate(G_hat_(3,3,3,3,cells1Red,cells(3),cells2),source=cmplx(.0_pREAL,.0_pREAL,pREAL)) !$OMP PARALLEL DO PRIVATE(l,m,n,o,xi_norm_2) do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red if (any([i,j+cells2Offset,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 - xi_norm_2 = cmplx(abs(dot_product(xi1st(:,i,k,j), xi1st(:,i,k,j))),0.0_pReal,pREAL) + xi_norm_2 = cmplx(abs(dot_product(xi1st(:,i,k,j), xi1st(:,i,k,j))),0.0_pREAL,pREAL) if (xi_norm_2%re > 1.e-16_pREAL) then #ifndef __INTEL_COMPILER do concurrent(l=1:3, m=1:3, n=1:3, o=1:3) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index f4e3f1bbf..67656fd1b 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -41,7 +41,7 @@ module discretization_mesh !!!! BEGIN DEPRECATED !!!!! integer, public, protected :: & mesh_maxNips !< max number of IPs in any CP element -!!!! BEGIN DEPRECATED !!!!! +!!!! END DEPRECATED !!!!! DM, public :: geomMesh @@ -158,11 +158,6 @@ subroutine discretization_mesh_init() CHKERRQ(err_PETSc) ! Get initial nodal coordinates - call DMGetCoordinatesLocal(geomMesh,coords_node0,err_PETSc) - CHKERRQ(err_PETSc) - call VecGetArrayF90(coords_node0, mesh_node0_temp,err_PETSc) - CHKERRQ(err_PETSc) - mesh_maxNips = FEM_nQuadrature(dimPlex,p_i) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,p_i)%p) @@ -174,9 +169,14 @@ subroutine discretization_mesh_init() CHKERRQ(err_PETSc) end do + call DMGetCoordinatesLocal(geomMesh,coords_node0,err_PETSc) + CHKERRQ(err_PETSc) + call VecGetArrayReadF90(coords_node0, mesh_node0_temp,err_PETSc) + CHKERRQ(err_PETSc) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pREAL) mesh_node0(1:dimPlex,:) = reshape(mesh_node0_temp,[dimPlex,mesh_Nnodes]) - + call VecRestoreArrayReadF90(coords_node0, mesh_node0_temp,err_PETSc) + CHKERRQ(err_PETSc) call discretization_init(int(materialAt),& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 12ac922b1..1220c9196 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -379,7 +379,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc real(pREAL), dimension(:), pointer :: x_scal, pf_scal real(pREAL), dimension(cellDof), target :: f_scal PetscReal :: IcellJMat(dimPlex,dimPlex) - PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer + PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, dev_null, basisFieldDer PetscInt :: cellStart, cellEnd, cell, component, face, & qPt, basis, comp, cidx, & numFields, & @@ -396,10 +396,6 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc call DMGetLocalSection(dm_local,section,err_PETSc) CHKERRQ(err_PETSc) - call DMGetDS(dm_local,prob,err_PETSc) - CHKERRQ(err_PETSc) - call PetscDSGetTabulation(prob,0_pPETSCINT,basisField,basisFieldDer,err_PETSc) - CHKERRQ(err_PETSc) call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc) CHKERRQ(err_PETSc) call DMGetLocalVector(dm_local,x_local,err_PETSc) @@ -422,6 +418,11 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc !-------------------------------------------------------------------------------------------------- ! evaluate field derivatives + call DMGetDS(dm_local,prob,err_PETSc) + CHKERRQ(err_PETSc) + call PetscDSGetTabulation(prob,0_pPETSCINT,dev_null,basisFieldDer,err_PETSc) + CHKERRQ(err_PETSc) + do cell = cellStart, cellEnd-1_pPETSCINT !< loop over all elements call PetscSectionGetNumFields(section,numFields,err_PETSc) @@ -496,6 +497,8 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,err_PETSc) CHKERRQ(err_PETSc) end do + call PetscDSRestoreTabulation(prob,0_pPETSCINT,dev_null,basisFieldDer,err_PETSc) + CHKERRQ(err_PETSc) call DMRestoreLocalVector(dm_local,x_local,err_PETSc) CHKERRQ(err_PETSc) @@ -521,8 +524,8 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA PetscReal, dimension(3,3) :: F, FAvg, FInv PetscReal :: detJ - PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & - pV0, pCellJ, pInvcellJ + PetscReal, dimension(:), pointer :: dev_null, basisFieldDer, & + pV0, pCellJ, pInvcellJ real(pREAL), dimension(:), pointer :: pK_e, x_scal @@ -546,7 +549,6 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P CHKERRQ(err_PETSc) call DMGetDS(dm_local,prob,err_PETSc) CHKERRQ(err_PETSc) - call PetscDSGetTabulation(prob,0_pPETSCINT,basisField,basisFieldDer,err_PETSc) call DMGetLocalSection(dm_local,section,err_PETSc) CHKERRQ(err_PETSc) call DMGetGlobalSection(dm_local,gSection,err_PETSc) @@ -569,6 +571,9 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P end if end if end do; end do + + call PetscDSGetTabulation(prob,0_pPETSCINT,dev_null,basisFieldDer,err_PETSc) + CHKERRQ(err_PETSc) call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc) CHKERRQ(err_PETSc) do cell = cellStart, cellEnd-1 !< loop over all elements @@ -642,6 +647,8 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P CHKERRQ(err_PETSc) call DMRestoreLocalVector(dm_local,x_local,err_PETSc) CHKERRQ(err_PETSc) + call PetscDSRestoreTabulation(prob,0_pPETSCINT,dev_null,basisFieldDer,err_PETSc) + CHKERRQ(err_PETSc) !-------------------------------------------------------------------------------------------------- ! apply boundary conditions @@ -781,7 +788,7 @@ subroutine FEM_mechanical_updateCoords() cellStart, cellEnd, c, n PetscSection :: section PetscQuadrature :: mechQuad - PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & + PetscReal, dimension(:), pointer :: basisField, dev_null, & nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes) real(pREAL), dimension(:), pointer :: x_scal @@ -815,7 +822,7 @@ subroutine FEM_mechanical_updateCoords() ! write ip displacements call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc) CHKERRQ(err_PETSc) - call PetscDSGetTabulation(mechQuad,0_pPETSCINT,basisField,basisFieldDer,err_PETSc) + call PetscDSGetTabulation(mechQuad,0_pPETSCINT,basisField,dev_null,err_PETSc) CHKERRQ(err_PETSc) allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pREAL) do c=cellStart,cellEnd-1_pPETSCINT @@ -842,6 +849,8 @@ subroutine FEM_mechanical_updateCoords() call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature])) call DMRestoreLocalVector(dm_local,x_local,err_PETSc) CHKERRQ(err_PETSc) + call PetscDSRestoreTabulation(mechQuad,0_pPETSCINT,basisField,dev_null,err_PETSc) + CHKERRQ(err_PETSc) end subroutine FEM_mechanical_updateCoords