2011-04-04 19:39:54 +05:30
! Copyright 2011 Max-Planck-Institut f<> r Eisenforschung GmbH
!
! This file is part of DAMASK,
! 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
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!##############################################################
2009-08-31 20:39:15 +05:30
!* $Id$
2007-03-20 19:25:22 +05:30
!##############################################################
2007-03-21 21:48:33 +05:30
MODULE mesh
2007-03-20 19:25:22 +05:30
!##############################################################
2007-03-21 21:48:33 +05:30
use prec , only : pReal , pInt
implicit none
! ---------------------------
2007-03-22 20:18:58 +05:30
! _Nelems : total number of elements in mesh
2007-04-25 13:03:24 +05:30
! _NcpElems : total number of CP elements in mesh
2007-03-22 20:18:58 +05:30
! _Nnodes : total number of nodes in mesh
2007-04-04 14:17:34 +05:30
! _maxNnodes : max number of nodes in any CP element
2007-04-25 13:03:24 +05:30
! _maxNips : max number of IPs in any CP element
! _maxNipNeighbors : max number of IP neighbors in any CP element
2007-04-04 14:17:34 +05:30
! _maxNsharedElems : max number of CP elements sharing a node
2007-03-22 20:18:58 +05:30
!
2008-03-25 18:22:27 +05:30
! _element : FEid, type(internal representation), material, texture, node indices
2011-05-24 21:27:59 +05:30
! _node0 : x,y,z coordinates (initially!)
! _node : x,y,z coordinates (after deformation!)
2007-03-26 14:20:57 +05:30
! _sharedElem : entryCount and list of elements containing node
2007-03-22 20:18:58 +05:30
!
2007-03-27 18:23:31 +05:30
! _mapFEtoCPelem : [sorted FEid, corresponding CPid]
! _mapFEtoCPnode : [sorted FEid, corresponding CPid]
2007-03-22 20:18:58 +05:30
!
! MISSING: these definitions should actually reside in the
! FE-solver specific part (different for MARC/ABAQUS)..!
! Hence, I suggest to prefix with "FE_"
!
2009-04-03 12:33:25 +05:30
! _Nnodes : # nodes in a specific type of element (how we use it)
! _NoriginalNodes : # nodes in a specific type of element (how it is originally defined by marc)
2007-03-26 14:20:57 +05:30
! _Nips : # IPs in a specific type of element
! _NipNeighbors : # IP neighbors in a specific type of element
! _ipNeighbor : +x,-x,+y,-y,+z,-z list of intra-element IPs and
2007-03-22 20:18:58 +05:30
! (negative) neighbor faces per own IP in a specific type of element
2009-01-16 23:06:37 +05:30
! _NfaceNodes : # nodes per face in a specific type of element
2007-03-26 14:20:57 +05:30
! _nodeOnFace : list of node indices on each face of a specific type of element
2010-05-10 20:24:59 +05:30
! _maxNnodesAtIP : max number of (equivalent) nodes attached to an IP
2009-03-30 20:20:19 +05:30
! _nodesAtIP : map IP index to two node indices in a specific type of element
2009-01-16 23:06:37 +05:30
! _NsubNodes : # subnodes required to fully define all IP volumes
2007-03-22 20:18:58 +05:30
! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype
2007-03-21 21:48:33 +05:30
! ---------------------------
2011-02-16 21:53:08 +05:30
integer ( pInt ) mesh_Nelems , mesh_NcpElems , mesh_NelemSets , mesh_maxNelemInSet
2009-10-12 21:31:49 +05:30
integer ( pInt ) mesh_Nmaterials
2011-02-16 21:53:08 +05:30
integer ( pInt ) mesh_Nnodes , mesh_maxNnodes , mesh_maxNips , mesh_maxNipNeighbors , mesh_maxNsharedElems , mesh_maxNsubNodes
2007-10-24 14:30:42 +05:30
integer ( pInt ) , dimension ( 2 ) :: mesh_maxValStateVar = 0_pInt
2009-10-12 21:31:49 +05:30
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
2009-08-11 22:01:57 +05:30
integer ( pInt ) , dimension ( : , : ) , allocatable , target :: mesh_mapFEtoCPelem , mesh_mapFEtoCPnode
2011-05-24 21:27:59 +05:30
integer ( pInt ) , dimension ( : , : ) , allocatable :: mesh_element , & ! FEid, type(internal representation), material, texture, node indices
mesh_sharedElem , & ! entryCount and list of elements containing node
2011-02-16 21:53:08 +05:30
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)
2011-05-24 21:27:59 +05:30
integer ( pInt ) , dimension ( : , : , : , : ) , allocatable :: mesh_ipNeighborhood ! 6 or less neighboring IPs as [element_num, IP_index]
2009-01-20 00:12:31 +05:30
2009-08-11 22:01:57 +05:30
real ( pReal ) , dimension ( : , : , : ) , allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element
2011-05-24 21:27:59 +05:30
real ( pReal ) , dimension ( : , : ) , allocatable :: mesh_node0 , & ! node coordinates (initially!)
mesh_node , & ! node coordinates (after deformation! ONLY FOR MARC!!!)
mesh_ipVolume ! volume associated with IP (initially!)
real ( pReal ) , dimension ( : , : , : ) , allocatable :: mesh_ipArea , & ! area of interface to neighboring IP (initially!)
mesh_ipCenterOfGravity ! center of gravity of IP (after deformation!)
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP (initially!)
integer ( pInt ) , dimension ( : , : , : ) , allocatable :: FE_nodesAtIP ! map IP index to two node indices in a specific type of element
2009-04-06 18:55:19 +05:30
integer ( pInt ) , dimension ( : , : , : ) , allocatable :: FE_ipNeighbor
integer ( pInt ) , dimension ( : , : , : ) , allocatable :: FE_subNodeParent
integer ( pInt ) , dimension ( : , : , : , : ) , allocatable :: FE_subNodeOnIPFace
2010-07-13 15:56:07 +05:30
logical :: noPart ! for cases where the ABAQUS input file does not use part/assembly information
2011-02-16 21:53:08 +05:30
logical , dimension ( 3 ) :: mesh_periodicSurface ! flag indicating periodic outer surfaces (used for fluxes)
2007-03-21 21:48:33 +05:30
2009-10-12 21:31:49 +05:30
integer ( pInt ) :: hypoelasticTableStyle
integer ( pInt ) :: initialcondTableStyle
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: FE_Nelemtypes = 10_pInt
integer ( pInt ) , parameter :: FE_maxNnodes = 8_pInt
integer ( pInt ) , parameter :: FE_maxNsubNodes = 56_pInt
integer ( pInt ) , parameter :: FE_maxNips = 27_pInt
integer ( pInt ) , parameter :: FE_maxNipNeighbors = 6_pInt
integer ( pInt ) , parameter :: FE_maxmaxNnodesAtIP = 8_pInt
integer ( pInt ) , parameter :: FE_NipFaceNodes = 4_pInt
2007-03-26 14:20:57 +05:30
integer ( pInt ) , dimension ( FE_Nelemtypes ) , parameter :: FE_Nnodes = &
2012-02-16 00:28:38 +05:30
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 )
2009-04-03 12:33:25 +05:30
integer ( pInt ) , dimension ( FE_Nelemtypes ) , parameter :: FE_NoriginalNodes = &
2012-02-16 00:28:38 +05:30
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
] , pInt )
2007-03-26 14:20:57 +05:30
integer ( pInt ) , dimension ( FE_Nelemtypes ) , parameter :: FE_Nips = &
2012-02-16 00:28:38 +05:30
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 )
2009-01-16 23:06:37 +05:30
integer ( pInt ) , dimension ( FE_Nelemtypes ) , parameter :: FE_NipNeighbors = &
2012-02-16 00:28:38 +05:30
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 )
2009-01-16 23:06:37 +05:30
integer ( pInt ) , dimension ( FE_Nelemtypes ) , parameter :: FE_NsubNodes = &
2012-02-16 00:28:38 +05:30
int ( [ 19 , & ! element 7
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
] , pInt )
2009-01-16 20:59:57 +05:30
integer ( pInt ) , dimension ( FE_maxNipNeighbors , FE_Nelemtypes ) , parameter :: FE_NfaceNodes = &
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2007-04-25 13:03:24 +05:30
4 , 4 , 4 , 4 , 4 , 4 , & ! element 7
2007-10-12 19:18:29 +05:30
3 , 3 , 3 , 3 , 0 , 0 , & ! element 134
2 , 2 , 2 , 2 , 0 , 0 , & ! element 11
2009-04-02 18:30:51 +05:30
2 , 2 , 2 , 2 , 0 , 0 , & ! element 27
2008-06-17 14:41:54 +05:30
3 , 3 , 3 , 3 , 0 , 0 , & ! element 157
2009-04-02 18:30:51 +05:30
3 , 4 , 4 , 4 , 3 , 0 , & ! element 136
2010-05-10 20:24:59 +05:30
4 , 4 , 4 , 4 , 4 , 4 , & ! element 21
2010-08-17 04:23:24 +05:30
4 , 4 , 4 , 4 , 4 , 4 , & ! element 117
2011-04-06 14:05:37 +05:30
4 , 4 , 4 , 4 , 4 , 4 , & ! element 57 (c3d20r == c3d8 --> copy of 7)
2 , 2 , 2 , 0 , 0 , 0 & ! element 155, 125, 128
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_maxNipNeighbors , FE_Nelemtypes / ) )
2010-05-10 20:24:59 +05:30
integer ( pInt ) , dimension ( FE_Nelemtypes ) , parameter :: FE_maxNnodesAtIP = &
2012-02-16 00:28:38 +05:30
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
] , pInt )
2009-01-16 20:59:57 +05:30
integer ( pInt ) , dimension ( FE_NipFaceNodes , FE_maxNipNeighbors , FE_Nelemtypes ) , parameter :: FE_nodeOnFace = &
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2007-04-25 13:03:24 +05:30
1 , 2 , 3 , 4 , & ! element 7
2 , 1 , 5 , 6 , &
3 , 2 , 6 , 7 , &
2008-06-17 14:41:54 +05:30
4 , 3 , 7 , 8 , &
2007-04-25 13:03:24 +05:30
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 , &
0 , 0 , 0 , 0 , &
2007-10-12 19:18:29 +05:30
0 , 0 , 0 , 0 , &
1 , 2 , 0 , 0 , & ! element 11
2 , 3 , 0 , 0 , &
3 , 4 , 0 , 0 , &
4 , 1 , 0 , 0 , &
0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , &
2009-04-02 18:30:51 +05:30
1 , 2 , 0 , 0 , & ! element 27
2 , 3 , 0 , 0 , &
3 , 4 , 0 , 0 , &
4 , 1 , 0 , 0 , &
2007-10-12 19:18:29 +05:30
0 , 0 , 0 , 0 , &
2008-06-17 14:41:54 +05:30
0 , 0 , 0 , 0 , &
1 , 2 , 3 , 0 , & ! element 157
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 , 4 , 5 , 2 , &
2 , 5 , 6 , 3 , &
1 , 3 , 6 , 4 , &
4 , 6 , 5 , 0 , &
2009-04-02 18:30:51 +05:30
0 , 0 , 0 , 0 , &
1 , 2 , 3 , 4 , & ! element 21
2 , 1 , 5 , 6 , &
3 , 2 , 6 , 7 , &
4 , 3 , 7 , 8 , &
4 , 1 , 5 , 8 , &
2010-05-10 20:24:59 +05:30
8 , 7 , 6 , 5 , &
1 , 2 , 3 , 4 , & ! element 117
2 , 1 , 5 , 6 , &
3 , 2 , 6 , 7 , &
4 , 3 , 7 , 8 , &
4 , 1 , 5 , 8 , &
2010-08-17 04:23:24 +05:30
8 , 7 , 6 , 5 , &
1 , 2 , 3 , 4 , & ! element 57 (c3d20r == c3d8 --> copy of 7)
2 , 1 , 5 , 6 , &
3 , 2 , 6 , 7 , &
4 , 3 , 7 , 8 , &
4 , 1 , 5 , 8 , &
2011-04-06 14:05:37 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipFaceNodes , FE_maxNipNeighbors , FE_Nelemtypes / ) )
2009-04-06 18:55:19 +05:30
2007-03-21 21:48:33 +05:30
CONTAINS
! ---------------------------
! subroutine mesh_init()
! function mesh_FEtoCPelement(FEid)
2007-03-26 14:20:57 +05:30
! function mesh_build_ipNeighorhood()
2007-03-21 21:48:33 +05:30
! ---------------------------
2009-01-20 00:12:31 +05:30
2007-04-25 13:03:24 +05:30
!***********************************************************
! initialization
!***********************************************************
2009-10-12 21:31:49 +05:30
subroutine mesh_init ( ip , element )
2009-01-20 00:12:31 +05:30
2011-05-11 22:31:03 +05:30
use DAMASK_interface
2012-02-16 00:28:38 +05:30
use , intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
2007-04-10 16:52:53 +05:30
use prec , only : pInt
2010-07-13 15:56:07 +05:30
use IO , only : IO_error , IO_open_InputFile , IO_abaqus_hasNoPart
2011-05-28 15:14:43 +05:30
use FEsolving , only : parallelExecution , FEsolving_execElem , FEsolving_execIP , calcMode , lastMode , FEmodelGeometry
2009-03-04 17:18:54 +05:30
2007-04-10 16:52:53 +05:30
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: fileUnit = 222_pInt
2009-10-12 21:31:49 +05:30
integer ( pInt ) e , element , ip
2009-06-18 19:58:02 +05:30
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2011-03-21 16:01:17 +05:30
write ( 6 , * )
write ( 6 , * ) '<<<+- mesh init -+>>>'
write ( 6 , * ) '$Id$'
2012-02-01 00:48:55 +05:30
#include "compilation_info.f90"
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-01-20 00:12:31 +05:30
2011-05-28 15:14:43 +05:30
call mesh_build_FEdata ( ) ! --- get properties of the different types of elements
2009-10-12 21:31:49 +05:30
2012-02-13 23:11:27 +05:30
call IO_open_inputFile ( fileUnit , FEmodelGeometry ) ! --- parse info from input file...
select case ( FEsolver )
case ( 'Spectral' )
call mesh_spectral_count_nodesAndElements ( fileUnit )
call mesh_spectral_count_cpElements ( )
call mesh_spectral_map_elements ( )
call mesh_spectral_map_nodes ( )
call mesh_spectral_count_cpSizes ( )
call mesh_spectral_build_nodes ( fileUnit )
call mesh_spectral_build_elements ( fileUnit )
2009-10-12 21:31:49 +05:30
2012-02-13 23:11:27 +05:30
case ( 'Marc' )
call mesh_marc_get_tableStyles ( fileUnit )
call mesh_marc_count_nodesAndElements ( fileUnit )
call mesh_marc_count_elementSets ( fileUnit )
call mesh_marc_map_elementSets ( fileUnit )
call mesh_marc_count_cpElements ( fileUnit )
call mesh_marc_map_elements ( fileUnit )
call mesh_marc_map_nodes ( fileUnit )
call mesh_marc_build_nodes ( fileUnit )
call mesh_marc_count_cpSizes ( fileunit )
call mesh_marc_build_elements ( fileUnit )
case ( 'Abaqus' )
noPart = IO_abaqus_hasNoPart ( fileUnit )
call mesh_abaqus_count_nodesAndElements ( fileUnit )
call mesh_abaqus_count_elementSets ( fileUnit )
call mesh_abaqus_count_materials ( fileUnit )
call mesh_abaqus_map_elementSets ( fileUnit )
call mesh_abaqus_map_materials ( fileUnit )
call mesh_abaqus_count_cpElements ( fileUnit )
call mesh_abaqus_map_elements ( fileUnit )
call mesh_abaqus_map_nodes ( fileUnit )
call mesh_abaqus_build_nodes ( fileUnit )
call mesh_abaqus_count_cpSizes ( fileunit )
call mesh_abaqus_build_elements ( fileUnit )
end select
call mesh_get_damaskOptions ( fileUnit )
close ( fileUnit )
call mesh_build_subNodeCoords ( )
call mesh_build_ipCoordinates ( )
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
2007-04-10 16:52:53 +05:30
2012-02-10 17:26:05 +05:30
FEsolving_execElem = [ 1_pInt , mesh_NcpElems ]
2012-02-16 00:28:38 +05:30
allocate ( FEsolving_execIP ( 2_pInt , mesh_NcpElems ) ) ; FEsolving_execIP = 1_pInt
2012-02-10 17:26:05 +05:30
forall ( e = 1_pInt : mesh_NcpElems ) FEsolving_execIP ( 2 , e ) = FE_Nips ( mesh_element ( 2 , e ) )
2009-06-15 18:41:21 +05:30
2009-10-12 21:31:49 +05:30
allocate ( calcMode ( mesh_maxNips , mesh_NcpElems ) )
calcMode = . false . ! pretend to have collected what first call is asking (F = I)
calcMode ( ip , mesh_FEasCP ( 'elem' , element ) ) = . true . ! first ip,el needs to be already pingponged to "calc"
lastMode = . true . ! and its mode is already known...
2009-06-15 18:41:21 +05:30
endsubroutine
2007-03-26 14:20:57 +05:30
2007-04-10 16:52:53 +05:30
2009-01-20 00:12:31 +05:30
2008-06-17 02:19:48 +05:30
!***********************************************************
! mapping of FE element types to internal representation
!***********************************************************
2009-06-15 18:41:21 +05:30
function FE_mapElemtype ( what )
2010-07-13 15:56:07 +05:30
use IO , only : IO_lc
2009-01-20 00:12:31 +05:30
2008-06-17 02:19:48 +05:30
implicit none
character ( len = * ) , intent ( in ) :: what
integer ( pInt ) FE_mapElemtype
2010-07-13 15:56:07 +05:30
select case ( IO_lc ( what ) )
2010-05-10 20:24:59 +05:30
case ( '7' , &
2010-07-13 15:56:07 +05:30
'c3d8' )
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 1_pInt ! Three-dimensional Arbitrarily Distorted Brick
2010-05-10 20:24:59 +05:30
case ( '134' , &
2010-07-13 15:56:07 +05:30
'c3d4' )
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 2_pInt ! Three-dimensional Four-node Tetrahedron
2010-05-10 20:24:59 +05:30
case ( '11' , &
2010-07-13 15:56:07 +05:30
'cpe4' )
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 3_pInt ! Arbitrary Quadrilateral Plane-strain
2010-05-10 20:24:59 +05:30
case ( '27' , &
2010-07-13 15:56:07 +05:30
'cpe8' )
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 4_pInt ! Plane Strain, Eight-node Distorted Quadrilateral
2008-06-17 14:41:54 +05:30
case ( '157' )
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 5_pInt ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations
2010-05-10 20:24:59 +05:30
case ( '136' , &
2010-07-13 15:56:07 +05:30
'c3d6' )
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 6_pInt ! Three-dimensional Arbitrarily Distorted Pentahedral
2010-05-10 20:24:59 +05:30
case ( '21' , &
2010-07-13 15:56:07 +05:30
'c3d20' )
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 7_pInt ! Three-dimensional Arbitrarily Distorted quadratic hexahedral
2010-05-10 20:24:59 +05:30
case ( '117' , &
'123' , &
2010-07-13 15:56:07 +05:30
'c3d8r' )
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 8_pInt ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration
2010-08-17 04:23:24 +05:30
case ( '57' , &
'c3d20r' )
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 9_pInt ! Three-dimensional Arbitrarily Distorted quad hexahedral with reduced integration
2011-04-06 14:05:37 +05:30
case ( '155' , &
'125' , &
'128' )
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 10_pInt ! Two-dimensional Plane Strain triangle (155: cubic shape function, 125/128: second order isoparametric)
2009-10-12 21:31:49 +05:30
case default
2012-02-10 17:26:05 +05:30
FE_mapElemtype = 0_pInt ! unknown element --> should raise an error upstream..!
2009-06-15 18:41:21 +05:30
endselect
2009-01-20 00:12:31 +05:30
2009-06-15 18:41:21 +05:30
endfunction
2008-06-17 02:19:48 +05:30
2009-01-20 00:12:31 +05:30
2007-03-26 14:20:57 +05:30
!***********************************************************
2007-04-10 16:52:53 +05:30
! FE to CP id mapping by binary search thru lookup array
2007-03-26 14:20:57 +05:30
!
2007-04-10 16:52:53 +05:30
! valid questions are 'elem', 'node'
!***********************************************************
2009-06-15 18:41:21 +05:30
function mesh_FEasCP ( what , id )
2007-04-10 16:52:53 +05:30
use prec , only : pInt
use IO , only : IO_lc
implicit none
character ( len = * ) , intent ( in ) :: what
integer ( pInt ) , intent ( in ) :: id
integer ( pInt ) , dimension ( : , : ) , pointer :: lookupMap
integer ( pInt ) mesh_FEasCP , lower , upper , center
mesh_FEasCP = 0_pInt
select case ( IO_lc ( what ( 1 : 4 ) ) )
case ( 'elem' )
lookupMap = > mesh_mapFEtoCPelem
case ( 'node' )
lookupMap = > mesh_mapFEtoCPnode
case default
return
2009-06-15 18:41:21 +05:30
endselect
2007-04-10 16:52:53 +05:30
lower = 1_pInt
2012-02-10 17:26:05 +05:30
upper = int ( size ( lookupMap , 2_pInt ) , pInt )
2007-04-10 16:52:53 +05:30
! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds?
2012-02-16 00:28:38 +05:30
if ( lookupMap ( 1_pInt , lower ) == id ) then
mesh_FEasCP = lookupMap ( 2_pInt , lower )
2007-04-10 16:52:53 +05:30
return
2012-02-16 00:28:38 +05:30
elseif ( lookupMap ( 1_pInt , upper ) == id ) then
mesh_FEasCP = lookupMap ( 2_pInt , upper )
2007-04-10 16:52:53 +05:30
return
endif
! binary search in between bounds
2012-02-16 00:28:38 +05:30
do while ( upper - lower > 1_pInt )
center = ( lower + upper ) / 2_pInt
if ( lookupMap ( 1_pInt , center ) < id ) then
2007-04-10 16:52:53 +05:30
lower = center
2012-02-16 00:28:38 +05:30
elseif ( lookupMap ( 1_pInt , center ) > id ) then
2007-04-10 16:52:53 +05:30
upper = center
else
2012-02-16 00:28:38 +05:30
mesh_FEasCP = lookupMap ( 2_pInt , center )
2007-04-10 16:52:53 +05:30
exit
2009-06-15 18:41:21 +05:30
endif
enddo
2007-04-10 16:52:53 +05:30
2009-06-15 18:41:21 +05:30
endfunction
2007-04-10 16:52:53 +05:30
2009-01-20 00:12:31 +05:30
2007-03-26 14:20:57 +05:30
!***********************************************************
2007-04-10 16:52:53 +05:30
! find face-matching element of same type
2009-10-12 21:31:49 +05:30
!***********************************************************
2011-02-16 21:53:08 +05:30
subroutine mesh_faceMatch ( elem , face , matchingElem , matchingFace )
use prec , only : pInt
implicit none
!*** output variables
integer ( pInt ) , intent ( out ) :: matchingElem , & ! matching CP element ID
matchingFace ! matching FE face ID
!*** input variables
integer ( pInt ) , intent ( in ) :: face , & ! FE face ID
elem ! FE elem ID
!*** local variables
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"
2012-02-16 00:28:38 +05:30
myType = mesh_element ( 2_pInt , elem ) ! figure elemType
2011-02-16 21:53:08 +05:30
2012-02-16 00:28:38 +05:30
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
NsharedElems = mesh_sharedElem ( 1_pInt , myFaceNodes ( n ) ) ! figure # shared elements for this node
2011-02-16 21:53:08 +05:30
if ( NsharedElems < minNsharedElems ) then
minNsharedElems = NsharedElems ! remember min # shared elems
lonelyNode = n ! remember most lonely node
endif
enddo
allocate ( element_seen ( minNsharedElems ) )
element_seen = 0_pInt
2012-02-16 00:28:38 +05:30
checkCandidate : do i = 1_pInt , minNsharedElems ! iterate over lonelyNode's shared elements
candidateElem = mesh_sharedElem ( 1_pInt + i , myFaceNodes ( lonelyNode ) ) ! present candidate elem
2011-02-16 21:53:08 +05:30
if ( all ( element_seen / = candidateElem ) ) then ! element seen for the first time?
element_seen ( i ) = candidateElem
2012-02-16 00:28:38 +05:30
candidateType = mesh_element ( 2_pInt , candidateElem ) ! figure elemType of candidate
checkCandidateFace : do candidateFace = 1_pInt , FE_maxNipNeighbors ! check each face of candidate
2011-02-16 21:53:08 +05:30
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 .
2012-02-16 00:28:38 +05:30
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 ) )
2011-02-16 21:53:08 +05:30
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
if ( checkTwins ) then
2012-02-16 00:28:38 +05:30
checkCandidateFaceTwins : do dir = 1_pInt , 3_pInt
do n = 1_pInt , FE_NfaceNodes ( candidateFace , candidateType ) ! loop through nodes on face
2011-02-16 21:53:08 +05:30
candidateFaceNode = mesh_FEasCP ( 'node' , mesh_element ( 4 + FE_nodeOnFace ( n , candidateFace , candidateType ) , candidateElem ) )
if ( all ( myFaceNodes / = mesh_nodeTwins ( dir , candidateFaceNode ) ) ) then ! node twin does not match either
2012-02-16 00:28:38 +05:30
if ( dir == 3_pInt ) then
2011-02-16 21:53:08 +05:30
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
deallocate ( element_seen )
endsubroutine
2007-03-26 14:20:57 +05:30
2009-04-06 18:55:19 +05:30
!********************************************************************
! get properties of different types of finite elements
!
! assign globals:
! FE_nodesAtIP, FE_ipNeighbor, FE_subNodeParent, FE_subNodeOnIPFace
!********************************************************************
2009-10-12 21:31:49 +05:30
subroutine mesh_build_FEdata ( )
2009-04-06 18:55:19 +05:30
use prec , only : pInt
implicit none
2010-05-10 20:24:59 +05:30
allocate ( FE_nodesAtIP ( FE_maxmaxNnodesAtIP , FE_maxNips , FE_Nelemtypes ) ) ; FE_nodesAtIP = 0_pInt
2009-04-06 18:55:19 +05:30
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
! fill FE_nodesAtIP with data
2010-05-10 20:24:59 +05:30
FE_nodesAtIP ( : , : FE_Nips ( 1 ) , 1 ) = & ! element 7
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-05-10 20:24:59 +05:30
1 , &
2 , &
4 , &
3 , &
5 , &
6 , &
8 , &
7 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_maxNnodesAtIP ( 1 ) , FE_Nips ( 1 ) / ) )
2010-05-10 20:24:59 +05:30
FE_nodesAtIP ( : , : FE_Nips ( 2 ) , 2 ) = & ! element 134
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-05-10 20:24:59 +05:30
1 , 2 , 3 , 4 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_maxNnodesAtIP ( 2 ) , FE_Nips ( 2 ) / ) )
2010-05-10 20:24:59 +05:30
FE_nodesAtIP ( : , : FE_Nips ( 3 ) , 3 ) = & ! element 11
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-05-10 20:24:59 +05:30
1 , &
2 , &
4 , &
3 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_maxNnodesAtIP ( 3 ) , FE_Nips ( 3 ) / ) )
2010-05-10 20:24:59 +05:30
FE_nodesAtIP ( : , : FE_Nips ( 4 ) , 4 ) = & ! element 27
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-05-10 20:24:59 +05:30
1 , 0 , &
1 , 2 , &
2 , 0 , &
1 , 4 , &
0 , 0 , &
2 , 3 , &
4 , 0 , &
3 , 4 , &
3 , 0 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_maxNnodesAtIP ( 4 ) , FE_Nips ( 4 ) / ) )
2010-05-10 20:24:59 +05:30
FE_nodesAtIP ( : , : FE_Nips ( 5 ) , 5 ) = & ! element 157
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-05-10 20:24:59 +05:30
1 , &
2 , &
3 , &
4 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_maxNnodesAtIP ( 5 ) , FE_Nips ( 5 ) / ) )
2010-05-10 20:24:59 +05:30
FE_nodesAtIP ( : , : FE_Nips ( 6 ) , 6 ) = & ! element 136
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-05-10 20:24:59 +05:30
1 , &
2 , &
3 , &
4 , &
5 , &
6 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_maxNnodesAtIP ( 6 ) , FE_Nips ( 6 ) / ) )
2010-08-17 04:23:24 +05:30
FE_nodesAtIP ( : , : FE_Nips ( 7 ) , 7 ) = & ! element 21
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
1 , 0 , 0 , 0 , &
1 , 2 , 0 , 0 , &
2 , 0 , 0 , 0 , &
1 , 4 , 0 , 0 , &
1 , 3 , 2 , 4 , &
2 , 3 , 0 , 0 , &
4 , 0 , 0 , 0 , &
3 , 4 , 0 , 0 , &
3 , 0 , 0 , 0 , &
1 , 5 , 0 , 0 , &
1 , 6 , 2 , 5 , &
2 , 6 , 0 , 0 , &
1 , 8 , 4 , 5 , &
0 , 0 , 0 , 0 , &
2 , 7 , 3 , 6 , &
4 , 8 , 0 , 0 , &
3 , 8 , 4 , 7 , &
3 , 7 , 0 , 0 , &
5 , 0 , 0 , 0 , &
5 , 6 , 0 , 0 , &
6 , 0 , 0 , 0 , &
5 , 8 , 0 , 0 , &
5 , 7 , 6 , 8 , &
6 , 7 , 0 , 0 , &
8 , 0 , 0 , 0 , &
7 , 8 , 0 , 0 , &
7 , 0 , 0 , 0 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_maxNnodesAtIP ( 7 ) , FE_Nips ( 7 ) / ) )
2011-02-16 21:53:08 +05:30
! 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)/))
2010-08-17 04:23:24 +05:30
FE_nodesAtIP ( : , : FE_Nips ( 9 ) , 9 ) = & ! element 57 (c3d20r == c3d8 --> copy of 7)
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-08-17 04:23:24 +05:30
1 , &
2 , &
4 , &
3 , &
5 , &
6 , &
8 , &
7 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_maxNnodesAtIP ( 9 ) , FE_Nips ( 9 ) / ) )
2011-04-06 14:05:37 +05:30
FE_nodesAtIP ( : , : FE_Nips ( 10 ) , 10 ) = & ! element 155, 125, 128
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2011-04-06 14:05:37 +05:30
1 , &
2 , &
3 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_maxNnodesAtIP ( 10 ) , FE_Nips ( 10 ) / ) )
2009-04-06 18:55:19 +05:30
2011-04-06 14:05:37 +05:30
! *** 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.
2009-04-06 18:55:19 +05:30
FE_ipNeighbor ( : FE_NipNeighbors ( 1 ) , : FE_Nips ( 1 ) , 1 ) = & ! element 7
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipNeighbors ( 1 ) , FE_Nips ( 1 ) / ) )
2009-04-06 18:55:19 +05:30
FE_ipNeighbor ( : FE_NipNeighbors ( 2 ) , : FE_Nips ( 2 ) , 2 ) = & ! element 134
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
- 1 , - 2 , - 3 , - 4 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipNeighbors ( 2 ) , FE_Nips ( 2 ) / ) )
2009-04-06 18:55:19 +05:30
FE_ipNeighbor ( : FE_NipNeighbors ( 3 ) , : FE_Nips ( 3 ) , 3 ) = & ! element 11
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
2 , - 4 , 3 , - 1 , &
- 2 , 1 , 4 , - 1 , &
4 , - 4 , - 3 , 1 , &
- 2 , 3 , - 3 , 2 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipNeighbors ( 3 ) , FE_Nips ( 3 ) / ) )
2009-04-06 18:55:19 +05:30
FE_ipNeighbor ( : FE_NipNeighbors ( 4 ) , : FE_Nips ( 4 ) , 4 ) = & ! element 27
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
2 , - 4 , 4 , - 1 , &
3 , 1 , 5 , - 1 , &
- 2 , 2 , 6 , - 1 , &
5 , - 4 , 7 , 1 , &
6 , 4 , 8 , 2 , &
- 2 , 5 , 9 , 3 , &
8 , - 4 , - 3 , 4 , &
9 , 7 , - 3 , 5 , &
- 2 , 8 , - 3 , 6 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipNeighbors ( 4 ) , FE_Nips ( 4 ) / ) )
2009-04-06 18:55:19 +05:30
FE_ipNeighbor ( : FE_NipNeighbors ( 5 ) , : FE_Nips ( 5 ) , 5 ) = & ! element 157
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
2 , - 4 , 3 , - 2 , 4 , - 1 , &
3 , - 2 , 1 , - 3 , 4 , - 1 , &
1 , - 3 , 2 , - 4 , 4 , - 1 , &
1 , - 3 , 2 , - 4 , 3 , - 2 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipNeighbors ( 5 ) , FE_Nips ( 5 ) / ) )
2009-04-06 18:55:19 +05:30
FE_ipNeighbor ( : FE_NipNeighbors ( 6 ) , : FE_Nips ( 6 ) , 6 ) = & ! element 136
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
2 , - 4 , 3 , - 2 , 4 , - 1 , &
- 3 , 1 , 3 , - 2 , 5 , - 1 , &
2 , - 4 , - 3 , 1 , 6 , - 1 , &
5 , - 4 , 6 , - 2 , - 5 , 1 , &
- 3 , 4 , 6 , - 2 , - 5 , 2 , &
5 , - 4 , - 3 , 4 , - 5 , 3 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipNeighbors ( 6 ) , FE_Nips ( 6 ) / ) )
2009-04-06 18:55:19 +05:30
FE_ipNeighbor ( : FE_NipNeighbors ( 7 ) , : FE_Nips ( 7 ) , 7 ) = & ! element 21
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-04-06 17:15:23 +05:30
2 , - 5 , 4 , - 2 , 10 , - 1 , &
2009-04-06 18:55:19 +05:30
3 , 1 , 5 , - 2 , 11 , - 1 , &
- 3 , 2 , 6 , - 2 , 12 , - 1 , &
5 , - 5 , 7 , 1 , 13 , - 1 , &
6 , 4 , 8 , 2 , 14 , - 1 , &
- 3 , 5 , 9 , 3 , 15 , - 1 , &
8 , - 5 , - 4 , 4 , 16 , - 1 , &
9 , 7 , - 4 , 5 , 17 , - 1 , &
- 3 , 8 , - 4 , 6 , 18 , - 1 , &
2010-04-06 17:15:23 +05:30
11 , - 5 , 13 , - 2 , 19 , 1 , &
12 , 10 , 14 , - 2 , 20 , 2 , &
- 3 , 11 , 15 , - 2 , 21 , 3 , &
14 , - 5 , 16 , 10 , 22 , 4 , &
15 , 13 , 17 , 11 , 23 , 5 , &
- 3 , 14 , 18 , 12 , 24 , 6 , &
17 , - 5 , - 4 , 13 , 25 , 7 , &
18 , 16 , - 4 , 14 , 26 , 8 , &
- 3 , 17 , - 4 , 15 , 27 , 9 , &
2009-04-06 18:55:19 +05:30
20 , - 5 , 22 , - 2 , - 6 , 10 , &
21 , 19 , 23 , - 2 , - 6 , 11 , &
- 3 , 20 , 24 , - 2 , - 6 , 12 , &
23 , - 5 , 25 , 19 , - 6 , 13 , &
24 , 22 , 26 , 20 , - 6 , 14 , &
- 3 , 23 , 27 , 21 , - 6 , 15 , &
26 , - 5 , - 4 , 22 , - 6 , 16 , &
27 , 25 , - 4 , 23 , - 6 , 17 , &
- 3 , 26 , - 4 , 24 , - 6 , 18 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipNeighbors ( 7 ) , FE_Nips ( 7 ) / ) )
2010-05-10 20:24:59 +05:30
FE_ipNeighbor ( : FE_NipNeighbors ( 8 ) , : FE_Nips ( 8 ) , 8 ) = & ! element 117
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-05-10 20:24:59 +05:30
- 3 , - 5 , - 4 , - 2 , - 6 , - 1 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipNeighbors ( 8 ) , FE_Nips ( 8 ) / ) )
2010-08-17 04:23:24 +05:30
FE_ipNeighbor ( : FE_NipNeighbors ( 9 ) , : FE_Nips ( 9 ) , 9 ) = & ! element 57 (c3d20r == c3d8 --> copy of 7)
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-08-17 04:23:24 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipNeighbors ( 9 ) , FE_Nips ( 9 ) / ) )
2011-04-06 14:05:37 +05:30
FE_ipNeighbor ( : FE_NipNeighbors ( 10 ) , : FE_Nips ( 10 ) , 10 ) = & ! element 155, 125, 128
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2011-04-06 14:05:37 +05:30
2 , - 3 , 3 , - 1 , &
- 2 , 1 , 3 , - 1 , &
2 , - 3 , - 2 , 1 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipNeighbors ( 10 ) , FE_Nips ( 10 ) / ) )
2011-04-06 14:05:37 +05:30
! *** FE_subNodeParent ***
2011-05-24 21:27:59 +05:30
! lists the group of nodes for which the center of gravity
2011-04-06 14:05:37 +05:30
! corresponds to the location of a each subnode.
! 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
2009-04-06 18:55:19 +05:30
FE_subNodeParent ( : FE_Nips ( 1 ) , : FE_NsubNodes ( 1 ) , 1 ) = & ! element 7
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_Nips ( 1 ) , FE_NsubNodes ( 1 ) / ) )
2010-05-10 20:24:59 +05:30
!FE_subNodeParent(:FE_Nips(2),:FE_NsubNodes(2),2) ! element 134 has no subnodes
2009-04-06 18:55:19 +05:30
FE_subNodeParent ( : FE_Nips ( 3 ) , : FE_NsubNodes ( 3 ) , 3 ) = & ! element 11
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
1 , 2 , 0 , 0 , &
2 , 3 , 0 , 0 , &
3 , 4 , 0 , 0 , &
4 , 1 , 0 , 0 , &
1 , 2 , 3 , 4 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_Nips ( 3 ) , FE_NsubNodes ( 3 ) / ) )
2009-04-06 18:55:19 +05:30
FE_subNodeParent ( : FE_Nips ( 4 ) , : FE_NsubNodes ( 4 ) , 4 ) = & ! element 27
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_Nips ( 4 ) , FE_NsubNodes ( 4 ) / ) )
2009-04-06 18:55:19 +05:30
!FE_subNodeParent(:FE_Nips(5),:FE_NsubNodes(5),5) = & ! element 157
! reshape((/&
! *still to be defined*
2012-02-16 00:28:38 +05:30
! ],pInt),(/FE_Nips(5),FE_NsubNodes(5)/))
2009-04-06 18:55:19 +05:30
FE_subNodeParent ( : FE_Nips ( 6 ) , : FE_NsubNodes ( 6 ) , 6 ) = & ! element 136
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_Nips ( 6 ) , FE_NsubNodes ( 6 ) / ) )
2009-04-06 18:55:19 +05:30
FE_subNodeParent ( : FE_Nips ( 7 ) , : FE_NsubNodes ( 7 ) , 7 ) = & ! element 21
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
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 , &
2 , 2 , 3 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
2 , 3 , 3 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
3 , 3 , 4 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
3 , 4 , 4 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
4 , 4 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
4 , 1 , 1 , 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 , 1 , 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
2 , 2 , 6 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
3 , 3 , 7 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
4 , 4 , 8 , 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 , 5 , 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
2 , 6 , 6 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
3 , 7 , 7 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
4 , 8 , 8 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
5 , 5 , 6 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
5 , 6 , 6 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
6 , 6 , 7 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
6 , 7 , 7 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
7 , 7 , 8 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
7 , 8 , 8 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
8 , 8 , 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
8 , 5 , 5 , 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 , 1 , 1 , 1 , 2 , 2 , 4 , 4 , 3 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
2 , 2 , 2 , 2 , 1 , 1 , 3 , 3 , 4 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
3 , 3 , 3 , 3 , 2 , 2 , 4 , 4 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
4 , 4 , 4 , 4 , 1 , 1 , 3 , 3 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
1 , 1 , 1 , 1 , 2 , 2 , 5 , 5 , 6 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
2 , 2 , 2 , 2 , 1 , 1 , 6 , 6 , 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
2 , 2 , 2 , 2 , 3 , 3 , 6 , 6 , 7 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
3 , 3 , 3 , 3 , 2 , 2 , 7 , 7 , 6 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
3 , 3 , 3 , 3 , 4 , 4 , 7 , 7 , 8 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
4 , 4 , 4 , 4 , 3 , 3 , 8 , 8 , 7 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
4 , 4 , 4 , 4 , 1 , 1 , 8 , 8 , 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
1 , 1 , 1 , 1 , 4 , 4 , 5 , 5 , 8 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
5 , 5 , 5 , 5 , 1 , 1 , 6 , 6 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
6 , 6 , 6 , 6 , 2 , 2 , 5 , 5 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
6 , 6 , 6 , 6 , 2 , 2 , 7 , 7 , 3 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
7 , 7 , 7 , 7 , 3 , 3 , 6 , 6 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
7 , 7 , 7 , 7 , 3 , 3 , 8 , 8 , 4 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
8 , 8 , 8 , 8 , 4 , 4 , 7 , 7 , 3 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
8 , 8 , 8 , 8 , 4 , 4 , 5 , 5 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
5 , 5 , 5 , 5 , 1 , 1 , 8 , 8 , 4 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
5 , 5 , 5 , 5 , 6 , 6 , 8 , 8 , 7 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
6 , 6 , 6 , 6 , 5 , 5 , 7 , 7 , 8 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
7 , 7 , 7 , 7 , 6 , 6 , 8 , 8 , 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
8 , 8 , 8 , 8 , 5 , 5 , 7 , 7 , 6 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 2 , 2 , 2 , 2 , 4 , 4 , 4 , 4 , 5 , 5 , 5 , 5 , 3 , 3 , 6 , 6 , 8 , 8 , 7 , &
2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 1 , 1 , 1 , 1 , 3 , 3 , 3 , 3 , 6 , 6 , 6 , 6 , 4 , 4 , 5 , 5 , 7 , 7 , 8 , &
3 , 3 , 3 , 3 , 3 , 3 , 3 , 3 , 2 , 2 , 2 , 2 , 4 , 4 , 4 , 4 , 7 , 7 , 7 , 7 , 1 , 1 , 6 , 6 , 8 , 8 , 5 , &
4 , 4 , 4 , 4 , 4 , 4 , 4 , 4 , 1 , 1 , 1 , 1 , 3 , 3 , 3 , 3 , 8 , 8 , 8 , 8 , 2 , 2 , 5 , 5 , 7 , 7 , 6 , &
5 , 5 , 5 , 5 , 5 , 5 , 5 , 5 , 1 , 1 , 1 , 1 , 6 , 6 , 6 , 6 , 8 , 8 , 8 , 8 , 2 , 2 , 4 , 4 , 7 , 7 , 3 , &
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_Nips ( 7 ) , FE_NsubNodes ( 7 ) / ) )
2010-05-10 20:24:59 +05:30
!FE_subNodeParent(:FE_Nips(8),:FE_NsubNodes(8),8) ! element 117 has no subnodes
2010-08-17 04:23:24 +05:30
FE_subNodeParent ( : FE_Nips ( 9 ) , : FE_NsubNodes ( 9 ) , 9 ) = & ! element 57 (c3d20r == c3d8 --> copy of 7)
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-08-17 04:23:24 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_Nips ( 9 ) , FE_NsubNodes ( 9 ) / ) )
2011-04-06 14:05:37 +05:30
FE_subNodeParent ( : FE_Nips ( 10 ) , : FE_NsubNodes ( 10 ) , 10 ) = & ! element 155, 125, 128
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2011-04-06 14:05:37 +05:30
1 , 2 , 0 , &
2 , 3 , 0 , &
3 , 1 , 0 , &
1 , 2 , 3 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_Nips ( 10 ) , FE_NsubNodes ( 10 ) / ) )
2011-04-06 14:05:37 +05:30
! *** FE_subNodeOnIPFace ***
! indicates which subnodes make up the interfaces enclosing the IP volume.
! The sorting convention is such that the outward pointing normal
! follows from a right-handed traversal of the face node list.
! For two-dimensional elements, which only have lines as "interface"
! one nevertheless has to specify each interface by a closed path,
! e.g., 1,2, 2,1, assuming the line connects nodes 1 and 2.
! 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.
2009-04-06 18:55:19 +05:30
FE_subNodeOnIPFace ( : FE_NipFaceNodes , : FE_NipNeighbors ( 1 ) , : FE_Nips ( 1 ) , 1 ) = & ! element 7
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipFaceNodes , FE_NipNeighbors ( 1 ) , FE_Nips ( 1 ) / ) )
2009-04-06 18:55:19 +05:30
FE_subNodeOnIPFace ( : FE_NipFaceNodes , : FE_NipNeighbors ( 2 ) , : FE_Nips ( 2 ) , 2 ) = & ! element 134
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
1 , 1 , 3 , 2 , & ! 1
1 , 1 , 2 , 4 , &
2 , 2 , 3 , 4 , &
1 , 1 , 4 , 3 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipFaceNodes , FE_NipNeighbors ( 2 ) , FE_Nips ( 2 ) / ) )
2009-04-06 18:55:19 +05:30
FE_subNodeOnIPFace ( : FE_NipFaceNodes , : FE_NipNeighbors ( 3 ) , : FE_Nips ( 3 ) , 3 ) = & ! element 11
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2011-04-06 14:05:37 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipFaceNodes , FE_NipNeighbors ( 3 ) , FE_Nips ( 3 ) / ) )
2009-04-06 18:55:19 +05:30
FE_subNodeOnIPFace ( : FE_NipFaceNodes , : FE_NipNeighbors ( 4 ) , : FE_Nips ( 4 ) , 4 ) = & ! element 27
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2011-04-06 14:05:37 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipFaceNodes , FE_NipNeighbors ( 4 ) , FE_Nips ( 4 ) / ) )
2009-04-06 18:55:19 +05:30
!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 ( : FE_NipFaceNodes , : FE_NipNeighbors ( 6 ) , : FE_Nips ( 6 ) , 6 ) = & ! element 136
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipFaceNodes , FE_NipNeighbors ( 6 ) , FE_Nips ( 6 ) / ) )
2009-04-06 18:55:19 +05:30
FE_subNodeOnIPFace ( : FE_NipFaceNodes , : FE_NipNeighbors ( 7 ) , : FE_Nips ( 7 ) , 7 ) = & ! element 21
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2009-04-06 18:55:19 +05:30
9 , 33 , 57 , 37 , & ! 1
1 , 17 , 44 , 16 , &
33 , 16 , 44 , 57 , &
1 , 9 , 37 , 17 , &
17 , 37 , 57 , 44 , &
1 , 16 , 33 , 9 , &
10 , 34 , 58 , 38 , & ! 2
9 , 37 , 57 , 33 , &
34 , 33 , 57 , 58 , &
9 , 10 , 38 , 37 , &
37 , 38 , 58 , 57 , &
9 , 33 , 34 , 10 , &
2 , 11 , 39 , 18 , & ! 3
10 , 38 , 58 , 34 , &
11 , 34 , 58 , 39 , &
10 , 2 , 18 , 38 , &
38 , 18 , 39 , 58 , &
10 , 34 , 11 , 2 , &
33 , 36 , 60 , 57 , & ! 4
16 , 44 , 43 , 15 , &
36 , 15 , 43 , 60 , &
16 , 33 , 57 , 44 , &
44 , 57 , 60 , 43 , &
16 , 15 , 36 , 33 , &
34 , 35 , 59 , 58 , & ! 5
33 , 57 , 60 , 36 , &
35 , 36 , 60 , 59 , &
33 , 34 , 58 , 57 , &
57 , 58 , 59 , 60 , &
33 , 36 , 35 , 34 , &
11 , 12 , 40 , 39 , & ! 6
34 , 58 , 59 , 35 , &
12 , 35 , 59 , 40 , &
34 , 11 , 39 , 58 , &
58 , 39 , 40 , 59 , &
34 , 35 , 12 , 11 , &
36 , 14 , 42 , 60 , & ! 7
15 , 43 , 20 , 4 , &
14 , 4 , 20 , 42 , &
15 , 36 , 60 , 43 , &
43 , 60 , 42 , 20 , &
15 , 4 , 14 , 36 , &
35 , 13 , 41 , 59 , & ! 8
36 , 60 , 42 , 14 , &
13 , 14 , 42 , 41 , &
36 , 35 , 59 , 60 , &
60 , 59 , 41 , 42 , &
36 , 14 , 13 , 35 , &
12 , 3 , 19 , 40 , & ! 9
35 , 59 , 41 , 13 , &
3 , 13 , 41 , 19 , &
35 , 12 , 40 , 59 , &
59 , 40 , 19 , 41 , &
35 , 13 , 3 , 12 , &
37 , 57 , 61 , 45 , & ! 10
17 , 21 , 52 , 44 , &
57 , 44 , 52 , 61 , &
17 , 37 , 45 , 21 , &
21 , 45 , 61 , 52 , &
17 , 44 , 57 , 37 , &
38 , 58 , 62 , 46 , & ! 11
37 , 45 , 61 , 57 , &
58 , 57 , 61 , 62 , &
37 , 38 , 46 , 45 , &
45 , 46 , 62 , 61 , &
37 , 57 , 58 , 38 , &
18 , 39 , 47 , 22 , & ! 12
38 , 46 , 62 , 58 , &
39 , 58 , 62 , 47 , &
38 , 18 , 22 , 46 , &
46 , 22 , 47 , 62 , &
38 , 58 , 39 , 18 , &
57 , 60 , 64 , 61 , & ! 13
44 , 52 , 51 , 43 , &
60 , 43 , 51 , 64 , &
44 , 57 , 61 , 52 , &
52 , 61 , 64 , 51 , &
44 , 43 , 60 , 57 , &
58 , 59 , 63 , 62 , & ! 14
57 , 61 , 64 , 60 , &
59 , 60 , 64 , 63 , &
57 , 58 , 62 , 61 , &
61 , 62 , 63 , 64 , &
57 , 60 , 59 , 58 , &
39 , 40 , 48 , 47 , & ! 15
58 , 62 , 63 , 59 , &
40 , 59 , 63 , 48 , &
58 , 39 , 47 , 62 , &
62 , 47 , 48 , 63 , &
58 , 59 , 40 , 39 , &
60 , 42 , 50 , 64 , & ! 16
43 , 51 , 24 , 20 , &
42 , 20 , 24 , 50 , &
43 , 60 , 64 , 51 , &
51 , 64 , 50 , 24 , &
43 , 20 , 42 , 60 , &
59 , 41 , 49 , 63 , & ! 17
60 , 64 , 50 , 42 , &
41 , 42 , 50 , 49 , &
60 , 59 , 63 , 64 , &
64 , 63 , 49 , 50 , &
60 , 42 , 41 , 59 , &
40 , 19 , 23 , 48 , & ! 18
59 , 63 , 49 , 41 , &
19 , 41 , 49 , 23 , &
59 , 40 , 48 , 63 , &
63 , 48 , 23 , 49 , &
59 , 41 , 19 , 40 , &
45 , 61 , 53 , 25 , & ! 19
21 , 5 , 32 , 52 , &
61 , 52 , 32 , 53 , &
21 , 45 , 25 , 5 , &
5 , 25 , 53 , 32 , &
21 , 52 , 61 , 45 , &
46 , 62 , 54 , 26 , & ! 20
45 , 25 , 53 , 61 , &
62 , 61 , 53 , 54 , &
45 , 46 , 26 , 25 , &
25 , 26 , 54 , 53 , &
45 , 61 , 62 , 46 , &
22 , 47 , 27 , 6 , & ! 21
46 , 26 , 54 , 62 , &
47 , 62 , 54 , 27 , &
46 , 22 , 6 , 26 , &
26 , 6 , 27 , 54 , &
46 , 62 , 47 , 22 , &
61 , 64 , 56 , 53 , & ! 22
52 , 32 , 31 , 51 , &
64 , 51 , 31 , 56 , &
52 , 61 , 53 , 32 , &
32 , 53 , 56 , 31 , &
52 , 51 , 64 , 61 , &
62 , 63 , 55 , 54 , & ! 23
61 , 53 , 56 , 64 , &
63 , 64 , 56 , 55 , &
61 , 62 , 54 , 53 , &
53 , 54 , 55 , 56 , &
61 , 64 , 63 , 62 , &
47 , 48 , 28 , 27 , & ! 24
62 , 54 , 55 , 63 , &
48 , 63 , 55 , 28 , &
62 , 47 , 27 , 54 , &
54 , 27 , 28 , 55 , &
62 , 63 , 48 , 47 , &
64 , 50 , 30 , 56 , & ! 25
51 , 31 , 8 , 24 , &
50 , 24 , 8 , 30 , &
51 , 64 , 56 , 31 , &
31 , 56 , 30 , 8 , &
51 , 24 , 50 , 64 , &
63 , 49 , 29 , 55 , & ! 26
64 , 56 , 30 , 50 , &
49 , 50 , 30 , 29 , &
64 , 63 , 55 , 56 , &
56 , 55 , 29 , 30 , &
64 , 50 , 49 , 63 , &
48 , 23 , 7 , 28 , & ! 27
63 , 55 , 29 , 49 , &
23 , 49 , 29 , 7 , &
63 , 48 , 28 , 55 , &
55 , 28 , 7 , 29 , &
63 , 49 , 23 , 48 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipFaceNodes , FE_NipNeighbors ( 7 ) , FE_Nips ( 7 ) / ) )
2010-05-10 20:24:59 +05:30
FE_subNodeOnIPFace ( : FE_NipFaceNodes , : FE_NipNeighbors ( 8 ) , : FE_Nips ( 8 ) , 8 ) = & ! element 117
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-05-10 20:24:59 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipFaceNodes , FE_NipNeighbors ( 8 ) , FE_Nips ( 8 ) / ) )
2010-08-17 04:23:24 +05:30
FE_subNodeOnIPFace ( : FE_NipFaceNodes , : FE_NipNeighbors ( 9 ) , : FE_Nips ( 9 ) , 9 ) = & ! element 57 (c3d20r == c3d8 --> copy of 7)
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2010-08-17 04:23:24 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipFaceNodes , FE_NipNeighbors ( 9 ) , FE_Nips ( 9 ) / ) )
2011-04-06 14:05:37 +05:30
FE_subNodeOnIPFace ( : FE_NipFaceNodes , : FE_NipNeighbors ( 10 ) , : FE_Nips ( 10 ) , 10 ) = & ! element 155, 125, 128
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2011-04-06 14:05:37 +05:30
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 &
2012-02-16 00:28:38 +05:30
] , pInt ) , ( / FE_NipFaceNodes , FE_NipNeighbors ( 10 ) , FE_Nips ( 10 ) / ) )
2009-04-06 18:55:19 +05:30
2009-06-15 18:41:21 +05:30
endsubroutine
2009-04-06 18:55:19 +05:30
2009-10-12 21:31:49 +05:30
2007-04-10 16:52:53 +05:30
!********************************************************************
2009-10-12 21:31:49 +05:30
! figure out table styles (Marc only)
2007-04-10 16:52:53 +05:30
!
2009-10-12 21:31:49 +05:30
! initialcondTableStyle, hypoelasticTableStyle
2007-04-10 16:52:53 +05:30
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_marc_get_tableStyles ( myUnit )
2009-01-20 00:12:31 +05:30
2007-03-26 14:20:57 +05:30
use prec , only : pInt
2007-04-10 16:52:53 +05:30
use IO
2007-03-26 14:20:57 +05:30
implicit none
2009-01-20 00:12:31 +05:30
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 6_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2009-01-20 00:12:31 +05:30
2009-10-12 21:31:49 +05:30
initialcondTableStyle = 0_pInt
hypoelasticTableStyle = 0_pInt
2007-04-10 16:52:53 +05:30
610 FORMAT ( A300 )
2009-01-20 00:12:31 +05:30
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2007-04-10 16:52:53 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2009-01-20 00:12:31 +05:30
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == 'table' . and . myPos ( 1_pInt ) . GT . 5 ) then
initialcondTableStyle = IO_intValue ( line , myPos , 4_pInt )
hypoelasticTableStyle = IO_intValue ( line , myPos , 5_pInt )
2009-10-12 21:31:49 +05:30
exit
endif
2009-06-15 18:41:21 +05:30
enddo
2009-01-20 00:12:31 +05:30
2011-09-13 21:24:06 +05:30
620 endsubroutine
2009-01-20 00:12:31 +05:30
2011-02-16 21:53:08 +05:30
!********************************************************************
2012-02-03 18:07:52 +05:30
! get any additional damask options from input file
2011-02-16 21:53:08 +05:30
!
! mesh_periodicSurface
!********************************************************************
2012-02-03 18:07:52 +05:30
subroutine mesh_get_damaskOptions ( myUnit )
2011-02-16 21:53:08 +05:30
2012-02-03 18:07:52 +05:30
use DAMASK_interface , only : FEsolver
2011-02-16 21:53:08 +05:30
use prec , only : pInt
use IO
implicit none
2011-10-18 14:50:29 +05:30
integer ( pInt ) , intent ( in ) :: myUnit
2011-02-16 21:53:08 +05:30
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 5_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: myPos
2012-02-03 18:07:52 +05:30
integer ( pInt ) chunk , Nchunks
2012-02-10 17:26:05 +05:30
character ( len = 300 ) line , keyword , damaskOption , v
2011-02-16 21:53:08 +05:30
mesh_periodicSurface = ( / . false . , . false . , . false . / )
610 FORMAT ( A300 )
2012-02-03 18:07:52 +05:30
select case ( FEsolver )
case ( 'Spectral' ) ! no special keyword needed, the damask option directly goes into the header
case ( 'Marc' )
keyword = '$damask'
case ( 'Abaqus' )
keyword = '**damask'
end select
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2011-02-16 21:53:08 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-03 18:07:52 +05:30
Nchunks = myPos ( 1 )
select case ( FEsolver )
case ( 'Marc' , 'Abaqus' )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == keyword . and . Nchunks > 1_pInt ) then ! found keyword for damask option and there is at least one more chunk to read
damaskOption = IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) )
2012-02-03 18:07:52 +05:30
select case ( damaskOption )
case ( 'periodic' ) ! damask Option that allows to specify periodic fluxes
2012-02-10 17:26:05 +05:30
do chunk = 3_pInt , Nchunks ! loop through chunks (skipping the keyword)
v = IO_lc ( IO_stringValue ( line , myPos , chunk ) ) ! chunk matches keyvalues x,y, or z?
mesh_periodicSurface ( 1 ) = mesh_periodicSurface ( 1 ) . or . v == 'x'
mesh_periodicSurface ( 2 ) = mesh_periodicSurface ( 2 ) . or . v == 'y'
mesh_periodicSurface ( 3 ) = mesh_periodicSurface ( 3 ) . or . v == 'z'
2012-02-03 18:07:52 +05:30
enddo
endselect
endif
case ( 'Spectral' )
2012-02-10 17:26:05 +05:30
damaskOption = IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) )
2012-02-03 18:07:52 +05:30
select case ( damaskOption )
case ( 'periodic' ) ! damask Option that allows to specify periodic fluxes
2012-02-10 17:26:05 +05:30
do chunk = 2_pInt , Nchunks ! loop through chunks (skipping the keyword)
v = IO_lc ( IO_stringValue ( line , myPos , chunk ) ) ! chunk matches keyvalues x,y, or z?
mesh_periodicSurface ( 1 ) = mesh_periodicSurface ( 1 ) . or . v == 'x'
mesh_periodicSurface ( 2 ) = mesh_periodicSurface ( 2 ) . or . v == 'y'
mesh_periodicSurface ( 3 ) = mesh_periodicSurface ( 3 ) . or . v == 'z'
2012-02-03 18:07:52 +05:30
enddo
endselect
endselect
2011-02-16 21:53:08 +05:30
enddo
2011-09-13 21:24:06 +05:30
620 endsubroutine
2011-02-16 21:53:08 +05:30
2007-03-26 14:20:57 +05:30
2010-05-06 22:10:47 +05:30
!********************************************************************
! count overall number of nodes and elements in mesh
!
! mesh_Nelems, mesh_Nnodes
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_spectral_count_nodesAndElements ( myUnit )
2010-05-06 22:10:47 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 7_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: myPos
integer ( pInt ) a , b , c , i , j , headerLength
2010-05-06 22:10:47 +05:30
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit
character ( len = 1024 ) line , keyword
2010-05-06 22:10:47 +05:30
mesh_Nnodes = 0_pInt
mesh_Nelems = 0_pInt
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
read ( myUnit , '(a1024)' ) line
2012-02-10 17:26:05 +05:30
myPos = IO_stringPos ( line , 2_pInt )
keyword = IO_lc ( IO_StringValue ( line , myPos , 2_pInt ) )
2011-10-18 14:50:29 +05:30
if ( keyword ( 1 : 4 ) == 'head' ) then
2012-02-10 17:26:05 +05:30
headerLength = IO_intValue ( line , myPos , 1_pInt ) + 1_pInt
2011-10-18 14:50:29 +05:30
else
2012-02-13 23:11:27 +05:30
call IO_error ( error_ID = 842_pInt )
2011-10-18 14:50:29 +05:30
endif
rewind ( myUnit )
2012-02-16 00:28:38 +05:30
do i = 1_pInt , headerLength
2011-10-18 14:50:29 +05:30
read ( myUnit , '(a1024)' ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_StringValue ( line , myPos , 1_pInt ) ) == 'resolution' ) then
2012-02-16 00:28:38 +05:30
do j = 2_pInt , 6_pInt , 2_pInt
2011-10-18 14:50:29 +05:30
select case ( IO_lc ( IO_stringValue ( line , myPos , j ) ) )
2010-05-10 20:24:59 +05:30
case ( 'a' )
2012-02-10 17:26:05 +05:30
a = IO_intValue ( line , myPos , j + 1_pInt )
2010-05-10 20:24:59 +05:30
case ( 'b' )
2012-02-10 17:26:05 +05:30
b = IO_intValue ( line , myPos , j + 1_pInt )
2010-05-10 20:24:59 +05:30
case ( 'c' )
2012-02-10 17:26:05 +05:30
c = IO_intValue ( line , myPos , j + 1_pInt )
2010-05-06 22:10:47 +05:30
end select
enddo
2010-10-01 16:12:15 +05:30
mesh_Nelems = a * b * c
2012-02-16 00:28:38 +05:30
mesh_Nnodes = ( 1_pInt + a ) * ( 1_pInt + b ) * ( 1_pInt + c )
2010-05-06 22:10:47 +05:30
endif
enddo
2011-10-18 14:50:29 +05:30
endsubroutine
2010-05-06 22:10:47 +05:30
2009-10-12 21:31:49 +05:30
!********************************************************************
! count overall number of nodes and elements in mesh
2007-03-27 18:23:31 +05:30
!
2009-10-12 21:31:49 +05:30
! mesh_Nelems, mesh_Nnodes
2007-04-10 16:52:53 +05:30
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_marc_count_nodesAndElements ( myUnit )
2009-10-12 21:31:49 +05:30
2007-03-21 21:48:33 +05:30
use prec , only : pInt
2007-04-10 16:52:53 +05:30
use IO
2007-03-21 21:48:33 +05:30
implicit none
2009-10-12 21:31:49 +05:30
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 4_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
mesh_Nnodes = 0_pInt
mesh_Nelems = 0_pInt
2007-04-10 16:52:53 +05:30
610 FORMAT ( A300 )
2009-10-12 21:31:49 +05:30
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2009-10-12 21:31:49 +05:30
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_StringValue ( line , myPos , 1_pInt ) ) == 'sizing' ) then
mesh_Nelems = IO_IntValue ( line , myPos , 3_pInt )
mesh_Nnodes = IO_IntValue ( line , myPos , 4_pInt )
2007-09-28 20:26:26 +05:30
exit
2009-06-15 18:41:21 +05:30
endif
enddo
2009-10-12 21:31:49 +05:30
2011-09-13 21:24:06 +05:30
620 endsubroutine
2009-04-03 12:33:25 +05:30
2007-10-23 18:39:46 +05:30
!********************************************************************
2009-10-12 21:31:49 +05:30
! count overall number of nodes and elements in mesh
2007-10-23 18:39:46 +05:30
!
2009-10-12 21:31:49 +05:30
! mesh_Nelems, mesh_Nnodes
2007-10-23 18:39:46 +05:30
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_count_nodesAndElements ( myUnit )
2009-01-20 00:12:31 +05:30
2007-10-23 18:39:46 +05:30
use prec , only : pInt
use IO
implicit none
2009-01-20 00:12:31 +05:30
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 2_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2009-01-20 00:12:31 +05:30
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit
2009-10-12 21:31:49 +05:30
logical inPart
2009-01-20 00:12:31 +05:30
2012-02-16 00:28:38 +05:30
mesh_Nnodes = 0_pInt
mesh_Nelems = 0_pInt
2009-10-12 21:31:49 +05:30
610 FORMAT ( A300 )
2009-01-20 00:12:31 +05:30
2009-10-12 21:31:49 +05:30
inPart = . false .
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'part' ) inPart = . false .
2009-10-12 21:31:49 +05:30
2010-07-13 15:56:07 +05:30
if ( inPart . or . noPart ) then
2012-02-10 17:26:05 +05:30
select case ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) )
2009-10-12 21:31:49 +05:30
case ( '*node' )
if ( &
2012-02-10 17:26:05 +05:30
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'print' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'file' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'response' &
2009-10-12 21:31:49 +05:30
) &
2011-10-18 14:50:29 +05:30
mesh_Nnodes = mesh_Nnodes + IO_countDataLines ( myUnit )
2009-10-12 21:31:49 +05:30
case ( '*element' )
if ( &
2012-02-10 17:26:05 +05:30
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'matrix' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'response' &
2009-10-12 21:31:49 +05:30
) then
2011-10-18 14:50:29 +05:30
mesh_Nelems = mesh_Nelems + IO_countDataLines ( myUnit )
2009-10-12 21:31:49 +05:30
endif
endselect
2009-06-15 18:41:21 +05:30
endif
enddo
2010-07-13 15:56:07 +05:30
2012-02-16 00:28:38 +05:30
620 if ( mesh_Nnodes < 2_pInt ) call IO_error ( error_ID = 900_pInt )
if ( mesh_Nelems == 0_pInt ) call IO_error ( error_ID = 901_pInt )
2009-10-12 21:31:49 +05:30
2009-06-15 18:41:21 +05:30
endsubroutine
2009-01-20 00:12:31 +05:30
2007-04-25 13:03:24 +05:30
!********************************************************************
2009-10-12 21:31:49 +05:30
! count overall number of element sets in mesh
2007-04-25 13:03:24 +05:30
!
2009-10-12 21:31:49 +05:30
! mesh_NelemSets, mesh_maxNelemInSet
2007-04-25 13:03:24 +05:30
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_marc_count_elementSets ( myUnit )
2009-01-20 00:12:31 +05:30
2007-04-25 13:03:24 +05:30
use prec , only : pInt
use IO
implicit none
2009-01-20 00:12:31 +05:30
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 2_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: myPos
2009-01-20 00:12:31 +05:30
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2009-01-20 00:12:31 +05:30
2009-10-12 21:31:49 +05:30
mesh_NelemSets = 0_pInt
mesh_maxNelemInSet = 0_pInt
610 FORMAT ( A300 )
2009-01-20 00:12:31 +05:30
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2009-10-12 21:31:49 +05:30
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_StringValue ( line , myPos , 1_pInt ) ) == 'define' . and . &
IO_lc ( IO_StringValue ( line , myPos , 2_pInt ) ) == 'element' ) then
2009-10-12 21:31:49 +05:30
mesh_NelemSets = mesh_NelemSets + 1_pInt
mesh_maxNelemInSet = max ( mesh_maxNelemInSet , &
2011-10-18 14:50:29 +05:30
IO_countContinousIntValues ( myUnit ) )
2009-06-15 18:41:21 +05:30
endif
enddo
2009-01-20 00:12:31 +05:30
2011-09-13 21:24:06 +05:30
620 endsubroutine
2007-04-25 13:03:24 +05:30
2009-01-20 00:12:31 +05:30
2007-04-10 16:52:53 +05:30
!********************************************************************
2009-10-12 21:31:49 +05:30
! count overall number of element sets in mesh
2007-04-10 16:52:53 +05:30
!
2009-10-12 21:31:49 +05:30
! mesh_NelemSets, mesh_maxNelemInSet
2007-04-10 16:52:53 +05:30
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_count_elementSets ( myUnit )
2009-01-20 00:12:31 +05:30
2007-04-10 16:52:53 +05:30
use prec , only : pInt
use IO
implicit none
2009-01-20 00:12:31 +05:30
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 2_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2009-01-20 00:12:31 +05:30
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit
2009-10-12 21:31:49 +05:30
logical inPart
mesh_NelemSets = 0_pInt
mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons
2007-04-10 16:52:53 +05:30
610 FORMAT ( A300 )
2009-01-20 00:12:31 +05:30
2009-10-12 21:31:49 +05:30
inPart = . false .
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'part' ) inPart = . false .
2009-10-12 21:31:49 +05:30
2012-02-10 17:26:05 +05:30
if ( ( inPart . or . noPart ) . and . IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*elset' ) &
2009-10-12 21:31:49 +05:30
mesh_NelemSets = mesh_NelemSets + 1_pInt
2009-06-15 18:41:21 +05:30
enddo
2009-01-20 00:12:31 +05:30
2010-07-13 15:56:07 +05:30
620 continue
2012-02-10 17:26:05 +05:30
if ( mesh_NelemSets == 0 ) call IO_error ( error_ID = 902_pInt )
2010-07-13 15:56:07 +05:30
2009-06-15 18:41:21 +05:30
endsubroutine
2007-04-10 16:52:53 +05:30
2009-01-20 00:12:31 +05:30
2007-04-25 13:03:24 +05:30
!********************************************************************
2009-10-12 21:31:49 +05:30
! count overall number of solid sections sets in mesh (Abaqus only)
2007-04-10 16:52:53 +05:30
!
2009-10-12 21:31:49 +05:30
! mesh_Nmaterials
2007-04-25 13:03:24 +05:30
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_count_materials ( myUnit )
2009-01-20 00:12:31 +05:30
2007-04-25 13:03:24 +05:30
use prec , only : pInt
use IO
implicit none
2009-01-20 00:12:31 +05:30
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 2_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2009-01-20 00:12:31 +05:30
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit
2009-10-12 21:31:49 +05:30
logical inPart
mesh_Nmaterials = 0_pInt
2007-04-25 13:03:24 +05:30
610 FORMAT ( A300 )
2009-01-20 00:12:31 +05:30
2009-10-12 21:31:49 +05:30
inPart = . false .
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'part' ) inPart = . false .
2009-10-12 21:31:49 +05:30
2010-07-13 15:56:07 +05:30
if ( ( inPart . or . noPart ) . and . &
2012-02-10 17:26:05 +05:30
IO_lc ( IO_StringValue ( line , myPos , 1_pInt ) ) == '*solid' . and . &
IO_lc ( IO_StringValue ( line , myPos , 2_pInt ) ) == 'section' ) &
2009-10-12 21:31:49 +05:30
mesh_Nmaterials = mesh_Nmaterials + 1_pInt
2009-06-15 18:41:21 +05:30
enddo
2009-01-20 00:12:31 +05:30
2012-02-16 00:28:38 +05:30
620 if ( mesh_Nmaterials == 0_pInt ) call IO_error ( error_ID = 903_pInt )
2010-07-13 15:56:07 +05:30
2009-06-15 18:41:21 +05:30
endsubroutine
2007-04-25 13:03:24 +05:30
2009-01-20 00:12:31 +05:30
2010-05-06 22:10:47 +05:30
!********************************************************************
! count overall number of cpElements in mesh
!
! mesh_NcpElems
!********************************************************************
subroutine mesh_spectral_count_cpElements ( )
use prec , only : pInt
implicit none
mesh_NcpElems = mesh_Nelems
endsubroutine
2007-04-25 13:03:24 +05:30
!********************************************************************
2009-10-12 21:31:49 +05:30
! count overall number of cpElements in mesh
2007-04-10 16:52:53 +05:30
!
2009-10-12 21:31:49 +05:30
! mesh_NcpElems
2007-04-25 13:03:24 +05:30
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_marc_count_cpElements ( myUnit )
2009-01-20 00:12:31 +05:30
2007-04-25 13:03:24 +05:30
use prec , only : pInt
use IO
implicit none
2009-01-20 00:12:31 +05:30
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 1_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: myPos
2009-04-03 12:33:25 +05:30
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit , i
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2009-01-20 00:12:31 +05:30
2009-10-12 21:31:49 +05:30
mesh_NcpElems = 0_pInt
2009-01-20 00:12:31 +05:30
2007-04-25 13:03:24 +05:30
610 FORMAT ( A300 )
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2009-04-03 12:33:25 +05:30
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == 'hypoelastic' ) then
2012-02-16 00:28:38 +05:30
do i = 1_pInt , 3_pInt + hypoelasticTableStyle ! Skip 3 or 4 lines
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
2009-04-03 12:33:25 +05:30
enddo
2011-10-18 14:50:29 +05:30
mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues ( myUnit )
2009-10-12 21:31:49 +05:30
exit
2009-04-03 12:33:25 +05:30
endif
enddo
2009-01-20 00:12:31 +05:30
2011-09-13 21:24:06 +05:30
620 endsubroutine
2007-04-25 13:03:24 +05:30
2009-10-12 21:31:49 +05:30
2007-04-10 16:52:53 +05:30
!********************************************************************
2009-10-12 21:31:49 +05:30
! count overall number of cpElements in mesh
2007-04-10 16:52:53 +05:30
!
2009-10-12 21:31:49 +05:30
! mesh_NcpElems
2007-04-10 16:52:53 +05:30
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_count_cpElements ( myUnit )
2009-01-20 00:12:31 +05:30
2007-04-10 16:52:53 +05:30
use prec , only : pInt
use IO
implicit none
2009-01-20 00:12:31 +05:30
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 2_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2009-01-20 00:12:31 +05:30
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit , i , k
2009-10-12 21:31:49 +05:30
logical materialFound
character ( len = 64 ) materialName , elemSetName
2012-02-16 00:28:38 +05:30
mesh_NcpElems = 0_pInt
2009-10-12 21:31:49 +05:30
610 FORMAT ( A300 )
2009-01-20 00:12:31 +05:30
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
select case ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) )
2009-10-12 21:31:49 +05:30
case ( '*material' )
2012-02-10 17:26:05 +05:30
materialName = trim ( IO_extractValue ( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) , 'name' ) ) ! extract name=value
2011-03-23 21:50:12 +05:30
materialFound = materialName / = '' ! valid name?
2009-10-12 21:31:49 +05:30
case ( '*user' )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_StringValue ( line , myPos , 2_pInt ) ) == 'material' . and . materialFound ) then
2012-02-16 00:28:38 +05:30
do i = 1_pInt , mesh_Nmaterials ! look thru material names
2011-03-23 21:50:12 +05:30
if ( materialName == mesh_nameMaterial ( i ) ) then ! found one
elemSetName = mesh_mapMaterial ( i ) ! take corresponding elemSet
2012-02-16 00:28:38 +05:30
do k = 1_pInt , mesh_NelemSets ! look thru all elemSet definitions
2011-03-23 21:50:12 +05:30
if ( elemSetName == mesh_nameElemSet ( k ) ) & ! matched?
mesh_NcpElems = mesh_NcpElems + mesh_mapElemSet ( 1 , k ) ! add those elem count
enddo
endif
2009-10-12 21:31:49 +05:30
enddo
materialFound = . false .
endif
endselect
2009-06-15 18:41:21 +05:30
enddo
2011-03-23 21:50:12 +05:30
2012-02-16 00:28:38 +05:30
620 if ( mesh_NcpElems == 0_pInt ) call IO_error ( error_ID = 906_pInt )
2010-07-13 15:56:07 +05:30
2009-10-12 21:31:49 +05:30
endsubroutine
2009-01-20 00:12:31 +05:30
2009-10-12 21:31:49 +05:30
!********************************************************************
! map element sets
2007-04-10 16:52:53 +05:30
!
2009-10-12 21:31:49 +05:30
! allocate globals: mesh_nameElemSet, mesh_mapElemSet
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_marc_map_elementSets ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 4_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit , elemSet
2009-10-12 21:31:49 +05:30
allocate ( mesh_nameElemSet ( mesh_NelemSets ) ) ; mesh_nameElemSet = ''
2012-02-10 17:26:05 +05:30
allocate ( mesh_mapElemSet ( 1_pInt + mesh_maxNelemInSet , mesh_NelemSets ) ) ; mesh_mapElemSet = 0_pInt
2009-10-12 21:31:49 +05:30
610 FORMAT ( A300 )
elemSet = 0_pInt
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 640 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == 'define' ) . and . &
( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'element' ) ) then
2012-02-16 00:28:38 +05:30
elemSet = elemSet + 1_pInt
2012-02-10 17:26:05 +05:30
mesh_nameElemSet ( elemSet ) = trim ( IO_stringValue ( line , myPos , 4_pInt ) )
2011-10-18 14:50:29 +05:30
mesh_mapElemSet ( : , elemSet ) = IO_continousIntValues ( myUnit , mesh_maxNelemInSet , mesh_nameElemSet , mesh_mapElemSet , mesh_NelemSets )
2009-10-12 21:31:49 +05:30
endif
enddo
2011-09-13 21:24:06 +05:30
640 endsubroutine
2009-10-12 21:31:49 +05:30
!********************************************************************
! Build element set mapping
!
! allocate globals: mesh_nameElemSet, mesh_mapElemSet
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_map_elementSets ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-10 17:26:05 +05:30
integer ( pInt ) , parameter :: maxNchunks = 4_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit , elemSet , i
2009-10-12 21:31:49 +05:30
logical inPart
allocate ( mesh_nameElemSet ( mesh_NelemSets ) ) ; mesh_nameElemSet = ''
2012-02-10 17:26:05 +05:30
allocate ( mesh_mapElemSet ( 1_pInt + mesh_maxNelemInSet , mesh_NelemSets ) ) ; mesh_mapElemSet = 0_pInt
2009-10-12 21:31:49 +05:30
610 FORMAT ( A300 )
elemSet = 0_pInt
inPart = . false .
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 640 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'part' ) inPart = . false .
2009-10-12 21:31:49 +05:30
2012-02-10 17:26:05 +05:30
if ( ( inPart . or . noPart ) . and . IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*elset' ) then
2009-10-12 21:31:49 +05:30
elemSet = elemSet + 1_pInt
2012-02-10 17:26:05 +05:30
mesh_nameElemSet ( elemSet ) = trim ( IO_extractValue ( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) , 'elset' ) )
mesh_mapElemSet ( : , elemSet ) = IO_continousIntValues ( myUnit , mesh_Nelems , mesh_nameElemSet , &
mesh_mapElemSet , elemSet - 1_pInt )
2009-10-12 21:31:49 +05:30
endif
enddo
2012-02-10 17:26:05 +05:30
640 do i = 1_pInt , elemSet
2010-07-13 15:56:07 +05:30
! write(6,*)'elemSetName: ',mesh_nameElemSet(i)
! write(6,*)'elems in Elset',mesh_mapElemSet(:,i)
2012-02-10 17:26:05 +05:30
if ( mesh_mapElemSet ( 1 , i ) == 0_pInt ) call IO_error ( error_ID = 904_pInt , ext_msg = mesh_nameElemSet ( i ) )
2010-07-13 15:56:07 +05:30
enddo
2009-10-12 21:31:49 +05:30
endsubroutine
!********************************************************************
! map solid section (Abaqus only)
!
! allocate globals: mesh_nameMaterial, mesh_mapMaterial
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_map_materials ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 20_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2012-02-16 00:28:38 +05:30
integer ( pInt ) myUnit , i , c
2011-04-13 19:46:22 +05:30
logical inPart
2009-10-12 21:31:49 +05:30
character ( len = 64 ) elemSetName , materialName
allocate ( mesh_nameMaterial ( mesh_Nmaterials ) ) ; mesh_nameMaterial = ''
allocate ( mesh_mapMaterial ( mesh_Nmaterials ) ) ; mesh_mapMaterial = ''
610 FORMAT ( A300 )
2012-02-16 00:28:38 +05:30
c = 0_pInt
2009-10-12 21:31:49 +05:30
inPart = . false .
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'part' ) inPart = . false .
2009-10-12 21:31:49 +05:30
2010-07-13 15:56:07 +05:30
if ( ( inPart . or . noPart ) . and . &
2012-02-10 17:26:05 +05:30
IO_lc ( IO_StringValue ( line , myPos , 1_pInt ) ) == '*solid' . and . &
IO_lc ( IO_StringValue ( line , myPos , 2_pInt ) ) == 'section' ) then
2009-10-12 21:31:49 +05:30
elemSetName = ''
materialName = ''
2012-02-16 00:28:38 +05:30
do i = 3_pInt , myPos ( 1_pInt )
2011-10-18 14:50:29 +05:30
if ( IO_extractValue ( IO_lc ( IO_stringValue ( line , myPos , i ) ) , 'elset' ) / = '' ) &
2012-02-10 17:26:05 +05:30
elemSetName = trim ( IO_extractValue ( IO_lc ( IO_stringValue ( line , myPos , i ) ) , 'elset' ) )
2011-10-18 14:50:29 +05:30
if ( IO_extractValue ( IO_lc ( IO_stringValue ( line , myPos , i ) ) , 'material' ) / = '' ) &
2012-02-10 17:26:05 +05:30
materialName = trim ( IO_extractValue ( IO_lc ( IO_stringValue ( line , myPos , i ) ) , 'material' ) )
2009-10-12 21:31:49 +05:30
enddo
if ( elemSetName / = '' . and . materialName / = '' ) then
2012-02-16 00:28:38 +05:30
c = c + 1_pInt
mesh_nameMaterial ( c ) = materialName ! name of material used for this section
mesh_mapMaterial ( c ) = elemSetName ! mapped to respective element set
2009-10-12 21:31:49 +05:30
endif
endif
enddo
2012-02-16 00:28:38 +05:30
620 if ( c == 0_pInt ) call IO_error ( error_ID = 905_pInt )
do i = 1_pInt , c
2010-07-13 15:56:07 +05:30
! write(6,*)'name of materials: ',i,mesh_nameMaterial(i)
! write(6,*)'name of elemSets: ',i,mesh_mapMaterial(i)
2012-02-10 17:26:05 +05:30
if ( mesh_nameMaterial ( i ) == '' . or . mesh_mapMaterial ( i ) == '' ) call IO_error ( error_ID = 905_pInt )
2010-07-13 15:56:07 +05:30
enddo
2009-10-12 21:31:49 +05:30
endsubroutine
2010-05-06 22:10:47 +05:30
!********************************************************************
! map nodes from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPnode
!********************************************************************
subroutine mesh_spectral_map_nodes ( )
use prec , only : pInt
implicit none
integer ( pInt ) i
2012-02-10 17:26:05 +05:30
allocate ( mesh_mapFEtoCPnode ( 2_pInt , mesh_Nnodes ) ) ; mesh_mapFEtoCPnode = 0_pInt
2010-05-06 22:10:47 +05:30
2012-02-16 00:28:38 +05:30
forall ( i = 1_pInt : mesh_Nnodes ) &
2011-09-13 21:24:06 +05:30
mesh_mapFEtoCPnode ( 1 : 2 , i ) = i
2010-07-13 15:56:07 +05:30
2010-05-06 22:10:47 +05:30
endsubroutine
!********************************************************************
2009-10-12 21:31:49 +05:30
! map nodes from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPnode
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_marc_map_nodes ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use math , only : qsort
use IO
implicit none
2012-02-10 17:26:05 +05:30
integer ( pInt ) , parameter :: maxNchunks = 1_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
integer ( pInt ) , dimension ( mesh_Nnodes ) :: node_count
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit , i
2009-10-12 21:31:49 +05:30
2012-02-10 17:26:05 +05:30
allocate ( mesh_mapFEtoCPnode ( 2_pInt , mesh_Nnodes ) ) ; mesh_mapFEtoCPnode = 0_pInt
2009-10-12 21:31:49 +05:30
610 FORMAT ( A300 )
node_count = 0_pInt
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 650 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == 'coordinates' ) then
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 650 ) line ! skip crap line
2012-02-16 00:28:38 +05:30
do i = 1_pInt , mesh_Nnodes
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 650 ) line
2012-02-10 17:26:05 +05:30
mesh_mapFEtoCPnode ( 1_pInt , i ) = IO_fixedIntValue ( line , [ 0_pInt , 10_pInt ] , 1_pInt )
mesh_mapFEtoCPnode ( 2_pInt , i ) = i
2009-10-12 21:31:49 +05:30
enddo
exit
endif
enddo
2012-02-10 17:26:05 +05:30
650 call qsort ( mesh_mapFEtoCPnode , 1_pInt , int ( size ( mesh_mapFEtoCPnode , 2_pInt ) , pInt ) )
2009-10-12 21:31:49 +05:30
endsubroutine
!********************************************************************
! map nodes from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPnode
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_map_nodes ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use math , only : qsort
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 2_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2012-02-16 00:28:38 +05:30
integer ( pInt ) myUnit , i , c , cpNode
2009-10-12 21:31:49 +05:30
logical inPart
2012-02-10 17:26:05 +05:30
allocate ( mesh_mapFEtoCPnode ( 2_pInt , mesh_Nnodes ) ) ; mesh_mapFEtoCPnode = 0_pInt
2009-10-12 21:31:49 +05:30
610 FORMAT ( A300 )
cpNode = 0_pInt
inPart = . false .
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 650 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'part' ) inPart = . false .
2009-10-12 21:31:49 +05:30
2010-07-13 15:56:07 +05:30
if ( ( inPart . or . noPart ) . and . &
2012-02-10 17:26:05 +05:30
IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*node' . and . &
( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'print' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'file' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'response' ) &
2009-10-12 21:31:49 +05:30
) then
2012-02-16 00:28:38 +05:30
c = IO_countDataLines ( myUnit )
do i = 1_pInt , c
2011-10-18 14:50:29 +05:30
backspace ( myUnit )
2009-10-12 21:31:49 +05:30
enddo
2012-02-16 00:28:38 +05:30
do i = 1_pInt , c
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 650 ) line
myPos = IO_stringPos ( line , maxNchunks )
2009-10-12 21:31:49 +05:30
cpNode = cpNode + 1_pInt
2012-02-10 17:26:05 +05:30
mesh_mapFEtoCPnode ( 1_pInt , cpNode ) = IO_intValue ( line , myPos , 1_pInt )
mesh_mapFEtoCPnode ( 2_pInt , cpNode ) = cpNode
2009-10-12 21:31:49 +05:30
enddo
endif
enddo
2012-02-10 17:26:05 +05:30
650 call qsort ( mesh_mapFEtoCPnode , 1_pInt , int ( size ( mesh_mapFEtoCPnode , 2_pInt ) , pInt ) )
2009-10-12 21:31:49 +05:30
2012-02-10 17:26:05 +05:30
if ( int ( size ( mesh_mapFEtoCPnode ) , pInt ) == 0_pInt ) call IO_error ( error_ID = 908_pInt )
2009-10-12 21:31:49 +05:30
endsubroutine
2010-05-06 22:10:47 +05:30
!********************************************************************
! map elements from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPelem
!********************************************************************
subroutine mesh_spectral_map_elements ( )
use prec , only : pInt
implicit none
integer ( pInt ) i
2012-02-10 17:26:05 +05:30
allocate ( mesh_mapFEtoCPelem ( 2_pInt , mesh_NcpElems ) ) ; mesh_mapFEtoCPelem = 0_pInt
2010-05-06 22:10:47 +05:30
2012-02-16 00:28:38 +05:30
forall ( i = 1_pInt : mesh_NcpElems ) &
2011-09-13 21:24:06 +05:30
mesh_mapFEtoCPelem ( 1 : 2 , i ) = i
2010-07-01 20:50:39 +05:30
2010-05-06 22:10:47 +05:30
endsubroutine
2009-10-12 21:31:49 +05:30
!********************************************************************
! map elements from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPelem
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_marc_map_elements ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use math , only : qsort
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 1_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2012-02-16 00:28:38 +05:30
integer ( pInt ) , dimension ( 1_pInt + mesh_NcpElems ) :: contInts
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit , i , cpElem
2009-10-12 21:31:49 +05:30
allocate ( mesh_mapFEtoCPelem ( 2 , mesh_NcpElems ) ) ; mesh_mapFEtoCPelem = 0_pInt
610 FORMAT ( A300 )
cpElem = 0_pInt
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 660 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == 'hypoelastic' ) then
2012-02-16 00:28:38 +05:30
do i = 1_pInt , 3_pInt + hypoelasticTableStyle ! skip three (or four if new table style!) lines
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 660 ) line
2009-10-12 21:31:49 +05:30
enddo
2011-10-18 14:50:29 +05:30
contInts = IO_continousIntValues ( myUnit , mesh_NcpElems , mesh_nameElemSet , mesh_mapElemSet , mesh_NelemSets )
2012-02-16 00:28:38 +05:30
do i = 1_pInt , contInts ( 1 )
cpElem = cpElem + 1_pInt
mesh_mapFEtoCPelem ( 1 , cpElem ) = contInts ( 1_pInt + i )
2009-10-12 21:31:49 +05:30
mesh_mapFEtoCPelem ( 2 , cpElem ) = cpElem
enddo
endif
enddo
2012-02-10 17:26:05 +05:30
660 call qsort ( mesh_mapFEtoCPelem , 1_pInt , int ( size ( mesh_mapFEtoCPelem , 2_pInt ) , pInt ) ) ! should be mesh_NcpElems
2009-10-12 21:31:49 +05:30
endsubroutine
!********************************************************************
! map elements from FE id to internal (consecutive) representation
!
! allocate globals: mesh_mapFEtoCPelem
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_map_elements ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use math , only : qsort
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 2_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit , i , j , k , cpElem
2009-10-12 21:31:49 +05:30
logical materialFound
2012-02-10 17:26:05 +05:30
character ( len = 64 ) materialName , elemSetName ! why limited to 64? ABAQUS?
2009-10-12 21:31:49 +05:30
allocate ( mesh_mapFEtoCPelem ( 2 , mesh_NcpElems ) ) ; mesh_mapFEtoCPelem = 0_pInt
610 FORMAT ( A300 )
cpElem = 0_pInt
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 660 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
select case ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) )
2009-10-12 21:31:49 +05:30
case ( '*material' )
2012-02-10 17:26:05 +05:30
materialName = trim ( IO_extractValue ( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) , 'name' ) ) ! extract name=value
2011-03-23 21:50:12 +05:30
materialFound = materialName / = '' ! valid name?
2009-10-12 21:31:49 +05:30
case ( '*user' )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'material' . and . materialFound ) then
2012-02-16 00:28:38 +05:30
do i = 1_pInt , mesh_Nmaterials ! look thru material names
2011-03-23 21:50:12 +05:30
if ( materialName == mesh_nameMaterial ( i ) ) then ! found one
elemSetName = mesh_mapMaterial ( i ) ! take corresponding elemSet
2012-02-16 00:28:38 +05:30
do k = 1_pInt , mesh_NelemSets ! look thru all elemSet definitions
2011-03-23 21:50:12 +05:30
if ( elemSetName == mesh_nameElemSet ( k ) ) then ! matched?
2012-02-16 00:28:38 +05:30
do j = 1_pInt , mesh_mapElemSet ( 1 , k )
2011-03-23 21:50:12 +05:30
cpElem = cpElem + 1_pInt
2012-02-16 00:28:38 +05:30
mesh_mapFEtoCPelem ( 1 , cpElem ) = mesh_mapElemSet ( 1_pInt + j , k ) ! store FE id
2011-03-23 21:50:12 +05:30
mesh_mapFEtoCPelem ( 2 , cpElem ) = cpElem ! store our id
enddo
endif
enddo
endif
2009-10-12 21:31:49 +05:30
enddo
materialFound = . false .
endif
endselect
enddo
2012-02-10 17:26:05 +05:30
660 call qsort ( mesh_mapFEtoCPelem , 1_pInt , int ( size ( mesh_mapFEtoCPelem , 2_pInt ) , pInt ) ) ! should be mesh_NcpElems
2009-10-12 21:31:49 +05:30
2012-02-10 17:26:05 +05:30
if ( int ( size ( mesh_mapFEtoCPelem ) , pInt ) < 2_pInt ) call IO_error ( error_ID = 907_pInt )
2011-09-13 21:24:06 +05:30
2009-10-12 21:31:49 +05:30
endsubroutine
2010-05-06 22:10:47 +05:30
!********************************************************************
! get maximum count of nodes, IPs, IP neighbors, and subNodes
! among cpElements
!
! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes
!********************************************************************
subroutine mesh_spectral_count_cpSizes ( )
use prec , only : pInt
implicit none
integer ( pInt ) t
t = FE_mapElemtype ( 'C3D8R' ) ! fake 3D hexahedral 8 node 1 IP element
mesh_maxNnodes = FE_Nnodes ( t )
mesh_maxNips = FE_Nips ( t )
mesh_maxNipNeighbors = FE_NipNeighbors ( t )
mesh_maxNsubNodes = FE_NsubNodes ( t )
2010-07-01 20:50:39 +05:30
2010-05-06 22:10:47 +05:30
endsubroutine
2009-10-12 21:31:49 +05:30
!********************************************************************
! get maximum count of nodes, IPs, IP neighbors, and subNodes
! among cpElements
!
! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_marc_count_cpSizes ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 2_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit , i , t , e
2009-10-12 21:31:49 +05:30
mesh_maxNnodes = 0_pInt
mesh_maxNips = 0_pInt
mesh_maxNipNeighbors = 0_pInt
mesh_maxNsubNodes = 0_pInt
610 FORMAT ( A300 )
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 630 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == 'connectivity' ) then
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 630 ) line ! Garbage line
2012-02-16 00:28:38 +05:30
do i = 1_pInt , mesh_Nelems ! read all elements
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 630 ) line
myPos = IO_stringPos ( line , maxNchunks ) ! limit to id and type
2012-02-10 17:26:05 +05:30
e = mesh_FEasCP ( 'elem' , IO_intValue ( line , myPos , 1_pInt ) )
2012-02-16 00:28:38 +05:30
if ( e / = 0_pInt ) then
2012-02-10 17:26:05 +05:30
t = FE_mapElemtype ( IO_stringValue ( line , myPos , 2_pInt ) )
2009-10-12 21:31:49 +05:30
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 ) )
2012-02-10 17:26:05 +05:30
call IO_skipChunks ( myUnit , FE_NoriginalNodes ( t ) - ( myPos ( 1_pInt ) - 2_pInt ) ) ! read on if FE_Nnodes exceeds node count present on current line
2009-10-12 21:31:49 +05:30
endif
enddo
exit
endif
enddo
2011-09-13 21:24:06 +05:30
630 endsubroutine
2009-10-12 21:31:49 +05:30
!********************************************************************
! get maximum count of nodes, IPs, IP neighbors, and subNodes
! among cpElements
!
! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_count_cpSizes ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 2_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2012-02-16 00:28:38 +05:30
integer ( pInt ) myUnit , i , c , t
2009-10-12 21:31:49 +05:30
logical inPart
mesh_maxNnodes = 0_pInt
mesh_maxNips = 0_pInt
mesh_maxNipNeighbors = 0_pInt
mesh_maxNsubNodes = 0_pInt
610 FORMAT ( A300 )
inPart = . false .
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'part' ) inPart = . false .
2009-10-12 21:31:49 +05:30
2010-07-13 15:56:07 +05:30
if ( ( inPart . or . noPart ) . and . &
2012-02-10 17:26:05 +05:30
IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*element' . and . &
( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'matrix' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'response' ) &
2009-10-12 21:31:49 +05:30
) then
2012-02-10 17:26:05 +05:30
t = FE_mapElemtype ( IO_extractValue ( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) , 'type' ) ) ! remember elem type
2012-02-16 00:28:38 +05:30
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
2011-10-18 14:50:29 +05:30
backspace ( myUnit )
2009-10-12 21:31:49 +05:30
enddo
2012-02-16 00:28:38 +05:30
do i = 1_pInt , c
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks ) ! limit to 64 nodes max
2012-02-16 00:28:38 +05:30
if ( mesh_FEasCP ( 'elem' , IO_intValue ( line , myPos , 1_pInt ) ) / = 0_pInt ) then ! disregard non CP elems
2009-10-12 21:31:49 +05:30
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
endif
enddo
2011-09-13 21:24:06 +05:30
620 endsubroutine
2009-10-12 21:31:49 +05:30
2010-05-06 22:10:47 +05:30
!********************************************************************
! store x,y,z coordinates of all nodes in mesh
!
! allocate globals:
! _node
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_spectral_build_nodes ( myUnit )
2010-05-06 22:10:47 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 7_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2011-10-18 14:50:29 +05:30
integer ( pInt ) a , b , c , n , i , j , headerLength
2010-05-06 22:10:47 +05:30
real ( pReal ) x , y , z
logical gotResolution , gotDimension
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit
character ( len = 1024 ) line , keyword
2010-05-06 22:10:47 +05:30
2011-07-31 21:12:59 +05:30
allocate ( mesh_node0 ( 3 , mesh_Nnodes ) ) ; mesh_node0 = 0.0_pReal
allocate ( mesh_node ( 3 , mesh_Nnodes ) ) ; mesh_node = 0.0_pReal
2010-07-01 20:50:39 +05:30
2010-05-06 22:10:47 +05:30
a = 1_pInt
b = 1_pInt
c = 1_pInt
x = 1.0_pReal
y = 1.0_pReal
z = 1.0_pReal
gotResolution = . false .
gotDimension = . false .
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
read ( myUnit , '(a1024)' ) line
2012-02-10 17:26:05 +05:30
myPos = IO_stringPos ( line , 2_pInt )
keyword = IO_lc ( IO_StringValue ( line , myPos , 2_pInt ) )
2011-10-18 14:50:29 +05:30
if ( keyword ( 1 : 4 ) == 'head' ) then
2012-02-10 17:26:05 +05:30
headerLength = IO_intValue ( line , myPos , 1_pInt ) + 1_pInt
2011-10-18 14:50:29 +05:30
else
2012-02-13 23:11:27 +05:30
call IO_error ( error_ID = 842_pInt )
2011-10-18 14:50:29 +05:30
endif
rewind ( myUnit )
2012-02-16 00:28:38 +05:30
do i = 1_pInt , headerLength
2011-10-18 14:50:29 +05:30
read ( myUnit , '(a1024)' ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
select case ( IO_lc ( IO_StringValue ( line , myPos , 1_pInt ) ) )
2010-05-06 22:10:47 +05:30
case ( 'dimension' )
2011-10-18 14:50:29 +05:30
gotDimension = . true .
2012-02-16 00:28:38 +05:30
do j = 2_pInt , 6_pInt , 2_pInt
2011-10-18 14:50:29 +05:30
select case ( IO_lc ( IO_stringValue ( line , myPos , j ) ) )
case ( 'x' )
2012-02-10 17:26:05 +05:30
x = IO_floatValue ( line , myPos , j + 1_pInt )
2011-10-18 14:50:29 +05:30
case ( 'y' )
2012-02-10 17:26:05 +05:30
y = IO_floatValue ( line , myPos , j + 1_pInt )
2011-10-18 14:50:29 +05:30
case ( 'z' )
2012-02-10 17:26:05 +05:30
z = IO_floatValue ( line , myPos , j + 1_pInt )
2011-10-18 14:50:29 +05:30
end select
enddo
case ( 'resolution' )
gotResolution = . true .
2012-02-16 00:28:38 +05:30
do j = 2_pInt , 6_pInt , 2_pInt
2011-10-18 14:50:29 +05:30
select case ( IO_lc ( IO_stringValue ( line , myPos , j ) ) )
case ( 'a' )
2012-02-10 17:26:05 +05:30
a = 1_pInt + IO_intValue ( line , myPos , j + 1_pInt )
2011-10-18 14:50:29 +05:30
case ( 'b' )
2012-02-10 17:26:05 +05:30
b = 1_pInt + IO_intValue ( line , myPos , j + 1_pInt )
2011-10-18 14:50:29 +05:30
case ( 'c' )
2012-02-10 17:26:05 +05:30
c = 1_pInt + IO_intValue ( line , myPos , j + 1_pInt )
2011-10-18 14:50:29 +05:30
end select
enddo
2010-05-10 20:24:59 +05:30
end select
2010-05-06 22:10:47 +05:30
enddo
! --- sanity checks ---
2012-02-13 23:11:27 +05:30
if ( ( . not . gotDimension ) . or . ( . not . gotResolution ) ) call IO_error ( error_ID = 842_pInt )
2012-02-16 00:28:38 +05:30
if ( ( a < 1_pInt ) . or . ( b < 1_pInt ) . or . ( c < 0_pInt ) ) call IO_error ( error_ID = 843_pInt ) ! 1_pInt is already added
2012-02-13 23:11:27 +05:30
if ( ( x < = 0.0_pReal ) . or . ( y < = 0.0_pReal ) . or . ( z < = 0.0_pReal ) ) call IO_error ( error_ID = 844_pInt )
2010-05-06 22:10:47 +05:30
2012-02-16 00:28:38 +05:30
forall ( n = 0_pInt : mesh_Nnodes - 1_pInt )
mesh_node0 ( 1 , n + 1_pInt ) = x * real ( mod ( n , a ) , pReal ) / real ( a - 1_pInt , pReal )
mesh_node0 ( 2 , n + 1_pInt ) = y * real ( mod ( n / a , b ) , pReal ) / real ( b - 1_pInt , pReal )
mesh_node0 ( 3 , n + 1_pInt ) = z * real ( mod ( n / a / b , c ) , pReal ) / real ( c - 1_pInt , pReal )
2012-01-17 22:23:56 +05:30
end forall
2011-10-18 14:50:29 +05:30
mesh_node = mesh_node0 !why?
2010-05-06 22:10:47 +05:30
2011-10-18 14:50:29 +05:30
endsubroutine
2010-05-06 22:10:47 +05:30
2009-10-12 21:31:49 +05:30
!********************************************************************
! store x,y,z coordinates of all nodes in mesh
!
! allocate globals:
! _node
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_marc_build_nodes ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , dimension ( 5 ) , parameter :: node_ends = int ( [ 0 , 10 , 30 , 50 , 70 ] , pInt )
integer ( pInt ) , parameter :: maxNchunks = 1_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit , i , j , m
2009-10-12 21:31:49 +05:30
2011-07-31 21:12:59 +05:30
allocate ( mesh_node0 ( 3 , mesh_Nnodes ) ) ; mesh_node0 = 0.0_pReal
allocate ( mesh_node ( 3 , mesh_Nnodes ) ) ; mesh_node = 0.0_pReal
2009-10-12 21:31:49 +05:30
610 FORMAT ( A300 )
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 670 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == 'coordinates' ) then
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 670 ) line ! skip crap line
2012-02-16 00:28:38 +05:30
do i = 1_pInt , mesh_Nnodes
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 670 ) line
2012-02-10 17:26:05 +05:30
m = mesh_FEasCP ( 'node' , IO_fixedIntValue ( line , node_ends , 1_pInt ) )
forall ( j = 1_pInt : 3_pInt ) mesh_node0 ( j , m ) = IO_fixedNoEFloatValue ( line , node_ends , j + 1_pInt )
2009-10-12 21:31:49 +05:30
enddo
exit
endif
enddo
2011-05-24 21:27:59 +05:30
670 mesh_node = mesh_node0
2009-10-12 21:31:49 +05:30
endsubroutine
!********************************************************************
! store x,y,z coordinates of all nodes in mesh
!
! allocate globals:
! _node
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_build_nodes ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 4_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2012-02-16 00:28:38 +05:30
integer ( pInt ) myUnit , i , j , m , c
2009-10-12 21:31:49 +05:30
logical inPart
2011-07-31 21:12:59 +05:30
allocate ( mesh_node0 ( 3 , mesh_Nnodes ) ) ; mesh_node0 = 0.0_pReal
allocate ( mesh_node ( 3 , mesh_Nnodes ) ) ; mesh_node = 0.0_pReal
2009-10-12 21:31:49 +05:30
610 FORMAT ( A300 )
inPart = . false .
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 670 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'part' ) inPart = . false .
2009-10-12 21:31:49 +05:30
2010-07-13 15:56:07 +05:30
if ( ( inPart . or . noPart ) . and . &
2012-02-10 17:26:05 +05:30
IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*node' . and . &
( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'print' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'file' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'response' ) &
2009-10-12 21:31:49 +05:30
) then
2012-02-16 00:28:38 +05:30
c = IO_countDataLines ( myUnit ) ! how many nodes are defined here?
do i = 1_pInt , c
2011-10-18 14:50:29 +05:30
backspace ( myUnit ) ! rewind to first entry
2009-10-12 21:31:49 +05:30
enddo
2012-02-16 00:28:38 +05:30
do i = 1_pInt , c
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 670 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
m = mesh_FEasCP ( 'node' , IO_intValue ( line , myPos , 1_pInt ) )
2012-02-16 00:28:38 +05:30
forall ( j = 1_pInt : 3_pInt ) mesh_node0 ( j , m ) = IO_floatValue ( line , myPos , j + 1_pInt )
2009-10-12 21:31:49 +05:30
enddo
endif
enddo
2012-02-10 17:26:05 +05:30
670 if ( int ( size ( mesh_node0 , 2_pInt ) , pInt ) / = mesh_Nnodes ) call IO_error ( error_ID = 909_pInt )
2011-05-24 21:27:59 +05:30
mesh_node = mesh_node0
2009-10-12 21:31:49 +05:30
endsubroutine
2010-05-06 22:10:47 +05:30
!********************************************************************
! store FEid, type, mat, tex, and node list per element
!
! allocate globals:
! _element
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_spectral_build_elements ( myUnit )
2010-05-06 22:10:47 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 7_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2012-01-12 22:31:24 +05:30
integer ( pInt ) a , b , c , e , i , j , homog , headerLength , maxIntCount
integer ( pInt ) , dimension ( : ) , allocatable :: microstructures
integer ( pInt ) , dimension ( 1 , 1 ) :: dummySet = 0_pInt
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit
2012-01-12 22:31:24 +05:30
character ( len = 65536 ) line , keyword
character ( len = 64 ) , dimension ( 1 ) :: dummyName = ''
2010-05-06 22:10:47 +05:30
a = 1_pInt
b = 1_pInt
c = 1_pInt
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2012-01-12 22:31:24 +05:30
read ( myUnit , '(a65536)' ) line
2012-02-10 17:26:05 +05:30
myPos = IO_stringPos ( line , 2_pInt )
keyword = IO_lc ( IO_StringValue ( line , myPos , 2_pInt ) )
2011-10-18 14:50:29 +05:30
if ( keyword ( 1 : 4 ) == 'head' ) then
2012-02-10 17:26:05 +05:30
headerLength = IO_intValue ( line , myPos , 1_pInt ) + 1_pInt
2011-10-18 14:50:29 +05:30
else
2012-02-13 23:11:27 +05:30
call IO_error ( error_ID = 842_pInt )
2011-10-18 14:50:29 +05:30
endif
rewind ( myUnit )
2012-02-16 00:28:38 +05:30
do i = 1_pInt , headerLength
2012-01-12 22:31:24 +05:30
read ( myUnit , '(a65536)' ) line
2011-10-18 14:50:29 +05:30
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
select case ( IO_lc ( IO_StringValue ( line , myPos , 1_pInt ) ) )
2010-05-06 22:10:47 +05:30
case ( 'resolution' )
2012-02-16 00:28:38 +05:30
do j = 2_pInt , 6_pInt , 2_pInt
2011-10-18 14:50:29 +05:30
select case ( IO_lc ( IO_stringValue ( line , myPos , j ) ) )
case ( 'a' )
2012-02-10 17:26:05 +05:30
a = 1_pInt + IO_intValue ( line , myPos , j + 1_pInt )
2011-10-18 14:50:29 +05:30
case ( 'b' )
2012-02-10 17:26:05 +05:30
b = 1_pInt + IO_intValue ( line , myPos , j + 1_pInt )
2011-10-18 14:50:29 +05:30
case ( 'c' )
2012-02-10 17:26:05 +05:30
c = 1_pInt + IO_intValue ( line , myPos , j + 1_pInt )
2011-10-18 14:50:29 +05:30
end select
enddo
case ( 'homogenization' )
2012-02-10 17:26:05 +05:30
homog = IO_intValue ( line , myPos , 2_pInt )
2010-05-10 20:24:59 +05:30
end select
2010-05-06 22:10:47 +05:30
enddo
2012-01-12 22:31:24 +05:30
maxIntCount = 0_pInt
i = 1_pInt
do while ( i > 0_pInt )
i = IO_countContinousIntValues ( myUnit )
maxIntCount = max ( maxIntCount , i )
enddo
rewind ( myUnit )
2012-02-16 00:28:38 +05:30
do i = 1_pInt , headerLength ! skip header
2012-01-12 22:31:24 +05:30
read ( myUnit , '(a65536)' ) line
enddo
2012-02-16 00:28:38 +05:30
100 allocate ( mesh_element ( 4_pInt + mesh_maxNnodes , mesh_NcpElems ) ) ; mesh_element = 0_pInt
2012-01-12 22:31:24 +05:30
allocate ( microstructures ( 1_pInt + maxIntCount ) ) ; microstructures = 2_pInt
2010-05-06 22:10:47 +05:30
e = 0_pInt
2012-01-12 22:31:24 +05:30
do while ( e < mesh_NcpElems . and . microstructures ( 1 ) > 0_pInt ) ! fill expected number of elements, stop at end of data (or blank line!)
2012-02-10 17:26:05 +05:30
microstructures = IO_continousIntValues ( myUnit , maxIntCount , dummyName , dummySet , 0_pInt ) ! get affected elements
2012-02-16 00:28:38 +05:30
do i = 1_pInt , microstructures ( 1_pInt )
e = e + 1_pInt ! valid element entry
2012-01-17 22:23:56 +05:30
mesh_element ( 1 , e ) = e ! FE id
mesh_element ( 2 , e ) = FE_mapElemtype ( 'C3D8R' ) ! elem type
mesh_element ( 3 , e ) = homog ! homogenization
mesh_element ( 4 , e ) = microstructures ( 1_pInt + i ) ! microstructure
2012-02-10 17:26:05 +05:30
mesh_element ( 5 , e ) = e + ( e - 1_pInt ) / ( a - 1_pInt ) + ( ( e - 1_pInt ) / ( ( a - 1_pInt ) * ( b - 1_pInt ) ) ) * a ! base node
mesh_element ( 6 , e ) = mesh_element ( 5 , e ) + 1_pInt
mesh_element ( 7 , e ) = mesh_element ( 5 , e ) + a + 1_pInt
2012-01-17 22:23:56 +05:30
mesh_element ( 8 , e ) = mesh_element ( 5 , e ) + a
mesh_element ( 9 , e ) = mesh_element ( 5 , e ) + a * b ! second floor base node
2012-02-10 17:26:05 +05:30
mesh_element ( 10 , e ) = mesh_element ( 9 , e ) + 1_pInt
mesh_element ( 11 , e ) = mesh_element ( 9 , e ) + a + 1_pInt
2012-01-17 22:23:56 +05:30
mesh_element ( 12 , e ) = mesh_element ( 9 , e ) + a
2012-01-12 22:31:24 +05:30
mesh_maxValStateVar ( 1 ) = max ( mesh_maxValStateVar ( 1 ) , mesh_element ( 3 , e ) ) !needed for statistics
mesh_maxValStateVar ( 2 ) = max ( mesh_maxValStateVar ( 2 ) , mesh_element ( 4 , e ) )
enddo
2010-05-06 22:10:47 +05:30
enddo
2010-07-01 20:50:39 +05:30
2012-01-12 22:31:24 +05:30
110 deallocate ( microstructures )
2012-02-13 23:11:27 +05:30
if ( e / = mesh_NcpElems ) call IO_error ( 880_pInt , e )
2012-01-12 22:31:24 +05:30
endsubroutine
2010-05-06 22:10:47 +05:30
2009-10-12 21:31:49 +05:30
!********************************************************************
! store FEid, type, mat, tex, and node list per element
!
! allocate globals:
! _element
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_marc_build_elements ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 66_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
character ( len = 300 ) line
2012-02-16 00:28:38 +05:30
integer ( pInt ) , dimension ( 1_pInt + mesh_NcpElems ) :: contInts
2011-10-18 14:50:29 +05:30
integer ( pInt ) myUnit , i , j , sv , val , e
2009-10-12 21:31:49 +05:30
2012-02-16 00:28:38 +05:30
allocate ( mesh_element ( 4_pInt + mesh_maxNnodes , mesh_NcpElems ) ) ; mesh_element = 0_pInt
2009-10-12 21:31:49 +05:30
610 FORMAT ( A300 )
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
2012-02-10 17:26:05 +05:30
myPos ( 1 : 1 + 2 * 1 ) = IO_stringPos ( line , 1_pInt )
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == 'connectivity' ) then
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line ! Garbage line
2012-02-16 00:28:38 +05:30
do i = 1_pInt , mesh_Nelems
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks ) ! limit to 64 nodes max (plus ID, type)
2012-02-10 17:26:05 +05:30
e = mesh_FEasCP ( 'elem' , IO_intValue ( line , myPos , 1_pInt ) )
2012-02-16 00:28:38 +05:30
if ( e / = 0_pInt ) then ! disregard non CP elems
2012-02-10 17:26:05 +05:30
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
2012-02-16 00:28:38 +05:30
forall ( j = 1_pInt : FE_Nnodes ( mesh_element ( 2 , e ) ) ) &
mesh_element ( j + 4_pInt , e ) = IO_IntValue ( line , myPos , j + 2_pInt ) ! copy FE ids of nodes
2012-02-10 17:26:05 +05:30
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
2009-10-12 21:31:49 +05:30
endif
enddo
exit
endif
enddo
2011-10-18 14:50:29 +05:30
620 rewind ( myUnit ) ! just in case "initial state" apears before "connectivity"
read ( myUnit , 610 , END = 620 ) line
2009-10-12 21:31:49 +05:30
do
2012-02-10 17:26:05 +05:30
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
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 630 ) line ! read line with index of state var
2012-02-10 17:26:05 +05:30
myPos ( 1 : 1 + 2 * 1 ) = IO_stringPos ( line , 1_pInt )
sv = IO_IntValue ( line , myPos , 1_pInt ) ! figure state variable index
if ( ( sv == 2_pInt ) . or . ( sv == 3_pInt ) ) then ! only state vars 2 and 3 of interest
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line ! read line with value of state var
2012-02-10 17:26:05 +05:30
myPos ( 1 : 1 + 2 * 1 ) = IO_stringPos ( line , 1_pInt )
2012-02-16 00:28:38 +05:30
do while ( scan ( IO_stringValue ( line , myPos , 1_pInt ) , '+-' , back = . true . ) > 1 ) ! is noEfloat value?
val = nint ( IO_fixedNoEFloatValue ( line , [ 0_pInt , 20_pInt ] , 1_pInt ) , pInt ) ! state var's value
mesh_maxValStateVar ( sv - 1_pInt ) = max ( val , mesh_maxValStateVar ( sv - 1_pInt ) ) ! remember max val of homogenization and microstructure index
if ( initialcondTableStyle == 2_pInt ) then
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 630 ) line ! read extra line
read ( myUnit , 610 , END = 630 ) line ! read extra line
2009-10-12 21:31:49 +05:30
endif
2011-10-18 14:50:29 +05:30
contInts = IO_continousIntValues ( myUnit , mesh_Nelems , mesh_nameElemSet , mesh_mapElemSet , mesh_NelemSets ) ! get affected elements
2012-02-16 00:28:38 +05:30
do i = 1_pInt , contInts ( 1 )
2012-02-10 17:26:05 +05:30
e = mesh_FEasCP ( 'elem' , contInts ( 1_pInt + i ) )
mesh_element ( 1_pInt + sv , e ) = val
2009-10-12 21:31:49 +05:30
enddo
2012-02-16 00:28:38 +05:30
if ( initialcondTableStyle == 0_pInt ) read ( myUnit , 610 , END = 620 ) line ! ignore IP range for old table style
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 630 ) line
2012-02-10 17:26:05 +05:30
myPos ( 1 : 1 + 2 * 1 ) = IO_stringPos ( line , 1_pInt )
2009-10-12 21:31:49 +05:30
enddo
endif
else
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 630 ) line
2009-10-12 21:31:49 +05:30
endif
enddo
2011-09-13 21:24:06 +05:30
630 endsubroutine
2009-10-12 21:31:49 +05:30
!********************************************************************
! store FEid, type, mat, tex, and node list per element
!
! allocate globals:
! _element
!********************************************************************
2011-10-18 14:50:29 +05:30
subroutine mesh_abaqus_build_elements ( myUnit )
2009-10-12 21:31:49 +05:30
use prec , only : pInt
use IO
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: maxNchunks = 65_pInt
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * maxNchunks ) :: myPos
2009-10-12 21:31:49 +05:30
2012-02-16 00:28:38 +05:30
integer ( pInt ) myUnit , i , j , k , c , e , t , homog , micro
2009-10-12 21:31:49 +05:30
logical inPart , materialFound
character ( len = 64 ) materialName , elemSetName
character ( len = 300 ) line
2012-02-16 00:28:38 +05:30
allocate ( mesh_element ( 4_pInt + mesh_maxNnodes , mesh_NcpElems ) ) ; mesh_element = 0_pInt
2009-10-12 21:31:49 +05:30
610 FORMAT ( A300 )
inPart = . false .
2011-10-18 14:50:29 +05:30
rewind ( myUnit )
2009-10-12 21:31:49 +05:30
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
2012-02-10 17:26:05 +05:30
myPos ( 1 : 1 + 2 * 2 ) = IO_stringPos ( line , 2_pInt )
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) == 'part' ) inPart = . false .
2009-10-12 21:31:49 +05:30
2010-07-13 15:56:07 +05:30
if ( ( inPart . or . noPart ) . and . &
2012-02-10 17:26:05 +05:30
IO_lc ( IO_stringValue ( line , myPos , 1_pInt ) ) == '*element' . and . &
( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'matrix' . and . &
IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) / = 'response' ) &
2009-10-12 21:31:49 +05:30
) then
2012-02-10 17:26:05 +05:30
t = FE_mapElemtype ( IO_extractValue ( IO_lc ( IO_stringValue ( line , myPos , 2_pInt ) ) , 'type' ) ) ! remember elem type
2012-02-16 00:28:38 +05:30
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
2011-10-18 14:50:29 +05:30
backspace ( myUnit )
2009-10-12 21:31:49 +05:30
enddo
2012-02-16 00:28:38 +05:30
do i = 1_pInt , c
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 620 ) line
myPos = IO_stringPos ( line , maxNchunks ) ! limit to 64 nodes max
2012-02-10 17:26:05 +05:30
e = mesh_FEasCP ( 'elem' , IO_intValue ( line , myPos , 1_pInt ) )
2012-02-16 00:28:38 +05:30
if ( e / = 0_pInt ) then ! disregard non CP elems
2012-02-10 17:26:05 +05:30
mesh_element ( 1 , e ) = IO_intValue ( line , myPos , 1_pInt ) ! FE id
2009-10-12 21:31:49 +05:30
mesh_element ( 2 , e ) = t ! elem type
2012-02-16 00:28:38 +05:30
forall ( j = 1_pInt : FE_Nnodes ( t ) ) &
2012-02-10 17:26:05 +05:30
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
2009-10-12 21:31:49 +05:30
endif
enddo
endif
enddo
2011-10-18 14:50:29 +05:30
620 rewind ( myUnit ) ! just in case "*material" definitions apear before "*element"
2009-10-12 21:31:49 +05:30
materialFound = . false .
do
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 630 ) line
myPos = IO_stringPos ( line , maxNchunks )
2012-02-10 17:26:05 +05:30
select case ( IO_lc ( IO_StringValue ( line , myPos , 1_pInt ) ) )
2009-10-12 21:31:49 +05:30
case ( '*material' )
2012-02-10 17:26:05 +05:30
materialName = trim ( IO_extractValue ( IO_lc ( IO_StringValue ( line , myPos , 2_pInt ) ) , 'name' ) ) ! extract name=value
2011-03-23 21:50:12 +05:30
materialFound = materialName / = '' ! valid name?
2009-10-12 21:31:49 +05:30
case ( '*user' )
2012-02-10 17:26:05 +05:30
if ( IO_lc ( IO_StringValue ( line , myPos , 2_pInt ) ) == 'material' . and . &
2009-10-12 21:31:49 +05:30
materialFound ) then
2011-10-18 14:50:29 +05:30
read ( myUnit , 610 , END = 630 ) line ! read homogenization and microstructure
2012-02-10 17:26:05 +05:30
myPos ( 1 : 1 + 2 * 2 ) = IO_stringPos ( line , 2_pInt )
2012-02-16 00:28:38 +05:30
homog = nint ( IO_floatValue ( line , myPos , 1_pInt ) , pInt )
micro = nint ( IO_floatValue ( line , myPos , 2_pInt ) , pInt )
do i = 1_pInt , mesh_Nmaterials ! look thru material names
2011-03-23 21:50:12 +05:30
if ( materialName == mesh_nameMaterial ( i ) ) then ! found one
elemSetName = mesh_mapMaterial ( i ) ! take corresponding elemSet
2012-02-16 00:28:38 +05:30
do k = 1_pInt , mesh_NelemSets ! look thru all elemSet definitions
2011-03-23 21:50:12 +05:30
if ( elemSetName == mesh_nameElemSet ( k ) ) then ! matched?
2012-02-16 00:28:38 +05:30
do j = 1_pInt , mesh_mapElemSet ( 1 , k )
2011-03-23 21:50:12 +05:30
e = mesh_FEasCP ( 'elem' , mesh_mapElemSet ( 1 + j , k ) )
mesh_element ( 3 , e ) = homog ! store homogenization
mesh_element ( 4 , e ) = micro ! store microstructure
mesh_maxValStateVar ( 1 ) = max ( mesh_maxValStateVar ( 1 ) , homog )
mesh_maxValStateVar ( 2 ) = max ( mesh_maxValStateVar ( 2 ) , micro )
enddo
endif
enddo
endif
2009-10-12 21:31:49 +05:30
enddo
materialFound = . false .
endif
endselect
enddo
2011-09-13 21:24:06 +05:30
630 endsubroutine
2009-10-12 21:31:49 +05:30
!********************************************************************
! get maximum count of shared elements among cpElements and
! build list of elements shared by each node in mesh
!
! _maxNsharedElems
! _sharedElem
!********************************************************************
2011-02-16 21:53:08 +05:30
subroutine mesh_build_sharedElems ( )
use prec , only : pInt
implicit none
integer ( pint ) e , & ! element index
t , & ! element type
node , & ! CP node index
j , & ! node index per element
2012-02-16 00:28:38 +05:30
myDim , & ! dimension index
2011-02-16 21:53:08 +05:30
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 ) ) )
2011-03-23 21:50:12 +05:30
2011-02-16 21:53:08 +05:30
node_count = 0_pInt
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems
2011-07-31 21:12:59 +05:30
t = mesh_element ( 2 , e ) ! get element type
2011-02-16 21:53:08 +05:30
node_seen = 0_pInt ! reset node duplicates
2012-02-16 00:28:38 +05:30
do j = 1_pInt , FE_Nnodes ( t ) ! check each node of element
2011-07-31 21:12:59 +05:30
node = mesh_FEasCP ( 'node' , mesh_element ( 4 + j , e ) ) ! translate to internal (consecutive) numbering
2011-02-16 21:53:08 +05:30
if ( all ( node_seen / = node ) ) then
node_count ( node ) = node_count ( node ) + 1_pInt ! if FE node not yet encountered -> count it
2012-02-16 00:28:38 +05:30
do myDim = 1_pInt , 3_pInt ! check in each dimension...
nodeTwin = mesh_nodeTwins ( myDim , node )
2011-07-31 21:12:59 +05:30
if ( nodeTwin > 0_pInt ) & ! if I am a twin of some node...
2011-02-16 21:53:08 +05:30
node_count ( nodeTwin ) = node_count ( nodeTwin ) + 1_pInt ! -> count me again for the twin node
enddo
endif
node_seen ( j ) = node ! remember this node to be counted already
enddo
enddo
2012-02-16 00:28:38 +05:30
mesh_maxNsharedElems = int ( maxval ( node_count ) , pInt ) ! most shared node
2011-02-16 21:53:08 +05:30
allocate ( mesh_sharedElem ( 1 + mesh_maxNsharedElems , mesh_Nnodes ) )
mesh_sharedElem = 0_pInt
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems
2011-02-16 21:53:08 +05:30
t = mesh_element ( 2 , e )
node_seen = 0_pInt
2012-02-16 00:28:38 +05:30
do j = 1_pInt , FE_Nnodes ( t )
node = mesh_FEasCP ( 'node' , mesh_element ( 4_pInt + j , e ) )
2011-02-16 21:53:08 +05:30
if ( all ( node_seen / = node ) ) then
mesh_sharedElem ( 1 , node ) = mesh_sharedElem ( 1 , node ) + 1_pInt ! count for each node the connected elements
2012-02-16 00:28:38 +05:30
mesh_sharedElem ( mesh_sharedElem ( 1 , node ) + 1_pInt , node ) = e ! store the respective element id
do myDim = 1_pInt , 3_pInt ! check in each dimension...
nodeTwin = mesh_nodeTwins ( myDim , node )
if ( nodeTwin > 0_pInt ) then ! if i am a twin of some node...
2011-02-16 21:53:08 +05:30
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
2009-10-12 21:31:49 +05:30
2011-02-16 21:53:08 +05:30
deallocate ( node_seen )
2009-10-12 21:31:49 +05:30
2011-02-16 21:53:08 +05:30
endsubroutine
2009-10-12 21:31:49 +05:30
!***********************************************************
! build up of IP neighborhood
!
! allocate globals
! _ipNeighborhood
!***********************************************************
2011-02-16 21:53:08 +05:30
subroutine mesh_build_ipNeighborhood ( )
use prec , only : pInt
implicit none
integer ( pInt ) myElem , & ! my CP element index
myIP , &
myType , & ! my element type
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
2011-04-13 19:46:22 +05:30
a , anchor
2011-02-16 21:53:08 +05:30
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
2012-02-16 00:28:38 +05:30
do myElem = 1_pInt , mesh_NcpElems ! loop over cpElems
2011-02-16 21:53:08 +05:30
myType = mesh_element ( 2 , myElem ) ! get elemType
2012-02-16 00:28:38 +05:30
do myIP = 1_pInt , FE_Nips ( myType ) ! loop over IPs of elem
2011-07-31 21:12:59 +05:30
2012-02-16 00:28:38 +05:30
do neighbor = 1_pInt , FE_NipNeighbors ( myType ) ! loop over neighbors of IP
2011-02-16 21:53:08 +05:30
neighboringIPkey = FE_ipNeighbor ( neighbor , myIP , myType )
2011-07-31 21:12:59 +05:30
2011-02-16 21:53:08 +05:30
!*** if the key is positive, the neighbor is inside the element
!*** that means, we have already found our neighboring IP
2012-02-16 00:28:38 +05:30
if ( neighboringIPkey > 0_pInt ) then
2011-02-16 21:53:08 +05:30
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
2012-02-16 00:28:38 +05:30
elseif ( neighboringIPkey < 0_pInt ) then ! neighboring element's IP
2011-02-16 21:53:08 +05:30
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
2010-05-10 20:24:59 +05:30
2011-02-16 21:53:08 +05:30
!*** find those nodes which build the link to the neighbor
NlinkedNodes = 0_pInt
linkedNodes = 0_pInt
2012-02-16 00:28:38 +05:30
do a = 1_pInt , FE_maxNnodesAtIP ( myType ) ! figure my anchor nodes on connecting face
2011-02-16 21:53:08 +05:30
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?
NlinkedNodes = NlinkedNodes + 1_pInt
2011-07-31 21:12:59 +05:30
linkedNodes ( NlinkedNodes ) = &
2012-02-16 00:28:38 +05:30
mesh_FEasCP ( 'node' , mesh_element ( 4_pInt + anchor , myElem ) ) ! CP id of anchor node
2011-02-16 21:53:08 +05:30
else ! something went wrong with the linkage, since not all anchors sit on my face
NlinkedNodes = 0_pInt
linkedNodes = 0_pInt
exit
endif
endif
enddo
!*** loop through the ips of my neighbor
!*** and try to find an ip with matching nodes
!*** also try to match with node twins
2011-07-31 21:12:59 +05:30
2012-02-16 00:28:38 +05:30
checkCandidateIP : do candidateIP = 1_pInt , FE_Nips ( neighboringType )
2011-02-16 21:53:08 +05:30
NmatchingNodes = 0_pInt
matchingNodes = 0_pInt
2012-02-16 00:28:38 +05:30
do a = 1_pInt , FE_maxNnodesAtIP ( neighboringType ) ! check each anchor node of that ip
2011-02-16 21:53:08 +05:30
anchor = FE_nodesAtIP ( a , candidateIP , neighboringType )
if ( anchor / = 0_pInt ) then ! valid anchor node
if ( any ( FE_nodeOnFace ( : , matchingFace , neighboringType ) == anchor ) ) then ! sits on matching face?
NmatchingNodes = NmatchingNodes + 1_pInt
2011-07-31 21:12:59 +05:30
matchingNodes ( NmatchingNodes ) = &
mesh_FEasCP ( 'node' , mesh_element ( 4 + anchor , matchingElem ) ) ! CP id of neighbor's anchor node
2011-02-16 21:53:08 +05:30
else ! no matching, because not all nodes sit on the matching face
NmatchingNodes = 0_pInt
matchingNodes = 0_pInt
exit
endif
endif
enddo
if ( NmatchingNodes / = NlinkedNodes ) & ! this ip has wrong count of anchors on face
cycle checkCandidateIP
!*** check "normal" nodes whether they match or not
checkTwins = . false .
2012-02-16 00:28:38 +05:30
do a = 1_pInt , NlinkedNodes
2011-02-16 21:53:08 +05:30
if ( all ( matchingNodes / = linkedNodes ( a ) ) ) then ! this linkedNode does not match any matchingNode
checkTwins = . true .
exit ! no need to search further
endif
enddo
!*** if no match found, then also check node twins
if ( checkTwins ) then
2012-02-16 00:28:38 +05:30
dir = int ( maxloc ( abs ( mesh_ipAreaNormal ( 1 : 3 , neighbor , myIP , myElem ) ) , 1 ) , pInt ) ! check for twins only in direction of the surface normal
do a = 1_pInt , NlinkedNodes
2011-02-24 14:56:30 +05:30
twin_of_linkedNode = mesh_nodeTwins ( dir , linkedNodes ( a ) )
2011-07-31 21:12:59 +05:30
if ( twin_of_linkedNode == 0_pInt . or . & ! twin of linkedNode does not exist...
all ( matchingNodes / = twin_of_linkedNode ) ) then ! ... or it does not match any matchingNode
cycle checkCandidateIP ! ... then check next candidateIP
2011-02-24 14:56:30 +05:30
endif
enddo
2011-02-16 21:53:08 +05:30
endif
!*** we found a match !!!
2011-07-31 21:12:59 +05:30
2011-02-16 21:53:08 +05:30
mesh_ipNeighborhood ( 1 , neighbor , myIP , myElem ) = matchingElem
mesh_ipNeighborhood ( 2 , neighbor , myIP , myElem ) = candidateIP
2011-02-24 14:56:30 +05:30
exit checkCandidateIP
2011-02-16 21:53:08 +05:30
enddo checkCandidateIP
2011-07-31 21:12:59 +05:30
endif ! end of valid external matching
endif ! end of internal/external matching
2011-02-16 21:53:08 +05:30
enddo
enddo
enddo
2009-01-20 00:12:31 +05:30
2011-02-16 21:53:08 +05:30
endsubroutine
2009-01-20 00:12:31 +05:30
2009-01-16 23:06:37 +05:30
!***********************************************************
! assignment of coordinates for subnodes in each cp element
!
! allocate globals
! _subNodeCoord
!***********************************************************
2009-06-15 18:41:21 +05:30
subroutine mesh_build_subNodeCoords ( )
2009-01-16 23:06:37 +05:30
use prec , only : pInt , pReal
implicit none
integer ( pInt ) e , t , n , p
2009-10-12 21:31:49 +05:30
2011-05-24 21:27:59 +05:30
if ( . not . allocated ( mesh_subNodeCoord ) ) then
allocate ( mesh_subNodeCoord ( 3 , mesh_maxNnodes + mesh_maxNsubNodes , mesh_NcpElems ) )
endif
mesh_subNodeCoord = 0.0_pReal
2009-01-16 23:06:37 +05:30
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems ! loop over cpElems
2009-01-16 23:06:37 +05:30
t = mesh_element ( 2 , e ) ! get elemType
2012-02-16 00:28:38 +05:30
do n = 1_pInt , FE_Nnodes ( t )
mesh_subNodeCoord ( 1 : 3 , n , e ) = mesh_node ( 1 : 3 , mesh_FEasCP ( 'node' , mesh_element ( 4_pInt + n , e ) ) ) ! loop over nodes of this element type
2009-01-16 23:06:37 +05:30
enddo
2012-02-16 00:28:38 +05:30
do n = 1_pInt , FE_NsubNodes ( t ) ! now for the true subnodes
do p = 1_pInt , FE_Nips ( t ) ! loop through possible parent nodes
if ( FE_subNodeParent ( p , n , t ) > 0_pInt ) & ! valid parent node
2011-05-24 21:27:59 +05:30
mesh_subNodeCoord ( 1 : 3 , FE_Nnodes ( t ) + n , e ) = mesh_subNodeCoord ( 1 : 3 , FE_Nnodes ( t ) + n , e ) &
2012-02-16 00:28:38 +05:30
+ mesh_node ( 1 : 3 , mesh_FEasCP ( 'node' , mesh_element ( 4_pInt + FE_subNodeParent ( p , n , t ) , e ) ) ) ! add up parents
2009-04-02 18:30:51 +05:30
enddo
2011-05-24 21:27:59 +05:30
mesh_subNodeCoord ( 1 : 3 , n + FE_Nnodes ( t ) , e ) = mesh_subNodeCoord ( 1 : 3 , n + FE_Nnodes ( t ) , e ) &
2012-02-16 00:28:38 +05:30
/ real ( count ( FE_subNodeParent ( : , n , t ) > 0_pInt ) , pReal )
2009-01-20 00:12:31 +05:30
enddo
2009-01-16 23:06:37 +05:30
enddo
2009-06-15 18:41:21 +05:30
endsubroutine
2009-01-16 23:06:37 +05:30
!***********************************************************
2011-08-10 22:07:17 +05:30
! calculation of IP coordinates
2009-01-16 23:06:37 +05:30
!
! allocate globals
2011-08-10 22:07:17 +05:30
! _ipCenterOfGravity
2009-01-16 23:06:37 +05:30
!***********************************************************
2011-08-10 22:07:17 +05:30
subroutine mesh_build_ipCoordinates ( )
2009-01-16 23:06:37 +05:30
2011-02-25 18:19:18 +05:30
use prec , only : pInt , tol_gravityNodePos
2009-01-16 23:06:37 +05:30
implicit none
2009-04-03 12:33:25 +05:30
integer ( pInt ) e , f , t , i , j , k , n
2012-02-16 00:28:38 +05:30
logical , dimension ( mesh_maxNnodes + mesh_maxNsubNodes ) :: gravityNode ! flagList to find subnodes determining center of grav
2009-01-16 23:06:37 +05:30
real ( pReal ) , dimension ( 3 , mesh_maxNnodes + mesh_maxNsubNodes ) :: gravityNodePos ! coordinates of subnodes determining center of grav
real ( pReal ) , dimension ( 3 ) :: centerOfGravity
2011-05-24 21:27:59 +05:30
2011-08-10 22:07:17 +05:30
if ( . not . allocated ( mesh_ipCenterOfGravity ) ) then
2011-05-24 21:27:59 +05:30
allocate ( mesh_ipCenterOfGravity ( 3 , mesh_maxNips , mesh_NcpElems ) )
endif
2009-08-11 22:01:57 +05:30
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems ! loop over cpElems
2009-10-12 21:31:49 +05:30
t = mesh_element ( 2 , e ) ! get elemType
2012-02-16 00:28:38 +05:30
do i = 1_pInt , FE_Nips ( t ) ! loop over IPs of elem
2009-12-22 18:53:20 +05:30
gravityNode = . false . ! reset flagList
2009-10-12 21:31:49 +05:30
gravityNodePos = 0.0_pReal ! reset coordinates
2012-02-16 00:28:38 +05:30
do f = 1_pInt , FE_NipNeighbors ( t ) ! loop over interfaces of IP
do n = 1_pInt , FE_NipFaceNodes ! loop over nodes on interface
2009-12-22 18:53:20 +05:30
gravityNode ( FE_subNodeOnIPFace ( n , f , i , t ) ) = . true .
2009-01-16 23:06:37 +05:30
gravityNodePos ( : , FE_subNodeOnIPFace ( n , f , i , t ) ) = mesh_subNodeCoord ( : , FE_subNodeOnIPFace ( n , f , i , t ) , e )
2009-06-15 18:41:21 +05:30
enddo
enddo
2009-12-22 18:53:20 +05:30
2012-02-16 00:28:38 +05:30
do j = 1_pInt , mesh_maxNnodes + mesh_maxNsubNodes - 1_pInt ! walk through entire flagList except last
2009-12-22 18:53:20 +05:30
if ( gravityNode ( j ) ) then ! valid node index
2012-02-16 00:28:38 +05:30
do k = j + 1_pInt , mesh_maxNnodes + mesh_maxNsubNodes ! walk through remainder of list
2011-02-25 18:19:18 +05:30
if ( gravityNode ( k ) . and . all ( abs ( gravityNodePos ( : , j ) - gravityNodePos ( : , k ) ) < tol_gravityNodePos ) ) then ! found duplicate
2009-12-22 18:53:20 +05:30
gravityNode ( j ) = . false . ! delete first instance
2009-04-02 18:30:51 +05:30
gravityNodePos ( : , j ) = 0.0_pReal
2009-10-12 21:31:49 +05:30
exit ! continue with next suspect
2009-06-15 18:41:21 +05:30
endif
enddo
endif
enddo
2012-02-10 17:26:05 +05:30
centerOfGravity = sum ( gravityNodePos , 2 ) / real ( count ( gravityNode ) , pReal )
2009-08-11 22:01:57 +05:30
mesh_ipCenterOfGravity ( : , i , e ) = centerOfGravity
2009-06-15 18:41:21 +05:30
enddo
enddo
2011-08-10 22:07:17 +05:30
endsubroutine
!***********************************************************
! calculation of IP volume
!
! allocate globals
! _ipVolume
!***********************************************************
subroutine mesh_build_ipVolumes ( )
use prec , only : pInt , tol_gravityNodePos
use math , only : math_volTetrahedron
implicit none
2009-01-16 23:06:37 +05:30
2011-08-10 22:07:17 +05:30
integer ( pInt ) e , f , t , i , j , n
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: Ntriangles = FE_NipFaceNodes - 2_pInt ! each interface is made up of this many triangles
2011-08-10 22:07:17 +05:30
real ( pReal ) , dimension ( 3 , FE_NipFaceNodes ) :: nPos ! coordinates of nodes on IP face
real ( pReal ) , dimension ( Ntriangles , FE_NipFaceNodes ) :: volume ! volumes of possible tetrahedra
if ( . not . allocated ( mesh_ipVolume ) ) then
allocate ( mesh_ipVolume ( mesh_maxNips , mesh_NcpElems ) )
endif
mesh_ipVolume = 0.0_pReal
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems ! loop over cpElems
2012-02-10 17:26:05 +05:30
t = mesh_element ( 2_pInt , e ) ! get elemType
2012-02-16 00:28:38 +05:30
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 ) &
2011-08-10 22:07:17 +05:30
nPos ( : , n ) = mesh_subNodeCoord ( : , FE_subNodeOnIPFace ( n , f , i , t ) , e )
2012-02-16 00:28:38 +05:30
forall ( n = 1_pInt : FE_NipFaceNodes , j = 1_pInt : Ntriangles ) & ! start at each interface node and build valid triangles to cover interface
2011-08-10 22:07:17 +05:30
volume ( j , n ) = math_volTetrahedron ( nPos ( : , n ) , & ! calc volume of respective tetrahedron to CoG
2012-02-10 17:26:05 +05:30
nPos ( : , 1_pInt + mod ( n - 1_pInt + j , FE_NipFaceNodes ) ) , & ! start at offset j
2012-02-16 00:28:38 +05:30
nPos ( : , 1_pInt + mod ( n - 1_pInt + j + 1_pInt , FE_NipFaceNodes ) ) , & ! and take j's neighbor
2011-08-10 22:07:17 +05:30
mesh_ipCenterOfGravity ( : , i , e ) )
mesh_ipVolume ( i , e ) = mesh_ipVolume ( i , e ) + sum ( volume ) ! add contribution from this interface
enddo
mesh_ipVolume ( i , e ) = mesh_ipVolume ( i , e ) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them
enddo
enddo
2009-06-15 18:41:21 +05:30
endsubroutine
2009-01-16 23:06:37 +05:30
!***********************************************************
2009-03-04 17:18:54 +05:30
! calculation of IP interface areas
2009-01-16 23:06:37 +05:30
!
! allocate globals
2009-03-04 17:18:54 +05:30
! _ipArea, _ipAreaNormal
2009-01-16 23:06:37 +05:30
!***********************************************************
2009-06-15 18:41:21 +05:30
subroutine mesh_build_ipAreas ( )
2009-01-16 23:06:37 +05:30
use prec , only : pInt , pReal
use math
implicit none
2009-03-05 19:48:50 +05:30
integer ( pInt ) e , f , t , i , j , n
2012-02-16 00:28:38 +05:30
integer ( pInt ) , parameter :: Ntriangles = FE_NipFaceNodes - 2_pInt ! each interface is made up of this many triangles
2009-10-12 21:31:49 +05:30
real ( pReal ) , dimension ( 3 , FE_NipFaceNodes ) :: nPos ! coordinates of nodes on IP face
2009-01-20 00:12:31 +05:30
real ( pReal ) , dimension ( 3 , Ntriangles , FE_NipFaceNodes ) :: normal
real ( pReal ) , dimension ( Ntriangles , FE_NipFaceNodes ) :: area
2009-01-16 23:06:37 +05:30
2009-01-20 00:12:31 +05:30
allocate ( mesh_ipArea ( mesh_maxNipNeighbors , mesh_maxNips , mesh_NcpElems ) ) ; mesh_ipArea = 0.0_pReal
2012-02-10 17:26:05 +05:30
allocate ( mesh_ipAreaNormal ( 3_pInt , mesh_maxNipNeighbors , mesh_maxNips , mesh_NcpElems ) ) ; mesh_ipAreaNormal = 0.0_pReal
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems ! loop over cpElems
2009-10-12 21:31:49 +05:30
t = mesh_element ( 2 , e ) ! get elemType
2012-02-16 00:28:38 +05:30
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 )
forall ( n = 1_pInt : FE_NipFaceNodes , j = 1_pInt : Ntriangles ) ! start at each interface node and build valid triangles to cover interface
2012-02-10 17:26:05 +05:30
normal ( : , j , n ) = math_vectorproduct ( nPos ( : , 1_pInt + mod ( n + j - 1_pInt , FE_NipFaceNodes ) ) - nPos ( : , n ) , & ! calc their normal vectors
nPos ( : , 1_pInt + mod ( n + j - 0_pInt , FE_NipFaceNodes ) ) - nPos ( : , n ) )
2011-02-25 14:55:53 +05:30
area ( j , n ) = sqrt ( sum ( normal ( : , j , n ) * normal ( : , j , n ) ) ) ! and area
2009-01-20 00:12:31 +05:30
end forall
2012-02-16 00:28:38 +05:30
forall ( n = 1_pInt : FE_NipFaceNodes , j = 1_pInt : Ntriangles , area ( j , n ) > 0.0_pReal ) &
2011-10-18 14:50:29 +05:30
normal ( 1 : 3 , j , n ) = normal ( 1 : 3 , j , n ) / area ( j , n ) ! make myUnit normal
2009-08-11 22:01:57 +05:30
2009-10-12 21:31:49 +05:30
mesh_ipArea ( f , i , e ) = sum ( area ) / ( FE_NipFaceNodes * 2.0_pReal ) ! area of parallelograms instead of triangles
2012-02-16 00:28:38 +05:30
mesh_ipAreaNormal ( : , f , i , e ) = sum ( sum ( normal , 3 ) , 2_pInt ) / & ! average of all valid normals
2012-02-10 17:26:05 +05:30
real ( count ( area > 0.0_pReal ) , pReal )
2009-01-16 23:06:37 +05:30
enddo
enddo
enddo
2007-04-10 16:52:53 +05:30
2009-06-15 18:41:21 +05:30
endsubroutine
2009-01-16 23:06:37 +05:30
2011-02-16 21:53:08 +05:30
!***********************************************************
! 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
2011-07-31 21:12:59 +05:30
logical , dimension ( mesh_Nnodes ) :: unpaired
2011-02-16 21:53:08 +05:30
allocate ( mesh_nodeTwins ( 3 , mesh_Nnodes ) )
mesh_nodeTwins = 0_pInt
tolerance = 0.001_pReal * minval ( mesh_ipVolume ) ** 0.333_pReal
2012-02-16 00:28:38 +05:30
do dir = 1_pInt , 3_pInt ! check periodicity in directions of x,y,z
2011-02-16 21:53:08 +05:30
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
2011-05-24 21:27:59 +05:30
minCoord = minval ( mesh_node0 ( dir , : ) )
maxCoord = maxval ( mesh_node0 ( dir , : ) )
2012-02-16 00:28:38 +05:30
do node = 1_pInt , mesh_Nnodes ! loop through all nodes and find surface nodes
2011-05-24 21:27:59 +05:30
if ( abs ( mesh_node0 ( dir , node ) - minCoord ) < = tolerance ) then
2011-02-16 21:53:08 +05:30
minimumNodes ( 1 ) = minimumNodes ( 1 ) + 1_pInt
2012-02-16 00:28:38 +05:30
minimumNodes ( minimumNodes ( 1 ) + 1_pInt ) = node
2011-05-24 21:27:59 +05:30
elseif ( abs ( mesh_node0 ( dir , node ) - maxCoord ) < = tolerance ) then
2011-02-16 21:53:08 +05:30
maximumNodes ( 1 ) = maximumNodes ( 1 ) + 1_pInt
2012-02-16 00:28:38 +05:30
maximumNodes ( maximumNodes ( 1 ) + 1_pInt ) = node
2011-02-16 21:53:08 +05:30
endif
enddo
!*** find the corresponding node on the other side with the same position in this dimension
2011-07-31 21:12:59 +05:30
unpaired = . true .
2012-02-16 00:28:38 +05:30
do n1 = 1_pInt , minimumNodes ( 1 )
minimumNode = minimumNodes ( n1 + 1_pInt )
2011-07-31 21:12:59 +05:30
if ( unpaired ( minimumNode ) ) then
2012-02-16 00:28:38 +05:30
do n2 = 1_pInt , maximumNodes ( 1 )
maximumNode = maximumNodes ( n2 + 1_pInt )
2011-07-31 21:12:59 +05:30
distance = abs ( mesh_node0 ( : , minimumNode ) - mesh_node0 ( : , maximumNode ) )
if ( sum ( distance ) - distance ( dir ) < = tolerance ) then ! minimum possible distance (within tolerance)
mesh_nodeTwins ( dir , minimumNode ) = maximumNode
mesh_nodeTwins ( dir , maximumNode ) = minimumNode
unpaired ( maximumNode ) = . false . ! remember this node, we don't have to look for his partner again
exit
endif
enddo
endif
2011-02-16 21:53:08 +05:30
enddo
endif
enddo
endsubroutine
2007-10-24 14:30:42 +05:30
!***********************************************************
! write statistics regarding input file parsing
! to the output file
!
!***********************************************************
2011-02-16 21:53:08 +05:30
subroutine mesh_tell_statistics ( )
2009-01-20 00:12:31 +05:30
2011-02-16 21:53:08 +05:30
use prec , only : pInt
use math , only : math_range
use IO , only : IO_error
2011-03-21 16:01:17 +05:30
use debug , only : debug_verbosity , &
debug_e , &
debug_i , &
debug_selectiveDebugger
2009-01-20 00:12:31 +05:30
2011-02-16 21:53:08 +05:30
implicit none
2009-01-20 00:12:31 +05:30
2011-02-16 21:53:08 +05:30
integer ( pInt ) , dimension ( : , : ) , allocatable :: mesh_HomogMicro
character ( len = 64 ) fmt
2009-01-20 00:12:31 +05:30
2011-02-16 21:53:08 +05:30
integer ( pInt ) i , e , n , f , t
2008-06-17 02:19:48 +05:30
2012-02-13 23:11:27 +05:30
if ( mesh_maxValStateVar ( 1 ) < 1_pInt ) call IO_error ( error_ID = 170_pInt ) ! no homogenization specified
if ( mesh_maxValStateVar ( 2 ) < 1_pInt ) call IO_error ( error_ID = 180_pInt ) ! no microstructure specified
2011-02-16 21:53:08 +05:30
allocate ( mesh_HomogMicro ( mesh_maxValStateVar ( 1 ) , mesh_maxValStateVar ( 2 ) ) ) ; mesh_HomogMicro = 0_pInt
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems
2012-02-13 23:11:27 +05:30
if ( mesh_element ( 3 , e ) < 1_pInt ) call IO_error ( error_ID = 170_pInt , e = e ) ! no homogenization specified
if ( mesh_element ( 4 , e ) < 1_pInt ) call IO_error ( error_ID = 180_pInt , e = e ) ! no microstructure specified
2011-02-16 21:53:08 +05:30
mesh_HomogMicro ( mesh_element ( 3 , e ) , mesh_element ( 4 , e ) ) = &
2012-02-16 00:28:38 +05:30
mesh_HomogMicro ( mesh_element ( 3 , e ) , mesh_element ( 4 , e ) ) + 1_pInt ! count combinations of homogenization and microstructure
2011-02-16 21:53:08 +05:30
enddo
2010-09-07 14:36:02 +05:30
2012-02-16 00:28:38 +05:30
if ( debug_verbosity > 0_pInt ) then
2011-03-21 16:01:17 +05:30
!$OMP CRITICAL (write2out)
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , * ) 'Input Parser: STATISTICS'
2011-03-21 16:01:17 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , * ) mesh_Nelems , ' : total number of elements in mesh'
write ( 6 , * ) mesh_NcpElems , ' : total number of CP elements in mesh'
write ( 6 , * ) mesh_Nnodes , ' : total number of nodes in mesh'
write ( 6 , * ) mesh_maxNnodes , ' : max number of nodes in any CP element'
write ( 6 , * ) mesh_maxNips , ' : max number of IPs in any CP element'
write ( 6 , * ) mesh_maxNipNeighbors , ' : max number of IP neighbors in any CP element'
write ( 6 , * ) mesh_maxNsubNodes , ' : max number of (additional) subnodes in any CP element'
write ( 6 , * ) mesh_maxNsharedElems , ' : max number of CP elements sharing a node'
2011-03-21 16:01:17 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , * ) 'Input Parser: HOMOGENIZATION/MICROSTRUCTURE'
2011-03-21 16:01:17 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , * ) mesh_maxValStateVar ( 1 ) , ' : maximum homogenization index'
write ( 6 , * ) mesh_maxValStateVar ( 2 ) , ' : maximum microstructure index'
2011-03-21 16:01:17 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( fmt , '(a,i32.32,a)' ) '(9x,a2,1x,' , mesh_maxValStateVar ( 2 ) , '(i8))'
write ( 6 , fmt ) '+-' , math_range ( mesh_maxValStateVar ( 2 ) )
write ( fmt , '(a,i32.32,a)' ) '(i8,1x,a2,1x,' , mesh_maxValStateVar ( 2 ) , '(i8))'
2012-02-16 00:28:38 +05:30
do i = 1_pInt , mesh_maxValStateVar ( 1 ) ! loop over all (possibly assigned) homogenizations
2012-02-01 00:48:55 +05:30
write ( 6 , fmt ) i , '| ' , mesh_HomogMicro ( i , : ) ! loop over all (possibly assigned) microstructures
2011-03-21 16:01:17 +05:30
enddo
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , * ) 'Input Parser: ADDITIONAL MPIE OPTIONS'
2011-03-21 16:01:17 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , * ) 'periodic surface : ' , mesh_periodicSurface
2011-03-21 16:01:17 +05:30
write ( 6 , * )
call flush ( 6 )
!$OMP END CRITICAL (write2out)
endif
if ( debug_verbosity > 1 ) then
2010-09-07 14:36:02 +05:30
!$OMP CRITICAL (write2out)
2011-02-16 21:53:08 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , * ) 'Input Parser: SUBNODE COORDINATES'
2011-02-16 21:53:08 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , '(a8,1x,a5,1x,a15,1x,a15,1x,a20,3(1x,a12))' ) 'elem' , 'IP' , 'IP neighbor' , 'IPFaceNodes' , 'subNodeOnIPFace' , 'x' , 'y' , 'z'
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems ! loop over cpElems
2011-03-21 16:01:17 +05:30
if ( debug_selectiveDebugger . and . debug_e / = e ) cycle
2011-02-16 21:53:08 +05:30
t = mesh_element ( 2 , e ) ! get elemType
2012-02-16 00:28:38 +05:30
do i = 1_pInt , FE_Nips ( t ) ! loop over IPs of elem
2011-03-21 16:01:17 +05:30
if ( debug_selectiveDebugger . and . debug_i / = i ) cycle
2012-02-16 00:28:38 +05:30
do f = 1_pInt , FE_NipNeighbors ( t ) ! loop over interfaces of IP
do n = 1_pInt , FE_NipFaceNodes ! loop over nodes on interface
2012-02-01 00:48:55 +05:30
write ( 6 , '(i8,1x,i5,1x,i15,1x,i15,1x,i20,3(1x,f12.8))' ) e , i , f , n , FE_subNodeOnIPFace ( n , f , i , t ) , &
2011-02-16 21:53:08 +05:30
mesh_subNodeCoord ( 1 , FE_subNodeOnIPFace ( n , f , i , t ) , e ) , &
mesh_subNodeCoord ( 2 , FE_subNodeOnIPFace ( n , f , i , t ) , e ) , &
mesh_subNodeCoord ( 3 , FE_subNodeOnIPFace ( n , f , i , t ) , e )
enddo
enddo
enddo
enddo
write ( 6 , * )
write ( 6 , * ) 'Input Parser: IP COORDINATES'
2012-02-01 00:48:55 +05:30
write ( 6 , '(a8,1x,a5,3(1x,a12))' ) 'elem' , 'IP' , 'x' , 'y' , 'z'
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems
2011-03-21 16:01:17 +05:30
if ( debug_selectiveDebugger . and . debug_e / = e ) cycle
2012-02-16 00:28:38 +05:30
do i = 1_pInt , FE_Nips ( mesh_element ( 2 , e ) )
2011-03-21 16:01:17 +05:30
if ( debug_selectiveDebugger . and . debug_i / = i ) cycle
2012-02-01 00:48:55 +05:30
write ( 6 , '(i8,1x,i5,3(1x,f12.8))' ) e , i , mesh_ipCenterOfGravity ( : , i , e )
2011-02-16 21:53:08 +05:30
enddo
enddo
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , * ) 'Input Parser: ELEMENT VOLUME'
2010-09-07 14:36:02 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , '(a13,1x,e15.8)' ) 'total volume' , sum ( mesh_ipVolume )
2010-09-07 14:36:02 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , '(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)' ) 'elem' , 'IP' , 'volume' , 'face' , 'area' , '-- normal --'
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems
2011-03-21 16:01:17 +05:30
if ( debug_selectiveDebugger . and . debug_e / = e ) cycle
2012-02-16 00:28:38 +05:30
do i = 1_pInt , FE_Nips ( mesh_element ( 2 , e ) )
2011-03-21 16:01:17 +05:30
if ( debug_selectiveDebugger . and . debug_i / = i ) cycle
2012-02-01 00:48:55 +05:30
write ( 6 , '(i8,1x,i5,1x,e15.8)' ) e , i , mesh_IPvolume ( i , e )
2012-02-16 00:28:38 +05:30
do f = 1_pInt , FE_NipNeighbors ( mesh_element ( 2 , e ) )
2012-02-01 00:48:55 +05:30
write ( 6 , '(i33,1x,e15.8,1x,3(f6.3,1x))' ) f , mesh_ipArea ( f , i , e ) , mesh_ipAreaNormal ( : , f , i , e )
2010-09-07 14:36:02 +05:30
enddo
enddo
enddo
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , * ) 'Input Parser: NODE TWINS'
2011-02-16 21:53:08 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , '(a6,3(3x,a6))' ) ' node' , 'twin_x' , 'twin_y' , 'twin_z'
2012-02-16 00:28:38 +05:30
do n = 1_pInt , mesh_Nnodes ! loop over cpNodes
2011-03-21 16:01:17 +05:30
if ( debug_e < = mesh_NcpElems ) then
if ( any ( mesh_element ( 5 : , debug_e ) == n ) ) then
2012-02-01 00:48:55 +05:30
write ( 6 , '(i6,3(3x,i6))' ) n , mesh_nodeTwins ( 1 : 3 , n )
2011-03-21 16:01:17 +05:30
endif
endif
2011-02-16 21:53:08 +05:30
enddo
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , * ) 'Input Parser: IP NEIGHBORHOOD'
2011-02-16 21:53:08 +05:30
write ( 6 , * )
2012-02-01 00:48:55 +05:30
write ( 6 , '(a8,1x,a10,1x,a10,1x,a3,1x,a13,1x,a13)' ) 'elem' , 'IP' , 'neighbor' , '' , 'elemNeighbor' , 'ipNeighbor'
2012-02-16 00:28:38 +05:30
do e = 1_pInt , mesh_NcpElems ! loop over cpElems
2011-03-21 16:01:17 +05:30
if ( debug_selectiveDebugger . and . debug_e / = e ) cycle
2011-02-16 21:53:08 +05:30
t = mesh_element ( 2 , e ) ! get elemType
2012-02-16 00:28:38 +05:30
do i = 1_pInt , FE_Nips ( t ) ! loop over IPs of elem
2011-03-21 16:01:17 +05:30
if ( debug_selectiveDebugger . and . debug_i / = i ) cycle
2012-02-16 00:28:38 +05:30
do n = 1_pInt , FE_NipNeighbors ( t ) ! loop over neighbors of IP
2012-02-01 00:48:55 +05:30
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 )
2011-02-16 21:53:08 +05:30
enddo
enddo
enddo
2010-09-30 13:01:53 +05:30
!$OMP END CRITICAL (write2out)
2011-02-16 21:53:08 +05:30
endif
deallocate ( mesh_HomogMicro )
endsubroutine
2012-02-13 23:11:27 +05:30
END MODULE mesh