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

- fixed bug in volume calculation of 2d elements 

- subroutine "mesh_build_ipVolumes" uses openmp parallelization
This commit is contained in:
Christoph Kords 2013-04-10 08:24:53 +00:00
parent 29758a3505
commit 338e338c9a
1 changed files with 23 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
code/mesh.f90
View File

@ -657,9 +657,12 @@ subroutine mesh_build_ipVolumes
math_areaTriangle math_areaTriangle
implicit none implicit none
integer(pInt) :: e,f,t,i,j,n,Ntriangles ! each interface is made up of this many triangles integer(pInt) :: e,f,t,i,j,n, &
integer(pInt), dimension(FE_maxNipFaceNodes) :: mySubNodes NmySubNodes, & !< number of subnodes already found for this (2d) ip cell
real(pReal), dimension(3,FE_maxNipFaceNodes) :: nPos ! coordinates of nodes on IP face Ntriangles !< each interface is made up of this many triangles
integer(pInt), dimension(FE_maxNipFaceNodes) :: mySubNodes !< list of subnodes that constitute this (2d) ip cell
real(pReal) :: myVolume
real(pReal), dimension(3,FE_maxNipFaceNodes) :: nPos !< coordinates of nodes on IP face
real(pReal), dimension(FE_maxNipFaceNodes-2_pInt,FE_maxNipFaceNodes) :: volume real(pReal), dimension(FE_maxNipFaceNodes-2_pInt,FE_maxNipFaceNodes) :: volume
if (.not. allocated(mesh_ipVolume)) then if (.not. allocated(mesh_ipVolume)) then
@ -667,36 +670,40 @@ subroutine mesh_build_ipVolumes
endif endif
mesh_ipVolume = 0.0_pReal mesh_ipVolume = 0.0_pReal
!$OMP PARALLEL DO PRIVATE(t,Ntriangles,nPos,volume,mySubNodes,NmySubNodes,myVolume)
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
select case (FE_dimension(t)) ! treat (1D), 2D, and 3D element types differently select case (FE_dimension(t)) ! treat (1D), 2D, and 3D element types differently
case (2) case (2_pInt)
Ntriangles = FE_NipNeighbors(t) - 2_pInt Ntriangles = FE_NipNeighbors(t) - 2_pInt
do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem
nPos = 0.0_pReal nPos = 0.0_pReal
volume = 0.0_pReal volume = 0.0_pReal
mySubNodes = 0_pInt mySubNodes = 0_pInt
do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG NmySubNodes = 0_pInt
faceLoop:do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG
do n = 1_pInt,FE_NipFaceNodes(t) do n = 1_pInt,FE_NipFaceNodes(t)
if (all(mySubNodes /= FE_subNodeOnIPFace(n,f,i,t))) then if (all(mySubNodes /= FE_subNodeOnIPFace(n,f,i,t))) then
mySubNodes(1) = mySubNodes(1) + 1_pInt NmySubNodes = NmySubNodes + 1_pInt
mySubNodes(mySubNodes(1)) = FE_subNodeOnIPFace(n,f,i,t) mySubNodes(NmySubNodes) = FE_subNodeOnIPFace(n,f,i,t)
nPos(1:3,f) = mesh_subNodeCoord(1:3,mySubNodes(mySubNodes(1)),e) nPos(1:3,NmySubNodes) = mesh_subNodeCoord(1:3,mySubNodes(NmySubNodes),e)
elseif (NmySubNodes == Fe_maxNipFaceNodes) then
exit faceLoop
endif endif
enddo enddo
enddo enddo faceLoop
forall (n = 1_pInt:FE_NipNeighbors(t), j = 1_pInt:Ntriangles) & ! start at each interface node and build valid triangles to cover interface forall (n = 1_pInt:FE_NipNeighbors(t), j = 1_pInt:Ntriangles) & ! start at each interface node and build valid triangles to cover interface
volume(j,n) = math_areaTriangle(nPos(1:3,n), & volume(j,n) = math_areaTriangle(nPos(1:3,n), &
nPos(1:3,1_pInt+mod(n-1_pInt +j ,FE_NipFaceNodes(t))), & nPos(1:3,1_pInt+mod(n+j-1_pInt,FE_NipNeighbors(t))), &
nPos(1:3,1_pInt+mod(n-1_pInt +j+1_pInt,FE_NipFaceNodes(t)))) nPos(1:3,1_pInt+mod(n+j ,FE_NipNeighbors(t)))) ! assuming element depth of 1
mesh_ipVolume(i,e) = sum(volume) / FE_NipFaceNodes(t) ! renormalize with interfaceNodeNum due to loop over them mesh_ipVolume(i,e) = sum(volume) / FE_NipFaceNodes(t) ! renormalize with interfaceNodeNum due to loop over them
enddo enddo
case (3) case (3_pInt)
Ntriangles = FE_NipFaceNodes(t) - 2_pInt Ntriangles = FE_NipFaceNodes(t) - 2_pInt
do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem
myVolume = 0.0_pReal
do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG
nPos = 0.0_pReal nPos = 0.0_pReal
volume = 0.0_pReal volume = 0.0_pReal
@ -707,12 +714,13 @@ subroutine mesh_build_ipVolumes
nPos(:,1_pInt+mod(n-1_pInt +j ,FE_NipFaceNodes(t))),& ! start at offset j nPos(:,1_pInt+mod(n-1_pInt +j ,FE_NipFaceNodes(t))),& ! start at offset j
nPos(:,1_pInt+mod(n-1_pInt +j+1_pInt,FE_NipFaceNodes(t))),& ! and take j's neighbor nPos(:,1_pInt+mod(n-1_pInt +j+1_pInt,FE_NipFaceNodes(t))),& ! and take j's neighbor
mesh_cellCenterCoordinates(i,e)) mesh_cellCenterCoordinates(i,e))
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface myVolume = myVolume + sum(volume) ! add contribution from this interface
enddo enddo
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes(t) ! renormalize with interfaceNodeNum due to loop over them mesh_ipVolume(i,e) = myVolume / FE_NipFaceNodes(t) ! renormalize with interfaceNodeNum due to loop over them
enddo enddo
endselect endselect
enddo enddo
!$OMP END PARALLEL DO
end subroutine mesh_build_ipVolumes end subroutine mesh_build_ipVolumes