From dbe15f88f2191c39ba7fca15edbc58de1bfc2390 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 17:19:19 +0200 Subject: [PATCH] bugfix: forgot to read file first draft of nonlocal functions for area normal, area, and volume --- src/mesh_marc.f90 | 151 +++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 144 insertions(+), 7 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 24284f889..517865718 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -176,7 +176,6 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni microstructureAt, & homogenizationAt - integer :: & fileFormatVersion, & hypoelasticTableStyle, & @@ -194,6 +193,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni integer, dimension(:,:), allocatable :: & mapElemSet !< list of elements in elementSet + inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension)) call inputRead_fileFormat(fileFormatVersion, & inputFile) call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, & @@ -209,17 +209,17 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni call inputRead_mapElemSets(nameElemSet,mapElemSet,& FILEUNIT) - allocate (mesh_mapFEtoCPelem(2,nElems), source = 0) + allocate (mesh_mapFEtoCPelem(2,nElems), source=0) call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,FILEUNIT) - allocate (mesh_mapFEtoCPnode(2,Nnodes), source = 0) + allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0) call inputRead_mapNodes(inputFile) call inputRead_elemType(elem, & nElems,FILEUNIT) call inputRead_elemNodes(node0_elem, & Nnodes,inputFile) - + connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,FILEUNIT) call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & @@ -423,22 +423,30 @@ subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersio character(len=300) :: line, & tmp - integer, dimension (1+size(mesh_mapFEtoCPelem,2)) :: contInts + integer, dimension(:), allocatable :: contInts integer :: i,cpElem + allocate(contInts(size(mesh_mapFEtoCPelem,2)+1)) + write(6,*) 'hallo';flush(6) + write(6,*) nameElemSet;flush(6) + write(6,*) mapElemSet;flush(6) + write(6,*) 'map', size(mesh_mapFEtoCPelem,2);flush(6) cpElem = 0 contInts = 0 rewind(fileUnit) do read (fileUnit,'(A300)',END=620) line + write(6,*) trim(line);flush(6) chunkPos = IO_stringPos(line) if (fileFormatVersion < 13) then ! Marc 2016 or earlier if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then do i=1,3+TableStyle ! skip three (or four if new table style!) lines read (fileUnit,'(A300)') line + write(6,*) i;flush(6) enddo contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,& - mapElemSet,size(nameElemSet)) + mapElemSet,size(nameElemSet)) + exit endif else ! Marc2017 and later @@ -739,7 +747,6 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza 630 end subroutine inputRead_microstructureAndHomogenization - !-------------------------------------------------------------------------------------------------- !> @brief Calculates cell node coordinates from element node coordinates !-------------------------------------------------------------------------------------------------- @@ -926,6 +933,7 @@ subroutine buildCellNodes(node_cell, & end subroutine buildCellNodes + !-------------------------------------------------------------------------------------------------- !> @brief Calculates IP coordinates as center of cell !-------------------------------------------------------------------------------------------------- @@ -950,6 +958,135 @@ subroutine buildIPcoordinates(IPcoordinates, & end subroutine buildIPcoordinates +!--------------------------------------------------------------------------------------------------- +!> @brief Calculates IP volume. +!> @details The IP volume is calculated differently depending on the cell type. +!> 2D cells assume an element depth of one in order to calculate the volume. +!--------------------------------------------------------------------------------------------------- +function IPvolume(elem,node,connectivity) + + type(tElement) :: elem + real(pReal), dimension(:,:), intent(in) :: node + integer, dimension(:,:,:), intent(in) :: connectivity + real(pReal), dimension(elem%nIPs,size(connectivity,3)) :: IPvolume + + integer :: e,i + + do e = 1,size(connectivity,3) + do i = 1,elem%nIPs + + select case (elem%cellType) + case (1) ! 2D 3node + IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e))) + + case (2) ! 2D 4node + IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), & ! assume planar shape, division in two triangles suffices + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e))) & + + math_areaTriangle(node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(4,i,e)), & + node(1:3,connectivity(1,i,e))) + case (3) ! 3D 4node + IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(4,i,e))) + case (4) ! 3D 8node + IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(6,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(4,i,e)), & + node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(8,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(6,i,e)), & + node(1:3,connectivity(7,i,e)), & + node(1:3,connectivity(8,i,e)), & + node(1:3,connectivity(3,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(8,i,e)), & + node(1:3,connectivity(5,i,e)), & + node(1:3,connectivity(6,i,e)), & + node(1:3,connectivity(1,i,e))) ! first triangulation + IPvolume(i,e) = IPvolume(i,e) & + + math_volTetrahedron(node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(4,i,e)), & + node(1:3,connectivity(7,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(4,i,e)), & + node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(5,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(7,i,e)), & + node(1:3,connectivity(8,i,e)), & + node(1:3,connectivity(5,i,e)), & + node(1:3,connectivity(4,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(5,i,e)), & + node(1:3,connectivity(6,i,e)), & + node(1:3,connectivity(7,i,e)), & + node(1:3,connectivity(2,i,e))) ! second triangulation + IPvolume(i,e) = IPvolume(i,e) * 0.5_pReal + end select + enddo + enddo + +end function IPvolume + +!-------------------------------------------------------------------------------------------------- +!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_build_ipAreas(ipArea,ipAreaNormal,elem,nElem,connectivity,node) + + type(tElement) :: elem + integer :: nElem + integer, dimension(:,:,:), intent(in) :: connectivity + real(pReal), dimension(:,:), intent(in) :: node + real(pReal), dimension(elem%nIPneighbors,elem%nIPs,nElem), intent(out) :: ipArea + real(pReal), dimension(3,elem%nIPneighbors,elem%nIPs,nElem), intent(out) :: ipAreaNormal + + real(pReal), dimension (3,size(elem%cellFace,2)) :: nodePos, normals + real(pReal), dimension(3) :: normal + integer :: e,i,f,n,m + + m = size(elem%cellFace,1) + + do e = 1,nElem + do i = 1,elem%nIPs + do f = 1,size(elem%cellFace,2) !???? + + forall(n = 1: size(elem%cellface,1)) & + nodePos(1:3,n) = node(1:3,connectivity(elem%cellface(n,f),i,e)) + + select case (elem%cellType) + case (1,2) ! 2D 3 or 4 node + normal(1) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector + normal(2) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector + normal(3) = 0.0_pReal + case (3) ! 3D 4node + normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & + nodePos(1:3,3) - nodePos(1:3,1)) + case (4) ! 3D 8node + ! 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 + normals(1:3,n) = 0.5_pReal & + * math_cross(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) + end select + + ipArea(f,i,e) = norm2(normal) + ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal + + enddo + enddo + enddo + +end subroutine mesh_build_ipAreas + !-------------------------------------------------------------------------------------------------- !> @brief Gives the FE to CP ID mapping by binary search through lookup array !! valid questions (what) are 'elem', 'node'