2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
2013-04-18 22:10:49 +05:30
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
2018-01-10 21:43:25 +05:30
!> @brief Sets up the mesh for the solvers MSC.Marc, Abaqus and the spectral solver
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
2014-04-29 23:20:59 +05:30
module mesh
2012-08-27 13:34:47 +05:30
use , intrinsic :: iso_c_binding
2012-03-09 01:55:28 +05:30
use prec , only : pReal , pInt
2019-01-28 18:29:54 +05:30
use mesh_base
2012-08-06 18:13:05 +05:30
2007-03-21 21:48:33 +05:30
implicit none
2012-03-09 01:55:28 +05:30
private
2013-01-16 16:10:53 +05:30
integer ( pInt ) , public , protected :: &
2015-03-25 21:36:19 +05:30
mesh_NcpElems , & !< total number of CP elements in local mesh
2018-09-24 00:23:35 +05:30
mesh_elemType , & !< Element type of the mesh (only support homogeneous meshes)
2012-06-15 21:40:21 +05:30
mesh_Nnodes , & !< total number of nodes in mesh
2013-04-15 13:43:20 +05:30
mesh_Ncellnodes , & !< total number of cell nodes in mesh (including duplicates)
2013-04-18 22:10:49 +05:30
mesh_Ncells , & !< total number of cells in mesh
2018-09-23 21:07:57 +05:30
mesh_NcellnodesPerElem , & !< number of cell nodes per element
2012-06-15 21:40:21 +05:30
mesh_maxNipNeighbors , & !< max number of IP neighbors in any CP element
2018-09-23 21:07:57 +05:30
mesh_maxNsharedElems !< max number of CP elements sharing a node
!!!! BEGIN DEPRECATED !!!!!
integer ( pInt ) , public , protected :: &
mesh_maxNips , & !< max number of IPs in any CP element
2018-09-23 20:35:01 +05:30
mesh_maxNcellnodes !< max number of cell nodes in any CP element
2018-09-23 21:07:57 +05:30
!!!! BEGIN DEPRECATED !!!!!
2012-03-09 01:55:28 +05:30
2018-09-23 22:20:54 +05:30
integer ( pInt ) , dimension ( : ) , allocatable , public , protected :: &
2018-10-04 09:33:48 +05:30
mesh_homogenizationAt , & !< homogenization ID of each element
mesh_microstructureAt !< microstructure ID of each element
2018-09-23 22:20:54 +05:30
2013-01-16 16:10:53 +05:30
integer ( pInt ) , dimension ( : , : ) , allocatable , public , protected :: &
2018-09-23 23:28:43 +05:30
mesh_element , & !DEPRECATED
2012-06-15 21:40:21 +05:30
mesh_sharedElem , & !< entryCount and list of elements containing node
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)
2018-01-10 21:43:25 +05:30
2013-01-16 16:10:53 +05:30
integer ( pInt ) , dimension ( : , : , : , : ) , allocatable , public , protected :: &
2012-10-24 19:33:02 +05:30
mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me]
2012-03-09 01:55:28 +05:30
2013-05-07 18:36:29 +05:30
real ( pReal ) , public , protected :: &
mesh_unitlength !< physical length of one unit in mesh
2012-03-09 01:55:28 +05:30
real ( pReal ) , dimension ( : , : ) , allocatable , public :: &
2013-04-22 00:18:59 +05:30
mesh_node , & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
mesh_cellnode !< cell node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
2018-01-10 21:43:25 +05:30
2013-01-16 16:10:53 +05:30
real ( pReal ) , dimension ( : , : ) , allocatable , public , protected :: &
mesh_ipVolume , & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!)
real ( pReal ) , dimension ( : , : , : ) , allocatable , public , protected :: &
2013-01-31 21:58:08 +05:30
mesh_ipArea !< area of interface to neighboring IP (initially!)
2018-01-10 21:43:25 +05:30
2013-01-16 16:10:53 +05:30
real ( pReal ) , dimension ( : , : , : ) , allocatable , public :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
2018-01-10 21:43:25 +05:30
real ( pReal ) , dimension ( : , : , : , : ) , allocatable , public , protected :: &
2012-06-15 21:40:21 +05:30
mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!)
2018-01-10 21:43:25 +05:30
2013-01-16 16:10:53 +05:30
logical , dimension ( 3 ) , public , protected :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes)
2013-02-28 02:11:14 +05:30
2018-10-10 03:22:54 +05:30
integer ( pInt ) , private :: &
mesh_maxNelemInSet , &
mesh_Nmaterials
2012-03-09 01:55:28 +05:30
integer ( pInt ) , dimension ( 2 ) , private :: &
mesh_maxValStateVar = 0_pInt
2018-01-10 21:43:25 +05:30
2018-09-23 18:57:51 +05:30
integer ( pInt ) , dimension ( : , : ) , allocatable , private :: &
2013-12-28 01:33:28 +05:30
mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( : , : , : ) , allocatable , private :: &
2013-04-22 00:18:59 +05:30
mesh_cell !< cell connectivity for each element,ip/cell
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( : , : , : ) , allocatable , private :: &
FE_nodesAtIP , & !< map IP index to node indices in a specific type of element
FE_ipNeighbor , & !< +x,-x,+y,-y,+z,-z list of intra-element IPs and(negative) neighbor faces per own IP in a specific type of element
FE_cell , & !< list of intra-element cell node IDs that constitute the cells in a specific type of element geometry
FE_cellface !< list of intra-cell cell node IDs that constitute the cell faces of a specific type of cell
2018-01-10 21:43:25 +05:30
2013-04-15 13:43:20 +05:30
real ( pReal ) , dimension ( : , : , : ) , allocatable , private :: &
FE_cellnodeParentnodeWeights !< list of node weights for the generation of cell nodes
2018-01-10 21:43:25 +05:30
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( : , : , : , : ) , allocatable , private :: &
FE_subNodeOnIPFace
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS)
2007-03-22 20:18:58 +05:30
! Hence, I suggest to prefix with "FE_"
2012-03-09 01:55:28 +05:30
integer ( pInt ) , parameter , public :: &
2013-04-15 13:43:20 +05:30
FE_Nelemtypes = 13_pInt , &
FE_Ngeomtypes = 10_pInt , &
FE_Ncelltypes = 4_pInt , &
FE_maxNnodes = 20_pInt , &
FE_maxNips = 27_pInt , &
FE_maxNipNeighbors = 6_pInt , &
2012-06-15 21:40:21 +05:30
FE_maxmaxNnodesAtIP = 8_pInt , & !< max number of (equivalent) nodes attached to an IP
2013-04-15 13:43:20 +05:30
FE_maxNmatchingNodesPerFace = 4_pInt , &
FE_maxNfaces = 6_pInt , &
FE_maxNcellnodes = 64_pInt , &
FE_maxNcellnodesPerCell = 8_pInt , &
FE_maxNcellfaces = 6_pInt , &
FE_maxNcellnodesPerCellface = 4_pInt
2018-01-10 21:43:25 +05:30
2012-11-16 04:15:20 +05:30
integer ( pInt ) , dimension ( FE_Nelemtypes ) , parameter , public :: FE_geomtype = & !< geometry type of particular element type
int ( [ &
1 , & ! element 6 (2D 3node 1ip)
2 , & ! element 125 (2D 6node 3ip)
3 , & ! element 11 (2D 4node 4ip)
4 , & ! element 27 (2D 8node 9ip)
2013-04-10 15:08:40 +05:30
3 , & ! element 54 (2D 8node 4ip)
2012-11-16 04:15:20 +05:30
5 , & ! element 134 (3D 4node 1ip)
6 , & ! element 157 (3D 5node 4ip)
6 , & ! element 127 (3D 10node 4ip)
7 , & ! element 136 (3D 6node 6ip)
8 , & ! element 117 (3D 8node 1ip)
9 , & ! element 7 (3D 8node 8ip)
9 , & ! element 57 (3D 20node 8ip)
10 & ! element 21 (3D 20node 27ip)
2012-02-16 00:28:38 +05:30
] , pInt )
2012-11-16 04:15:20 +05:30
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( FE_Ngeomtypes ) , parameter , public :: FE_celltype = & !< cell type that is used by each geometry type
2012-11-16 04:15:20 +05:30
int ( [ &
2013-04-15 13:43:20 +05:30
1 , & ! element 6 (2D 3node 1ip)
2 , & ! element 125 (2D 6node 3ip)
2 , & ! element 11 (2D 4node 4ip)
2 , & ! element 27 (2D 8node 9ip)
3 , & ! element 134 (3D 4node 1ip)
4 , & ! element 127 (3D 10node 4ip)
4 , & ! element 136 (3D 6node 6ip)
4 , & ! element 117 (3D 8node 1ip)
4 , & ! element 7 (3D 8node 8ip)
4 & ! element 21 (3D 20node 27ip)
2012-02-16 00:28:38 +05:30
] , pInt )
2012-11-16 04:15:20 +05:30
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( FE_Ngeomtypes ) , parameter , public :: FE_dimension = & !< dimension of geometry type
2013-02-05 18:57:37 +05:30
int ( [ &
2 , & ! element 6 (2D 3node 1ip)
2 , & ! element 125 (2D 6node 3ip)
2 , & ! element 11 (2D 4node 4ip)
2 , & ! element 27 (2D 8node 9ip)
3 , & ! element 134 (3D 4node 1ip)
3 , & ! element 127 (3D 10node 4ip)
3 , & ! element 136 (3D 6node 6ip)
3 , & ! element 117 (3D 8node 1ip)
3 , & ! element 7 (3D 8node 8ip)
3 & ! element 21 (3D 20node 27ip)
] , pInt )
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( FE_Nelemtypes ) , parameter , public :: FE_Nnodes = & !< number of nodes that constitute a specific type of element
2012-11-16 04:15:20 +05:30
int ( [ &
3 , & ! element 6 (2D 3node 1ip)
2013-04-15 13:43:20 +05:30
6 , & ! element 125 (2D 6node 3ip)
2012-11-16 04:15:20 +05:30
4 , & ! element 11 (2D 4node 4ip)
2013-04-15 13:43:20 +05:30
8 , & ! element 27 (2D 8node 9ip)
8 , & ! element 54 (2D 8node 4ip)
2012-11-16 04:15:20 +05:30
4 , & ! element 134 (3D 4node 1ip)
2013-04-15 13:43:20 +05:30
5 , & ! element 157 (3D 5node 4ip)
10 , & ! element 127 (3D 10node 4ip)
2012-11-16 04:15:20 +05:30
6 , & ! element 136 (3D 6node 6ip)
8 , & ! element 117 (3D 8node 1ip)
8 , & ! element 7 (3D 8node 8ip)
2013-04-15 13:43:20 +05:30
20 , & ! element 57 (3D 20node 8ip)
20 & ! element 21 (3D 20node 27ip)
2012-02-16 00:28:38 +05:30
] , pInt )
2012-11-16 04:15:20 +05:30
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( FE_Ngeomtypes ) , parameter , public :: FE_Nfaces = & !< number of faces of a specific type of element geometry
2012-11-16 04:15:20 +05:30
int ( [ &
2013-04-15 13:43:20 +05:30
3 , & ! element 6 (2D 3node 1ip)
2012-11-16 04:15:20 +05:30
3 , & ! element 125 (2D 6node 3ip)
4 , & ! element 11 (2D 4node 4ip)
2013-04-15 13:43:20 +05:30
4 , & ! element 27 (2D 8node 9ip)
4 , & ! element 134 (3D 4node 1ip)
2012-11-16 04:15:20 +05:30
4 , & ! element 127 (3D 10node 4ip)
2013-04-15 13:43:20 +05:30
5 , & ! element 136 (3D 6node 6ip)
6 , & ! element 117 (3D 8node 1ip)
6 , & ! element 7 (3D 8node 8ip)
6 & ! element 21 (3D 20node 27ip)
2012-02-16 00:28:38 +05:30
] , pInt )
2012-11-16 04:15:20 +05:30
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( FE_Ngeomtypes ) , parameter , private :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry
2012-11-16 04:15:20 +05:30
int ( [ &
3 , & ! element 6 (2D 3node 1ip)
2013-04-15 13:43:20 +05:30
3 , & ! element 125 (2D 6node 3ip)
2012-11-16 04:15:20 +05:30
4 , & ! element 11 (2D 4node 4ip)
4 , & ! element 27 (2D 8node 9ip)
4 , & ! element 134 (3D 4node 1ip)
2013-04-15 13:43:20 +05:30
4 , & ! element 127 (3D 10node 4ip)
2012-11-16 04:15:20 +05:30
6 , & ! element 136 (3D 6node 6ip)
2013-04-15 13:43:20 +05:30
8 , & ! element 117 (3D 8node 1ip)
8 , & ! element 7 (3D 8node 8ip)
8 & ! element 21 (3D 20node 27ip)
2012-11-16 04:15:20 +05:30
] , pInt )
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( FE_maxNfaces , FE_Ngeomtypes ) , parameter , private :: &
FE_NmatchingNodesPerFace = & !< number of matching nodes per face in a specific type of element geometry
2012-11-16 04:15:20 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
2 , 2 , 2 , 0 , 0 , 0 , & ! element 6 (2D 3node 1ip)
2012-11-16 04:15:20 +05:30
2 , 2 , 2 , 0 , 0 , 0 , & ! element 125 (2D 6node 3ip)
2 , 2 , 2 , 2 , 0 , 0 , & ! element 11 (2D 4node 4ip)
2 , 2 , 2 , 2 , 0 , 0 , & ! element 27 (2D 8node 9ip)
3 , 3 , 3 , 3 , 0 , 0 , & ! element 134 (3D 4node 1ip)
3 , 3 , 3 , 3 , 0 , 0 , & ! element 127 (3D 10node 4ip)
3 , 4 , 4 , 4 , 3 , 0 , & ! element 136 (3D 6node 6ip)
4 , 4 , 4 , 4 , 4 , 4 , & ! element 117 (3D 8node 1ip)
4 , 4 , 4 , 4 , 4 , 4 , & ! element 7 (3D 8node 8ip)
4 , 4 , 4 , 4 , 4 , 4 & ! element 21 (3D 20node 27ip)
2018-01-10 21:43:25 +05:30
] , pInt ) , [ FE_maxNipNeighbors , FE_Ngeomtypes ] )
2012-11-16 04:15:20 +05:30
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( FE_maxNmatchingNodesPerFace , FE_maxNfaces , FE_Ngeomtypes ) , &
parameter , private :: FE_face = & !< List of node indices on each face of a specific type of element geometry
2012-02-16 00:28:38 +05:30
reshape ( int ( [ &
2012-11-16 04:15:20 +05:30
1 , 2 , 0 , 0 , & ! element 6 (2D 3node 1ip)
2 , 3 , 0 , 0 , &
3 , 1 , 0 , 0 , &
2007-04-25 13:03:24 +05:30
0 , 0 , 0 , 0 , &
2007-10-12 19:18:29 +05:30
0 , 0 , 0 , 0 , &
2012-11-16 04:15:20 +05:30
0 , 0 , 0 , 0 , &
1 , 2 , 0 , 0 , & ! element 125 (2D 6node 3ip)
2 , 3 , 0 , 0 , &
3 , 1 , 0 , 0 , &
0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , &
1 , 2 , 0 , 0 , & ! element 11 (2D 4node 4ip)
2007-10-12 19:18:29 +05:30
2 , 3 , 0 , 0 , &
3 , 4 , 0 , 0 , &
4 , 1 , 0 , 0 , &
0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , &
2012-11-16 04:15:20 +05:30
1 , 2 , 0 , 0 , & ! element 27 (2D 8node 9ip)
2009-04-02 18:30:51 +05:30
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 , &
2012-11-16 04:15:20 +05:30
1 , 2 , 3 , 0 , & ! element 134 (3D 4node 1ip)
2008-06-17 14:41:54 +05:30
1 , 4 , 2 , 0 , &
2 , 3 , 4 , 0 , &
1 , 3 , 4 , 0 , &
0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , &
2012-11-16 04:15:20 +05:30
1 , 2 , 3 , 0 , & ! element 127 (3D 10node 4ip)
1 , 4 , 2 , 0 , &
2 , 4 , 3 , 0 , &
1 , 3 , 4 , 0 , &
0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , &
1 , 2 , 3 , 0 , & ! element 136 (3D 6node 6ip)
2008-06-17 14:41:54 +05:30
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 , &
2012-11-16 04:15:20 +05:30
1 , 2 , 3 , 4 , & ! element 117 (3D 8node 1ip)
2009-04-02 18:30:51 +05:30
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 , &
2012-11-16 04:15:20 +05:30
1 , 2 , 3 , 4 , & ! element 7 (3D 8node 8ip)
2010-05-10 20:24:59 +05:30
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 , &
2012-11-16 04:15:20 +05:30
1 , 2 , 3 , 4 , & ! element 21 (3D 20node 27ip)
2010-08-17 04:23:24 +05:30
2 , 1 , 5 , 6 , &
3 , 2 , 6 , 7 , &
4 , 3 , 7 , 8 , &
4 , 1 , 5 , 8 , &
2012-11-16 04:15:20 +05:30
8 , 7 , 6 , 5 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_maxNmatchingNodesPerFace , FE_maxNfaces , FE_Ngeomtypes ] )
integer ( pInt ) , dimension ( FE_Ngeomtypes ) , parameter , private :: FE_Ncellnodes = & !< number of cell nodes in a specific geometry type
int ( [ &
3 , & ! element 6 (2D 3node 1ip)
7 , & ! element 125 (2D 6node 3ip)
9 , & ! element 11 (2D 4node 4ip)
2013-04-18 23:26:14 +05:30
16 , & ! element 27 (2D 8node 9ip)
2013-04-15 13:43:20 +05:30
4 , & ! element 134 (3D 4node 1ip)
15 , & ! element 127 (3D 10node 4ip)
2013-04-18 23:26:14 +05:30
21 , & ! element 136 (3D 6node 6ip)
2013-04-15 13:43:20 +05:30
8 , & ! element 117 (3D 8node 1ip)
27 , & ! element 7 (3D 8node 8ip)
64 & ! element 21 (3D 20node 27ip)
] , pInt )
integer ( pInt ) , dimension ( FE_Ncelltypes ) , parameter , private :: FE_NcellnodesPerCell = & !< number of cell nodes in a specific cell type
int ( [ &
3 , & ! (2D 3node)
4 , & ! (2D 4node)
4 , & ! (3D 4node)
8 & ! (3D 8node)
] , pInt )
integer ( pInt ) , dimension ( FE_Ncelltypes ) , parameter , private :: FE_NcellnodesPerCellface = & !< number of cell nodes per cell face in a specific cell type
int ( [ &
2 , & ! (2D 3node)
2 , & ! (2D 4node)
3 , & ! (3D 4node)
4 & ! (3D 8node)
] , pInt )
integer ( pInt ) , dimension ( FE_Ngeomtypes ) , parameter , public :: FE_Nips = & !< number of IPs in a specific type of element
int ( [ &
1 , & ! element 6 (2D 3node 1ip)
3 , & ! element 125 (2D 6node 3ip)
4 , & ! element 11 (2D 4node 4ip)
9 , & ! element 27 (2D 8node 9ip)
1 , & ! element 134 (3D 4node 1ip)
4 , & ! element 127 (3D 10node 4ip)
6 , & ! element 136 (3D 6node 6ip)
1 , & ! element 117 (3D 8node 1ip)
8 , & ! element 7 (3D 8node 8ip)
27 & ! element 21 (3D 20node 27ip)
] , pInt )
integer ( pInt ) , dimension ( FE_Ncelltypes ) , parameter , public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type
int ( [ &
3 , & ! (2D 3node)
4 , & ! (2D 4node)
4 , & ! (3D 4node)
6 & ! (3D 8node)
] , pInt )
integer ( pInt ) , dimension ( FE_Ngeomtypes ) , parameter , private :: FE_maxNnodesAtIP = & !< maximum number of parent nodes that belong to an IP for a specific type of element
int ( [ &
3 , & ! element 6 (2D 3node 1ip)
1 , & ! element 125 (2D 6node 3ip)
1 , & ! element 11 (2D 4node 4ip)
2 , & ! element 27 (2D 8node 9ip)
4 , & ! element 134 (3D 4node 1ip)
1 , & ! element 127 (3D 10node 4ip)
1 , & ! element 136 (3D 6node 6ip)
8 , & ! element 117 (3D 8node 1ip)
1 , & ! element 7 (3D 8node 8ip)
4 & ! element 21 (3D 20node 27ip)
] , pInt )
2018-09-23 19:01:19 +05:30
integer ( pInt ) , private :: &
2018-09-23 20:35:01 +05:30
mesh_Nelems , & !< total number of elements in mesh (including non-DAMASK elements)
2018-09-23 19:27:21 +05:30
mesh_maxNnodes , & !< max number of nodes in any CP element
2018-09-23 19:01:19 +05:30
mesh_NelemSets
2018-09-23 18:57:51 +05:30
character ( len = 64 ) , dimension ( : ) , allocatable , private :: &
mesh_nameElemSet , & !< names of elementSet
mesh_nameMaterial , & !< names of material in solid section
mesh_mapMaterial !< name of elementSet for material
integer ( pInt ) , dimension ( : , : ) , allocatable , private :: &
mesh_mapElemSet !< list of elements in elementSet
integer ( pInt ) , dimension ( : , : ) , allocatable , target , private :: &
mesh_mapFEtoCPelem , & !< [sorted FEid, corresponding CPid]
mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid]
logical , private :: noPart !< for cases where the ABAQUS input file does not use part/assembly information
2013-04-22 00:18:59 +05:30
2013-04-18 22:10:49 +05:30
public :: &
mesh_init , &
2013-04-22 00:18:59 +05:30
mesh_build_cellnodes , &
2013-04-18 22:10:49 +05:30
mesh_build_ipVolumes , &
mesh_build_ipCoordinates , &
2013-04-22 00:18:59 +05:30
mesh_cellCenterCoordinates , &
2018-09-23 21:27:48 +05:30
mesh_FEasCP
2012-06-15 21:40:21 +05:30
private :: &
2018-09-23 21:27:48 +05:30
mesh_get_damaskOptions , &
mesh_build_cellconnectivity , &
mesh_build_ipAreas , &
FE_mapElemtype , &
mesh_faceMatch , &
mesh_build_FEdata , &
mesh_build_nodeTwins , &
mesh_build_sharedElems , &
mesh_build_ipNeighborhood , &
2013-04-18 22:10:49 +05:30
mesh_abaqus_count_nodesAndElements , &
mesh_abaqus_count_elementSets , &
mesh_abaqus_count_materials , &
mesh_abaqus_map_elementSets , &
mesh_abaqus_map_materials , &
mesh_abaqus_count_cpElements , &
mesh_abaqus_map_elements , &
mesh_abaqus_map_nodes , &
mesh_abaqus_build_nodes , &
mesh_abaqus_count_cpSizes , &
2018-09-23 21:27:48 +05:30
mesh_abaqus_build_elements
2013-04-18 22:10:49 +05:30
2019-01-28 18:29:54 +05:30
2019-02-01 16:54:23 +05:30
type , public , extends ( tMesh ) :: tMesh_abaqus
2019-01-28 18:29:54 +05:30
integer ( pInt ) :: &
mesh_Nelems , & !< total number of elements in mesh (including non-DAMASK elements)
mesh_maxNnodes , & !< max number of nodes in any CP element
mesh_NelemSets , &
mesh_maxNelemInSet , &
mesh_Nmaterials
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
logical :: noPart !< for cases where the ABAQUS input file does not use part/assembly information
contains
2019-02-01 16:54:23 +05:30
procedure , pass ( self ) :: tMesh_abaqus_init
generic , public :: init = > tMesh_abaqus_init
end type tMesh_abaqus
2019-01-28 18:29:54 +05:30
2019-02-01 16:54:23 +05:30
type ( tMesh_abaqus ) , public , protected :: theMesh
2019-01-28 18:29:54 +05:30
2012-03-09 01:55:28 +05:30
contains
2007-03-21 21:48:33 +05:30
2019-02-01 16:54:23 +05:30
subroutine tMesh_abaqus_init ( self , elemType , nodes )
2019-01-28 18:29:54 +05:30
implicit none
class ( tMesh_abaqus ) :: self
2019-02-01 16:54:23 +05:30
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: nodes
integer ( pInt ) , intent ( in ) :: elemType
call self % tMesh % init ( 'mesh' , elemType , nodes )
2019-02-02 16:16:38 +05:30
call theMesh % setNelems ( mesh_NcpElems )
2019-01-28 18:29:54 +05:30
end subroutine tMesh_abaqus_init
2009-01-20 00:12:31 +05:30
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief initializes the mesh by calling all necessary private routines the mesh module
!! Order and routines strongly depend on type of solver
!--------------------------------------------------------------------------------------------------
2013-02-05 18:57:37 +05:30
subroutine mesh_init ( ip , el )
2018-02-02 17:06:09 +05:30
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
2017-10-05 20:05:34 +05:30
use , intrinsic :: iso_fortran_env , only : &
compiler_version , &
compiler_options
#endif
2017-02-05 04:23:22 +05:30
use DAMASK_interface
2012-06-15 21:40:21 +05:30
use IO , only : &
2015-04-11 10:39:15 +05:30
IO_open_InputFile , &
IO_timeStamp , &
IO_error , &
IO_write_jobFile
2013-10-16 18:08:00 +05:30
use debug , only : &
debug_e , &
2014-12-19 00:11:02 +05:30
debug_i , &
debug_level , &
debug_mesh , &
debug_levelBasic
2013-04-16 22:37:27 +05:30
use numerics , only : &
2013-05-07 18:36:29 +05:30
usePingPong , &
2014-12-19 00:11:02 +05:30
numerics_unitlength , &
worldrank
2012-06-15 21:40:21 +05:30
use FEsolving , only : &
2019-01-28 19:06:44 +05:30
modelName , &
calcMode , & FEsolving_execElem , &
2018-09-23 22:34:17 +05:30
FEsolving_execIP
2018-01-10 21:43:25 +05:30
2007-04-10 16:52:53 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , parameter :: FILEUNIT = 222_pInt
2018-09-24 01:01:30 +05:30
integer ( pInt ) , intent ( in ) , optional :: el , ip
2013-02-08 21:25:53 +05:30
integer ( pInt ) :: j
2014-12-19 00:11:02 +05:30
logical :: myDebug
2018-01-10 21:43:25 +05:30
2016-06-29 20:29:22 +05:30
write ( 6 , '(/,a)' ) ' <<<+- mesh init -+>>>'
write ( 6 , '(a15,a)' ) ' Current time: ' , IO_timeStamp ( )
2012-02-01 00:48:55 +05:30
#include "compilation_info.f90"
2013-01-08 16:39:20 +05:30
2012-03-09 01:55:28 +05:30
call mesh_build_FEdata ! get properties of the different types of elements
2013-05-07 18:36:29 +05:30
mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh
2018-01-10 21:43:25 +05:30
2014-12-19 00:11:02 +05:30
myDebug = ( iand ( debug_level ( debug_mesh ) , debug_levelBasic ) / = 0_pInt )
2013-05-07 18:36:29 +05:30
2013-12-11 22:19:20 +05:30
call IO_open_inputFile ( FILEUNIT , modelName ) ! parse info from input file...
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Opened input file' ; flush ( 6 )
2019-01-31 15:59:56 +05:30
noPart = hasNoPart ( FILEUNIT )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_count_nodesAndElements ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Counted nodes/elements' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_count_elementSets ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Counted element sets' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_count_materials ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Counted materials' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_map_elementSets ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Mapped element sets' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_map_materials ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Mapped materials' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_count_cpElements ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Counted CP elements' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_map_elements ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Mapped elements' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_map_nodes ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Mapped nodes' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_build_nodes ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Built nodes' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_count_cpSizes ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Counted CP sizes' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_abaqus_build_elements ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Built elements' ; flush ( 6 )
2013-12-11 22:19:20 +05:30
call mesh_get_damaskOptions ( FILEUNIT )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Got DAMASK options' ; flush ( 6 )
2013-04-22 00:18:59 +05:30
call mesh_build_cellconnectivity
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Built cell connectivity' ; flush ( 6 )
2013-04-24 00:00:56 +05:30
mesh_cellnode = mesh_build_cellnodes ( mesh_node , mesh_Ncellnodes )
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Built cell nodes' ; flush ( 6 )
2012-03-09 01:55:28 +05:30
call mesh_build_ipCoordinates
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Built IP coordinates' ; flush ( 6 )
2012-03-09 01:55:28 +05:30
call mesh_build_ipVolumes
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Built IP volumes' ; flush ( 6 )
2012-03-09 01:55:28 +05:30
call mesh_build_ipAreas
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Built IP areas' ; flush ( 6 )
2015-11-26 02:25:17 +05:30
close ( FILEUNIT )
2012-03-09 01:55:28 +05:30
call mesh_build_nodeTwins
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Built node twins' ; flush ( 6 )
2012-03-09 01:55:28 +05:30
call mesh_build_sharedElems
2014-12-19 00:11:02 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Built shared elements' ; flush ( 6 )
2012-03-09 01:55:28 +05:30
call mesh_build_ipNeighborhood
2015-11-26 02:25:17 +05:30
if ( myDebug ) write ( 6 , '(a)' ) ' Built IP neighborhood' ; flush ( 6 )
2012-02-13 23:11:27 +05:30
2013-10-16 18:08:00 +05:30
if ( usePingPong . and . ( mesh_Nelems / = mesh_NcpElems ) ) &
2014-01-17 07:08:35 +05:30
call IO_error ( 600_pInt ) ! ping-pong must be disabled when having non-DAMASK elements
2013-10-16 18:08:00 +05:30
if ( debug_e < 1 . or . debug_e > mesh_NcpElems ) &
call IO_error ( 602_pInt , ext_msg = 'element' ) ! selected element does not exist
if ( debug_i < 1 . or . debug_i > FE_Nips ( FE_geomtype ( mesh_element ( 2_pInt , debug_e ) ) ) ) &
call IO_error ( 602_pInt , ext_msg = 'IP' ) ! selected element does not have requested IP
2014-01-17 07:08:35 +05:30
FEsolving_execElem = [ 1_pInt , mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements
2018-09-23 20:56:13 +05:30
allocate ( FEsolving_execIP ( 2_pInt , mesh_NcpElems ) , source = 1_pInt ) ! parallel loop bounds set to comprise from first IP...
2014-01-17 07:08:35 +05:30
forall ( j = 1_pInt : mesh_NcpElems ) FEsolving_execIP ( 2 , j ) = FE_Nips ( FE_geomtype ( mesh_element ( 2 , j ) ) ) ! ...up to own IP count for each element
2009-10-12 21:31:49 +05:30
allocate ( calcMode ( mesh_maxNips , mesh_NcpElems ) )
2012-06-15 21:40:21 +05:30
calcMode = . false . ! pretend to have collected what first call is asking (F = I)
2013-02-05 18:57:37 +05:30
calcMode ( ip , mesh_FEasCP ( 'elem' , el ) ) = . true . ! first ip,el needs to be already pingponged to "calc"
2013-04-19 18:11:06 +05:30
2018-09-23 21:07:57 +05:30
!!!! COMPATIBILITY HACK !!!!
! for a homogeneous mesh, all elements have the same number of IPs and and cell nodes.
! hence, xxPerElem instead of maxXX
mesh_NcellnodesPerElem = mesh_maxNcellnodes
2018-09-23 22:20:54 +05:30
! better name
2018-10-04 09:33:48 +05:30
mesh_homogenizationAt = mesh_element ( 3 , : )
mesh_microstructureAt = mesh_element ( 4 , : )
2018-09-23 21:07:57 +05:30
!!!!!!!!!!!!!!!!!!!!!!!!
2019-02-01 16:54:23 +05:30
call theMesh % init ( mesh_element ( 2 , 1 ) , mesh_node0 )
2019-01-31 15:59:56 +05:30
contains
2019-02-02 13:59:58 +05:30
2019-01-31 15:59:56 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief check if the input file for Abaqus contains part info
!--------------------------------------------------------------------------------------------------
logical function hasNoPart ( fileUnit )
2019-01-31 21:31:26 +05:30
use IO , only : &
IO_stringPos , &
IO_stringValue , &
IO_lc
2019-01-31 15:59:56 +05:30
implicit none
integer ( pInt ) , intent ( in ) :: fileUnit
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
character ( len = 65536 ) :: line
hasNoPart = . true .
rewind ( fileUnit )
do
2019-02-02 13:59:58 +05:30
read ( fileUnit , '(a65536)' , END = 620 ) line
2019-01-31 15:59:56 +05:30
chunkPos = IO_stringPos ( line )
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*part' ) then
hasNoPart = . false .
exit
endif
enddo
620 end function hasNoPart
2012-03-09 01:55:28 +05:30
end subroutine mesh_init
2009-01-20 00:12:31 +05:30
2014-04-04 20:10:30 +05:30
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Gives the FE to CP ID mapping by binary search through lookup array
!! valid questions (what) are 'elem', 'node'
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
integer ( pInt ) function mesh_FEasCP ( what , myID )
2013-04-10 15:49:16 +05:30
use IO , only : &
IO_lc
2012-03-09 01:55:28 +05:30
2007-04-10 16:52:53 +05:30
implicit none
character ( len = * ) , intent ( in ) :: what
2012-03-09 01:55:28 +05:30
integer ( pInt ) , intent ( in ) :: myID
2018-01-10 21:43:25 +05:30
2007-04-10 16:52:53 +05:30
integer ( pInt ) , dimension ( : , : ) , pointer :: lookupMap
2012-03-09 01:55:28 +05:30
integer ( pInt ) :: lower , upper , center
2018-01-10 21:43:25 +05:30
2007-04-10 16:52:53 +05:30
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
2018-01-10 21:43:25 +05:30
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 )
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
if ( lookupMap ( 1_pInt , lower ) == myID ) then ! 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
mesh_FEasCP = lookupMap ( 2_pInt , lower )
2007-04-10 16:52:53 +05:30
return
2012-03-09 01:55:28 +05:30
elseif ( lookupMap ( 1_pInt , upper ) == myID ) then
2012-02-16 00:28:38 +05:30
mesh_FEasCP = lookupMap ( 2_pInt , upper )
2007-04-10 16:52:53 +05:30
return
endif
2013-10-16 18:34:59 +05:30
binarySearch : do while ( upper - lower > 1_pInt )
2012-02-16 00:28:38 +05:30
center = ( lower + upper ) / 2_pInt
2012-03-09 01:55:28 +05:30
if ( lookupMap ( 1_pInt , center ) < myID ) then
2007-04-10 16:52:53 +05:30
lower = center
2012-03-09 01:55:28 +05:30
elseif ( lookupMap ( 1_pInt , center ) > myID ) 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
2013-10-16 18:34:59 +05:30
enddo binarySearch
2018-01-10 21:43:25 +05:30
2012-03-09 01:55:28 +05:30
end function mesh_FEasCP
2019-01-28 18:29:54 +05:30
2009-01-20 00:12:31 +05:30
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
2013-04-15 13:43:20 +05:30
!> @brief Split CP elements into cells.
2018-01-10 21:43:25 +05:30
!> @details Build a mapping between cells and the corresponding cell nodes ('mesh_cell').
!> Cell nodes that are also matching nodes are unique in the list of cell nodes,
2013-04-22 00:18:59 +05:30
!> all others (currently) might be stored more than once.
!> Also allocates the 'mesh_node' array.
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
2013-04-22 00:18:59 +05:30
subroutine mesh_build_cellconnectivity
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
implicit none
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( : ) , allocatable :: &
matchingNode2cellnode
integer ( pInt ) , dimension ( : , : ) , allocatable :: &
cellnodeParent
integer ( pInt ) , dimension ( mesh_maxNcellnodes ) :: &
localCellnode2globalCellnode
2013-10-16 18:08:00 +05:30
integer ( pInt ) :: &
2018-01-10 21:43:25 +05:30
e , t , g , c , n , i , &
2013-04-15 13:43:20 +05:30
matchingNodeID , &
localCellnodeID
2018-01-10 21:43:25 +05:30
2013-12-30 15:36:01 +05:30
allocate ( mesh_cell ( FE_maxNcellnodesPerCell , mesh_maxNips , mesh_NcpElems ) , source = 0_pInt )
allocate ( matchingNode2cellnode ( mesh_Nnodes ) , source = 0_pInt )
allocate ( cellnodeParent ( 2_pInt , mesh_maxNcellnodes * mesh_NcpElems ) , source = 0_pInt )
2018-01-10 21:43:25 +05:30
2013-10-16 18:08:00 +05:30
!--------------------------------------------------------------------------------------------------
! Count cell nodes (including duplicates) and generate cell connectivity list
2013-04-22 00:18:59 +05:30
mesh_Ncellnodes = 0_pInt
mesh_Ncells = 0_pInt
do e = 1_pInt , mesh_NcpElems ! loop over cpElems
t = mesh_element ( 2_pInt , e ) ! get element type
g = FE_geomtype ( t ) ! get geometry type
c = FE_celltype ( g ) ! get cell type
localCellnode2globalCellnode = 0_pInt
mesh_Ncells = mesh_Ncells + FE_Nips ( g )
do i = 1_pInt , FE_Nips ( g ) ! loop over ips=cells in this element
do n = 1_pInt , FE_NcellnodesPerCell ( c ) ! loop over cell nodes in this cell
localCellnodeID = FE_cell ( n , i , g )
if ( localCellnodeID < = FE_NmatchingNodes ( g ) ) then ! this cell node is a matching node
matchingNodeID = mesh_element ( 4_pInt + localCellnodeID , e )
if ( matchingNode2cellnode ( matchingNodeID ) == 0_pInt ) then ! if this matching node does not yet exist in the glbal cell node list ...
mesh_Ncellnodes = mesh_Ncellnodes + 1_pInt ! ... count it as cell node ...
matchingNode2cellnode ( matchingNodeID ) = mesh_Ncellnodes ! ... and remember its global ID
cellnodeParent ( 1_pInt , mesh_Ncellnodes ) = e ! ... and where it belongs to
cellnodeParent ( 2_pInt , mesh_Ncellnodes ) = localCellnodeID
2013-04-15 13:43:20 +05:30
endif
2013-04-22 00:18:59 +05:30
mesh_cell ( n , i , e ) = matchingNode2cellnode ( matchingNodeID )
else ! this cell node is no matching node
if ( localCellnode2globalCellnode ( localCellnodeID ) == 0_pInt ) then ! if this local cell node does not yet exist in the global cell node list ...
mesh_Ncellnodes = mesh_Ncellnodes + 1_pInt ! ... count it as cell node ...
localCellnode2globalCellnode ( localCellnodeID ) = mesh_Ncellnodes ! ... and remember its global ID ...
cellnodeParent ( 1_pInt , mesh_Ncellnodes ) = e ! ... and it belongs to
cellnodeParent ( 2_pInt , mesh_Ncellnodes ) = localCellnodeID
endif
mesh_cell ( n , i , e ) = localCellnode2globalCellnode ( localCellnodeID )
endif
2012-06-15 21:40:21 +05:30
enddo
enddo
2013-04-22 00:18:59 +05:30
enddo
2013-04-15 13:43:20 +05:30
2013-04-22 00:18:59 +05:30
allocate ( mesh_cellnodeParent ( 2_pInt , mesh_Ncellnodes ) )
allocate ( mesh_cellnode ( 3_pInt , mesh_Ncellnodes ) )
forall ( n = 1_pInt : mesh_Ncellnodes )
mesh_cellnodeParent ( 1 , n ) = cellnodeParent ( 1 , n )
mesh_cellnodeParent ( 2 , n ) = cellnodeParent ( 2 , n )
endforall
2013-04-15 13:43:20 +05:30
2013-04-22 00:18:59 +05:30
end subroutine mesh_build_cellconnectivity
!--------------------------------------------------------------------------------------------------
!> @brief Calculate position of cellnodes from the given position of nodes
2018-01-10 21:43:25 +05:30
!> Build list of cellnodes' coordinates.
2013-04-22 00:18:59 +05:30
!> Cellnode coordinates are calculated from a weighted sum of node coordinates.
!--------------------------------------------------------------------------------------------------
2013-04-24 00:00:56 +05:30
function mesh_build_cellnodes ( nodes , Ncellnodes )
2018-01-10 21:43:25 +05:30
2013-04-22 00:18:59 +05:30
implicit none
2013-10-16 18:08:00 +05:30
integer ( pInt ) , intent ( in ) :: Ncellnodes !< requested number of cellnodes
2013-04-22 00:18:59 +05:30
real ( pReal ) , dimension ( 3 , mesh_Nnodes ) , intent ( in ) :: nodes
2013-04-24 00:00:56 +05:30
real ( pReal ) , dimension ( 3 , Ncellnodes ) :: mesh_build_cellnodes
2013-04-22 00:18:59 +05:30
2013-10-16 18:08:00 +05:30
integer ( pInt ) :: &
2018-01-10 21:43:25 +05:30
e , t , n , m , &
2013-04-22 00:18:59 +05:30
localCellnodeID
real ( pReal ) , dimension ( 3 ) :: &
myCoords
2018-01-10 21:43:25 +05:30
2013-04-22 00:18:59 +05:30
mesh_build_cellnodes = 0.0_pReal
2013-04-18 23:26:14 +05:30
!$OMP PARALLEL DO PRIVATE(e,localCellnodeID,t,myCoords)
2013-04-24 00:00:56 +05:30
do n = 1_pInt , Ncellnodes ! loop over cell nodes
2013-04-22 00:18:59 +05:30
e = mesh_cellnodeParent ( 1 , n )
localCellnodeID = mesh_cellnodeParent ( 2 , n )
2013-04-15 13:43:20 +05:30
t = mesh_element ( 2 , e ) ! get element type
myCoords = 0.0_pReal
do m = 1_pInt , FE_Nnodes ( t )
2013-04-22 00:18:59 +05:30
myCoords = myCoords + nodes ( 1 : 3 , mesh_element ( 4_pInt + m , e ) ) &
2013-04-15 13:43:20 +05:30
* FE_cellnodeParentnodeWeights ( m , localCellnodeID , t )
enddo
2013-04-22 00:18:59 +05:30
mesh_build_cellnodes ( 1 : 3 , n ) = myCoords / sum ( FE_cellnodeParentnodeWeights ( : , localCellnodeID , t ) )
2013-04-15 13:43:20 +05:30
enddo
2013-04-18 23:26:14 +05:30
!$OMP END PARALLEL DO
2013-04-15 13:43:20 +05:30
2013-04-22 00:18:59 +05:30
end function mesh_build_cellnodes
2011-02-16 21:53:08 +05:30
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
2013-04-15 13:43:20 +05:30
!> @details The IP volume is calculated differently depending on the cell type.
!> 2D cells assume an element depth of one in order to calculate the volume.
!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal
!> shape with a cell face as basis and the central ip at the tip. This subvolume is
2018-01-10 21:43:25 +05:30
!> calculated as an average of four tetrahedals with three corners on the cell face
2013-04-15 13:43:20 +05:30
!> and one corner at the central ip.
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipVolumes
2013-04-10 15:49:16 +05:30
use math , only : &
math_volTetrahedron , &
math_areaTriangle
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
implicit none
2013-04-15 13:43:20 +05:30
integer ( pInt ) :: e , t , g , c , i , m , f , n
real ( pReal ) , dimension ( FE_maxNcellnodesPerCellface , FE_maxNcellfaces ) :: subvolume
2011-02-16 21:53:08 +05:30
2019-02-01 15:52:06 +05:30
allocate ( mesh_ipVolume ( mesh_maxNips , mesh_NcpElems ) , source = 0.0_pReal )
2018-01-10 21:43:25 +05:30
2013-04-15 13:43:20 +05:30
!$OMP PARALLEL DO PRIVATE(t,g,c,m,subvolume)
do e = 1_pInt , mesh_NcpElems ! loop over cpElems
t = mesh_element ( 2_pInt , e ) ! get element type
g = FE_geomtype ( t ) ! get geometry type
c = FE_celltype ( g ) ! get cell type
select case ( c )
case ( 1_pInt ) ! 2D 3node
forall ( i = 1_pInt : FE_Nips ( g ) ) & ! loop over ips=cells in this element
2013-04-22 00:18:59 +05:30
mesh_ipVolume ( i , e ) = math_areaTriangle ( mesh_cellnode ( 1 : 3 , mesh_cell ( 1 , i , e ) ) , &
mesh_cellnode ( 1 : 3 , mesh_cell ( 2 , i , e ) ) , &
mesh_cellnode ( 1 : 3 , mesh_cell ( 3 , i , e ) ) )
2018-01-10 21:43:25 +05:30
2013-04-15 13:43:20 +05:30
case ( 2_pInt ) ! 2D 4node
forall ( i = 1_pInt : FE_Nips ( g ) ) & ! loop over ips=cells in this element
2013-04-22 00:18:59 +05:30
mesh_ipVolume ( i , e ) = math_areaTriangle ( mesh_cellnode ( 1 : 3 , mesh_cell ( 1 , i , e ) ) , & ! here we assume a planar shape, so division in two triangles suffices
mesh_cellnode ( 1 : 3 , mesh_cell ( 2 , i , e ) ) , &
mesh_cellnode ( 1 : 3 , mesh_cell ( 3 , i , e ) ) ) &
+ math_areaTriangle ( mesh_cellnode ( 1 : 3 , mesh_cell ( 3 , i , e ) ) , &
mesh_cellnode ( 1 : 3 , mesh_cell ( 4 , i , e ) ) , &
mesh_cellnode ( 1 : 3 , mesh_cell ( 1 , i , e ) ) )
2013-04-15 13:43:20 +05:30
case ( 3_pInt ) ! 3D 4node
forall ( i = 1_pInt : FE_Nips ( g ) ) & ! loop over ips=cells in this element
2013-04-22 00:18:59 +05:30
mesh_ipVolume ( i , e ) = math_volTetrahedron ( mesh_cellnode ( 1 : 3 , mesh_cell ( 1 , i , e ) ) , &
mesh_cellnode ( 1 : 3 , mesh_cell ( 2 , i , e ) ) , &
mesh_cellnode ( 1 : 3 , mesh_cell ( 3 , i , e ) ) , &
mesh_cellnode ( 1 : 3 , mesh_cell ( 4 , i , e ) ) )
2013-04-15 13:43:20 +05:30
2019-01-28 18:29:54 +05:30
case ( 4_pInt ) ! 3D 8node
m = FE_NcellnodesPerCellface ( c )
do i = 1_pInt , FE_Nips ( g ) ! loop over ips=cells in this element
subvolume = 0.0_pReal
forall ( f = 1_pInt : FE_NipNeighbors ( c ) , n = 1_pInt : FE_NcellnodesPerCellface ( c ) ) &
subvolume ( n , f ) = math_volTetrahedron ( &
mesh_cellnode ( 1 : 3 , mesh_cell ( FE_cellface ( n , f , c ) , i , e ) ) , &
mesh_cellnode ( 1 : 3 , mesh_cell ( FE_cellface ( 1 + mod ( n , m ) , f , c ) , i , e ) ) , &
mesh_cellnode ( 1 : 3 , mesh_cell ( FE_cellface ( 1 + mod ( n + 1 , m ) , f , c ) , i , e ) ) , &
mesh_ipCoordinates ( 1 : 3 , i , e ) )
mesh_ipVolume ( i , e ) = 0.5_pReal * sum ( subvolume ) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two
enddo
end select
enddo
!$OMP END PARALLEL DO
2018-01-10 21:43:25 +05:30
2019-01-28 18:29:54 +05:30
end subroutine mesh_build_ipVolumes
2009-10-12 21:31:49 +05:30
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-28 18:29:54 +05:30
!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates'
! Called by all solvers in mesh_init in order to initialize the ip coordinates.
! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus,
! so no need to use this subroutine anymore; Marc however only provides nodal displacements,
! so in this case the ip coordinates are always calculated on the basis of this subroutine.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES,
! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME.
! HAS TO BE CHANGED IN A LATER VERSION.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-28 18:29:54 +05:30
subroutine mesh_build_ipCoordinates
2009-10-12 21:31:49 +05:30
implicit none
2019-01-28 18:29:54 +05:30
integer ( pInt ) :: e , t , g , c , i , n
real ( pReal ) , dimension ( 3 ) :: myCoords
2009-10-12 21:31:49 +05:30
2019-01-28 18:29:54 +05:30
if ( . not . allocated ( mesh_ipCoordinates ) ) &
allocate ( mesh_ipCoordinates ( 3 , mesh_maxNips , mesh_NcpElems ) , source = 0.0_pReal )
2009-10-12 21:31:49 +05:30
2019-01-28 18:29:54 +05:30
!$OMP PARALLEL DO PRIVATE(t,g,c,myCoords)
do e = 1_pInt , mesh_NcpElems ! loop over cpElems
t = mesh_element ( 2_pInt , e ) ! get element type
g = FE_geomtype ( t ) ! get geometry type
c = FE_celltype ( g ) ! get cell type
do i = 1_pInt , FE_Nips ( g ) ! loop over ips=cells in this element
myCoords = 0.0_pReal
do n = 1_pInt , FE_NcellnodesPerCell ( c ) ! loop over cell nodes in this cell
myCoords = myCoords + mesh_cellnode ( 1 : 3 , mesh_cell ( n , i , e ) )
2012-06-15 21:40:21 +05:30
enddo
2019-01-28 18:29:54 +05:30
mesh_ipCoordinates ( 1 : 3 , i , e ) = myCoords / real ( FE_NcellnodesPerCell ( c ) , pReal )
enddo
2009-10-12 21:31:49 +05:30
enddo
2019-01-28 18:29:54 +05:30
!$OMP END PARALLEL DO
2009-10-12 21:31:49 +05:30
2019-01-28 18:29:54 +05:30
end subroutine mesh_build_ipCoordinates
2009-10-12 21:31:49 +05:30
2012-03-09 01:55:28 +05:30
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-28 18:29:54 +05:30
!> @brief Calculates cell center coordinates.
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-28 18:29:54 +05:30
pure function mesh_cellCenterCoordinates ( ip , el )
2009-10-12 21:31:49 +05:30
implicit none
2019-01-28 18:29:54 +05:30
integer ( pInt ) , intent ( in ) :: el , & !< element number
ip !< integration point number
real ( pReal ) , dimension ( 3 ) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell
integer ( pInt ) :: t , g , c , n
2018-01-10 21:43:25 +05:30
2019-01-28 18:29:54 +05:30
t = mesh_element ( 2_pInt , el ) ! get element type
g = FE_geomtype ( t ) ! get geometry type
c = FE_celltype ( g ) ! get cell type
mesh_cellCenterCoordinates = 0.0_pReal
do n = 1_pInt , FE_NcellnodesPerCell ( c ) ! loop over cell nodes in this cell
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode ( 1 : 3 , mesh_cell ( n , ip , el ) )
2009-10-12 21:31:49 +05:30
enddo
2019-01-28 18:29:54 +05:30
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / real ( FE_NcellnodesPerCell ( c ) , pReal )
end function mesh_cellCenterCoordinates
2018-01-10 21:43:25 +05:30
2010-05-06 22:10:47 +05:30
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Count overall number of nodes and elements in mesh and stores them in
!! 'mesh_Nelems' and 'mesh_Nnodes'
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_count_nodesAndElements ( fileUnit )
2009-10-12 21:31:49 +05:30
2012-03-09 01:55:28 +05:30
use IO , only : IO_lc , &
IO_stringValue , &
IO_stringPos , &
IO_countDataLines , &
IO_error
2018-01-10 21:43:25 +05:30
2009-10-12 21:31:49 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2018-01-10 21:43:25 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2012-06-15 21:40:21 +05:30
character ( len = 300 ) :: line
2019-02-02 14:26:11 +05:30
integer :: myStat
2012-06-15 21:40:21 +05:30
logical :: inPart
2009-10-12 21:31:49 +05:30
2012-06-15 21:40:21 +05:30
mesh_Nnodes = 0_pInt
mesh_Nelems = 0_pInt
2019-02-02 14:26:11 +05:30
2012-06-15 21:40:21 +05:30
inPart = . false .
2019-02-02 14:26:11 +05:30
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) == 'part' ) inPart = . false .
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
if ( inPart . or . noPart ) then
2015-08-28 13:08:48 +05:30
select case ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) )
2012-06-15 21:40:21 +05:30
case ( '*node' )
if ( &
2015-08-28 13:08:48 +05:30
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'print' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'file' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'response' &
2012-06-15 21:40:21 +05:30
) &
2013-12-11 22:19:20 +05:30
mesh_Nnodes = mesh_Nnodes + IO_countDataLines ( fileUnit )
2012-06-15 21:40:21 +05:30
case ( '*element' )
if ( &
2015-08-28 13:08:48 +05:30
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'matrix' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'response' &
2012-06-15 21:40:21 +05:30
) then
2013-12-11 22:19:20 +05:30
mesh_Nelems = mesh_Nelems + IO_countDataLines ( fileUnit )
2012-06-15 21:40:21 +05:30
endif
endselect
2009-10-12 21:31:49 +05:30
endif
enddo
2018-01-10 21:43:25 +05:30
2019-02-02 14:26:11 +05:30
if ( mesh_Nnodes < 2_pInt ) call IO_error ( error_ID = 900_pInt )
2012-06-15 21:40:21 +05:30
if ( mesh_Nelems == 0_pInt ) call IO_error ( error_ID = 901_pInt )
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
end subroutine mesh_abaqus_count_nodesAndElements
2009-10-12 21:31:49 +05:30
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief count overall number of element sets in mesh and write 'mesh_NelemSets' and
!! 'mesh_maxNelemInSet'
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_count_elementSets ( fileUnit )
2012-06-15 21:40:21 +05:30
use IO , only : IO_lc , &
IO_stringValue , &
IO_stringPos , &
IO_error
2010-05-06 22:10:47 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2010-05-06 22:10:47 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2012-06-15 21:40:21 +05:30
character ( len = 300 ) :: line
2019-02-02 14:26:11 +05:30
integer :: myStat
2012-06-15 21:40:21 +05:30
logical :: inPart
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
mesh_NelemSets = 0_pInt
2013-04-18 22:10:49 +05:30
mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
inPart = . false .
2019-02-02 14:26:11 +05:30
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) == 'part' ) inPart = . false .
2018-01-10 21:43:25 +05:30
2015-08-28 13:08:48 +05:30
if ( ( inPart . or . noPart ) . and . IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*elset' ) &
2012-06-15 21:40:21 +05:30
mesh_NelemSets = mesh_NelemSets + 1_pInt
enddo
2010-07-01 20:50:39 +05:30
2012-06-15 21:40:21 +05:30
if ( mesh_NelemSets == 0 ) call IO_error ( error_ID = 902_pInt )
2010-05-06 22:10:47 +05:30
2012-06-15 21:40:21 +05:30
end subroutine mesh_abaqus_count_elementSets
2010-05-06 22:10:47 +05:30
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
2012-06-15 21:40:21 +05:30
! count overall number of solid sections sets in mesh (Abaqus only)
2009-10-12 21:31:49 +05:30
!
2012-06-15 21:40:21 +05:30
! mesh_Nmaterials
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_count_materials ( fileUnit )
2009-10-12 21:31:49 +05:30
2012-03-09 01:55:28 +05:30
use IO , only : IO_lc , &
IO_stringValue , &
IO_stringPos , &
2012-06-15 21:40:21 +05:30
IO_error
2009-10-12 21:31:49 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2018-01-10 21:43:25 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2012-06-15 21:40:21 +05:30
character ( len = 300 ) :: line
2019-02-02 14:26:11 +05:30
integer :: myStat
logical :: inPart
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
mesh_Nmaterials = 0_pInt
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
inPart = . false .
2019-02-02 14:26:11 +05:30
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) == 'part' ) inPart = . false .
2009-10-12 21:31:49 +05:30
2012-06-15 21:40:21 +05:30
if ( ( inPart . or . noPart ) . and . &
2015-08-28 13:08:48 +05:30
IO_lc ( IO_StringValue ( line , chunkPos , 1_pInt ) ) == '*solid' . and . &
IO_lc ( IO_StringValue ( line , chunkPos , 2_pInt ) ) == 'section' ) &
2012-06-15 21:40:21 +05:30
mesh_Nmaterials = mesh_Nmaterials + 1_pInt
enddo
2019-02-02 14:26:11 +05:30
if ( mesh_Nmaterials == 0_pInt ) call IO_error ( error_ID = 903_pInt )
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
end subroutine mesh_abaqus_count_materials
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
2018-01-10 21:43:25 +05:30
! Build element set mapping
2012-06-15 21:40:21 +05:30
!
! allocate globals: mesh_nameElemSet, mesh_mapElemSet
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_map_elementSets ( fileUnit )
2012-06-15 21:40:21 +05:30
use IO , only : IO_lc , &
IO_stringValue , &
IO_stringPos , &
IO_extractValue , &
IO_continuousIntValues , &
IO_error
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2012-06-15 21:40:21 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2012-06-15 21:40:21 +05:30
character ( len = 300 ) :: line
2019-02-02 14:26:11 +05:30
integer :: myStat
logical :: inPart
integer ( pInt ) :: elemSet , i
2012-06-15 21:40:21 +05:30
2018-09-23 20:56:13 +05:30
allocate ( mesh_nameElemSet ( mesh_NelemSets ) ) ; mesh_nameElemSet = ''
allocate ( mesh_mapElemSet ( 1_pInt + mesh_maxNelemInSet , mesh_NelemSets ) , source = 0_pInt )
2009-10-12 21:31:49 +05:30
2012-06-15 21:40:21 +05:30
2019-02-02 14:26:11 +05:30
elemSet = 0_pInt
inPart = . false .
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) == 'part' ) inPart = . false .
2018-01-10 21:43:25 +05:30
2015-08-28 13:08:48 +05:30
if ( ( inPart . or . noPart ) . and . IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*elset' ) then
2012-06-15 21:40:21 +05:30
elemSet = elemSet + 1_pInt
2015-08-28 13:08:48 +05:30
mesh_nameElemSet ( elemSet ) = trim ( IO_extractValue ( IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) , 'elset' ) )
2013-12-11 22:19:20 +05:30
mesh_mapElemSet ( : , elemSet ) = IO_continuousIntValues ( fileUnit , mesh_Nelems , mesh_nameElemSet , &
2012-06-15 21:40:21 +05:30
mesh_mapElemSet , elemSet - 1_pInt )
2009-10-12 21:31:49 +05:30
endif
enddo
2019-02-02 14:26:11 +05:30
do i = 1_pInt , elemSet
2012-06-15 21:40:21 +05:30
if ( mesh_mapElemSet ( 1 , i ) == 0_pInt ) call IO_error ( error_ID = 904_pInt , ext_msg = mesh_nameElemSet ( i ) )
enddo
2009-10-12 21:31:49 +05:30
2012-06-15 21:40:21 +05:30
end subroutine mesh_abaqus_map_elementSets
2009-10-12 21:31:49 +05:30
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
2012-06-15 21:40:21 +05:30
! map solid section (Abaqus only)
2009-10-12 21:31:49 +05:30
!
2012-06-15 21:40:21 +05:30
! allocate globals: mesh_nameMaterial, mesh_mapMaterial
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_map_materials ( fileUnit )
2012-06-15 21:40:21 +05:30
use IO , only : IO_lc , &
IO_stringValue , &
IO_stringPos , &
IO_extractValue , &
IO_error
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2012-06-15 21:40:21 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2019-02-02 14:26:11 +05:30
character ( len = 300 ) :: line
integer :: myStat
logical :: inPart
integer ( pInt ) :: i , c
2012-06-15 21:40:21 +05:30
character ( len = 64 ) :: elemSetName , materialName
2018-01-10 21:43:25 +05:30
2018-09-23 20:56:13 +05:30
allocate ( mesh_nameMaterial ( mesh_Nmaterials ) ) ; mesh_nameMaterial = ''
allocate ( mesh_mapMaterial ( mesh_Nmaterials ) ) ; mesh_mapMaterial = ''
2012-06-15 21:40:21 +05:30
2019-02-02 14:26:11 +05:30
c = 0_pInt
inPart = . false .
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) == 'part' ) inPart = . false .
2012-06-15 21:40:21 +05:30
if ( ( inPart . or . noPart ) . and . &
2015-08-28 13:08:48 +05:30
IO_lc ( IO_StringValue ( line , chunkPos , 1_pInt ) ) == '*solid' . and . &
IO_lc ( IO_StringValue ( line , chunkPos , 2_pInt ) ) == 'section' ) then
2012-06-15 21:40:21 +05:30
elemSetName = ''
materialName = ''
2015-08-28 13:08:48 +05:30
do i = 3_pInt , chunkPos ( 1_pInt )
if ( IO_extractValue ( IO_lc ( IO_stringValue ( line , chunkPos , i ) ) , 'elset' ) / = '' ) &
elemSetName = trim ( IO_extractValue ( IO_lc ( IO_stringValue ( line , chunkPos , i ) ) , 'elset' ) )
if ( IO_extractValue ( IO_lc ( IO_stringValue ( line , chunkPos , i ) ) , 'material' ) / = '' ) &
materialName = trim ( IO_extractValue ( IO_lc ( IO_stringValue ( line , chunkPos , i ) ) , 'material' ) )
2012-06-15 21:40:21 +05:30
enddo
if ( elemSetName / = '' . and . materialName / = '' ) then
c = c + 1_pInt
mesh_nameMaterial ( c ) = materialName ! name of material used for this section
mesh_mapMaterial ( c ) = elemSetName ! mapped to respective element set
2018-01-10 21:43:25 +05:30
endif
2012-06-15 21:40:21 +05:30
endif
enddo
2019-02-02 14:26:11 +05:30
if ( c == 0_pInt ) call IO_error ( error_ID = 905_pInt )
2012-06-15 21:40:21 +05:30
do i = 1_pInt , c
if ( mesh_nameMaterial ( i ) == '' . or . mesh_mapMaterial ( i ) == '' ) call IO_error ( error_ID = 905_pInt )
enddo
end subroutine mesh_abaqus_map_materials
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Count overall number of CP elements in mesh and stores them in 'mesh_NcpElems'
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_count_cpElements ( fileUnit )
2012-06-15 21:40:21 +05:30
use IO , only : IO_lc , &
IO_stringValue , &
IO_stringPos , &
IO_error , &
IO_extractValue
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2018-01-10 21:43:25 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2019-02-02 14:26:11 +05:30
character ( len = 300 ) :: line
integer :: myStat
logical :: materialFound
2012-06-15 21:40:21 +05:30
integer ( pInt ) :: i , k
character ( len = 64 ) :: materialName , elemSetName
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
mesh_NcpElems = 0_pInt
2019-02-02 14:26:11 +05:30
materialFound = . false .
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
select case ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) )
2012-06-15 21:40:21 +05:30
case ( '*material' )
2015-08-28 13:08:48 +05:30
materialName = trim ( IO_extractValue ( IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) , 'name' ) ) ! extract name=value
2012-06-15 21:40:21 +05:30
materialFound = materialName / = '' ! valid name?
case ( '*user' )
2015-08-28 13:08:48 +05:30
if ( IO_lc ( IO_StringValue ( line , chunkPos , 2_pInt ) ) == 'material' . and . materialFound ) then
2012-06-15 21:40:21 +05:30
do i = 1_pInt , mesh_Nmaterials ! look thru material names
if ( materialName == mesh_nameMaterial ( i ) ) then ! found one
elemSetName = mesh_mapMaterial ( i ) ! take corresponding elemSet
do k = 1_pInt , mesh_NelemSets ! look thru all elemSet definitions
if ( elemSetName == mesh_nameElemSet ( k ) ) & ! matched?
mesh_NcpElems = mesh_NcpElems + mesh_mapElemSet ( 1 , k ) ! add those elem count
enddo
endif
enddo
materialFound = . false .
endif
endselect
enddo
2018-01-10 21:43:25 +05:30
2019-02-02 14:26:11 +05:30
if ( mesh_NcpElems == 0_pInt ) call IO_error ( error_ID = 906_pInt )
2012-06-15 21:40:21 +05:30
end subroutine mesh_abaqus_count_cpElements
!--------------------------------------------------------------------------------------------------
!> @brief Maps elements from FE ID to internal (consecutive) representation.
!! Allocates global array 'mesh_mapFEtoCPelem'
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_map_elements ( fileUnit )
2009-10-12 21:31:49 +05:30
2013-02-01 12:56:21 +05:30
use math , only : math_qsort
2012-03-09 01:55:28 +05:30
use IO , only : IO_lc , &
IO_stringValue , &
IO_stringPos , &
IO_extractValue , &
IO_error
2018-01-10 21:43:25 +05:30
2009-10-12 21:31:49 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2009-10-12 21:31:49 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2012-03-09 01:55:28 +05:30
character ( len = 300 ) :: line
2019-02-02 14:26:11 +05:30
integer :: myStat
logical :: materialFound
integer ( pInt ) :: i , j , k , cpElem
2012-06-15 21:40:21 +05:30
character ( len = 64 ) materialName , elemSetName ! why limited to 64? ABAQUS?
2009-10-12 21:31:49 +05:30
2018-09-23 20:56:13 +05:30
allocate ( mesh_mapFEtoCPelem ( 2 , mesh_NcpElems ) , source = 0_pInt )
2009-10-12 21:31:49 +05:30
2019-02-02 14:26:11 +05:30
cpElem = 0_pInt
materialFound = . false .
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
select case ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) )
2009-10-12 21:31:49 +05:30
case ( '*material' )
2015-08-28 13:08:48 +05:30
materialName = trim ( IO_extractValue ( IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) , 'name' ) ) ! extract name=value
2012-06-15 21:40:21 +05:30
materialFound = materialName / = '' ! valid name?
2009-10-12 21:31:49 +05:30
case ( '*user' )
2015-08-28 13:08:48 +05:30
if ( IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) == 'material' . and . materialFound ) then
2012-03-09 01:55:28 +05:30
do i = 1_pInt , mesh_Nmaterials ! look thru material names
if ( materialName == mesh_nameMaterial ( i ) ) then ! found one
elemSetName = mesh_mapMaterial ( i ) ! take corresponding elemSet
do k = 1_pInt , mesh_NelemSets ! look thru all elemSet definitions
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-03-09 01:55:28 +05:30
mesh_mapFEtoCPelem ( 1 , cpElem ) = mesh_mapElemSet ( 1_pInt + j , k ) ! store FE id
mesh_mapFEtoCPelem ( 2 , cpElem ) = cpElem ! store our id
2011-03-23 21:50:12 +05:30
enddo
endif
enddo
endif
2009-10-12 21:31:49 +05:30
enddo
materialFound = . false .
endif
endselect
enddo
2019-02-02 14:26:11 +05:30
call math_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
2012-03-09 01:55:28 +05:30
end subroutine mesh_abaqus_map_elements
2009-10-12 21:31:49 +05:30
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Maps node from FE ID to internal (consecutive) representation.
!! Allocates global array 'mesh_mapFEtoCPnode'
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_map_nodes ( fileUnit )
2010-05-06 22:10:47 +05:30
2013-02-01 12:56:21 +05:30
use math , only : math_qsort
2012-03-09 01:55:28 +05:30
use IO , only : IO_lc , &
IO_stringValue , &
IO_stringPos , &
2012-06-15 21:40:21 +05:30
IO_countDataLines , &
2012-03-09 01:55:28 +05:30
IO_intValue , &
2012-06-15 21:40:21 +05:30
IO_error
2012-03-09 01:55:28 +05:30
2009-10-12 21:31:49 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2012-06-15 21:40:21 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2019-02-02 14:26:11 +05:30
character ( len = 300 ) :: line
integer :: myStat
logical :: inPart
integer ( pInt ) :: i , c , cpNode
2012-06-15 21:40:21 +05:30
2018-09-23 20:56:13 +05:30
allocate ( mesh_mapFEtoCPnode ( 2_pInt , mesh_Nnodes ) , source = 0_pInt )
2019-02-02 14:26:11 +05:30
cpNode = 0_pInt
inPart = . false .
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) == 'part' ) inPart = . false .
2012-06-15 21:40:21 +05:30
if ( ( inPart . or . noPart ) . and . &
2015-08-28 13:08:48 +05:30
IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*node' . and . &
( IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'print' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'file' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'response' ) &
2012-06-15 21:40:21 +05:30
) then
2013-12-11 22:19:20 +05:30
c = IO_countDataLines ( fileUnit )
2012-06-15 21:40:21 +05:30
do i = 1_pInt , c
2013-12-11 22:19:20 +05:30
backspace ( fileUnit )
2012-06-15 21:40:21 +05:30
enddo
do i = 1_pInt , c
2019-02-02 21:49:12 +05:30
read ( fileUnit , '(a300)' ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
2012-06-15 21:40:21 +05:30
cpNode = cpNode + 1_pInt
2015-08-28 13:08:48 +05:30
mesh_mapFEtoCPnode ( 1_pInt , cpNode ) = IO_intValue ( line , chunkPos , 1_pInt )
2012-06-15 21:40:21 +05:30
mesh_mapFEtoCPnode ( 2_pInt , cpNode ) = cpNode
2009-10-12 21:31:49 +05:30
enddo
endif
enddo
2012-06-15 21:40:21 +05:30
2019-02-02 14:26:11 +05:30
call math_qsort ( mesh_mapFEtoCPnode , 1_pInt , int ( size ( mesh_mapFEtoCPnode , 2_pInt ) , pInt ) )
2012-06-15 21:40:21 +05:30
if ( int ( size ( mesh_mapFEtoCPnode ) , pInt ) == 0_pInt ) call IO_error ( error_ID = 908_pInt )
end subroutine mesh_abaqus_map_nodes
!--------------------------------------------------------------------------------------------------
!> @brief store x,y,z coordinates of all nodes in mesh.
!! Allocates global arrays 'mesh_node0' and 'mesh_node'
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_build_nodes ( fileUnit )
2013-05-17 22:22:19 +05:30
use IO , only : &
IO_lc , &
IO_stringValue , &
IO_floatValue , &
IO_stringPos , &
IO_error , &
IO_countDataLines , &
IO_intValue
2012-06-15 21:40:21 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2012-06-15 21:40:21 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2012-06-15 21:40:21 +05:30
character ( len = 300 ) :: line
2019-02-02 14:26:11 +05:30
integer :: myStat
2012-06-15 21:40:21 +05:30
logical :: inPart
2019-02-02 14:26:11 +05:30
integer ( pInt ) :: i , j , m , c
2018-01-10 21:43:25 +05:30
2018-09-23 20:56:13 +05:30
allocate ( mesh_node0 ( 3 , mesh_Nnodes ) , source = 0.0_pReal )
allocate ( mesh_node ( 3 , mesh_Nnodes ) , source = 0.0_pReal )
2009-10-12 21:31:49 +05:30
2012-06-15 21:40:21 +05:30
inPart = . false .
2019-02-02 14:26:11 +05:30
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) == 'part' ) inPart = . false .
2012-06-15 21:40:21 +05:30
if ( ( inPart . or . noPart ) . and . &
2015-08-28 13:08:48 +05:30
IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*node' . and . &
( IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'print' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'file' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'response' ) &
2012-06-15 21:40:21 +05:30
) then
2013-12-11 22:19:20 +05:30
c = IO_countDataLines ( fileUnit ) ! how many nodes are defined here?
2012-06-15 21:40:21 +05:30
do i = 1_pInt , c
2013-12-11 22:19:20 +05:30
backspace ( fileUnit ) ! rewind to first entry
2012-06-15 21:40:21 +05:30
enddo
do i = 1_pInt , c
2019-02-02 21:49:12 +05:30
read ( fileUnit , '(a300)' ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
m = mesh_FEasCP ( 'node' , IO_intValue ( line , chunkPos , 1_pInt ) )
2013-02-09 13:53:47 +05:30
do j = 1_pInt , 3_pInt
2015-08-28 13:08:48 +05:30
mesh_node0 ( j , m ) = mesh_unitlength * IO_floatValue ( line , chunkPos , j + 1_pInt )
2018-01-10 21:43:25 +05:30
enddo
2012-06-15 21:40:21 +05:30
enddo
endif
enddo
2019-02-02 14:26:11 +05:30
if ( int ( size ( mesh_node0 , 2_pInt ) , pInt ) / = mesh_Nnodes ) call IO_error ( error_ID = 909_pInt )
2012-06-15 21:40:21 +05:30
mesh_node = mesh_node0
end subroutine mesh_abaqus_build_nodes
!--------------------------------------------------------------------------------------------------
!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements.
2018-01-10 21:43:25 +05:30
!! Sets global values 'mesh_maxNnodes', 'mesh_maxNips', 'mesh_maxNipNeighbors',
2014-12-15 17:21:32 +05:30
!! and 'mesh_maxNcellnodes'
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_count_cpSizes ( fileUnit )
2012-03-09 01:55:28 +05:30
use IO , only : IO_lc , &
IO_stringValue , &
IO_stringPos , &
IO_extractValue , &
IO_error , &
IO_countDataLines , &
IO_intValue
2009-10-12 21:31:49 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2009-10-12 21:31:49 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2012-03-09 01:55:28 +05:30
character ( len = 300 ) :: line
2019-02-02 14:26:11 +05:30
integer :: myStat
2012-03-09 01:55:28 +05:30
logical :: inPart
2019-02-02 14:26:11 +05:30
integer ( pInt ) :: i , c , t , g
2009-10-12 21:31:49 +05:30
mesh_maxNnodes = 0_pInt
mesh_maxNips = 0_pInt
mesh_maxNipNeighbors = 0_pInt
2013-04-15 13:43:20 +05:30
mesh_maxNcellnodes = 0_pInt
2009-10-12 21:31:49 +05:30
inPart = . false .
2019-02-02 14:26:11 +05:30
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 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 . &
2015-08-28 13:08:48 +05:30
IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*element' . and . &
( IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'matrix' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'response' ) &
2009-10-12 21:31:49 +05:30
) then
2015-08-28 13:08:48 +05:30
t = FE_mapElemtype ( IO_extractValue ( IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) , 'type' ) ) ! remember elem type
2012-11-16 04:15:20 +05:30
g = FE_geomtype ( t )
2013-04-15 13:43:20 +05:30
c = FE_celltype ( g )
mesh_maxNnodes = max ( mesh_maxNnodes , FE_Nnodes ( t ) )
2012-11-16 04:15:20 +05:30
mesh_maxNips = max ( mesh_maxNips , FE_Nips ( g ) )
2013-04-15 13:43:20 +05:30
mesh_maxNipNeighbors = max ( mesh_maxNipNeighbors , FE_NipNeighbors ( c ) )
mesh_maxNcellnodes = max ( mesh_maxNcellnodes , FE_Ncellnodes ( g ) )
2009-10-12 21:31:49 +05:30
endif
enddo
2018-01-10 21:43:25 +05:30
2019-02-02 14:26:11 +05:30
end subroutine mesh_abaqus_count_cpSizes
2009-10-12 21:31:49 +05:30
2012-06-15 21:40:21 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Store FEid, type, mat, tex, and node list per elemen.
!! Allocates global array 'mesh_element'
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_abaqus_build_elements ( fileUnit )
2012-03-09 01:55:28 +05:30
use IO , only : IO_lc , &
IO_stringValue , &
IO_stringPos , &
2012-06-15 21:40:21 +05:30
IO_intValue , &
IO_extractValue , &
2012-03-09 01:55:28 +05:30
IO_floatValue , &
2018-09-23 18:49:47 +05:30
IO_countDataLines , &
IO_error
2009-10-12 21:31:49 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2009-10-12 21:31:49 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2019-02-02 14:26:11 +05:30
character ( len = 300 ) :: line
integer :: myStat
2019-02-02 21:49:12 +05:30
logical :: inPart , materialFound
2013-04-18 22:10:49 +05:30
integer ( pInt ) :: i , j , k , c , e , t , homog , micro , nNodesAlreadyRead
2012-06-15 21:40:21 +05:30
character ( len = 64 ) :: materialName , elemSetName
2018-09-23 23:28:43 +05:30
allocate ( mesh_element ( 4_pInt + mesh_maxNnodes , mesh_NcpElems ) , source = 0_pInt )
2018-09-23 18:49:47 +05:30
mesh_elemType = - 1_pInt
2009-10-12 21:31:49 +05:30
inPart = . false .
2019-02-02 14:26:11 +05:30
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-31 22:00:04 +05:30
chunkPos = IO_stringPos ( line )
2015-08-28 13:08:48 +05:30
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*part' ) inPart = . true .
if ( IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*end' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 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 . &
2015-08-28 13:08:48 +05:30
IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) == '*element' . and . &
( IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'output' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'matrix' . and . &
IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) / = 'response' ) &
2012-06-15 21:40:21 +05:30
) then
2018-09-23 18:49:47 +05:30
t = FE_mapElemtype ( IO_extractValue ( IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) ) , 'type' ) ) ! remember elem type
2013-12-11 22:19:20 +05:30
c = IO_countDataLines ( fileUnit )
2012-02-16 00:28:38 +05:30
do i = 1_pInt , c
2013-12-11 22:19:20 +05:30
backspace ( fileUnit )
2009-10-12 21:31:49 +05:30
enddo
2012-02-16 00:28:38 +05:30
do i = 1_pInt , c
2019-02-02 21:49:12 +05:30
read ( fileUnit , '(a300)' ) line
2019-02-02 13:59:58 +05:30
chunkPos = IO_stringPos ( line ) ! limit to 64 nodes max
2015-08-28 13:08:48 +05:30
e = mesh_FEasCP ( 'elem' , IO_intValue ( line , chunkPos , 1_pInt ) )
2018-09-23 18:49:47 +05:30
if ( e / = 0_pInt ) then ! disregard non CP elems
2018-09-23 22:10:14 +05:30
mesh_element ( 1 , e ) = - 1_pInt ! DEPRECATED
2018-09-23 18:49:47 +05:30
if ( mesh_elemType / = t . and . mesh_elemType / = - 1_pInt ) &
call IO_error ( 191 , el = t , ip = mesh_elemType )
mesh_elemType = t
2012-11-16 04:15:20 +05:30
mesh_element ( 2 , e ) = t ! elem type
2013-04-15 13:43:20 +05:30
nNodesAlreadyRead = 0_pInt
2015-08-28 13:08:48 +05:30
do j = 1_pInt , chunkPos ( 1 ) - 1_pInt
mesh_element ( 4_pInt + j , e ) = mesh_FEasCP ( 'node' , IO_intValue ( line , chunkPos , 1_pInt + j ) ) ! put CP ids of nodes to position 5:
2013-02-07 16:15:10 +05:30
enddo
2015-08-28 13:08:48 +05:30
nNodesAlreadyRead = chunkPos ( 1 ) - 1_pInt
2013-04-15 13:43:20 +05:30
do while ( nNodesAlreadyRead < FE_Nnodes ( t ) ) ! read on if not all nodes in one line
2019-02-02 21:49:12 +05:30
read ( fileUnit , '(a300)' ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
do j = 1_pInt , chunkPos ( 1 )
2013-04-22 00:18:59 +05:30
mesh_element ( 4_pInt + nNodesAlreadyRead + j , e ) &
2015-08-28 13:08:48 +05:30
= mesh_FEasCP ( 'node' , IO_IntValue ( line , chunkPos , j ) ) ! CP ids of nodes
2013-04-15 13:43:20 +05:30
enddo
2015-08-28 13:08:48 +05:30
nNodesAlreadyRead = nNodesAlreadyRead + chunkPos ( 1 )
2013-04-15 13:43:20 +05:30
enddo
2012-06-15 21:40:21 +05:30
endif
2009-10-12 21:31:49 +05:30
enddo
endif
enddo
2018-01-10 21:43:25 +05:30
2019-02-02 14:26:11 +05:30
rewind ( fileUnit ) ! just in case "*material" definitions apear before "*element"
2009-10-12 21:31:49 +05:30
materialFound = . false .
2019-02-02 14:26:11 +05:30
myStat = 0
rewind ( fileUnit )
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
select case ( IO_lc ( IO_StringValue ( line , chunkPos , 1_pInt ) ) )
2009-10-12 21:31:49 +05:30
case ( '*material' )
2015-08-28 13:08:48 +05:30
materialName = trim ( IO_extractValue ( IO_lc ( IO_StringValue ( line , chunkPos , 2_pInt ) ) , 'name' ) ) ! extract name=value
2012-03-09 01:55:28 +05:30
materialFound = materialName / = '' ! valid name?
2009-10-12 21:31:49 +05:30
case ( '*user' )
2015-08-28 13:08:48 +05:30
if ( IO_lc ( IO_StringValue ( line , chunkPos , 2_pInt ) ) == 'material' . and . &
2009-10-12 21:31:49 +05:30
materialFound ) then
2019-02-02 21:49:12 +05:30
read ( fileUnit , '(a300)' ) line ! read homogenization and microstructure
2015-08-31 22:00:04 +05:30
chunkPos = IO_stringPos ( line )
2015-08-28 13:08:48 +05:30
homog = nint ( IO_floatValue ( line , chunkPos , 1_pInt ) , pInt )
micro = nint ( IO_floatValue ( line , chunkPos , 2_pInt ) , pInt )
2012-03-09 01:55:28 +05:30
do i = 1_pInt , mesh_Nmaterials ! look thru material names
if ( materialName == mesh_nameMaterial ( i ) ) then ! found one
elemSetName = mesh_mapMaterial ( i ) ! take corresponding elemSet
do k = 1_pInt , mesh_NelemSets ! look thru all elemSet definitions
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 ) )
2012-03-09 01:55:28 +05:30
mesh_element ( 3 , e ) = homog ! store homogenization
mesh_element ( 4 , e ) = micro ! store microstructure
2011-03-23 21:50:12 +05:30
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
2019-02-02 14:26:11 +05:30
end subroutine mesh_abaqus_build_elements
2019-01-28 18:29:54 +05:30
2009-10-12 21:31:49 +05:30
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief get any additional damask options from input file, sets mesh_periodicSurface
!--------------------------------------------------------------------------------------------------
2013-12-11 22:19:20 +05:30
subroutine mesh_get_damaskOptions ( fileUnit )
2012-06-15 21:40:21 +05:30
use IO , only : &
IO_lc , &
IO_stringValue , &
IO_stringPos
2011-02-16 21:53:08 +05:30
2012-06-19 20:03:24 +05:30
implicit none
2013-12-11 22:19:20 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2012-06-15 21:40:21 +05:30
2015-08-28 13:08:48 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
2019-02-02 14:26:11 +05:30
character ( len = 300 ) :: line
integer :: myStat
logical :: inPart
2012-06-19 20:03:24 +05:30
integer ( pInt ) chunk , Nchunks
2019-02-02 14:26:11 +05:30
character ( len = 300 ) :: damaskOption , v
character ( len = * ) , parameter :: keyword = '**damask'
2012-06-15 21:40:21 +05:30
2013-04-22 14:02:51 +05:30
mesh_periodicSurface = . false .
2019-02-02 14:26:11 +05:30
myStat = 0
2013-12-11 22:19:20 +05:30
rewind ( fileUnit )
2019-02-02 14:26:11 +05:30
do while ( myStat == 0 )
read ( fileUnit , '(a300)' , iostat = myStat ) line
2015-08-28 13:08:48 +05:30
chunkPos = IO_stringPos ( line )
Nchunks = chunkPos ( 1 )
if ( IO_lc ( IO_stringValue ( line , chunkPos , 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 , chunkPos , 2_pInt ) )
2012-06-19 20:03:24 +05:30
select case ( damaskOption )
2013-04-22 14:02:51 +05:30
case ( 'periodic' ) ! damask Option that allows to specify periodic fluxes
do chunk = 3_pInt , Nchunks ! loop through chunks (skipping the keyword)
2015-08-28 13:08:48 +05:30
v = IO_lc ( IO_stringValue ( line , chunkPos , chunk ) ) ! chunk matches keyvalues x,y, or z?
2012-06-19 20:03:24 +05:30
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'
enddo
endselect
endif
enddo
2012-06-15 21:40:21 +05:30
2019-02-02 14:26:11 +05:30
end subroutine mesh_get_damaskOptions
2019-01-28 18:29:54 +05:30
2012-06-15 21:40:21 +05:30
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal'
!--------------------------------------------------------------------------------------------------
2012-06-15 21:40:21 +05:30
subroutine mesh_build_ipAreas
2013-04-15 13:43:20 +05:30
use math , only : &
2016-01-10 19:04:26 +05:30
math_crossproduct
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
implicit none
2013-04-15 13:43:20 +05:30
integer ( pInt ) :: e , t , g , c , i , f , n , m
real ( pReal ) , dimension ( 3 , FE_maxNcellnodesPerCellface ) :: nodePos , normals
real ( pReal ) , dimension ( 3 ) :: normal
2012-06-15 21:40:21 +05:30
2013-12-28 01:33:28 +05:30
allocate ( mesh_ipArea ( mesh_maxNipNeighbors , mesh_maxNips , mesh_NcpElems ) , source = 0.0_pReal )
allocate ( mesh_ipAreaNormal ( 3_pInt , mesh_maxNipNeighbors , mesh_maxNips , mesh_NcpElems ) , source = 0.0_pReal )
2013-04-15 13:43:20 +05:30
!$OMP PARALLEL DO PRIVATE(t,g,c,nodePos,normal,normals)
do e = 1_pInt , mesh_NcpElems ! loop over cpElems
t = mesh_element ( 2_pInt , e ) ! get element type
g = FE_geomtype ( t ) ! get geometry type
c = FE_celltype ( g ) ! get cell type
select case ( c )
case ( 1_pInt , 2_pInt ) ! 2D 3 or 4 node
do i = 1_pInt , FE_Nips ( g ) ! loop over ips=cells in this element
do f = 1_pInt , FE_NipNeighbors ( c ) ! loop over cell faces
forall ( n = 1_pInt : FE_NcellnodesPerCellface ( c ) ) &
2013-04-22 00:18:59 +05:30
nodePos ( 1 : 3 , n ) = mesh_cellnode ( 1 : 3 , mesh_cell ( FE_cellface ( n , f , c ) , i , e ) )
2013-04-15 13:43:20 +05:30
normal ( 1 ) = nodePos ( 2 , 2 ) - nodePos ( 2 , 1 ) ! x_normal = y_connectingVector
normal ( 2 ) = - ( nodePos ( 1 , 2 ) - nodePos ( 1 , 1 ) ) ! y_normal = -x_connectingVector
normal ( 3 ) = 0.0_pReal
2016-01-10 19:04:26 +05:30
mesh_ipArea ( f , i , e ) = norm2 ( normal )
mesh_ipAreaNormal ( 1 : 3 , f , i , e ) = normal / norm2 ( normal ) ! ensure unit length of area normal
2013-04-15 13:43:20 +05:30
enddo
2013-04-09 23:37:30 +05:30
enddo
2013-04-15 13:43:20 +05:30
case ( 3_pInt ) ! 3D 4node
do i = 1_pInt , FE_Nips ( g ) ! loop over ips=cells in this element
do f = 1_pInt , FE_NipNeighbors ( c ) ! loop over cell faces
forall ( n = 1_pInt : FE_NcellnodesPerCellface ( c ) ) &
2013-04-22 00:18:59 +05:30
nodePos ( 1 : 3 , n ) = mesh_cellnode ( 1 : 3 , mesh_cell ( FE_cellface ( n , f , c ) , i , e ) )
2016-01-10 19:04:26 +05:30
normal = math_crossproduct ( nodePos ( 1 : 3 , 2 ) - nodePos ( 1 : 3 , 1 ) , &
2013-04-15 13:43:20 +05:30
nodePos ( 1 : 3 , 3 ) - nodePos ( 1 : 3 , 1 ) )
2016-01-10 19:04:26 +05:30
mesh_ipArea ( f , i , e ) = norm2 ( normal )
mesh_ipAreaNormal ( 1 : 3 , f , i , e ) = normal / norm2 ( normal ) ! ensure unit length of area normal
2013-04-15 13:43:20 +05:30
enddo
2013-04-09 23:37:30 +05:30
enddo
2013-02-05 18:57:37 +05:30
2013-04-15 13:43:20 +05:30
case ( 4_pInt ) ! 3D 8node
2018-01-10 21:43:25 +05:30
! for this cell type we get the normal of the quadrilateral face as an average of
2013-04-15 13:43:20 +05:30
! four normals of triangular subfaces; since the face consists only of two triangles,
2018-01-10 21:43:25 +05:30
! the sum has to be divided by two; this whole prcedure tries to compensate for
! probable non-planar cell surfaces
2013-04-15 13:43:20 +05:30
m = FE_NcellnodesPerCellface ( c )
do i = 1_pInt , FE_Nips ( g ) ! loop over ips=cells in this element
do f = 1_pInt , FE_NipNeighbors ( c ) ! loop over cell faces
forall ( n = 1_pInt : FE_NcellnodesPerCellface ( c ) ) &
2013-04-22 00:18:59 +05:30
nodePos ( 1 : 3 , n ) = mesh_cellnode ( 1 : 3 , mesh_cell ( FE_cellface ( n , f , c ) , i , e ) )
2013-04-15 13:43:20 +05:30
forall ( n = 1_pInt : FE_NcellnodesPerCellface ( c ) ) &
normals ( 1 : 3 , n ) = 0.5_pReal &
2016-01-10 19:04:26 +05:30
* math_crossproduct ( nodePos ( 1 : 3 , 1 + mod ( n , m ) ) - nodePos ( 1 : 3 , n ) , &
2013-04-15 13:43:20 +05:30
nodePos ( 1 : 3 , 1 + mod ( n + 1 , m ) ) - nodePos ( 1 : 3 , n ) )
normal = 0.5_pReal * sum ( normals , 2 )
2016-01-10 19:04:26 +05:30
mesh_ipArea ( f , i , e ) = norm2 ( normal )
mesh_ipAreaNormal ( 1 : 3 , f , i , e ) = normal / norm2 ( normal )
2013-04-15 13:43:20 +05:30
enddo
enddo
end select
enddo
!$OMP END PARALLEL DO
2018-01-10 21:43:25 +05:30
2015-04-11 17:17:33 +05:30
end subroutine mesh_build_ipAreas
2018-01-10 21:43:25 +05:30
2019-01-28 18:29:54 +05:30
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief assignment of twin nodes for each cp node, allocate globals '_nodeTwins'
!--------------------------------------------------------------------------------------------------
2012-06-15 21:40:21 +05:30
subroutine mesh_build_nodeTwins
2013-04-19 18:11:06 +05:30
implicit none
integer ( pInt ) dir , & ! direction of periodicity
node , &
minimumNode , &
maximumNode , &
n1 , &
n2
2015-03-25 21:36:19 +05:30
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
2013-04-19 18:11:06 +05:30
logical , dimension ( mesh_Nnodes ) :: unpaired
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
allocate ( mesh_nodeTwins ( 3 , mesh_Nnodes ) )
mesh_nodeTwins = 0_pInt
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
tolerance = 0.001_pReal * minval ( mesh_ipVolume ) ** 0.333_pReal
2018-01-10 21:43:25 +05:30
2015-03-25 21:36:19 +05:30
do dir = 1_pInt , 3_pInt ! check periodicity in directions of x,y,z
if ( mesh_periodicSurface ( dir ) ) then ! only if periodicity is requested
2018-01-10 21:43:25 +05:30
!*** find out which nodes sit on the surface
2013-04-19 18:11:06 +05:30
!*** and have a minimum or maximum position in this dimension
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
minimumNodes = 0_pInt
maximumNodes = 0_pInt
minCoord = minval ( mesh_node0 ( dir , : ) )
maxCoord = maxval ( mesh_node0 ( dir , : ) )
2015-03-25 21:36:19 +05:30
do node = 1_pInt , mesh_Nnodes ! loop through all nodes and find surface nodes
2013-04-19 18:11:06 +05:30
if ( abs ( mesh_node0 ( dir , node ) - minCoord ) < = tolerance ) then
minimumNodes ( 1 ) = minimumNodes ( 1 ) + 1_pInt
minimumNodes ( minimumNodes ( 1 ) + 1_pInt ) = node
elseif ( abs ( mesh_node0 ( dir , node ) - maxCoord ) < = tolerance ) then
maximumNodes ( 1 ) = maximumNodes ( 1 ) + 1_pInt
maximumNodes ( maximumNodes ( 1 ) + 1_pInt ) = node
endif
enddo
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
!*** find the corresponding node on the other side with the same position in this dimension
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
unpaired = . true .
do n1 = 1_pInt , minimumNodes ( 1 )
minimumNode = minimumNodes ( n1 + 1_pInt )
if ( unpaired ( minimumNode ) ) then
do n2 = 1_pInt , maximumNodes ( 1 )
maximumNode = maximumNodes ( n2 + 1_pInt )
distance = abs ( mesh_node0 ( : , minimumNode ) - mesh_node0 ( : , maximumNode ) )
2015-03-25 21:36:19 +05:30
if ( sum ( distance ) - distance ( dir ) < = tolerance ) then ! minimum possible distance (within tolerance)
2013-04-19 18:11:06 +05:30
mesh_nodeTwins ( dir , minimumNode ) = maximumNode
mesh_nodeTwins ( dir , maximumNode ) = minimumNode
2015-03-25 21:36:19 +05:30
unpaired ( maximumNode ) = . false . ! remember this node, we don't have to look for his partner again
2013-04-19 18:11:06 +05:30
exit
endif
enddo
endif
enddo
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
endif
enddo
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
end subroutine mesh_build_nodeTwins
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
2018-01-10 21:43:25 +05:30
!> @brief get maximum count of shared elements among cpElements and build list of elements shared
2013-04-18 22:10:49 +05:30
!! by each node in mesh. Allocate globals '_maxNsharedElems' and '_sharedElem'
!--------------------------------------------------------------------------------------------------
2012-06-15 21:40:21 +05:30
subroutine mesh_build_sharedElems
2013-04-19 18:11:06 +05:30
implicit none
integer ( pint ) e , & ! element index
g , & ! element type
node , & ! CP node index
2018-01-10 21:43:25 +05:30
n , & ! node index per element
myDim , & ! dimension index
2013-04-19 18:11:06 +05:30
nodeTwin ! node twin in the specified dimension
integer ( pInt ) , dimension ( mesh_Nnodes ) :: node_count
2018-09-23 20:56:13 +05:30
integer ( pInt ) , dimension ( : ) , allocatable :: node_seen
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
allocate ( node_seen ( maxval ( FE_NmatchingNodes ) ) )
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
node_count = 0_pInt
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
do e = 1_pInt , mesh_NcpElems
g = FE_geomtype ( mesh_element ( 2 , e ) ) ! get elemGeomType
node_seen = 0_pInt ! reset node duplicates
2013-04-22 00:18:59 +05:30
do n = 1_pInt , FE_NmatchingNodes ( g ) ! check each node of element
node = mesh_element ( 4 + n , e )
2013-04-19 18:11:06 +05:30
if ( all ( node_seen / = node ) ) then
node_count ( node ) = node_count ( node ) + 1_pInt ! if FE node not yet encountered -> count it
2013-04-22 00:18:59 +05:30
do myDim = 1_pInt , 3_pInt ! check in each dimension...
2013-04-19 18:11:06 +05:30
nodeTwin = mesh_nodeTwins ( myDim , node )
if ( nodeTwin > 0_pInt ) & ! if I am a twin of some node...
node_count ( nodeTwin ) = node_count ( nodeTwin ) + 1_pInt ! -> count me again for the twin node
enddo
endif
node_seen ( n ) = node ! remember this node to be counted already
enddo
enddo
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
mesh_maxNsharedElems = int ( maxval ( node_count ) , pInt ) ! most shared node
2018-01-10 21:43:25 +05:30
2018-09-23 20:56:13 +05:30
allocate ( mesh_sharedElem ( 1 + mesh_maxNsharedElems , mesh_Nnodes ) , source = 0_pInt )
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
do e = 1_pInt , mesh_NcpElems
g = FE_geomtype ( mesh_element ( 2 , e ) ) ! get elemGeomType
node_seen = 0_pInt
do n = 1_pInt , FE_NmatchingNodes ( g )
2013-04-22 00:18:59 +05:30
node = mesh_element ( 4_pInt + n , e )
2013-04-19 18:11:06 +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
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...
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 ( n ) = node
enddo
enddo
2018-01-10 21:43:25 +05:30
2012-03-09 01:55:28 +05:30
end subroutine mesh_build_sharedElems
2009-10-12 21:31:49 +05:30
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
2013-12-28 01:33:28 +05:30
!> @brief build up of IP neighborhood, allocate globals '_ipNeighborhood'
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
subroutine mesh_build_ipNeighborhood
2013-04-19 18:11:06 +05:30
use math , only : &
math_mul3x3
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
implicit none
integer ( pInt ) :: myElem , & ! my CP element index
myIP , &
myType , & ! my element type
myFace , &
neighbor , & ! neighor index
2018-01-10 21:43:25 +05:30
neighboringIPkey , & ! positive integer indicating the neighboring IP (for intra-element) and negative integer indicating the face towards neighbor (for neighboring element)
2013-04-19 18:11:06 +05:30
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
a , anchor , &
2018-01-10 21:43:25 +05:30
neighboringIP , &
2013-04-19 18:11:06 +05:30
neighboringElem , &
pointingToMe
integer ( pInt ) , dimension ( FE_maxmaxNnodesAtIP ) :: &
linkedNodes = 0_pInt , &
matchingNodes
logical checkTwins
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
allocate ( mesh_ipNeighborhood ( 3 , mesh_maxNipNeighbors , mesh_maxNips , mesh_NcpElems ) )
mesh_ipNeighborhood = 0_pInt
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
do myElem = 1_pInt , mesh_NcpElems ! loop over cpElems
myType = FE_geomtype ( mesh_element ( 2 , myElem ) ) ! get elemGeomType
do myIP = 1_pInt , FE_Nips ( myType ) ! loop over IPs of elem
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
do neighbor = 1_pInt , FE_NipNeighbors ( FE_celltype ( myType ) ) ! loop over neighbors of IP
neighboringIPkey = FE_ipNeighbor ( neighbor , myIP , myType )
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
!*** if the key is positive, the neighbor is inside the element
!*** that means, we have already found our neighboring IP
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
if ( neighboringIPkey > 0_pInt ) then
mesh_ipNeighborhood ( 1 , neighbor , myIP , myElem ) = myElem
mesh_ipNeighborhood ( 2 , neighbor , myIP , myElem ) = neighboringIPkey
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
!*** 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
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
elseif ( neighboringIPkey < 0_pInt ) then ! neighboring element's IP
myFace = - neighboringIPkey
call mesh_faceMatch ( myElem , myFace , matchingElem , matchingFace ) ! get face and CP elem id of face match
if ( matchingElem > 0_pInt ) then ! found match?
neighboringType = FE_geomtype ( mesh_element ( 2 , matchingElem ) )
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
!*** trivial solution if neighbor has only one IP
2018-01-10 21:43:25 +05:30
if ( FE_Nips ( neighboringType ) == 1_pInt ) then
2013-04-19 18:11:06 +05:30
mesh_ipNeighborhood ( 1 , neighbor , myIP , myElem ) = matchingElem
mesh_ipNeighborhood ( 2 , neighbor , myIP , myElem ) = 1_pInt
cycle
endif
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
!*** find those nodes which build the link to the neighbor
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
NlinkedNodes = 0_pInt
linkedNodes = 0_pInt
do a = 1_pInt , FE_maxNnodesAtIP ( myType ) ! figure my anchor nodes on connecting face
anchor = FE_nodesAtIP ( a , myIP , myType )
if ( anchor / = 0_pInt ) then ! valid anchor node
if ( any ( FE_face ( : , myFace , myType ) == anchor ) ) then ! ip anchor sits on face?
NlinkedNodes = NlinkedNodes + 1_pInt
2013-04-22 00:18:59 +05:30
linkedNodes ( NlinkedNodes ) = mesh_element ( 4_pInt + anchor , myElem ) ! CP id of anchor node
2013-04-19 18:11:06 +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
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
!*** loop through the ips of my neighbor
!*** and try to find an ip with matching nodes
!*** also try to match with node twins
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
checkCandidateIP : do candidateIP = 1_pInt , FE_Nips ( neighboringType )
NmatchingNodes = 0_pInt
matchingNodes = 0_pInt
do a = 1_pInt , FE_maxNnodesAtIP ( neighboringType ) ! check each anchor node of that ip
anchor = FE_nodesAtIP ( a , candidateIP , neighboringType )
if ( anchor / = 0_pInt ) then ! valid anchor node
if ( any ( FE_face ( : , matchingFace , neighboringType ) == anchor ) ) then ! sits on matching face?
NmatchingNodes = NmatchingNodes + 1_pInt
2013-04-22 00:18:59 +05:30
matchingNodes ( NmatchingNodes ) = mesh_element ( 4 + anchor , matchingElem ) ! CP id of neighbor's anchor node
2013-04-19 18:11:06 +05:30
else ! no matching, because not all nodes sit on the matching face
NmatchingNodes = 0_pInt
matchingNodes = 0_pInt
exit
endif
endif
enddo
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
if ( NmatchingNodes / = NlinkedNodes ) & ! this ip has wrong count of anchors on face
cycle checkCandidateIP
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
!*** check "normal" nodes whether they match or not
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
checkTwins = . false .
do a = 1_pInt , NlinkedNodes
if ( all ( matchingNodes / = linkedNodes ( a ) ) ) then ! this linkedNode does not match any matchingNode
checkTwins = . true .
exit ! no need to search further
endif
enddo
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
!*** if no match found, then also check node twins
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
if ( checkTwins ) then
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
twin_of_linkedNode = mesh_nodeTwins ( dir , linkedNodes ( a ) )
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
endif
enddo
endif
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
!*** we found a match !!!
2018-01-10 21:43:25 +05:30
2013-04-19 18:11:06 +05:30
mesh_ipNeighborhood ( 1 , neighbor , myIP , myElem ) = matchingElem
mesh_ipNeighborhood ( 2 , neighbor , myIP , myElem ) = candidateIP
2018-01-10 21:43:25 +05:30
exit checkCandidateIP
2013-04-19 18:11:06 +05:30
enddo checkCandidateIP
endif ! end of valid external matching
endif ! end of internal/external matching
enddo
enddo
enddo
do myElem = 1_pInt , mesh_NcpElems ! loop over cpElems
myType = FE_geomtype ( mesh_element ( 2 , myElem ) ) ! get elemGeomType
do myIP = 1_pInt , FE_Nips ( myType ) ! loop over IPs of elem
do neighbor = 1_pInt , FE_NipNeighbors ( FE_celltype ( myType ) ) ! loop over neighbors of IP
neighboringElem = mesh_ipNeighborhood ( 1 , neighbor , myIP , myElem )
neighboringIP = mesh_ipNeighborhood ( 2 , neighbor , myIP , myElem )
if ( neighboringElem > 0_pInt . and . neighboringIP > 0_pInt ) then ! if neighbor exists ...
neighboringType = FE_geomtype ( mesh_element ( 2 , neighboringElem ) )
do pointingToMe = 1_pInt , FE_NipNeighbors ( FE_celltype ( neighboringType ) ) ! find neighboring index that points from my neighbor to myself
if ( myElem == mesh_ipNeighborhood ( 1 , pointingToMe , neighboringIP , neighboringElem ) &
. and . myIP == mesh_ipNeighborhood ( 2 , pointingToMe , neighboringIP , neighboringElem ) ) then ! possible candidate
if ( math_mul3x3 ( mesh_ipAreaNormal ( 1 : 3 , neighbor , myIP , myElem ) , &
mesh_ipAreaNormal ( 1 : 3 , pointingToMe , neighboringIP , neighboringElem ) ) < 0.0_pReal ) then ! area normals have opposite orientation (we have to check that because of special case for single element with two ips and periodicity. In this case the neighbor is identical in two different directions.)
mesh_ipNeighborhood ( 3 , neighbor , myIP , myElem ) = pointingToMe ! found match
exit ! so no need to search further
endif
endif
enddo
endif
enddo
enddo
enddo
2018-01-10 21:43:25 +05:30
2012-03-09 01:55:28 +05:30
end subroutine mesh_build_ipNeighborhood
2019-01-28 18:29:54 +05:30
2009-01-20 00:12:31 +05:30
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief mapping of FE element types to internal representation
!--------------------------------------------------------------------------------------------------
2012-06-15 21:40:21 +05:30
integer ( pInt ) function FE_mapElemtype ( what )
2013-04-10 15:08:40 +05:30
use IO , only : IO_lc , IO_error
2012-04-24 22:32:27 +05:30
2012-06-15 21:40:21 +05:30
implicit none
character ( len = * ) , intent ( in ) :: what
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +05:30
select case ( IO_lc ( what ) )
2019-01-31 15:59:56 +05:30
case ( 'cpe4' , &
2014-10-23 13:47:37 +05:30
'cpe4t' )
2012-06-15 21:40:21 +05:30
FE_mapElemtype = 3_pInt ! Arbitrary Quadrilateral Plane-strain
2019-01-31 15:59:56 +05:30
case ( 'cpe8' , &
2014-10-23 13:47:37 +05:30
'cpe8t' )
2012-11-16 04:15:20 +05:30
FE_mapElemtype = 4_pInt ! Plane Strain, Eight-node Distorted Quadrilateral
2019-01-31 15:59:56 +05:30
case ( 'c3d4' , &
2014-10-23 13:47:37 +05:30
'c3d4t' )
2013-04-10 18:52:02 +05:30
FE_mapElemtype = 6_pInt ! Three-dimensional Four-node Tetrahedron
2019-01-31 15:59:56 +05:30
case ( 'c3d6' , &
2014-10-23 13:47:37 +05:30
'c3d6t' )
2013-04-10 18:52:02 +05:30
FE_mapElemtype = 9_pInt ! Three-dimensional Arbitrarily Distorted Pentahedral
2019-01-31 15:59:56 +05:30
case ( 'c3d8r' , &
2014-10-23 13:47:37 +05:30
'c3d8rt' )
2013-04-10 18:52:02 +05:30
FE_mapElemtype = 10_pInt ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration
2019-01-31 15:59:56 +05:30
case ( 'c3d8' , &
2014-10-23 13:47:37 +05:30
'c3d8t' )
2013-04-10 18:52:02 +05:30
FE_mapElemtype = 11_pInt ! Three-dimensional Arbitrarily Distorted Brick
2019-01-31 15:59:56 +05:30
case ( 'c3d20r' , &
2014-10-23 13:47:37 +05:30
'c3d20rt' )
2013-04-10 18:52:02 +05:30
FE_mapElemtype = 12_pInt ! Three-dimensional Arbitrarily Distorted quad hexahedral with reduced integration
2019-01-31 15:59:56 +05:30
case ( 'c3d20' , &
2014-10-23 13:47:37 +05:30
'c3d20t' )
2013-04-10 18:52:02 +05:30
FE_mapElemtype = 13_pInt ! Three-dimensional Arbitrarily Distorted quadratic hexahedral
2018-01-10 21:43:25 +05:30
case default
2013-04-10 15:08:40 +05:30
call IO_error ( error_ID = 190_pInt , ext_msg = IO_lc ( what ) )
2012-11-16 04:15:20 +05:30
end select
2012-04-24 22:32:27 +05:30
2012-06-15 21:40:21 +05:30
end function FE_mapElemtype
2012-05-08 20:27:06 +05:30
2013-04-18 22:10:49 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief find face-matching element of same type
!--------------------------------------------------------------------------------------------------
2012-06-15 21:40:21 +05:30
subroutine mesh_faceMatch ( elem , face , matchingElem , matchingFace )
2012-05-08 20:27:06 +05:30
2012-06-15 21:40:21 +05:30
implicit none
integer ( pInt ) , intent ( out ) :: matchingElem , & ! matching CP element ID
2018-01-10 21:43:25 +05:30
matchingFace ! matching face ID
2013-04-22 00:18:59 +05:30
integer ( pInt ) , intent ( in ) :: face , & ! face ID
elem ! CP elem ID
2013-04-15 13:43:20 +05:30
integer ( pInt ) , dimension ( FE_NmatchingNodesPerFace ( face , FE_geomtype ( mesh_element ( 2 , elem ) ) ) ) :: &
2012-06-15 21:40:21 +05:30
myFaceNodes ! global node ids on my face
integer ( pInt ) :: myType , &
candidateType , &
candidateElem , &
candidateFace , &
candidateFaceNode , &
minNsharedElems , &
NsharedElems , &
lonelyNode = 0_pInt , &
i , &
n , &
dir ! periodicity direction
integer ( pInt ) , dimension ( : ) , allocatable :: element_seen
logical checkTwins
2012-04-24 22:32:27 +05:30
2012-06-15 21:40:21 +05:30
matchingElem = 0_pInt
matchingFace = 0_pInt
minNsharedElems = mesh_maxNsharedElems + 1_pInt ! init to worst case
2012-11-16 04:15:20 +05:30
myType = FE_geomtype ( mesh_element ( 2_pInt , elem ) ) ! figure elemGeomType
2012-05-08 20:27:06 +05:30
2013-04-15 13:43:20 +05:30
do n = 1_pInt , FE_NmatchingNodesPerFace ( face , myType ) ! loop over nodes on face
2013-04-22 00:18:59 +05:30
myFaceNodes ( n ) = mesh_element ( 4_pInt + FE_face ( n , face , myType ) , elem ) ! CP id of face node
2012-06-15 21:40:21 +05:30
NsharedElems = mesh_sharedElem ( 1_pInt , myFaceNodes ( n ) ) ! figure # shared elements for this node
if ( NsharedElems < minNsharedElems ) then
minNsharedElems = NsharedElems ! remember min # shared elems
lonelyNode = n ! remember most lonely node
endif
enddo
2012-05-08 20:27:06 +05:30
2012-06-15 21:40:21 +05:30
allocate ( element_seen ( minNsharedElems ) )
element_seen = 0_pInt
2012-05-08 20:27:06 +05:30
2012-06-15 21:40:21 +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
if ( all ( element_seen / = candidateElem ) ) then ! element seen for the first time?
element_seen ( i ) = candidateElem
2012-11-16 04:15:20 +05:30
candidateType = FE_geomtype ( mesh_element ( 2_pInt , candidateElem ) ) ! figure elemGeomType of candidate
2012-06-15 21:40:21 +05:30
checkCandidateFace : do candidateFace = 1_pInt , FE_maxNipNeighbors ! check each face of candidate
2013-04-15 13:43:20 +05:30
if ( FE_NmatchingNodesPerFace ( candidateFace , candidateType ) &
/ = FE_NmatchingNodesPerFace ( face , myType ) & ! incompatible face
2012-11-16 04:15:20 +05:30
. or . ( candidateElem == elem . and . candidateFace == face ) ) then ! this is my face
2012-06-15 21:40:21 +05:30
cycle checkCandidateFace
endif
checkTwins = . false .
2013-04-15 13:43:20 +05:30
do n = 1_pInt , FE_NmatchingNodesPerFace ( candidateFace , candidateType ) ! loop through nodes on face
2013-04-22 00:18:59 +05:30
candidateFaceNode = mesh_element ( 4_pInt + FE_face ( n , candidateFace , candidateType ) , candidateElem )
2013-04-15 13:43:20 +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
2012-06-15 21:40:21 +05:30
exit
endif
enddo
if ( checkTwins ) then
checkCandidateFaceTwins : do dir = 1_pInt , 3_pInt
2013-04-15 13:43:20 +05:30
do n = 1_pInt , FE_NmatchingNodesPerFace ( candidateFace , candidateType ) ! loop through nodes on face
2013-04-22 00:18:59 +05:30
candidateFaceNode = mesh_element ( 4 + FE_face ( n , candidateFace , candidateType ) , candidateElem )
2013-04-15 13:43:20 +05:30
if ( all ( myFaceNodes / = mesh_nodeTwins ( dir , candidateFaceNode ) ) ) then ! node twin does not match either
2012-06-15 21:40:21 +05:30
if ( dir == 3_pInt ) then
cycle checkCandidateFace
else
2013-04-15 13:43:20 +05:30
cycle checkCandidateFaceTwins ! try twins in next dimension
2012-06-15 21:40:21 +05:30
endif
endif
enddo
exit checkCandidateFaceTwins
enddo checkCandidateFaceTwins
endif
matchingFace = candidateFace
matchingElem = candidateElem
2013-04-15 13:43:20 +05:30
exit checkCandidate ! found my matching candidate
2012-06-15 21:40:21 +05:30
enddo checkCandidateFace
endif
enddo checkCandidate
2012-05-08 20:27:06 +05:30
2012-06-15 21:40:21 +05:30
end subroutine mesh_faceMatch
2012-05-08 20:27:06 +05:30
2012-04-10 20:45:46 +05:30
2013-04-15 13:43:20 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief get properties of different types of finite elements
!> @details assign globals: FE_nodesAtIP, FE_ipNeighbor, FE_cellnodeParentnodeWeights, FE_subNodeOnIPFace
!--------------------------------------------------------------------------------------------------
2012-06-15 21:40:21 +05:30
subroutine mesh_build_FEdata
2012-04-11 22:58:08 +05:30
implicit none
2012-11-16 04:15:20 +05:30
integer ( pInt ) :: me
2018-09-23 20:56:13 +05:30
allocate ( FE_nodesAtIP ( FE_maxmaxNnodesAtIP , FE_maxNips , FE_Ngeomtypes ) , source = 0_pInt )
allocate ( FE_ipNeighbor ( FE_maxNipNeighbors , FE_maxNips , FE_Ngeomtypes ) , source = 0_pInt )
allocate ( FE_cell ( FE_maxNcellnodesPerCell , FE_maxNips , FE_Ngeomtypes ) , source = 0_pInt )
allocate ( FE_cellnodeParentnodeWeights ( FE_maxNnodes , FE_maxNcellnodes , FE_Nelemtypes ) , source = 0.0_pReal )
allocate ( FE_cellface ( FE_maxNcellnodesPerCellface , FE_maxNcellfaces , FE_Ncelltypes ) , source = 0_pInt )
2013-04-15 13:43:20 +05:30
!*** fill FE_nodesAtIP with data ***
2012-11-16 04:15:20 +05:30
me = 0_pInt
me = me + 1_pInt
FE_nodesAtIP ( 1 : FE_maxNnodesAtIP ( me ) , 1 : FE_Nips ( me ) , me ) = & ! element 6 (2D 3node 1ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2012-11-16 04:15:20 +05:30
1 , 2 , 3 &
] , pInt ) , [ FE_maxNnodesAtIP ( me ) , FE_Nips ( me ) ] )
2012-04-20 17:28:41 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
FE_nodesAtIP ( 1 : FE_maxNnodesAtIP ( me ) , 1 : FE_Nips ( me ) , me ) = & ! element 125 (2D 6node 3ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2012-11-16 04:15:20 +05:30
1 , &
2 , &
3 &
] , pInt ) , [ FE_maxNnodesAtIP ( me ) , FE_Nips ( me ) ] )
2012-04-12 00:16:36 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
FE_nodesAtIP ( 1 : FE_maxNnodesAtIP ( me ) , 1 : FE_Nips ( me ) , me ) = & ! element 11 (2D 4node 4ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
1 , &
2 , &
4 , &
3 &
2012-11-16 04:15:20 +05:30
] , pInt ) , [ FE_maxNnodesAtIP ( me ) , FE_Nips ( me ) ] )
2012-04-11 22:58:08 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
FE_nodesAtIP ( 1 : FE_maxNnodesAtIP ( me ) , 1 : FE_Nips ( me ) , me ) = & ! element 27 (2D 8node 9ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
1 , 0 , &
1 , 2 , &
2 , 0 , &
1 , 4 , &
0 , 0 , &
2 , 3 , &
4 , 0 , &
3 , 4 , &
3 , 0 &
2012-11-16 04:15:20 +05:30
] , pInt ) , [ FE_maxNnodesAtIP ( me ) , FE_Nips ( me ) ] )
2012-04-11 22:58:08 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
FE_nodesAtIP ( 1 : FE_maxNnodesAtIP ( me ) , 1 : FE_Nips ( me ) , me ) = & ! element 134 (3D 4node 1ip)
reshape ( int ( [ &
1 , 2 , 3 , 4 &
] , pInt ) , [ FE_maxNnodesAtIP ( me ) , FE_Nips ( me ) ] )
me = me + 1_pInt
FE_nodesAtIP ( 1 : FE_maxNnodesAtIP ( me ) , 1 : FE_Nips ( me ) , me ) = & ! element 127 (3D 10node 4ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
1 , &
2 , &
3 , &
4 &
2012-11-16 04:15:20 +05:30
] , pInt ) , [ FE_maxNnodesAtIP ( me ) , FE_Nips ( me ) ] )
2012-04-11 22:58:08 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
FE_nodesAtIP ( 1 : FE_maxNnodesAtIP ( me ) , 1 : FE_Nips ( me ) , me ) = & ! element 136 (3D 6node 6ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
1 , &
2 , &
3 , &
4 , &
5 , &
6 &
2012-11-16 04:15:20 +05:30
] , pInt ) , [ FE_maxNnodesAtIP ( me ) , FE_Nips ( me ) ] )
2012-04-20 17:28:41 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
FE_nodesAtIP ( 1 : FE_maxNnodesAtIP ( me ) , 1 : FE_Nips ( me ) , me ) = & ! element 117 (3D 8node 1ip)
reshape ( int ( [ &
1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 &
] , pInt ) , [ FE_maxNnodesAtIP ( me ) , FE_Nips ( me ) ] )
me = me + 1_pInt
FE_nodesAtIP ( 1 : FE_maxNnodesAtIP ( me ) , 1 : FE_Nips ( me ) , me ) = & ! element 7 (3D 8node 8ip)
reshape ( int ( [ &
1 , &
2 , &
4 , &
3 , &
5 , &
6 , &
8 , &
7 &
] , pInt ) , [ FE_maxNnodesAtIP ( me ) , FE_Nips ( me ) ] )
me = me + 1_pInt
FE_nodesAtIP ( 1 : FE_maxNnodesAtIP ( me ) , 1 : FE_Nips ( me ) , me ) = & ! element 21 (3D 20node 27ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
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-11-16 04:15:20 +05:30
] , pInt ) , [ FE_maxNnodesAtIP ( me ) , FE_Nips ( me ) ] )
2012-04-11 22:58:08 +05:30
2018-01-10 21:43:25 +05:30
2012-06-15 21:40:21 +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.
2012-11-16 04:15:20 +05:30
me = 0_pInt
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_ipNeighbor ( 1 : FE_NipNeighbors ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 6 (2D 3node 1ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2012-11-16 04:15:20 +05:30
- 2 , - 3 , - 1 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_NipNeighbors ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2018-01-10 21:43:25 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_ipNeighbor ( 1 : FE_NipNeighbors ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 125 (2D 6node 3ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2012-11-16 04:15:20 +05:30
2 , - 3 , 3 , - 1 , &
- 2 , 1 , 3 , - 1 , &
2 , - 3 , - 2 , 1 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_NipNeighbors ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2018-01-10 21:43:25 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_ipNeighbor ( 1 : FE_NipNeighbors ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 11 (2D 4node 4ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2 , - 4 , 3 , - 1 , &
- 2 , 1 , 4 , - 1 , &
4 , - 4 , - 3 , 1 , &
- 2 , 3 , - 3 , 2 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_NipNeighbors ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_ipNeighbor ( 1 : FE_NipNeighbors ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 27 (2D 8node 9ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
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 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_NipNeighbors ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_ipNeighbor ( 1 : FE_NipNeighbors ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 134 (3D 4node 1ip)
2012-11-16 04:15:20 +05:30
reshape ( int ( [ &
- 1 , - 2 , - 3 , - 4 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_NipNeighbors ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_ipNeighbor ( 1 : FE_NipNeighbors ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 127 (3D 10node 4ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2 , - 4 , 3 , - 2 , 4 , - 1 , &
2012-11-16 04:15:20 +05:30
- 2 , 1 , 3 , - 2 , 4 , - 1 , &
2 , - 4 , - 3 , 1 , 4 , - 1 , &
2 , - 4 , 3 , - 2 , - 3 , 1 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_NipNeighbors ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_ipNeighbor ( 1 : FE_NipNeighbors ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 136 (3D 6node 6ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
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 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_NipNeighbors ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_ipNeighbor ( 1 : FE_NipNeighbors ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 117 (3D 8node 1ip)
2012-11-16 04:15:20 +05:30
reshape ( int ( [ &
- 3 , - 5 , - 4 , - 2 , - 6 , - 1 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_NipNeighbors ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_ipNeighbor ( 1 : FE_NipNeighbors ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 7 (3D 8node 8ip)
2012-11-16 04:15:20 +05:30
reshape ( int ( [ &
2 , - 5 , 3 , - 2 , 5 , - 1 , &
- 3 , 1 , 4 , - 2 , 6 , - 1 , &
4 , - 5 , - 4 , 1 , 7 , - 1 , &
- 3 , 3 , - 4 , 2 , 8 , - 1 , &
6 , - 5 , 7 , - 2 , - 6 , 1 , &
- 3 , 5 , 8 , - 2 , - 6 , 2 , &
8 , - 5 , - 4 , 5 , - 6 , 3 , &
- 3 , 7 , - 4 , 6 , - 6 , 4 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_NipNeighbors ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_ipNeighbor ( 1 : FE_NipNeighbors ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 21 (3D 20node 27ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2 , - 5 , 4 , - 2 , 10 , - 1 , &
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 , &
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 , &
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 &
2013-04-15 13:43:20 +05:30
] , pInt ) , [ FE_NipNeighbors ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-06-15 21:40:21 +05:30
2013-04-15 13:43:20 +05:30
! *** FE_cell ***
2012-11-16 04:15:20 +05:30
me = 0_pInt
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cell ( 1 : FE_NcellnodesPerCell ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 6 (2D 3node 1ip)
reshape ( int ( [ &
1 , 2 , 3 &
] , pInt ) , [ FE_NcellnodesPerCell ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cell ( 1 : FE_NcellnodesPerCell ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 125 (2D 6node 3ip)
2012-11-16 04:15:20 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
1 , 4 , 7 , 6 , &
2 , 5 , 7 , 4 , &
3 , 6 , 7 , 5 &
] , pInt ) , [ FE_NcellnodesPerCell ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cell ( 1 : FE_NcellnodesPerCell ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 11 (2D 4node 4ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
1 , 5 , 9 , 8 , &
5 , 2 , 6 , 9 , &
8 , 9 , 7 , 4 , &
9 , 6 , 3 , 7 &
] , pInt ) , [ FE_NcellnodesPerCell ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cell ( 1 : FE_NcellnodesPerCell ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 27 (2D 8node 9ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
1 , 5 , 13 , 12 , &
5 , 6 , 14 , 13 , &
6 , 2 , 7 , 14 , &
12 , 13 , 16 , 11 , &
13 , 14 , 15 , 16 , &
14 , 7 , 8 , 15 , &
11 , 16 , 10 , 4 , &
16 , 15 , 9 , 10 , &
15 , 8 , 3 , 9 &
] , pInt ) , [ FE_NcellnodesPerCell ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cell ( 1 : FE_NcellnodesPerCell ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 134 (3D 4node 1ip)
reshape ( int ( [ &
1 , 2 , 3 , 4 &
] , pInt ) , [ FE_NcellnodesPerCell ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cell ( 1 : FE_NcellnodesPerCell ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 127 (3D 10node 4ip)
2012-11-16 04:15:20 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
1 , 5 , 11 , 7 , 8 , 12 , 15 , 14 , &
5 , 2 , 6 , 11 , 12 , 9 , 13 , 15 , &
7 , 11 , 6 , 3 , 14 , 15 , 13 , 10 , &
8 , 12 , 15 , 4 , 4 , 9 , 13 , 10 &
] , pInt ) , [ FE_NcellnodesPerCell ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cell ( 1 : FE_NcellnodesPerCell ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 136 (3D 6node 6ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
1 , 7 , 16 , 9 , 10 , 17 , 21 , 19 , &
7 , 2 , 8 , 16 , 17 , 11 , 18 , 21 , &
9 , 16 , 8 , 3 , 19 , 21 , 18 , 12 , &
10 , 17 , 21 , 19 , 4 , 13 , 20 , 15 , &
17 , 11 , 18 , 21 , 13 , 5 , 14 , 20 , &
19 , 21 , 18 , 12 , 15 , 20 , 14 , 6 &
] , pInt ) , [ FE_NcellnodesPerCell ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cell ( 1 : FE_NcellnodesPerCell ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 117 (3D 8node 1ip)
reshape ( int ( [ &
1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 &
] , pInt ) , [ FE_NcellnodesPerCell ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cell ( 1 : FE_NcellnodesPerCell ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 7 (3D 8node 8ip)
2012-11-16 04:15:20 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
1 , 9 , 21 , 12 , 13 , 22 , 27 , 25 , &
9 , 2 , 10 , 21 , 22 , 14 , 23 , 27 , &
12 , 21 , 11 , 4 , 25 , 27 , 24 , 16 , &
21 , 10 , 3 , 11 , 27 , 23 , 15 , 24 , &
13 , 22 , 27 , 25 , 5 , 17 , 26 , 20 , &
22 , 14 , 23 , 27 , 17 , 6 , 18 , 26 , &
25 , 27 , 24 , 16 , 20 , 26 , 19 , 8 , &
27 , 23 , 15 , 24 , 26 , 18 , 7 , 19 &
] , pInt ) , [ FE_NcellnodesPerCell ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cell ( 1 : FE_NcellnodesPerCell ( FE_celltype ( me ) ) , 1 : FE_Nips ( me ) , me ) = & ! element 21 (3D 20node 27ip)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
1 , 9 , 33 , 16 , 17 , 37 , 57 , 44 , &
9 , 10 , 34 , 33 , 37 , 38 , 58 , 57 , &
10 , 2 , 11 , 34 , 38 , 18 , 39 , 58 , &
16 , 33 , 36 , 15 , 44 , 57 , 60 , 43 , &
33 , 34 , 35 , 36 , 57 , 58 , 59 , 60 , &
34 , 11 , 12 , 35 , 58 , 39 , 40 , 59 , &
15 , 36 , 14 , 4 , 43 , 60 , 42 , 20 , &
36 , 35 , 13 , 14 , 60 , 59 , 41 , 42 , &
35 , 12 , 3 , 13 , 59 , 40 , 19 , 41 , &
17 , 37 , 57 , 44 , 21 , 45 , 61 , 52 , &
37 , 38 , 58 , 57 , 45 , 46 , 62 , 61 , &
38 , 18 , 39 , 58 , 46 , 22 , 47 , 62 , &
44 , 57 , 60 , 43 , 52 , 61 , 64 , 51 , &
57 , 58 , 59 , 60 , 61 , 62 , 63 , 64 , &
58 , 39 , 40 , 59 , 62 , 47 , 48 , 63 , &
43 , 60 , 42 , 20 , 51 , 64 , 50 , 24 , &
60 , 59 , 41 , 42 , 64 , 63 , 49 , 50 , &
59 , 40 , 19 , 41 , 63 , 48 , 23 , 49 , &
21 , 45 , 61 , 52 , 5 , 25 , 53 , 32 , &
45 , 46 , 62 , 61 , 25 , 26 , 54 , 53 , &
46 , 22 , 47 , 62 , 26 , 6 , 27 , 54 , &
52 , 61 , 64 , 51 , 32 , 53 , 56 , 31 , &
61 , 62 , 63 , 64 , 53 , 54 , 55 , 56 , &
62 , 47 , 48 , 63 , 54 , 27 , 28 , 55 , &
51 , 64 , 50 , 24 , 31 , 56 , 30 , 8 , &
64 , 63 , 49 , 50 , 56 , 55 , 29 , 30 , &
63 , 48 , 23 , 49 , 55 , 28 , 7 , 29 &
] , pInt ) , [ FE_NcellnodesPerCell ( FE_celltype ( me ) ) , FE_Nips ( me ) ] )
! *** FE_cellnodeParentnodeWeights ***
! center of gravity of the weighted nodes gives the position of the cell node.
! fill with 0.
! example: face-centered cell node with face nodes 1,2,5,6 to be used in,
! e.g., an 8 node element, would be encoded:
! 1, 1, 0, 0, 1, 1, 0, 0
2012-11-16 04:15:20 +05:30
me = 0_pInt
2013-04-15 13:43:20 +05:30
me = me + 1_pInt
2018-01-10 21:43:25 +05:30
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 6 (2D 3node 1ip)
2013-04-15 13:43:20 +05:30
reshape ( real ( [ &
2018-01-10 21:43:25 +05:30
1 , 0 , 0 , &
0 , 1 , 0 , &
2013-04-15 13:43:20 +05:30
0 , 0 , 1 &
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
me = me + 1_pInt
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 125 (2D 6node 3ip)
reshape ( real ( [ &
2018-01-10 21:43:25 +05:30
1 , 0 , 0 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , 0 , 0 , &
2013-04-15 13:43:20 +05:30
0 , 0 , 1 , 0 , 0 , 0 , &
0 , 0 , 0 , 1 , 0 , 0 , &
0 , 0 , 0 , 0 , 1 , 0 , &
0 , 0 , 0 , 0 , 0 , 1 , &
1 , 1 , 1 , 2 , 2 , 2 &
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
2018-01-10 21:43:25 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 11 (2D 4node 4ip)
reshape ( real ( [ &
2018-01-10 21:43:25 +05:30
1 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , &
2013-04-15 13:43:20 +05:30
0 , 0 , 1 , 0 , &
2018-01-10 21:43:25 +05:30
0 , 0 , 0 , 1 , &
2013-04-15 13:43:20 +05:30
1 , 1 , 0 , 0 , &
0 , 1 , 1 , 0 , &
0 , 0 , 1 , 1 , &
1 , 0 , 0 , 1 , &
1 , 1 , 1 , 1 &
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 27 (2D 8node 9ip)
reshape ( real ( [ &
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , &
0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , &
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , &
1 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , 2 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , 0 , 2 , 0 , 0 , &
0 , 0 , 1 , 0 , 0 , 2 , 0 , 0 , &
0 , 0 , 1 , 0 , 0 , 0 , 2 , 0 , &
0 , 0 , 0 , 1 , 0 , 0 , 2 , 0 , &
0 , 0 , 0 , 1 , 0 , 0 , 0 , 2 , &
1 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , &
4 , 1 , 1 , 1 , 8 , 2 , 2 , 8 , &
1 , 4 , 1 , 1 , 8 , 8 , 2 , 2 , &
1 , 1 , 4 , 1 , 2 , 8 , 8 , 2 , &
1 , 1 , 1 , 4 , 2 , 2 , 8 , 8 &
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 54 (2D 8node 4ip)
reshape ( real ( [ &
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , &
0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , &
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , &
0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , &
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , &
1 , 1 , 1 , 1 , 2 , 2 , 2 , 2 &
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2018-01-10 21:43:25 +05:30
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 134 (3D 4node 1ip)
2013-04-15 13:43:20 +05:30
reshape ( real ( [ &
2018-01-10 21:43:25 +05:30
1 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , &
2013-04-15 13:43:20 +05:30
0 , 0 , 1 , 0 , &
2018-01-10 21:43:25 +05:30
0 , 0 , 0 , 1 &
2013-04-15 13:43:20 +05:30
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2018-01-10 21:43:25 +05:30
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 157 (3D 5node 4ip)
2013-04-15 13:43:20 +05:30
reshape ( real ( [ &
1 , 0 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , 0 , &
0 , 0 , 1 , 0 , 0 , &
0 , 0 , 0 , 1 , 0 , &
1 , 1 , 0 , 0 , 0 , &
0 , 1 , 1 , 0 , 0 , &
1 , 0 , 1 , 0 , 0 , &
1 , 0 , 0 , 1 , 0 , &
0 , 1 , 0 , 1 , 0 , &
0 , 0 , 1 , 1 , 0 , &
1 , 1 , 1 , 0 , 0 , &
1 , 1 , 0 , 1 , 0 , &
0 , 1 , 1 , 1 , 0 , &
1 , 0 , 1 , 1 , 0 , &
0 , 0 , 0 , 0 , 1 &
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 127 (3D 10node 4ip)
reshape ( real ( [ &
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , &
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , &
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , &
1 , 1 , 1 , 0 , 2 , 2 , 2 , 0 , 0 , 0 , &
1 , 1 , 0 , 1 , 2 , 0 , 0 , 2 , 2 , 0 , &
0 , 1 , 1 , 1 , 0 , 2 , 0 , 0 , 2 , 2 , &
1 , 0 , 1 , 1 , 0 , 0 , 2 , 2 , 0 , 2 , &
3 , 3 , 3 , 3 , 4 , 4 , 4 , 4 , 4 , 4 &
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
me = me + 1_pInt
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 136 (3D 6node 6ip)
reshape ( real ( [ &
1 , 0 , 0 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , 0 , 0 , &
0 , 0 , 1 , 0 , 0 , 0 , &
0 , 0 , 0 , 1 , 0 , 0 , &
0 , 0 , 0 , 0 , 1 , 0 , &
0 , 0 , 0 , 0 , 0 , 1 , &
1 , 1 , 0 , 0 , 0 , 0 , &
0 , 1 , 1 , 0 , 0 , 0 , &
1 , 0 , 1 , 0 , 0 , 0 , &
1 , 0 , 0 , 1 , 0 , 0 , &
0 , 1 , 0 , 0 , 1 , 0 , &
0 , 0 , 1 , 0 , 0 , 1 , &
0 , 0 , 0 , 1 , 1 , 0 , &
0 , 0 , 0 , 0 , 1 , 1 , &
0 , 0 , 0 , 1 , 0 , 1 , &
1 , 1 , 1 , 0 , 0 , 0 , &
1 , 1 , 0 , 1 , 1 , 0 , &
0 , 1 , 1 , 0 , 1 , 1 , &
1 , 0 , 1 , 1 , 0 , 1 , &
0 , 0 , 0 , 1 , 1 , 1 , &
1 , 1 , 1 , 1 , 1 , 1 &
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
me = me + 1_pInt
2018-01-10 21:43:25 +05:30
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 117 (3D 8node 1ip)
2013-04-15 13:43:20 +05:30
reshape ( real ( [ &
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , &
0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , &
0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , &
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , &
0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , &
0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , &
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 &
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
me = me + 1_pInt
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 7 (3D 8node 8ip)
reshape ( real ( [ &
2018-01-10 21:43:25 +05:30
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , & !
2013-04-15 13:43:20 +05:30
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , & ! 5
2018-01-10 21:43:25 +05:30
0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , & !
1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , & !
2013-04-15 13:43:20 +05:30
0 , 1 , 1 , 0 , 0 , 0 , 0 , 0 , & ! 10
2018-01-10 21:43:25 +05:30
0 , 0 , 1 , 1 , 0 , 0 , 0 , 0 , & !
1 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , & !
1 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , & !
0 , 1 , 0 , 0 , 0 , 1 , 0 , 0 , & !
2013-04-15 13:43:20 +05:30
0 , 0 , 1 , 0 , 0 , 0 , 1 , 0 , & ! 15
2018-01-10 21:43:25 +05:30
0 , 0 , 0 , 1 , 0 , 0 , 0 , 1 , & !
0 , 0 , 0 , 0 , 1 , 1 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 1 , 1 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 1 , 1 , & !
2013-04-15 13:43:20 +05:30
0 , 0 , 0 , 0 , 1 , 0 , 0 , 1 , & ! 20
2018-01-10 21:43:25 +05:30
1 , 1 , 1 , 1 , 0 , 0 , 0 , 0 , & !
1 , 1 , 0 , 0 , 1 , 1 , 0 , 0 , & !
0 , 1 , 1 , 0 , 0 , 1 , 1 , 0 , & !
0 , 0 , 1 , 1 , 0 , 0 , 1 , 1 , & !
2013-04-15 13:43:20 +05:30
1 , 0 , 0 , 1 , 1 , 0 , 0 , 1 , & ! 25
2018-01-10 21:43:25 +05:30
0 , 0 , 0 , 0 , 1 , 1 , 1 , 1 , & !
1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 & !
2013-04-15 13:43:20 +05:30
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
me = me + 1_pInt
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 57 (3D 20node 8ip)
reshape ( real ( [ &
2018-01-10 21:43:25 +05:30
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & ! 5
0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & ! 10
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 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 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , & ! 15
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , & ! 20
1 , 1 , 1 , 1 , 0 , 0 , 0 , 0 , 2 , 2 , 2 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
1 , 1 , 0 , 0 , 1 , 1 , 0 , 0 , 2 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 2 , 2 , 0 , 0 , & !
0 , 1 , 1 , 0 , 0 , 1 , 1 , 0 , 0 , 2 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 2 , 2 , 0 , & !
0 , 0 , 1 , 1 , 0 , 0 , 1 , 1 , 0 , 0 , 2 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 2 , 2 , & !
1 , 0 , 0 , 1 , 1 , 0 , 0 , 1 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 2 , 2 , 0 , 0 , 2 , & ! 25
0 , 0 , 0 , 0 , 1 , 1 , 1 , 1 , 0 , 0 , 0 , 0 , 2 , 2 , 2 , 2 , 0 , 0 , 0 , 0 , & !
3 , 3 , 3 , 3 , 3 , 3 , 3 , 3 , 4 , 4 , 4 , 4 , 4 , 4 , 4 , 4 , 4 , 4 , 4 , 4 & !
2013-04-15 13:43:20 +05:30
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-15 13:43:20 +05:30
FE_cellnodeParentnodeWeights ( 1 : FE_Nnodes ( me ) , 1 : FE_Ncellnodes ( FE_geomtype ( me ) ) , me ) = & ! element 21 (3D 20node 27ip)
reshape ( real ( [ &
2018-01-10 21:43:25 +05:30
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & ! 5
0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & ! 10
0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & ! 15
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , & !
0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , & !
0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , & !
0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , & ! 20
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , & !
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & ! 25
0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , & ! 30
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , & !
4 , 1 , 1 , 1 , 0 , 0 , 0 , 0 , 8 , 2 , 2 , 8 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
1 , 4 , 1 , 1 , 0 , 0 , 0 , 0 , 8 , 8 , 2 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
1 , 1 , 4 , 1 , 0 , 0 , 0 , 0 , 2 , 8 , 8 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & ! 35
1 , 1 , 1 , 4 , 0 , 0 , 0 , 0 , 2 , 2 , 8 , 8 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , & !
4 , 1 , 0 , 0 , 1 , 1 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 2 , 0 , 0 , & !
1 , 4 , 0 , 0 , 1 , 1 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 2 , 8 , 0 , 0 , & !
0 , 4 , 1 , 0 , 0 , 1 , 1 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 2 , 0 , & !
0 , 1 , 4 , 0 , 0 , 1 , 1 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 2 , 8 , 0 , & ! 40
0 , 0 , 4 , 1 , 0 , 0 , 1 , 1 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 2 , & !
0 , 0 , 1 , 4 , 0 , 0 , 1 , 1 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 2 , 8 , & !
1 , 0 , 0 , 4 , 1 , 0 , 0 , 1 , 0 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 2 , 0 , 0 , 8 , & !
4 , 0 , 0 , 1 , 1 , 0 , 0 , 1 , 0 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 8 , 0 , 0 , 2 , & !
1 , 1 , 0 , 0 , 4 , 1 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 0 , 0 , 0 , 8 , 2 , 0 , 0 , & ! 45
1 , 1 , 0 , 0 , 1 , 4 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 8 , 0 , 0 , & !
0 , 1 , 1 , 0 , 0 , 4 , 1 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 0 , 0 , 0 , 8 , 2 , 0 , & !
0 , 1 , 1 , 0 , 0 , 1 , 4 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 8 , 0 , & !
0 , 0 , 1 , 1 , 0 , 0 , 4 , 1 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 0 , 0 , 0 , 8 , 2 , & !
0 , 0 , 1 , 1 , 0 , 0 , 1 , 4 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 0 , 0 , 0 , 2 , 8 , & ! 50
1 , 0 , 0 , 1 , 1 , 0 , 0 , 4 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 2 , 0 , 0 , 8 , & !
1 , 0 , 0 , 1 , 4 , 0 , 0 , 1 , 0 , 0 , 0 , 2 , 0 , 0 , 0 , 8 , 8 , 0 , 0 , 2 , & !
0 , 0 , 0 , 0 , 4 , 1 , 1 , 1 , 0 , 0 , 0 , 0 , 8 , 2 , 2 , 8 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 1 , 4 , 1 , 1 , 0 , 0 , 0 , 0 , 8 , 8 , 2 , 2 , 0 , 0 , 0 , 0 , & !
0 , 0 , 0 , 0 , 1 , 1 , 4 , 1 , 0 , 0 , 0 , 0 , 2 , 8 , 8 , 2 , 0 , 0 , 0 , 0 , & ! 55
0 , 0 , 0 , 0 , 1 , 1 , 1 , 4 , 0 , 0 , 0 , 0 , 2 , 2 , 8 , 8 , 0 , 0 , 0 , 0 , & !
24 , 8 , 4 , 8 , 8 , 4 , 3 , 4 , 32 , 12 , 12 , 32 , 12 , 4 , 4 , 12 , 32 , 12 , 4 , 12 , & !
8 , 24 , 8 , 4 , 4 , 8 , 4 , 3 , 32 , 32 , 12 , 12 , 12 , 12 , 4 , 4 , 12 , 32 , 12 , 4 , & !
4 , 8 , 24 , 8 , 3 , 4 , 8 , 4 , 12 , 32 , 32 , 12 , 4 , 12 , 12 , 4 , 4 , 12 , 32 , 12 , & !
8 , 4 , 8 , 24 , 4 , 3 , 4 , 8 , 12 , 12 , 32 , 32 , 4 , 4 , 12 , 12 , 12 , 4 , 12 , 32 , & ! 60
8 , 4 , 3 , 4 , 24 , 8 , 4 , 8 , 12 , 4 , 4 , 12 , 32 , 12 , 12 , 32 , 32 , 12 , 4 , 12 , & !
4 , 8 , 4 , 3 , 8 , 24 , 8 , 4 , 12 , 12 , 4 , 4 , 32 , 32 , 12 , 12 , 12 , 32 , 12 , 4 , & !
3 , 4 , 8 , 4 , 4 , 8 , 24 , 8 , 4 , 12 , 12 , 4 , 12 , 32 , 32 , 12 , 4 , 12 , 32 , 12 , & !
4 , 3 , 4 , 8 , 8 , 4 , 8 , 24 , 4 , 4 , 12 , 12 , 12 , 12 , 32 , 32 , 12 , 4 , 12 , 32 & !
2013-04-15 13:43:20 +05:30
] , pReal ) , [ FE_Nnodes ( me ) , FE_Ncellnodes ( FE_geomtype ( me ) ) ] )
! *** FE_cellface ***
me = 0_pInt
me = me + 1_pInt
2013-04-18 22:10:49 +05:30
FE_cellface ( 1 : FE_NcellnodesPerCellface ( me ) , 1 : FE_NipNeighbors ( me ) , me ) = & ! 2D 3node, VTK_TRIANGLE (5)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
2 , 3 , &
3 , 1 , &
1 , 2 &
] , pInt ) , [ FE_NcellnodesPerCellface ( me ) , FE_NipNeighbors ( me ) ] )
2012-06-15 21:40:21 +05:30
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-18 22:10:49 +05:30
FE_cellface ( 1 : FE_NcellnodesPerCellface ( me ) , 1 : FE_NipNeighbors ( me ) , me ) = & ! 2D 4node, VTK_QUAD (9)
2012-11-16 04:15:20 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
2 , 3 , &
4 , 1 , &
3 , 4 , &
1 , 2 &
] , pInt ) , [ FE_NcellnodesPerCellface ( me ) , FE_NipNeighbors ( me ) ] )
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-18 22:10:49 +05:30
FE_cellface ( 1 : FE_NcellnodesPerCellface ( me ) , 1 : FE_NipNeighbors ( me ) , me ) = & ! 3D 4node, VTK_TETRA (10)
2012-11-16 04:15:20 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
1 , 3 , 2 , &
1 , 2 , 4 , &
2 , 3 , 4 , &
1 , 4 , 3 &
] , pInt ) , [ FE_NcellnodesPerCellface ( me ) , FE_NipNeighbors ( me ) ] )
2012-11-16 04:15:20 +05:30
me = me + 1_pInt
2013-04-18 22:10:49 +05:30
FE_cellface ( 1 : FE_NcellnodesPerCellface ( me ) , 1 : FE_NipNeighbors ( me ) , me ) = & ! 3D 8node, VTK_HEXAHEDRON (12)
2012-06-15 21:40:21 +05:30
reshape ( int ( [ &
2013-04-15 13:43:20 +05:30
2 , 3 , 7 , 6 , &
4 , 1 , 5 , 8 , &
3 , 4 , 8 , 7 , &
1 , 2 , 6 , 5 , &
5 , 6 , 7 , 8 , &
1 , 4 , 3 , 2 &
] , pInt ) , [ FE_NcellnodesPerCellface ( me ) , FE_NipNeighbors ( me ) ] )
2012-06-15 21:40:21 +05:30
end subroutine mesh_build_FEdata
2012-08-16 20:25:23 +05:30
end module mesh