diff --git a/PRIVATE b/PRIVATE index 8b9836aaf..4b5636231 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 8b9836aaf5798727d96a316a2eb27df03bbd2d07 +Subproject commit 4b563623162f033b11d9d0ac95fc421de279f690 diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 8a12bef76..e441339cb 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -72,7 +72,6 @@ contains !-------------------------------------------------------------------------------------------------- subroutine CPFEM_initAll - call parallelization_init call DAMASK_interface_init call prec_init call IO_init diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 360c911ff..7abb0e746 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -7,20 +7,6 @@ !> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief Interfaces DAMASK with MSC.Marc -!> @details Usage: -!> @details - choose material as hypela2 -!> @details - set statevariable 2 to index of homogenization -!> @details - set statevariable 3 to index of microstructure -!> @details - use nonsymmetric option for solver (e.g. direct profile or multifrontal sparse, the latter seems to be faster!) -!> @details - in case of ddm (domain decomposition) a SYMMETRIC solver has to be used, i.e uncheck "non-symmetric" -!> @details Marc subroutines used: -!> @details - hypela2 -!> @details - uedinc -!> @details - flux -!> @details - quit -!> @details Marc common blocks included: -!> @details - concom: lovl, inc -!> @details - creeps: timinc !-------------------------------------------------------------------------------------------------- #define QUOTE(x) #x #define PASTE(x,y) x ## y @@ -65,14 +51,8 @@ subroutine DAMASK_interface_init print'(/,a)', ' Version: '//DAMASKVERSION - ! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md -#if __INTEL_COMPILER >= 1800 - print'(/,a)', ' Compiled with: '//compiler_version() - print'(a)', ' Compiler options: '//compiler_options() -#else - print'(/,a,i4.4,a,i8.8)', ' Compiled with Intel fortran version :', __INTEL_COMPILER,& - ', build date :', __INTEL_COMPILER_BUILD_DATE -#endif + print'(/,a)', ' Compiled with: '//compiler_version() + print'(a)', ' Compiler options: '//compiler_options() print'(/,a)', ' Compiled on: '//__DATE__//' at '//__TIME__ @@ -239,7 +219,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & real(pReal), dimension(6) :: stress real(pReal), dimension(6,6) :: ddsdde integer :: computationMode, i, cp_en, node, CPnodeID - integer(4) :: defaultNumThreadsInt !< default value set by Marc + integer(pI32) :: defaultNumThreadsInt !< default value set by Marc integer(pInt), save :: & theInc = -1_pInt, & !< needs description @@ -250,13 +230,13 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & logical, save :: & lastIncConverged = .false., & !< needs description outdatedByNewInc = .false., & !< needs description - CPFEM_init_done = .false., & !< remember whether init has been done already + CPFEM_init_done = .false., & !< remember whether init has been done already debug_basic = .true. class(tNode), pointer :: & debug_Marc ! pointer to Marc debug options if(debug_basic) then - print'(a,/,i8,i8,i2)', ' MSC.MARC information on shape of element(2), IP:', m, nn + print'(a,/,i8,i8,i2)', ' MSC.Marc information on shape of element(2), IP:', m, nn print'(a,2(i1))', ' Jacobian: ', ngens,ngens print'(a,i1)', ' Direct stress: ', ndi print'(a,i1)', ' Shear stress: ', nshear @@ -271,7 +251,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & endif defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc - call omp_set_num_threads(1) ! no openMP + call omp_set_num_threads(1_pI32) ! no openMP if (.not. CPFEM_init_done) then CPFEM_init_done = .true. @@ -379,14 +359,28 @@ subroutine flux(f,ts,n,time) subroutine uedinc(inc,incsub) use prec use CPFEM + use discretization_marc implicit none integer, intent(in) :: inc, incsub + integer :: n, nqncomp, nqdatatype integer, save :: inc_written + real(pReal), allocatable, dimension(:,:) :: d_n #include QUOTE(PASTE(./marc/include/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment) + if (inc > inc_written) then + allocate(d_n(3,size(mesh_FEM2DAMASK_node))) + do n = 1, size(d_n,2) + if (mesh_FEM2DAMASK_node(n) /= -1) then + call nodvar(1,n,d_n(1:3,mesh_FEM2DAMASK_node(n)),nqncomp,nqdatatype) + if(nqncomp == 2) d_n(3,mesh_FEM2DAMASK_node(n)) = 0.0_pReal + endif + enddo + + call discretization_marc_UpdateNodeAndIpCoords(d_n) call CPFEM_results(inc,cptim) + inc_written = inc endif diff --git a/src/element.f90 b/src/element.f90 index 722a7fd96..5b7af36bd 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -686,7 +686,7 @@ module element 1, 5,11, 7, 8,12,15,14, & 5, 2, 6,11,12, 9,13,15, & 7,11, 6, 3,14,15,13,10, & - 8,12,15, 4, 4, 9,13,10 & + 8,12,15,14, 4, 9,13,10 & #if !defined(__GFORTRAN__) ],shape(CELL6)) #else diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index 675e57bd3..d92623215 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -20,6 +20,14 @@ module discretization_marc implicit none private + real(pReal), public, protected :: & + mesh_unitlength !< physical length of one unit in mesh MD: needs systematic_name + + integer, dimension(:), allocatable, public, protected :: & + mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID MD: Needs systematic name + mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID MD: needs systematic_name + + type tCellNodeDefinition integer, dimension(:,:), allocatable :: parents integer, dimension(:,:), allocatable :: weights @@ -27,15 +35,12 @@ module discretization_marc type(tCellNodeDefinition), dimension(:), allocatable :: cellNodeDefinition - real(pReal), public, protected :: & - mesh_unitlength !< physical length of one unit in mesh - - integer, dimension(:), allocatable, public :: & - mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID - mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID + integer, dimension(:,:,:), allocatable :: & + connectivity_cell !< cell connectivity for each element,ip/cell public :: & - discretization_marc_init + discretization_marc_init, & + discretization_marc_UpdateNodeAndIpCoords contains @@ -45,25 +50,22 @@ contains !-------------------------------------------------------------------------------------------------- subroutine discretization_marc_init - real(pReal), dimension(:,:), allocatable :: & + real(pReal), dimension(:,:), allocatable :: & node0_elem, & !< node x,y,z coordinates (initially!) node0_cell type(tElement) :: elem - integer, dimension(:), allocatable :: & + integer, dimension(:), allocatable :: & materialAt integer:: & - Nnodes, & !< total number of nodes in the mesh Nelems, & !< total number of elements in the mesh debug_e, debug_i - real(pReal), dimension(:,:), allocatable :: & + real(pReal), dimension(:,:), allocatable :: & IP_reshaped - integer,dimension(:,:,:), allocatable :: & - connectivity_cell !< cell connectivity for each element,ip/cell - integer, dimension(:,:), allocatable :: & + integer, dimension(:,:), allocatable :: & connectivity_elem - real(pReal), dimension(:,:,:,:),allocatable :: & + real(pReal), dimension(:,:,:,:), allocatable :: & unscaledNormals class(tNode), pointer :: & @@ -90,18 +92,14 @@ subroutine discretization_marc_init allocate(cellNodeDefinition(elem%nNodes-1)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) + call buildCells(connectivity_cell,cellNodeDefinition,& elem,connectivity_elem) - allocate(node0_cell(3,maxval(connectivity_cell))) - call buildCellNodes(node0_cell,& - cellNodeDefinition,node0_elem) - allocate(IP_reshaped(3,elem%nIPs*nElems),source=0.0_pReal) - call buildIPcoordinates(IP_reshaped,reshape(connectivity_cell,[elem%NcellNodesPerCell,& - elem%nIPs*nElems]),node0_cell) + node0_cell = buildCellNodes(node0_elem) + + IP_reshaped = buildIPcoordinates(node0_cell) - call discretization_init(materialAt,& - IP_reshaped,& - node0_cell) + call discretization_init(materialAt, IP_reshaped, node0_cell) call writeGeometry(elem,connectivity_elem,& reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),& @@ -109,55 +107,65 @@ subroutine discretization_marc_init !-------------------------------------------------------------------------------------------------- ! geometry information required by the nonlocal CP model - call geometry_plastic_nonlocal_setIPvolume(IPvolume(elem,node0_cell,connectivity_cell)) - unscaledNormals = IPareaNormal(elem,nElems,connectivity_cell,node0_cell) + call geometry_plastic_nonlocal_setIPvolume(IPvolume(elem,node0_cell)) + unscaledNormals = IPareaNormal(elem,nElems,node0_cell) call geometry_plastic_nonlocal_setIParea(norm2(unscaledNormals,1)) call geometry_plastic_nonlocal_setIPareaNormal(unscaledNormals/spread(norm2(unscaledNormals,1),1,3)) - call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(elem,connectivity_cell)) + call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(elem)) call geometry_plastic_nonlocal_results end subroutine discretization_marc_init +!-------------------------------------------------------------------------------------------------- +!> @brief Calculate and set current nodal and IP positions (including cell nodes) +!-------------------------------------------------------------------------------------------------- +subroutine discretization_marc_UpdateNodeAndIpCoords(d_n) + + real(pReal), dimension(:,:), intent(in) :: d_n + + real(pReal), dimension(:,:), allocatable :: node_cell + + + node_cell = buildCellNodes(discretization_NodeCoords0(1:3,1:maxval(mesh_FEM2DAMASK_node)) + d_n) + + call discretization_setNodeCoords(node_cell) + call discretization_setIPcoords(buildIPcoordinates(node_cell)) + +end subroutine discretization_marc_UpdateNodeAndIpCoords + + !-------------------------------------------------------------------------------------------------- !> @brief Write all information needed for the DADF5 geometry !-------------------------------------------------------------------------------------------------- subroutine writeGeometry(elem, & - connectivity_elem,connectivity_cell, & + connectivity_elem,connectivity_cell_reshaped, & coordinates_nodes,coordinates_points) type(tElement), intent(in) :: & elem integer, dimension(:,:), intent(in) :: & connectivity_elem, & - connectivity_cell + connectivity_cell_reshaped real(pReal), dimension(:,:), intent(in) :: & coordinates_nodes, & coordinates_points - integer, dimension(:,:), allocatable :: & - connectivity_temp - real(pReal), dimension(:,:), allocatable :: & - coordinates_temp call results_openJobFile call results_closeGroup(results_addGroup('geometry')) - connectivity_temp = connectivity_elem - call results_writeDataset('geometry',connectivity_temp,'T_e',& + call results_writeDataset('geometry',connectivity_elem,'T_e',& 'connectivity of the elements','-') - connectivity_temp = connectivity_cell - call results_writeDataset('geometry',connectivity_temp,'T_c', & + call results_writeDataset('geometry',connectivity_cell_reshaped,'T_c', & 'connectivity of the cells','-') call results_addAttribute('VTK_TYPE',elem%vtkType,'geometry/T_c') - coordinates_temp = coordinates_nodes - call results_writeDataset('geometry',coordinates_temp,'x_n', & + call results_writeDataset('geometry',coordinates_nodes,'x_n', & 'initial coordinates of the nodes','m') - coordinates_temp = coordinates_points - call results_writeDataset('geometry',coordinates_temp,'x_p', & + call results_writeDataset('geometry',coordinates_points,'x_p', & 'initial coordinates of the materialpoints (cell centers)','m') call results_closeJobFile @@ -186,13 +194,13 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,materialAt) nElems integer, dimension(:), allocatable :: & matNumber !< material numbers for hypoelastic material - character(len=pStringLen), dimension(:), allocatable :: inputFile !< file content, separated per lines - character(len=pStringLen), dimension(:), allocatable :: & + inputFile, & !< file content, separated per lines nameElemSet integer, dimension(:,:), allocatable :: & mapElemSet !< list of elements in elementSet + inputFile = IO_readlines(trim(getSolverJobName())//trim(InputFileExtension)) call inputRead_fileFormat(fileFormatVersion, & inputFile) @@ -588,20 +596,20 @@ subroutine inputRead_elemType(elem, & mapElemtype = 1 ! Two-dimensional Plane Strain Triangle case ( '125') ! 155, 128 (need test) mapElemtype = 2 ! Two-dimensional Plane Strain triangle (155: cubic shape function, 125/128: second order isoparametric) - !case ( '11') ! need test - ! mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain + case ( '11') + mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain case ( '27') mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral case ( '54') mapElemtype = 5 ! Plane Strain, Eight-node Distorted Quadrilateral with reduced integration - !case ( '134') ! need test - ! mapElemtype = 6 ! Three-dimensional Four-node Tetrahedron + case ( '134') + mapElemtype = 6 ! Three-dimensional Four-node Tetrahedron !case ( '157') ! need test ! mapElemtype = 7 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations - !case ( '127') ! need test - ! mapElemtype = 8 ! Three-dimensional Ten-node Tetrahedron - !case ( '136') ! need test - ! mapElemtype = 9 ! Three-dimensional Arbitrarily Distorted Pentahedral + case ( '127') + mapElemtype = 8 ! Three-dimensional Ten-node Tetrahedron + case ( '136') + mapElemtype = 9 ! Three-dimensional Arbitrarily Distorted Pentahedral case ( '117') ! 123 (need test) mapElemtype = 10 ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration case ( '7') @@ -727,11 +735,11 @@ end subroutine inputRead_material !-------------------------------------------------------------------------------------------------- !> @brief Calculates cell node coordinates from element node coordinates !-------------------------------------------------------------------------------------------------- -pure subroutine buildCells(connectivity_cell,cellNodeDefinition, & +pure subroutine buildCells(connectivity,definition, & elem,connectivity_elem) - type(tCellNodeDefinition), dimension(:), intent(out) :: cellNodeDefinition ! definition of cell nodes for increasing number of parents - integer, dimension(:,:,:),intent(out) :: connectivity_cell + type(tCellNodeDefinition), dimension(:), intent(out) :: definition ! definition of cell nodes for increasing number of parents + integer, dimension(:,:,:),intent(out) :: connectivity type(tElement), intent(in) :: elem ! element definition integer, dimension(:,:), intent(in) :: connectivity_elem ! connectivity of the elements @@ -745,7 +753,7 @@ pure subroutine buildCells(connectivity_cell,cellNodeDefinition, & !--------------------------------------------------------------------------------------------------- ! initialize global connectivity to negative local connectivity - connectivity_cell = -spread(elem%cell,3,Nelem) ! local cell node ID + connectivity = -spread(elem%cell,3,Nelem) ! local cell node ID !--------------------------------------------------------------------------------------------------- ! set connectivity of cell nodes that coincide with FE nodes (defined by 1 parent node) @@ -753,9 +761,7 @@ pure subroutine buildCells(connectivity_cell,cellNodeDefinition, & do e = 1, Nelem do c = 1, elem%NcellNodes realNode: if (count(elem%cellNodeParentNodeWeights(:,c) /= 0) == 1) then - where(connectivity_cell(:,:,e) == -c) - connectivity_cell(:,:,e) = connectivity_elem(c,e) - end where + where(connectivity(:,:,e) == -c) connectivity(:,:,e) = connectivity_elem(c,e) endif realNode enddo enddo @@ -821,8 +827,8 @@ pure subroutine buildCells(connectivity_cell,cellNodeDefinition, & enddo i = uniqueRows(candidates_global(1:2*nParentNodes,:)) - allocate(cellNodeDefinition(nParentNodes-1)%parents(i,nParentNodes)) - allocate(cellNodeDefinition(nParentNodes-1)%weights(i,nParentNodes)) + allocate(definition(nParentNodes-1)%parents(i,nParentNodes)) + allocate(definition(nParentNodes-1)%weights(i,nParentNodes)) i = 1 n = 1 @@ -836,15 +842,15 @@ pure subroutine buildCells(connectivity_cell,cellNodeDefinition, & do while (n+j<= size(candidates_local)*Nelem) if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit - where (connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) == -candidates_global(nParentNodes*2+2,n+j)) ! still locally defined - connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) = nCellNode + 1 ! gets current new cell node id + where (connectivity(:,:,candidates_global(nParentNodes*2+1,n+j)) == -candidates_global(nParentNodes*2+2,n+j)) ! still locally defined + connectivity(:,:,candidates_global(nParentNodes*2+1,n+j)) = nCellNode + 1 ! gets current new cell node id end where j = j+1 enddo nCellNode = nCellNode + 1 - cellNodeDefinition(nParentNodes-1)%parents(i,:) = parentsAndWeights(:,1) - cellNodeDefinition(nParentNodes-1)%weights(i,:) = parentsAndWeights(:,2) + definition(nParentNodes-1)%parents(i,:) = parentsAndWeights(:,1) + definition(nParentNodes-1)%weights(i,:) = parentsAndWeights(:,2) i = i + 1 n = n+j enddo @@ -884,55 +890,62 @@ end subroutine buildCells !-------------------------------------------------------------------------------------------------- !> @brief Calculates cell node coordinates from element node coordinates !-------------------------------------------------------------------------------------------------- -pure subroutine buildCellNodes(node_cell, & - definition,node_elem) +pure function buildCellNodes(node_elem) - real(pReal), dimension(:,:), intent(out) :: node_cell !< cell node coordinates - type(tCellNodeDefinition), dimension(:), intent(in) :: definition !< cell node definition (weights and parents) real(pReal), dimension(:,:), intent(in) :: node_elem !< element nodes + real(pReal), dimension(:,:), allocatable :: buildCellNodes !< cell node coordinates integer :: i, j, k, n - n = size(node_elem,2) - node_cell(:,1:n) = node_elem !< initial nodes coincide with element nodes - do i = 1, size(cellNodeDefinition,1) + allocate(buildCellNodes(3,maxval(connectivity_cell))) + n = size(node_elem,2) + buildCellNodes(:,1:n) = node_elem !< initial nodes coincide with element nodes + + do i = 1, size(cellNodeDefinition) do j = 1, size(cellNodeDefinition(i)%parents,1) n = n+1 - node_cell(:,n) = 0.0_pReal + buildCellNodes(:,n) = 0.0_pReal do k = 1, size(cellNodeDefinition(i)%parents,2) - node_cell(:,n) = node_cell(:,n) & - + node_cell(:,definition(i)%parents(j,k)) * real(definition(i)%weights(j,k),pReal) + buildCellNodes(:,n) = buildCellNodes(:,n) & + + buildCellNodes(:,cellNodeDefinition(i)%parents(j,k)) & + * real(cellNodeDefinition(i)%weights(j,k),pReal) enddo - node_cell(:,n) = node_cell(:,n)/real(sum(definition(i)%weights(j,:)),pReal) + buildCellNodes(:,n) = buildCellNodes(:,n)/real(sum(cellNodeDefinition(i)%weights(j,:)),pReal) enddo enddo -end subroutine buildCellNodes +end function buildCellNodes !-------------------------------------------------------------------------------------------------- !> @brief Calculates IP coordinates as center of cell !-------------------------------------------------------------------------------------------------- -pure subroutine buildIPcoordinates(IPcoordinates, & - connectivity_cell,node_cell) +pure function buildIPcoordinates(node_cell) - real(pReal), dimension(:,:), intent(out):: IPcoordinates !< cell-center/IP coordinates - integer, dimension(:,:), intent(in) :: connectivity_cell !< connectivity for each cell - real(pReal), dimension(:,:), intent(in) :: node_cell !< cell node coordinates + real(pReal), dimension(:,:), intent(in) :: node_cell !< cell node coordinates + real(pReal), dimension(:,:), allocatable :: buildIPcoordinates !< cell-center/IP coordinates - integer :: i, n + integer, dimension(:,:), allocatable :: connectivity_cell_reshaped + integer :: i, n, NcellNodesPerCell,Ncells - do i = 1, size(connectivity_cell,2) - IPcoordinates(:,i) = 0.0_pReal - do n = 1, size(connectivity_cell,1) - IPcoordinates(:,i) = IPcoordinates(:,i) & - + node_cell(:,connectivity_cell(n,i)) + + NcellNodesPerCell = size(connectivity_cell,1) + Ncells = size(connectivity_cell,2)*size(connectivity_cell,3) + connectivity_cell_reshaped = reshape(connectivity_cell,[NcellNodesPerCell,Ncells]) + + allocate(buildIPcoordinates(3,Ncells)) + + do i = 1, size(connectivity_cell_reshaped,2) + buildIPcoordinates(:,i) = 0.0_pReal + do n = 1, size(connectivity_cell_reshaped,1) + buildIPcoordinates(:,i) = buildIPcoordinates(:,i) & + + node_cell(:,connectivity_cell_reshaped(n,i)) enddo - IPcoordinates(:,i) = IPcoordinates(:,i)/real(size(connectivity_cell,1),pReal) + buildIPcoordinates(:,i) = buildIPcoordinates(:,i)/real(size(connectivity_cell_reshaped,1),pReal) enddo -end subroutine buildIPcoordinates +end function buildIPcoordinates !--------------------------------------------------------------------------------------------------- @@ -940,50 +953,50 @@ end subroutine buildIPcoordinates !> @details The IP volume is calculated differently depending on the cell type. !> 2D cells assume an element depth of 1.0 !--------------------------------------------------------------------------------------------------- -pure function IPvolume(elem,node,connectivity) +pure function IPvolume(elem,node) type(tElement), intent(in) :: elem real(pReal), dimension(:,:), intent(in) :: node - integer, dimension(:,:,:), intent(in) :: connectivity - real(pReal), dimension(elem%nIPs,size(connectivity,3)) :: IPvolume + real(pReal), dimension(elem%nIPs,size(connectivity_cell,3)) :: IPvolume real(pReal), dimension(3) :: x0,x1,x2,x3,x4,x5,x6,x7 integer :: e,i - do e = 1,size(connectivity,3) + + do e = 1,size(connectivity_cell,3) do i = 1,elem%nIPs select case (elem%cellType) case (1) ! 2D 3node - IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), & - node(1:3,connectivity(2,i,e)), & - node(1:3,connectivity(3,i,e))) + IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity_cell(1,i,e)), & + node(1:3,connectivity_cell(2,i,e)), & + node(1:3,connectivity_cell(3,i,e))) case (2) ! 2D 4node - IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), & ! assume planar shape, division in two triangles suffices - node(1:3,connectivity(2,i,e)), & - node(1:3,connectivity(3,i,e))) & - + math_areaTriangle(node(1:3,connectivity(3,i,e)), & - node(1:3,connectivity(4,i,e)), & - node(1:3,connectivity(1,i,e))) + IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity_cell(1,i,e)), & ! assume planar shape, division in two triangles suffices + node(1:3,connectivity_cell(2,i,e)), & + node(1:3,connectivity_cell(3,i,e))) & + + math_areaTriangle(node(1:3,connectivity_cell(3,i,e)), & + node(1:3,connectivity_cell(4,i,e)), & + node(1:3,connectivity_cell(1,i,e))) case (3) ! 3D 4node - IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity(1,i,e)), & - node(1:3,connectivity(2,i,e)), & - node(1:3,connectivity(3,i,e)), & - node(1:3,connectivity(4,i,e))) + IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity_cell(1,i,e)), & + node(1:3,connectivity_cell(2,i,e)), & + node(1:3,connectivity_cell(3,i,e)), & + node(1:3,connectivity_cell(4,i,e))) case (4) ! 3D 8node ! J. Grandy, Efficient Calculation of Volume of Hexahedral Cells ! Lawrence Livermore National Laboratory ! https://www.osti.gov/servlets/purl/632793 - x0 = node(1:3,connectivity(1,i,e)) - x1 = node(1:3,connectivity(2,i,e)) - x2 = node(1:3,connectivity(4,i,e)) - x3 = node(1:3,connectivity(3,i,e)) - x4 = node(1:3,connectivity(5,i,e)) - x5 = node(1:3,connectivity(6,i,e)) - x6 = node(1:3,connectivity(8,i,e)) - x7 = node(1:3,connectivity(7,i,e)) + x0 = node(1:3,connectivity_cell(1,i,e)) + x1 = node(1:3,connectivity_cell(2,i,e)) + x2 = node(1:3,connectivity_cell(4,i,e)) + x3 = node(1:3,connectivity_cell(3,i,e)) + x4 = node(1:3,connectivity_cell(5,i,e)) + x5 = node(1:3,connectivity_cell(6,i,e)) + x6 = node(1:3,connectivity_cell(8,i,e)) + x7 = node(1:3,connectivity_cell(7,i,e)) IPvolume(i,e) = dot_product((x7-x1)+(x6-x0),math_cross((x7-x2), (x3-x0))) & + dot_product((x6-x0), math_cross((x7-x2)+(x5-x0),(x7-x4))) & + dot_product((x7-x1), math_cross((x5-x0), (x7-x4)+(x3-x0))) @@ -998,11 +1011,10 @@ end function IPvolume !-------------------------------------------------------------------------------------------------- !> @brief calculation of IP interface areas !-------------------------------------------------------------------------------------------------- -pure function IPareaNormal(elem,nElem,connectivity,node) +pure function IPareaNormal(elem,nElem,node) type(tElement), intent(in) :: elem integer, intent(in) :: nElem - integer, dimension(:,:,:), intent(in) :: connectivity real(pReal), dimension(:,:), intent(in) :: node real(pReal), dimension(3,elem%nIPneighbors,elem%nIPs,nElem) :: ipAreaNormal @@ -1015,7 +1027,7 @@ pure function IPareaNormal(elem,nElem,connectivity,node) do e = 1,nElem do i = 1,elem%nIPs do f = 1,elem%nIPneighbors - nodePos = node(1:3,connectivity(elem%cellface(1:m,f),i,e)) + nodePos = node(1:3,connectivity_cell(elem%cellface(1:m,f),i,e)) select case (elem%cellType) case (1,2) ! 2D 3 or 4 node @@ -1046,23 +1058,22 @@ end function IPareaNormal !-------------------------------------------------------------------------------------------------- !> @brief IP neighborhood !-------------------------------------------------------------------------------------------------- -function IPneighborhood(elem,connectivity) +function IPneighborhood(elem) type(tElement), intent(in) :: elem ! definition of the element in use - integer, dimension(:,:,:), intent(in) :: connectivity ! cell connectivity integer, dimension(3,size(elem%cellFace,2), & - size(connectivity,2),size(connectivity,3)) :: IPneighborhood ! neighboring IPs as [element ID, IP ID, face ID] + size(connectivity_cell,2),size(connectivity_cell,3)) :: IPneighborhood ! neighboring IPs as [element ID, IP ID, face ID] integer, dimension(size(elem%cellFace,1)+3,& - size(elem%cellFace,2)*size(connectivity,2)*size(connectivity,3)) :: face - integer, dimension(size(connectivity,1)) :: myConnectivity - integer, dimension(size(elem%cellFace,1)) :: face_unordered + size(elem%cellFace,2)*size(connectivity_cell,2)*size(connectivity_cell,3)) :: face + integer, dimension(size(connectivity_cell,1)) :: myConnectivity + integer, dimension(size(elem%cellFace,1)) :: face_unordered integer :: e,i,f,n,c,s c = 0 - do e = 1, size(connectivity,3) - do i = 1, size(connectivity,2) - myConnectivity = connectivity(:,i,e) + do e = 1, size(connectivity_cell,3) + do i = 1, size(connectivity_cell,2) + myConnectivity = connectivity_cell(:,i,e) do f = 1, size(elem%cellFace,2) c = c + 1 face_unordered = myConnectivity(elem%cellFace(:,f)) diff --git a/src/parallelization.f90 b/src/parallelization.f90 index edf56407b..8413ba825 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -9,9 +9,8 @@ module parallelization #ifdef PETSc #include use petscsys -#endif !$ use OMP_LIB - +#endif use prec implicit none @@ -21,6 +20,7 @@ module parallelization worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only) worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only) +#ifdef PETSc public :: & parallelization_init @@ -32,16 +32,12 @@ contains subroutine parallelization_init integer :: err, typeSize -!$ integer :: got_env, OMP_NUM_THREADS, threadLevel +!$ integer :: got_env, threadLevel +!$ integer(pI32) :: OMP_NUM_THREADS !$ character(len=6) NumThreadsString -#ifdef PETSc + + PetscErrorCode :: petsc_err - -#else - print'(/,a)', ' <<<+- parallelization init -+>>>'; flush(OUTPUT_UNIT) -#endif - -#ifdef PETSc #ifdef _OPENMP ! If openMP is enabled, check if the MPI libary supports it and initialize accordingly. ! Otherwise, the first call to PETSc will do the initialization. @@ -64,7 +60,7 @@ subroutine parallelization_init #endif CHKERRQ(petsc_err) -call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,err) + call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,err) if (err /= 0) error stop 'Could not determine worldrank' if (worldrank == 0) print'(/,a)', ' <<<+- parallelization init -+>>>' @@ -80,7 +76,6 @@ call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,err) call MPI_Type_size(MPI_DOUBLE,typeSize,err) if (err /= 0) error stop 'Could not determine MPI real size' if (typeSize*8 /= storage_size(0.0_pReal)) error stop 'Mismatch between MPI and DAMASK real' -#endif if (worldrank /= 0) then close(OUTPUT_UNIT) ! disable output @@ -102,5 +97,6 @@ call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,err) !$ call omp_set_num_threads(OMP_NUM_THREADS) end subroutine parallelization_init +#endif end module parallelization