From d9a98417caf7a442e4b0b103b54368cd008fc28d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 15 Nov 2012 22:45:20 +0000 Subject: [PATCH] switched element library to geomType based. saves to copy same geometry description for different elements that are essentially similar regarding the IP number but differ in total node count. introduced quadratic tetrahedron (Marc element 127 -- element 157 might also work, but did not perform well in fully elastic calc so far) --- code/CPFEM.f90 | 5 +- code/DAMASK_marc.f90 | 5 +- code/constitutive.f90 | 11 +- code/constitutive_nonlocal.f90 | 47 +- code/crystallite.f90 | 9 +- code/homogenization.f90 | 4 +- code/homogenization_RGC.f90 | 6 +- code/material.f90 | 20 +- code/mesh.f90 | 1282 +++++++++++++++++--------------- 9 files changed, 729 insertions(+), 660 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 3d5aedeb8..e93cc8a70 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -277,7 +277,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt mesh_maxNips, & mesh_element, & FE_Nips, & - FE_Nnodes + FE_Nnodes, & + FE_geomtype use material, only: homogenization_maxNgrains, & microstructure_elemhomo, & material_phase @@ -532,7 +533,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element? - forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others + forall (i = 2:FE_Nips(FE_geomtype(mesh_element(2,e)))) ! copy results of first IP to all others materialpoint_P(1:3,1:3,i,e) = materialpoint_P(1:3,1:3,1,e) materialpoint_F(1:3,1:3,i,e) = materialpoint_F(1:3,1:3,1,e) materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) diff --git a/code/DAMASK_marc.f90 b/code/DAMASK_marc.f90 index 73afbce8f..24235ff56 100644 --- a/code/DAMASK_marc.f90 +++ b/code/DAMASK_marc.f90 @@ -242,7 +242,8 @@ subroutine hypela2(& mesh_node, & mesh_build_subNodeCoords, & mesh_build_ipCoordinates, & - FE_Nnodes + FE_Nnodes, & + FE_geomtype use CPFEM, only: CPFEM_initAll,CPFEM_general,CPFEM_init_done !$ use numerics, only: DAMASK_NumThreadsInt ! number of threads set by DAMASK_NUM_THREADS @@ -361,7 +362,7 @@ subroutine hypela2(& else computationMode = 3 ! plain collect endif - do node = 1,FE_Nnodes(mesh_element(2,cp_en)) + do node = 1,FE_Nnodes(FE_geomtype(mesh_element(2,cp_en))) FEnodeID = mesh_FEasCP('node',mesh_element(4+node,cp_en)) mesh_node(1:3,FEnodeID) = mesh_node0(1:3,FEnodeID) + numerics_unitlength * dispt(1:3,node) enddo diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 9d005ad37..c48571197 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -93,9 +93,10 @@ subroutine constitutive_init use mesh, only: mesh_maxNips, & mesh_NcpElems, & mesh_element, & + mesh_ipNeighborhood, & FE_Nips, & FE_NipNeighbors, & - mesh_ipNeighborhood + FE_geomtype use material, only: material_phase, & material_Nphase, & material_localFileExt, & @@ -221,7 +222,7 @@ endif do e = 1_pInt,mesh_NcpElems ! loop over elements myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = 1_pInt,FE_Nips(mesh_element(2,e)) ! loop over IPs + do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) ! loop over IPs do g = 1_pInt,myNgrains ! loop over grains select case(phase_elasticity(material_phase(g,i,e))) @@ -375,8 +376,8 @@ do e = 1_pInt,mesh_NcpElems ! loop over element constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(myInstance) case (constitutive_nonlocal_label) - select case(mesh_element(2,e)) - case (1_pInt,6_pInt,7_pInt,8_pInt,9_pInt) + select case(FE_geomtype(mesh_element(2,e))) + case (7_pInt,8_pInt,9_pInt,10_pInt) ! all fine case default call IO_error(253_pInt,e,i,g) @@ -417,7 +418,7 @@ enddo call constitutive_nonlocal_stateInit(constitutive_state0(1,1:iMax,1:eMax)) do e = 1_pInt,mesh_NcpElems ! loop over elements myNgrains = homogenization_Ngrains(mesh_element(3,e)) - forall(i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:myNgrains) + forall(i = 1_pInt:FE_Nips(FE_geomtype(mesh_element(2,e))), g = 1_pInt:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p constitutive_state(g,i,e)%p = constitutive_state0(g,i,e)%p ! need to be defined for first call of constitutive_microstructure in crystallite_init endforall diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 18dc55800..fad1e4623 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -952,8 +952,9 @@ use math, only: math_sampleGaussVar use mesh, only: mesh_ipVolume, & mesh_NcpElems, & mesh_maxNips, & + mesh_element, & FE_Nips, & - mesh_element + FE_geomtype use material, only: material_phase, & phase_plasticityInstance, & phase_plasticity @@ -1012,7 +1013,7 @@ do myInstance = 1_pInt,maxNinstance minimumIpVolume = 1e99_pReal do el = 1_pInt,mesh_NcpElems - do ip = 1_pInt,FE_Nips(mesh_element(2,el)) + do ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,el))) if (constitutive_nonlocal_label == phase_plasticity(material_phase(1,ip,el)) & .and. myInstance == phase_plasticityInstance(material_phase(1,ip,el))) then totalVolume = totalVolume + mesh_ipVolume(ip,el) @@ -1029,7 +1030,7 @@ do myInstance = 1_pInt,maxNinstance do while(meanDensity < constitutive_nonlocal_rhoSglRandom(myInstance)) call random_number(rnd) el = nint(rnd(1)*real(mesh_NcpElems,pReal)+0.5_pReal,pInt) - ip = nint(rnd(2)*real(FE_Nips(mesh_element(2,el)),pReal)+0.5_pReal,pInt) + ip = nint(rnd(2)*real(FE_Nips(FE_geomtype(mesh_element(2,el))),pReal)+0.5_pReal,pInt) if (constitutive_nonlocal_label == phase_plasticity(material_phase(1,ip,el)) & .and. myInstance == phase_plasticityInstance(material_phase(1,ip,el))) then s = nint(rnd(3)*real(ns,pReal)+0.5_pReal,pInt) @@ -1042,7 +1043,7 @@ do myInstance = 1_pInt,maxNinstance ! homogeneous distribution of density with some noise else do el = 1_pInt,mesh_NcpElems - do ip = 1_pInt,FE_Nips(mesh_element(2,el)) + do ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,el))) if (constitutive_nonlocal_label == phase_plasticity(material_phase(1,ip,el)) & .and. myInstance == phase_plasticityInstance(material_phase(1,ip,el))) then do f = 1_pInt,lattice_maxNslipFamily @@ -1177,12 +1178,13 @@ use debug, only: debug_level, & use mesh, only: mesh_NcpElems, & mesh_maxNips, & mesh_element, & - FE_NipNeighbors, & - FE_maxNipNeighbors, & mesh_ipNeighborhood, & mesh_ipCoordinates, & mesh_ipVolume, & - mesh_ipAreaNormal + mesh_ipAreaNormal, & + FE_NipNeighbors, & + FE_maxNipNeighbors, & + FE_geomtype use material, only: homogenization_maxNgrains, & material_phase, & phase_localPlasticity, & @@ -1347,7 +1349,7 @@ if (.not. phase_localPlasticity(phase) .and. constitutive_nonlocal_shortRangeStr !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities - do n = 1_pInt,FE_NipNeighbors(mesh_element(2,el)) + do n = 1_pInt,FE_NipNeighbors(FE_geomtype(mesh_element(2,el))) neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) if (neighboring_el > 0 .and. neighboring_ip > 0) then @@ -2031,11 +2033,12 @@ use math, only: math_norm3, & use mesh, only: mesh_NcpElems, & mesh_maxNips, & mesh_element, & - FE_NipNeighbors, & mesh_ipNeighborhood, & mesh_ipVolume, & mesh_ipArea, & - mesh_ipAreaNormal + mesh_ipAreaNormal, & + FE_NipNeighbors, & + FE_geomtype use material, only: homogenization_maxNgrains, & material_phase, & phase_plasticityInstance, & @@ -2300,7 +2303,7 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then my_Fe = Fe(1:3,1:3,g,ip,el) my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,g,ip,el)) - do n = 1_pInt,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors + do n = 1_pInt,FE_NipNeighbors(FE_geomtype(mesh_element(2,el))) ! loop through my neighbors neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) neighboring_n = mesh_ipNeighborhood(3,n,ip,el) @@ -2561,9 +2564,10 @@ use material, only: material_phase, & homogenization_maxNgrains use mesh, only: mesh_element, & mesh_ipNeighborhood, & - FE_NipNeighbors, & mesh_maxNips, & - mesh_NcpElems + mesh_NcpElems, & + FE_NipNeighbors, & + FE_geomtype use lattice, only: lattice_sn, & lattice_sd @@ -2594,7 +2598,7 @@ integer(pInt) Nneighbors, & ! real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor real(pReal), dimension(2,constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(1,i,e))),& constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(1,i,e))),& - FE_NipNeighbors(mesh_element(2,e))) :: & + FE_NipNeighbors(FE_geomtype(mesh_element(2,e)))) :: & compatibility ! compatibility for current element and ip real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: & slipNormal, & @@ -2606,10 +2610,10 @@ logical, dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance(mat belowThreshold -Nneighbors = FE_NipNeighbors(mesh_element(2,e)) -my_phase = material_phase(1,i,e) -my_texture = material_texture(1,i,e) -my_instance = phase_plasticityInstance(my_phase) +Nneighbors = FE_NipNeighbors(FE_geomtype(mesh_element(2,e))) +my_phase = material_phase(1,i,e) +my_texture = material_texture(1,i,e) +my_instance = phase_plasticityInstance(my_phase) my_structure = constitutive_nonlocal_structure(my_instance) ns = constitutive_nonlocal_totalNslip(my_instance) slipNormal(1:3,1:ns) = lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,my_instance), my_structure) @@ -2766,10 +2770,11 @@ use mesh, only: mesh_NcpElems, & mesh_maxNips, & mesh_element, & mesh_node0, & - FE_Nips, & mesh_cellCenterCoordinates, & mesh_ipVolume, & - mesh_periodicSurface + mesh_periodicSurface, & + FE_Nips, & + FE_geomtype use material, only: homogenization_maxNgrains, & material_phase, & phase_localPlasticity, & @@ -2890,7 +2895,7 @@ if (.not. phase_localPlasticity(phase)) then !* but only consider nonlocal neighbors within a certain cutoff radius R do neighboring_el = 1_pInt,mesh_NcpElems -ipLoop: do neighboring_ip = 1_pInt,FE_Nips(mesh_element(2,neighboring_el)) +ipLoop: do neighboring_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighboring_el))) neighboring_phase = material_phase(g,neighboring_ip,neighboring_el) if (phase_localPlasticity(neighboring_phase)) then cycle diff --git a/code/crystallite.f90 b/code/crystallite.f90 index e69deb57c..6003b3387 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -3038,7 +3038,8 @@ use material, only: material_phase, & phase_plasticityInstance use mesh, only: mesh_element, & mesh_ipNeighborhood, & - FE_NipNeighbors + FE_NipNeighbors, & + FE_geomtype use constitutive_nonlocal, only: constitutive_nonlocal_structure, & constitutive_nonlocal_updateCompatibility @@ -3104,19 +3105,19 @@ logical error ! --- calculate disorientation between me and my neighbor --- - do n = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors + do n = 1_pInt,FE_NipNeighbors(FE_geomtype(mesh_element(2,e))) ! loop through my neighbors neighboring_e = mesh_ipNeighborhood(1,n,i,e) neighboring_i = mesh_ipNeighborhood(2,n,i,e) if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase if (.not. phase_localPlasticity(neighboringPhase)) then ! neighbor got also nonlocal plasticity neighboringInstance = phase_plasticityInstance(neighboringPhase) - neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure + neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me crystallite_disorientation(:,n,1,i,e) = & math_QuaternionDisorientation( crystallite_orientation(1:4,1,i,e), & crystallite_orientation(1:4,1,neighboring_i,neighboring_e), & - crystallite_symmetryID(1,i,e)) ! calculate disorientation + crystallite_symmetryID(1,i,e)) ! calculate disorientation else ! for neighbor with different phase crystallite_disorientation(:,n,1,i,e) = (/0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal/) ! 180 degree rotation about 100 axis endif diff --git a/code/homogenization.f90 b/code/homogenization.f90 index b647f0dda..af039415d 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -79,7 +79,7 @@ subroutine homogenization_init(Temperature) use debug, only: debug_level, debug_homogenization, debug_levelBasic use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat, IO_write_jobFile, & IO_write_jobBinaryIntFile - use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips + use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips,FE_geomtype use material use constitutive, only: constitutive_maxSizePostResults use crystallite, only: crystallite_maxSizePostResults @@ -178,7 +178,7 @@ subroutine homogenization_init(Temperature) !$OMP PARALLEL DO PRIVATE(myInstance) do e = 1,mesh_NcpElems ! loop over elements myInstance = homogenization_typeInstance(mesh_element(3,e)) - do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs + do i = 1,FE_Nips(FE_geomtype(mesh_element(2,e))) ! loop over IPs select case(homogenization_type(mesh_element(3,e))) case (homogenization_isostrain_label) if (homogenization_isostrain_sizeState(myInstance) > 0_pInt) then diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 77cd981fe..3ee3c3b9b 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -77,7 +77,7 @@ subroutine homogenization_RGC_init(& math_sampleRandomOri,& math_EulerToR,& INRAD - use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips + use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips,FE_geomtype use IO use material @@ -164,7 +164,7 @@ subroutine homogenization_RGC_init(& myInstance = homogenization_typeInstance(mesh_element(3,e)) if (all (homogenization_RGC_angles(:,myInstance) >= 399.9_pReal)) then homogenization_RGC_orientation(:,:,1,e) = math_EulerToR(math_sampleRandomOri()) - do i = 1_pInt,FE_Nips(mesh_element(2,e)) + do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) if (microstructure_elemhomo(mesh_element(4,e))) then homogenization_RGC_orientation(:,:,i,e) = homogenization_RGC_orientation(:,:,1,e) else @@ -172,7 +172,7 @@ subroutine homogenization_RGC_init(& endif enddo else - do i = 1_pInt,FE_Nips(mesh_element(2,e)) + do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) homogenization_RGC_orientation(:,:,i,e) = math_EulerToR(homogenization_RGC_angles(:,myInstance)*inRad) enddo endif diff --git a/code/material.f90 b/code/material.f90 index b972eab20..00319c49d 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -672,7 +672,8 @@ subroutine material_populateGrains mesh_maxNips, & mesh_NcpElems, & mesh_ipVolume, & - FE_Nips + FE_Nips, & + FE_geomtype use IO, only: IO_error, & IO_hybridIA use FEsolving, only: FEsolving_execIP @@ -720,6 +721,7 @@ subroutine material_populateGrains ! identify maximum grain count per IP (from element) and find grains per homog/micro pair do e = 1_pInt,mesh_NcpElems + t = FE_geomtype(mesh_element(2,e)) homog = mesh_element(3,e) micro = mesh_element(4,e) if (homog < 1_pInt .or. homog > material_Nhomogenization) & ! out of bounds @@ -729,7 +731,7 @@ subroutine material_populateGrains if (microstructure_elemhomo(micro)) then dGrains = homogenization_Ngrains(homog) else - dGrains = homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e)) + dGrains = homogenization_Ngrains(homog) * FE_Nips(t) endif Ngrains(homog,micro) = Ngrains(homog,micro) + dGrains Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt @@ -767,15 +769,16 @@ subroutine material_populateGrains grain = 0_pInt do hme = 1_pInt, Nelems(homog,micro) e = elemsOfHomogMicro(hme,homog,micro) ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex + t = FE_geomtype(mesh_element(2,e)) if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs - volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(mesh_element(2,e)),e))/& + volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(t),e))/& real(dGrains,pReal) grain = grain + dGrains ! wind forward by NgrainsPerIP else - forall (i = 1_pInt:FE_Nips(mesh_element(2,e))) & ! loop over IPs + forall (i = 1_pInt:FE_Nips(t)) & ! loop over IPs volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = & mesh_ipVolume(i,e)/dGrains ! assign IPvolume/Ngrains to all grains of IP - grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP + grain = grain + FE_Nips(t) * dGrains ! wind forward by Nips*NgrainsPerIP endif enddo @@ -886,8 +889,9 @@ subroutine material_populateGrains grain = 0_pInt do hme = 1_pInt, Nelems(homog,micro) e = elemsOfHomogMicro(hme,homog,micro) ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex + t = FE_geomtype(mesh_element(2,e)) if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs - forall (i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:dGrains) ! loop over IPs and grains + forall (i = 1_pInt:FE_Nips(t), g = 1_pInt:dGrains) ! loop over IPs and grains material_volume(g,i,e) = volumeOfGrain(grain+g) material_phase(g,i,e) = phaseOfGrain(grain+g) material_texture(g,i,e) = textureOfGrain(grain+g) @@ -896,13 +900,13 @@ subroutine material_populateGrains FEsolving_execIP(2,e) = 1_pInt ! restrict calculation to first IP only, since all other results are to be copied from this grain = grain + dGrains ! wind forward by NgrainsPerIP else - forall (i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:dGrains) ! loop over IPs and grains + forall (i = 1_pInt:FE_Nips(t), g = 1_pInt:dGrains) ! loop over IPs and grains material_volume(g,i,e) = volumeOfGrain(grain+(i-1_pInt)*dGrains+g) material_phase(g,i,e) = phaseOfGrain(grain+(i-1_pInt)*dGrains+g) material_texture(g,i,e) = textureOfGrain(grain+(i-1_pInt)*dGrains+g) material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1_pInt)*dGrains+g) end forall - grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP + grain = grain + FE_Nips(t) * dGrains ! wind forward by Nips*NgrainsPerIP endif enddo endif ! active homog,micro pair diff --git a/code/mesh.f90 b/code/mesh.f90 index d70b62be1..8703ddef8 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -1,7 +1,7 @@ -! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH +! Copyright 2011 Max-Planck-Institut fŸr Eisenforschung GmbH ! ! This file is part of DAMASK, -! the Düsseldorf Advanced Material Simulation Kit. +! the DŸsseldorf Advanced Material Simulation Kit. ! ! DAMASK is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by @@ -19,11 +19,11 @@ !-------------------------------------------------------------------------------------------------- !* $Id$ !-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!! Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!! Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH -!! Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!! Krishna Komerla, Max-Planck-Institut für Eisenforschung GmbH +!> @author Franz Roters, Max-Planck-Institut fŸr Eisenforschung GmbH +!! Philip Eisenlohr, Max-Planck-Institut fŸr Eisenforschung GmbH +!! Christoph Koords, Max-Planck-Institut fŸr Eisenforschung GmbH +!! Martin Diehl, Max-Planck-Institut fŸr Eisenforschung GmbH +!! Krishna Komerla, Max-Planck-Institut fŸr Eisenforschung GmbH !> @brief Sets up the mesh for the solvers MSC.Marc, Abaqus and the spectral solver !-------------------------------------------------------------------------------------------------- @@ -109,7 +109,8 @@ module mesh ! Hence, I suggest to prefix with "FE_" integer(pInt), parameter, public :: & - FE_Nelemtypes = 10_pInt, & + FE_Nelemtypes = 12_pInt, & + FE_Ngeomtypes = 10_pInt, & FE_maxNnodes = 8_pInt, & FE_maxNsubNodes = 56_pInt, & FE_maxNips = 27_pInt, & @@ -117,155 +118,186 @@ module mesh FE_maxmaxNnodesAtIP = 8_pInt, & !< max number of (equivalent) nodes attached to an IP FE_NipFaceNodes = 4_pInt - integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nnodes = & !< nodes in a specific type of element (how we use it) - int([8, & ! element 7 - 4, & ! element 134 - 4, & ! element 11 - 4, & ! element 27 - 4, & ! element 157 - 6, & ! element 136 - 8, & ! element 21 - 8, & ! element 117 - 8, & ! element 57 (c3d20r == c3d8 --> copy of 7) - 3 & ! element 155, 125, 128 - ],pInt) - integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nips = & !< IPs in a specific type of element - int([8, & ! element 7 - 1, & ! element 134 - 4, & ! element 11 - 9, & ! element 27 - 4, & ! element 157 - 6, & ! element 136 - 27,& ! element 21 - 1, & ! element 117 - 8, & ! element 57 (c3d20r == c3d8 --> copy of 7) - 3 & ! element 155, 125, 128 - ],pInt) - integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_NipNeighbors = & !< IP neighbors in a specific type of element - int([6, & ! element 7 - 4, & ! element 134 - 4, & ! element 11 - 4, & ! element 27 - 6, & ! element 157 - 6, & ! element 136 - 6, & ! element 21 - 6, & ! element 117 - 6, & ! element 57 (c3d20r == c3d8 --> copy of 7) - 4 & ! element 155, 125, 128 - ],pInt) - integer(pInt), dimension(FE_Nelemtypes), parameter, private :: FE_NsubNodes = & !< subnodes required to fully define all IP volumes - int([19,& ! element 7 ! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype - 0, & ! element 134 - 5, & ! element 11 - 12,& ! element 27 - 0, & ! element 157 - 15,& ! element 136 - 56,& ! element 21 - 0, & ! element 117 - 19,& ! element 57 (c3d20r == c3d8 --> copy of 7) - 4 & ! element 155, 125, 128 + integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_geomtype = & !< geometry type of particular element type + int([ & + 1, & ! element 6 (2D 3node 1ip) + 2, & ! element 125 (2D 6node 3ip) + 3, & ! element 11 (2D 4node 4ip) + 4, & ! element 27 (2D 8node 9ip) + 5, & ! element 134 (3D 4node 1ip) + 6, & ! element 157 (3D 5node 4ip) + 6, & ! element 127 (3D 10node 4ip) + 7, & ! element 136 (3D 6node 6ip) + 8, & ! element 117 (3D 8node 1ip) + 9, & ! element 7 (3D 8node 8ip) + 9, & ! element 57 (3D 20node 8ip) + 10 & ! element 21 (3D 20node 27ip) ],pInt) + integer(pInt), dimension(FE_Nelemtypes), parameter, private :: FE_NoriginalNodes = & !< nodes in a specific type of element (how it is originally defined by marc) - int([8, & ! element 7 - 4, & ! element 134 - 4, & ! element 11 - 8, & ! element 27 - 4, & ! element 157 - 6, & ! element 136 - 20,& ! element 21 - 8, & ! element 117 - 20,& ! element 57 (c3d20r == c3d8 --> copy of 7) - 6 & ! element 155, 125, 128 + int([ & + 3, & ! element 6 (2D 3node 1ip) + 6, & ! element 125 (2D 6node 3ip) + 4, & ! element 11 (2D 4node 4ip) + 8, & ! element 27 (2D 8node 9ip) + 4, & ! element 134 (3D 4node 1ip) + 5, & ! element 157 (3D 5node 4ip) + 10, & ! element 127 (3D 10node 4ip) + 6, & ! element 136 (3D 6node 6ip) + 8, & ! element 117 (3D 8node 1ip) + 8, & ! element 7 (3D 8node 8ip) + 20, & ! element 57 (3D 20node 8ip) + 20 & ! element 21 (3D 20node 27ip) ],pInt) - integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter, private :: FE_NfaceNodes = &!< nodes per face in a specific type of element - reshape(int([& - 4,4,4,4,4,4, & ! element 7 - 3,3,3,3,0,0, & ! element 134 - 2,2,2,2,0,0, & ! element 11 - 2,2,2,2,0,0, & ! element 27 - 3,3,3,3,0,0, & ! element 157 - 3,4,4,4,3,0, & ! element 136 - 4,4,4,4,4,4, & ! element 21 - 4,4,4,4,4,4, & ! element 117 - 4,4,4,4,4,4, & ! element 57 (c3d20r == c3d8 --> copy of 7) - 2,2,2,0,0,0 & ! element 155, 125, 128 - ],pInt),[FE_maxNipNeighbors,FE_Nelemtypes]) - integer(pInt), dimension(FE_Nelemtypes), parameter, private :: FE_maxNnodesAtIP = & !< map IP index to two node indices in a specific type of element - int([1, & ! element 7 - 4, & ! element 134 - 1, & ! element 11 - 2, & ! element 27 - 1, & ! element 157 - 1, & ! element 136 - 4, & ! element 21 - 8, & ! element 117 - 1, & ! element 57 (c3d20r == c3d8 --> copy of 7) - 1 & ! element 155, 125, 128 + + integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_Nnodes = & !< nodes in a specific type of element (how we use it) + 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) - integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes), parameter, private :: & + + integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_Nips = & !< IPs in a specific type of element + int([ & + 1, & ! element 6 (2D 3node 1ip) + 3, & ! element 125 (2D 6node 3ip) + 4, & ! element 11 (2D 4node 4ip) + 9, & ! element 27 (2D 8node 9ip) + 1, & ! element 134 (3D 4node 1ip) + 4, & ! element 127 (3D 10node 4ip) + 6, & ! element 136 (3D 6node 6ip) + 1, & ! element 117 (3D 8node 1ip) + 8, & ! element 7 (3D 8node 8ip) + 27 & ! element 21 (3D 20node 27ip) + ],pInt) + + integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_NipNeighbors = & !< IP neighbors in a specific type of element + int([ & + 3, & ! element 6 (2D 3node 1ip) + 4, & ! element 125 (2D 6node 3ip) + 4, & ! element 11 (2D 4node 4ip) + 4, & ! element 27 (2D 8node 9ip) + 4, & ! element 134 (3D 4node 1ip) + 6, & ! element 127 (3D 10node 4ip) + 6, & ! element 136 (3D 6node 6ip) + 6, & ! element 117 (3D 8node 1ip) + 6, & ! element 7 (3D 8node 8ip) + 6 & ! element 21 (3D 20node 27ip) + ],pInt) + + integer(pInt), dimension(FE_Ngeomtypes), parameter, private :: FE_NsubNodes = & !< subnodes required to fully define all IP volumes + int([ & + 0, & ! element 6 (2D 3node 1ip) + 4, & ! element 125 (2D 6node 3ip) + 5, & ! element 11 (2D 4node 4ip) + 12, & ! element 27 (2D 8node 9ip) + 0, & ! element 134 (3D 4node 1ip) + 11, & ! element 127 (3D 10node 4ip) + 15, & ! element 136 (3D 6node 6ip) + 0, & ! element 117 (3D 8node 1ip) + 19, & ! element 7 (3D 8node 8ip) + 56 & ! element 21 (3D 20node 27ip) + ],pInt) + + integer(pInt), dimension(FE_maxNipNeighbors,FE_Ngeomtypes), parameter, private :: FE_NfaceNodes = &!< nodes per face in a specific type of element + reshape(int([ & + 2,2,2,0,0,0, & ! element 6 (2D 6node 1ip) + 2,2,2,0,0,0, & ! element 125 (2D 6node 3ip) + 2,2,2,2,0,0, & ! element 11 (2D 4node 4ip) + 2,2,2,2,0,0, & ! element 27 (2D 8node 9ip) + 3,3,3,3,0,0, & ! element 134 (3D 4node 1ip) + 3,3,3,3,0,0, & ! element 127 (3D 10node 4ip) + 3,4,4,4,3,0, & ! element 136 (3D 6node 6ip) + 4,4,4,4,4,4, & ! element 117 (3D 8node 1ip) + 4,4,4,4,4,4, & ! element 7 (3D 8node 8ip) + 4,4,4,4,4,4 & ! element 21 (3D 20node 27ip) + ],pInt),[FE_maxNipNeighbors,FE_Ngeomtypes]) + + integer(pInt), dimension(FE_Ngeomtypes), parameter, private :: FE_maxNnodesAtIP = & !< map IP index to two node indices in a specific type of element + int([ & + 3, & ! element 6 (2D 6node 1ip) + 1, & ! element 125 (2D 6node 3ip) + 1, & ! element 11 (2D 4node 4ip) + 2, & ! element 27 (2D 8node 9ip) + 4, & ! element 134 (3D 4node 1ip) + 1, & ! element 127 (3D 10node 4ip) + 1, & ! element 136 (3D 6node 6ip) + 8, & ! element 117 (3D 8node 1ip) + 1, & ! element 7 (3D 8node 8ip) + 4 & ! element 21 (3D 20node 27ip) + ],pInt) + + integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_Ngeomtypes), parameter, private :: & FE_nodeOnFace = & !< List of node indices on each face of a specific type of element reshape(int([& - 1,2,3,4 , & ! element 7 - 2,1,5,6 , & - 3,2,6,7 , & - 4,3,7,8 , & - 4,1,5,8 , & - 8,7,6,5 , & - 1,2,3,0 , & ! element 134 - 1,4,2,0 , & - 2,3,4,0 , & - 1,3,4,0 , & + 1,2,0,0 , & ! element 6 (2D 3node 1ip) + 2,3,0,0 , & + 3,1,0,0 , & 0,0,0,0 , & 0,0,0,0 , & - 1,2,0,0 , & ! element 11 + 0,0,0,0 , & + 1,2,0,0 , & ! element 125 (2D 6node 3ip) + 2,3,0,0 , & + 3,1,0,0 , & + 0,0,0,0 , & + 0,0,0,0 , & + 0,0,0,0 , & + 1,2,0,0 , & ! element 11 (2D 4node 4ip) 2,3,0,0 , & 3,4,0,0 , & 4,1,0,0 , & 0,0,0,0 , & 0,0,0,0 , & - 1,2,0,0 , & ! element 27 + 1,2,0,0 , & ! element 27 (2D 8node 9ip) 2,3,0,0 , & 3,4,0,0 , & 4,1,0,0 , & 0,0,0,0 , & 0,0,0,0 , & - 1,2,3,0 , & ! element 157 + 1,2,3,0 , & ! element 134 (3D 4node 1ip) 1,4,2,0 , & 2,3,4,0 , & 1,3,4,0 , & 0,0,0,0 , & 0,0,0,0 , & - 1,2,3,0 , & ! element 136 + 1,2,3,0 , & ! element 127 (3D 10node 4ip) + 1,4,2,0 , & + 2,4,3,0 , & + 1,3,4,0 , & + 0,0,0,0 , & + 0,0,0,0 , & + 1,2,3,0 , & ! element 136 (3D 6node 6ip) 1,4,5,2 , & 2,5,6,3 , & 1,3,6,4 , & 4,6,5,0 , & 0,0,0,0 , & - 1,2,3,4 , & ! element 21 + 1,2,3,4 , & ! element 117 (3D 8node 1ip) 2,1,5,6 , & 3,2,6,7 , & 4,3,7,8 , & 4,1,5,8 , & 8,7,6,5 , & - 1,2,3,4 , & ! element 117 + 1,2,3,4 , & ! element 7 (3D 8node 8ip) 2,1,5,6 , & 3,2,6,7 , & 4,3,7,8 , & 4,1,5,8 , & 8,7,6,5 , & - 1,2,3,4 , & ! element 57 (c3d20r == c3d8 --> copy of 7) + 1,2,3,4 , & ! element 21 (3D 20node 27ip) 2,1,5,6 , & 3,2,6,7 , & 4,3,7,8 , & 4,1,5,8 , & - 8,7,6,5 , & - 1,2,0,0 , & ! element 155,125,128 - 2,3,0,0 , & - 3,1,0,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 0,0,0,0 & - ],pInt),[FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes]) + 8,7,6,5 & + ],pInt),[FE_NipFaceNodes,FE_maxNipNeighbors,FE_Ngeomtypes]) integer(pInt), dimension(:,:,:), allocatable, private :: & FE_nodesAtIP, & !< map IP index to two node indices in a specific type of element @@ -380,23 +412,23 @@ subroutine mesh_init(ip,element) write(6,*) '$Id$' #include "compilation_info.f90" !$OMP END CRITICAL (write2out) - if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem) - if (allocated(mesh_mapFEtoCPnode)) deallocate(mesh_mapFEtoCPnode) - if (allocated(mesh_node0)) deallocate(mesh_node0) - if (allocated(mesh_node)) deallocate(mesh_node) - if (allocated(mesh_element)) deallocate(mesh_element) - if (allocated(mesh_subNodeCoord)) deallocate(mesh_subNodeCoord) - if (allocated(mesh_ipCoordinates)) deallocate(mesh_ipCoordinates) - if (allocated(mesh_ipArea)) deallocate(mesh_ipArea) - if (allocated(mesh_ipAreaNormal)) deallocate(mesh_ipAreaNormal) - if (allocated(mesh_sharedElem)) deallocate(mesh_sharedElem) + if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem) + if (allocated(mesh_mapFEtoCPnode)) deallocate(mesh_mapFEtoCPnode) + if (allocated(mesh_node0)) deallocate(mesh_node0) + if (allocated(mesh_node)) deallocate(mesh_node) + if (allocated(mesh_element)) deallocate(mesh_element) + if (allocated(mesh_subNodeCoord)) deallocate(mesh_subNodeCoord) + if (allocated(mesh_ipCoordinates)) deallocate(mesh_ipCoordinates) + if (allocated(mesh_ipArea)) deallocate(mesh_ipArea) + if (allocated(mesh_ipAreaNormal)) deallocate(mesh_ipAreaNormal) + if (allocated(mesh_sharedElem)) deallocate(mesh_sharedElem) if (allocated(mesh_ipNeighborhood)) deallocate(mesh_ipNeighborhood) - if (allocated(mesh_ipVolume)) deallocate(mesh_ipVolume) - if (allocated(mesh_nodeTwins)) deallocate(mesh_nodeTwins) - if (allocated(FE_nodesAtIP)) deallocate(FE_nodesAtIP) - if (allocated(FE_ipNeighbor))deallocate(FE_ipNeighbor) - if (allocated(FE_subNodeParent)) deallocate(FE_subNodeParent) - if (allocated(FE_subNodeOnIPFace)) deallocate(FE_subNodeOnIPFace) + if (allocated(mesh_ipVolume)) deallocate(mesh_ipVolume) + if (allocated(mesh_nodeTwins)) deallocate(mesh_nodeTwins) + if (allocated(FE_nodesAtIP)) deallocate(FE_nodesAtIP) + if (allocated(FE_ipNeighbor)) deallocate(FE_ipNeighbor) + if (allocated(FE_subNodeParent)) deallocate(FE_subNodeParent) + if (allocated(FE_subNodeOnIPFace)) deallocate(FE_subNodeOnIPFace) call mesh_build_FEdata ! get properties of the different types of elements #ifdef Spectral call IO_open_file(fileUnit,geometryFile) ! parse info from geometry file... @@ -425,15 +457,25 @@ subroutine mesh_init(ip,element) #endif #ifdef Marc call IO_open_inputFile(fileUnit,modelName) ! parse info from input file... + write(6,*) '1<<<+- mesh init -+>>>' call mesh_marc_get_tableStyles(fileUnit) + write(6,*) '2<<<+- mesh init -+>>>' call mesh_marc_count_nodesAndElements(fileUnit) + write(6,*) '3<<<+- mesh init -+>>>' call mesh_marc_count_elementSets(fileUnit) + write(6,*) '4<<<+- mesh init -+>>>' call mesh_marc_map_elementSets(fileUnit) + write(6,*) '5<<<+- mesh init -+>>>' call mesh_marc_count_cpElements(fileUnit) + write(6,*) '6<<<+- mesh init -+>>>' call mesh_marc_map_elements(fileUnit) + write(6,*) '7<<<+- mesh init -+>>>' call mesh_marc_map_nodes(fileUnit) + write(6,*) '8<<<+- mesh init -+>>>' call mesh_marc_build_nodes(fileUnit) + write(6,*) '9<<<+- mesh init -+>>>' call mesh_marc_count_cpSizes(fileunit) + write(6,*) '10<<<+- mesh init -+>>>' call mesh_marc_build_elements(fileUnit) #endif #ifdef Abaqus @@ -455,13 +497,21 @@ subroutine mesh_init(ip,element) call mesh_get_damaskOptions(fileUnit) close (fileUnit) + write(6,*) '11<<<+- mesh init -+>>>' call mesh_build_subNodeCoords + write(6,*) '12<<<+- mesh init -+>>>' call mesh_build_ipCoordinates + write(6,*) '13<<<+- mesh init -+>>>' call mesh_build_ipVolumes + write(6,*) '14<<<+- mesh init -+>>>' call mesh_build_ipAreas + write(6,*) '15<<<+- mesh init -+>>>' call mesh_build_nodeTwins + write(6,*) '16<<<+- mesh init -+>>>' call mesh_build_sharedElems + write(6,*) '17<<<+- mesh init -+>>>' call mesh_build_ipNeighborhood + write(6,*) 'stat<<<+- mesh init -+>>>' call mesh_tell_statistics parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive @@ -469,7 +519,7 @@ subroutine mesh_init(ip,element) FEsolving_execElem = [ 1_pInt,mesh_NcpElems] if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP) allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt - forall (e = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e)) + forall (e = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(FE_geomtype(mesh_element(2,e))) if (allocated(calcMode)) deallocate(calcMode) allocate(calcMode(mesh_maxNips,mesh_NcpElems)) @@ -548,7 +598,7 @@ subroutine mesh_build_subNodeCoords !$OMP PARALLEL DO PRIVATE(mySubNodeCoord,t,Nparents) do e = 1_pInt,mesh_NcpElems ! loop over cpElems mySubNodeCoord = 0.0_pReal - t = mesh_element(2,e) ! get elemType + t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType do n = 1_pInt,FE_Nnodes(t) mySubNodeCoord(1:3,n) = mesh_node(1:3,mesh_FEasCP('node',mesh_element(4_pInt+n,e))) ! loop over nodes of this element type enddo @@ -587,7 +637,7 @@ subroutine mesh_build_ipVolumes mesh_ipVolume = 0.0_pReal do e = 1_pInt,mesh_NcpElems ! loop over cpElems - t = mesh_element(2_pInt,e) ! get elemType + t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG forall (n = 1_pInt:FE_NipFaceNodes) & @@ -631,7 +681,7 @@ subroutine mesh_build_ipCoordinates !$OMP PARALLEL DO PRIVATE(t,gravityNode,gravityNodePos) do e = 1_pInt,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType + t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem gravityNode = .false. ! reset flagList gravityNodePos = 0.0_pReal ! reset coordinates @@ -683,7 +733,7 @@ logical, dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav -t = mesh_element(2,e) ! get elemType +t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType gravityNode = .false. ! reset flagList gravityNodePos = 0.0_pReal ! reset coordinates do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP @@ -997,7 +1047,7 @@ subroutine mesh_spectral_count_cpSizes implicit none integer(pInt) :: t - t = FE_mapElemtype('C3D8R') ! fake 3D hexahedral 8 node 1 IP element + t = FE_geomtype(FE_mapElemtype('C3D8R')) ! fake 3D hexahedral 8 node 1 IP element mesh_maxNnodes = FE_Nnodes(t) mesh_maxNips = FE_Nips(t) @@ -2282,10 +2332,13 @@ subroutine mesh_marc_count_nodesAndElements(myUnit) read (myUnit,610,END=620) line myPos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'sizing') then + if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'sizing') & mesh_Nelems = IO_IntValue (line,myPos,3_pInt) - mesh_Nnodes = IO_IntValue (line,myPos,4_pInt) - exit + if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'coordinates') then + read (myUnit,610,END=620) line + myPos = IO_stringPos(line,maxNchunks) + mesh_Nnodes = IO_IntValue (line,myPos,2_pInt) + exit ! assumes that "coordinates" comes later in file endif enddo @@ -2573,7 +2626,7 @@ subroutine mesh_marc_count_cpSizes(myUnit) integer(pInt), parameter :: maxNchunks = 2_pInt integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) :: line - integer(pInt) :: i,t,e + integer(pInt) :: i,t,g,e mesh_maxNnodes = 0_pInt mesh_maxNips = 0_pInt @@ -2593,10 +2646,11 @@ subroutine mesh_marc_count_cpSizes(myUnit) e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) if (e /= 0_pInt) then t = FE_mapElemtype(IO_stringValue(line,myPos,2_pInt)) - mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) - mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) - mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) - mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) + g = FE_geomtype(t) + mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(g)) + mesh_maxNips = max(mesh_maxNips,FE_Nips(g)) + mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(g)) + mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(g)) call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-2_pInt)) ! read on if FE_Nnodes exceeds node count present on current line endif enddo @@ -2624,12 +2678,12 @@ subroutine mesh_marc_build_elements(myUnit) implicit none integer(pInt), intent(in) :: myUnit - integer(pInt), parameter :: maxNchunks = 66_pInt + integer(pInt), parameter :: maxNchunks = 66_pInt ! limit to 64 nodes max (plus ID, type) integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) line integer(pInt), dimension(1_pInt+mesh_NcpElems) :: contInts - integer(pInt) :: i,j,sv,myVal,e + integer(pInt) :: i,j,t,sv,myVal,e allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt @@ -2640,17 +2694,18 @@ subroutine mesh_marc_build_elements(myUnit) read (myUnit,610,END=620) line myPos(1:1+2*1) = IO_stringPos(line,1_pInt) if( IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'connectivity' ) then - read (myUnit,610,END=620) line ! Garbage line + read (myUnit,610,END=620) line ! garbage line do i = 1_pInt,mesh_Nelems read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type) + myPos = IO_stringPos(line,maxNchunks) e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) if (e /= 0_pInt) then ! disregard non CP elems + t = FE_mapElemtype(IO_StringValue(line,myPos,2_pInt)) ! elem type + mesh_element(2,e) = t mesh_element(1,e) = IO_IntValue (line,myPos,1_pInt) ! FE id - mesh_element(2,e) = FE_mapElemtype(IO_StringValue(line,myPos,2_pInt)) ! elem type - forall (j = 1_pInt:FE_Nnodes(mesh_element(2,e))) & + forall (j = 1_pInt:FE_Nnodes(FE_geomtype(t))) & mesh_element(j+4_pInt,e) = IO_IntValue(line,myPos,j+2_pInt) ! copy FE ids of nodes - call IO_skipChunks(myUnit,FE_NoriginalNodes(mesh_element(2_pInt,e))-(myPos(1_pInt)-2_pInt)) ! read on if FE_Nnodes exceeds node count present on current line + call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-2_pInt)) ! read on if FE_Nnodes exceeds node count present on current line endif enddo exit @@ -2663,7 +2718,7 @@ subroutine mesh_marc_build_elements(myUnit) myPos(1:1+2*2) = IO_stringPos(line,2_pInt) if( (IO_lc(IO_stringValue(line,myPos,1_pInt)) == 'initial') .and. & (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'state') ) then - if (initialcondTableStyle == 2_pInt) read (myUnit,610,END=620) line ! read extra line for new style + if (initialcondTableStyle == 2_pInt) read (myUnit,610,END=620) line ! read extra line for new style read (myUnit,610,END=630) line ! read line with index of state var myPos(1:1+2*1) = IO_stringPos(line,1_pInt) sv = IO_IntValue(line,myPos,1_pInt) ! figure state variable index @@ -2674,8 +2729,8 @@ subroutine mesh_marc_build_elements(myUnit) myVal = nint(IO_fixedNoEFloatValue(line,[0_pInt,20_pInt],1_pInt),pInt) ! state var's value mesh_maxValStateVar(sv-1_pInt) = max(myVal,mesh_maxValStateVar(sv-1_pInt)) ! remember max val of homogenization and microstructure index if (initialcondTableStyle == 2_pInt) then - read (myUnit,610,END=630) line ! read extra line - read (myUnit,610,END=630) line ! read extra line + read (myUnit,610,END=630) line ! read extra line + read (myUnit,610,END=630) line ! read extra line endif contInts = IO_continuousIntValues& ! get affected elements (myUnit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) @@ -3232,7 +3287,7 @@ subroutine mesh_abaqus_count_cpSizes(myUnit) integer(pInt), parameter :: maxNchunks = 2_pInt integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos character(len=300) :: line - integer(pInt) :: i,c,t + integer(pInt) :: i,c,t,g logical :: inPart mesh_maxNnodes = 0_pInt @@ -3258,21 +3313,12 @@ subroutine mesh_abaqus_count_cpSizes(myUnit) IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & ) then t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type - if (t==0_pInt) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_count_cpSizes') - c = IO_countDataLines(myUnit) - do i = 1_pInt,c - backspace(myUnit) - enddo - do i = 1_pInt,c - read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max - if (mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) /= 0_pInt) then ! disregard non CP elems - mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) - mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) - mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) - mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) - endif - enddo + if (t == 0_pInt) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_count_cpSizes') + g = FE_geomtype(t) + mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(g)) + mesh_maxNips = max(mesh_maxNips,FE_Nips(g)) + mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(g)) + mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(g)) endif enddo @@ -3325,29 +3371,29 @@ subroutine mesh_abaqus_build_elements(myUnit) IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'matrix' .and. & IO_lc(IO_stringValue(line,myPos,2_pInt)) /= 'response' ) & ) then - t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type - if (t==0_pInt) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_build_elements') + t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'type')) ! remember elem type + if (t == 0_pInt) call IO_error(error_ID=910_pInt,ext_msg='mesh_abaqus_build_elements') c = IO_countDataLines(myUnit) do i = 1_pInt,c backspace(myUnit) enddo do i = 1_pInt,c read (myUnit,610,END=620) line - myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max + myPos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max e = mesh_FEasCP('elem',IO_intValue(line,myPos,1_pInt)) - if (e /= 0_pInt) then ! disregard non CP elems - mesh_element(1,e) = IO_intValue(line,myPos,1_pInt) ! FE id - mesh_element(2,e) = t ! elem type - forall (j=1_pInt:FE_Nnodes(t)) & - mesh_element(4_pInt+j,e) = IO_intValue(line,myPos,1_pInt+j) ! copy FE ids of nodes to position 5: - call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-1_pInt)) ! read on (even multiple lines) if FE_NoriginalNodes exceeds required node count + if (e /= 0_pInt) then ! disregard non CP elems + mesh_element(1,e) = IO_intValue(line,myPos,1_pInt) ! FE id + mesh_element(2,e) = t ! elem type + forall (j=1_pInt:FE_Nnodes(FE_geomtype(t))) & + mesh_element(4_pInt+j,e) = IO_intValue(line,myPos,1_pInt+j) ! copy FE ids of nodes to position 5: + call IO_skipChunks(myUnit,FE_NoriginalNodes(t)-(myPos(1_pInt)-1_pInt)) ! read on (even multiple lines) if FE_NoriginalNodes exceeds required node count endif enddo endif enddo -620 rewind(myUnit) ! just in case "*material" definitions apear before "*element" +620 rewind(myUnit) ! just in case "*material" definitions apear before "*element" materialFound = .false. do @@ -3477,7 +3523,7 @@ subroutine mesh_build_ipAreas 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 do e = 1_pInt,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType + t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem 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) @@ -3605,8 +3651,7 @@ allocate(node_seen(maxval(FE_Nnodes))) node_count = 0_pInt do e = 1_pInt,mesh_NcpElems - t = mesh_element(2,e) ! get element type - + t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType node_seen = 0_pInt ! reset node duplicates do j = 1_pInt,FE_Nnodes(t) ! check each node of element node = mesh_FEasCP('node',mesh_element(4+j,e)) ! translate to internal (consecutive) numbering @@ -3628,7 +3673,7 @@ allocate(mesh_sharedElem(1+mesh_maxNsharedElems,mesh_Nnodes)) mesh_sharedElem = 0_pInt do e = 1_pInt,mesh_NcpElems - t = mesh_element(2,e) + t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType node_seen = 0_pInt do j = 1_pInt,FE_Nnodes(t) node = mesh_FEasCP('node',mesh_element(4_pInt+j,e)) @@ -3691,7 +3736,7 @@ mesh_ipNeighborhood = 0_pInt do myElem = 1_pInt,mesh_NcpElems ! loop over cpElems - myType = mesh_element(2,myElem) ! get elemType + myType = FE_geomtype(mesh_element(2,myElem)) ! get elemGeomType do myIP = 1_pInt,FE_Nips(myType) ! loop over IPs of elem do neighbor = 1_pInt,FE_NipNeighbors(myType) ! loop over neighbors of IP @@ -3710,9 +3755,9 @@ do myElem = 1_pInt,mesh_NcpElems elseif (neighboringIPkey < 0_pInt) 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) + call mesh_faceMatch(myElem, myFace, matchingElem, matchingFace) ! get face and CP elem id of face match + if (matchingElem > 0_pInt) then ! found match? + neighboringType = FE_geomtype(mesh_element(2,matchingElem)) !*** trivial solution if neighbor has only one IP @@ -3728,12 +3773,12 @@ do myElem = 1_pInt,mesh_NcpElems linkedNodes = 0_pInt do a = 1_pInt,FE_maxNnodesAtIP(myType) ! figure my anchor nodes on connecting face anchor = FE_nodesAtIP(a,myIP,myType) - if (anchor /= 0_pInt) then ! valid anchor node - if (any(FE_nodeOnFace(:,myFace,myType) == anchor)) then ! ip anchor sits on face? + if (anchor /= 0_pInt) then ! valid anchor node + if (any(FE_nodeOnFace(:,myFace,myType) == anchor)) then ! ip anchor sits on face? NlinkedNodes = NlinkedNodes + 1_pInt linkedNodes(NlinkedNodes) = & mesh_FEasCP('node',mesh_element(4_pInt+anchor,myElem)) ! CP id of anchor node - else ! something went wrong with the linkage, since not all anchors sit on my face + else ! something went wrong with the linkage, since not all anchors sit on my face NlinkedNodes = 0_pInt linkedNodes = 0_pInt exit @@ -3801,13 +3846,13 @@ checkCandidateIP: do candidateIP = 1_pInt,FE_Nips(neighboringType) enddo enddo do myElem = 1_pInt,mesh_NcpElems ! loop over cpElems - myType = mesh_element(2,myElem) ! get elemType + myType = FE_geomtype(mesh_element(2,myElem)) ! get elemGeomType do myIP = 1_pInt,FE_Nips(myType) ! loop over IPs of elem do neighbor = 1_pInt,FE_NipNeighbors(myType) ! loop over neighbors of IP neighboringElem = mesh_ipNeighborhood(1,neighbor,myIP,myElem) neighboringIP = mesh_ipNeighborhood(2,neighbor,myIP,myElem) if (neighboringElem > 0_pInt .and. neighboringIP > 0_pInt) then ! if neighbor exists ... - neighboringType = mesh_element(2,neighboringElem) + neighboringType = FE_geomtype(mesh_element(2,neighboringElem)) do pointingToMe = 1_pInt,FE_NipNeighbors(neighboringType) ! find neighboring index that points from my neighbor to myself if ( myElem == mesh_ipNeighborhood(1,pointingToMe,neighboringIP,neighboringElem) & .and. myIP == mesh_ipNeighborhood(2,pointingToMe,neighboringIP,neighboringElem)) then ! possible candidate @@ -3899,13 +3944,13 @@ enddo write(6,*) write(6,'(a8,1x,a5,1x,2(a15,1x),a20,3(1x,a12))')& 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z' - do e = 1_pInt,mesh_NcpElems ! loop over cpElems + do e = 1_pInt,mesh_NcpElems ! loop over cpElems if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle - t = mesh_element(2,e) ! get elemType - do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem + t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle - do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP - do n = 1_pInt,FE_NipFaceNodes ! loop over nodes on interface + do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP + do n = 1_pInt,FE_NipFaceNodes ! loop over nodes on interface write(6,'(i8,1x,i5,2(1x,i15),1x,i20,3(1x,f12.8))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),& mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),& mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),& @@ -3919,7 +3964,7 @@ enddo write(6,'(a8,1x,a5,3(1x,a12))') 'elem','IP','x','y','z' do e = 1_pInt,mesh_NcpElems if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle - do i = 1_pInt,FE_Nips(mesh_element(2,e)) + do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle write(6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCoordinates(:,i,e) enddo @@ -3932,10 +3977,11 @@ enddo write(6,'(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)') 'elem','IP','volume','face','area','-- normal --' do e = 1_pInt,mesh_NcpElems if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle - do i = 1_pInt,FE_Nips(mesh_element(2,e)) + t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType + do i = 1_pInt,FE_Nips(t) if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle write(6,'(i8,1x,i5,1x,e15.8)') e,i,mesh_IPvolume(i,e) - do f = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) + do f = 1_pInt,FE_NipNeighbors(t) write(6,'(i33,1x,e15.8,1x,3(f6.3,1x))') f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) enddo enddo @@ -3955,12 +4001,12 @@ enddo write(6,*) 'Input Parser: IP NEIGHBORHOOD' write(6,*) write(6,'(a8,1x,a10,1x,a10,1x,a3,1x,a13,1x,a13)') 'elem','IP','neighbor','','elemNeighbor','ipNeighbor' - do e = 1_pInt,mesh_NcpElems ! loop over cpElems + do e = 1_pInt,mesh_NcpElems ! loop over cpElems if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle - t = mesh_element(2,e) ! get elemType - do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem + t = FE_geomtype(mesh_element(2,e)) ! get elemGeomType + do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle - do n = 1_pInt,FE_NipNeighbors(t) ! loop over neighbors of IP + do n = 1_pInt,FE_NipNeighbors(t) ! loop over neighbors of IP write(6,'(i8,1x,i10,1x,i10,1x,a3,1x,i13,1x,i13)') e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) enddo enddo @@ -3984,40 +4030,44 @@ integer(pInt) function FE_mapElemtype(what) character(len=*), intent(in) :: what select case (IO_lc(what)) - case ( '7', & - 'c3d8') - FE_mapElemtype = 1_pInt ! Three-dimensional Arbitrarily Distorted Brick - case ('134', & - 'c3d4') - FE_mapElemtype = 2_pInt ! Three-dimensional Four-node Tetrahedron + case ( '6') + FE_mapElemtype = 1_pInt ! Two-dimensional Plane Strain Triangle + case ( '155', & + '125', & + '128') + FE_mapElemtype = 2_pInt ! Two-dimensional Plane Strain triangle (155: cubic shape function, 125/128: second order isoparametric) case ( '11', & 'cpe4') FE_mapElemtype = 3_pInt ! Arbitrary Quadrilateral Plane-strain case ( '27', & 'cpe8') - FE_mapElemtype = 4_pInt ! Plane Strain, Eight-node Distorted Quadrilateral + FE_mapElemtype = 4_pInt ! Plane Strain, Eight-node Distorted Quadrilateral + case ('134', & + 'c3d4') + FE_mapElemtype = 5_pInt ! Three-dimensional Four-node Tetrahedron case ('157') - FE_mapElemtype = 5_pInt ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations + FE_mapElemtype = 6_pInt ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations + case ('127') + FE_mapElemtype = 7_pInt ! Three-dimensional Ten-node Tetrahedron case ('136', & 'c3d6') - FE_mapElemtype = 6_pInt ! Three-dimensional Arbitrarily Distorted Pentahedral - case ( '21', & - 'c3d20') - FE_mapElemtype = 7_pInt ! Three-dimensional Arbitrarily Distorted quadratic hexahedral + FE_mapElemtype = 8_pInt ! Three-dimensional Arbitrarily Distorted Pentahedral case ( '117', & '123', & 'c3d8r') - FE_mapElemtype = 8_pInt ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration + FE_mapElemtype = 9_pInt ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration + case ( '7', & + 'c3d8') + FE_mapElemtype = 10_pInt ! Three-dimensional Arbitrarily Distorted Brick case ( '57', & 'c3d20r') - FE_mapElemtype = 9_pInt ! Three-dimensional Arbitrarily Distorted quad hexahedral with reduced integration - case ( '155', & - '125', & - '128') - FE_mapElemtype = 10_pInt ! Two-dimensional Plane Strain triangle (155: cubic shape function, 125/128: second order isoparametric) + FE_mapElemtype = 11_pInt ! Three-dimensional Arbitrarily Distorted quad hexahedral with reduced integration + case ( '21', & + 'c3d20') + FE_mapElemtype = 12_pInt ! Three-dimensional Arbitrarily Distorted quadratic hexahedral case default FE_mapElemtype = 0_pInt ! unknown element --> should raise an error upstream..! - endselect + end select end function FE_mapElemtype @@ -4037,7 +4087,7 @@ integer(pInt), intent(in) :: face, & elem ! FE elem ID !*** local variables -integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: & +integer(pInt), dimension(FE_NfaceNodes(face,FE_geomtype(mesh_element(2,elem)))) :: & myFaceNodes ! global node ids on my face integer(pInt) :: myType, & candidateType, & @@ -4056,7 +4106,7 @@ logical checkTwins matchingElem = 0_pInt matchingFace = 0_pInt minNsharedElems = mesh_maxNsharedElems + 1_pInt ! init to worst case -myType = mesh_element(2_pInt,elem) ! figure elemType +myType = FE_geomtype(mesh_element(2_pInt,elem)) ! figure elemGeomType do n = 1_pInt,FE_NfaceNodes(face,myType) ! loop over nodes on face myFaceNodes(n) = mesh_FEasCP('node',mesh_element(4_pInt+FE_nodeOnFace(n,face,myType),elem)) ! CP id of face node @@ -4074,14 +4124,14 @@ checkCandidate: do i = 1_pInt,minNsharedElems candidateElem = mesh_sharedElem(1_pInt+i,myFaceNodes(lonelyNode)) ! present candidate elem if (all(element_seen /= candidateElem)) then ! element seen for the first time? element_seen(i) = candidateElem - candidateType = mesh_element(2_pInt,candidateElem) ! figure elemType of candidate + candidateType = FE_geomtype(mesh_element(2_pInt,candidateElem)) ! figure elemGeomType of candidate checkCandidateFace: do candidateFace = 1_pInt,FE_maxNipNeighbors ! check each face of candidate - if (FE_NfaceNodes(candidateFace,candidateType) /= FE_NfaceNodes(face,myType) & ! incompatible face - .or. (candidateElem == elem .and. candidateFace == face)) then ! this is my face + if (FE_NfaceNodes(candidateFace,candidateType) /= FE_NfaceNodes(face,myType) & ! incompatible face + .or. (candidateElem == elem .and. candidateFace == face)) then ! this is my face cycle checkCandidateFace endif checkTwins = .false. - do n = 1_pInt,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face + do n = 1_pInt,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face candidateFaceNode = mesh_FEasCP('node', mesh_element(4_pInt+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem)) if (all(myFaceNodes /= candidateFaceNode)) then ! candidate node does not match any of my face nodes checkTwins = .true. ! perhaps the twin nodes do match @@ -4124,38 +4174,40 @@ end subroutine mesh_faceMatch subroutine mesh_build_FEdata implicit none - allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Nelemtypes)) ; FE_nodesAtIP = 0_pInt - allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_ipNeighbor = 0_pInt - allocate(FE_subNodeParent(FE_maxNips,FE_maxNsubNodes,FE_Nelemtypes)) ; FE_subNodeParent = 0_pInt - allocate(FE_subNodeOnIPFace(FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_subNodeOnIPFace = 0_pInt + integer(pInt) :: me + allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Ngeomtypes)) ; FE_nodesAtIP = 0_pInt + allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Ngeomtypes)) ; FE_ipNeighbor = 0_pInt + allocate(FE_subNodeParent(FE_maxNips,FE_maxNsubNodes,FE_Ngeomtypes)) ; FE_subNodeParent = 0_pInt + allocate(FE_subNodeOnIPFace(FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Ngeomtypes)) ; FE_subNodeOnIPFace = 0_pInt ! fill FE_nodesAtIP with data - FE_nodesAtIP(1:FE_maxNnodesAtIP(1),1:FE_Nips(1),1) = & ! element 7 + me = 0_pInt + + me = me + 1_pInt + FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 6 (2D 3node 1ip) + reshape(int([& + 1,2,3 & + ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 125 (2D 6node 3ip) reshape(int([& 1, & 2, & - 4, & - 3, & - 5, & - 6, & - 8, & - 7 & - ],pInt),[FE_maxNnodesAtIP(1),FE_Nips(1)]) + 3 & + ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - FE_nodesAtIP(1:FE_maxNnodesAtIP(2),1:FE_Nips(2),2) = & ! element 134 - reshape(int([& - 1,2,3,4 & - ],pInt),[FE_maxNnodesAtIP(2),FE_Nips(2)]) - - FE_nodesAtIP(1:FE_maxNnodesAtIP(3),1:FE_Nips(3),3) = & ! element 11 + me = me + 1_pInt + FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 11 (2D 4node 4ip) reshape(int([& 1, & 2, & 4, & 3 & - ],pInt),[FE_maxNnodesAtIP(3),FE_Nips(3)]) + ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - FE_nodesAtIP(1:FE_maxNnodesAtIP(4),1:FE_Nips(4),4) = & ! element 27 + me = me + 1_pInt + FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 27 (2D 8node 9ip) reshape(int([& 1,0, & 1,2, & @@ -4166,17 +4218,25 @@ subroutine mesh_build_FEdata 4,0, & 3,4, & 3,0 & - ],pInt),[FE_maxNnodesAtIP(4),FE_Nips(4)]) + ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - FE_nodesAtIP(1:FE_maxNnodesAtIP(5),1:FE_Nips(5),5) = & ! element 157 + me = me + 1_pInt + FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 134 (3D 4node 1ip) + reshape(int([& + 1,2,3,4 & + ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 127 (3D 10node 4ip) reshape(int([& 1, & 2, & 3, & 4 & - ],pInt),[FE_maxNnodesAtIP(5),FE_Nips(5)]) + ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - FE_nodesAtIP(1:FE_maxNnodesAtIP(6),1:FE_Nips(6),6) = & ! element 136 + me = me + 1_pInt + FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 136 (3D 6node 6ip) reshape(int([& 1, & 2, & @@ -4184,9 +4244,29 @@ subroutine mesh_build_FEdata 4, & 5, & 6 & - ],pInt),[FE_maxNnodesAtIP(6),FE_Nips(6)]) + ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) - FE_nodesAtIP(1:FE_maxNnodesAtIP(7),1:FE_Nips(7),7) = & ! element 21 + me = me + 1_pInt + FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 117 (3D 8node 1ip) + reshape(int([& + 1,2,3,4,5,6,7,8 & + ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 7 (3D 8node 8ip) + reshape(int([& + 1, & + 2, & + 4, & + 3, & + 5, & + 6, & + 8, & + 7 & + ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_nodesAtIP(1:FE_maxNnodesAtIP(me),1:FE_Nips(me),me) = & ! element 21 (3D 20node 27ip) reshape(int([& 1,0, 0,0, & 1,2, 0,0, & @@ -4215,64 +4295,41 @@ subroutine mesh_build_FEdata 8,0, 0,0, & 7,8, 0,0, & 7,0, 0,0 & - ],pInt),[FE_maxNnodesAtIP(7),FE_Nips(7)]) + ],pInt),[FE_maxNnodesAtIP(me),FE_Nips(me)]) -! FE_nodesAtIP(:,:FE_Nips(8),8) = & ! element 117 (c3d8r --> single IP per element, so no need for this mapping) -! reshape((/& -! 1,2,3,4,5,6,7,8 & -! /),(/FE_maxNnodesAtIP(8),FE_Nips(8)/)) - - FE_nodesAtIP(1:FE_maxNnodesAtIP(9),1:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape(int([& - 1, & - 2, & - 4, & - 3, & - 5, & - 6, & - 8, & - 7 & - ],pInt),[FE_maxNnodesAtIP(9),FE_Nips(9)]) - - FE_nodesAtIP(1:FE_maxNnodesAtIP(10),1:FE_Nips(10),10) = & ! element 155, 125, 128 - reshape(int([& - 1, & - 2, & - 3 & - ],pInt),[FE_maxNnodesAtIP(10),FE_Nips(10)]) ! *** FE_ipNeighbor *** ! is a list of the neighborhood of each IP. ! It is sorted in (local) +x,-x, +y,-y, +z,-z direction. ! Positive integers denote an intra-FE IP identifier. ! Negative integers denote the interface behind which the neighboring (extra-FE) IP will be located. + me = 0_pInt - FE_ipNeighbor(1:FE_NipNeighbors(1),1:FE_Nips(1),1) = & ! element 7 + me = me + 1_pInt + FE_ipNeighbor(1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 6 (2D 3node 1ip) reshape(int([& - 2,-5, 3,-2, 5,-1, & - -3, 1, 4,-2, 6,-1, & - 4,-5,-4, 1, 7,-1, & - -3, 3,-4, 2, 8,-1, & - 6,-5, 7,-2,-6, 1, & - -3, 5, 8,-2,-6, 2, & - 8,-5,-4, 5,-6, 3, & - -3, 7,-4, 6,-6, 4 & - ],pInt),[FE_NipNeighbors(1),FE_Nips(1)]) - - FE_ipNeighbor(1:FE_NipNeighbors(2),1:FE_Nips(2),2) = & ! element 134 + -2,-3,-1 & + ],pInt),[FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_ipNeighbor(1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 125 (2D 6node 3ip) reshape(int([& - -1,-2,-3,-4 & - ],pInt),[FE_NipNeighbors(2),FE_Nips(2)]) - - FE_ipNeighbor(1:FE_NipNeighbors(3),1:FE_Nips(3),3) = & ! element 11 + 2,-3, 3,-1, & + -2, 1, 3,-1, & + 2,-3,-2, 1 & + ],pInt),[FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_ipNeighbor(1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 11 (2D 4node 4ip) reshape(int([& 2,-4, 3,-1, & -2, 1, 4,-1, & 4,-4,-3, 1, & -2, 3,-3, 2 & - ],pInt),[FE_NipNeighbors(3),FE_Nips(3)]) + ],pInt),[FE_NipNeighbors(me),FE_Nips(me)]) - FE_ipNeighbor(1:FE_NipNeighbors(4),1:FE_Nips(4),4) = & ! element 27 + me = me + 1_pInt + FE_ipNeighbor(1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 27 (2D 8node 9ip) reshape(int([& 2,-4, 4,-1, & 3, 1, 5,-1, & @@ -4283,17 +4340,25 @@ subroutine mesh_build_FEdata 8,-4,-3, 4, & 9, 7,-3, 5, & -2, 8,-3, 6 & - ],pInt),[FE_NipNeighbors(4),FE_Nips(4)]) + ],pInt),[FE_NipNeighbors(me),FE_Nips(me)]) - FE_ipNeighbor(1:FE_NipNeighbors(5),1:FE_Nips(5),5) = & ! element 157 + me = me + 1_pInt + FE_ipNeighbor(1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 134 (3D 4node 1ip) + reshape(int([& + -1,-2,-3,-4 & + ],pInt),[FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_ipNeighbor(1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 127 (3D 10node 4ip) reshape(int([& 2,-4, 3,-2, 4,-1, & - 3,-2, 1,-3, 4,-1, & - 1,-3, 2,-4, 4,-1, & - 1,-3, 2,-4, 3,-2 & - ],pInt),[FE_NipNeighbors(5),FE_Nips(5)]) + -2, 1, 3,-2, 4,-1, & + 2,-4,-3, 1, 4,-1, & + 2,-4, 3,-2,-3, 1 & + ],pInt),[FE_NipNeighbors(me),FE_Nips(me)]) - FE_ipNeighbor(1:FE_NipNeighbors(6),1:FE_Nips(6),6) = & ! element 136 + me = me + 1_pInt + FE_ipNeighbor(1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 136 (3D 6node 6ip) reshape(int([& 2,-4, 3,-2, 4,-1, & -3, 1, 3,-2, 5,-1, & @@ -4301,9 +4366,29 @@ subroutine mesh_build_FEdata 5,-4, 6,-2,-5, 1, & -3, 4, 6,-2,-5, 2, & 5,-4,-3, 4,-5, 3 & - ],pInt),[FE_NipNeighbors(6),FE_Nips(6)]) + ],pInt),[FE_NipNeighbors(me),FE_Nips(me)]) - FE_ipNeighbor(1:FE_NipNeighbors(7),1:FE_Nips(7),7) = & ! element 21 + me = me + 1_pInt + FE_ipNeighbor(1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 117 (3D 8node 1ip) + reshape(int([& + -3,-5,-4,-2,-6,-1 & + ],pInt),[FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_ipNeighbor(1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 7 (3D 8node 8ip) + reshape(int([& + 2,-5, 3,-2, 5,-1, & + -3, 1, 4,-2, 6,-1, & + 4,-5,-4, 1, 7,-1, & + -3, 3,-4, 2, 8,-1, & + 6,-5, 7,-2,-6, 1, & + -3, 5, 8,-2,-6, 2, & + 8,-5,-4, 5,-6, 3, & + -3, 7,-4, 6,-6, 4 & + ],pInt),[FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_ipNeighbor(1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 21 (3D 20node 27ip) reshape(int([& 2,-5, 4,-2,10,-1, & 3, 1, 5,-2,11,-1, & @@ -4332,32 +4417,9 @@ subroutine mesh_build_FEdata 26,-5,-4,22,-6,16, & 27,25,-4,23,-6,17, & -3,26,-4,24,-6,18 & - ],pInt),[FE_NipNeighbors(7),FE_Nips(7)]) + ],pInt),[FE_NipNeighbors(me),FE_Nips(me)]) -FE_ipNeighbor(1:FE_NipNeighbors(8),1:FE_Nips(8),8) = & ! element 117 - reshape(int([& - -3,-5,-4,-2,-6,-1 & - ],pInt),[FE_NipNeighbors(8),FE_Nips(8)]) - FE_ipNeighbor(1:FE_NipNeighbors(9),1:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape(int([& - 2,-5, 3,-2, 5,-1, & - -3, 1, 4,-2, 6,-1, & - 4,-5,-4, 1, 7,-1, & - -3, 3,-4, 2, 8,-1, & - 6,-5, 7,-2,-6, 1, & - -3, 5, 8,-2,-6, 2, & - 8,-5,-4, 5,-6, 3, & - -3, 7,-4, 6,-6, 4 & - ],pInt),[FE_NipNeighbors(9),FE_Nips(9)]) - - FE_ipNeighbor(1:FE_NipNeighbors(10),1:FE_Nips(10),10) = & ! element 155, 125, 128 - reshape(int([& - 2,-3, 3,-1, & - -2, 1, 3,-1, & - 2,-3,-2, 1 & - ],pInt),[FE_NipNeighbors(10),FE_Nips(10)]) - ! *** FE_subNodeParent *** ! lists the group of nodes for which the center of gravity ! corresponds to the location of a each subnode. @@ -4365,8 +4427,94 @@ FE_ipNeighbor(1:FE_NipNeighbors(8),1:FE_Nips(8),8) = & ! element 117 ! example: face-centered subnode with faceNodes 1,2,3,4 to be used in, ! e.g., a 8 IP grid, would be encoded: ! 1, 2, 3, 4, 0, 0, 0, 0 + me = 0_pInt - FE_subNodeParent(1:FE_Nips(1),1:FE_NsubNodes(1),1) = & ! element 7 + me = me + 1_pInt + FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 6 (2D 3node 1ip) has no subnodes + 0_pInt + + me = me + 1_pInt + FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 125 (2D 6node 3ip) + reshape(int([& + 1, 2, 0, & + 2, 3, 0, & + 3, 1, 0, & + 1, 2, 3 & + ],pInt),[FE_Nips(me),FE_NsubNodes(me)]) + + me = me + 1_pInt + FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 11 (2D 4node 4ip) + reshape(int([& + 1, 2, 0, 0, & + 2, 3, 0, 0, & + 3, 4, 0, 0, & + 4, 1, 0, 0, & + 1, 2, 3, 4 & + ],pInt),[FE_Nips(me),FE_NsubNodes(me)]) + + me = me + 1_pInt + FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 27 (2D 8node 9ip) + reshape(int([& + 1, 1, 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, 3, 3, 0, 0, 0, 0, 0, 0, & + 3, 3, 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, 1, 1, 0, 0, 0, 0, 0, 0, & + 1, 1, 1, 1, 2, 2, 4, 4, 3, & + 2, 2, 2, 2, 1, 1, 3, 3, 4, & + 3, 3, 3, 3, 2, 2, 4, 4, 1, & + 4, 4, 4, 4, 1, 1, 3, 3, 2 & + ],pInt),[FE_Nips(me),FE_NsubNodes(me)]) + + me = me + 1_pInt + FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 134 (3D 4node 1ip) has no subnodes + 0_pInt + + me = me + 1_pInt + FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 127 (3D 10node 4ip) + reshape(int([& + 1, 2, 0, 0, & + 2, 3, 0, 0, & + 3, 1, 0, 0, & + 1, 4, 0, 0, & + 2, 4, 0, 0, & + 3, 4, 0, 0, & + 1, 2, 3, 0, & + 1, 2, 4, 0, & + 2, 3, 4, 0, & + 1, 3, 4, 0, & + 1, 2, 3, 4 & + ],pInt),[FE_Nips(me),FE_NsubNodes(me)]) + + me = me + 1_pInt + FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 136 (3D 6node 6ip) + reshape(int([& + 1, 2, 0, 0, 0, 0, & + 2, 3, 0, 0, 0, 0, & + 3, 1, 0, 0, 0, 0, & + 1, 4, 0, 0, 0, 0, & + 2, 5, 0, 0, 0, 0, & + 3, 6, 0, 0, 0, 0, & + 4, 5, 0, 0, 0, 0, & + 5, 6, 0, 0, 0, 0, & + 6, 4, 0, 0, 0, 0, & + 1, 2, 3, 0, 0, 0, & + 1, 2, 4, 5, 0, 0, & + 2, 3, 5, 6, 0, 0, & + 1, 3, 4, 6, 0, 0, & + 4, 5, 6, 0, 0, 0, & + 1, 2, 3, 4, 5, 6 & + ],pInt),[FE_Nips(me),FE_NsubNodes(me)]) + + me = me + 1_pInt + FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 117 (3D 8node 1ip) has no subnodes + 0_pInt + + me = me + 1_pInt + FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 7 (3D 8node 8ip) reshape(int([& 1, 2, 0, 0, 0, 0, 0, 0, & 2, 3, 0, 0, 0, 0, 0, 0, & @@ -4387,60 +4535,10 @@ FE_ipNeighbor(1:FE_NipNeighbors(8),1:FE_Nips(8),8) = & ! element 117 1, 4, 8, 5, 0, 0, 0, 0, & 5, 6, 7, 8, 0, 0, 0, 0, & 1, 2, 3, 4, 5, 6, 7, 8 & - ],pInt),(/FE_Nips(1),FE_NsubNodes(1)/)) + ],pInt),[FE_Nips(me),FE_NsubNodes(me)]) -!FE_subNodeParent(:FE_Nips(2),:FE_NsubNodes(2),2) ! element 134 has no subnodes - - FE_subNodeParent(1:FE_Nips(3),1:FE_NsubNodes(3),3) = & ! element 11 - reshape(int([& - 1, 2, 0, 0, & - 2, 3, 0, 0, & - 3, 4, 0, 0, & - 4, 1, 0, 0, & - 1, 2, 3, 4 & - ],pInt),[FE_Nips(3),FE_NsubNodes(3)]) - - FE_subNodeParent(1:FE_Nips(4),1:FE_NsubNodes(4),4) = & ! element 27 - reshape(int([& - 1, 1, 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, 3, 3, 0, 0, 0, 0, 0, 0, & - 3, 3, 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, 1, 1, 0, 0, 0, 0, 0, 0, & - 1, 1, 1, 1, 2, 2, 4, 4, 3, & - 2, 2, 2, 2, 1, 1, 3, 3, 4, & - 3, 3, 3, 3, 2, 2, 4, 4, 1, & - 4, 4, 4, 4, 1, 1, 3, 3, 2 & - ],pInt),[FE_Nips(4),FE_NsubNodes(4)]) - - !FE_subNodeParent(:FE_Nips(5),:FE_NsubNodes(5),5) = & ! element 157 - ! reshape((/& - ! *still to be defined* - ! ],pInt),(/FE_Nips(5),FE_NsubNodes(5)/)) - - FE_subNodeParent(1:FE_Nips(6),1:FE_NsubNodes(6),6) = & ! element 136 - reshape(int([& - 1, 2, 0, 0, 0, 0, & - 2, 3, 0, 0, 0, 0, & - 3, 1, 0, 0, 0, 0, & - 1, 4, 0, 0, 0, 0, & - 2, 5, 0, 0, 0, 0, & - 3, 6, 0, 0, 0, 0, & - 4, 5, 0, 0, 0, 0, & - 5, 6, 0, 0, 0, 0, & - 6, 4, 0, 0, 0, 0, & - 1, 2, 3, 0, 0, 0, & - 1, 2, 4, 5, 0, 0, & - 2, 3, 5, 6, 0, 0, & - 1, 3, 4, 6, 0, 0, & - 4, 5, 6, 0, 0, 0, & - 1, 2, 3, 4, 5, 6 & - ],pInt),[FE_Nips(6),FE_NsubNodes(6)]) - - FE_subNodeParent(1:FE_Nips(7),1:FE_NsubNodes(7),7) = & ! element 21 + me = me + 1_pInt + FE_subNodeParent(1:FE_Nips(me),1:FE_NsubNodes(me),me) = & ! element 21 (3D 20node 27ip) reshape(int([& 1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 1, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & @@ -4498,41 +4596,9 @@ FE_ipNeighbor(1:FE_NipNeighbors(8),1:FE_Nips(8),8) = & ! element 117 6, 6, 6, 6, 6, 6, 6, 6, 2, 2, 2, 2, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 8, 8, 4, & 7, 7, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 5, 5, 1, & 8, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 6, 6, 2 & - ],pInt),[FE_Nips(7),FE_NsubNodes(7)]) + ],pInt),[FE_Nips(me),FE_NsubNodes(me)]) -!FE_subNodeParent(:FE_Nips(8),:FE_NsubNodes(8),8) ! element 117 has no subnodes - FE_subNodeParent(1:FE_Nips(9),1:FE_NsubNodes(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape(int([& - 1, 2, 0, 0, 0, 0, 0, 0, & - 2, 3, 0, 0, 0, 0, 0, 0, & - 3, 4, 0, 0, 0, 0, 0, 0, & - 4, 1, 0, 0, 0, 0, 0, 0, & - 1, 5, 0, 0, 0, 0, 0, 0, & - 2, 6, 0, 0, 0, 0, 0, 0, & - 3, 7, 0, 0, 0, 0, 0, 0, & - 4, 8, 0, 0, 0, 0, 0, 0, & - 5, 6, 0, 0, 0, 0, 0, 0, & - 6, 7, 0, 0, 0, 0, 0, 0, & - 7, 8, 0, 0, 0, 0, 0, 0, & - 8, 5, 0, 0, 0, 0, 0, 0, & - 1, 2, 3, 4, 0, 0, 0, 0, & - 1, 2, 6, 5, 0, 0, 0, 0, & - 2, 3, 7, 6, 0, 0, 0, 0, & - 3, 4, 8, 7, 0, 0, 0, 0, & - 1, 4, 8, 5, 0, 0, 0, 0, & - 5, 6, 7, 8, 0, 0, 0, 0, & - 1, 2, 3, 4, 5, 6, 7, 8 & - ],pInt),[FE_Nips(9),FE_NsubNodes(9)]) - - FE_subNodeParent(1:FE_Nips(10),1:FE_NsubNodes(10),10) = & ! element 155, 125, 128 - reshape(int([& - 1, 2, 0, & - 2, 3, 0, & - 3, 1, 0, & - 1, 2, 3 & - ],pInt),[FE_Nips(10),FE_NsubNodes(10)]) - ! *** FE_subNodeOnIPFace *** ! indicates which subnodes make up the interfaces enclosing the IP volume. ! The sorting convention is such that the outward pointing normal @@ -4543,8 +4609,187 @@ FE_ipNeighbor(1:FE_NipNeighbors(8),1:FE_Nips(8),8) = & ! element 117 ! This will result in zero ipVolume and interfaceArea, but is not ! detrimental at the moment since non-local constitutive laws are ! currently not foreseen in 2D cases. + me = 0_pInt - FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(1),1:FE_Nips(1),1) = & ! element 7 + me = me + 1_pInt + FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 6 (2D 3node 1ip) + reshape(int([& + 2, 3, 3, 2 , & ! 1 + 3, 1, 1, 3 , & + 1, 2, 2, 1 & + ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 125 (2D 6node 3ip) + reshape(int([& + 4, 7, 7, 4 , & ! 1 + 1, 6, 6, 1 , & + 6, 7, 7, 6 , & + 1, 4, 4, 1 , & + 2, 5, 5, 2 , & ! 2 + 4, 7, 7, 4 , & + 5, 7, 7, 5 , & + 2, 4, 4, 2 , & + 5, 7, 7, 5 , & ! 3 + 3, 6, 6, 3 , & + 3, 5, 5, 3 , & + 6, 7, 7, 6 & + ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 11 (2D 4node 4ip) + reshape(int([& + 5, 9, 9, 5 , & ! 1 + 1, 8, 8, 1 , & + 8, 9, 9, 8 , & + 1, 5, 5, 1 , & + 2, 6, 6, 2 , & ! 2 + 5, 9, 9, 5 , & + 6, 9, 9, 6 , & + 2, 5, 5, 2 , & + 3, 6, 6, 3 , & ! 3 + 7, 9, 9, 7 , & + 3, 7, 7, 3 , & + 6, 9, 9, 6 , & + 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)]) + + me = me + 1_pInt + FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 27 (2D 8node 9ip) + reshape(int([& + 9,17,17, 9 , & ! 1 + 1,16,16, 1 , & + 16,17,17,16 , & + 1, 9, 9, 1 , & + 10,18,18,10 , & ! 2 + 9,17,17, 9 , & + 17,18,18,17 , & + 9,10,10, 9 , & + 2,11,11, 2 , & ! 3 + 10,18,18,10 , & + 11,18,18,11 , & + 2,10,10, 2 , & + 17,20,20,17 , & ! 4 + 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 , & + 19,20,20,19 , & + 3,12,12, 3 , & ! 9 + 13,19,19,13 , & + 3,13,13, 3 , & + 12,19,19,12 & + ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 134 (3D 4node 1ip) + reshape(int([& + 1, 1, 3, 2, & ! 1 + 1, 1, 2, 4, & + 2, 2, 3, 4, & + 1, 1, 4, 3 & + ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 127 (3D 10node 4ip) + reshape(int([& + 5,11,15,12 , & ! 1 + 1, 8,14, 7 , & + 7,14,15,11 , & + 1, 5,12, 8 , & + 8,12,15,14 , & + 1, 7,11, 5 , & + 2, 6,13, 9 , & ! 2 + 5,12,15,11 , & + 6,11,15,13 , & + 2, 9,12, 5 , & + 9,13,15,12 , & + 2,13,11, 6 , & + 6,13,15,11 , & ! 3 + 3, 7,14,10 , & + 3,10,13, 6 , & + 7,11,15,14 , & + 13,10,14,15 , & + 3, 6,11, 7 , & + 9,12,15,13 , & ! 4 + 4,10,14, 8 , & + 10,13,15,14 , & + 4, 8,12, 9 , & + 4, 9,13,10 , & + 8,10,15,12 & + ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 136 (3D 6node 6ip) + reshape(int([& + 7,16,21,17, & ! 1 + 1,10,19, 9, & + 9,19,21,16, & + 1, 7,17,10, & + 10,17,21,19, & + 1, 9,16, 7, & + 2, 8,18,11, & ! 2 + 7,17,21,16, & + 8,16,21,18, & + 2,11,17, 7, & + 11,18,21,17, & + 2, 7,16, 8, & + 8,18,21,16, & ! 3 + 3, 9,19,12, & + 3,12,18, 8, & + 9,16,21,19, & + 12,19,21,18, & + 3, 8,16, 9, & + 13,17,21,20, & ! 4 + 4,15,19,10, & + 15,20,21,19, & + 4,10,17,13, & + 4,13,20,15, & + 10,19,21,17, & + 5,11,18,14, & ! 5 + 13,20,21,17, & + 14,18,21,20, & + 5,13,17,11, & + 5,14,20,13, & + 11,17,21,18, & + 14,20,21,18, & ! 6 + 6,12,19,15, & + 6,14,18,12, & + 15,19,21,20, & + 6,15,20,14, & + 12,18,21,19 & + ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 117 (3D 8node 1ip) + reshape(int([& + 2, 3, 7, 6, & ! 1 + 1, 5, 8, 4, & + 3, 4, 8, 7, & + 1, 2, 6, 5, & + 5, 6, 7, 8, & + 1, 4, 3, 2 & + ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) + + me = me + 1_pInt + FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 7 (3D 8node 8ip) reshape(int([& 9,21,27,22, & ! 1 1,13,25,12, & @@ -4594,122 +4839,10 @@ FE_ipNeighbor(1:FE_NipNeighbors(8),1:FE_Nips(8),8) = & ! element 117 18,26,27,23, & 7,19,26,18, & 15,23,27,24 & - ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(1),FE_Nips(1)]) + ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) - FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(2),1:FE_Nips(2),2) = & ! element 134 - reshape(int([& - 1, 1, 3, 2, & ! 1 - 1, 1, 2, 4, & - 2, 2, 3, 4, & - 1, 1, 4, 3 & - ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(2),FE_Nips(2)]) - - FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(3),1:FE_Nips(3),3) = & ! element 11 - reshape(int([& - 5, 9, 9, 5 , & ! 1 - 1, 8, 8, 1 , & - 8, 9, 9, 8 , & - 1, 5, 5, 1 , & - 2, 6, 6, 2 , & ! 2 - 5, 9, 9, 5 , & - 6, 9, 9, 6 , & - 2, 5, 5, 2 , & - 3, 6, 6, 3 , & ! 3 - 7, 9, 9, 7 , & - 3, 7, 7, 3 , & - 6, 9, 9, 6 , & - 7, 9, 9, 7 , & ! 4 - 4, 8, 8, 4 , & - 4, 7, 7, 4 , & - 8, 9, 9, 8 & - ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(3),FE_Nips(3)]) - - FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(4),1:FE_Nips(4),4) = & ! element 27 - reshape(int([& - 9,17,17, 9 , & ! 1 - 1,16,16, 1 , & - 16,17,17,16 , & - 1, 9, 9, 1 , & - 10,18,18,10 , & ! 2 - 9,17,17, 9 , & - 17,18,18,17 , & - 9,10,10, 9 , & - 2,11,11, 2 , & ! 3 - 10,18,18,10 , & - 11,18,18,11 , & - 2,10,10, 2 , & - 17,20,20,17 , & ! 4 - 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 , & - 19,20,20,19 , & - 3,12,12, 3 , & ! 9 - 13,19,19,13 , & - 3,13,13, 3 , & - 12,19,19,12 & - ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(4),FE_Nips(4)]) - - !FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(5),:FE_Nips(5),5) = & ! element 157 - ! reshape((/& - ! *still to be defined* - ! /),(/FE_NipFaceNodes,FE_NipNeighbors(5),FE_Nips(5)/)) - - FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(6),1:FE_Nips(6),6) = & ! element 136 - reshape(int([& - 7,16,21,17, & ! 1 - 1,10,19, 9, & - 9,19,21,16, & - 1, 7,17,10, & - 10,17,21,19, & - 1, 9,16, 7, & - 2, 8,18,11, & ! 2 - 7,17,21,16, & - 8,16,21,18, & - 2,11,17, 7, & - 11,18,21,17, & - 2, 7,16, 8, & - 8,18,21,16, & ! 3 - 3, 9,19,12, & - 3,12,18, 8, & - 9,16,21,19, & - 12,19,21,18, & - 3, 8,16, 9, & - 13,17,21,20, & ! 4 - 4,15,19,10, & - 15,20,21,19, & - 4,10,17,13, & - 4,13,20,15, & - 10,19,21,17, & - 5,11,18,14, & ! 5 - 13,20,21,17, & - 14,18,21,20, & - 5,13,17,11, & - 5,14,20,13, & - 11,17,21,18, & - 14,20,21,18, & ! 6 - 6,12,19,15, & - 6,14,18,12, & - 15,19,21,20, & - 6,15,20,14, & - 12,18,21,19 & - ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(6),FE_Nips(6)]) - - FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(7),1:FE_Nips(7),7) = & ! element 21 + me = me + 1_pInt + FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(me),1:FE_Nips(me),me) = & ! element 21 (3D 20node 27ip) reshape(int([& 9,33,57,37, & ! 1 1,17,44,16, & @@ -4873,85 +5006,8 @@ FE_ipNeighbor(1:FE_NipNeighbors(8),1:FE_Nips(8),8) = & ! element 117 63,48,28,55, & 55,28, 7,29, & 63,49,23,48 & - ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(7),FE_Nips(7)]) + ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(me),FE_Nips(me)]) - FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(8),1:FE_Nips(8),8) = & ! element 117 - reshape(int([& - 2, 3, 7, 6, & ! 1 - 1, 5, 8, 4, & - 3, 4, 8, 7, & - 1, 2, 6, 5, & - 5, 6, 7, 8, & - 1, 4, 3, 2 & - ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(8),FE_Nips(8)]) - - FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(9),1:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape(int([& - 9,21,27,22, & ! 1 - 1,13,25,12, & - 12,25,27,21, & - 1, 9,22,13, & - 13,22,27,25, & - 1,12,21, 9, & - 2,10,23,14, & ! 2 - 9,22,27,21, & - 10,21,27,23, & - 2,14,22, 9, & - 14,23,27,22, & - 2, 9,21,10, & - 11,24,27,21, & ! 3 - 4,12,25,16, & - 4,16,24,11, & - 12,21,27,25, & - 16,25,27,24, & - 4,11,21,12, & - 3,15,23,10, & ! 4 - 11,21,27,24, & - 3,11,24,15, & - 10,23,27,21, & - 15,24,27,23, & - 3,10,21,11, & - 17,22,27,26, & ! 5 - 5,20,25,13, & - 20,26,27,25, & - 5,13,22,17, & - 5,17,26,20, & - 13,25,27,22, & - 6,14,23,18, & ! 6 - 17,26,27,22, & - 18,23,27,26, & - 6,17,22,14, & - 6,18,26,17, & - 14,22,27,23, & - 19,26,27,24, & ! 7 - 8,16,25,20, & - 8,19,24,16, & - 20,25,27,26, & - 8,20,26,19, & - 16,24,27,25, & - 7,18,23,15, & ! 8 - 19,24,27,26, & - 7,15,24,19, & - 18,26,27,23, & - 7,19,26,18, & - 15,23,27,24 & - ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(9),FE_Nips(9)]) - - FE_subNodeOnIPFace(1:FE_NipFaceNodes,1:FE_NipNeighbors(10),1:FE_Nips(10),10) = & ! element 155, 125, 128 - reshape(int([& - 4, 7, 7, 4 , & ! 1 - 1, 6, 6, 1 , & - 6, 7, 7, 6 , & - 1, 4, 4, 1 , & - 2, 5, 5, 2 , & ! 2 - 4, 7, 7, 4 , & - 5, 7, 7, 5 , & - 2, 4, 4, 2 , & - 5, 7, 7, 5 , & ! 3 - 3, 6, 6, 3 , & - 3, 5, 5, 3 , & - 6, 7, 7, 6 & - ],pInt),[FE_NipFaceNodes,FE_NipNeighbors(10),FE_Nips(10)]) end subroutine mesh_build_FEdata