From cac472c506b135f5bb2bae8ddd8f792e5ef5ab63 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 May 2019 00:13:43 +0200 Subject: [PATCH] cleaning more complex initialization will be used by FEM solvers only --- src/mesh_grid.f90 | 295 ++++------------------------------------------ 1 file changed, 23 insertions(+), 272 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index e752f5d41..a7da77bf0 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -35,8 +35,8 @@ module mesh mesh_unitlength !< physical length of one unit in mesh real(pReal), dimension(:,:), allocatable, private :: & - mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) - mesh_cellnode !< cell node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) + mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) + real(pReal), dimension(:,:), allocatable, public, protected :: & mesh_ipVolume, & !< volume associated with IP (initially!) @@ -53,53 +53,14 @@ module mesh logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) -integer(pInt), dimension(:,:), allocatable, private :: & - mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID integer(pInt),dimension(:,:,:), allocatable, private :: & mesh_cell !< cell connectivity for each element,ip/cell - integer(pInt), dimension(:,:,:), allocatable, private :: & - FE_cellface !< list of intra-cell cell node IDs that constitute the cell faces of a specific type of cell - - -! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) -! Hence, I suggest to prefix with "FE_" - integer(pInt), parameter, private :: & FE_Ngeomtypes = 10_pInt, & FE_Ncelltypes = 4_pInt, & - FE_maxNcellnodesPerCell = 8_pInt, & - FE_maxNcellfaces = 6_pInt, & - FE_maxNcellnodesPerCellface = 4_pInt - - - - integer(pInt), dimension(FE_Ncelltypes), parameter, private :: FE_NcellnodesPerCell = & !< number of cell nodes in a specific cell type - int([ & - 3, & ! (2D 3node) - 4, & ! (2D 4node) - 4, & ! (3D 4node) - 8 & ! (3D 8node) - ],pInt) - - integer(pInt), dimension(FE_Ncelltypes), parameter, private :: FE_NcellnodesPerCellface = & !< number of cell nodes per cell face in a specific cell type - int([& - 2, & ! (2D 3node) - 2, & ! (2D 4node) - 3, & ! (3D 4node) - 4 & ! (3D 8node) - ],pInt) - - - integer(pInt), dimension(FE_Ncelltypes), parameter, private :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type - int([& - 3, & ! (2D 3node) - 4, & ! (2D 4node) - 4, & ! (3D 4node) - 6 & ! (3D 8node) - ],pInt) - + FE_maxNcellnodesPerCell = 8_pInt integer(pInt), dimension(3), public, protected :: & grid !< (global) grid @@ -117,13 +78,10 @@ integer(pInt), dimension(:,:), allocatable, private :: & mesh_init private :: & - mesh_build_cellconnectivity, & mesh_build_ipAreas, & - mesh_build_FEdata, & mesh_spectral_build_nodes, & mesh_spectral_build_elements, & mesh_spectral_build_ipNeighborhood, & - mesh_build_cellnodes, & mesh_build_ipCoordinates type, public, extends(tMesh) :: tMesh_grid @@ -234,16 +192,14 @@ subroutine mesh_init(ip,el) if (myDebug) write(6,'(a)') ' Built elements'; flush(6) - call mesh_build_FEdata ! get properties of the different types of elements - call mesh_build_cellconnectivity - if (myDebug) write(6,'(a)') ' Built cell connectivity'; flush(6) - mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) + if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) mesh_ipCoordinates = mesh_build_ipCoordinates() if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) allocate(mesh_ipVolume(1,theMesh%nElems),source=product([geomSize(1:2),size3]/real([grid(1:2),grid3]))) if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) - call mesh_build_ipAreas + mesh_ipArea = mesh_build_ipAreas() + call mesh_build_ipAreas2 if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) call mesh_spectral_build_ipNeighborhood @@ -261,9 +217,6 @@ subroutine mesh_init(ip,el) !!!! COMPATIBILITY HACK !!!! -! for a homogeneous mesh, all elements have the same number of IPs and and cell nodes. -! hence, xxPerElem instead of maxXX -! better name theMesh%homogenizationAt = mesh_element(3,:) theMesh%microstructureAt = mesh_element(4,:) !!!!!!!!!!!!!!!!!!!!!!!! @@ -643,236 +596,34 @@ end function mesh_nodesAroundCentres !-------------------------------------------------------------------------------------------------- -!> @brief Split CP elements into cells. -!> @details Build a mapping between cells and the corresponding cell nodes ('mesh_cell'). -!> Cell nodes that are also matching nodes are unique in the list of cell nodes, -!> all others (currently) might be stored more than once. -!> Also allocates the 'mesh_node' array. +!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_cellconnectivity +pure function mesh_build_ipAreas() - implicit none - integer(pInt), dimension(:), allocatable :: & - matchingNode2cellnode - integer(pInt), dimension(:,:), allocatable :: & - cellnodeParent - integer(pInt), dimension(theMesh%elem%Ncellnodes) :: & - localCellnode2globalCellnode - integer(pInt) :: & - e,n,i, & - matchingNodeID, & - localCellnodeID - - integer(pInt), dimension(FE_Ngeomtypes), parameter :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry - int([ & - 3, & ! element 6 (2D 3node 1ip) - 3, & ! element 125 (2D 6node 3ip) - 4, & ! element 11 (2D 4node 4ip) - 4, & ! element 27 (2D 8node 9ip) - 4, & ! element 134 (3D 4node 1ip) - 4, & ! element 127 (3D 10node 4ip) - 6, & ! element 136 (3D 6node 6ip) - 8, & ! element 117 (3D 8node 1ip) - 8, & ! element 7 (3D 8node 8ip) - 8 & ! element 21 (3D 20node 27ip) - ],pInt) + real(pReal), dimension(6,1,theMesh%nElems) :: mesh_build_ipAreas - allocate(mesh_cell(FE_maxNcellnodesPerCell,theMesh%elem%nIPs,theMesh%nElems), source=0_pInt) - allocate(matchingNode2cellnode(theMesh%nNodes), source=0_pInt) - allocate(cellnodeParent(2_pInt,theMesh%elem%Ncellnodes*theMesh%nElems), source=0_pInt) - - mesh_Ncells = theMesh%nElems*theMesh%elem%nIPs -!-------------------------------------------------------------------------------------------------- -! Count cell nodes (including duplicates) and generate cell connectivity list - mesh_Ncellnodes = 0_pInt - - do e = 1_pInt,theMesh%nElems - localCellnode2globalCellnode = 0_pInt - do i = 1_pInt,theMesh%elem%nIPs - do n = 1_pInt,theMesh%elem%NcellnodesPerCell - localCellnodeID = theMesh%elem%cell(n,i) - if (localCellnodeID <= FE_NmatchingNodes(theMesh%elem%geomType)) then ! this cell node is a matching node - matchingNodeID = mesh_element(4_pInt+localCellnodeID,e) - if (matchingNode2cellnode(matchingNodeID) == 0_pInt) then ! if this matching node does not yet exist in the glbal cell node list ... - mesh_Ncellnodes = mesh_Ncellnodes + 1_pInt ! ... count it as cell node ... - matchingNode2cellnode(matchingNodeID) = mesh_Ncellnodes ! ... and remember its global ID - cellnodeParent(1_pInt,mesh_Ncellnodes) = e ! ... and where it belongs to - cellnodeParent(2_pInt,mesh_Ncellnodes) = localCellnodeID - endif - mesh_cell(n,i,e) = matchingNode2cellnode(matchingNodeID) - else ! this cell node is no matching node - if (localCellnode2globalCellnode(localCellnodeID) == 0_pInt) then ! if this local cell node does not yet exist in the global cell node list ... - mesh_Ncellnodes = mesh_Ncellnodes + 1_pInt ! ... count it as cell node ... - localCellnode2globalCellnode(localCellnodeID) = mesh_Ncellnodes ! ... and remember its global ID ... - cellnodeParent(1_pInt,mesh_Ncellnodes) = e ! ... and it belongs to - cellnodeParent(2_pInt,mesh_Ncellnodes) = localCellnodeID - endif - mesh_cell(n,i,e) = localCellnode2globalCellnode(localCellnodeID) - endif - enddo - enddo - enddo - - allocate(mesh_cellnodeParent(2_pInt,mesh_Ncellnodes)) - allocate(mesh_cellnode(3_pInt,mesh_Ncellnodes)) - - forall(n = 1_pInt:mesh_Ncellnodes) - mesh_cellnodeParent(1,n) = cellnodeParent(1,n) - mesh_cellnodeParent(2,n) = cellnodeParent(2,n) - endforall - -end subroutine mesh_build_cellconnectivity - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculate position of cellnodes from the given position of nodes -!> Build list of cellnodes' coordinates. -!> Cellnode coordinates are calculated from a weighted sum of node coordinates. -!-------------------------------------------------------------------------------------------------- -function mesh_build_cellnodes(nodes,Ncellnodes) - - implicit none - integer(pInt), intent(in) :: Ncellnodes !< requested number of cellnodes - real(pReal), dimension(3,mesh_Nnodes), intent(in) :: nodes - real(pReal), dimension(3,Ncellnodes) :: mesh_build_cellnodes - - integer(pInt) :: & - e,n,m, & - localCellnodeID - real(pReal), dimension(3) :: & - myCoords - - mesh_build_cellnodes = 0.0_pReal -!$OMP PARALLEL DO PRIVATE(e,localCellnodeID,myCoords) - do n = 1_pInt,Ncellnodes ! loop over cell nodes - e = mesh_cellnodeParent(1,n) - localCellnodeID = mesh_cellnodeParent(2,n) - myCoords = 0.0_pReal - do m = 1_pInt,theMesh%elem%nNodes - myCoords = myCoords + nodes(1:3,mesh_element(4_pInt+m,e)) & - * theMesh%elem%cellNodeParentNodeWeights(m,localCellnodeID) - enddo - mesh_build_cellnodes(1:3,n) = myCoords / sum(theMesh%elem%cellNodeParentNodeWeights(:,localCellnodeID)) - enddo -!$OMP END PARALLEL DO - -end function mesh_build_cellnodes - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates' -! Called by all solvers in mesh_init in order to initialize the ip coordinates. -! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus, -! so no need to use this subroutine anymore; Marc however only provides nodal displacements, -! so in this case the ip coordinates are always calculated on the basis of this subroutine. -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES, -! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME. -! HAS TO BE CHANGED IN A LATER VERSION. -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!-------------------------------------------------------------------------------------------------- + mesh_build_ipAreas(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3)) + mesh_build_ipAreas(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1)) + mesh_build_ipAreas(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2)) + +end function mesh_build_ipAreas !-------------------------------------------------------------------------------------------------- !> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipAreas - use math, only: & - math_crossproduct +subroutine mesh_build_ipAreas2 - implicit none - integer(pInt) :: e,t,g,c,i,f,n,m - real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals - real(pReal), dimension(3) :: normal + allocate(mesh_ipAreaNormal(3,6,1,theMesh%nElems), source=0.0_pReal) - allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(mesh_ipAreaNormal(3_pInt,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - - !$OMP PARALLEL DO PRIVATE(t,g,c,nodePos,normal,normals) - do e = 1_pInt,theMesh%nElems ! loop over cpElems - c = theMesh%elem%cellType - select case (c) - - case (4_pInt) - ! for this cell type we get the normal of the quadrilateral face as an average of - ! four normals of triangular subfaces; since the face consists only of two triangles, - ! the sum has to be divided by two; this whole prcedure tries to compensate for - ! probable non-planar cell surfaces - m = FE_NcellnodesPerCellface(c) - do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element - do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e)) - forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) & - normals(1:3,n) = 0.5_pReal & - * math_crossproduct(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & - nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n)) - normal = 0.5_pReal * sum(normals,2) - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) - enddo - enddo - - end select - enddo - !$OMP END PARALLEL DO - -end subroutine mesh_build_ipAreas - - -!-------------------------------------------------------------------------------------------------- -!> @brief get properties of different types of finite elements -!> @details assign globals: FE_nodesAtIP, FE_ipNeighbor, FE_subNodeOnIPFace -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_FEdata - - implicit none - integer(pInt) :: me - allocate(FE_cellface(FE_maxNcellnodesPerCellface,FE_maxNcellfaces,FE_Ncelltypes), source=0_pInt) - - - ! *** FE_cellface *** - me = 0_pInt - - me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 3node, VTK_TRIANGLE (5) - reshape(int([& - 2,3, & - 3,1, & - 1,2 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) - - me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 4node, VTK_QUAD (9) - reshape(int([& - 2,3, & - 4,1, & - 3,4, & - 1,2 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) - - me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 4node, VTK_TETRA (10) - reshape(int([& - 1,3,2, & - 1,2,4, & - 2,3,4, & - 1,4,3 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) - - me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 8node, VTK_HEXAHEDRON (12) - reshape(int([& - 2,3,7,6, & - 4,1,5,8, & - 3,4,8,7, & - 1,2,6,5, & - 5,6,7,8, & - 1,4,3,2 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) + mesh_ipAreaNormal(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_ipAreaNormal(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_ipAreaNormal(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_ipAreaNormal(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_ipAreaNormal(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,theMesh%nElems) + mesh_ipAreaNormal(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,theMesh%nElems) - -end subroutine mesh_build_FEdata +end subroutine mesh_build_ipAreas2 end module mesh