223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

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.
This commit is contained in:
Philip Eisenlohr 2013-02-05 13:27:37 +00:00
parent cc9eb685fe
commit 79e7deca55
1 changed files with 108 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

191
code/mesh.f90
View File

@ -112,7 +112,6 @@ module mesh
integer(pInt), public, protected :: & integer(pInt), public, protected :: &
res1_red, & res1_red, &
homog homog
integer(pInt), private :: i
#endif #endif
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) ! 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) 20 & ! element 21 (3D 20node 27ip)
],pInt) ],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) integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_Nnodes = & !< nodes in a specific type of element (how we use it)
int([ & int([ &
3, & ! element 6 (2D 3node 1ip) 3, & ! element 6 (2D 3node 1ip)
@ -386,7 +399,7 @@ contains
!> @brief initializes the mesh by calling all necessary private routines the mesh module !> @brief initializes the mesh by calling all necessary private routines the mesh module
!! Order and routines strongly depend on type of solver !! Order and routines strongly depend on type of solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_init(ip,element) subroutine mesh_init(ip,el)
use DAMASK_interface use DAMASK_interface
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) 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 implicit none
integer(pInt), parameter :: fileUnit = 222_pInt integer(pInt), parameter :: fileUnit = 222_pInt
integer(pInt) :: e, element, ip integer(pInt), intent(in) :: el, ip
integer(pInt) :: e
write(6,*) write(6,*)
write(6,*) '<<<+- mesh init -+>>>' 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- ! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and reso-
! lution-independent divergence ! lution-independent divergence
if (divergence_correction == 1_pInt) then if (divergence_correction == 1_pInt) then
do i = 1_pInt, 3_pInt do e = 1_pInt, 3_pInt
if (i/=minloc(geomdim,1) .and. i/=maxloc(geomdim,1)) scaledDim=geomdim/geomdim(i) if (e /= minloc(geomdim,1) .and. e /= maxloc(geomdim,1)) scaledDim = geomdim/geomdim(e)
enddo enddo
elseif (divergence_correction == 2_pInt) then elseif (divergence_correction == 2_pInt) then
do i = 1_pInt, 3_pInt do e = 1_pInt, 3_pInt
if (i/=minloc(geomdim/res,1) .and. i/=maxloc(geomdim/res,1)) scaledDim=geomdim/geomdim(i)*res(i) if (e /= minloc(geomdim/res,1) .and. e /= maxloc(geomdim/res,1)) scaledDim = geomdim/geomdim(e)*res(e)
enddo enddo
else else
scaledDim = geomdim 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 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) if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt 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))) 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) if (allocated(calcMode)) deallocate(calcMode)
allocate(calcMode(mesh_maxNips,mesh_NcpElems)) allocate(calcMode(mesh_maxNips,mesh_NcpElems))
calcMode = .false. ! pretend to have collected what first call is asking (F = I) 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... lastMode = .true. ! and its mode is already known...
end subroutine mesh_init end subroutine mesh_init
@ -3434,24 +3448,35 @@ subroutine mesh_build_ipAreas
real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal
real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: area 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 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 do e = 1_pInt,mesh_NcpElems ! loop over cpElems
t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType
do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem
do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP 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) select case (FE_dimension(t)) ! treat (1D), 2D, and 3D element types differently
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 case (2)
nPos(:,1_pInt+mod(n+j-0_pInt,FE_NipFaceNodes)) - nPos(:,n)) forall (n = 1:2) nPos(1:3,n) = mesh_subNodeCoord(1:3,FE_subNodeOnIPFace(n,f,i,t),e)
area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area mesh_ipAreaNormal(1,f,i,e) = nPos(2,2) - nPos(2,1) ! x_normal = y_connectingVector
end forall mesh_ipAreaNormal(2,f,i,e) = -(nPos(1,2) - nPos(1,1)) ! y_normal = -x_connectingVector
forall (n = 1_pInt:FE_NipFaceNodes, j = 1_pInt:Ntriangles, area(j,n) > 0.0_pReal) & mesh_ipArea(f,i,e) = sqrt(sum(mesh_ipAreaNormal(1:3,f,i,e)*mesh_ipAreaNormal(1:3,f,i,e))) ! and area
normal(1:3,j,n) = normal(1:3,j,n) / area(j,n) ! make myUnit normal mesh_ipAreaNormal(1:3,f,i,e) = mesh_ipAreaNormal(1:3,f,i,e) / mesh_ipArea(f,i,e) ! have unit normal
mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles case (3)
mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2_pInt)/& ! average of all valid normals forall (n = 1_pInt:FE_NipFaceNodes) nPos(1:3,n) = mesh_subNodeCoord(1:3,FE_subNodeOnIPFace(n,f,i,t),e)
real(count(area > 0.0_pReal),pReal) 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 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) FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 27 (2D 8node 9ip)
reshape(int([& reshape(int([&
1, 1, 2, 0, 0, 0, 0, 0, 0, & 1, 1, 2, 0, 0, 0, 0, 0, 0, &
1, 2, 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, 2, 3, 0, 0, 0, 0, 0, 0, &
2, 3, 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, 3, 4, 0, 0, 0, 0, 0, 0, &
3, 4, 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, 4, 1, 0, 0, 0, 0, 0, 0, &
4, 1, 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, & 1, 1, 1, 1, 2, 2, 4, 4, 3, &
2, 2, 2, 2, 1, 1, 3, 3, 4, & 2, 2, 2, 2, 1, 1, 3, 3, 4, &
3, 3, 3, 3, 2, 2, 4, 4, 1, & 3, 3, 3, 3, 2, 2, 4, 4, 1, &
4, 4, 4, 4, 1, 1, 3, 3, 2 & 4, 4, 4, 4, 1, 1, 3, 3, 2 &
],pInt),[FE_Nips(me),FE_NsubNodes(me)]) ],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) FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 125 (2D 6node 3ip)
reshape(int([& reshape(int([&
4, 7, 7, 4 , & ! 1 4, 7, 7, 4 , & ! 1
1, 6, 6, 1 , & 6, 1, 1, 6 , &
6, 7, 7, 6 , & 7, 6, 6, 7 , &
1, 4, 4, 1 , & 1, 4, 4, 1 , &
2, 5, 5, 2 , & ! 2 2, 5, 5, 2 , & ! 2
4, 7, 7, 4 , & 7, 4, 4, 7 , &
5, 7, 7, 5 , & 5, 7, 7, 5 , &
2, 4, 4, 2 , & 4, 2, 2, 4 , &
5, 7, 7, 5 , & ! 3 7, 5, 5, 7 , & ! 3
3, 6, 6, 3 , & 3, 6, 6, 3 , &
3, 5, 5, 3 , & 5, 3, 3, 5 , &
6, 7, 7, 6 & 6, 7, 7, 6 &
],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) ],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) FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 11 (2D 4node 4ip)
reshape(int([& reshape(int([&
5, 9, 9, 5 , & ! 1 5, 9, 9, 5 , & ! 1
1, 8, 8, 1 , & 8, 1, 1, 8 , &
8, 9, 9, 8 , & 9, 8, 8, 9 , &
1, 5, 5, 1 , & 1, 5, 5, 1 , &
2, 6, 6, 2 , & ! 2 2, 6, 6, 2 , & ! 2
5, 9, 9, 5 , & 9, 5, 5, 9 , &
6, 9, 9, 6 , & 6, 9, 9, 6 , &
2, 5, 5, 2 , & 5, 2, 2, 5 , &
3, 6, 6, 3 , & ! 3 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 , & 7, 9, 9, 7 , &
3, 7, 7, 3 , & 3, 7, 7, 3 , &
6, 9, 9, 6 , & 9, 6, 6, 9 &
7, 9, 9, 7 , & ! 4
4, 8, 8, 4 , &
4, 7, 7, 4 , &
8, 9, 9, 8 &
],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)])
me = me + 1_pInt me = me + 1_pInt
FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 27 (2D 8node 9ip) FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 27 (2D 8node 9ip)
reshape(int([& reshape(int([&
9,17,17, 9 , & ! 1 5,13,13, 5 , & ! 1
1,16,16, 1 , & 12, 1, 1,12 , &
16,17,17,16 , & 13,12,12,13 , &
1, 9, 9, 1 , & 1, 5, 5, 1 , &
10,18,18,10 , & ! 2 6,14,14, 6 , & ! 2
9,17,17, 9 , & 13, 5, 5,13 , &
17,18,18,17 , & 14,13,13,14 , &
9,10,10, 9 , & 5, 6, 6, 5 , &
2,11,11, 2 , & ! 3 2, 7, 7, 2 , & ! 3
10,18,18,10 , & 14, 6, 6,14 , &
11,18,18,11 , & 7,14,14, 7 , &
2,10,10, 2 , & 6, 2, 2, 6 , &
17,20,20,17 , & ! 4 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,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 , & 13,14,14,13 , &
19,20,20,19 , & 7, 8, 8, 7 , & ! 6
3,12,12, 3 , & ! 9 15,14,14,15 , &
13,19,19,13 , & 8,15,15, 8 , &
3,13,13, 3 , & 14, 7, 7,14 , &
12,19,19,12 & 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)]) ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)])
me = me + 1_pInt me = me + 1_pInt