From e0c2dbc948d39790b127980ec0defcd88bc09a90 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 6 Jul 2020 16:33:48 +0200 Subject: [PATCH 1/6] look for new name --- env/DAMASK.csh | 6 +++--- env/DAMASK.sh | 4 ++-- env/DAMASK.zsh | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/env/DAMASK.csh b/env/DAMASK.csh index 98693d6b2..5529eac23 100644 --- a/env/DAMASK.csh +++ b/env/DAMASK.csh @@ -9,8 +9,8 @@ source $ENV_ROOT/CONFIG set path = ($DAMASK_ROOT/bin $path) -set SOLVER=`which DAMASK_spectral` -if ( "x$DAMASK_NUM_THREADS" == "x" ) then +set SOLVER=`which DAMASK_grid` +if ( "x$DAMASK_NUM_THREADS" == "x" ) then set DAMASK_NUM_THREADS=1 endif @@ -30,7 +30,7 @@ if ( $?prompt ) then echo echo Using environment with ... echo "DAMASK $DAMASK_ROOT" - echo "Grid Solver $SOLVER" + echo "Grid Solver $SOLVER" if ( $?PETSC_DIR) then echo "PETSc location $PETSC_DIR" endif diff --git a/env/DAMASK.sh b/env/DAMASK.sh index 5c3d2ba85..2151e842b 100644 --- a/env/DAMASK.sh +++ b/env/DAMASK.sh @@ -35,7 +35,7 @@ cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd PATH=${DAMASK_ROOT}/bin:$PATH -SOLVER=$(type -p DAMASK_spectral || true 2>/dev/null) +SOLVER=$(type -p DAMASK_grid || true 2>/dev/null) [ "x$SOLVER" == "x" ] && SOLVER=$(blink 'Not found!') [ "x$DAMASK_NUM_THREADS" == "x" ] && DAMASK_NUM_THREADS=1 @@ -56,7 +56,7 @@ if [ ! -z "$PS1" ]; then echo echo Using environment with ... echo "DAMASK $DAMASK_ROOT $BRANCH" - echo "Grid Solver $SOLVER" + echo "Grid Solver $SOLVER" 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 831268a7e..377aa5304 100644 --- a/env/DAMASK.zsh +++ b/env/DAMASK.zsh @@ -27,7 +27,7 @@ cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd PATH=${DAMASK_ROOT}/bin:$PATH -SOLVER=$(which DAMASK_spectral || true 2>/dev/null) +SOLVER=$(which DAMASK_grid || true 2>/dev/null) [[ "x$SOLVER" == "x" ]] && SOLVER=$(blink 'Not found!') [[ "x$DAMASK_NUM_THREADS" == "x" ]] && DAMASK_NUM_THREADS=1 @@ -48,7 +48,7 @@ if [ ! -z "$PS1" ]; then echo echo "Using environment with ..." echo "DAMASK $DAMASK_ROOT $BRANCH" - echo "Grid Solver $SOLVER" + echo "Grid Solver $SOLVER" if [ "x$PETSC_DIR" != "x" ]; then echo -n "PETSc location " [ -d $PETSC_DIR ] && echo $PETSC_DIR || blink $PETSC_DIR From dbca47f113b64ad87ba1391a8a43bdd7bdf669e0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 Jul 2020 12:38:21 +0200 Subject: [PATCH 2/6] better handling of corner cases copy and paste from python code --- src/quaternions.f90 | 2 ++ src/rotations.f90 | 6 +++--- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index c337bd681..a396f7f67 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -464,6 +464,8 @@ subroutine selfTest real(pReal), dimension(4) :: qu type(quaternion) :: q, q_2 + if(dNeq(abs(P),1.0_pReal)) call IO_error(0,ext_msg='P not in {-1,+1}') + call random_number(qu) qu = (qu-0.5_pReal) * 2.0_pReal q = quaternion(qu) diff --git a/src/rotations.f90 b/src/rotations.f90 index f95c54e98..7879523d2 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -581,7 +581,7 @@ pure function om2eu(om) result(eu) real(pReal), dimension(3) :: eu real(pReal) :: zeta - if (abs(om(3,3)) < 1.0_pReal) then + if(dNeq(abs(om(3,3)),1.0_pReal,1.e-9_pReal)) then zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal) eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), & acos(om(3,3)), & @@ -589,8 +589,8 @@ pure function om2eu(om) result(eu) else eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ] end if - - where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) + where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal + where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) end function om2eu From 9f79faf819026061cb16e2e2aa210bc1bcbb5e68 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 Jul 2020 22:48:08 +0200 Subject: [PATCH 3/6] using tolerances as in python results in invalid operations --- src/rotations.f90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/rotations.f90 b/src/rotations.f90 index 7879523d2..5d865c2c3 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -574,6 +574,7 @@ end function om2qu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief orientation matrix to Euler angles +!> @details Two step check for special cases to avoid invalid operations (not needed for python) !--------------------------------------------------------------------------------------------------- pure function om2eu(om) result(eu) @@ -581,11 +582,13 @@ pure function om2eu(om) result(eu) real(pReal), dimension(3) :: eu real(pReal) :: zeta - if(dNeq(abs(om(3,3)),1.0_pReal,1.e-9_pReal)) then + if (dNeq(abs(om(3,3)),1.0_pReal,1.e-5_pReal)) then zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal) eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), & acos(om(3,3)), & atan2(om(1,3)*zeta, om(2,3)*zeta)] + elseif(dNeq(abs(om(3,3)),1.0_pReal,1.e-8_pReal)) then + eu = [PI*.75_pReal,acos(math_clip(om(3,3),-1.0_pReal,1.0_pReal)),PI*.25_pReal] else eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ] end if From b0e06f180a18941bf47d00b86987825393c1ee15 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 Jul 2020 22:48:29 +0200 Subject: [PATCH 4/6] reading numpy manual ... https://numpy.org/doc/stable/reference/generated/numpy.isclose.html rtol=1.e-8, b=1 gives atol+btol=1.e-8+1e-9, simplified to 1e-8 --- python/damask/_rotation.py | 2 +- python/tests/test_Rotation.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 9bb2d300a..4c01bd430 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -805,7 +805,7 @@ class Rotation: """Rotation matrix to Bunge-Euler angles.""" with np.errstate(invalid='ignore',divide='ignore'): zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2) - eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-9), + eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,0.0), np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]), np.pi*0.5*(1-om[...,2,2:3]), np.zeros(om.shape[:-2]+(1,)), diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 21d0b5fae..46c10d8c2 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -150,7 +150,7 @@ def om2qu(a): def om2eu(om): """Rotation matrix to Bunge-Euler angles.""" - if not np.isclose(np.abs(om[2,2]),1.0,1.e-9): + if not np.isclose(np.abs(om[2,2]),1.0,0.0): zeta = 1.0/np.sqrt(1.0-om[2,2]**2) eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), np.arccos(om[2,2]), From 2b58b3df971887b92396eb28ee7a533c7d8063c9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 Jul 2020 06:55:40 +0200 Subject: [PATCH 5/6] closer to thoroughly tested python code but with extra safety measures to prevent invalid operations --- src/rotations.f90 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/rotations.f90 b/src/rotations.f90 index 5d865c2c3..07053545d 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -582,13 +582,11 @@ pure function om2eu(om) result(eu) real(pReal), dimension(3) :: eu real(pReal) :: zeta - if (dNeq(abs(om(3,3)),1.0_pReal,1.e-5_pReal)) then - zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal) + if (dNeq(abs(om(3,3)),1.0_pReal,1.e-8_pReal)) then + zeta = 1.0_pReal/sqrt(math_clip(1.0_pReal-om(3,3)**2.0_pReal,1e-64_pReal,1.0_pReal)) eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), & - acos(om(3,3)), & + acos(math_clip(om(3,3),-1.0_pReal,1.0_pReal)), & atan2(om(1,3)*zeta, om(2,3)*zeta)] - elseif(dNeq(abs(om(3,3)),1.0_pReal,1.e-8_pReal)) then - eu = [PI*.75_pReal,acos(math_clip(om(3,3),-1.0_pReal,1.0_pReal)),PI*.25_pReal] else eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ] end if From 133aa9111c86f3cca0aa911b337b23c4c73dc744 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 Jul 2020 07:18:32 +0200 Subject: [PATCH 6/6] small polishing --- src/marc/discretization_marc.f90 | 159 +++++++++++++++---------------- 1 file changed, 79 insertions(+), 80 deletions(-) diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index be16f5fc0..f0b77543b 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -21,24 +21,24 @@ module discretization_marc implicit none private - + type tCellNodeDefinition integer, dimension(:,:), allocatable :: parents integer, dimension(:,:), allocatable :: weights end type tCellNodeDefinition - + 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 public :: & discretization_marc_init - + contains !-------------------------------------------------------------------------------------------------- @@ -68,16 +68,16 @@ subroutine discretization_marc_init connectivity_elem real(pReal), dimension(:,:,:,:),allocatable :: & unscaledNormals - + class(tNode), pointer :: & num_commercialFEM - + write(6,'(/,a)') ' <<<+- discretization_marc init -+>>>'; flush(6) !--------------------------------------------------------------------------------- ! read debug parameters - debug_e = debug_root%get_asInt('element',defaultVal=1) - debug_i = debug_root%get_asInt('integrationpoint',defaultVal=1) + debug_e = debug_root%get_asInt('element',defaultVal=1) + debug_i = debug_root%get_asInt('integrationpoint',defaultVal=1) !-------------------------------------------------------------------------------- ! read numerics parameter and do sanity check @@ -90,10 +90,10 @@ subroutine discretization_marc_init if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP') - + FEsolving_execElem = [1,nElems] FEsolving_execIP = [1,elem%nIPs] - + allocate(cellNodeDefinition(elem%nNodes-1)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& @@ -108,7 +108,7 @@ subroutine discretization_marc_init call discretization_init(microstructureAt,homogenizationAt,& IP_reshaped,& node0_cell) - + call writeGeometry(elem,connectivity_elem,& reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),& node0_cell,IP_reshaped) @@ -120,7 +120,7 @@ subroutine discretization_marc_init call geometry_plastic_nonlocal_setIParea(norm2(unscaledNormals,1)) call geometry_plastic_nonlocal_setIPareaNormal(unscaledNormals/spread(norm2(unscaledNormals,1),1,3)) call geometry_plastic_nonlocal_results - + end subroutine discretization_marc_init @@ -139,7 +139,7 @@ subroutine writeGeometry(elem, & real(pReal), dimension(:,:), intent(in) :: & coordinates_nodes, & coordinates_points - + integer, dimension(:,:), allocatable :: & connectivity_temp real(pReal), dimension(:,:), allocatable :: & @@ -147,24 +147,24 @@ subroutine writeGeometry(elem, & call results_openJobFile call results_closeGroup(results_addGroup('geometry')) - + connectivity_temp = connectivity_elem call results_writeDataset('geometry',connectivity_temp,'T_e',& 'connectivity of the elements','-') - + connectivity_temp = connectivity_cell call results_writeDataset('geometry',connectivity_temp,'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', & 'initial coordinates of the nodes','m') - + coordinates_temp = coordinates_points call results_writeDataset('geometry',coordinates_temp,'x_p', & 'initial coordinates of the materialpoints','m') - + call results_closeJobFile end subroutine writeGeometry @@ -183,7 +183,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni integer, dimension(:), allocatable, intent(out) :: & microstructureAt, & homogenizationAt - + integer :: & fileFormatVersion, & hypoelasticTableStyle, & @@ -193,7 +193,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni 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 :: & nameElemSet integer, dimension(:,:), allocatable :: & @@ -209,8 +209,8 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni hypoelasticTableStyle,inputFile) call inputRead_NnodesAndElements(nNodes,nElems,& inputFile) - - + + call inputRead_mapElemSets(nameElemSet,mapElemSet,& inputFile) @@ -239,13 +239,13 @@ end subroutine inputRead !> @brief Figures out version of Marc input file format !-------------------------------------------------------------------------------------------------- subroutine inputRead_fileFormat(fileFormat,fileContent) - + integer, intent(out) :: fileFormat character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines - + integer, allocatable, dimension(:) :: chunkPos integer :: l - + do l = 1, size(fileContent) chunkPos = IO_stringPos(fileContent(l)) if(chunkPos(1) < 2) cycle @@ -262,13 +262,13 @@ end subroutine inputRead_fileFormat !> @brief Figures out table styles for initial cond and hypoelastic !-------------------------------------------------------------------------------------------------- subroutine inputRead_tableStyles(initialcond,hypoelastic,fileContent) - + integer, intent(out) :: initialcond, hypoelastic character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos integer :: l - + initialcond = 0 hypoelastic = 0 @@ -290,7 +290,7 @@ end subroutine inputRead_tableStyles !-------------------------------------------------------------------------------------------------- subroutine inputRead_matNumber(matNumber, & tableStyle,fileContent) - + integer, allocatable, dimension(:), intent(out) :: matNumber integer, intent(in) :: tableStyle character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines @@ -326,7 +326,7 @@ end subroutine inputRead_matNumber !-------------------------------------------------------------------------------------------------- subroutine inputRead_NnodesAndElements(nNodes,nElems,& fileContent) - + integer, intent(out) :: nNodes, nElems character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines @@ -355,7 +355,7 @@ end subroutine inputRead_NnodesAndElements !-------------------------------------------------------------------------------------------------- subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,& fileContent) - + integer, intent(out) :: nElemSets, maxNelemInSet character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines @@ -401,7 +401,7 @@ end subroutine inputRead_NelemSets !-------------------------------------------------------------------------------------------------- subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,& fileContent) - + character(len=pStringLen), dimension(:), allocatable, intent(out) :: nameElemSet integer, dimension(:,:), allocatable, intent(out) :: mapElemSet character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines @@ -434,7 +434,7 @@ end subroutine inputRead_mapElemSets !-------------------------------------------------------------------------------------------------- subroutine inputRead_mapElems(FEM2DAMASK, & nElems,nNodesPerElem,fileContent) - + integer, allocatable, dimension(:), intent(out) :: FEM2DAMASK integer, intent(in) :: nElems, & !< number of elements @@ -515,10 +515,10 @@ end subroutine inputRead_mapNodes subroutine inputRead_elemNodes(nodes, & nNode,fileContent) - real(pReal), allocatable, dimension(:,:), intent(out) :: nodes + real(pReal), allocatable, dimension(:,:), intent(out) :: nodes integer, intent(in) :: nNode character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines - + integer, allocatable, dimension(:) :: chunkPos integer :: i,j,m,l @@ -580,15 +580,15 @@ subroutine inputRead_elemType(elem, & endif enddo - contains + contains !-------------------------------------------------------------------------------------------------- !> @brief mapping of Marc element types to internal representation !-------------------------------------------------------------------------------------------------- integer function mapElemtype(what) - + character(len=*), intent(in) :: what - + select case (IO_lc(what)) case ( '6') mapElemtype = 1 ! Two-dimensional Plane Strain Triangle @@ -630,20 +630,20 @@ end subroutine inputRead_elemType !> @brief Stores node IDs !-------------------------------------------------------------------------------------------------- function inputRead_connectivityElem(nElem,nNodes,fileContent) - + integer, intent(in) :: & nElem, & nNodes !< number of nodes per element character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines - + integer, dimension(nNodes,nElem) :: & inputRead_connectivityElem - + integer, allocatable, dimension(:) :: chunkPos - + integer, dimension(1+nElem) :: contInts integer :: i,k,j,t,e,l,nNodesAlreadyRead - + do l = 1, size(fileContent) chunkPos = IO_stringPos(fileContent(l)) if(chunkPos(1) < 1) cycle @@ -694,7 +694,7 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos - + integer, dimension(1+nElem) :: contInts integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead,l,k,m @@ -727,42 +727,41 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza endif endif enddo - + end subroutine inputRead_microstructureAndHomogenization !-------------------------------------------------------------------------------------------------- -!> @brief Calculates cell node coordinates from element node coordinates +!> @brief Determine cell connectivity and definition of cell nodes !-------------------------------------------------------------------------------------------------- subroutine buildCells(connectivity_cell,cellNodeDefinition, & 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(tElement), intent(in) :: elem ! element definition integer, dimension(:,:), intent(in) :: connectivity_elem ! connectivity of the elements - + integer,dimension(:), allocatable :: candidates_local integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global - - integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID - + + integer :: e,n,c,p,s,i,m,j,& + nParentNodes,nCellNode,Nelem,candidateID + Nelem = size(connectivity_elem,2) !--------------------------------------------------------------------------------------------------- ! initialize global connectivity to negative local connectivity connectivity_cell = -spread(elem%cell,3,Nelem) ! local cell node ID - + !--------------------------------------------------------------------------------------------------- ! set connectivity of cell nodes that coincide with FE nodes (defined by 1 parent node) ! and renumber local (negative) to global (positive) node ID 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_cell(:,:,e) == -c) connectivity_cell(:,:,e) = connectivity_elem(c,e) endif realNode enddo enddo @@ -772,7 +771,7 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & !--------------------------------------------------------------------------------------------------- ! set connectivity of cell nodes that are defined by 2,...,nNodes real nodes do nParentNodes = 2, elem%nNodes - + ! get IDs of local cell nodes that are defined by the current number of parent nodes candidates_local = [integer::] do c = 1, elem%NcellNodes @@ -780,11 +779,11 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & candidates_local = [candidates_local,c] enddo s = size(candidates_local) - + if (allocated(candidates_global)) deallocate(candidates_global) allocate(candidates_global(nParentNodes*2+2,s*Nelem)) ! stores parent node ID + weight together with element ID and cellnode id (local) parentsAndWeights = reshape([(0, i = 1,2*nParentNodes)],[nParentNodes,2]) ! (re)allocate - + do e = 1, Nelem do i = 1, size(candidates_local) candidateID = (e-1)*size(candidates_local)+i ! including duplicates, runs to (Nelem*size(candidates_local)) @@ -792,18 +791,18 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & p = 0 do j = 1, size(elem%cellNodeParentNodeWeights(:,c)) if (elem%cellNodeParentNodeWeights(j,c) /= 0) then ! real node 'j' partly defines cell node 'c' - p = p + 1 + p = p + 1 parentsAndWeights(p,1:2) = [connectivity_elem(j,e),elem%cellNodeParentNodeWeights(j,c)] endif enddo ! store (and order) real node IDs and their weights together with the element number and local ID do p = 1, nParentNodes m = maxloc(parentsAndWeights(:,1),1) - + candidates_global(p, candidateID) = parentsAndWeights(m,1) candidates_global(p+nParentNodes, candidateID) = parentsAndWeights(m,2) candidates_global(nParentNodes*2+1:nParentNodes*2+2,candidateID) = [e,c] - + parentsAndWeights(m,1) = -huge(parentsAndWeights(m,1)) ! out of the competition enddo enddo @@ -811,7 +810,7 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & ! sort according to real node IDs + weight (from left to right) call math_sort(candidates_global,sortDim=1) ! sort according to first column - + do p = 2, nParentNodes*2 n = 1 do while(n <= size(candidates_local)*Nelem) @@ -826,11 +825,11 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & n = e+1 enddo enddo - + i = uniqueRows(candidates_global(1:2*nParentNodes,:)) allocate(cellNodeDefinition(nParentNodes-1)%parents(i,nParentNodes)) allocate(cellNodeDefinition(nParentNodes-1)%weights(i,nParentNodes)) - + i = 1 n = 1 do while(n <= size(candidates_local)*Nelem) @@ -846,7 +845,7 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & 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 end where - + j = j+1 enddo nCellNode = nCellNode + 1 @@ -863,7 +862,7 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & !> @brief count unique rows (same rows need to be stored consecutively) !------------------------------------------------------------------------------------------------ pure function uniqueRows(A) result(u) - + integer, dimension(:,:), intent(in) :: A !< array, rows need to be sorted integer :: & @@ -882,14 +881,14 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & u = u+1 r = r+d enddo - + end function uniqueRows end subroutine buildCells !-------------------------------------------------------------------------------------------------- -!> @brief Calculates cell node coordinates from element node coordinates +!> @brief Calculate cell node coordinates from element node coordinates !-------------------------------------------------------------------------------------------------- subroutine buildCellNodes(node_cell, & definition,node_elem) @@ -897,9 +896,9 @@ subroutine buildCellNodes(node_cell, & 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 - + integer :: i, j, k, n - + n = size(node_elem,2) node_cell(:,1:n) = node_elem !< initial nodes coincide with element nodes @@ -914,12 +913,12 @@ subroutine buildCellNodes(node_cell, & node_cell(:,n) = node_cell(:,n)/real(sum(definition(i)%weights(j,:)),pReal) enddo enddo - + end subroutine buildCellNodes !-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP coordinates as center of cell +!> @brief Calculate IP coordinates as center of cell !-------------------------------------------------------------------------------------------------- subroutine buildIPcoordinates(IPcoordinates, & connectivity_cell,node_cell) @@ -927,7 +926,7 @@ subroutine buildIPcoordinates(IPcoordinates, & 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 - + integer :: i, n do i = 1, size(connectivity_cell,2) @@ -938,7 +937,7 @@ subroutine buildIPcoordinates(IPcoordinates, & enddo IPcoordinates(:,i) = IPcoordinates(:,i)/real(size(connectivity_cell,1),pReal) enddo - + end subroutine buildIPcoordinates @@ -948,11 +947,11 @@ end subroutine buildIPcoordinates !> 2D cells assume an element depth of 1.0 !--------------------------------------------------------------------------------------------------- function IPvolume(elem,node,connectivity) - + 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(3) :: x0,x1,x2,x3,x4,x5,x6,x7 @@ -1011,19 +1010,19 @@ function IPareaNormal(elem,nElem,connectivity,node) 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 real(pReal), dimension (3,size(elem%cellFace,1)) :: nodePos integer :: e,i,f,n,m - + m = size(elem%cellFace,1) 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)) - + select case (elem%cellType) case (1,2) ! 2D 3 or 4 node IPareaNormal(1,f,i,e) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector @@ -1038,7 +1037,7 @@ function IPareaNormal(elem,nElem,connectivity,node) ! the sum has to be divided by two; this whole prcedure tries to compensate for ! probable non-planar cell surfaces IPareaNormal(1:3,f,i,e) = 0.0_pReal - do n = 1, m + do n = 1, m IPareaNormal(1:3,f,i,e) = IPareaNormal(1:3,f,i,e) & + math_cross(nodePos(1:3,mod(n+0,m)+1) - nodePos(1:3,n), & nodePos(1:3,mod(n+1,m)+1) - nodePos(1:3,n)) * 0.5_pReal