From 79e7deca5556e20a8d73c9acd1120e807bfedafc Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 5 Feb 2013 13:27:37 +0000 Subject: [PATCH] introduced dedicated calculation of interface normals in IP neighborhood. (2D elements did not work so far..!) fixed definitions of some (2D) interface planes to have all normals pointing outward. removed "int, private :: i" defined for spectral compilation. --- code/mesh.f90 | 191 ++++++++++++++++++++++++++++---------------------- 1 file changed, 108 insertions(+), 83 deletions(-) diff --git a/code/mesh.f90 b/code/mesh.f90 index 15fccacec..7b5a1ae5d 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -112,7 +112,6 @@ module mesh integer(pInt), public, protected :: & res1_red, & homog - integer(pInt), private :: i #endif ! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) @@ -160,6 +159,20 @@ module mesh 20 & ! element 21 (3D 20node 27ip) ],pInt) + integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_dimension = & !< dimension of element type + int([ & + 2, & ! element 6 (2D 3node 1ip) + 2, & ! element 125 (2D 6node 3ip) + 2, & ! element 11 (2D 4node 4ip) + 2, & ! element 27 (2D 8node 9ip) + 3, & ! element 134 (3D 4node 1ip) + 3, & ! element 127 (3D 10node 4ip) + 3, & ! element 136 (3D 6node 6ip) + 3, & ! element 117 (3D 8node 1ip) + 3, & ! element 7 (3D 8node 8ip) + 3 & ! element 21 (3D 20node 27ip) + ],pInt) + integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_Nnodes = & !< nodes in a specific type of element (how we use it) int([ & 3, & ! element 6 (2D 3node 1ip) @@ -386,7 +399,7 @@ contains !> @brief initializes the mesh by calling all necessary private routines the mesh module !! Order and routines strongly depend on type of solver !-------------------------------------------------------------------------------------------------- -subroutine mesh_init(ip,element) +subroutine mesh_init(ip,el) use DAMASK_interface use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) @@ -412,7 +425,8 @@ subroutine mesh_init(ip,element) implicit none integer(pInt), parameter :: fileUnit = 222_pInt - integer(pInt) :: e, element, ip + integer(pInt), intent(in) :: el, ip + integer(pInt) :: e write(6,*) write(6,*) '<<<+- mesh init -+>>>' @@ -449,12 +463,12 @@ subroutine mesh_init(ip,element) ! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and reso- ! lution-independent divergence if (divergence_correction == 1_pInt) then - do i = 1_pInt, 3_pInt - if (i/=minloc(geomdim,1) .and. i/=maxloc(geomdim,1)) scaledDim=geomdim/geomdim(i) + do e = 1_pInt, 3_pInt + if (e /= minloc(geomdim,1) .and. e /= maxloc(geomdim,1)) scaledDim = geomdim/geomdim(e) enddo elseif (divergence_correction == 2_pInt) then - do i = 1_pInt, 3_pInt - if (i/=minloc(geomdim/res,1) .and. i/=maxloc(geomdim/res,1)) scaledDim=geomdim/geomdim(i)*res(i) + do e = 1_pInt, 3_pInt + if (e /= minloc(geomdim/res,1) .and. e /= maxloc(geomdim/res,1)) scaledDim = geomdim/geomdim(e)*res(e) enddo else scaledDim = geomdim @@ -513,7 +527,7 @@ subroutine mesh_init(ip,element) parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive - FEsolving_execElem = [ 1_pInt,mesh_NcpElems] + FEsolving_execElem = [ 1_pInt,mesh_NcpElems ] if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP) allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt forall (e = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(FE_geomtype(mesh_element(2,e))) @@ -521,7 +535,7 @@ subroutine mesh_init(ip,element) if (allocated(calcMode)) deallocate(calcMode) allocate(calcMode(mesh_maxNips,mesh_NcpElems)) calcMode = .false. ! pretend to have collected what first call is asking (F = I) - calcMode(ip,mesh_FEasCP('elem',element)) = .true. ! first ip,el needs to be already pingponged to "calc" + calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" lastMode = .true. ! and its mode is already known... end subroutine mesh_init @@ -3434,24 +3448,35 @@ subroutine mesh_build_ipAreas real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: area - allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal + allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal allocate(mesh_ipAreaNormal(3_pInt,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal do e = 1_pInt,mesh_NcpElems ! loop over cpElems t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem - do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP - forall (n = 1_pInt:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles) ! start at each interface node and build valid triangles to cover interface - normal(:,j,n) = math_vectorproduct(nPos(:,1_pInt+mod(n+j-1_pInt,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors - nPos(:,1_pInt+mod(n+j-0_pInt,FE_NipFaceNodes)) - nPos(:,n)) - area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area - end forall - forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles, area(j,n) > 0.0_pReal) & - normal(1:3,j,n) = normal(1:3,j,n) / area(j,n) ! make myUnit normal - - mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles - mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2_pInt)/& ! average of all valid normals - real(count(area > 0.0_pReal),pReal) + do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP + select case (FE_dimension(t)) ! treat (1D), 2D, and 3D element types differently + + case (2) + forall (n = 1:2) nPos(1:3,n) = mesh_subNodeCoord(1:3,FE_subNodeOnIPFace(n,f,i,t),e) + mesh_ipAreaNormal(1,f,i,e) = nPos(2,2) - nPos(2,1) ! x_normal = y_connectingVector + mesh_ipAreaNormal(2,f,i,e) = -(nPos(1,2) - nPos(1,1)) ! y_normal = -x_connectingVector + mesh_ipArea(f,i,e) = sqrt(sum(mesh_ipAreaNormal(1:3,f,i,e)*mesh_ipAreaNormal(1:3,f,i,e))) ! and area + mesh_ipAreaNormal(1:3,f,i,e) = mesh_ipAreaNormal(1:3,f,i,e) / mesh_ipArea(f,i,e) ! have unit normal + + case (3) + forall (n = 1_pInt:FE_NipFaceNodes) nPos(1:3,n) = mesh_subNodeCoord(1:3,FE_subNodeOnIPFace(n,f,i,t),e) + forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles) ! start at each interface node and build valid triangles to cover interface + normal(1:3,j,n) = math_vectorproduct(nPos(1:3,1_pInt+mod(n+j-1_pInt,FE_NipFaceNodes)) - nPos(1:3,n), & ! calc their normal vectors + nPos(1:3,1_pInt+mod(n+j-0_pInt,FE_NipFaceNodes)) - nPos(1:3,n)) + area(j,n) = sqrt(sum(normal(1:3,j,n)*normal(1:3,j,n))) ! and area + end forall + forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles, area(j,n) > 0.0_pReal) & + normal(1:3,j,n) = normal(1:3,j,n) / area(j,n) ! have each normal of unit length + mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! vector product gives area of parallelograms instead of triangles + mesh_ipAreaNormal(1:3,f,i,e) = sum(sum(normal,3),2)/ & + max(1.0_pReal,real(count(area > 0.0_pReal),pReal)) ! average of all valid normals (not necessarily unit length---beware!!/fix?) + + end select enddo enddo enddo @@ -4364,16 +4389,16 @@ subroutine mesh_build_FEdata FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 27 (2D 8node 9ip) reshape(int([& 1, 1, 2, 0, 0, 0, 0, 0, 0, & - 1, 2, 2, 0, 0, 0, 0, 0, 0, & - 2, 2, 3, 0, 0, 0, 0, 0, 0, & - 2, 3, 3, 0, 0, 0, 0, 0, 0, & - 3, 3, 4, 0, 0, 0, 0, 0, 0, & - 3, 4, 4, 0, 0, 0, 0, 0, 0, & - 4, 4, 1, 0, 0, 0, 0, 0, 0, & - 4, 1, 1, 0, 0, 0, 0, 0, 0, & - 1, 1, 1, 1, 2, 2, 4, 4, 3, & - 2, 2, 2, 2, 1, 1, 3, 3, 4, & - 3, 3, 3, 3, 2, 2, 4, 4, 1, & + 1, 2, 2, 0, 0, 0, 0, 0, 0, & + 2, 2, 3, 0, 0, 0, 0, 0, 0, & + 2, 3, 3, 0, 0, 0, 0, 0, 0, & + 3, 3, 4, 0, 0, 0, 0, 0, 0, & + 3, 4, 4, 0, 0, 0, 0, 0, 0, & + 4, 4, 1, 0, 0, 0, 0, 0, 0, & + 4, 1, 1, 0, 0, 0, 0, 0, 0, & + 1, 1, 1, 1, 2, 2, 4, 4, 3, & + 2, 2, 2, 2, 1, 1, 3, 3, 4, & + 3, 3, 3, 3, 2, 2, 4, 4, 1, & 4, 4, 4, 4, 1, 1, 3, 3, 2 & ],pInt),[FE_Nips(me),FE_NsubNodes(me)]) @@ -4531,16 +4556,16 @@ subroutine mesh_build_FEdata FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 125 (2D 6node 3ip) reshape(int([& 4, 7, 7, 4 , & ! 1 - 1, 6, 6, 1 , & - 6, 7, 7, 6 , & + 6, 1, 1, 6 , & + 7, 6, 6, 7 , & 1, 4, 4, 1 , & 2, 5, 5, 2 , & ! 2 - 4, 7, 7, 4 , & + 7, 4, 4, 7 , & 5, 7, 7, 5 , & - 2, 4, 4, 2 , & - 5, 7, 7, 5 , & ! 3 + 4, 2, 2, 4 , & + 7, 5, 5, 7 , & ! 3 3, 6, 6, 3 , & - 3, 5, 5, 3 , & + 5, 3, 3, 5 , & 6, 7, 7, 6 & ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) @@ -4548,62 +4573,62 @@ subroutine mesh_build_FEdata FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 11 (2D 4node 4ip) reshape(int([& 5, 9, 9, 5 , & ! 1 - 1, 8, 8, 1 , & - 8, 9, 9, 8 , & + 8, 1, 1, 8 , & + 9, 8, 8, 9 , & 1, 5, 5, 1 , & 2, 6, 6, 2 , & ! 2 - 5, 9, 9, 5 , & + 9, 5, 5, 9 , & 6, 9, 9, 6 , & - 2, 5, 5, 2 , & - 3, 6, 6, 3 , & ! 3 + 5, 2, 2, 5 , & + 9, 7, 7, 9 , & ! 3 + 4, 8, 8, 4 , & + 7, 4, 4, 7 , & + 8, 9, 9, 8 , & + 6, 3, 3, 6 , & ! 4 7, 9, 9, 7 , & 3, 7, 7, 3 , & - 6, 9, 9, 6 , & - 7, 9, 9, 7 , & ! 4 - 4, 8, 8, 4 , & - 4, 7, 7, 4 , & - 8, 9, 9, 8 & + 9, 6, 6, 9 & ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) me = me + 1_pInt FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 27 (2D 8node 9ip) reshape(int([& - 9,17,17, 9 , & ! 1 - 1,16,16, 1 , & - 16,17,17,16 , & - 1, 9, 9, 1 , & - 10,18,18,10 , & ! 2 - 9,17,17, 9 , & - 17,18,18,17 , & - 9,10,10, 9 , & - 2,11,11, 2 , & ! 3 - 10,18,18,10 , & - 11,18,18,11 , & - 2,10,10, 2 , & - 17,20,20,17 , & ! 4 + 5,13,13, 5 , & ! 1 + 12, 1, 1,12 , & + 13,12,12,13 , & + 1, 5, 5, 1 , & + 6,14,14, 6 , & ! 2 + 13, 5, 5,13 , & + 14,13,13,14 , & + 5, 6, 6, 5 , & + 2, 7, 7, 2 , & ! 3 + 14, 6, 6,14 , & + 7,14,14, 7 , & + 6, 2, 2, 6 , & + 13,16,16,13 , & ! 4 + 11,12,12,11 , & + 16,11,11,16 , & + 12,13,13,12 , & + 14,15,15,14 , & ! 5 + 16,13,13,16 , & 15,16,16,15 , & - 15,20,20,15 , & - 16,17,17,16 , & - 18,19,19,18 , & ! 5 - 17,20,20,17 , & - 19,20,20,19 , & - 17,18,18,17 , & - 11,12,12,11 , & ! 6 - 18,19,19,18 , & - 12,19,19,12 , & - 11,18,18,11 , & - 14,20,20,14 , & ! 7 - 4,15,15, 4 , & - 4,14,14, 4 , & - 15,20,20,15 , & - 13,19,19,13 , & ! 8 - 14,20,20,14 , & 13,14,14,13 , & - 19,20,20,19 , & - 3,12,12, 3 , & ! 9 - 13,19,19,13 , & - 3,13,13, 3 , & - 12,19,19,12 & + 7, 8, 8, 7 , & ! 6 + 15,14,14,15 , & + 8,15,15, 8 , & + 14, 7, 7,14 , & + 16,10,10,16 , & ! 7 + 4,11,11, 4 , & + 10, 4, 4,10 , & + 11,16,16,11 , & + 15, 9, 9,15 , & ! 8 + 10,16,16,10 , & + 9,10,10, 9 , & + 16,15,15,16 , & + 8, 3, 3, 8 , & ! 9 + 9,15,15, 9 , & + 3, 9, 9, 3 , & + 15, 8, 8,15 & ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) me = me + 1_pInt