Now able to have real periodicity for fluxes (fluxes leaving the model on one side enter on the other side).

To enable this feature one has to add the following somewhere in the marc input file:
   $mpie periodic x y z
      for having periodicity in all directions
   $mpie periodic z x
      for having periodicity in x and z direction
   etc.
Note that this only works for regular meshes!!!

code/mesh.f90

@@ -42,24 +42,26 @@
 
 !     order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype
 ! ---------------------------
- integer(pInt) mesh_Nelems,mesh_NcpElems,mesh_NelemSets,mesh_maxNelemInSet
+ integer(pInt) mesh_Nelems, mesh_NcpElems, mesh_NelemSets, mesh_maxNelemInSet
  integer(pInt) mesh_Nmaterials
- integer(pInt) mesh_Nnodes,mesh_maxNnodes,mesh_maxNips,mesh_maxNipNeighbors,mesh_maxNsharedElems,mesh_maxNsubNodes
+ integer(pInt) mesh_Nnodes, mesh_maxNnodes, mesh_maxNips, mesh_maxNipNeighbors, mesh_maxNsharedElems, mesh_maxNsubNodes
  integer(pInt), dimension(2) :: mesh_maxValStateVar = 0_pInt
  character(len=64), dimension(:),   allocatable :: mesh_nameElemSet, &    ! names of elementSet
                                                    mesh_nameMaterial, &   ! names of material in solid section
                                                    mesh_mapMaterial       ! name of elementSet for material
  integer(pInt), dimension(:,:),     allocatable :: mesh_mapElemSet        ! list of elements in elementSet
  integer(pInt), dimension(:,:),     allocatable, target :: mesh_mapFEtoCPelem, mesh_mapFEtoCPnode
- integer(pInt), dimension(:,:),     allocatable :: mesh_element, mesh_sharedElem
+ integer(pInt), dimension(:,:),     allocatable :: mesh_element, &
+                                                   mesh_sharedElem, &
+                                                   mesh_nodeTwins         ! node twins are surface nodes that lie exactly on opposite sides of the mesh (surfaces nodes with equal coordinate values in two dimensions)
  integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood
 
  real(pReal),   dimension(:,:,:),   allocatable :: mesh_subNodeCoord      ! coordinates of subnodes per element
- real(pReal),   dimension(:,:),     allocatable :: mesh_ipVolume          ! volume associated with IP
+ real(pReal),   dimension(:,:),     allocatable :: mesh_node, & 
+                                                   mesh_ipVolume          ! volume associated with IP
  real(pReal),   dimension(:,:,:),   allocatable :: mesh_ipArea, &         ! area of interface to neighboring IP
                                                    mesh_ipCenterOfGravity ! center of gravity of IP
  real(pReal),   dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal      ! area normal of interface to neighboring IP
- real(pReal),                       allocatable :: mesh_node (:,:)
  
  integer(pInt), dimension(:,:,:),   allocatable :: FE_nodesAtIP
  integer(pInt), dimension(:,:,:),   allocatable :: FE_ipNeighbor
@@ -67,6 +69,7 @@
  integer(pInt), dimension(:,:,:,:), allocatable :: FE_subNodeOnIPFace
  
  logical :: noPart                                                        ! for cases where the ABAQUS input file does not use part/assembly information
+ logical, dimension(3) :: mesh_periodicSurface                            ! flag indicating periodic outer surfaces (used for fluxes)
 
  integer(pInt) :: hypoelasticTableStyle
  integer(pInt) :: initialcondTableStyle
@@ -266,6 +269,7 @@
                          call mesh_marc_build_nodes(fileUnit)
                          call mesh_marc_count_cpSizes(fileunit)
                          call mesh_marc_build_elements(fileUnit)
+                         call mesh_marc_get_mpieOptions(fileUnit)
      case ('Abaqus')
                          noPart = IO_abaqus_hasNoPart(fileUnit)
                          call mesh_abaqus_count_nodesAndElements(fileUnit)
@@ -282,11 +286,12 @@
    end select
    close (fileUnit)
    
-   call mesh_build_sharedElems()
-   call mesh_build_ipNeighborhood()
    call mesh_build_subNodeCoords()
    call mesh_build_ipVolumes()
    call mesh_build_ipAreas()
+   call mesh_build_nodeTwins()
+   call mesh_build_sharedElems()
+   call mesh_build_ipNeighborhood()
    call mesh_tell_statistics()
 
    parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems))      ! plus potential killer from non-local constitutive
@@ -413,56 +418,97 @@
 !***********************************************************
 ! find face-matching element of same type
 !***********************************************************
- function mesh_faceMatch(elem,face)
+subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace)
 
- use prec, only: pInt
+use prec, only: pInt
- implicit none
+implicit none
 
- integer(pInt) face,elem
+!*** output variables
- integer(pInt), dimension(2) :: mesh_faceMatch                                     ! matching element's ID and corresponding face ID 
+integer(pInt), intent(out) ::     matchingElem, &       ! matching CP element ID
- integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: myFaceNodes ! global node ids on my face
+                                  matchingFace          ! matching FE face ID 
- integer(pInt) minN,NsharedElems, &
-               t, lonelyNode, &
-               candidateType,candidateElem, &
-               i,f,n
 
- minN = mesh_maxNsharedElems+1         ! init to worst case
+!*** input variables
- mesh_faceMatch = 0_pInt               ! intialize to "no match found"
+integer(pInt), intent(in) ::      face, &               ! FE face ID
- t = mesh_element(2,elem)              ! figure elemType
+                                  elem                  ! FE elem ID
 
- do n = 1,FE_NfaceNodes(face,t)        ! loop over nodes on face
+!*** local variables
-   myFaceNodes(n) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(n,face,t),elem)) ! CP id of face node
+integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: &
+                                  myFaceNodes           ! global node ids on my face
+integer(pInt)                     myType, &
+                                  candidateType, &
+                                  candidateElem, &
+                                  candidateFace, &
+                                  candidateFaceNode, &
+                                  minNsharedElems, &
+                                  NsharedElems, &
+                                  lonelyNode, &
+                                  i, &
+                                  n, &
+                                  dir                   ! periodicity direction
+integer(pInt), dimension(:), allocatable :: element_seen
+logical checkTwins
+
+
+minNsharedElems = mesh_maxNsharedElems + 1_pInt                   ! init to worst case
+matchingFace = 0_pInt
+matchingElem = 0_pInt                                             ! intialize to "no match found"
+myType = mesh_element(2,elem)                                     ! figure elemType
+
+do n = 1,FE_NfaceNodes(face,myType)                               ! loop over nodes on face
+  myFaceNodes(n) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(n,face,myType),elem)) ! CP id of face node
   NsharedElems = mesh_sharedElem(1,myFaceNodes(n))                ! figure # shared elements for this node
-   if (NsharedElems < minN) then
+  if (NsharedElems < minNsharedElems) then
-     minN = NsharedElems               ! remember min # shared elems
+    minNsharedElems = NsharedElems                                ! remember min # shared elems
     lonelyNode = n                                                ! remember most lonely node
   endif
- enddo
+enddo
 
+allocate(element_seen(minNsharedElems))
+element_seen = 0_pInt
 
-candidate: do i = 1,minN               ! iterate over lonelyNode's shared elements
+checkCandidate: do i = 1,minNsharedElems                            ! iterate over lonelyNode's shared elements
   candidateElem = mesh_sharedElem(1+i,myFaceNodes(lonelyNode))      ! present candidate elem
-   if (candidateElem == elem) cycle candidate         ! my own element ?
+  if (all(element_seen /= candidateElem)) then                      ! element seen for the first time?
+    element_seen(i) = candidateElem
     candidateType = mesh_element(2,candidateElem)                   ! figure elemType of candidate
-
+checkCandidateFace: do candidateFace = 1,FE_maxNipNeighbors         ! check each face of candidate
-candidateFace: do f = 1,FE_maxNipNeighbors            ! check each face of candidate
+      if (FE_NfaceNodes(candidateFace,candidateType) /= FE_NfaceNodes(face,myType) & ! incompatible face
-        if (FE_NfaceNodes(f,candidateType) /= FE_NfaceNodes(face,t)) &
+          .or. (candidateElem == elem .and. candidateFace == face)) then  ! this is my face
-          cycle candidateFace                         ! incompatible face
+        cycle checkCandidateFace
-        do n = 1,FE_NfaceNodes(f,candidateType)       ! loop through nodes on face
+      endif
-          if (all(myFaceNodes /= &
+      checkTwins = .false.
-              mesh_FEasCP('node', &
+      do n = 1,FE_NfaceNodes(candidateFace,candidateType)           ! loop through nodes on face
-                mesh_element(4+FE_nodeOnFace(n,f,candidateType),candidateElem)))) &
+        candidateFaceNode = mesh_FEasCP('node', mesh_element(4+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem))
-            cycle candidateFace  ! other face node not one of my face nodes
+        if (all(myFaceNodes /= candidateFaceNode)) then             ! candidate node does not match any of my face nodes
+          checkTwins = .true.                                       ! perhaps the twin nodes do match
+          exit
+        endif
       enddo
-        mesh_faceMatch(1) = f
+      if(checkTwins) then
-        mesh_faceMatch(2) = candidateElem
+checkCandidateFaceTwins: do dir = 1,3
-        exit candidate                                ! found my matching candidate
+          do n = 1,FE_NfaceNodes(candidateFace,candidateType)       ! loop through nodes on face
-      enddo candidateFace
+            candidateFaceNode = mesh_FEasCP('node', mesh_element(4+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem))
- enddo candidate
+            if (all(myFaceNodes /= mesh_nodeTwins(dir,candidateFaceNode))) then ! node twin does not match either
+              if (dir == 3) then
+                cycle checkCandidateFace
+              else
+                cycle checkCandidateFaceTwins                       ! try twins in next dimension
+              endif
+            endif
+          enddo
+          exit checkCandidateFaceTwins
+        enddo checkCandidateFaceTwins
+      endif
+      matchingFace = candidateFace
+      matchingElem = candidateElem
+      exit checkCandidate                                           ! found my matching candidate
+    enddo checkCandidateFace
+  endif
+enddo checkCandidate
 
- return
+deallocate(element_seen)
 
- endfunction
+endsubroutine
 
  
 !********************************************************************
@@ -562,10 +608,10 @@ candidateFace: do f = 1,FE_maxNipNeighbors            ! check each face of candi
     7,8, 0,0,  &
     7,0, 0,0   &
     /),(/FE_maxNnodesAtIP(7),FE_Nips(7)/))
-!! FE_nodesAtIP(:,:FE_Nips(8),8) = &  ! element 117 WHY IS THIS NOT DEFINED (AS FOR 134)??
+! FE_nodesAtIP(:,:FE_Nips(8),8) = &  ! element 117 (c3d8r --> single IP per element, so no need for this mapping)
-!!    reshape((/&
+!    reshape((/&
-!!    1,2,3,4,5,6,7,8   &
+!    1,2,3,4,5,6,7,8   &
-!!    /),(/FE_maxNnodesAtIP(8),FE_Nips(8)/))
+!    /),(/FE_maxNnodesAtIP(8),FE_Nips(8)/))
  FE_nodesAtIP(:,:FE_Nips(9),9) = &  ! element 57 (c3d20r == c3d8 --> copy of 7)
     reshape((/&
     1,  &
@@ -1257,6 +1303,55 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = &  ! element 117
  endsubroutine
 
 
+!********************************************************************
+! get any additional mpie options from input file (Marc only)
+!
+! mesh_periodicSurface
+!********************************************************************
+subroutine mesh_marc_get_mpieOptions(unit)
+
+use prec, only: pInt
+use IO
+implicit none
+
+integer(pInt), intent(in) :: unit
+
+integer(pInt), parameter :: maxNchunks = 5
+integer(pInt), dimension (1+2*maxNchunks) :: pos
+integer(pInt) chunk
+character(len=300) line
+
+mesh_periodicSurface = (/.false., .false., .false./)
+
+610 FORMAT(A300)
+
+rewind(unit)
+do 
+  read (unit,610,END=620) line
+  pos = IO_stringPos(line,maxNchunks)
+
+  if (IO_lc(IO_stringValue(line,pos,1)) == '$mpie') then              ! found keyword for user defined input
+    if (IO_lc(IO_stringValue(line,pos,2)) == 'periodic' &             ! found keyword 'periodic'
+        .and. pos(1) > 3) then                                        ! and there is at least one more chunk to read
+      do chunk = 2,pos(1)                                             ! loop through chunks (skipping the keyword)
+        select case(IO_stringValue(line,pos,chunk))                   ! chunk matches keyvalues x,y or z?
+          case('x')
+            mesh_periodicSurface(1) = .true.
+          case('y')
+            mesh_periodicSurface(2) = .true.
+          case('z')
+            mesh_periodicSurface(3) = .true.
+        end select
+      enddo
+    endif
+  endif
+enddo
+
+620 return
+
+endsubroutine
+
+ 
 !********************************************************************
 ! count overall number of nodes and elements in mesh
 !
@@ -2686,53 +2781,70 @@ subroutine mesh_marc_count_cpSizes (unit)
 ! _maxNsharedElems
 ! _sharedElem
 !********************************************************************
- subroutine mesh_build_sharedElems ()
+subroutine mesh_build_sharedElems()
 
- use prec, only: pInt
+use prec, only: pInt
- use IO
+implicit none
- implicit none
 
- integer(pint) e,t,n,j
+integer(pint)   e, &      ! element index
- integer(pInt), dimension (mesh_Nnodes) :: node_count
+                t, &      ! element type
- integer(pInt), dimension (:), allocatable :: node_seen
+                node, &   ! CP node index
+                j, &      ! node index per element 
+                dim, &    ! dimension index 
+                nodeTwin  ! node twin in the specified dimension
+integer(pInt), dimension (mesh_Nnodes) :: node_count
+integer(pInt), dimension (:), allocatable :: node_seen
 
- allocate(node_seen(maxval(FE_Nnodes)))
+allocate(node_seen(maxval(FE_Nnodes)))
- node_count = 0_pInt
 
- do e = 1,mesh_NcpElems
+node_count = 0_pInt
+
+do e = 1,mesh_NcpElems
+  t = mesh_element(2,e)
 
-   t = mesh_element(2,e)                                                  ! my type
   node_seen = 0_pInt                                                      ! reset node duplicates
   do j = 1,FE_Nnodes(t)                                                   ! check each node of element
-     n = mesh_FEasCP('node',mesh_element(4+j,e))
+    node = mesh_FEasCP('node',mesh_element(4+j,e))
-     if (all(node_seen /= n)) node_count(n) = node_count(n) + 1_pInt      ! if FE node not yet encountered -> count it
+    if (all(node_seen /= node)) then
-     node_seen(j) = n                                                     ! remember this node to be counted already
+      node_count(node) = node_count(node) + 1_pInt                        ! if FE node not yet encountered -> count it
+      do dim = 1,3                                                        ! check in each dimension...
+        nodeTwin = mesh_nodeTwins(dim,node)
+        if (nodeTwin > 0) &                                               ! if i am a twin of some node...
+          node_count(nodeTwin) = node_count(nodeTwin) + 1_pInt            ! -> count me again for the twin node
       enddo
-
- enddo
-
- mesh_maxNsharedElems = maxval(node_count)                                ! most shared node
- 
- allocate ( mesh_sharedElem( 1+mesh_maxNsharedElems,mesh_Nnodes) )
- mesh_sharedElem = 0_pInt
-
- do e = 1,mesh_NcpElems
-   node_seen = 0_pInt
-   do j = 1,FE_Nnodes(mesh_element(2,e))
-     n = mesh_FEasCP('node',mesh_element(4+j,e))
-     if (all(node_seen /= n)) then
-       mesh_sharedElem(1,n) = mesh_sharedElem(1,n) + 1                    ! count for each node the connected elements
-       mesh_sharedElem(1+mesh_sharedElem(1,n),n) = e                      ! store the respective element id
     endif
-     node_seen(j) = n
+    node_seen(j) = node                                                   ! remember this node to be counted already
   enddo
+enddo
+
+mesh_maxNsharedElems = maxval(node_count)                                 ! most shared node
+
+allocate(mesh_sharedElem(1+mesh_maxNsharedElems,mesh_Nnodes))
+mesh_sharedElem = 0_pInt
+
+do e = 1,mesh_NcpElems
+  t = mesh_element(2,e)
+  node_seen = 0_pInt
+  do j = 1,FE_Nnodes(t)
+    node = mesh_FEasCP('node',mesh_element(4+j,e))
+    if (all(node_seen /= node)) then
+      mesh_sharedElem(1,node) = mesh_sharedElem(1,node) + 1_pInt          ! count for each node the connected elements
+      mesh_sharedElem(mesh_sharedElem(1,node)+1,node) = e                 ! store the respective element id
+      do dim = 1,3                                                        ! check in each dimension...
+        nodeTwin = mesh_nodeTwins(dim,node)
+        if (nodeTwin > 0) then                                            ! if i am a twin of some node...
+          mesh_sharedElem(1,nodeTwin) = mesh_sharedElem(1,nodeTwin) + 1_pInt ! ...count me again for the twin
+          mesh_sharedElem(mesh_sharedElem(1,nodeTwin)+1,nodeTwin) = e     ! store the respective element id
+        endif
       enddo
+    endif
+    node_seen(j) = node
+  enddo
+enddo
 
- deallocate (node_seen)
+deallocate(node_seen)
 
- return
+endsubroutine
-
- endsubroutine
 
 
 
@@ -2742,74 +2854,152 @@ subroutine mesh_marc_count_cpSizes (unit)
 ! allocate globals
 ! _ipNeighborhood
 !***********************************************************
- subroutine mesh_build_ipNeighborhood()
+subroutine mesh_build_ipNeighborhood()
 
- use prec, only: pInt
+use prec, only: pInt
- implicit none
+implicit none
 
- integer(pInt) e,t,i,j,k,l,m,n,a,anchor, neighborType
+integer(pInt)                   myElem, &               ! my CP element index
- integer(pInt) neighbor,neighboringElem,neighboringIP
+                                myIP, &
- integer(pInt), dimension(2) :: matchingElem                              ! face and id of matching elem
+                                myType, &               ! my element type
- integer(pInt), dimension(FE_maxmaxNnodesAtIP) :: linkedNodes, &
+                                myFace, &
+                                neighbor, &             ! neighor index
+                                neighboringIPkey, &     ! positive integer indicating the neighboring IP (for intra-element) and negative integer indicating the face towards neighbor (for neighboring element) 
+                                candidateIP, &
+                                neighboringType, &      ! element type of neighbor
+                                NlinkedNodes, &         ! number of linked nodes
+                                twin_of_linkedNode, &   ! node twin of a specific linkedNode
+                                NmatchingNodes, &       ! number of matching nodes
+                                dir, &                  ! direction of periodicity
+                                matchingElem, &         ! CP elem number of matching element
+                                matchingFace, &         ! face ID of matching element
+                                k, a, anchor
+integer(pInt), dimension(FE_maxmaxNnodesAtIP) :: &
+                                linkedNodes, &
                                 matchingNodes
+logical checkTwins
 
+allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems))
+mesh_ipNeighborhood = 0_pInt
+
+linkedNodes = 0_pInt
+
+do myElem = 1,mesh_NcpElems                                                         ! loop over cpElems
+  myType = mesh_element(2,myElem)                                                   ! get elemType
+  do myIP = 1,FE_Nips(myType)                                                       ! loop over IPs of elem
+    do neighbor = 1,FE_NipNeighbors(myType)                                         ! loop over neighbors of IP
+      neighboringIPkey = FE_ipNeighbor(neighbor,myIP,myType)
+      
+      !*** if the key is positive, the neighbor is inside the element
+      !*** that means, we have already found our neighboring IP
+      
+      if (neighboringIPkey > 0) then
+        mesh_ipNeighborhood(1,neighbor,myIP,myElem) = myElem
+        mesh_ipNeighborhood(2,neighbor,myIP,myElem) = neighboringIPkey
+
+
+      !*** if the key is negative, the neighbor resides in a neighboring element
+      !*** that means, we have to look through the face indicated by the key and see which element is behind that face
+      
+      elseif (neighboringIPkey < 0) then                                            ! neighboring element's IP
+        myFace = -neighboringIPkey
+        call mesh_faceMatch(myElem, myFace, matchingElem, matchingFace)             ! get face and CP elem id of face match
+        if (matchingElem > 0_pInt) then                                             ! found match?
+          neighboringType = mesh_element(2,matchingElem)
+
+          !*** trivial solution if neighbor has only one IP
+          
+          if (FE_Nips(neighboringType) == 1_pInt) then            
+            mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem
+            mesh_ipNeighborhood(2,neighbor,myIP,myElem) = 1_pInt
+            cycle
+          endif
+
+          !*** find those nodes which build the link to the neighbor
+          
+          NlinkedNodes = 0_pInt
           linkedNodes = 0_pInt
-
+          do a = 1,FE_maxNnodesAtIP(myType)                                         ! figure my anchor nodes on connecting face
- allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipNeighborhood = 0_pInt
+            anchor = FE_nodesAtIP(a,myIP,myType)
- 
+            if (anchor /= 0_pInt) then                                              ! valid anchor node
- do e = 1,mesh_NcpElems                                                   ! loop over cpElems
+              if (any(FE_nodeOnFace(:,myFace,myType) == anchor)) then               ! ip anchor sits on face?
-   t = mesh_element(2,e)                                                  ! get elemType
+                NlinkedNodes = NlinkedNodes + 1_pInt
-   do i = 1,FE_Nips(t)                                                    ! loop over IPs of elem
+                linkedNodes(NlinkedNodes) = mesh_element(4+anchor,myElem)           ! FE id of anchor node
-     do n = 1,FE_NipNeighbors(t)                                          ! loop over neighbors of IP
+              else                                                                  ! something went wrong with the linkage, since not all anchors sit on my face
-       neighbor = FE_ipNeighbor(n,i,t)
+                NlinkedNodes = 0_pInt
-       if (neighbor > 0) then                                             ! intra-element IP
+                linkedNodes = 0_pInt
-         neighboringElem = e
+                exit
-         neighboringIP   = neighbor
+              endif
-       else                                                               ! neighboring element's IP
-         neighboringElem = 0_pInt
-         neighboringIP = 0_pInt
-         matchingElem = mesh_faceMatch(e,-neighbor)                       ! get face and CP elem id of face match
-         if (matchingElem(2) > 0_pInt) then                               ! found match?
-           neighborType = mesh_element(2,matchingElem(2))
-           l = 0_pInt
-           do a = 1,FE_maxNnodesAtIP(t)                                   ! figure my anchor nodes on connecting face
-             anchor = FE_nodesAtIP(a,i,t)
-             if (any(FE_nodeOnFace(:,-neighbor,t) == anchor)) then        ! ip anchor sits on face?
-               l = l+1_pInt
-               linkedNodes(l) = mesh_element(4+anchor,e)                  ! FE id of anchor node
             endif
           enddo
 
-  neighborIP: do j = 1,FE_Nips(neighborType)                              ! loop over neighboring ips
+          !*** loop through the ips of my neighbor
-             m = 0_pInt
+          !*** and try to find an ip with matching nodes
+          !*** also try to match with node twins
+          
+checkCandidateIP: do candidateIP = 1,FE_Nips(neighboringType)
+            NmatchingNodes = 0_pInt
             matchingNodes = 0_pInt
-             do a = 1,FE_maxNnodesAtIP(neighborType)                      ! check each anchor node of that ip
+            do a = 1,FE_maxNnodesAtIP(neighboringType)                              ! check each anchor node of that ip
-               anchor = FE_nodesAtIP(a,j,neighborType)
+              anchor = FE_nodesAtIP(a,candidateIP,neighboringType)
-               if (    anchor /= 0_pInt .and. &    ! valid anchor node
+              if (anchor /= 0_pInt) then                                            ! valid anchor node
-                   any(FE_nodeOnFace(:,matchingElem(1),neighborType) == anchor) ) then
+                if (any(FE_nodeOnFace(:,matchingFace,neighboringType) == anchor)) then  ! sits on matching face?
-                 m = m+1_pInt
+                  NmatchingNodes = NmatchingNodes + 1_pInt
-                 matchingNodes(m) = mesh_element(4+anchor,matchingElem(2)) ! FE id of neighbor's anchor node
+                  matchingNodes(NmatchingNodes) = mesh_element(4+anchor,matchingElem)   ! FE id of neighbor's anchor node
+                else                                                                ! no matching, because not all nodes sit on the matching face
+                  NmatchingNodes = 0_pInt
+                  matchingNodes = 0_pInt
+                  exit
+                endif
               endif
             enddo
-             if (m /= l) cycle neighborIP                                 ! this ip has wrong count of anchors on face
+
-             do a = 1,l
+            if (NmatchingNodes /= NlinkedNodes) &                                   ! this ip has wrong count of anchors on face
-               if (all(matchingNodes /= linkedNodes(a))) cycle neighborIP
+              cycle checkCandidateIP
+            
+            !*** check "normal" nodes whether they match or not
+            
+            checkTwins = .false.
+            do a = 1,NlinkedNodes
+              if (all(matchingNodes /= linkedNodes(a))) then                        ! this linkedNode does not match any matchingNode
+                checkTwins = .true.
+                exit                                                                ! no need to search further
+              endif
             enddo
-             neighboringElem = matchingElem(2)                            ! survival of the fittest
+            
-             neighboringIP   = j
+            !*** if no match found, then also check node twins
-             exit neighborIP
+            
-           enddo neighborIP
+            if(checkTwins) then
+periodicityDirection: do dir = 1,3
+                do a = 1,NlinkedNodes
+                  twin_of_linkedNode = mesh_nodeTwins(dir,linkedNodes(a))
+                  if (twin_of_linkedNode == 0_pInt &                                ! twin of linkedNode does not exist...
+                      .or. all(matchingNodes /= twin_of_linkedNode)) then           ! ... or it does not match any matchingNode
+                    if (dir < 3) then                                               ! no match in this direction...
+                      cycle periodicityDirection                                    ! ... so try in different direction
+                    else                                                            ! no matching in any direction...
+                      cycle checkCandidateIP                                        ! ... so check next candidateIP
+                    endif
+                  endif
+                enddo
+                exit periodicityDirection
+              enddo periodicityDirection
+            endif
+
+            !*** we found a match !!!
+            
+            mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem
+            mesh_ipNeighborhood(2,neighbor,myIP,myElem) = candidateIP
+            exit checkCandidateIP
+            
+          enddo checkCandidateIP
         endif                                                                     ! end of valid external matching
       endif                                                                       ! end of internal/external matching
-       mesh_ipNeighborhood(1,n,i,e) = neighboringElem
-       mesh_ipNeighborhood(2,n,i,e) = neighboringIP
-     enddo
     enddo
   enddo
+enddo
 
- return
+endsubroutine
-
- endsubroutine
 
 
 
@@ -2959,72 +3149,114 @@ subroutine mesh_marc_count_cpSizes (unit)
  endsubroutine
 
 
+!***********************************************************
+! assignment of twin nodes for each cp node
+!
+! allocate globals
+! _nodeTwins
+!***********************************************************
+subroutine mesh_build_nodeTwins()
+
+use prec, only: pInt, pReal
+implicit none
+
+integer(pInt) dir, &      ! direction of periodicity
+              node, &
+              minimumNode, &
+              maximumNode, &
+              n1, &
+              n2
+integer(pInt), dimension(mesh_Nnodes+1) :: minimumNodes, maximumNodes ! list of surface nodes (minimum and maximum coordinate value) with first entry giving the number of nodes
+real(pReal)   minCoord, maxCoord, &     ! extreme positions in one dimension
+              tolerance                 ! tolerance below which positions are assumed identical
+real(pReal), dimension(3) ::  distance  ! distance between two nodes in all three coordinates
+logical, dimension(mesh_Nnodes) ::  node_seen
+
+allocate(mesh_nodeTwins(3,mesh_Nnodes))
+mesh_nodeTwins = 0_pInt
+
+tolerance = 0.001_pReal * minval(mesh_ipVolume) ** 0.333_pReal
+
+do dir = 1,3                                    ! check periodicity in directions of x,y,z
+  if (mesh_periodicSurface(dir)) then           ! only if periodicity is requested
+
+    
+    !*** find out which nodes sit on the surface 
+    !*** and have a minimum or maximum position in this dimension
+    
+    minimumNodes = 0_pInt
+    maximumNodes = 0_pInt
+    minCoord = minval(mesh_node(dir,:))
+    maxCoord = maxval(mesh_node(dir,:))
+    do node = 1,mesh_Nnodes                     ! loop through all nodes and find surface nodes
+      if (abs(mesh_node(dir,node) - minCoord) <= tolerance) then
+        minimumNodes(1) = minimumNodes(1) + 1_pInt
+        minimumNodes(minimumNodes(1)+1) = node
+      elseif (abs(mesh_node(dir,node) - maxCoord) <= tolerance) then
+        maximumNodes(1) = maximumNodes(1) + 1_pInt
+        maximumNodes(maximumNodes(1)+1) = node
+      endif
+    enddo
+    
+    
+    !*** find the corresponding node on the other side with the same position in this dimension
+    
+    node_seen = .false.
+    do n1 = 1,minimumNodes(1)
+      minimumNode = minimumNodes(n1+1)
+      if (node_seen(minimumNode)) &
+        cycle
+      do n2 = 1,maximumNodes(1)
+        maximumNode = maximumNodes(n2+1)
+        distance = abs(mesh_node(:,minimumNode) - mesh_node(:,maximumNode))
+        if (sum(distance) - distance(dir) <= tolerance) then        ! minimum possible distance (within tolerance)
+          mesh_nodeTwins(dir,minimumNode) = maximumNode
+          mesh_nodeTwins(dir,maximumNode) = minimumNode
+          node_seen(maximumNode) = .true.                           ! remember this node, we don't have to look for his partner again
+          exit
+        endif
+      enddo
+    enddo
+
+  endif
+enddo
+
+endsubroutine
+
+
+
 !***********************************************************
 ! write statistics regarding input file parsing
 ! to the output file
 ! 
 !***********************************************************
- subroutine mesh_tell_statistics()
+subroutine mesh_tell_statistics()
 
- use prec, only: pInt
+use prec, only: pInt
- use math, only: math_range
+use math, only: math_range
- use IO, only: IO_error
+use IO, only: IO_error
- use debug, only: verboseDebugger
+use debug, only: verboseDebugger
 
- implicit none
+implicit none
 
- integer(pInt), dimension (:,:), allocatable :: mesh_HomogMicro
+integer(pInt), dimension (:,:), allocatable :: mesh_HomogMicro
- character(len=64) fmt
+character(len=64) fmt
 
- integer(pInt) i,e,n,f,t
+integer(pInt) i,e,n,f,t
 
- if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(110) ! no homogenization specified
+if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(110) ! no homogenization specified
- if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(120) ! no microstructure specified
+if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(120) ! no microstructure specified
  
- allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt
+allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt
- do e = 1,mesh_NcpElems
+do e = 1,mesh_NcpElems
   if (mesh_element(3,e) < 1_pInt) call IO_error(110,e) ! no homogenization specified
   if (mesh_element(4,e) < 1_pInt) call IO_error(120,e) ! no homogenization specified
   mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) = &
   mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1 ! count combinations of homogenization and microstructure
- enddo
+enddo
 
- if (verboseDebugger) then
+if (verboseDebugger) then
   !$OMP CRITICAL (write2out)
-   write(6,*)
-   write(6,*) 'Input Parser: IP COORDINATES'
-   write(6,'(a5,x,a5,3(x,a12))') 'elem','IP','x','y','z'
-   do e = 1,mesh_NcpElems
-     do i = 1,FE_Nips(mesh_element(2,e))
-       write (6,'(i5,x,i5,3(x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e)
-     enddo
-   enddo 
-   write(6,*)
-   write(6,*) "Input Parser: IP NEIGHBORHOOD"
-   write(6,*)
-   write(6,"(a10,x,a10,x,a10,x,a3,x,a13,x,a13)") "elem","IP","neighbor","","elemNeighbor","ipNeighbor"
-   do e = 1,mesh_NcpElems                  ! loop over cpElems
-     t = mesh_element(2,e)                 ! get elemType
-     do i = 1,FE_Nips(t)                   ! loop over IPs of elem
-       do n = 1,FE_NipNeighbors(t)         ! loop over neighbors of IP
-         write (6,"(i10,x,i10,x,i10,x,a3,x,i13,x,i13)") e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e)
-       enddo
-     enddo
-   enddo
-   write (6,*)
-    write (6,*) "Input Parser: ELEMENT VOLUME"
-    write (6,*)
-    write (6,"(a13,x,e15.8)") "total volume", sum(mesh_ipVolume)
-    write (6,*)
-    write (6,"(a5,x,a5,x,a15,x,a5,x,a15,x,a16)") "elem","IP","volume","face","area","-- normal --"
-    do e = 1,mesh_NcpElems
-      do i = 1,FE_Nips(mesh_element(2,e))
-        write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e)
-        do f = 1,FE_NipNeighbors(mesh_element(2,e))
-          write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e)
-        enddo
-      enddo
-    enddo
     write (6,*)
     write (6,*) "Input Parser: SUBNODE COORDINATES"
     write (6,*)
@@ -3042,11 +3274,52 @@ subroutine mesh_marc_count_cpSizes (unit)
         enddo
       enddo
     enddo
+    write(6,*)
+    write(6,*) 'Input Parser: IP COORDINATES'
+    write(6,'(a5,x,a5,3(x,a12))') 'elem','IP','x','y','z'
+    do e = 1,mesh_NcpElems
+      do i = 1,FE_Nips(mesh_element(2,e))
+        write (6,'(i5,x,i5,3(x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e)
+      enddo
+    enddo 
+    write (6,*)
+    write (6,*) "Input Parser: ELEMENT VOLUME"
+    write (6,*)
+    write (6,"(a13,x,e15.8)") "total volume", sum(mesh_ipVolume)
+    write (6,*)
+    write (6,"(a5,x,a5,x,a15,x,a5,x,a15,x,a16)") "elem","IP","volume","face","area","-- normal --"
+    do e = 1,mesh_NcpElems
+      do i = 1,FE_Nips(mesh_element(2,e))
+        write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e)
+        do f = 1,FE_NipNeighbors(mesh_element(2,e))
+          write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e)
+        enddo
+      enddo
+    enddo
+    write (6,*)
+    write (6,*) "Input Parser: NODE TWINS"
+    write (6,*)
+    write(6,'(a6,3(3(x),a6))') '  node','twin_x','twin_y','twin_z'
+    do n = 1,mesh_Nnodes                    ! loop over cpNodes
+      write(6,'(i6,3(3(x),i6))') n, mesh_nodeTwins(1:3,n)
+    enddo
+    write(6,*)
+    write(6,*) "Input Parser: IP NEIGHBORHOOD"
+    write(6,*)
+    write(6,"(a10,x,a10,x,a10,x,a3,x,a13,x,a13)") "elem","IP","neighbor","","elemNeighbor","ipNeighbor"
+    do e = 1,mesh_NcpElems                  ! loop over cpElems
+      t = mesh_element(2,e)                 ! get elemType
+      do i = 1,FE_Nips(t)                   ! loop over IPs of elem
+        do n = 1,FE_NipNeighbors(t)         ! loop over neighbors of IP
+          write (6,"(i10,x,i10,x,i10,x,a3,x,i13,x,i13)") e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e)
+        enddo
+      enddo
+    enddo
   !$OMP END CRITICAL (write2out)
- endif
+endif
 
 
- !$OMP CRITICAL (write2out)
+!$OMP CRITICAL (write2out)
   write (6,*)
   write (6,*) "Input Parser: STATISTICS"
   write (6,*)
@@ -3070,15 +3343,18 @@ subroutine mesh_marc_count_cpSizes (unit)
   do i=1,mesh_maxValStateVar(1)      ! loop over all (possibly assigned) homogenizations
     write (6,fmt) i,"| ",mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstrcutures
   enddo
- write (6,*)
+  write(6,*)
+  write(6,*) "Input Parser: ADDITIONAL MPIE OPTIONS"
+  write(6,*)
+  write(6,*) "periodic surface : ", mesh_periodicSurface
+  write(6,*)
   call flush(6)
- !$OMP END CRITICAL (write2out)
+!$OMP END CRITICAL (write2out)
 
- deallocate(mesh_HomogMicro)
+deallocate(mesh_HomogMicro)
- return
  
- endsubroutine
+endsubroutine
 
  
- END MODULE mesh
+END MODULE mesh