diff --git a/docs/examples/Mesh/MeshHex8-FineLayer-nodes_paraview.py b/docs/examples/Mesh/MeshHex8-FineLayer-nodes_paraview.py index 47893db..9bce4d7 100644 --- a/docs/examples/Mesh/MeshHex8-FineLayer-nodes_paraview.py +++ b/docs/examples/Mesh/MeshHex8-FineLayer-nodes_paraview.py @@ -1,362 +1,364 @@ import numpy as np import h5py import GooseFEM as gf # ====================== create fictitious configuration + store to HDF5-file ====================== # create mesh object mesh = gf.Mesh.Hex8.FineLayer(9,17,27) # set node identifier for all boundary sets Z = lambda mesh: np.zeros((mesh.nnode()),dtype='int') nodesFront = Z(mesh); nodesFront [mesh.nodesFront() ] = 1 nodesBack = Z(mesh); nodesBack [mesh.nodesBack() ] = 1 nodesLeft = Z(mesh); nodesLeft [mesh.nodesLeft() ] = 1 nodesRight = Z(mesh); nodesRight [mesh.nodesRight() ] = 1 nodesBottom = Z(mesh); nodesBottom [mesh.nodesBottom() ] = 1 nodesTop = Z(mesh); nodesTop [mesh.nodesTop() ] = 1 nodesFrontFace = Z(mesh); nodesFrontFace [mesh.nodesFrontFace() ] = 1 nodesBackFace = Z(mesh); nodesBackFace [mesh.nodesBackFace() ] = 1 nodesLeftFace = Z(mesh); nodesLeftFace [mesh.nodesLeftFace() ] = 1 nodesRightFace = Z(mesh); nodesRightFace [mesh.nodesRightFace() ] = 1 nodesBottomFace = Z(mesh); nodesBottomFace [mesh.nodesBottomFace() ] = 1 nodesTopFace = Z(mesh); nodesTopFace [mesh.nodesTopFace() ] = 1 nodesFrontBottomEdge = Z(mesh); nodesFrontBottomEdge [mesh.nodesFrontBottomEdge() ] = 1 nodesFrontTopEdge = Z(mesh); nodesFrontTopEdge [mesh.nodesFrontTopEdge() ] = 1 nodesFrontLeftEdge = Z(mesh); nodesFrontLeftEdge [mesh.nodesFrontLeftEdge() ] = 1 nodesFrontRightEdge = Z(mesh); nodesFrontRightEdge [mesh.nodesFrontRightEdge() ] = 1 nodesBackBottomEdge = Z(mesh); nodesBackBottomEdge [mesh.nodesBackBottomEdge() ] = 1 nodesBackTopEdge = Z(mesh); nodesBackTopEdge [mesh.nodesBackTopEdge() ] = 1 nodesBackLeftEdge = Z(mesh); nodesBackLeftEdge [mesh.nodesBackLeftEdge() ] = 1 nodesBackRightEdge = Z(mesh); nodesBackRightEdge [mesh.nodesBackRightEdge() ] = 1 nodesBottomLeftEdge = Z(mesh); nodesBottomLeftEdge [mesh.nodesBottomLeftEdge() ] = 1 nodesBottomRightEdge = Z(mesh); nodesBottomRightEdge [mesh.nodesBottomRightEdge() ] = 1 nodesTopLeftEdge = Z(mesh); nodesTopLeftEdge [mesh.nodesTopLeftEdge() ] = 1 nodesTopRightEdge = Z(mesh); nodesTopRightEdge [mesh.nodesTopRightEdge() ] = 1 nodesFrontBottomOpenEdge = Z(mesh); nodesFrontBottomOpenEdge [mesh.nodesFrontBottomOpenEdge() ] = 1 nodesFrontTopOpenEdge = Z(mesh); nodesFrontTopOpenEdge [mesh.nodesFrontTopOpenEdge() ] = 1 nodesFrontLeftOpenEdge = Z(mesh); nodesFrontLeftOpenEdge [mesh.nodesFrontLeftOpenEdge() ] = 1 nodesFrontRightOpenEdge = Z(mesh); nodesFrontRightOpenEdge [mesh.nodesFrontRightOpenEdge() ] = 1 nodesBackBottomOpenEdge = Z(mesh); nodesBackBottomOpenEdge [mesh.nodesBackBottomOpenEdge() ] = 1 nodesBackTopOpenEdge = Z(mesh); nodesBackTopOpenEdge [mesh.nodesBackTopOpenEdge() ] = 1 nodesBackLeftOpenEdge = Z(mesh); nodesBackLeftOpenEdge [mesh.nodesBackLeftOpenEdge() ] = 1 nodesBackRightOpenEdge = Z(mesh); nodesBackRightOpenEdge [mesh.nodesBackRightOpenEdge() ] = 1 nodesBottomLeftOpenEdge = Z(mesh); nodesBottomLeftOpenEdge [mesh.nodesBottomLeftOpenEdge() ] = 1 nodesBottomRightOpenEdge = Z(mesh); nodesBottomRightOpenEdge [mesh.nodesBottomRightOpenEdge() ] = 1 nodesTopLeftOpenEdge = Z(mesh); nodesTopLeftOpenEdge [mesh.nodesTopLeftOpenEdge() ] = 1 nodesTopRightOpenEdge = Z(mesh); nodesTopRightOpenEdge [mesh.nodesTopRightOpenEdge() ] = 1 nodesFrontBottomLeftCorner = Z(mesh); nodesFrontBottomLeftCorner [mesh.nodesFrontBottomLeftCorner() ] = 1 nodesFrontBottomRightCorner = Z(mesh); nodesFrontBottomRightCorner[mesh.nodesFrontBottomRightCorner()] = 1 nodesFrontTopLeftCorner = Z(mesh); nodesFrontTopLeftCorner [mesh.nodesFrontTopLeftCorner() ] = 1 nodesFrontTopRightCorner = Z(mesh); nodesFrontTopRightCorner [mesh.nodesFrontTopRightCorner() ] = 1 nodesBackBottomLeftCorner = Z(mesh); nodesBackBottomLeftCorner [mesh.nodesBackBottomLeftCorner() ] = 1 nodesBackBottomRightCorner = Z(mesh); nodesBackBottomRightCorner [mesh.nodesBackBottomRightCorner() ] = 1 nodesBackTopLeftCorner = Z(mesh); nodesBackTopLeftCorner [mesh.nodesBackTopLeftCorner() ] = 1 nodesBackTopRightCorner = Z(mesh); nodesBackTopRightCorner [mesh.nodesBackTopRightCorner() ] = 1 # DOFs after eliminating periodicity dofsPeriodic = mesh.dofsPeriodic()[:,0] # open data file f = h5py.File('MeshHex8-FineLayer-nodes_paraview.hdf5','w') -# write particle positions, and a dummy connectivity +# write nodal positions, and a dummy connectivity f.create_dataset('/coor' ,data=mesh.coor() ) f.create_dataset('/conn' ,data=np.arange(mesh.nnode()) ) f.create_dataset('/nodesFront' ,data=nodesFront ) f.create_dataset('/nodesBack' ,data=nodesBack ) f.create_dataset('/nodesLeft' ,data=nodesLeft ) f.create_dataset('/nodesRight' ,data=nodesRight ) f.create_dataset('/nodesBottom' ,data=nodesBottom ) f.create_dataset('/nodesTop' ,data=nodesTop ) f.create_dataset('/nodesFrontFace' ,data=nodesFrontFace ) f.create_dataset('/nodesBackFace' ,data=nodesBackFace ) f.create_dataset('/nodesLeftFace' ,data=nodesLeftFace ) f.create_dataset('/nodesRightFace' ,data=nodesRightFace ) f.create_dataset('/nodesBottomFace' ,data=nodesBottomFace ) f.create_dataset('/nodesTopFace' ,data=nodesTopFace ) f.create_dataset('/nodesFrontBottomEdge' ,data=nodesFrontBottomEdge ) f.create_dataset('/nodesFrontTopEdge' ,data=nodesFrontTopEdge ) f.create_dataset('/nodesFrontLeftEdge' ,data=nodesFrontLeftEdge ) f.create_dataset('/nodesFrontRightEdge' ,data=nodesFrontRightEdge ) f.create_dataset('/nodesBackBottomEdge' ,data=nodesBackBottomEdge ) f.create_dataset('/nodesBackTopEdge' ,data=nodesBackTopEdge ) f.create_dataset('/nodesBackLeftEdge' ,data=nodesBackLeftEdge ) f.create_dataset('/nodesBackRightEdge' ,data=nodesBackRightEdge ) f.create_dataset('/nodesBottomLeftEdge' ,data=nodesBottomLeftEdge ) f.create_dataset('/nodesBottomRightEdge' ,data=nodesBottomRightEdge ) f.create_dataset('/nodesTopLeftEdge' ,data=nodesTopLeftEdge ) f.create_dataset('/nodesTopRightEdge' ,data=nodesTopRightEdge ) f.create_dataset('/nodesFrontBottomOpenEdge' ,data=nodesFrontBottomOpenEdge ) f.create_dataset('/nodesFrontTopOpenEdge' ,data=nodesFrontTopOpenEdge ) f.create_dataset('/nodesFrontLeftOpenEdge' ,data=nodesFrontLeftOpenEdge ) f.create_dataset('/nodesFrontRightOpenEdge' ,data=nodesFrontRightOpenEdge ) f.create_dataset('/nodesBackBottomOpenEdge' ,data=nodesBackBottomOpenEdge ) f.create_dataset('/nodesBackTopOpenEdge' ,data=nodesBackTopOpenEdge ) f.create_dataset('/nodesBackLeftOpenEdge' ,data=nodesBackLeftOpenEdge ) f.create_dataset('/nodesBackRightOpenEdge' ,data=nodesBackRightOpenEdge ) f.create_dataset('/nodesBottomLeftOpenEdge' ,data=nodesBottomLeftOpenEdge ) f.create_dataset('/nodesBottomRightOpenEdge' ,data=nodesBottomRightOpenEdge ) f.create_dataset('/nodesTopLeftOpenEdge' ,data=nodesTopLeftOpenEdge ) f.create_dataset('/nodesTopRightOpenEdge' ,data=nodesTopRightOpenEdge ) f.create_dataset('/nodesFrontBottomLeftCorner' ,data=nodesFrontBottomLeftCorner ) f.create_dataset('/nodesFrontBottomRightCorner',data=nodesFrontBottomRightCorner) f.create_dataset('/nodesFrontTopLeftCorner' ,data=nodesFrontTopLeftCorner ) f.create_dataset('/nodesFrontTopRightCorner' ,data=nodesFrontTopRightCorner ) f.create_dataset('/nodesBackBottomLeftCorner' ,data=nodesBackBottomLeftCorner ) f.create_dataset('/nodesBackBottomRightCorner' ,data=nodesBackBottomRightCorner ) f.create_dataset('/nodesBackTopLeftCorner' ,data=nodesBackTopLeftCorner ) f.create_dataset('/nodesBackTopRightCorner' ,data=nodesBackTopRightCorner ) f.create_dataset('/dofsPeriodic' ,data=dofsPeriodic ) # ======================================== write XDMF-file ========================================= +# basic file-format xmf = ''' MeshHex8-FineLayer-nodes_paraview.hdf5:/conn MeshHex8-FineLayer-nodes_paraview.hdf5:/coor MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFront MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBack MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesLeft MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesRight MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottom MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTop MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontFace MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackFace MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesLeftFace MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesRightFace MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottomFace MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTopFace MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontBottomEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontTopEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontLeftEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontRightEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackBottomEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackTopEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackLeftEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackRightEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottomLeftEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottomRightEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTopLeftEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTopRightEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontBottomOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontTopOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontLeftOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontRightOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackBottomOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackTopOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackLeftOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackRightOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottomLeftOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBottomRightOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTopLeftOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesTopRightOpenEdge MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontBottomLeftCorner MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontBottomRightCorner MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontTopLeftCorner MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesFrontTopRightCorner MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackBottomLeftCorner MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackBottomRightCorner MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackTopLeftCorner MeshHex8-FineLayer-nodes_paraview.hdf5:/nodesBackTopRightCorner MeshHex8-FineLayer-nodes_paraview.hdf5:/dofsPeriodic ''' +# write to file, fill mesh dimensions open('MeshHex8-FineLayer-nodes_paraview.xmf','w').write(xmf.format( nnode = mesh.nnode(), )) diff --git a/docs/examples/Mesh/MeshHex8-FineLayer_paraview.py b/docs/examples/Mesh/MeshHex8-FineLayer_paraview.py index c04f6a9..5198fbb 100644 --- a/docs/examples/Mesh/MeshHex8-FineLayer_paraview.py +++ b/docs/examples/Mesh/MeshHex8-FineLayer_paraview.py @@ -1,44 +1,56 @@ import numpy as np import h5py import GooseFEM as gf # ====================== create fictitious configuration + store to HDF5-file ====================== # create mesh object mesh = gf.Mesh.Hex8.FineLayer(9,17,27) # open data file f = h5py.File('MeshHex8-FineLayer_paraview.hdf5','w') -# write particle positions, and a dummy connectivity -f.create_dataset('/coor',data=mesh.coor()) -f.create_dataset('/conn',data=mesh.conn()) +# element set +elementsMiddleLayer = np.zeros((mesh.nelem()),dtype='int') +elementsMiddleLayer[mesh.elementsMiddleLayer()] = 1 + +# write nodal coordinates and connectivity +f.create_dataset('/coor' ,data=mesh.coor() ) +f.create_dataset('/conn' ,data=mesh.conn() ) +f.create_dataset('/elementsMiddleLayer',data=elementsMiddleLayer) # ======================================== write XDMF-file ========================================= +# basic file-format xmf = ''' MeshHex8-FineLayer_paraview.hdf5:/conn MeshHex8-FineLayer_paraview.hdf5:/coor + + + MeshHex8-FineLayer_paraview.hdf5:/elementsMiddleLayer + + ''' +# write to file, fill mesh dimensions open('MeshHex8-FineLayer_paraview.xmf','w').write(xmf.format( nnode = mesh.nnode(), nelem = mesh.nelem(), nne = mesh.nne(), )) diff --git a/docs/examples/Mesh/MeshHex8-Regular-nodes_paraview.py b/docs/examples/Mesh/MeshHex8-Regular-nodes_paraview.py index c03ef85..403e2f7 100644 --- a/docs/examples/Mesh/MeshHex8-Regular-nodes_paraview.py +++ b/docs/examples/Mesh/MeshHex8-Regular-nodes_paraview.py @@ -1,362 +1,364 @@ import numpy as np import h5py import GooseFEM as gf # ====================== create fictitious configuration + store to HDF5-file ====================== # create mesh object mesh = gf.Mesh.Hex8.Regular(6,8,12) # set node identifier for all boundary sets Z = lambda mesh: np.zeros((mesh.nnode()),dtype='int') nodesFront = Z(mesh); nodesFront [mesh.nodesFront() ] = 1 nodesBack = Z(mesh); nodesBack [mesh.nodesBack() ] = 1 nodesLeft = Z(mesh); nodesLeft [mesh.nodesLeft() ] = 1 nodesRight = Z(mesh); nodesRight [mesh.nodesRight() ] = 1 nodesBottom = Z(mesh); nodesBottom [mesh.nodesBottom() ] = 1 nodesTop = Z(mesh); nodesTop [mesh.nodesTop() ] = 1 nodesFrontFace = Z(mesh); nodesFrontFace [mesh.nodesFrontFace() ] = 1 nodesBackFace = Z(mesh); nodesBackFace [mesh.nodesBackFace() ] = 1 nodesLeftFace = Z(mesh); nodesLeftFace [mesh.nodesLeftFace() ] = 1 nodesRightFace = Z(mesh); nodesRightFace [mesh.nodesRightFace() ] = 1 nodesBottomFace = Z(mesh); nodesBottomFace [mesh.nodesBottomFace() ] = 1 nodesTopFace = Z(mesh); nodesTopFace [mesh.nodesTopFace() ] = 1 nodesFrontBottomEdge = Z(mesh); nodesFrontBottomEdge [mesh.nodesFrontBottomEdge() ] = 1 nodesFrontTopEdge = Z(mesh); nodesFrontTopEdge [mesh.nodesFrontTopEdge() ] = 1 nodesFrontLeftEdge = Z(mesh); nodesFrontLeftEdge [mesh.nodesFrontLeftEdge() ] = 1 nodesFrontRightEdge = Z(mesh); nodesFrontRightEdge [mesh.nodesFrontRightEdge() ] = 1 nodesBackBottomEdge = Z(mesh); nodesBackBottomEdge [mesh.nodesBackBottomEdge() ] = 1 nodesBackTopEdge = Z(mesh); nodesBackTopEdge [mesh.nodesBackTopEdge() ] = 1 nodesBackLeftEdge = Z(mesh); nodesBackLeftEdge [mesh.nodesBackLeftEdge() ] = 1 nodesBackRightEdge = Z(mesh); nodesBackRightEdge [mesh.nodesBackRightEdge() ] = 1 nodesBottomLeftEdge = Z(mesh); nodesBottomLeftEdge [mesh.nodesBottomLeftEdge() ] = 1 nodesBottomRightEdge = Z(mesh); nodesBottomRightEdge [mesh.nodesBottomRightEdge() ] = 1 nodesTopLeftEdge = Z(mesh); nodesTopLeftEdge [mesh.nodesTopLeftEdge() ] = 1 nodesTopRightEdge = Z(mesh); nodesTopRightEdge [mesh.nodesTopRightEdge() ] = 1 nodesFrontBottomOpenEdge = Z(mesh); nodesFrontBottomOpenEdge [mesh.nodesFrontBottomOpenEdge() ] = 1 nodesFrontTopOpenEdge = Z(mesh); nodesFrontTopOpenEdge [mesh.nodesFrontTopOpenEdge() ] = 1 nodesFrontLeftOpenEdge = Z(mesh); nodesFrontLeftOpenEdge [mesh.nodesFrontLeftOpenEdge() ] = 1 nodesFrontRightOpenEdge = Z(mesh); nodesFrontRightOpenEdge [mesh.nodesFrontRightOpenEdge() ] = 1 nodesBackBottomOpenEdge = Z(mesh); nodesBackBottomOpenEdge [mesh.nodesBackBottomOpenEdge() ] = 1 nodesBackTopOpenEdge = Z(mesh); nodesBackTopOpenEdge [mesh.nodesBackTopOpenEdge() ] = 1 nodesBackLeftOpenEdge = Z(mesh); nodesBackLeftOpenEdge [mesh.nodesBackLeftOpenEdge() ] = 1 nodesBackRightOpenEdge = Z(mesh); nodesBackRightOpenEdge [mesh.nodesBackRightOpenEdge() ] = 1 nodesBottomLeftOpenEdge = Z(mesh); nodesBottomLeftOpenEdge [mesh.nodesBottomLeftOpenEdge() ] = 1 nodesBottomRightOpenEdge = Z(mesh); nodesBottomRightOpenEdge [mesh.nodesBottomRightOpenEdge() ] = 1 nodesTopLeftOpenEdge = Z(mesh); nodesTopLeftOpenEdge [mesh.nodesTopLeftOpenEdge() ] = 1 nodesTopRightOpenEdge = Z(mesh); nodesTopRightOpenEdge [mesh.nodesTopRightOpenEdge() ] = 1 nodesFrontBottomLeftCorner = Z(mesh); nodesFrontBottomLeftCorner [mesh.nodesFrontBottomLeftCorner() ] = 1 nodesFrontBottomRightCorner = Z(mesh); nodesFrontBottomRightCorner[mesh.nodesFrontBottomRightCorner()] = 1 nodesFrontTopLeftCorner = Z(mesh); nodesFrontTopLeftCorner [mesh.nodesFrontTopLeftCorner() ] = 1 nodesFrontTopRightCorner = Z(mesh); nodesFrontTopRightCorner [mesh.nodesFrontTopRightCorner() ] = 1 nodesBackBottomLeftCorner = Z(mesh); nodesBackBottomLeftCorner [mesh.nodesBackBottomLeftCorner() ] = 1 nodesBackBottomRightCorner = Z(mesh); nodesBackBottomRightCorner [mesh.nodesBackBottomRightCorner() ] = 1 nodesBackTopLeftCorner = Z(mesh); nodesBackTopLeftCorner [mesh.nodesBackTopLeftCorner() ] = 1 nodesBackTopRightCorner = Z(mesh); nodesBackTopRightCorner [mesh.nodesBackTopRightCorner() ] = 1 # DOFs after eliminating periodicity dofsPeriodic = mesh.dofsPeriodic()[:,0] # open data file f = h5py.File('MeshHex8-Regular-nodes_paraview.hdf5','w') -# write particle positions, and a dummy connectivity +# write nodal positions, and a dummy connectivity f.create_dataset('/coor' ,data=mesh.coor() ) f.create_dataset('/conn' ,data=np.arange(mesh.nnode()) ) f.create_dataset('/nodesFront' ,data=nodesFront ) f.create_dataset('/nodesBack' ,data=nodesBack ) f.create_dataset('/nodesLeft' ,data=nodesLeft ) f.create_dataset('/nodesRight' ,data=nodesRight ) f.create_dataset('/nodesBottom' ,data=nodesBottom ) f.create_dataset('/nodesTop' ,data=nodesTop ) f.create_dataset('/nodesFrontFace' ,data=nodesFrontFace ) f.create_dataset('/nodesBackFace' ,data=nodesBackFace ) f.create_dataset('/nodesLeftFace' ,data=nodesLeftFace ) f.create_dataset('/nodesRightFace' ,data=nodesRightFace ) f.create_dataset('/nodesBottomFace' ,data=nodesBottomFace ) f.create_dataset('/nodesTopFace' ,data=nodesTopFace ) f.create_dataset('/nodesFrontBottomEdge' ,data=nodesFrontBottomEdge ) f.create_dataset('/nodesFrontTopEdge' ,data=nodesFrontTopEdge ) f.create_dataset('/nodesFrontLeftEdge' ,data=nodesFrontLeftEdge ) f.create_dataset('/nodesFrontRightEdge' ,data=nodesFrontRightEdge ) f.create_dataset('/nodesBackBottomEdge' ,data=nodesBackBottomEdge ) f.create_dataset('/nodesBackTopEdge' ,data=nodesBackTopEdge ) f.create_dataset('/nodesBackLeftEdge' ,data=nodesBackLeftEdge ) f.create_dataset('/nodesBackRightEdge' ,data=nodesBackRightEdge ) f.create_dataset('/nodesBottomLeftEdge' ,data=nodesBottomLeftEdge ) f.create_dataset('/nodesBottomRightEdge' ,data=nodesBottomRightEdge ) f.create_dataset('/nodesTopLeftEdge' ,data=nodesTopLeftEdge ) f.create_dataset('/nodesTopRightEdge' ,data=nodesTopRightEdge ) f.create_dataset('/nodesFrontBottomOpenEdge' ,data=nodesFrontBottomOpenEdge ) f.create_dataset('/nodesFrontTopOpenEdge' ,data=nodesFrontTopOpenEdge ) f.create_dataset('/nodesFrontLeftOpenEdge' ,data=nodesFrontLeftOpenEdge ) f.create_dataset('/nodesFrontRightOpenEdge' ,data=nodesFrontRightOpenEdge ) f.create_dataset('/nodesBackBottomOpenEdge' ,data=nodesBackBottomOpenEdge ) f.create_dataset('/nodesBackTopOpenEdge' ,data=nodesBackTopOpenEdge ) f.create_dataset('/nodesBackLeftOpenEdge' ,data=nodesBackLeftOpenEdge ) f.create_dataset('/nodesBackRightOpenEdge' ,data=nodesBackRightOpenEdge ) f.create_dataset('/nodesBottomLeftOpenEdge' ,data=nodesBottomLeftOpenEdge ) f.create_dataset('/nodesBottomRightOpenEdge' ,data=nodesBottomRightOpenEdge ) f.create_dataset('/nodesTopLeftOpenEdge' ,data=nodesTopLeftOpenEdge ) f.create_dataset('/nodesTopRightOpenEdge' ,data=nodesTopRightOpenEdge ) f.create_dataset('/nodesFrontBottomLeftCorner' ,data=nodesFrontBottomLeftCorner ) f.create_dataset('/nodesFrontBottomRightCorner',data=nodesFrontBottomRightCorner) f.create_dataset('/nodesFrontTopLeftCorner' ,data=nodesFrontTopLeftCorner ) f.create_dataset('/nodesFrontTopRightCorner' ,data=nodesFrontTopRightCorner ) f.create_dataset('/nodesBackBottomLeftCorner' ,data=nodesBackBottomLeftCorner ) f.create_dataset('/nodesBackBottomRightCorner' ,data=nodesBackBottomRightCorner ) f.create_dataset('/nodesBackTopLeftCorner' ,data=nodesBackTopLeftCorner ) f.create_dataset('/nodesBackTopRightCorner' ,data=nodesBackTopRightCorner ) f.create_dataset('/dofsPeriodic' ,data=dofsPeriodic ) # ======================================== write XDMF-file ========================================= +# basic file-format xmf = ''' MeshHex8-Regular-nodes_paraview.hdf5:/conn MeshHex8-Regular-nodes_paraview.hdf5:/coor MeshHex8-Regular-nodes_paraview.hdf5:/nodesFront MeshHex8-Regular-nodes_paraview.hdf5:/nodesBack MeshHex8-Regular-nodes_paraview.hdf5:/nodesLeft MeshHex8-Regular-nodes_paraview.hdf5:/nodesRight MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottom MeshHex8-Regular-nodes_paraview.hdf5:/nodesTop MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontFace MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackFace MeshHex8-Regular-nodes_paraview.hdf5:/nodesLeftFace MeshHex8-Regular-nodes_paraview.hdf5:/nodesRightFace MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottomFace MeshHex8-Regular-nodes_paraview.hdf5:/nodesTopFace MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontBottomEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontTopEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontLeftEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontRightEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackBottomEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackTopEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackLeftEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackRightEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottomLeftEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottomRightEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesTopLeftEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesTopRightEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontBottomOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontTopOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontLeftOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontRightOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackBottomOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackTopOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackLeftOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackRightOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottomLeftOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesBottomRightOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesTopLeftOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesTopRightOpenEdge MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontBottomLeftCorner MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontBottomRightCorner MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontTopLeftCorner MeshHex8-Regular-nodes_paraview.hdf5:/nodesFrontTopRightCorner MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackBottomLeftCorner MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackBottomRightCorner MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackTopLeftCorner MeshHex8-Regular-nodes_paraview.hdf5:/nodesBackTopRightCorner MeshHex8-Regular-nodes_paraview.hdf5:/dofsPeriodic ''' +# write to file, fill mesh dimensions open('MeshHex8-Regular-nodes_paraview.xmf','w').write(xmf.format( nnode = mesh.nnode(), )) diff --git a/docs/examples/Mesh/MeshHex8-Regular_paraview.py b/docs/examples/Mesh/MeshHex8-Regular_paraview.py index 3b32502..8a65f2d 100644 --- a/docs/examples/Mesh/MeshHex8-Regular_paraview.py +++ b/docs/examples/Mesh/MeshHex8-Regular_paraview.py @@ -1,44 +1,46 @@ import numpy as np import h5py import GooseFEM as gf # ====================== create fictitious configuration + store to HDF5-file ====================== # create mesh object mesh = gf.Mesh.Hex8.Regular(6,8,12) # open data file f = h5py.File('MeshHex8-Regular_paraview.hdf5','w') -# write particle positions, and a dummy connectivity +# write nodal coordinates and connectivity f.create_dataset('/coor',data=mesh.coor()) f.create_dataset('/conn',data=mesh.conn()) # ======================================== write XDMF-file ========================================= +# basic file-format xmf = ''' MeshHex8-Regular_paraview.hdf5:/conn MeshHex8-Regular_paraview.hdf5:/coor ''' +# write to file, fill mesh dimensions open('MeshHex8-Regular_paraview.xmf','w').write(xmf.format( nnode = mesh.nnode(), nelem = mesh.nelem(), nne = mesh.nne(), )) diff --git a/docs/examples/Mesh/MeshQuad4-FineLayer-nodes_paraview.py b/docs/examples/Mesh/MeshQuad4-FineLayer-nodes_paraview.py new file mode 100644 index 0000000..2398c14 --- /dev/null +++ b/docs/examples/Mesh/MeshQuad4-FineLayer-nodes_paraview.py @@ -0,0 +1,140 @@ + +import numpy as np +import h5py +import GooseFEM as gf + +# ====================== create fictitious configuration + store to HDF5-file ====================== + +# create mesh object +mesh = gf.Mesh.Quad4.FineLayer(9,17) + +# set node identifier for all boundary sets +Z = lambda mesh: np.zeros((mesh.nnode()),dtype='int') +nodesBottomEdge = Z(mesh); nodesBottomEdge [mesh.nodesBottomEdge() ] = 1 +nodesTopEdge = Z(mesh); nodesTopEdge [mesh.nodesTopEdge() ] = 1 +nodesLeftEdge = Z(mesh); nodesLeftEdge [mesh.nodesLeftEdge() ] = 1 +nodesRightEdge = Z(mesh); nodesRightEdge [mesh.nodesRightEdge() ] = 1 +nodesBottomOpenEdge = Z(mesh); nodesBottomOpenEdge [mesh.nodesBottomOpenEdge() ] = 1 +nodesTopOpenEdge = Z(mesh); nodesTopOpenEdge [mesh.nodesTopOpenEdge() ] = 1 +nodesLeftOpenEdge = Z(mesh); nodesLeftOpenEdge [mesh.nodesLeftOpenEdge() ] = 1 +nodesRightOpenEdge = Z(mesh); nodesRightOpenEdge [mesh.nodesRightOpenEdge() ] = 1 +nodesBottomLeftCorner = Z(mesh); nodesBottomLeftCorner [mesh.nodesBottomLeftCorner() ] = 1 +nodesBottomRightCorner = Z(mesh); nodesBottomRightCorner[mesh.nodesBottomRightCorner()] = 1 +nodesTopLeftCorner = Z(mesh); nodesTopLeftCorner [mesh.nodesTopLeftCorner() ] = 1 +nodesTopRightCorner = Z(mesh); nodesTopRightCorner [mesh.nodesTopRightCorner() ] = 1 + +# DOFs after eliminating periodicity +dofsPeriodic = mesh.dofsPeriodic()[:,0] + +# open data file +f = h5py.File('MeshQuad4-FineLayer-nodes_paraview.hdf5','w') + +# write nodal positions, and a dummy connectivity +f.create_dataset('/coor' ,data=mesh.coor() ) +f.create_dataset('/conn' ,data=np.arange(mesh.nnode())) +f.create_dataset('/nodesBottomEdge' ,data=nodesBottomEdge ) +f.create_dataset('/nodesTopEdge' ,data=nodesTopEdge ) +f.create_dataset('/nodesLeftEdge' ,data=nodesLeftEdge ) +f.create_dataset('/nodesRightEdge' ,data=nodesRightEdge ) +f.create_dataset('/nodesBottomOpenEdge' ,data=nodesBottomOpenEdge ) +f.create_dataset('/nodesTopOpenEdge' ,data=nodesTopOpenEdge ) +f.create_dataset('/nodesLeftOpenEdge' ,data=nodesLeftOpenEdge ) +f.create_dataset('/nodesRightOpenEdge' ,data=nodesRightOpenEdge ) +f.create_dataset('/nodesBottomLeftCorner' ,data=nodesBottomLeftCorner ) +f.create_dataset('/nodesBottomRightCorner',data=nodesBottomRightCorner ) +f.create_dataset('/nodesTopLeftCorner' ,data=nodesTopLeftCorner ) +f.create_dataset('/nodesTopRightCorner' ,data=nodesTopRightCorner ) +f.create_dataset('/dofsPeriodic' ,data=dofsPeriodic ) + +# ======================================== write XDMF-file ========================================= + +# basic file-format +xmf = ''' + + + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/conn + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/coor + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesBottomEdge + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesTopEdge + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesLeftEdge + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesRightEdge + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesBottomOpenEdge + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesTopOpenEdge + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesLeftOpenEdge + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesRightOpenEdge + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesBottomLeftCorner + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesBottomRightCorner + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesTopLeftCorner + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/nodesTopRightCorner + + + + + MeshQuad4-FineLayer-nodes_paraview.hdf5:/dofsPeriodic + + + + + +''' + +# write to file, fill mesh dimensions +open('MeshQuad4-FineLayer-nodes_paraview.xmf','w').write(xmf.format( + nnode = mesh.nnode(), +)) diff --git a/docs/examples/Mesh/MeshQuad4-FineLayer_paraview.py b/docs/examples/Mesh/MeshQuad4-FineLayer_paraview.py index 3a027ef..4ec21e1 100644 --- a/docs/examples/Mesh/MeshQuad4-FineLayer_paraview.py +++ b/docs/examples/Mesh/MeshQuad4-FineLayer_paraview.py @@ -1,115 +1,56 @@ -# ========================================= load libraries ========================================= - -import GooseFEM as gf import numpy as np import h5py +import GooseFEM as gf -# ========================================== define mesh =========================================== - -mesh = gf.Mesh.Quad4.FineLayer(36,36) -coor = mesh.coor() -conn = mesh.conn() - -# ======================================== write HDF5=file ========================================= +# ====================== create fictitious configuration + store to HDF5-file ====================== -# set file-name base -name = 'MeshQuad4-FineLayer_paraview' +# create mesh object +mesh = gf.Mesh.Quad4.FineLayer(9,17) # open data file -f = h5py.File('{name:s}.hdf5'.format(name=name),'w') +f = h5py.File('MeshQuad4-FineLayer_paraview.hdf5','w') -# write particle positions, and a dummy connectivity -f.create_dataset('/coor',data=coor) -f.create_dataset('/conn',data=conn) +# element set +elementsMiddleLayer = np.zeros((mesh.nelem()),dtype='int') +elementsMiddleLayer[mesh.elementsMiddleLayer()] = 1 -# create a sample deformation: simple shear -for inc,gamma in enumerate(np.linspace(0,1,100)): - - # - initialize displacement (should be 3-d for ParaView) - U = np.zeros((mesh.nnode(),3),dtype='float64') - # - set - U[:,0] += gamma * coor[:,1] - # - store - f.create_dataset('/disp/%d'%inc,data=U) +# write nodal coordinates and connectivity +f.create_dataset('/coor' ,data=mesh.coor() ) +f.create_dataset('/conn' ,data=mesh.conn() ) +f.create_dataset('/elementsMiddleLayer',data=elementsMiddleLayer) # ======================================== write XDMF-file ========================================= -# --------------------------------- format of the main structure ---------------------------------- - +# basic file-format xmf = ''' - -{series:s} - + + + + MeshQuad4-FineLayer_paraview.hdf5:/conn + + + + + MeshQuad4-FineLayer_paraview.hdf5:/coor + + + + + MeshQuad4-FineLayer_paraview.hdf5:/elementsMiddleLayer + + + ''' -# ----------------------------- format of an increment in time-series ------------------------------ - -grid = ''' - -''' - -# ------------------------------------------- write file ------------------------------------------- - -# initialize string that will contain the full time series -txt = '' - -# loop over all increments, append the time series -for inc in range(100): - txt += grid.format( - name = name, - inc = inc, - nnode = mesh.nnode(), - ndim = mesh.ndim(), - nelem = mesh.nelem(), - nne = mesh.nne() - ) - -# write xmf-file, fix the indentation -fname = '{name:s}.xmf'.format(name=name) -open(fname,'w').write(xmf.format(series=' '+txt.replace('\n','\n '))) - -# ======================================== print node=sets ========================================= - -np.set_printoptions(linewidth=200) - -print('nodesBottomEdge =', mesh.nodesBottomEdge()) -print('nodesTopEdge =', mesh.nodesTopEdge()) -print('nodesLeftEdge =', mesh.nodesLeftEdge()) -print('nodesRightEdge =', mesh.nodesRightEdge()) -print('nodesLeftBottomCorner =', mesh.nodesLeftBottomCorner()) -print('nodesRightBottomCorner =', mesh.nodesRightBottomCorner()) -print('nodesRightTopCorner =', mesh.nodesRightTopCorner()) -print('nodesLeftTopCorner =', mesh.nodesLeftTopCorner()) - -print('nodesPeriodic =') -print(mesh.nodesPeriodic()) -print('nodesOrigin =') -print(mesh.nodesOrigin()) - -# check uniqueness of tied nodes -ties = mesh.nodesPeriodic()[:,1] - -if len(np.unique(ties)) != len(ties): - raise IOError('Nodal ties not unique') +# write to file, fill mesh dimensions +open('MeshQuad4-FineLayer_paraview.xmf','w').write(xmf.format( + nnode = mesh.nnode(), + nelem = mesh.nelem(), + nne = mesh.nne(), +)) diff --git a/docs/examples/Mesh/MeshQuad4-Regular-nodes_paraview.py b/docs/examples/Mesh/MeshQuad4-Regular-nodes_paraview.py new file mode 100644 index 0000000..ed1e6d6 --- /dev/null +++ b/docs/examples/Mesh/MeshQuad4-Regular-nodes_paraview.py @@ -0,0 +1,140 @@ + +import numpy as np +import h5py +import GooseFEM as gf + +# ====================== create fictitious configuration + store to HDF5-file ====================== + +# create mesh object +mesh = gf.Mesh.Quad4.Regular(9,11) + +# set node identifier for all boundary sets +Z = lambda mesh: np.zeros((mesh.nnode()),dtype='int') +nodesBottomEdge = Z(mesh); nodesBottomEdge [mesh.nodesBottomEdge() ] = 1 +nodesTopEdge = Z(mesh); nodesTopEdge [mesh.nodesTopEdge() ] = 1 +nodesLeftEdge = Z(mesh); nodesLeftEdge [mesh.nodesLeftEdge() ] = 1 +nodesRightEdge = Z(mesh); nodesRightEdge [mesh.nodesRightEdge() ] = 1 +nodesBottomOpenEdge = Z(mesh); nodesBottomOpenEdge [mesh.nodesBottomOpenEdge() ] = 1 +nodesTopOpenEdge = Z(mesh); nodesTopOpenEdge [mesh.nodesTopOpenEdge() ] = 1 +nodesLeftOpenEdge = Z(mesh); nodesLeftOpenEdge [mesh.nodesLeftOpenEdge() ] = 1 +nodesRightOpenEdge = Z(mesh); nodesRightOpenEdge [mesh.nodesRightOpenEdge() ] = 1 +nodesBottomLeftCorner = Z(mesh); nodesBottomLeftCorner [mesh.nodesBottomLeftCorner() ] = 1 +nodesBottomRightCorner = Z(mesh); nodesBottomRightCorner[mesh.nodesBottomRightCorner()] = 1 +nodesTopLeftCorner = Z(mesh); nodesTopLeftCorner [mesh.nodesTopLeftCorner() ] = 1 +nodesTopRightCorner = Z(mesh); nodesTopRightCorner [mesh.nodesTopRightCorner() ] = 1 + +# DOFs after eliminating periodicity +dofsPeriodic = mesh.dofsPeriodic()[:,0] + +# open data file +f = h5py.File('MeshQuad4-Regular-nodes_paraview.hdf5','w') + +# write nodal positions, and a dummy connectivity +f.create_dataset('/coor' ,data=mesh.coor() ) +f.create_dataset('/conn' ,data=np.arange(mesh.nnode())) +f.create_dataset('/nodesBottomEdge' ,data=nodesBottomEdge ) +f.create_dataset('/nodesTopEdge' ,data=nodesTopEdge ) +f.create_dataset('/nodesLeftEdge' ,data=nodesLeftEdge ) +f.create_dataset('/nodesRightEdge' ,data=nodesRightEdge ) +f.create_dataset('/nodesBottomOpenEdge' ,data=nodesBottomOpenEdge ) +f.create_dataset('/nodesTopOpenEdge' ,data=nodesTopOpenEdge ) +f.create_dataset('/nodesLeftOpenEdge' ,data=nodesLeftOpenEdge ) +f.create_dataset('/nodesRightOpenEdge' ,data=nodesRightOpenEdge ) +f.create_dataset('/nodesBottomLeftCorner' ,data=nodesBottomLeftCorner ) +f.create_dataset('/nodesBottomRightCorner',data=nodesBottomRightCorner ) +f.create_dataset('/nodesTopLeftCorner' ,data=nodesTopLeftCorner ) +f.create_dataset('/nodesTopRightCorner' ,data=nodesTopRightCorner ) +f.create_dataset('/dofsPeriodic' ,data=dofsPeriodic ) + +# ======================================== write XDMF-file ========================================= + +# basic file-format +xmf = ''' + + + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/conn + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/coor + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesBottomEdge + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesTopEdge + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesLeftEdge + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesRightEdge + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesBottomOpenEdge + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesTopOpenEdge + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesLeftOpenEdge + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesRightOpenEdge + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesBottomLeftCorner + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesBottomRightCorner + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesTopLeftCorner + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/nodesTopRightCorner + + + + + MeshQuad4-Regular-nodes_paraview.hdf5:/dofsPeriodic + + + + + +''' + +# write to file, fill mesh dimensions +open('MeshQuad4-Regular-nodes_paraview.xmf','w').write(xmf.format( + nnode = mesh.nnode(), +)) diff --git a/docs/examples/Mesh/MeshQuad4-Regular_paraview.py b/docs/examples/Mesh/MeshQuad4-Regular_paraview.py index 78ba252..3c58990 100644 --- a/docs/examples/Mesh/MeshQuad4-Regular_paraview.py +++ b/docs/examples/Mesh/MeshQuad4-Regular_paraview.py @@ -1,115 +1,46 @@ -# ========================================= load libraries ========================================= - -import GooseFEM as gf import numpy as np import h5py +import GooseFEM as gf -# ========================================== define mesh =========================================== - -mesh = gf.Mesh.Quad4.Regular(4,5) -coor = mesh.coor() -conn = mesh.conn() - -# ======================================== write HDF5=file ========================================= +# ====================== create fictitious configuration + store to HDF5-file ====================== -# set file-name base -name = 'MeshQuad4-Regular_paraview' +# create mesh object +mesh = gf.Mesh.Quad4.Regular(9,11) # open data file -f = h5py.File('{name:s}.hdf5'.format(name=name),'w') +f = h5py.File('MeshQuad4-Regular_paraview.hdf5','w') -# write particle positions, and a dummy connectivity -f.create_dataset('/coor',data=coor) -f.create_dataset('/conn',data=conn) - -# create a sample deformation: simple shear -for inc,gamma in enumerate(np.linspace(0,1,100)): - - # - initialize displacement (should be 3-d for ParaView) - U = np.zeros((mesh.nnode(),3),dtype='float64') - # - set - U[:,0] += gamma * coor[:,1] - # - store - f.create_dataset('/disp/%d'%inc,data=U) +# write nodal coordinates and connectivity +f.create_dataset('/coor',data=mesh.coor()) +f.create_dataset('/conn',data=mesh.conn()) # ======================================== write XDMF-file ========================================= -# --------------------------------- format of the main structure ---------------------------------- - +# basic file-format xmf = ''' - -{series:s} - + + + + MeshQuad4-Regular_paraview.hdf5:/conn + + + + + MeshQuad4-Regular_paraview.hdf5:/coor + + + ''' -# ----------------------------- format of an increment in time-series ------------------------------ - -grid = ''' - -''' - -# ------------------------------------------- write file ------------------------------------------- - -# initialize string that will contain the full time series -txt = '' - -# loop over all increments, append the time series -for inc in range(100): - txt += grid.format( - name = name, - inc = inc, - nnode = mesh.nnode(), - ndim = mesh.ndim(), - nelem = mesh.nelem(), - nne = mesh.nne() - ) - -# write xmf-file, fix the indentation -fname = '{name:s}.xmf'.format(name=name) -open(fname,'w').write(xmf.format(series=' '+txt.replace('\n','\n '))) - -# ======================================== print node=sets ========================================= - -np.set_printoptions(linewidth=200) - -print('nodesBottomEdge =', mesh.nodesBottomEdge()) -print('nodesTopEdge =', mesh.nodesTopEdge()) -print('nodesLeftEdge =', mesh.nodesLeftEdge()) -print('nodesRightEdge =', mesh.nodesRightEdge()) -print('nodesLeftBottomCorner =', mesh.nodesLeftBottomCorner()) -print('nodesRightBottomCorner =', mesh.nodesRightBottomCorner()) -print('nodesRightTopCorner =', mesh.nodesRightTopCorner()) -print('nodesLeftTopCorner =', mesh.nodesLeftTopCorner()) - -print('nodesPeriodic =') -print(mesh.nodesPeriodic()) -print('nodesOrigin =') -print(mesh.nodesOrigin()) - -# check uniqueness of tied nodes -ties = mesh.nodesPeriodic()[:,1] - -if len(np.unique(ties)) != len(ties): - raise IOError('Nodal ties not unique') +# write to file, fill mesh dimensions +open('MeshQuad4-Regular_paraview.xmf','w').write(xmf.format( + nnode = mesh.nnode(), + nelem = mesh.nelem(), + nne = mesh.nne(), +)) diff --git a/docs/examples/Mesh/MeshTri3-Regular-nodes_paraview.py b/docs/examples/Mesh/MeshTri3-Regular-nodes_paraview.py new file mode 100644 index 0000000..b320616 --- /dev/null +++ b/docs/examples/Mesh/MeshTri3-Regular-nodes_paraview.py @@ -0,0 +1,140 @@ + +import numpy as np +import h5py +import GooseFEM as gf + +# ====================== create fictitious configuration + store to HDF5-file ====================== + +# create mesh object +mesh = gf.Mesh.Tri3.Regular(9,11) + +# set node identifier for all boundary sets +Z = lambda mesh: np.zeros((mesh.nnode()),dtype='int') +nodesBottomEdge = Z(mesh); nodesBottomEdge [mesh.nodesBottomEdge() ] = 1 +nodesTopEdge = Z(mesh); nodesTopEdge [mesh.nodesTopEdge() ] = 1 +nodesLeftEdge = Z(mesh); nodesLeftEdge [mesh.nodesLeftEdge() ] = 1 +nodesRightEdge = Z(mesh); nodesRightEdge [mesh.nodesRightEdge() ] = 1 +nodesBottomOpenEdge = Z(mesh); nodesBottomOpenEdge [mesh.nodesBottomOpenEdge() ] = 1 +nodesTopOpenEdge = Z(mesh); nodesTopOpenEdge [mesh.nodesTopOpenEdge() ] = 1 +nodesLeftOpenEdge = Z(mesh); nodesLeftOpenEdge [mesh.nodesLeftOpenEdge() ] = 1 +nodesRightOpenEdge = Z(mesh); nodesRightOpenEdge [mesh.nodesRightOpenEdge() ] = 1 +nodesBottomLeftCorner = Z(mesh); nodesBottomLeftCorner [mesh.nodesBottomLeftCorner() ] = 1 +nodesBottomRightCorner = Z(mesh); nodesBottomRightCorner[mesh.nodesBottomRightCorner()] = 1 +nodesTopLeftCorner = Z(mesh); nodesTopLeftCorner [mesh.nodesTopLeftCorner() ] = 1 +nodesTopRightCorner = Z(mesh); nodesTopRightCorner [mesh.nodesTopRightCorner() ] = 1 + +# DOFs after eliminating periodicity +dofsPeriodic = mesh.dofsPeriodic()[:,0] + +# open data file +f = h5py.File('MeshTri3-Regular-nodes_paraview.hdf5','w') + +# write nodal positions, and a dummy connectivity +f.create_dataset('/coor' ,data=mesh.coor() ) +f.create_dataset('/conn' ,data=np.arange(mesh.nnode())) +f.create_dataset('/nodesBottomEdge' ,data=nodesBottomEdge ) +f.create_dataset('/nodesTopEdge' ,data=nodesTopEdge ) +f.create_dataset('/nodesLeftEdge' ,data=nodesLeftEdge ) +f.create_dataset('/nodesRightEdge' ,data=nodesRightEdge ) +f.create_dataset('/nodesBottomOpenEdge' ,data=nodesBottomOpenEdge ) +f.create_dataset('/nodesTopOpenEdge' ,data=nodesTopOpenEdge ) +f.create_dataset('/nodesLeftOpenEdge' ,data=nodesLeftOpenEdge ) +f.create_dataset('/nodesRightOpenEdge' ,data=nodesRightOpenEdge ) +f.create_dataset('/nodesBottomLeftCorner' ,data=nodesBottomLeftCorner ) +f.create_dataset('/nodesBottomRightCorner',data=nodesBottomRightCorner ) +f.create_dataset('/nodesTopLeftCorner' ,data=nodesTopLeftCorner ) +f.create_dataset('/nodesTopRightCorner' ,data=nodesTopRightCorner ) +f.create_dataset('/dofsPeriodic' ,data=dofsPeriodic ) + +# ======================================== write XDMF-file ========================================= + +# basic file-format +xmf = ''' + + + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/conn + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/coor + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesBottomEdge + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesTopEdge + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesLeftEdge + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesRightEdge + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesBottomOpenEdge + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesTopOpenEdge + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesLeftOpenEdge + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesRightOpenEdge + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesBottomLeftCorner + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesBottomRightCorner + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesTopLeftCorner + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/nodesTopRightCorner + + + + + MeshTri3-Regular-nodes_paraview.hdf5:/dofsPeriodic + + + + + +''' + +# write to file, fill mesh dimensions +open('MeshTri3-Regular-nodes_paraview.xmf','w').write(xmf.format( + nnode = mesh.nnode(), +)) diff --git a/docs/examples/Mesh/MeshTri3-Regular_paraview.py b/docs/examples/Mesh/MeshTri3-Regular_paraview.py index 1e95a2f..66cc261 100644 --- a/docs/examples/Mesh/MeshTri3-Regular_paraview.py +++ b/docs/examples/Mesh/MeshTri3-Regular_paraview.py @@ -1,115 +1,46 @@ -# ========================================= load libraries ========================================= - -import GooseFEM as gf import numpy as np import h5py +import GooseFEM as gf -# ========================================== define mesh =========================================== - -mesh = gf.Mesh.Tri3.Regular(4,5) -coor = mesh.coor() -conn = mesh.conn() - -# ======================================== write HDF5=file ========================================= +# ====================== create fictitious configuration + store to HDF5-file ====================== -# set file-name base -name = 'MeshTri3-Regular_paraview' +# create mesh object +mesh = gf.Mesh.Tri3.Regular(9,11) # open data file -f = h5py.File('{name:s}.hdf5'.format(name=name),'w') +f = h5py.File('MeshTri3-Regular_paraview.hdf5','w') -# write particle positions, and a dummy connectivity -f.create_dataset('/coor',data=coor) -f.create_dataset('/conn',data=conn) - -# create a sample deformation: simple shear -for inc,gamma in enumerate(np.linspace(0,1,100)): - - # - initialize displacement (should be 3-d for ParaView) - U = np.zeros((mesh.nnode(),3),dtype='float64') - # - set - U[:,0] += gamma * coor[:,1] - # - store - f.create_dataset('/disp/%d'%inc,data=U) +# write nodal coordinates and connectivity +f.create_dataset('/coor',data=mesh.coor()) +f.create_dataset('/conn',data=mesh.conn()) # ======================================== write XDMF-file ========================================= -# --------------------------------- format of the main structure ---------------------------------- - +# basic file-format xmf = ''' - -{series:s} - + + + + MeshTri3-Regular_paraview.hdf5:/conn + + + + + MeshTri3-Regular_paraview.hdf5:/coor + + + ''' -# ----------------------------- format of an increment in time-series ------------------------------ - -grid = ''' - -''' - -# ------------------------------------------- write file ------------------------------------------- - -# initialize string that will contain the full time series -txt = '' - -# loop over all increments, append the time series -for inc in range(100): - txt += grid.format( - name = name, - inc = inc, - nnode = mesh.nnode(), - ndim = mesh.ndim(), - nelem = mesh.nelem(), - nne = mesh.nne() - ) - -# write xmf-file, fix the indentation -fname = '{name:s}.xmf'.format(name=name) -open(fname,'w').write(xmf.format(series=' '+txt.replace('\n','\n '))) - -# ======================================== print node=sets ========================================= - -np.set_printoptions(linewidth=200) - -print('nodesBottomEdge =', mesh.nodesBottomEdge()) -print('nodesTopEdge =', mesh.nodesTopEdge()) -print('nodesLeftEdge =', mesh.nodesLeftEdge()) -print('nodesRightEdge =', mesh.nodesRightEdge()) -print('nodesLeftBottomCorner =', mesh.nodesLeftBottomCorner()) -print('nodesRightBottomCorner =', mesh.nodesRightBottomCorner()) -print('nodesRightTopCorner =', mesh.nodesRightTopCorner()) -print('nodesLeftTopCorner =', mesh.nodesLeftTopCorner()) - -print('nodesPeriodic =') -print(mesh.nodesPeriodic()) -print('nodesOrigin =') -print(mesh.nodesOrigin()) - -# check uniqueness of tied nodes -ties = mesh.nodesPeriodic()[:,1] - -if len(np.unique(ties)) != len(ties): - raise IOError('Nodal ties not unique') +# write to file, fill mesh dimensions +open('MeshTri3-Regular_paraview.xmf','w').write(xmf.format( + nnode = mesh.nnode(), + nelem = mesh.nelem(), + nne = mesh.nne(), +)) diff --git a/python/python_interface.cpp b/python/python_interface.cpp index 3bc3025..a01ceea 100644 --- a/python/python_interface.cpp +++ b/python/python_interface.cpp @@ -1,492 +1,505 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include #include "../src/GooseFEM/GooseFEM.h" // alias for short-hand notation below namespace py = pybind11; PYBIND11_MODULE(GooseFEM, m) { // ================================================================================================= m.doc() = "Some simple finite element meshes and operations"; // define submodules "mXXX" py::module mMesh = m .def_submodule("Mesh" , "Generic mesh routines" ); py::module mMeshTri3 = mMesh.def_submodule("Tri3" , "Linear triangular elements (2D)" ); py::module mMeshQuad4 = mMesh.def_submodule("Quad4" , "Linear quadrilateral elements (2D)" ); py::module mMeshHex8 = mMesh.def_submodule("Hex8" , "Linear hexahedron (brick) elements (3D)"); // ======================================= GooseFEM/Mesh.h ======================================== mMesh.def("elem2node", &GooseFEM::Mesh::elem2node, "Elements connect to each node: [ number of elements , element numbers ]", py::arg("conn") ); mMesh.def("dofs", &GooseFEM::Mesh::dofs, "List with DOF-numbers (in sequential order)", py::arg("nnode"), py::arg("ndim") ); mMesh.def("renumber", py::overload_cast(&GooseFEM::Mesh::renumber), "Renumber DOF-list to use the lowest possible index", py::arg("dofs") ); mMesh.def("renumber", py::overload_cast(&GooseFEM::Mesh::renumber), "Renumber DOF-list to begin or end with 'idx'", py::arg("dofs"), py::arg("idx"), py::arg("location")="end" ); // ====================================== GooseFEM/MeshHex8.h ====================================== py::class_(mMeshHex8, "Regular") // constructor .def( py::init(), "mesh with nx*ny*nz 'pixels' and edge size h", py::arg("nx"), py::arg("ny"), py::arg("nz"), py::arg("h")=1. ) // sizes .def("nelem" , &GooseFEM::Mesh::Hex8::Regular::nelem ) .def("nnode" , &GooseFEM::Mesh::Hex8::Regular::nnode ) .def("nne" , &GooseFEM::Mesh::Hex8::Regular::nne ) .def("ndim" , &GooseFEM::Mesh::Hex8::Regular::ndim ) // mesh .def("coor" , &GooseFEM::Mesh::Hex8::Regular::coor ) .def("conn" , &GooseFEM::Mesh::Hex8::Regular::conn ) // boundary nodes: planes .def("nodesFront" , &GooseFEM::Mesh::Hex8::Regular::nodesFront ) .def("nodesBack" , &GooseFEM::Mesh::Hex8::Regular::nodesBack ) .def("nodesLeft" , &GooseFEM::Mesh::Hex8::Regular::nodesLeft ) .def("nodesRight" , &GooseFEM::Mesh::Hex8::Regular::nodesRight ) .def("nodesBottom" , &GooseFEM::Mesh::Hex8::Regular::nodesBottom ) .def("nodesTop" , &GooseFEM::Mesh::Hex8::Regular::nodesTop ) // boundary nodes: faces .def("nodesFrontFace" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontFace ) .def("nodesBackFace" , &GooseFEM::Mesh::Hex8::Regular::nodesBackFace ) .def("nodesLeftFace" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFace ) .def("nodesRightFace" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFace ) .def("nodesBottomFace" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFace ) .def("nodesTopFace" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFace ) // boundary nodes: edges .def("nodesFrontBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomEdge ) .def("nodesFrontTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopEdge ) .def("nodesFrontLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftEdge ) .def("nodesFrontRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightEdge ) .def("nodesBackBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomEdge ) .def("nodesBackTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopEdge ) .def("nodesBackLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftEdge ) .def("nodesBackRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightEdge ) .def("nodesBottomLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftEdge ) .def("nodesBottomRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightEdge ) .def("nodesTopLeftEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftEdge ) .def("nodesTopRightEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightEdge ) // boundary nodes: faces (aliases) .def("nodesBottomFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontEdge ) .def("nodesBottomBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackEdge ) .def("nodesTopFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontEdge ) .def("nodesTopBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackEdge ) .def("nodesLeftBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomEdge ) .def("nodesLeftFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontEdge ) .def("nodesLeftBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackEdge ) .def("nodesLeftTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopEdge ) .def("nodesRightBottomEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomEdge ) .def("nodesRightTopEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopEdge ) .def("nodesRightFrontEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontEdge ) .def("nodesRightBackEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackEdge ) // boundary nodes: edges, without corners .def("nodesFrontBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomOpenEdge ) .def("nodesFrontTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopOpenEdge ) .def("nodesFrontLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftOpenEdge ) .def("nodesFrontRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightOpenEdge ) .def("nodesBackBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomOpenEdge ) .def("nodesBackTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopOpenEdge ) .def("nodesBackLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftOpenEdge ) .def("nodesBackRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightOpenEdge ) .def("nodesBottomLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftOpenEdge ) .def("nodesBottomRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightOpenEdge ) .def("nodesTopLeftOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftOpenEdge ) .def("nodesTopRightOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightOpenEdge ) // boundary nodes: edges, without corners (aliases) .def("nodesBottomFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontOpenEdge ) .def("nodesBottomBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackOpenEdge ) .def("nodesTopFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontOpenEdge ) .def("nodesTopBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackOpenEdge ) .def("nodesLeftBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomOpenEdge ) .def("nodesLeftFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontOpenEdge ) .def("nodesLeftBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackOpenEdge ) .def("nodesLeftTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopOpenEdge ) .def("nodesRightBottomOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomOpenEdge ) .def("nodesRightTopOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopOpenEdge ) .def("nodesRightFrontOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontOpenEdge ) .def("nodesRightBackOpenEdge" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackOpenEdge ) // boundary nodes: corners .def("nodesFrontBottomLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomLeftCorner ) .def("nodesFrontBottomRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomRightCorner) .def("nodesFrontTopLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopLeftCorner ) .def("nodesFrontTopRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopRightCorner ) .def("nodesBackBottomLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomLeftCorner ) .def("nodesBackBottomRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomRightCorner ) .def("nodesBackTopLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopLeftCorner ) .def("nodesBackTopRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackTopRightCorner ) // boundary nodes: corners (aliases) .def("nodesFrontLeftBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftBottomCorner ) .def("nodesBottomFrontLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontLeftCorner ) .def("nodesBottomLeftFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftFrontCorner ) .def("nodesLeftFrontBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontBottomCorner ) .def("nodesLeftBottomFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomFrontCorner ) .def("nodesFrontRightBottomCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightBottomCorner) .def("nodesBottomFrontRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontRightCorner) .def("nodesBottomRightFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightFrontCorner) .def("nodesRightFrontBottomCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontBottomCorner) .def("nodesRightBottomFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomFrontCorner) .def("nodesFrontLeftTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftTopCorner ) .def("nodesTopFrontLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontLeftCorner ) .def("nodesTopLeftFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftFrontCorner ) .def("nodesLeftFrontTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontTopCorner ) .def("nodesLeftTopFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopFrontCorner ) .def("nodesFrontRightTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightTopCorner ) .def("nodesTopFrontRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontRightCorner ) .def("nodesTopRightFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightFrontCorner ) .def("nodesRightFrontTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontTopCorner ) .def("nodesRightTopFrontCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopFrontCorner ) .def("nodesBackLeftBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftBottomCorner ) .def("nodesBottomBackLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackLeftCorner ) .def("nodesBottomLeftBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftBackCorner ) .def("nodesLeftBackBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackBottomCorner ) .def("nodesLeftBottomBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomBackCorner ) .def("nodesBackRightBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightBottomCorner ) .def("nodesBottomBackRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackRightCorner ) .def("nodesBottomRightBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightBackCorner ) .def("nodesRightBackBottomCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackBottomCorner ) .def("nodesRightBottomBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomBackCorner ) .def("nodesBackLeftTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftTopCorner ) .def("nodesTopBackLeftCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackLeftCorner ) .def("nodesTopLeftBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftBackCorner ) .def("nodesLeftBackTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackTopCorner ) .def("nodesLeftTopBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopBackCorner ) .def("nodesBackRightTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesBackRightTopCorner ) .def("nodesTopBackRightCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopBackRightCorner ) .def("nodesTopRightBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesTopRightBackCorner ) .def("nodesRightBackTopCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightBackTopCorner ) .def("nodesRightTopBackCorner" , &GooseFEM::Mesh::Hex8::Regular::nodesRightTopBackCorner ) // periodicity .def("nodesPeriodic" , &GooseFEM::Mesh::Hex8::Regular::nodesPeriodic ) .def("nodesOrigin" , &GooseFEM::Mesh::Hex8::Regular::nodesOrigin ) .def("dofs" , &GooseFEM::Mesh::Hex8::Regular::dofs ) .def("dofsPeriodic" , &GooseFEM::Mesh::Hex8::Regular::dofsPeriodic ) // print to screen .def("__repr__", [](const GooseFEM::Mesh::Hex8::Regular &a){ return ""; } ); // ------------------------------------------------------------------------------------------------- py::class_(mMeshHex8, "FineLayer") // constructor .def( py::init(), "mesh with nx*ny*nz 'pixels' and edge size h", py::arg("nx"), py::arg("ny"), py::arg("nz"), py::arg("h")=1., py::arg("nfine")=1 ) // sizes .def("nelem" , &GooseFEM::Mesh::Hex8::FineLayer::nelem ) .def("nnode" , &GooseFEM::Mesh::Hex8::FineLayer::nnode ) .def("nne" , &GooseFEM::Mesh::Hex8::FineLayer::nne ) .def("ndim" , &GooseFEM::Mesh::Hex8::FineLayer::ndim ) .def("shape" , &GooseFEM::Mesh::Hex8::FineLayer::shape ) // mesh .def("coor" , &GooseFEM::Mesh::Hex8::FineLayer::coor ) .def("conn" , &GooseFEM::Mesh::Hex8::FineLayer::conn ) + // element sets + .def("elementsMiddleLayer" , &GooseFEM::Mesh::Hex8::FineLayer::elementsMiddleLayer ) // boundary nodes: planes .def("nodesFront" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFront ) .def("nodesBack" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBack ) .def("nodesLeft" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeft ) .def("nodesRight" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRight ) .def("nodesBottom" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottom ) .def("nodesTop" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTop ) // boundary nodes: faces .def("nodesFrontFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontFace ) .def("nodesBackFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackFace ) .def("nodesLeftFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFace ) .def("nodesRightFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFace ) .def("nodesBottomFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFace ) .def("nodesTopFace" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFace ) // boundary nodes: edges .def("nodesFrontBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomEdge ) .def("nodesFrontTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopEdge ) .def("nodesFrontLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftEdge ) .def("nodesFrontRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightEdge ) .def("nodesBackBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomEdge ) .def("nodesBackTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopEdge ) .def("nodesBackLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftEdge ) .def("nodesBackRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightEdge ) .def("nodesBottomLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftEdge ) .def("nodesBottomRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightEdge ) .def("nodesTopLeftEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftEdge ) .def("nodesTopRightEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightEdge ) // boundary nodes: faces (aliases) .def("nodesBottomFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontEdge ) .def("nodesBottomBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackEdge ) .def("nodesTopFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontEdge ) .def("nodesTopBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackEdge ) .def("nodesLeftBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomEdge ) .def("nodesLeftFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontEdge ) .def("nodesLeftBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackEdge ) .def("nodesLeftTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopEdge ) .def("nodesRightBottomEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomEdge ) .def("nodesRightTopEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopEdge ) .def("nodesRightFrontEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontEdge ) .def("nodesRightBackEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackEdge ) // boundary nodes: edges, without corners .def("nodesFrontBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomOpenEdge ) .def("nodesFrontTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopOpenEdge ) .def("nodesFrontLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftOpenEdge ) .def("nodesFrontRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightOpenEdge ) .def("nodesBackBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomOpenEdge ) .def("nodesBackTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopOpenEdge ) .def("nodesBackLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftOpenEdge ) .def("nodesBackRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightOpenEdge ) .def("nodesBottomLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftOpenEdge ) .def("nodesBottomRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightOpenEdge ) .def("nodesTopLeftOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftOpenEdge ) .def("nodesTopRightOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightOpenEdge ) // boundary nodes: edges, without corners (aliases) .def("nodesBottomFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontOpenEdge ) .def("nodesBottomBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackOpenEdge ) .def("nodesTopFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontOpenEdge ) .def("nodesTopBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackOpenEdge ) .def("nodesLeftBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomOpenEdge ) .def("nodesLeftFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontOpenEdge ) .def("nodesLeftBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackOpenEdge ) .def("nodesLeftTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopOpenEdge ) .def("nodesRightBottomOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomOpenEdge ) .def("nodesRightTopOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopOpenEdge ) .def("nodesRightFrontOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontOpenEdge ) .def("nodesRightBackOpenEdge" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackOpenEdge ) // boundary nodes: corners .def("nodesFrontBottomLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomLeftCorner ) .def("nodesFrontBottomRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomRightCorner) .def("nodesFrontTopLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopLeftCorner ) .def("nodesFrontTopRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopRightCorner ) .def("nodesBackBottomLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomLeftCorner ) .def("nodesBackBottomRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomRightCorner ) .def("nodesBackTopLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopLeftCorner ) .def("nodesBackTopRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopRightCorner ) // boundary nodes: corners (aliases) .def("nodesFrontLeftBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftBottomCorner ) .def("nodesBottomFrontLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontLeftCorner ) .def("nodesBottomLeftFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftFrontCorner ) .def("nodesLeftFrontBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontBottomCorner ) .def("nodesLeftBottomFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomFrontCorner ) .def("nodesFrontRightBottomCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightBottomCorner) .def("nodesBottomFrontRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontRightCorner) .def("nodesBottomRightFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightFrontCorner) .def("nodesRightFrontBottomCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontBottomCorner) .def("nodesRightBottomFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomFrontCorner) .def("nodesFrontLeftTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftTopCorner ) .def("nodesTopFrontLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontLeftCorner ) .def("nodesTopLeftFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftFrontCorner ) .def("nodesLeftFrontTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontTopCorner ) .def("nodesLeftTopFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopFrontCorner ) .def("nodesFrontRightTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightTopCorner ) .def("nodesTopFrontRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontRightCorner ) .def("nodesTopRightFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightFrontCorner ) .def("nodesRightFrontTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontTopCorner ) .def("nodesRightTopFrontCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopFrontCorner ) .def("nodesBackLeftBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftBottomCorner ) .def("nodesBottomBackLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackLeftCorner ) .def("nodesBottomLeftBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftBackCorner ) .def("nodesLeftBackBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackBottomCorner ) .def("nodesLeftBottomBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomBackCorner ) .def("nodesBackRightBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightBottomCorner ) .def("nodesBottomBackRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackRightCorner ) .def("nodesBottomRightBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightBackCorner ) .def("nodesRightBackBottomCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackBottomCorner ) .def("nodesRightBottomBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomBackCorner ) .def("nodesBackLeftTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftTopCorner ) .def("nodesTopBackLeftCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackLeftCorner ) .def("nodesTopLeftBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftBackCorner ) .def("nodesLeftBackTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackTopCorner ) .def("nodesLeftTopBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopBackCorner ) .def("nodesBackRightTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightTopCorner ) .def("nodesTopBackRightCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackRightCorner ) .def("nodesTopRightBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightBackCorner ) .def("nodesRightBackTopCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackTopCorner ) .def("nodesRightTopBackCorner" , &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopBackCorner ) // periodicity .def("nodesPeriodic" , &GooseFEM::Mesh::Hex8::FineLayer::nodesPeriodic ) .def("nodesOrigin" , &GooseFEM::Mesh::Hex8::FineLayer::nodesOrigin ) .def("dofs" , &GooseFEM::Mesh::Hex8::FineLayer::dofs ) .def("dofsPeriodic" , &GooseFEM::Mesh::Hex8::FineLayer::dofsPeriodic ) // print to screen .def("__repr__", [](const GooseFEM::Mesh::Hex8::FineLayer &a){ return ""; } ); // ===================================== GooseFEM/MeshQuad4.h ===================================== py::class_(mMeshQuad4, "Regular") .def( py::init(), "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge size 'h'", py::arg("nx"), py::arg("ny"), py::arg("h")=1. ) .def("coor" , &GooseFEM::Mesh::Quad4::Regular::coor ) .def("conn" , &GooseFEM::Mesh::Quad4::Regular::conn ) .def("nelem" , &GooseFEM::Mesh::Quad4::Regular::nelem ) .def("nnode" , &GooseFEM::Mesh::Quad4::Regular::nnode ) .def("nne" , &GooseFEM::Mesh::Quad4::Regular::nne ) .def("ndim" , &GooseFEM::Mesh::Quad4::Regular::ndim ) .def("nodesBottomEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesBottomEdge ) .def("nodesTopEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesTopEdge ) .def("nodesLeftEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesLeftEdge ) .def("nodesRightEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesRightEdge ) + .def("nodesBottomOpenEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesBottomOpenEdge ) + .def("nodesTopOpenEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesTopOpenEdge ) + .def("nodesLeftOpenEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesLeftOpenEdge ) + .def("nodesRightOpenEdge" , &GooseFEM::Mesh::Quad4::Regular::nodesRightOpenEdge ) .def("nodesBottomLeftCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesBottomLeftCorner ) .def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::Regular::nodesBottomRightCorner) .def("nodesTopLeftCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesTopLeftCorner ) .def("nodesTopRightCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesTopRightCorner ) .def("nodesLeftBottomCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesLeftBottomCorner ) .def("nodesLeftTopCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesLeftTopCorner ) .def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::Regular::nodesRightBottomCorner) .def("nodesRightTopCorner" , &GooseFEM::Mesh::Quad4::Regular::nodesRightTopCorner ) .def("nodesPeriodic" , &GooseFEM::Mesh::Quad4::Regular::nodesPeriodic ) .def("nodesOrigin" , &GooseFEM::Mesh::Quad4::Regular::nodesOrigin ) .def("dofs" , &GooseFEM::Mesh::Quad4::Regular::dofs ) .def("dofsPeriodic" , &GooseFEM::Mesh::Quad4::Regular::dofsPeriodic ) .def("__repr__", [](const GooseFEM::Mesh::Quad4::Regular &a){ return ""; } ); // ------------------------------------------------------------------------------------------------- py::class_(mMeshQuad4, "FineLayer") .def( - py::init(), + py::init(), "FineLayer mesh: 'nx' pixels in horizontal direction (length 'Lx'), idem in vertical direction", py::arg("nx"), py::arg("ny"), py::arg("h")=1., - py::arg("nfine")=1, - py::arg("nskip")=0 + py::arg("nfine")=1 ) .def("shape" , &GooseFEM::Mesh::Quad4::FineLayer::shape ) .def("coor" , &GooseFEM::Mesh::Quad4::FineLayer::coor ) .def("conn" , &GooseFEM::Mesh::Quad4::FineLayer::conn ) .def("nelem" , &GooseFEM::Mesh::Quad4::FineLayer::nelem ) .def("nnode" , &GooseFEM::Mesh::Quad4::FineLayer::nnode ) .def("nne" , &GooseFEM::Mesh::Quad4::FineLayer::nne ) .def("ndim" , &GooseFEM::Mesh::Quad4::FineLayer::ndim ) .def("elementsMiddleLayer" , &GooseFEM::Mesh::Quad4::FineLayer::elementsMiddleLayer ) .def("nodesBottomEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomEdge ) .def("nodesTopEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesTopEdge ) .def("nodesLeftEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftEdge ) .def("nodesRightEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesRightEdge ) + .def("nodesBottomOpenEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomOpenEdge ) + .def("nodesTopOpenEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesTopOpenEdge ) + .def("nodesLeftOpenEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftOpenEdge ) + .def("nodesRightOpenEdge" , &GooseFEM::Mesh::Quad4::FineLayer::nodesRightOpenEdge ) .def("nodesBottomLeftCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomLeftCorner ) .def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomRightCorner) .def("nodesTopLeftCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesTopLeftCorner ) .def("nodesTopRightCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesTopRightCorner ) .def("nodesLeftBottomCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftBottomCorner ) .def("nodesLeftTopCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftTopCorner ) .def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightBottomCorner) .def("nodesRightTopCorner" , &GooseFEM::Mesh::Quad4::FineLayer::nodesRightTopCorner ) .def("nodesPeriodic" , &GooseFEM::Mesh::Quad4::FineLayer::nodesPeriodic ) .def("nodesOrigin" , &GooseFEM::Mesh::Quad4::FineLayer::nodesOrigin ) .def("dofs" , &GooseFEM::Mesh::Quad4::FineLayer::dofs ) .def("dofsPeriodic" , &GooseFEM::Mesh::Quad4::FineLayer::dofsPeriodic ) .def("__repr__", [](const GooseFEM::Mesh::Quad4::FineLayer &a){ return ""; } ); // ===================================== GooseFEM/MeshTri3.h ====================================== py::class_(mMeshTri3, "Regular") .def( py::init(), "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge size 'h'", py::arg("nx"), py::arg("ny"), py::arg("h")=1. ) .def("coor" , &GooseFEM::Mesh::Tri3::Regular::coor ) .def("conn" , &GooseFEM::Mesh::Tri3::Regular::conn ) .def("nelem" , &GooseFEM::Mesh::Tri3::Regular::nelem ) .def("nnode" , &GooseFEM::Mesh::Tri3::Regular::nnode ) .def("nne" , &GooseFEM::Mesh::Tri3::Regular::nne ) .def("ndim" , &GooseFEM::Mesh::Tri3::Regular::ndim ) .def("nodesBottomEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesBottomEdge ) .def("nodesTopEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesTopEdge ) .def("nodesLeftEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftEdge ) .def("nodesRightEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesRightEdge ) + .def("nodesBottomOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesBottomOpenEdge ) + .def("nodesTopOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesTopOpenEdge ) + .def("nodesLeftOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftOpenEdge ) + .def("nodesRightOpenEdge" , &GooseFEM::Mesh::Tri3::Regular::nodesRightOpenEdge ) .def("nodesBottomLeftCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesBottomLeftCorner ) .def("nodesBottomRightCorner", &GooseFEM::Mesh::Tri3::Regular::nodesBottomRightCorner) .def("nodesTopLeftCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesTopLeftCorner ) .def("nodesTopRightCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesTopRightCorner ) .def("nodesLeftBottomCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftBottomCorner ) .def("nodesLeftTopCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesLeftTopCorner ) .def("nodesRightBottomCorner", &GooseFEM::Mesh::Tri3::Regular::nodesRightBottomCorner) .def("nodesRightTopCorner" , &GooseFEM::Mesh::Tri3::Regular::nodesRightTopCorner ) .def("nodesPeriodic" , &GooseFEM::Mesh::Tri3::Regular::nodesPeriodic ) .def("nodesOrigin" , &GooseFEM::Mesh::Tri3::Regular::nodesOrigin ) .def("dofs" , &GooseFEM::Mesh::Tri3::Regular::dofs ) .def("dofsPeriodic" , &GooseFEM::Mesh::Tri3::Regular::dofsPeriodic ) .def("__repr__", [](const GooseFEM::Mesh::Tri3::Regular &a){ return ""; } ); // ------------------------------------------------------------------------------------------------- mMeshTri3.def("getOrientation",&GooseFEM::Mesh::Tri3::getOrientation, "Get the orientation of each element", py::arg("coor"), py::arg("conn") ); // ------------------------------------------------------------------------------------------------- mMeshTri3.def("retriangulate",&GooseFEM::Mesh::Tri3::retriangulate, "Re-triangulate existing mesh", py::arg("coor"), py::arg("conn"), py::arg("orientation")=-1 ); // ================================================================================================= } diff --git a/src/GooseFEM/MeshHex8.cpp b/src/GooseFEM/MeshHex8.cpp index f27a375..9d181f1 100644 --- a/src/GooseFEM/MeshHex8.cpp +++ b/src/GooseFEM/MeshHex8.cpp @@ -1,2608 +1,2625 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESHHEX8_CPP #define GOOSEFEM_MESHHEX8_CPP // ------------------------------------------------------------------------------------------------- #include "MeshHex8.h" // ===================================== GooseFEM::Mesh::Hex8 ====================================== namespace GooseFEM { namespace Mesh { namespace Hex8 { // ===================================== CLASS - REGULAR MESH ====================================== // ------------------------------------------ constructor ------------------------------------------ inline Regular::Regular(size_t nelx, size_t nely, size_t nelz, double h): m_h(h), m_nelx(nelx), m_nely(nely), m_nelz(nelz) { assert( m_nelx >= 1 ); assert( m_nely >= 1 ); assert( m_nelz >= 1 ); m_nnode = (m_nelx+1) * (m_nely+1) * (m_nelz+1); m_nelem = m_nelx * m_nely * m_nelz ; } // -------------------------------------- number of elements --------------------------------------- inline size_t Regular::nelem() { return m_nelem; } // ---------------------------------------- number of nodes ---------------------------------------- inline size_t Regular::nnode() { return m_nnode; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t Regular::nne() { return m_nne; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t Regular::ndim() { return m_ndim; } // --------------------------------- coordinates (nodal positions) --------------------------------- inline MatD Regular::coor() { MatD out(m_nnode,m_ndim); ColD x = ColD::LinSpaced(m_nelx+1, 0.0, m_h*static_cast(m_nelx)); ColD y = ColD::LinSpaced(m_nely+1, 0.0, m_h*static_cast(m_nely)); ColD z = ColD::LinSpaced(m_nelz+1, 0.0, m_h*static_cast(m_nelz)); size_t inode = 0; for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) { for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) { for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) { out(inode,0) = x(ix); out(inode,1) = y(iy); out(inode,2) = z(iz); ++inode; } } } return out; } // ---------------------------- connectivity (node-numbers per element) ---------------------------- inline MatS Regular::conn() { MatS out(m_nelem,m_nne); size_t ielem = 0; for ( size_t iz = 0 ; iz < m_nelz ; ++iz ) { for ( size_t iy = 0 ; iy < m_nely ; ++iy ) { for ( size_t ix = 0 ; ix < m_nelx ; ++ix ) { - out(ielem,0) = (iz+0)*(m_nely+1)*(m_nelx+1) + (iy+0)*(m_nelx+1) + (ix+0); - out(ielem,1) = (iz+0)*(m_nely+1)*(m_nelx+1) + (iy+0)*(m_nelx+1) + (ix+1); - out(ielem,3) = (iz+0)*(m_nely+1)*(m_nelx+1) + (iy+1)*(m_nelx+1) + (ix+0); - out(ielem,2) = (iz+0)*(m_nely+1)*(m_nelx+1) + (iy+1)*(m_nelx+1) + (ix+1); - out(ielem,4) = (iz+1)*(m_nely+1)*(m_nelx+1) + (iy+0)*(m_nelx+1) + (ix+0); - out(ielem,5) = (iz+1)*(m_nely+1)*(m_nelx+1) + (iy+0)*(m_nelx+1) + (ix+1); - out(ielem,7) = (iz+1)*(m_nely+1)*(m_nelx+1) + (iy+1)*(m_nelx+1) + (ix+0); - out(ielem,6) = (iz+1)*(m_nely+1)*(m_nelx+1) + (iy+1)*(m_nelx+1) + (ix+1); + out(ielem,0) = (iy )*(m_nelx+1) + (ix ) + (iz )*(m_nely+1)*(m_nelx+1); + out(ielem,1) = (iy )*(m_nelx+1) + (ix+1) + (iz )*(m_nely+1)*(m_nelx+1); + out(ielem,3) = (iy+1)*(m_nelx+1) + (ix ) + (iz )*(m_nely+1)*(m_nelx+1); + out(ielem,2) = (iy+1)*(m_nelx+1) + (ix+1) + (iz )*(m_nely+1)*(m_nelx+1); + out(ielem,4) = (iy )*(m_nelx+1) + (ix ) + (iz+1)*(m_nely+1)*(m_nelx+1); + out(ielem,5) = (iy )*(m_nelx+1) + (ix+1) + (iz+1)*(m_nely+1)*(m_nelx+1); + out(ielem,7) = (iy+1)*(m_nelx+1) + (ix ) + (iz+1)*(m_nely+1)*(m_nelx+1); + out(ielem,6) = (iy+1)*(m_nelx+1) + (ix+1) + (iz+1)*(m_nely+1)*(m_nelx+1); ++ielem; } } } return out; } // ------------------------------ node-numbers along the front plane ------------------------------- inline ColS Regular::nodesFront() { ColS out((m_nelx+1)*(m_nely+1)); for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(iy*(m_nelx+1)+ix) = iy*(m_nelx+1) + ix; return out; } // ------------------------------- node-numbers along the back plane ------------------------------- inline ColS Regular::nodesBack() { ColS out((m_nelx+1)*(m_nely+1)); for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(iy*(m_nelx+1)+ix) = iy*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1); return out; } // ------------------------------- node-numbers along the left plane ------------------------------- inline ColS Regular::nodesLeft() { ColS out((m_nely+1)*(m_nelz+1)); for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) out(iz*(m_nely+1)+iy) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1); return out; } // ------------------------------ node-numbers along the right plane ------------------------------- inline ColS Regular::nodesRight() { ColS out((m_nely+1)*(m_nelz+1)); for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) out(iz*(m_nely+1)+iy) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx; return out; } // ------------------------------ node-numbers along the bottom plane ------------------------------ inline ColS Regular::nodesBottom() { ColS out((m_nelx+1)*(m_nelz+1)); for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(iz*(m_nelx+1)+ix) = ix + iz*(m_nelx+1)*(m_nely+1); return out; } // ------------------------------- node-numbers along the top plane -------------------------------- inline ColS Regular::nodesTop() { ColS out((m_nelx+1)*(m_nelz+1)); for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(iz*(m_nelx+1)+ix) = ix + m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1); return out; } // ------------------------------- node-numbers along the front face ------------------------------- inline ColS Regular::nodesFrontFace() { ColS out((m_nelx-1)*(m_nely-1)); for ( size_t iy = 1 ; iy < m_nely ; ++iy ) for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) out((iy-1)*(m_nelx-1)+(ix-1)) = iy*(m_nelx+1) + ix; return out; } // ------------------------------- node-numbers along the back face -------------------------------- inline ColS Regular::nodesBackFace() { ColS out((m_nelx-1)*(m_nely-1)); for ( size_t iy = 1 ; iy < m_nely ; ++iy ) { for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) { out((iy-1)*(m_nelx-1)+(ix-1)) = iy*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1); } } return out; } // ------------------------------- node-numbers along the left face -------------------------------- inline ColS Regular::nodesLeftFace() { ColS out((m_nely-1)*(m_nelz-1)); for ( size_t iz = 1 ; iz < m_nelz ; ++iz ) for ( size_t iy = 1 ; iy < m_nely ; ++iy ) out((iz-1)*(m_nely-1)+(iy-1)) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1); return out; } // ------------------------------- node-numbers along the right face ------------------------------- inline ColS Regular::nodesRightFace() { ColS out((m_nely-1)*(m_nelz-1)); for ( size_t iz = 1 ; iz < m_nelz ; ++iz ) for ( size_t iy = 1 ; iy < m_nely ; ++iy ) out((iz-1)*(m_nely-1)+(iy-1)) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx; return out; } // ------------------------------ node-numbers along the bottom face ------------------------------- inline ColS Regular::nodesBottomFace() { ColS out((m_nelx-1)*(m_nelz-1)); for ( size_t iz = 1 ; iz < m_nelz ; ++iz ) for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) out((iz-1)*(m_nelx-1)+(ix-1)) = ix + iz*(m_nelx+1)*(m_nely+1); return out; } // -------------------------------- node-numbers along the top face -------------------------------- inline ColS Regular::nodesTopFace() { ColS out((m_nelx-1)*(m_nelz-1)); for ( size_t iz = 1 ; iz < m_nelz ; ++iz ) for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) out((iz-1)*(m_nelx-1)+(ix-1)) = ix + m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1); return out; } // --------------------------- node-numbers along the front-bottom edge ---------------------------- inline ColS Regular::nodesFrontBottomEdge() { ColS out(m_nelx+1); for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(ix) = ix; return out; } // ----------------------------- node-numbers along the front-top edge ----------------------------- inline ColS Regular::nodesFrontTopEdge() { ColS out(m_nelx+1); for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) - out(ix) = m_nely*(m_nelx+1) + ix; + out(ix) = ix + m_nely*(m_nelx+1); return out; } // ---------------------------- node-numbers along the front-left edge ----------------------------- inline ColS Regular::nodesFrontLeftEdge() { ColS out(m_nely+1); for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) out(iy) = iy*(m_nelx+1); return out; } // ---------------------------- node-numbers along the front-right edge ---------------------------- inline ColS Regular::nodesFrontRightEdge() { ColS out(m_nely+1); for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) out(iy) = iy*(m_nelx+1) + m_nelx; return out; } // ---------------------------- node-numbers along the back-bottom edge ---------------------------- inline ColS Regular::nodesBackBottomEdge() { ColS out(m_nelx+1); for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(ix) = ix + m_nelz*(m_nely+1)*(m_nelx+1); return out; } // ----------------------------- node-numbers along the back-top edge ------------------------------ inline ColS Regular::nodesBackTopEdge() { ColS out(m_nelx+1); for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(ix) = m_nely*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1); return out; } // ----------------------------- node-numbers along the back-left edge ----------------------------- inline ColS Regular::nodesBackLeftEdge() { ColS out(m_nely+1); for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) out(iy) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1); return out; } // ---------------------------- node-numbers along the back-right edge ----------------------------- inline ColS Regular::nodesBackRightEdge() { ColS out(m_nely+1); for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) out(iy) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1) + m_nelx; return out; } // ---------------------------- node-numbers along the bottom-left edge ---------------------------- inline ColS Regular::nodesBottomLeftEdge() { ColS out(m_nelz+1); for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) out(iz) = iz*(m_nelx+1)*(m_nely+1); return out; } // --------------------------- node-numbers along the bottom-right edge ---------------------------- inline ColS Regular::nodesBottomRightEdge() { ColS out(m_nelz+1); for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) out(iz) = iz*(m_nelx+1)*(m_nely+1) + m_nelx; return out; } // ----------------------------- node-numbers along the top-left edge ------------------------------ inline ColS Regular::nodesTopLeftEdge() { ColS out(m_nelz+1); for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) out(iz) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1); return out; } // ----------------------------- node-numbers along the top-right edge ----------------------------- inline ColS Regular::nodesTopRightEdge() { ColS out(m_nelz+1); for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) out(iz) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx; return out; } // -------------------------------------------- aliases -------------------------------------------- inline ColS Regular::nodesBottomFrontEdge() { return nodesFrontBottomEdge(); } inline ColS Regular::nodesBottomBackEdge() { return nodesBackBottomEdge(); } inline ColS Regular::nodesTopFrontEdge() { return nodesFrontTopEdge(); } inline ColS Regular::nodesTopBackEdge() { return nodesBackTopEdge(); } inline ColS Regular::nodesLeftBottomEdge() { return nodesBottomLeftEdge(); } inline ColS Regular::nodesLeftFrontEdge() { return nodesFrontLeftEdge(); } inline ColS Regular::nodesLeftBackEdge() { return nodesBackLeftEdge(); } inline ColS Regular::nodesLeftTopEdge() { return nodesTopLeftEdge(); } inline ColS Regular::nodesRightBottomEdge() { return nodesBottomRightEdge(); } inline ColS Regular::nodesRightTopEdge() { return nodesTopRightEdge(); } inline ColS Regular::nodesRightFrontEdge() { return nodesFrontRightEdge(); } inline ColS Regular::nodesRightBackEdge() { return nodesBackRightEdge(); } // ------------------- node-numbers along the front-bottom edge, without corners ------------------- inline ColS Regular::nodesFrontBottomOpenEdge() { ColS out(m_nelx-1); for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) out(ix-1) = ix; return out; } // -------------------- node-numbers along the front-top edge, without corners --------------------- inline ColS Regular::nodesFrontTopOpenEdge() { ColS out(m_nelx-1); for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) - out(ix-1) = m_nely*(m_nelx+1) + ix; + out(ix-1) = ix + m_nely*(m_nelx+1); return out; } // -------------------- node-numbers along the front-left edge, without corners -------------------- inline ColS Regular::nodesFrontLeftOpenEdge() { ColS out(m_nely-1); for ( size_t iy = 1 ; iy < m_nely ; ++iy ) out(iy-1) = iy*(m_nelx+1); return out; } // ------------------- node-numbers along the front-right edge, without corners -------------------- inline ColS Regular::nodesFrontRightOpenEdge() { ColS out(m_nely-1); for ( size_t iy = 1 ; iy < m_nely ; ++iy ) out(iy-1) = iy*(m_nelx+1) + m_nelx; return out; } // ------------------- node-numbers along the back-bottom edge, without corners -------------------- inline ColS Regular::nodesBackBottomOpenEdge() { ColS out(m_nelx-1); for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) out(ix-1) = ix + m_nelz*(m_nely+1)*(m_nelx+1); return out; } // --------------------- node-numbers along the back-top edge, without corners --------------------- inline ColS Regular::nodesBackTopOpenEdge() { ColS out(m_nelx-1); for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) out(ix-1) = m_nely*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1); return out; } // -------------------- node-numbers along the back-left edge, without corners --------------------- inline ColS Regular::nodesBackLeftOpenEdge() { ColS out(m_nely-1); for ( size_t iy = 1 ; iy < m_nely ; ++iy ) out(iy-1) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1); return out; } // -------------------- node-numbers along the back-right edge, without corners -------------------- inline ColS Regular::nodesBackRightOpenEdge() { ColS out(m_nely-1); for ( size_t iy = 1 ; iy < m_nely ; ++iy ) out(iy-1) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1) + m_nelx; return out; } // ------------------- node-numbers along the bottom-left edge, without corners -------------------- inline ColS Regular::nodesBottomLeftOpenEdge() { ColS out(m_nelz-1); for ( size_t iz = 1 ; iz < m_nelz ; ++iz ) out(iz-1) = iz*(m_nelx+1)*(m_nely+1); return out; } // ------------------- node-numbers along the bottom-right edge, without corners ------------------- inline ColS Regular::nodesBottomRightOpenEdge() { ColS out(m_nelz-1); for ( size_t iz = 1 ; iz < m_nelz ; ++iz ) out(iz-1) = iz*(m_nelx+1)*(m_nely+1) + m_nelx; return out; } // --------------------- node-numbers along the top-left edge, without corners --------------------- inline ColS Regular::nodesTopLeftOpenEdge() { ColS out(m_nelz-1); for ( size_t iz = 1 ; iz < m_nelz ; ++iz ) out(iz-1) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1); return out; } // -------------------- node-numbers along the top-right edge, without corners --------------------- inline ColS Regular::nodesTopRightOpenEdge() { ColS out(m_nelz-1); for ( size_t iz = 1 ; iz < m_nelz ; ++iz ) out(iz-1) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx; return out; } // -------------------------------------------- aliases -------------------------------------------- inline ColS Regular::nodesBottomFrontOpenEdge() { return nodesFrontBottomOpenEdge(); } inline ColS Regular::nodesBottomBackOpenEdge() { return nodesBackBottomOpenEdge(); } inline ColS Regular::nodesTopFrontOpenEdge() { return nodesFrontTopOpenEdge(); } inline ColS Regular::nodesTopBackOpenEdge() { return nodesBackTopOpenEdge(); } inline ColS Regular::nodesLeftBottomOpenEdge() { return nodesBottomLeftOpenEdge(); } inline ColS Regular::nodesLeftFrontOpenEdge() { return nodesFrontLeftOpenEdge(); } inline ColS Regular::nodesLeftBackOpenEdge() { return nodesBackLeftOpenEdge(); } inline ColS Regular::nodesLeftTopOpenEdge() { return nodesTopLeftOpenEdge(); } inline ColS Regular::nodesRightBottomOpenEdge() { return nodesBottomRightOpenEdge(); } inline ColS Regular::nodesRightTopOpenEdge() { return nodesTopRightOpenEdge(); } inline ColS Regular::nodesRightFrontOpenEdge() { return nodesFrontRightOpenEdge(); } inline ColS Regular::nodesRightBackOpenEdge() { return nodesBackRightOpenEdge(); } // -------------------------- node-number of the front-bottom-left corner -------------------------- inline size_t Regular::nodesFrontBottomLeftCorner() { return 0; } // ------------------------- node-number of the front-bottom-right corner -------------------------- inline size_t Regular::nodesFrontBottomRightCorner() { return m_nelx; } // --------------------------- node-number of the front-top-left corner ---------------------------- inline size_t Regular::nodesFrontTopLeftCorner() { return m_nely*(m_nelx+1); } // --------------------------- node-number of the front-top-right corner --------------------------- inline size_t Regular::nodesFrontTopRightCorner() { return m_nely*(m_nelx+1) + m_nelx; } // -------------------------- node-number of the back-bottom-left corner --------------------------- inline size_t Regular::nodesBackBottomLeftCorner() { return m_nelz*(m_nely+1)*(m_nelx+1); } // -------------------------- node-number of the back-bottom-right corner -------------------------- inline size_t Regular::nodesBackBottomRightCorner() { return m_nelx + m_nelz*(m_nely+1)*(m_nelx+1); } // ---------------------------- node-number of the back-top-left corner ---------------------------- inline size_t Regular::nodesBackTopLeftCorner() { return m_nely*(m_nelx+1) + m_nelz*(m_nely+1)*(m_nelx+1); } // --------------------------- node-number of the back-top-right corner ---------------------------- inline size_t Regular::nodesBackTopRightCorner() { return m_nely*(m_nelx+1) + m_nelx + m_nelz*(m_nely+1)*(m_nelx+1); } // -------------------------------------------- aliases -------------------------------------------- inline size_t Regular::nodesFrontLeftBottomCorner() { return nodesFrontBottomLeftCorner(); } inline size_t Regular::nodesBottomFrontLeftCorner() { return nodesFrontBottomLeftCorner(); } inline size_t Regular::nodesBottomLeftFrontCorner() { return nodesFrontBottomLeftCorner(); } inline size_t Regular::nodesLeftFrontBottomCorner() { return nodesFrontBottomLeftCorner(); } inline size_t Regular::nodesLeftBottomFrontCorner() { return nodesFrontBottomLeftCorner(); } inline size_t Regular::nodesFrontRightBottomCorner() { return nodesFrontBottomRightCorner(); } inline size_t Regular::nodesBottomFrontRightCorner() { return nodesFrontBottomRightCorner(); } inline size_t Regular::nodesBottomRightFrontCorner() { return nodesFrontBottomRightCorner(); } inline size_t Regular::nodesRightFrontBottomCorner() { return nodesFrontBottomRightCorner(); } inline size_t Regular::nodesRightBottomFrontCorner() { return nodesFrontBottomRightCorner(); } inline size_t Regular::nodesFrontLeftTopCorner() { return nodesFrontTopLeftCorner(); } inline size_t Regular::nodesTopFrontLeftCorner() { return nodesFrontTopLeftCorner(); } inline size_t Regular::nodesTopLeftFrontCorner() { return nodesFrontTopLeftCorner(); } inline size_t Regular::nodesLeftFrontTopCorner() { return nodesFrontTopLeftCorner(); } inline size_t Regular::nodesLeftTopFrontCorner() { return nodesFrontTopLeftCorner(); } inline size_t Regular::nodesFrontRightTopCorner() { return nodesFrontTopRightCorner(); } inline size_t Regular::nodesTopFrontRightCorner() { return nodesFrontTopRightCorner(); } inline size_t Regular::nodesTopRightFrontCorner() { return nodesFrontTopRightCorner(); } inline size_t Regular::nodesRightFrontTopCorner() { return nodesFrontTopRightCorner(); } inline size_t Regular::nodesRightTopFrontCorner() { return nodesFrontTopRightCorner(); } inline size_t Regular::nodesBackLeftBottomCorner() { return nodesBackBottomLeftCorner(); } inline size_t Regular::nodesBottomBackLeftCorner() { return nodesBackBottomLeftCorner(); } inline size_t Regular::nodesBottomLeftBackCorner() { return nodesBackBottomLeftCorner(); } inline size_t Regular::nodesLeftBackBottomCorner() { return nodesBackBottomLeftCorner(); } inline size_t Regular::nodesLeftBottomBackCorner() { return nodesBackBottomLeftCorner(); } inline size_t Regular::nodesBackRightBottomCorner() { return nodesBackBottomRightCorner(); } inline size_t Regular::nodesBottomBackRightCorner() { return nodesBackBottomRightCorner(); } inline size_t Regular::nodesBottomRightBackCorner() { return nodesBackBottomRightCorner(); } inline size_t Regular::nodesRightBackBottomCorner() { return nodesBackBottomRightCorner(); } inline size_t Regular::nodesRightBottomBackCorner() { return nodesBackBottomRightCorner(); } inline size_t Regular::nodesBackLeftTopCorner() { return nodesBackTopLeftCorner(); } inline size_t Regular::nodesTopBackLeftCorner() { return nodesBackTopLeftCorner(); } inline size_t Regular::nodesTopLeftBackCorner() { return nodesBackTopLeftCorner(); } inline size_t Regular::nodesLeftBackTopCorner() { return nodesBackTopLeftCorner(); } inline size_t Regular::nodesLeftTopBackCorner() { return nodesBackTopLeftCorner(); } inline size_t Regular::nodesBackRightTopCorner() { return nodesBackTopRightCorner(); } inline size_t Regular::nodesTopBackRightCorner() { return nodesBackTopRightCorner(); } inline size_t Regular::nodesTopRightBackCorner() { return nodesBackTopRightCorner(); } inline size_t Regular::nodesRightBackTopCorner() { return nodesBackTopRightCorner(); } inline size_t Regular::nodesRightTopBackCorner() { return nodesBackTopRightCorner(); } // ------------------------------ node-numbers of periodic node-pairs ------------------------------ inline MatS Regular::nodesPeriodic() { // faces ColS fro = nodesFrontFace(); ColS bck = nodesBackFace(); ColS lft = nodesLeftFace(); ColS rgt = nodesRightFace(); ColS bot = nodesBottomFace(); ColS top = nodesTopFace(); // edges ColS froBot = nodesFrontBottomOpenEdge(); ColS froTop = nodesFrontTopOpenEdge(); ColS froLft = nodesFrontLeftOpenEdge(); ColS froRgt = nodesFrontRightOpenEdge(); ColS bckBot = nodesBackBottomOpenEdge(); ColS bckTop = nodesBackTopOpenEdge(); ColS bckLft = nodesBackLeftOpenEdge(); ColS bckRgt = nodesBackRightOpenEdge(); ColS botLft = nodesBottomLeftOpenEdge(); ColS botRgt = nodesBottomRightOpenEdge(); ColS topLft = nodesTopLeftOpenEdge(); ColS topRgt = nodesTopRightOpenEdge(); // allocate nodal ties // - number of tying per category size_t tface = fro.size() + lft.size() + bot.size(); size_t tedge = 3*froBot.size() + 3*froLft.size() + 3*botLft.size(); size_t tnode = 7; // - allocate MatS out(tface+tedge+tnode, 2); // counter size_t i = 0; // tie all corners to one corner out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontBottomRightCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomRightCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomLeftCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopLeftCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopRightCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopRightCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopLeftCorner(); ++i; // tie all corresponding edges to each other (exclude corners) for ( auto j = 0 ; j(nodePer.rows()); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs for ( size_t i = 0 ; i < nper ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) out(nodePer(i,1),j) = out(nodePer(i,0),j); // renumber "out" to be sequential return GooseFEM::Mesh::renumber(out); } // ==================================== CLASS - FINELAYER MESH ===================================== // ------------------------------------------ constructor ------------------------------------------ inline FineLayer::FineLayer(size_t nelx, size_t nely, size_t nelz, double h, size_t nfine): m_h(h), m_nelx(nelx), m_nelz(nelz) { // basic assumptions assert( nelx >= 1 ); assert( nely >= 1 ); assert( nelz >= 1 ); // store basic info m_Lx = m_h * static_cast(nelx); m_Lz = m_h * static_cast(nelz); // compute element size in y-direction (use symmetry, compute upper half) // ---------------------------------------------------------------------- // temporary variables size_t nmin, ntot; ColS nhx(nely), nhy(nely), nhz(nely); ColI refine(nely); // minimum height in y-direction (half of the height because of symmetry) if ( nely % 2 == 0 ) nmin = nely /2; else nmin = (nely +1)/2; // minimum number of fine layers in y-direction (minimum 1, middle layer part of this half) if ( nfine % 2 == 0 ) nfine = nfine /2 + 1; else nfine = (nfine+1)/2; if ( nfine < 1 ) nfine = 1; if ( nfine > nmin ) nfine = nmin; // initialize to state with only fine elements nhx .setOnes(); nhy .setOnes(); nhz .setOnes(); refine.setConstant(-1); // loop over element layers in y-direction, try to coarsen using these rules: // (1) element size in y-direction <= distance to origin in y-direction // (2) element size in x-(z-)direction should fit the total number of elements in x-(z-)direction // (3) a certain number of layers have the minimum size "1" (are fine) for ( size_t iy = nfine ; ; ) { // initialize current size in y-direction if ( iy == nfine ) ntot = nfine; // check to stop if ( iy >= nely or ntot >= nmin ) { nely = iy; break; } // rules (1,2) satisfied: coarsen in x-direction (and z-direction) if ( 3*nhy(iy) <= ntot and nelx%(3*nhx(iy)) == 0 ) { // - process refinement in x-direction refine (iy ) = 0; nhy (iy ) *= 2; nhy.segment(iy+1,nely-iy-1) *= 3; nhx.segment(iy ,nely-iy ) *= 3; // - rule (2) satisfied: coarsen next element layer in z-direction if ( iy+1 < nely and ntot+nhy(iy) < nmin ) { if ( nelz%(3*nhz(iy+1)) == 0 ) { // - update the number of elements in y-direction ntot += nhy(iy); // - proceed to next element layer in y-direction ++iy; // - process refinement in z-direction refine (iy ) = 2; nhy (iy ) = nhy(iy-1); nhz.segment(iy,nely-iy) *= 3; } } } // rules (1,2) satisfied: coarse in z-direction else if ( 3*nhy(iy) <= ntot and nelz%(3*nhz(iy)) == 0 ) { // - process refinement in z-direction refine (iy ) = 2; nhy (iy ) *= 2; nhy.segment(iy+1,nely-iy-1) *= 3; nhz.segment(iy ,nely-iy ) *= 3; } // update the number of elements in y-direction ntot += nhy(iy); // proceed to next element layer in y-direction ++iy; // check to stop if ( iy >= nely or ntot >= nmin ) { nely = iy; break; } } // symmetrize, compute full information // ------------------------------------ // allocate mesh constructor parameters m_nhx .conservativeResize(nely*2-1); m_nhy .conservativeResize(nely*2-1); m_nhz .conservativeResize(nely*2-1); m_refine .conservativeResize(nely*2-1); m_nelx .conservativeResize(nely*2-1); m_nelz .conservativeResize(nely*2-1); m_nnd .conservativeResize(nely*2 ); m_startElem.conservativeResize(nely*2-1); m_startNode.conservativeResize(nely*2 ); // fill // - lower half for ( size_t iy = 0 ; iy < nely ; ++iy ) { m_nhx (iy) = nhx (nely-iy-1); m_nhy (iy) = nhy (nely-iy-1); m_nhz (iy) = nhz (nely-iy-1); m_refine(iy) = refine(nely-iy-1); } // - upper half for ( size_t iy = 0 ; iy < nely-1 ; ++iy ) { m_nhx (iy+nely) = nhx (iy+1); m_nhy (iy+nely) = nhy (iy+1); m_nhz (iy+nely) = nhz (iy+1); m_refine(iy+nely) = refine(iy+1); } // update size nely = m_nhx.size(); // compute the number of elements per element layer in y-direction for ( size_t iy = 0 ; iy < nely ; ++iy ) { m_nelx(iy) = nelx / m_nhx(iy); m_nelz(iy) = nelz / m_nhz(iy); } // compute the number of nodes per node layer in y-direction // - bottom half for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) m_nnd(iy) = (m_nelx(iy)+1) * (m_nelz(iy)+1); // - top half for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) m_nnd(iy+1) = (m_nelx(iy)+1) * (m_nelz(iy)+1); // compute mesh dimensions // ----------------------- // initialize m_nnode = 0; m_nelem = 0; m_startNode(0) = 0; // loop over element layers (bottom -> middle, elements become finer) for ( size_t i = 0 ; i < (nely-1)/2 ; ++i ) { // - store the first element of the layer m_startElem(i) = m_nelem; // - add the nodes of this layer if ( m_refine(i) == 0 ) { m_nnode += (3*m_nelx(i)+1) * ( m_nelz(i)+1); } else if ( m_refine(i) == 2 ) { m_nnode += ( m_nelx(i)+1) * (3*m_nelz(i)+1); } else { m_nnode += ( m_nelx(i)+1) * ( m_nelz(i)+1); } // - add the elements of this layer if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ) * ( m_nelz(i) ); } else if ( m_refine(i) == 2 ) { m_nelem += ( m_nelx(i) ) * (4*m_nelz(i) ); } else { m_nelem += ( m_nelx(i) ) * ( m_nelz(i) ); } // - store the starting node of the next layer m_startNode(i+1) = m_nnode; } // loop over element layers (middle -> top, elements become coarser) for ( size_t i = (nely-1)/2 ; i < nely ; ++i ) { // - store the first element of the layer m_startElem(i) = m_nelem; // - add the nodes of this layer if ( m_refine(i) == 0 ) { m_nnode += (5*m_nelx(i)+1) * ( m_nelz(i)+1); } else if ( m_refine(i) == 2 ) { m_nnode += ( m_nelx(i)+1) * (5*m_nelz(i)+1); } else { m_nnode += ( m_nelx(i)+1) * ( m_nelz(i)+1); } // - add the elements of this layer if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ) * ( m_nelz(i) ); } else if ( m_refine(i) == 2 ) { m_nelem += ( m_nelx(i) ) * (4*m_nelz(i) ); } else { m_nelem += ( m_nelx(i) ) * ( m_nelz(i) ); } // - store the starting node of the next layer m_startNode(i+1) = m_nnode; } // - add the top row of nodes m_nnode += (m_nelx(nely-1)+1) * (m_nelz(nely-1)+1); } // -------------------------------------- number of elements --------------------------------------- inline size_t FineLayer::nelem() { return m_nelem; } // ---------------------------------------- number of nodes ---------------------------------------- inline size_t FineLayer::nnode() { return m_nnode; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t FineLayer::nne() { return m_nne; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t FineLayer::ndim() { return m_ndim; } // ---------------------------- actual number of nodes in one direction ---------------------------- inline size_t FineLayer::shape(size_t i) { assert( i >= 0 and i <= 2 ); if ( i == 0 ) return m_nelx.maxCoeff(); else if ( i == 2 ) return m_nelz.maxCoeff(); else return m_nhy .sum(); } // --------------------------------- coordinates (nodal positions) --------------------------------- inline MatD FineLayer::coor() { // allocate output MatD out(m_nnode, m_ndim); // current node, number of element layers size_t inode = 0; size_t nely = static_cast(m_nhy.size()); // y-position of each main node layer (i.e. excluding node layers for refinement/coarsening) // - allocate ColD y(nely+1); // - initialize y(0) = 0.0; // - compute for ( size_t iy = 1 ; iy < nely+1 ; ++iy ) y(iy) = y(iy-1) + m_nhy(iy-1) * m_h; // loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of nodes // -------------------------------------------------------------------------------------------- for ( size_t iy = 0 ; ; ++iy ) { // get positions along the x- and z-axis ColD x = ColD::LinSpaced(m_nelx(iy)+1, 0.0, m_Lx); ColD z = ColD::LinSpaced(m_nelz(iy)+1, 0.0, m_Lz); // add nodes of the bottom layer of this element for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) { for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { out(inode,0) = x(ix); out(inode,1) = y(iy); out(inode,2) = z(iz); ++inode; } } // stop at middle layer if ( iy == (nely-1)/2 ) break; // add extra nodes of the intermediate layer, for refinement in x-direction if ( m_refine(iy) == 0 ) { // - get position offset in x- and y-direction double dx = m_h * static_cast(m_nhx(iy)/3); double dy = m_h * static_cast(m_nhy(iy)/2); // - add nodes of the intermediate layer for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) { for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { for ( size_t j = 0 ; j < 2 ; ++j ) { out(inode,0) = x(ix) + dx * static_cast(j+1); out(inode,1) = y(iy) + dy; out(inode,2) = z(iz); ++inode; } } } } // add extra nodes of the intermediate layer, for refinement in z-direction else if ( m_refine(iy) == 2 ) { // - get position offset in y- and z-direction double dz = m_h * static_cast(m_nhz(iy)/3); double dy = m_h * static_cast(m_nhy(iy)/2); // - add nodes of the intermediate layer for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) { for ( size_t j = 0 ; j < 2 ; ++j ) { for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { out(inode,0) = x(ix); out(inode,1) = y(iy) + dy; out(inode,2) = z(iz) + dz * static_cast(j+1); ++inode; } } } } } // loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes // -------------------------------------------------------------------------------------- for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { // get positions along the x- and z-axis ColD x = ColD::LinSpaced(m_nelx(iy)+1, 0.0, m_Lx); ColD z = ColD::LinSpaced(m_nelz(iy)+1, 0.0, m_Lz); // add extra nodes of the intermediate layer, for refinement in x-direction if ( m_refine(iy) == 0 ) { // - get position offset in x- and y-direction double dx = m_h * static_cast(m_nhx(iy)/3); double dy = m_h * static_cast(m_nhy(iy)/2); // - add nodes of the intermediate layer for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) { for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { for ( size_t j = 0 ; j < 2 ; ++j ) { out(inode,0) = x(ix) + dx * static_cast(j+1); out(inode,1) = y(iy) + dy; out(inode,2) = z(iz); ++inode; } } } } // add extra nodes of the intermediate layer, for refinement in z-direction else if ( m_refine(iy) == 2 ) { // - get position offset in y- and z-direction double dz = m_h * static_cast(m_nhz(iy)/3); double dy = m_h * static_cast(m_nhy(iy)/2); // - add nodes of the intermediate layer for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) { for ( size_t j = 0 ; j < 2 ; ++j ) { for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { out(inode,0) = x(ix); out(inode,1) = y(iy) + dy; out(inode,2) = z(iz) + dz * static_cast(j+1); ++inode; } } } } // add nodes of the top layer of this element for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) { for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { out(inode,0) = x(ix ); out(inode,1) = y(iy+1); out(inode,2) = z(iz ); ++inode; } } } return out; } // ---------------------------- connectivity (node-numbers per element) ---------------------------- inline MatS FineLayer::conn() { // allocate output MatS out(m_nelem, m_nne); // current element, number of element layers, starting nodes of each node layer size_t ielem = 0; size_t nely = static_cast(m_nhy.size()); size_t bot,mid,top; // loop over all element layers for ( size_t iy = 0 ; iy < nely ; ++iy ) { // - get: starting nodes of bottom(, middle) and top layer bot = m_startNode(iy ); mid = m_startNode(iy ) + m_nnd(iy); top = m_startNode(iy+1); // - define connectivity: no coarsening/refinement if ( m_refine(iy) == -1 ) { for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) { for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { out(ielem,0) = bot + (ix ) + (iz ) * (m_nelx(iy)+1); out(ielem,1) = bot + (ix+1) + (iz ) * (m_nelx(iy)+1); out(ielem,2) = top + (ix+1) + (iz ) * (m_nelx(iy)+1); out(ielem,3) = top + (ix ) + (iz ) * (m_nelx(iy)+1); out(ielem,4) = bot + (ix ) + (iz+1) * (m_nelx(iy)+1); out(ielem,5) = bot + (ix+1) + (iz+1) * (m_nelx(iy)+1); out(ielem,6) = top + (ix+1) + (iz+1) * (m_nelx(iy)+1); out(ielem,7) = top + (ix ) + (iz+1) * (m_nelx(iy)+1); ielem++; } } } // - define connectivity: refinement along the x-direction (below the middle layer) else if ( m_refine(iy) == 0 and iy <= (nely-1)/2 ) { for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) { for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { // -- bottom element out(ielem,0) = bot + ( ix ) + (iz ) * ( m_nelx(iy)+1); out(ielem,1) = bot + ( ix+1) + (iz ) * ( m_nelx(iy)+1); out(ielem,2) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) ); out(ielem,3) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) ); out(ielem,4) = bot + ( ix ) + (iz+1) * ( m_nelx(iy)+1); out(ielem,5) = bot + ( ix+1) + (iz+1) * ( m_nelx(iy)+1); out(ielem,6) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) ); out(ielem,7) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) ); ielem++; // -- top-right element out(ielem,0) = bot + ( ix+1) + (iz ) * ( m_nelx(iy)+1); out(ielem,1) = top + (3*ix+3) + (iz ) * (3*m_nelx(iy)+1); out(ielem,2) = top + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1); out(ielem,3) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) ); out(ielem,4) = bot + ( ix+1) + (iz+1) * ( m_nelx(iy)+1); out(ielem,5) = top + (3*ix+3) + (iz+1) * (3*m_nelx(iy)+1); out(ielem,6) = top + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1); out(ielem,7) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) ); ielem++; // -- top-center element out(ielem,0) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) ); out(ielem,1) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) ); out(ielem,2) = top + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1); out(ielem,3) = top + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1); out(ielem,4) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) ); out(ielem,5) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) ); out(ielem,6) = top + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1); out(ielem,7) = top + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1); ielem++; // -- top-left element out(ielem,0) = bot + ( ix ) + (iz ) * ( m_nelx(iy)+1); out(ielem,1) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) ); out(ielem,2) = top + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1); out(ielem,3) = top + (3*ix ) + (iz ) * (3*m_nelx(iy)+1); out(ielem,4) = bot + ( ix ) + (iz+1) * ( m_nelx(iy)+1); out(ielem,5) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) ); out(ielem,6) = top + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1); out(ielem,7) = top + (3*ix ) + (iz+1) * (3*m_nelx(iy)+1); ielem++; } } } // - define connectivity: coarsening along the x-direction (above the middle layer) else if ( m_refine(iy) == 0 and iy > (nely-1)/2 ) { for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) { for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { // -- lower-left element out(ielem,0) = bot + (3*ix ) + (iz ) * (3*m_nelx(iy)+1); out(ielem,1) = bot + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1); out(ielem,2) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) ); out(ielem,3) = top + ( ix ) + (iz ) * ( m_nelx(iy)+1); out(ielem,4) = bot + (3*ix ) + (iz+1) * (3*m_nelx(iy)+1); out(ielem,5) = bot + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1); out(ielem,6) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) ); out(ielem,7) = top + ( ix ) + (iz+1) * ( m_nelx(iy)+1); ielem++; // -- lower-center element out(ielem,0) = bot + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1); out(ielem,1) = bot + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1); out(ielem,2) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) ); out(ielem,3) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) ); out(ielem,4) = bot + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1); out(ielem,5) = bot + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1); out(ielem,6) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) ); out(ielem,7) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) ); ielem++; // -- lower-right element out(ielem,0) = bot + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1); out(ielem,1) = bot + (3*ix+3) + (iz ) * (3*m_nelx(iy)+1); out(ielem,2) = top + ( ix+1) + (iz ) * ( m_nelx(iy)+1); out(ielem,3) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) ); out(ielem,4) = bot + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1); out(ielem,5) = bot + (3*ix+3) + (iz+1) * (3*m_nelx(iy)+1); out(ielem,6) = top + ( ix+1) + (iz+1) * ( m_nelx(iy)+1); out(ielem,7) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) ); ielem++; // -- upper element out(ielem,0) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) ); out(ielem,1) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) ); out(ielem,2) = top + ( ix+1) + (iz ) * ( m_nelx(iy)+1); out(ielem,3) = top + ( ix ) + (iz ) * ( m_nelx(iy)+1); out(ielem,4) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) ); out(ielem,5) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) ); out(ielem,6) = top + ( ix+1) + (iz+1) * ( m_nelx(iy)+1); out(ielem,7) = top + ( ix ) + (iz+1) * ( m_nelx(iy)+1); ielem++; } } } // - define connectivity: refinement along the z-direction (below the middle layer) else if ( m_refine(iy) == 2 and iy <= (nely-1)/2 ) { for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) { for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { // -- bottom element out(ielem,0) = bot + (ix ) + iz * (m_nelx(iy)+1); out(ielem,1) = bot + (ix+1) + iz * (m_nelx(iy)+1); out(ielem,2) = bot + (ix+1) + ( iz+1) * (m_nelx(iy)+1); out(ielem,3) = bot + (ix ) + ( iz+1) * (m_nelx(iy)+1); out(ielem,4) = mid + (ix ) + 2*iz * (m_nelx(iy)+1); out(ielem,5) = mid + (ix+1) + 2*iz * (m_nelx(iy)+1); out(ielem,6) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1); out(ielem,7) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1); ielem++; // -- top-back element out(ielem,0) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1); out(ielem,1) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1); out(ielem,2) = top + (ix+1) + (3*iz+2) * (m_nelx(iy)+1); out(ielem,3) = top + (ix ) + (3*iz+2) * (m_nelx(iy)+1); out(ielem,4) = bot + (ix ) + ( iz+1) * (m_nelx(iy)+1); out(ielem,5) = bot + (ix+1) + ( iz+1) * (m_nelx(iy)+1); out(ielem,6) = top + (ix+1) + (3*iz+3) * (m_nelx(iy)+1); out(ielem,7) = top + (ix ) + (3*iz+3) * (m_nelx(iy)+1); ielem++; // -- top-center element out(ielem,0) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1); out(ielem,1) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1); out(ielem,2) = top + (ix+1) + (3*iz+1) * (m_nelx(iy)+1); out(ielem,3) = top + (ix ) + (3*iz+1) * (m_nelx(iy)+1); out(ielem,4) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1); out(ielem,5) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1); out(ielem,6) = top + (ix+1) + (3*iz+2) * (m_nelx(iy)+1); out(ielem,7) = top + (ix ) + (3*iz+2) * (m_nelx(iy)+1); ielem++; // -- top-front element out(ielem,0) = bot + (ix ) + ( iz ) * (m_nelx(iy)+1); out(ielem,1) = bot + (ix+1) + ( iz ) * (m_nelx(iy)+1); out(ielem,2) = top + (ix+1) + (3*iz ) * (m_nelx(iy)+1); out(ielem,3) = top + (ix ) + (3*iz ) * (m_nelx(iy)+1); out(ielem,4) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1); out(ielem,5) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1); out(ielem,6) = top + (ix+1) + (3*iz+1) * (m_nelx(iy)+1); out(ielem,7) = top + (ix ) + (3*iz+1) * (m_nelx(iy)+1); ielem++; } } } // - define connectivity: coarsening along the z-direction (above the middle layer) else if ( m_refine(iy) == 2 and iy > (nely-1)/2 ) { for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) { for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { // -- bottom-front element out(ielem,0) = bot + (ix ) + (3*iz ) * (m_nelx(iy)+1); out(ielem,1) = bot + (ix+1) + (3*iz ) * (m_nelx(iy)+1); out(ielem,2) = top + (ix+1) + ( iz ) * (m_nelx(iy)+1); out(ielem,3) = top + (ix ) + ( iz ) * (m_nelx(iy)+1); out(ielem,4) = bot + (ix ) + (3*iz+1) * (m_nelx(iy)+1); out(ielem,5) = bot + (ix+1) + (3*iz+1) * (m_nelx(iy)+1); out(ielem,6) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1); out(ielem,7) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1); ielem++; // -- bottom-center element out(ielem,0) = bot + (ix ) + (3*iz+1) * (m_nelx(iy)+1); out(ielem,1) = bot + (ix+1) + (3*iz+1) * (m_nelx(iy)+1); out(ielem,2) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1); out(ielem,3) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1); out(ielem,4) = bot + (ix ) + (3*iz+2) * (m_nelx(iy)+1); out(ielem,5) = bot + (ix+1) + (3*iz+2) * (m_nelx(iy)+1); out(ielem,6) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1); out(ielem,7) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1); ielem++; // -- bottom-back element out(ielem,0) = bot + (ix ) + (3*iz+2) * (m_nelx(iy)+1); out(ielem,1) = bot + (ix+1) + (3*iz+2) * (m_nelx(iy)+1); out(ielem,2) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1); out(ielem,3) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1); out(ielem,4) = bot + (ix ) + (3*iz+3) * (m_nelx(iy)+1); out(ielem,5) = bot + (ix+1) + (3*iz+3) * (m_nelx(iy)+1); out(ielem,6) = top + (ix+1) + ( iz+1) * (m_nelx(iy)+1); out(ielem,7) = top + (ix ) + ( iz+1) * (m_nelx(iy)+1); ielem++; // -- top element out(ielem,0) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1); out(ielem,1) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1); out(ielem,2) = top + (ix+1) + ( iz ) * (m_nelx(iy)+1); out(ielem,3) = top + (ix ) + ( iz ) * (m_nelx(iy)+1); out(ielem,4) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1); out(ielem,5) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1); out(ielem,6) = top + (ix+1) + ( iz+1) * (m_nelx(iy)+1); out(ielem,7) = top + (ix ) + ( iz+1) * (m_nelx(iy)+1); ielem++; } } } } return out; } +// ------------------------------ element numbers of the middle layer ------------------------------ + +inline ColS FineLayer::elementsMiddleLayer() +{ + // number of element layers in y-direction, the index of the middle layer + size_t nely = static_cast(m_nhy.size()); + size_t iy = (nely-1)/2; + + ColS out(m_nelx(iy)*m_nelz(iy)); + + for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) + for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) + out(ix+iz*m_nelx(iy)) = m_startElem(iy) + ix + iz*m_nelx(iy); + + return out; +} + // ------------------------------ node-numbers along the front plane ------------------------------- inline ColS FineLayer::nodesFront() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) { if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1; else n += m_nelx(iy) + 1; } // - top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1; else n += m_nelx(iy) + 1; } // allocate node-list ColS out(n); // initialize counter: current index in the node-list "out" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) { // -- bottom node layer for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { out(j) = m_startNode(iy) + ix; ++j; } // -- refinement layer if ( m_refine(iy) == 0 ) { for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy) + ix + m_nnd(iy); ++j; } } } // top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { // -- refinement layer if ( m_refine(iy) == 0 ) { for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy) + ix + m_nnd(iy); ++j; } } // -- top node layer for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { out(j) = m_startNode(iy+1) + ix; ++j; } } return out; } // ------------------------------- node-numbers along the back plane ------------------------------- inline ColS FineLayer::nodesBack() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) { if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1; else n += m_nelx(iy) + 1; } // - top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1; else n += m_nelx(iy) + 1; } // allocate node-list ColS out(n); // initialize counter: current index in the node-list "out" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) { // -- bottom node layer for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { out(j) = m_startNode(iy) + ix + (m_nelx(iy)+1)*m_nelz(iy); ++j; } // -- refinement layer if ( m_refine(iy) == 0 ) { for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy); ++j; } } } // top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { // -- refinement layer if ( m_refine(iy) == 0 ) { for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy); ++j; } } // -- top node layer for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { out(j) = m_startNode(iy+1) + ix + (m_nelx(iy)+1)*m_nelz(iy); ++j; } } return out; } // ------------------------------- node-numbers along the left plane ------------------------------- inline ColS FineLayer::nodesLeft() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) { if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1; else n += m_nelz(iy) + 1; } // - top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1; else n += m_nelz(iy) + 1; } // allocate node-list ColS out(n); // initialize counter: current index in the node-list "out" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) { // -- bottom node layer for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) { out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1); ++j; } // -- refinement layer if ( m_refine(iy) == 2 ) { for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy); ++j; } } } // top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { // -- refinement layer if ( m_refine(iy) == 2 ) { for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy); ++j; } } // -- top node layer for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) { out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1); ++j; } } return out; } // ------------------------------ node-numbers along the right plane ------------------------------- inline ColS FineLayer::nodesRight() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) { if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1; else n += m_nelz(iy) + 1; } // - top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1; else n += m_nelz(iy) + 1; } // allocate node-list ColS out(n); // initialize counter: current index in the node-list "out" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) { // -- bottom node layer for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) { out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy); ++j; } // -- refinement layer if ( m_refine(iy) == 2 ) { for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy); ++j; } } } // top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { // -- refinement layer if ( m_refine(iy) == 2 ) { for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy); ++j; } } // -- top node layer for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) { out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1) + m_nelx(iy); ++j; } } return out; } // ------------------------------ node-numbers along the bottom plane ------------------------------ inline ColS FineLayer::nodesBottom() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // allocate node list ColS out(m_nnd(nely)); // counter size_t j = 0; // fill node list for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix ) { for ( size_t iz = 0 ; iz < m_nelz(0)+1 ; ++iz ) { out(j) = m_startNode(0) + ix + iz * (m_nelx(0)+1); ++j; } } return out; } // ------------------------------- node-numbers along the top plane -------------------------------- inline ColS FineLayer::nodesTop() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // allocate node list ColS out(m_nnd(nely)); // counter size_t j = 0; // fill node list for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix ) { for ( size_t iz = 0 ; iz < m_nelz(nely-1)+1 ; ++iz ) { out(j) = m_startNode(nely) + ix + iz * (m_nelx(nely-1)+1); ++j; } } return out; } // ------------------------------- node-numbers along the front face ------------------------------- inline ColS FineLayer::nodesFrontFace() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) { if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1; else n += m_nelx(iy) - 1; } // - top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) { if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1; else n += m_nelx(iy) - 1; } // allocate node-list ColS out(n); // initialize counter: current index in the node-list "out" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) { // -- bottom node layer for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy) + ix; ++j; } // -- refinement layer if ( m_refine(iy) == 0 ) { for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy) + ix + m_nnd(iy); ++j; } } } // top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) { // -- refinement layer if ( m_refine(iy) == 0 ) { for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy) + ix + m_nnd(iy); ++j; } } // -- top node layer for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy+1) + ix; ++j; } } return out; } // ------------------------------- node-numbers along the back face -------------------------------- inline ColS FineLayer::nodesBackFace() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) { if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1; else n += m_nelx(iy) - 1; } // - top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) { if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1; else n += m_nelx(iy) - 1; } // allocate node-list ColS out(n); // initialize counter: current index in the node-list "out" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) { // -- bottom node layer for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy) + ix + (m_nelx(iy)+1)*m_nelz(iy); ++j; } // -- refinement layer if ( m_refine(iy) == 0 ) { for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy); ++j; } } } // top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) { // -- refinement layer if ( m_refine(iy) == 0 ) { for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy); ++j; } } // -- top node layer for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) { out(j) = m_startNode(iy+1) + ix + (m_nelx(iy)+1)*m_nelz(iy); ++j; } } return out; } // ------------------------------- node-numbers along the left face -------------------------------- inline ColS FineLayer::nodesLeftFace() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) { if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1; else n += m_nelz(iy) - 1; } // - top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) { if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1; else n += m_nelz(iy) - 1; } // allocate node-list ColS out(n); // initialize counter: current index in the node-list "out" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) { // -- bottom node layer for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1); ++j; } // -- refinement layer if ( m_refine(iy) == 2 ) { for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy); ++j; } } } // top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) { // -- refinement layer if ( m_refine(iy) == 2 ) { for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy); ++j; } } // -- top node layer for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1); ++j; } } return out; } // ------------------------------- node-numbers along the right face ------------------------------- inline ColS FineLayer::nodesRightFace() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // number of boundary nodes // - initialize size_t n = 0; // - bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) { if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1; else n += m_nelz(iy) - 1; } // - top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) { if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1; else n += m_nelz(iy) - 1; } // allocate node-list ColS out(n); // initialize counter: current index in the node-list "out" size_t j = 0; // bottom half: bottom node layer (+ middle node layer) for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) { // -- bottom node layer for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy); ++j; } // -- refinement layer if ( m_refine(iy) == 2 ) { for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy); ++j; } } } // top half: (middle node layer +) top node layer for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) { // -- refinement layer if ( m_refine(iy) == 2 ) { for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy); ++j; } } // -- top node layer for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) { out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1) + m_nelx(iy); ++j; } } return out; } // ------------------------------ node-numbers along the bottom face ------------------------------- inline ColS FineLayer::nodesBottomFace() { // allocate node list ColS out((m_nelx(0)-1)*(m_nelz(0)-1)); // counter size_t j = 0; // fill node list for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix ) { for ( size_t iz = 1 ; iz < m_nelz(0) ; ++iz ) { out(j) = m_startNode(0) + ix + iz * (m_nelx(0)+1); ++j; } } return out; } // -------------------------------- node-numbers along the top face -------------------------------- inline ColS FineLayer::nodesTopFace() { // number of element layers in y-direction size_t nely = static_cast(m_nhy.size()); // allocate node list ColS out((m_nelx(nely-1)-1)*(m_nelz(nely-1)-1)); // counter size_t j = 0; // fill node list for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix ) { for ( size_t iz = 1 ; iz < m_nelz(nely-1) ; ++iz ) { out(j) = m_startNode(nely) + ix + iz * (m_nelx(nely-1)+1); ++j; } } return out; } // --------------------------- node-numbers along the front-bottom edge ---------------------------- inline ColS FineLayer::nodesFrontBottomEdge() { ColS out(m_nelx(0)+1); for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix ) out(ix) = m_startNode(0) + ix; return out; } // ----------------------------- node-numbers along the front-top edge ----------------------------- inline ColS FineLayer::nodesFrontTopEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(m_nelx(nely-1)+1); for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix ) out(ix) = m_startNode(nely) + ix; return out; } // ---------------------------- node-numbers along the front-left edge ----------------------------- inline ColS FineLayer::nodesFrontLeftEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(nely+1); for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) out(iy) = m_startNode(iy); for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) out(iy+1) = m_startNode(iy+1); return out; } // ---------------------------- node-numbers along the front-right edge ---------------------------- inline ColS FineLayer::nodesFrontRightEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(nely+1); for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) out(iy) = m_startNode(iy) + m_nelx(iy); for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) out(iy+1) = m_startNode(iy+1) + m_nelx(iy); return out; } // ---------------------------- node-numbers along the back-bottom edge ---------------------------- inline ColS FineLayer::nodesBackBottomEdge() { ColS out(m_nelx(0)+1); for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix ) out(ix) = m_startNode(0) + ix + (m_nelx(0)+1)*(m_nelz(0)); return out; } // ----------------------------- node-numbers along the back-top edge ------------------------------ inline ColS FineLayer::nodesBackTopEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(m_nelx(nely-1)+1); for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix ) out(ix) = m_startNode(nely) + ix + (m_nelx(nely-1)+1)*(m_nelz(nely-1)); return out; } // ----------------------------- node-numbers along the back-left edge ----------------------------- inline ColS FineLayer::nodesBackLeftEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(nely+1); for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) out(iy) = m_startNode(iy) + (m_nelx(iy)+1)*(m_nelz(iy)); for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) out(iy+1) = m_startNode(iy+1) + (m_nelx(iy)+1)*(m_nelz(iy)); return out; } // ---------------------------- node-numbers along the back-right edge ----------------------------- inline ColS FineLayer::nodesBackRightEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(nely+1); for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) out(iy) = m_startNode(iy) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy)); for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) out(iy+1) = m_startNode(iy+1) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy)); return out; } // ---------------------------- node-numbers along the bottom-left edge ---------------------------- inline ColS FineLayer::nodesBottomLeftEdge() { ColS out(m_nelz(0)+1); for ( size_t iz = 0 ; iz < m_nelz(0)+1 ; ++iz ) out(iz) = m_startNode(0) + iz * (m_nelx(0)+1); return out; } // --------------------------- node-numbers along the bottom-right edge ---------------------------- inline ColS FineLayer::nodesBottomRightEdge() { ColS out(m_nelz(0)+1); for ( size_t iz = 0 ; iz < m_nelz(0)+1 ; ++iz ) out(iz) = m_startNode(0) + m_nelx(0) + iz * (m_nelx(0)+1); return out; } // ----------------------------- node-numbers along the top-left edge ------------------------------ inline ColS FineLayer::nodesTopLeftEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(m_nelz(nely-1)+1); for ( size_t iz = 0 ; iz < m_nelz(nely-1)+1 ; ++iz ) out(iz) = m_startNode(nely) + iz * (m_nelx(nely-1)+1); return out; } // ----------------------------- node-numbers along the top-right edge ----------------------------- inline ColS FineLayer::nodesTopRightEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(m_nelz(nely-1)+1); for ( size_t iz = 0 ; iz < m_nelz(nely-1)+1 ; ++iz ) out(iz) = m_startNode(nely) + m_nelx(nely-1) + iz * (m_nelx(nely-1)+1); return out; } // -------------------------------------------- aliases -------------------------------------------- inline ColS FineLayer::nodesBottomFrontEdge() { return nodesFrontBottomEdge(); } inline ColS FineLayer::nodesBottomBackEdge() { return nodesBackBottomEdge(); } inline ColS FineLayer::nodesTopFrontEdge() { return nodesFrontTopEdge(); } inline ColS FineLayer::nodesTopBackEdge() { return nodesBackTopEdge(); } inline ColS FineLayer::nodesLeftBottomEdge() { return nodesBottomLeftEdge(); } inline ColS FineLayer::nodesLeftFrontEdge() { return nodesFrontLeftEdge(); } inline ColS FineLayer::nodesLeftBackEdge() { return nodesBackLeftEdge(); } inline ColS FineLayer::nodesLeftTopEdge() { return nodesTopLeftEdge(); } inline ColS FineLayer::nodesRightBottomEdge() { return nodesBottomRightEdge(); } inline ColS FineLayer::nodesRightTopEdge() { return nodesTopRightEdge(); } inline ColS FineLayer::nodesRightFrontEdge() { return nodesFrontRightEdge(); } inline ColS FineLayer::nodesRightBackEdge() { return nodesBackRightEdge(); } // ------------------- node-numbers along the front-bottom edge, without corners ------------------- inline ColS FineLayer::nodesFrontBottomOpenEdge() { ColS out(m_nelx(0)-1); for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix ) out(ix-1) = m_startNode(0) + ix; return out; } // -------------------- node-numbers along the front-top edge, without corners --------------------- inline ColS FineLayer::nodesFrontTopOpenEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(m_nelx(nely-1)-1); for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix ) out(ix-1) = m_startNode(nely) + ix; return out; } // -------------------- node-numbers along the front-left edge, without corners -------------------- inline ColS FineLayer::nodesFrontLeftOpenEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(nely-1); for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) out(iy-1) = m_startNode(iy); for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) out(iy) = m_startNode(iy+1); return out; } // ------------------- node-numbers along the front-right edge, without corners -------------------- inline ColS FineLayer::nodesFrontRightOpenEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(nely-1); for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) out(iy-1) = m_startNode(iy) + m_nelx(iy); for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) out(iy) = m_startNode(iy+1) + m_nelx(iy); return out; } // ------------------- node-numbers along the back-bottom edge, without corners -------------------- inline ColS FineLayer::nodesBackBottomOpenEdge() { ColS out(m_nelx(0)-1); for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix ) out(ix-1) = m_startNode(0) + ix + (m_nelx(0)+1)*(m_nelz(0)); return out; } // --------------------- node-numbers along the back-top edge, without corners --------------------- inline ColS FineLayer::nodesBackTopOpenEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(m_nelx(nely-1)-1); for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix ) out(ix-1) = m_startNode(nely) + ix + (m_nelx(nely-1)+1)*(m_nelz(nely-1)); return out; } // -------------------- node-numbers along the back-left edge, without corners --------------------- inline ColS FineLayer::nodesBackLeftOpenEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(nely-1); for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) out(iy-1) = m_startNode(iy) + (m_nelx(iy)+1)*(m_nelz(iy)); for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) out(iy) = m_startNode(iy+1) + (m_nelx(iy)+1)*(m_nelz(iy)); return out; } // -------------------- node-numbers along the back-right edge, without corners -------------------- inline ColS FineLayer::nodesBackRightOpenEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(nely-1); for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) out(iy-1) = m_startNode(iy) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy)); for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) out(iy) = m_startNode(iy+1) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy)); return out; } // ------------------- node-numbers along the bottom-left edge, without corners -------------------- inline ColS FineLayer::nodesBottomLeftOpenEdge() { ColS out(m_nelz(0)-1); for ( size_t iz = 1 ; iz < m_nelz(0) ; ++iz ) out(iz-1) = m_startNode(0) + iz * (m_nelx(0)+1); return out; } // ------------------- node-numbers along the bottom-right edge, without corners ------------------- inline ColS FineLayer::nodesBottomRightOpenEdge() { ColS out(m_nelz(0)-1); for ( size_t iz = 1 ; iz < m_nelz(0) ; ++iz ) out(iz-1) = m_startNode(0) + m_nelx(0) + iz * (m_nelx(0)+1); return out; } // --------------------- node-numbers along the top-left edge, without corners --------------------- inline ColS FineLayer::nodesTopLeftOpenEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(m_nelz(nely-1)-1); for ( size_t iz = 1 ; iz < m_nelz(nely-1) ; ++iz ) out(iz-1) = m_startNode(nely) + iz * (m_nelx(nely-1)+1); return out; } // -------------------- node-numbers along the top-right edge, without corners --------------------- inline ColS FineLayer::nodesTopRightOpenEdge() { size_t nely = static_cast(m_nhy.size()); ColS out(m_nelz(nely-1)-1); for ( size_t iz = 1 ; iz < m_nelz(nely-1) ; ++iz ) out(iz-1) = m_startNode(nely) + m_nelx(nely-1) + iz * (m_nelx(nely-1)+1); return out; } // -------------------------------------------- aliases -------------------------------------------- inline ColS FineLayer::nodesBottomFrontOpenEdge() { return nodesFrontBottomOpenEdge(); } inline ColS FineLayer::nodesBottomBackOpenEdge() { return nodesBackBottomOpenEdge(); } inline ColS FineLayer::nodesTopFrontOpenEdge() { return nodesFrontTopOpenEdge(); } inline ColS FineLayer::nodesTopBackOpenEdge() { return nodesBackTopOpenEdge(); } inline ColS FineLayer::nodesLeftBottomOpenEdge() { return nodesBottomLeftOpenEdge(); } inline ColS FineLayer::nodesLeftFrontOpenEdge() { return nodesFrontLeftOpenEdge(); } inline ColS FineLayer::nodesLeftBackOpenEdge() { return nodesBackLeftOpenEdge(); } inline ColS FineLayer::nodesLeftTopOpenEdge() { return nodesTopLeftOpenEdge(); } inline ColS FineLayer::nodesRightBottomOpenEdge() { return nodesBottomRightOpenEdge(); } inline ColS FineLayer::nodesRightTopOpenEdge() { return nodesTopRightOpenEdge(); } inline ColS FineLayer::nodesRightFrontOpenEdge() { return nodesFrontRightOpenEdge(); } inline ColS FineLayer::nodesRightBackOpenEdge() { return nodesBackRightOpenEdge(); } // -------------------------- node-number of the front-bottom-left corner -------------------------- inline size_t FineLayer::nodesFrontBottomLeftCorner() { return m_startNode(0); } // ------------------------- node-number of the front-bottom-right corner -------------------------- inline size_t FineLayer::nodesFrontBottomRightCorner() { return m_startNode(0) + m_nelx(0); } // --------------------------- node-number of the front-top-left corner ---------------------------- inline size_t FineLayer::nodesFrontTopLeftCorner() { size_t nely = static_cast(m_nhy.size()); return m_startNode(nely); } // --------------------------- node-number of the front-top-right corner --------------------------- inline size_t FineLayer::nodesFrontTopRightCorner() { size_t nely = static_cast(m_nhy.size()); return m_startNode(nely) + m_nelx(nely-1); } // -------------------------- node-number of the back-bottom-left corner --------------------------- inline size_t FineLayer::nodesBackBottomLeftCorner() { return m_startNode(0) + (m_nelx(0)+1)*(m_nelz(0)); } // -------------------------- node-number of the back-bottom-right corner -------------------------- inline size_t FineLayer::nodesBackBottomRightCorner() { return m_startNode(0) + m_nelx(0) + (m_nelx(0)+1)*(m_nelz(0)); } // ---------------------------- node-number of the back-top-left corner ---------------------------- inline size_t FineLayer::nodesBackTopLeftCorner() { size_t nely = static_cast(m_nhy.size()); return m_startNode(nely) + (m_nelx(nely-1)+1)*(m_nelz(nely-1)); } // --------------------------- node-number of the back-top-right corner ---------------------------- inline size_t FineLayer::nodesBackTopRightCorner() { size_t nely = static_cast(m_nhy.size()); return m_startNode(nely) + m_nelx(nely-1) + (m_nelx(nely-1)+1)*(m_nelz(nely-1)); } // -------------------------------------------- aliases -------------------------------------------- inline size_t FineLayer::nodesFrontLeftBottomCorner() { return nodesFrontBottomLeftCorner(); } inline size_t FineLayer::nodesBottomFrontLeftCorner() { return nodesFrontBottomLeftCorner(); } inline size_t FineLayer::nodesBottomLeftFrontCorner() { return nodesFrontBottomLeftCorner(); } inline size_t FineLayer::nodesLeftFrontBottomCorner() { return nodesFrontBottomLeftCorner(); } inline size_t FineLayer::nodesLeftBottomFrontCorner() { return nodesFrontBottomLeftCorner(); } inline size_t FineLayer::nodesFrontRightBottomCorner() { return nodesFrontBottomRightCorner(); } inline size_t FineLayer::nodesBottomFrontRightCorner() { return nodesFrontBottomRightCorner(); } inline size_t FineLayer::nodesBottomRightFrontCorner() { return nodesFrontBottomRightCorner(); } inline size_t FineLayer::nodesRightFrontBottomCorner() { return nodesFrontBottomRightCorner(); } inline size_t FineLayer::nodesRightBottomFrontCorner() { return nodesFrontBottomRightCorner(); } inline size_t FineLayer::nodesFrontLeftTopCorner() { return nodesFrontTopLeftCorner(); } inline size_t FineLayer::nodesTopFrontLeftCorner() { return nodesFrontTopLeftCorner(); } inline size_t FineLayer::nodesTopLeftFrontCorner() { return nodesFrontTopLeftCorner(); } inline size_t FineLayer::nodesLeftFrontTopCorner() { return nodesFrontTopLeftCorner(); } inline size_t FineLayer::nodesLeftTopFrontCorner() { return nodesFrontTopLeftCorner(); } inline size_t FineLayer::nodesFrontRightTopCorner() { return nodesFrontTopRightCorner(); } inline size_t FineLayer::nodesTopFrontRightCorner() { return nodesFrontTopRightCorner(); } inline size_t FineLayer::nodesTopRightFrontCorner() { return nodesFrontTopRightCorner(); } inline size_t FineLayer::nodesRightFrontTopCorner() { return nodesFrontTopRightCorner(); } inline size_t FineLayer::nodesRightTopFrontCorner() { return nodesFrontTopRightCorner(); } inline size_t FineLayer::nodesBackLeftBottomCorner() { return nodesBackBottomLeftCorner(); } inline size_t FineLayer::nodesBottomBackLeftCorner() { return nodesBackBottomLeftCorner(); } inline size_t FineLayer::nodesBottomLeftBackCorner() { return nodesBackBottomLeftCorner(); } inline size_t FineLayer::nodesLeftBackBottomCorner() { return nodesBackBottomLeftCorner(); } inline size_t FineLayer::nodesLeftBottomBackCorner() { return nodesBackBottomLeftCorner(); } inline size_t FineLayer::nodesBackRightBottomCorner() { return nodesBackBottomRightCorner(); } inline size_t FineLayer::nodesBottomBackRightCorner() { return nodesBackBottomRightCorner(); } inline size_t FineLayer::nodesBottomRightBackCorner() { return nodesBackBottomRightCorner(); } inline size_t FineLayer::nodesRightBackBottomCorner() { return nodesBackBottomRightCorner(); } inline size_t FineLayer::nodesRightBottomBackCorner() { return nodesBackBottomRightCorner(); } inline size_t FineLayer::nodesBackLeftTopCorner() { return nodesBackTopLeftCorner(); } inline size_t FineLayer::nodesTopBackLeftCorner() { return nodesBackTopLeftCorner(); } inline size_t FineLayer::nodesTopLeftBackCorner() { return nodesBackTopLeftCorner(); } inline size_t FineLayer::nodesLeftBackTopCorner() { return nodesBackTopLeftCorner(); } inline size_t FineLayer::nodesLeftTopBackCorner() { return nodesBackTopLeftCorner(); } inline size_t FineLayer::nodesBackRightTopCorner() { return nodesBackTopRightCorner(); } inline size_t FineLayer::nodesTopBackRightCorner() { return nodesBackTopRightCorner(); } inline size_t FineLayer::nodesTopRightBackCorner() { return nodesBackTopRightCorner(); } inline size_t FineLayer::nodesRightBackTopCorner() { return nodesBackTopRightCorner(); } inline size_t FineLayer::nodesRightTopBackCorner() { return nodesBackTopRightCorner(); } // ------------------------------ node-numbers of periodic node-pairs ------------------------------ inline MatS FineLayer::nodesPeriodic() { // faces ColS fro = nodesFrontFace(); ColS bck = nodesBackFace(); ColS lft = nodesLeftFace(); ColS rgt = nodesRightFace(); ColS bot = nodesBottomFace(); ColS top = nodesTopFace(); // edges ColS froBot = nodesFrontBottomOpenEdge(); ColS froTop = nodesFrontTopOpenEdge(); ColS froLft = nodesFrontLeftOpenEdge(); ColS froRgt = nodesFrontRightOpenEdge(); ColS bckBot = nodesBackBottomOpenEdge(); ColS bckTop = nodesBackTopOpenEdge(); ColS bckLft = nodesBackLeftOpenEdge(); ColS bckRgt = nodesBackRightOpenEdge(); ColS botLft = nodesBottomLeftOpenEdge(); ColS botRgt = nodesBottomRightOpenEdge(); ColS topLft = nodesTopLeftOpenEdge(); ColS topRgt = nodesTopRightOpenEdge(); // allocate nodal ties // - number of tying per category size_t tface = fro.size() + lft.size() + bot.size(); size_t tedge = 3*froBot.size() + 3*froLft.size() + 3*botLft.size(); size_t tnode = 7; // - allocate MatS out(tface+tedge+tnode, 2); // counter size_t i = 0; // tie all corners to one corner out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontBottomRightCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomRightCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomLeftCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopLeftCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopRightCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopRightCorner(); ++i; out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopLeftCorner(); ++i; // tie all corresponding edges to each other (exclude corners) for ( auto j = 0 ; j(nodePer.rows()); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs for ( size_t i = 0 ; i < nper ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) out(nodePer(i,1),j) = out(nodePer(i,0),j); // renumber "out" to be sequential return GooseFEM::Mesh::renumber(out); } // ================================================================================================= }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/MeshHex8.h b/src/GooseFEM/MeshHex8.h index 2a7b3f6..8290448 100644 --- a/src/GooseFEM/MeshHex8.h +++ b/src/GooseFEM/MeshHex8.h @@ -1,327 +1,329 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESHHEX8_H #define GOOSEFEM_MESHHEX8_H // ------------------------------------------------------------------------------------------------- #include "Mesh.h" // ===================================== GooseFEM::Mesh::Hex8 ====================================== namespace GooseFEM { namespace Mesh { namespace Hex8 { // ========================================== REGULAR MESH ========================================= class Regular { private: double m_h; // elementary element edge-size (in all directions) - size_t m_nelx; // number of element x-direction (length == "m_nelx * m_h") - size_t m_nely; // number of element y-direction (length == "m_nely * m_h") - size_t m_nelz; // number of element z-direction (length == "m_nely * m_h") + size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h") + size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h") + size_t m_nelz; // number of elements in z-direction (length == "m_nely * m_h") size_t m_nelem; // number of elements size_t m_nnode; // number of nodes size_t m_nne=8; // number of nodes-per-element size_t m_ndim=3; // number of dimensions public: - // mesh with "nelx*nely*nelz" 'pixels' of edge size "h" + // mesh with "nelx*nely*nelz" 'elements' of edge size "h" Regular(size_t nelx, size_t nely, size_t nelz, double h=1.); // sizes size_t nelem(); // number of elements size_t nnode(); // number of nodes size_t nne(); // number of nodes-per-element size_t ndim(); // number of dimensions // mesh MatD coor(); // nodal positions [nnode ,ndim] MatS conn(); // connectivity [nelem ,nne ] // boundary nodes: planes ColS nodesFront(); // node-numbers along the front plane ColS nodesBack(); // node-numbers along the back plane ColS nodesLeft(); // node-numbers along the left plane ColS nodesRight(); // node-numbers along the right plane ColS nodesBottom(); // node-numbers along the bottom plane ColS nodesTop(); // node-numbers along the top plane // boundary nodes: faces ColS nodesFrontFace(); // node-numbers along the front face ColS nodesBackFace(); // node-numbers along the back face ColS nodesLeftFace(); // node-numbers along the left face ColS nodesRightFace(); // node-numbers along the right face ColS nodesBottomFace(); // node-numbers along the bottom face ColS nodesTopFace(); // node-numbers along the top face // boundary nodes: edges ColS nodesFrontBottomEdge(); // node-numbers along the front - bottom edge ColS nodesFrontTopEdge(); // node-numbers along the front - top edge ColS nodesFrontLeftEdge(); // node-numbers along the front - left edge ColS nodesFrontRightEdge(); // node-numbers along the front - right edge ColS nodesBackBottomEdge(); // node-numbers along the back - bottom edge ColS nodesBackTopEdge(); // node-numbers along the back - top edge ColS nodesBackLeftEdge(); // node-numbers along the back - left edge ColS nodesBackRightEdge(); // node-numbers along the back - right edge ColS nodesBottomLeftEdge(); // node-numbers along the bottom - left edge ColS nodesBottomRightEdge(); // node-numbers along the bottom - right edge ColS nodesTopLeftEdge(); // node-numbers along the top - left edge ColS nodesTopRightEdge(); // node-numbers along the top - right edge // boundary nodes: edges (aliases) ColS nodesBottomFrontEdge(); // alias, see above: nodesFrontBottomEdge ColS nodesBottomBackEdge(); // alias, see above: nodesBackBottomEdge ColS nodesTopFrontEdge(); // alias, see above: nodesFrontTopEdge ColS nodesTopBackEdge(); // alias, see above: nodesBackTopEdge ColS nodesLeftBottomEdge(); // alias, see above: nodesBottomLeftEdge ColS nodesLeftFrontEdge(); // alias, see above: nodesFrontLeftEdge ColS nodesLeftBackEdge(); // alias, see above: nodesBackLeftEdge ColS nodesLeftTopEdge(); // alias, see above: nodesTopLeftEdge ColS nodesRightBottomEdge(); // alias, see above: nodesBottomRightEdge ColS nodesRightTopEdge(); // alias, see above: nodesTopRightEdge ColS nodesRightFrontEdge(); // alias, see above: nodesFrontRightEdge ColS nodesRightBackEdge(); // alias, see above: nodesBackRightEdge // boundary nodes: edges, without corners ColS nodesFrontBottomOpenEdge(); // node-numbers along the front - bottom edge ColS nodesFrontTopOpenEdge(); // node-numbers along the front - top edge ColS nodesFrontLeftOpenEdge(); // node-numbers along the front - left edge ColS nodesFrontRightOpenEdge(); // node-numbers along the front - right edge ColS nodesBackBottomOpenEdge(); // node-numbers along the back - bottom edge ColS nodesBackTopOpenEdge(); // node-numbers along the back - top edge ColS nodesBackLeftOpenEdge(); // node-numbers along the back - left edge ColS nodesBackRightOpenEdge(); // node-numbers along the back - right edge ColS nodesBottomLeftOpenEdge(); // node-numbers along the bottom - left edge ColS nodesBottomRightOpenEdge(); // node-numbers along the bottom - right edge ColS nodesTopLeftOpenEdge(); // node-numbers along the top - left edge ColS nodesTopRightOpenEdge(); // node-numbers along the top - right edge // boundary nodes: edges, without corners (aliases) ColS nodesBottomFrontOpenEdge(); // alias, see above: nodesFrontBottomOpenEdge ColS nodesBottomBackOpenEdge(); // alias, see above: nodesBackBottomOpenEdge ColS nodesTopFrontOpenEdge(); // alias, see above: nodesFrontTopOpenEdge ColS nodesTopBackOpenEdge(); // alias, see above: nodesBackTopOpenEdge ColS nodesLeftBottomOpenEdge(); // alias, see above: nodesBottomLeftOpenEdge ColS nodesLeftFrontOpenEdge(); // alias, see above: nodesFrontLeftOpenEdge ColS nodesLeftBackOpenEdge(); // alias, see above: nodesBackLeftOpenEdge ColS nodesLeftTopOpenEdge(); // alias, see above: nodesTopLeftOpenEdge ColS nodesRightBottomOpenEdge(); // alias, see above: nodesBottomRightOpenEdge ColS nodesRightTopOpenEdge(); // alias, see above: nodesTopRightOpenEdge ColS nodesRightFrontOpenEdge(); // alias, see above: nodesFrontRightOpenEdge ColS nodesRightBackOpenEdge(); // alias, see above: nodesBackRightOpenEdge // boundary nodes: corners size_t nodesFrontBottomLeftCorner(); // node-number of the front - bottom - left corner size_t nodesFrontBottomRightCorner(); // node-number of the front - bottom - right corner size_t nodesFrontTopLeftCorner(); // node-number of the front - top - left corner size_t nodesFrontTopRightCorner(); // node-number of the front - top - right corner size_t nodesBackBottomLeftCorner(); // node-number of the back - bottom - left corner size_t nodesBackBottomRightCorner(); // node-number of the back - bottom - right corner size_t nodesBackTopLeftCorner(); // node-number of the back - top - left corner size_t nodesBackTopRightCorner(); // node-number of the back - top - right corner // boundary nodes: corners (aliases) size_t nodesFrontLeftBottomCorner(); // alias, see above: nodesFrontBottomLeftCorner size_t nodesBottomFrontLeftCorner(); // alias, see above: nodesFrontBottomLeftCorner size_t nodesBottomLeftFrontCorner(); // alias, see above: nodesFrontBottomLeftCorner size_t nodesLeftFrontBottomCorner(); // alias, see above: nodesFrontBottomLeftCorner size_t nodesLeftBottomFrontCorner(); // alias, see above: nodesFrontBottomLeftCorner size_t nodesFrontRightBottomCorner(); // alias, see above: nodesFrontBottomRightCorner size_t nodesBottomFrontRightCorner(); // alias, see above: nodesFrontBottomRightCorner size_t nodesBottomRightFrontCorner(); // alias, see above: nodesFrontBottomRightCorner size_t nodesRightFrontBottomCorner(); // alias, see above: nodesFrontBottomRightCorner size_t nodesRightBottomFrontCorner(); // alias, see above: nodesFrontBottomRightCorner size_t nodesFrontLeftTopCorner(); // alias, see above: nodesFrontTopLeftCorner size_t nodesTopFrontLeftCorner(); // alias, see above: nodesFrontTopLeftCorner size_t nodesTopLeftFrontCorner(); // alias, see above: nodesFrontTopLeftCorner size_t nodesLeftFrontTopCorner(); // alias, see above: nodesFrontTopLeftCorner size_t nodesLeftTopFrontCorner(); // alias, see above: nodesFrontTopLeftCorner size_t nodesFrontRightTopCorner(); // alias, see above: nodesFrontTopRightCorner size_t nodesTopFrontRightCorner(); // alias, see above: nodesFrontTopRightCorner size_t nodesTopRightFrontCorner(); // alias, see above: nodesFrontTopRightCorner size_t nodesRightFrontTopCorner(); // alias, see above: nodesFrontTopRightCorner size_t nodesRightTopFrontCorner(); // alias, see above: nodesFrontTopRightCorner size_t nodesBackLeftBottomCorner(); // alias, see above: nodesBackBottomLeftCorner size_t nodesBottomBackLeftCorner(); // alias, see above: nodesBackBottomLeftCorner size_t nodesBottomLeftBackCorner(); // alias, see above: nodesBackBottomLeftCorner size_t nodesLeftBackBottomCorner(); // alias, see above: nodesBackBottomLeftCorner size_t nodesLeftBottomBackCorner(); // alias, see above: nodesBackBottomLeftCorner size_t nodesBackRightBottomCorner(); // alias, see above: nodesBackBottomRightCorner size_t nodesBottomBackRightCorner(); // alias, see above: nodesBackBottomRightCorner size_t nodesBottomRightBackCorner(); // alias, see above: nodesBackBottomRightCorner size_t nodesRightBackBottomCorner(); // alias, see above: nodesBackBottomRightCorner size_t nodesRightBottomBackCorner(); // alias, see above: nodesBackBottomRightCorner size_t nodesBackLeftTopCorner(); // alias, see above: nodesBackTopLeftCorner size_t nodesTopBackLeftCorner(); // alias, see above: nodesBackTopLeftCorner size_t nodesTopLeftBackCorner(); // alias, see above: nodesBackTopLeftCorner size_t nodesLeftBackTopCorner(); // alias, see above: nodesBackTopLeftCorner size_t nodesLeftTopBackCorner(); // alias, see above: nodesBackTopLeftCorner size_t nodesBackRightTopCorner(); // alias, see above: nodesBackTopRightCorner size_t nodesTopBackRightCorner(); // alias, see above: nodesBackTopRightCorner size_t nodesTopRightBackCorner(); // alias, see above: nodesBackTopRightCorner size_t nodesRightBackTopCorner(); // alias, see above: nodesBackTopRightCorner size_t nodesRightTopBackCorner(); // alias, see above: nodesBackTopRightCorner // periodicity MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent) size_t nodesOrigin(); // front-bottom-left node, used as reference for periodicity MatS dofs(); // DOF-numbers for each component of each node (sequential) MatS dofsPeriodic(); // ,, for the case that the periodicity if fully eliminated }; // ====================== MESH WITH A FINE LAYER THAT EXPONENTIALLY COARSENS ======================= class FineLayer { private: double m_h; // elementary element edge-size (in all directions) double m_Lx, m_Lz; // mesh size in "x" and "z" ColS m_nelx, m_nelz; // number of elements in "x" and "z" (per el.layer in "y") ColS m_nnd; // total number of nodes in the main node layer (per nd.layer in "y") ColS m_nhx, m_nhy, m_nhz; // element size in each direction (per el.layer in "y") - ColI m_refine; // refine direction (-1: no refine, 0: "x", ...) (per el.layer in "y") + ColI m_refine; // refine direction (-1:no refine, 0:"x", 2:"z") (per el.layer in "y") ColS m_startElem; // start element (per el.layer in "y") ColS m_startNode; // start node (per nd.layer in "y") size_t m_nelem; // number of elements size_t m_nnode; // number of nodes size_t m_nne=8; // number of nodes-per-element size_t m_ndim=3; // number of dimensions public: - // mesh with "nelx*nely*nelz" 'pixels' of edge size "h"; elements are coarsened in "y"-direction + // mesh with "nelx*nely*nelz" elements of edge size "h"; elements are coarsened in "y"-direction FineLayer(size_t nelx, size_t nely, size_t nelz, double h=1., size_t nfine=1); // sizes size_t nelem(); // number of elements size_t nnode(); // number of nodes size_t nne(); // number of nodes-per-element size_t ndim(); // number of dimensions size_t shape(size_t i); // actual shape in a certain direction // mesh MatD coor(); // nodal positions [nnode ,ndim] MatS conn(); // connectivity [nelem ,nne ] + // element sets + ColS elementsMiddleLayer(); // elements in the middle, fine, layer // boundary nodes: planes ColS nodesFront(); // node-numbers along the front plane ColS nodesBack(); // node-numbers along the back plane ColS nodesLeft(); // node-numbers along the left plane ColS nodesRight(); // node-numbers along the right plane ColS nodesBottom(); // node-numbers along the bottom plane ColS nodesTop(); // node-numbers along the top plane // boundary nodes: faces ColS nodesFrontFace(); // node-numbers along the front face ColS nodesBackFace(); // node-numbers along the back face ColS nodesLeftFace(); // node-numbers along the left face ColS nodesRightFace(); // node-numbers along the right face ColS nodesBottomFace(); // node-numbers along the bottom face ColS nodesTopFace(); // node-numbers along the top face // boundary nodes: edges ColS nodesFrontBottomEdge(); // node-numbers along the front - bottom edge ColS nodesFrontTopEdge(); // node-numbers along the front - top edge ColS nodesFrontLeftEdge(); // node-numbers along the front - left edge ColS nodesFrontRightEdge(); // node-numbers along the front - right edge ColS nodesBackBottomEdge(); // node-numbers along the back - bottom edge ColS nodesBackTopEdge(); // node-numbers along the back - top edge ColS nodesBackLeftEdge(); // node-numbers along the back - left edge ColS nodesBackRightEdge(); // node-numbers along the back - right edge ColS nodesBottomLeftEdge(); // node-numbers along the bottom - left edge ColS nodesBottomRightEdge(); // node-numbers along the bottom - right edge ColS nodesTopLeftEdge(); // node-numbers along the top - left edge ColS nodesTopRightEdge(); // node-numbers along the top - right edge // boundary nodes: edges (aliases) ColS nodesBottomFrontEdge(); // alias, see above: nodesFrontBottomEdge ColS nodesBottomBackEdge(); // alias, see above: nodesBackBottomEdge ColS nodesTopFrontEdge(); // alias, see above: nodesFrontTopEdge ColS nodesTopBackEdge(); // alias, see above: nodesBackTopEdge ColS nodesLeftBottomEdge(); // alias, see above: nodesBottomLeftEdge ColS nodesLeftFrontEdge(); // alias, see above: nodesFrontLeftEdge ColS nodesLeftBackEdge(); // alias, see above: nodesBackLeftEdge ColS nodesLeftTopEdge(); // alias, see above: nodesTopLeftEdge ColS nodesRightBottomEdge(); // alias, see above: nodesBottomRightEdge ColS nodesRightTopEdge(); // alias, see above: nodesTopRightEdge ColS nodesRightFrontEdge(); // alias, see above: nodesFrontRightEdge ColS nodesRightBackEdge(); // alias, see above: nodesBackRightEdge // boundary nodes: edges, without corners ColS nodesFrontBottomOpenEdge(); // node-numbers along the front - bottom edge ColS nodesFrontTopOpenEdge(); // node-numbers along the front - top edge ColS nodesFrontLeftOpenEdge(); // node-numbers along the front - left edge ColS nodesFrontRightOpenEdge(); // node-numbers along the front - right edge ColS nodesBackBottomOpenEdge(); // node-numbers along the back - bottom edge ColS nodesBackTopOpenEdge(); // node-numbers along the back - top edge ColS nodesBackLeftOpenEdge(); // node-numbers along the back - left edge ColS nodesBackRightOpenEdge(); // node-numbers along the back - right edge ColS nodesBottomLeftOpenEdge(); // node-numbers along the bottom - left edge ColS nodesBottomRightOpenEdge(); // node-numbers along the bottom - right edge ColS nodesTopLeftOpenEdge(); // node-numbers along the top - left edge ColS nodesTopRightOpenEdge(); // node-numbers along the top - right edge // boundary nodes: edges, without corners (aliases) ColS nodesBottomFrontOpenEdge(); // alias, see above: nodesFrontBottomOpenEdge ColS nodesBottomBackOpenEdge(); // alias, see above: nodesBackBottomOpenEdge ColS nodesTopFrontOpenEdge(); // alias, see above: nodesFrontTopOpenEdge ColS nodesTopBackOpenEdge(); // alias, see above: nodesBackTopOpenEdge ColS nodesLeftBottomOpenEdge(); // alias, see above: nodesBottomLeftOpenEdge ColS nodesLeftFrontOpenEdge(); // alias, see above: nodesFrontLeftOpenEdge ColS nodesLeftBackOpenEdge(); // alias, see above: nodesBackLeftOpenEdge ColS nodesLeftTopOpenEdge(); // alias, see above: nodesTopLeftOpenEdge ColS nodesRightBottomOpenEdge(); // alias, see above: nodesBottomRightOpenEdge ColS nodesRightTopOpenEdge(); // alias, see above: nodesTopRightOpenEdge ColS nodesRightFrontOpenEdge(); // alias, see above: nodesFrontRightOpenEdge ColS nodesRightBackOpenEdge(); // alias, see above: nodesBackRightOpenEdge // boundary nodes: corners size_t nodesFrontBottomLeftCorner(); // node-number of the front - bottom - left corner size_t nodesFrontBottomRightCorner(); // node-number of the front - bottom - right corner size_t nodesFrontTopLeftCorner(); // node-number of the front - top - left corner size_t nodesFrontTopRightCorner(); // node-number of the front - top - right corner size_t nodesBackBottomLeftCorner(); // node-number of the back - bottom - left corner size_t nodesBackBottomRightCorner(); // node-number of the back - bottom - right corner size_t nodesBackTopLeftCorner(); // node-number of the back - top - left corner size_t nodesBackTopRightCorner(); // node-number of the back - top - right corner // boundary nodes: corners (aliases) size_t nodesFrontLeftBottomCorner(); // alias, see above: nodesFrontBottomLeftCorner size_t nodesBottomFrontLeftCorner(); // alias, see above: nodesFrontBottomLeftCorner size_t nodesBottomLeftFrontCorner(); // alias, see above: nodesFrontBottomLeftCorner size_t nodesLeftFrontBottomCorner(); // alias, see above: nodesFrontBottomLeftCorner size_t nodesLeftBottomFrontCorner(); // alias, see above: nodesFrontBottomLeftCorner size_t nodesFrontRightBottomCorner(); // alias, see above: nodesFrontBottomRightCorner size_t nodesBottomFrontRightCorner(); // alias, see above: nodesFrontBottomRightCorner size_t nodesBottomRightFrontCorner(); // alias, see above: nodesFrontBottomRightCorner size_t nodesRightFrontBottomCorner(); // alias, see above: nodesFrontBottomRightCorner size_t nodesRightBottomFrontCorner(); // alias, see above: nodesFrontBottomRightCorner size_t nodesFrontLeftTopCorner(); // alias, see above: nodesFrontTopLeftCorner size_t nodesTopFrontLeftCorner(); // alias, see above: nodesFrontTopLeftCorner size_t nodesTopLeftFrontCorner(); // alias, see above: nodesFrontTopLeftCorner size_t nodesLeftFrontTopCorner(); // alias, see above: nodesFrontTopLeftCorner size_t nodesLeftTopFrontCorner(); // alias, see above: nodesFrontTopLeftCorner size_t nodesFrontRightTopCorner(); // alias, see above: nodesFrontTopRightCorner size_t nodesTopFrontRightCorner(); // alias, see above: nodesFrontTopRightCorner size_t nodesTopRightFrontCorner(); // alias, see above: nodesFrontTopRightCorner size_t nodesRightFrontTopCorner(); // alias, see above: nodesFrontTopRightCorner size_t nodesRightTopFrontCorner(); // alias, see above: nodesFrontTopRightCorner size_t nodesBackLeftBottomCorner(); // alias, see above: nodesBackBottomLeftCorner size_t nodesBottomBackLeftCorner(); // alias, see above: nodesBackBottomLeftCorner size_t nodesBottomLeftBackCorner(); // alias, see above: nodesBackBottomLeftCorner size_t nodesLeftBackBottomCorner(); // alias, see above: nodesBackBottomLeftCorner size_t nodesLeftBottomBackCorner(); // alias, see above: nodesBackBottomLeftCorner size_t nodesBackRightBottomCorner(); // alias, see above: nodesBackBottomRightCorner size_t nodesBottomBackRightCorner(); // alias, see above: nodesBackBottomRightCorner size_t nodesBottomRightBackCorner(); // alias, see above: nodesBackBottomRightCorner size_t nodesRightBackBottomCorner(); // alias, see above: nodesBackBottomRightCorner size_t nodesRightBottomBackCorner(); // alias, see above: nodesBackBottomRightCorner size_t nodesBackLeftTopCorner(); // alias, see above: nodesBackTopLeftCorner size_t nodesTopBackLeftCorner(); // alias, see above: nodesBackTopLeftCorner size_t nodesTopLeftBackCorner(); // alias, see above: nodesBackTopLeftCorner size_t nodesLeftBackTopCorner(); // alias, see above: nodesBackTopLeftCorner size_t nodesLeftTopBackCorner(); // alias, see above: nodesBackTopLeftCorner size_t nodesBackRightTopCorner(); // alias, see above: nodesBackTopRightCorner size_t nodesTopBackRightCorner(); // alias, see above: nodesBackTopRightCorner size_t nodesTopRightBackCorner(); // alias, see above: nodesBackTopRightCorner size_t nodesRightBackTopCorner(); // alias, see above: nodesBackTopRightCorner size_t nodesRightTopBackCorner(); // alias, see above: nodesBackTopRightCorner // periodicity MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent) size_t nodesOrigin(); // front-bottom-left node, used as reference for periodicity MatS dofs(); // DOF-numbers for each component of each node (sequential) MatS dofsPeriodic(); // ,, for the case that the periodicity if fully eliminated }; // ================================================================================================= }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/MeshQuad4.cpp b/src/GooseFEM/MeshQuad4.cpp index 9e23c3b..58e52e5 100644 --- a/src/GooseFEM/MeshQuad4.cpp +++ b/src/GooseFEM/MeshQuad4.cpp @@ -1,876 +1,922 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESHQUAD4_CPP #define GOOSEFEM_MESHQUAD4_CPP // ------------------------------------------------------------------------------------------------- #include "MeshQuad4.h" // ===================================== GooseFEM::Mesh::Quad4 ===================================== namespace GooseFEM { namespace Mesh { namespace Quad4 { // ===================================== CLASS - REGULAR MESH ====================================== // ------------------------------------------ constructor ------------------------------------------ -inline Regular::Regular(size_t nx, size_t ny, double h): -m_nx(nx), m_ny(ny), m_h(h) +inline Regular::Regular(size_t nelx, size_t nely, double h): +m_h(h), m_nelx(nelx), m_nely(nely) { - assert( m_nx >= 1 ); - assert( m_ny >= 1 ); + assert( m_nelx >= 1 ); + assert( m_nely >= 1 ); - m_nnode = (m_nx+1) * (m_ny+1); - m_nelem = m_nx * m_ny ; + m_nnode = (m_nelx+1) * (m_nely+1); + m_nelem = m_nelx * m_nely ; } // -------------------------------------- number of elements --------------------------------------- inline size_t Regular::nelem() { return m_nelem; } // ---------------------------------------- number of nodes ---------------------------------------- inline size_t Regular::nnode() { return m_nnode; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t Regular::nne() { return m_nne; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t Regular::ndim() { return m_ndim; } // --------------------------------- coordinates (nodal positions) --------------------------------- inline MatD Regular::coor() { - MatD coor(m_nnode,m_ndim); + MatD out(m_nnode, m_ndim); - ColD x = ColD::LinSpaced(m_nx+1, 0.0, m_h*static_cast(m_nx)); - ColD y = ColD::LinSpaced(m_ny+1, 0.0, m_h*static_cast(m_ny)); + ColD x = ColD::LinSpaced(m_nelx+1, 0.0, m_h*static_cast(m_nelx)); + ColD y = ColD::LinSpaced(m_nely+1, 0.0, m_h*static_cast(m_nely)); size_t inode = 0; - for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy ) { - for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix ) { - coor(inode,0) = x(ix); - coor(inode,1) = y(iy); + for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) { + for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) { + out(inode,0) = x(ix); + out(inode,1) = y(iy); ++inode; } } - return coor; + return out; } // ---------------------------- connectivity (node-numbers per element) ---------------------------- inline MatS Regular::conn() { - MatS conn(m_nelem,m_nne); + MatS out(m_nelem,m_nne); size_t ielem = 0; - for ( size_t iy = 0 ; iy < m_ny ; ++iy ) { - for ( size_t ix = 0 ; ix < m_nx ; ++ix ) { - conn(ielem,0) = (iy+0)*(m_nx+1) + (ix+0); - conn(ielem,1) = (iy+0)*(m_nx+1) + (ix+1); - conn(ielem,3) = (iy+1)*(m_nx+1) + (ix+0); - conn(ielem,2) = (iy+1)*(m_nx+1) + (ix+1); + for ( size_t iy = 0 ; iy < m_nely ; ++iy ) { + for ( size_t ix = 0 ; ix < m_nelx ; ++ix ) { + out(ielem,0) = (iy )*(m_nelx+1) + (ix ); + out(ielem,1) = (iy )*(m_nelx+1) + (ix+1); + out(ielem,3) = (iy+1)*(m_nelx+1) + (ix ); + out(ielem,2) = (iy+1)*(m_nelx+1) + (ix+1); ++ielem; } } - return conn; + return out; } // ------------------------------ node-numbers along the bottom edge ------------------------------- inline ColS Regular::nodesBottomEdge() { - ColS nodes(m_nx+1); + ColS out(m_nelx+1); - for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix ) - nodes(ix) = ix; + for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) + out(ix) = ix; - return nodes; + return out; } // -------------------------------- node-numbers along the top edge -------------------------------- inline ColS Regular::nodesTopEdge() { - ColS nodes(m_nx+1); + ColS out(m_nelx+1); - for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix ) - nodes(ix) = ix + m_ny*(m_nx+1); + for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) + out(ix) = ix + m_nely*(m_nelx+1); - return nodes; + return out; } -// ----------------------------- node-numbers along the node left edge ----------------------------- +// ------------------------------- node-numbers along the left edge -------------------------------- inline ColS Regular::nodesLeftEdge() { - ColS nodes(m_ny+1); + ColS out(m_nely+1); - for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy ) - nodes(iy) = iy*(m_nx+1); + for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) + out(iy) = iy*(m_nelx+1); - return nodes; + return out; } // ------------------------------- node-numbers along the right edge ------------------------------- inline ColS Regular::nodesRightEdge() { - ColS nodes(m_ny+1); + ColS out(m_nely+1); + + for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) + out(iy) = iy*(m_nelx+1) + m_nelx; + + return out; +} + +// ---------------------- node-numbers along the bottom edge, without corners ---------------------- + +inline ColS Regular::nodesBottomOpenEdge() +{ + ColS out(m_nelx-1); + + for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) + out(ix-1) = ix; + + return out; +} + +// ----------------------- node-numbers along the top edge, without corners ------------------------ + +inline ColS Regular::nodesTopOpenEdge() +{ + ColS out(m_nelx-1); + + for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) + out(ix-1) = ix + m_nely*(m_nelx+1); + + return out; +} + +// ----------------------- node-numbers along the left edge, without corners ----------------------- + +inline ColS Regular::nodesLeftOpenEdge() +{ + ColS out(m_nely-1); - for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy ) - nodes(iy) = iy*(m_nx+1) + m_nx; + for ( size_t iy = 1 ; iy < m_nely ; ++iy ) + out(iy-1) = iy*(m_nelx+1); - return nodes; + return out; +} + +// ---------------------- node-numbers along the right edge, without corners ----------------------- + +inline ColS Regular::nodesRightOpenEdge() +{ + ColS out(m_nely-1); + + for ( size_t iy = 1 ; iy < m_nely ; ++iy ) + out(iy-1) = iy*(m_nelx+1) + m_nelx; + + return out; } // ----------------------------- node-number of the bottom-left corner ----------------------------- inline size_t Regular::nodesBottomLeftCorner() { return 0; } // ---------------------------- node-number of the bottom-right corner ----------------------------- inline size_t Regular::nodesBottomRightCorner() { - return m_nx; + return m_nelx; } // ------------------------------ node-number of the top-left corner ------------------------------- inline size_t Regular::nodesTopLeftCorner() { - return m_ny*(m_nx+1); + return m_nely*(m_nelx+1); } // ------------------------------ node-number of the top-right corner ------------------------------ inline size_t Regular::nodesTopRightCorner() { - return (m_ny+1)*(m_nx+1) - 1; + return m_nely*(m_nelx+1) + m_nelx; } // ----------------------------- node-number of the corners (aliases) ------------------------------ inline size_t Regular::nodesLeftBottomCorner() { return nodesBottomLeftCorner(); } inline size_t Regular::nodesLeftTopCorner() { return nodesTopLeftCorner(); } inline size_t Regular::nodesRightBottomCorner() { return nodesBottomRightCorner(); } inline size_t Regular::nodesRightTopCorner() { return nodesTopRightCorner(); } // ------------------------------ node-numbers of periodic node-pairs ------------------------------ inline MatS Regular::nodesPeriodic() { - // edges - ColS bot = nodesBottomEdge(); - ColS top = nodesTopEdge(); - ColS rgt = nodesRightEdge(); - ColS lft = nodesLeftEdge(); + // edges (without corners) + ColS bot = nodesBottomOpenEdge(); + ColS top = nodesTopOpenEdge(); + ColS lft = nodesLeftOpenEdge(); + ColS rgt = nodesRightOpenEdge(); // allocate nodal ties - MatS nodes(bot.size()-2 + lft.size()-2 + 3, 2); + // - number of tying per category + size_t tedge = bot.size() + lft.size(); + size_t tnode = 3; + // - allocate + MatS out(tedge+tnode, 2); // counter size_t i = 0; // tie all corners to one corner - nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesBottomRightCorner(); ++i; - nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopRightCorner(); ++i; - nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopLeftCorner(); ++i; + out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesBottomRightCorner(); ++i; + out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopRightCorner(); ++i; + out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopLeftCorner(); ++i; - // tie edges to each other (exclude corners) - for ( auto j = 1 ; j < bot.size()-1 ; ++j ) { nodes(i,0) = bot(j); nodes(i,1) = top(j); ++i; } - for ( auto j = 1 ; j < lft.size()-1 ; ++j ) { nodes(i,0) = lft(j); nodes(i,1) = rgt(j); ++i; } + // tie all corresponding edges to each other + for ( auto j = 0 ; j(nodePer.rows()); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs for ( size_t i = 0 ; i < nper ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) out(nodePer(i,1),j) = out(nodePer(i,0),j); // renumber "out" to be sequential return GooseFEM::Mesh::renumber(out); } // ==================================== CLASS - FINELAYER MESH ===================================== // ------------------------------------------ constructor ------------------------------------------ -inline FineLayer::FineLayer(size_t nx, size_t ny, double h, size_t nfine, size_t nskip): -m_h(h), m_nx(nx) +inline FineLayer::FineLayer(size_t nelx, size_t nely, double h, size_t nfine): +m_h(h), m_nelx(nelx) { - assert( nx >= 1 ); - assert( ny >= 1 ); + // basic assumptions + assert( nelx >= 1 ); + assert( nely >= 1 ); + + // store basic info + m_Lx = m_h * static_cast(nelx); - // local counters - size_t N; + // compute element size in y-direction (use symmetry, compute upper half) + // ---------------------------------------------------------------------- - // compute the element size : based on symmetric half - // -------------------------------------------------- + // temporary variables + size_t nmin, ntot; + ColS nhx(nely), nhy(nely); + ColI refine(nely); - // convert height to half of the height - if ( ny % 2 == 0 ) ny = ny /2; - else ny = (ny+1)/2; + // minimum height in y-direction (half of the height because of symmetry) + if ( nely % 2 == 0 ) nmin = nely /2; + else nmin = (nely +1)/2; - // convert to half the number of fine layer (minimum 1) + // minimum number of fine layers in y-direction (minimum 1, middle layer part of this half) if ( nfine % 2 == 0 ) nfine = nfine /2 + 1; else nfine = (nfine+1)/2; if ( nfine < 1 ) nfine = 1; + if ( nfine > nmin ) nfine = nmin; + + // initialize to state with only fine elements + nhx .setOnes(); + nhy .setOnes(); + refine.setConstant(-1); + + // loop over element layers in y-direction, try to coarsen using these rules: + // (1) element size in y-direction <= distance to origin in y-direction + // (2) element size in x-direction should fit the total number of elements in x-direction + // (3) a certain number of layers have the minimum size "1" (are fine) + for ( size_t iy = nfine ; ; ) + { + // initialize current size in y-direction + if ( iy == nfine ) ntot = nfine; + // check to stop + if ( iy >= nely or ntot >= nmin ) { nely = iy; break; } - // check the number of fine layers from the center - assert( nfine <= ny ); - - // define arrays to determine to coarsen - ColS n (ny+1); // size of the element (number of times "h"): even sizes -> coarsen - ColS nsum(ny+1); // current 'size' of the mesh in y-direction - - // initialize all elements of size "1" - n.setOnes(); + // rules (1,2) satisfied: coarsen in x-direction + if ( 3*nhy(iy) <= ntot and nelx%(3*nhx(iy)) == 0 ) + { + refine (iy ) = 0; + nhy (iy ) *= 2; + nhy.segment(iy+1,nely-iy-1) *= 3; + nhx.segment(iy ,nely-iy ) *= 3; + } - // size of the mesh in y-direction after "i" elements - // - initialize - nsum(0) = 1; - // - compute based on the current value of "n" - for ( size_t i = 1 ; i < ny+1 ; ++i ) - nsum(i) = nsum(i-1) + n(i-1); - - // loop in vertical direction and check to coarsen; rules: - // * the size of the element cannot be greater than the distance - // * the size of the element should match the horizontal direction - // * skip a certain number of layers that have the minimum size "1" (are fine) - // - initialize next element size - size_t next = 1; - // - loop to coarsen - for ( size_t i = nfine ; i < ny+1 ; ++i ) - { - // -- set element size, based on whether to coarsen or not - if ( 3*n(i-1) <= nsum(i-1) and m_nx % ( 3*next ) == 0 ) { n(i) = next * 2; next *= 3; } - else { n(i) = next; } - // -- compute the total number of elements - nsum(i) = nsum(i-1) + n(i); + // update the number of elements in y-direction + ntot += nhy(iy); + // proceed to next element layer in y-direction + ++iy; + // check to stop + if ( iy >= nely or ntot >= nmin ) { nely = iy; break; } } - // truncate to size "N" such that the height does not exceed "ny" - for ( N = 0 ; N < ny+1 ; ++N ) - if ( nsum(N) > ny+1 ) - break; - - // fiddle a little bit to approximate the size as best as possible - if ( N > 1 and N < ny ) + // symmetrize, compute full information + // ------------------------------------ + + // allocate mesh constructor parameters + m_nhx .conservativeResize(nely*2-1); + m_nhy .conservativeResize(nely*2-1); + m_refine .conservativeResize(nely*2-1); + m_nelx .conservativeResize(nely*2-1); + m_nnd .conservativeResize(nely*2 ); + m_startElem.conservativeResize(nely*2-1); + m_startNode.conservativeResize(nely*2 ); + + // fill + // - lower half + for ( size_t iy = 0 ; iy < nely ; ++iy ) { - // - set all sequential steps to last value, make sure not to use a coarsening step - if ( n(N-1) % 2 == 0 ) for ( size_t i = N ; i < ny+1 ; ++i ) n(i) = n(i-1)/2*3; - else for ( size_t i = N ; i < ny+1 ; ++i ) n(i) = n(i-1); - // - fix total size - for ( size_t i = N ; i < ny+1 ; ++i ) nsum(i) = nsum(i-1) + n(i); - // - truncate such that the height does not exceed "ny" - for ( N = 0 ; N < ny ; ++N ) - if ( nsum(N) > ny+1 ) - break; - // - check to extend one - if ( nsum(N) - ny < ny - nsum(N-1) ) - ++N; + m_nhx (iy) = nhx (nely-iy-1); + m_nhy (iy) = nhy (nely-iy-1); + m_refine(iy) = refine(nely-iy-1); } - - // skip the last "nskip" coarsening steps - if ( nskip > 0 ) + // - upper half + for ( size_t iy = 0 ; iy < nely-1 ; ++iy ) { - // - counter : number of steps skipped - size_t iskip = 0; - // - loop to skip - for ( size_t i = N ; i-- > 1 ; ) - { - // -- check -> quit if the number of steps is reached - if ( iskip >= nskip ) break; - // -- check -> quit if the finest layer has been reached - if ( n(i) == 1 ) break; - // -- if coarsening step found -> remove - if ( n(i) % 2 == 0 ) - { - size_t k = n(i) / 2; - for ( size_t j = i ; j < ny+1 ; ++j ) n(j) = k; - iskip++; - } - } - // - update the height after "i" elements - for ( size_t i = 1 ; i < ny+1 ; ++i ) - nsum(i) = nsum(i-1) + n(i); - // - truncate such that the height does not exceed "ny" - for ( N = 0 ; N < ny+1 ; ++N ) - if ( nsum(N) > ny+1 ) - break; + m_nhx (iy+nely) = nhx (iy+1); + m_nhy (iy+nely) = nhy (iy+1); + m_refine(iy+nely) = refine(iy+1); } - // truncate - n.conservativeResize(N); + // update size + nely = m_nhx.size(); - // compute mesh dimensions : based on the entire mesh - // -------------------------------------------------- + // compute the number of elements per element layer in y-direction + for ( size_t iy = 0 ; iy < nely ; ++iy ) + m_nelx(iy) = nelx / m_nhx(iy); - // symmetrize - // - get size - N = static_cast(n.size()); - // - allocate - m_nh.conservativeResize(2*N-1); - // - store symmetrized - for ( size_t i = 0 ; i < N ; i++ ) m_nh( i ) = n(N-1-i); - for ( size_t i = 1 ; i < N ; i++ ) m_nh(N+i-1) = n( i); + // compute the number of nodes per node layer in y-direction + for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) m_nnd(iy ) = m_nelx(iy)+1; + for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) m_nnd(iy+1) = m_nelx(iy)+1; - // allocate counters - m_startElem.conservativeResize(m_nh.size() ); - m_startNode.conservativeResize(m_nh.size()+1); + // compute mesh dimensions + // ----------------------- - // compute the dimensions of the mesh - // - initialize + // initialize m_nnode = 0; m_nelem = 0; - N = static_cast(m_nh.size()); m_startNode(0) = 0; - // loop from bottom to middle : elements become finer - for ( size_t i = 0 ; i < (N-1)/2 ; ++i ) + // loop over element layers (bottom -> middle, elements become finer) + for ( size_t i = 0 ; i < (nely-1)/2 ; ++i ) { - // - store the first element of the row + // - store the first element of the layer m_startElem(i) = m_nelem; - // - add the nodes of this row - if ( m_nh(i)%2 == 0 ) { m_nnode += 3 * m_nx/(3*m_nh(i)/2) + 1; } - else { m_nnode += m_nx/ m_nh(i) + 1; } - // - add the elements of this row - if ( m_nh(i)%2 == 0 ) { m_nelem += 4 * m_nx/(3*m_nh(i)/2); } - else { m_nelem += m_nx/ m_nh(i) ; } - // - store the starting node of the next row + // - add the nodes of this layer + if ( m_refine(i) == 0 ) { m_nnode += (3*m_nelx(i)+1); } + else { m_nnode += ( m_nelx(i)+1); } + // - add the elements of this layer + if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ); } + else { m_nelem += ( m_nelx(i) ); } + // - store the starting node of the next layer m_startNode(i+1) = m_nnode; } - // loop from middle to top : elements become coarser - for ( size_t i = (N-1)/2 ; i < N ; ++i ) + // loop over element layers (middle -> top, elements become coarser) + for ( size_t i = (nely-1)/2 ; i < nely ; ++i ) { - // - store the first element of the row + // - store the first element of the layer m_startElem(i) = m_nelem; - // - add the nodes of this row - if ( m_nh(i)%2 == 0 ) { m_nnode += 1 + m_nx/(m_nh(i)/2) + 2 * m_nx/(3*m_nh(i)/2); } - else { m_nnode += 1 + m_nx/ m_nh(i); } - // - add the elements of this row - if ( m_nh(i)%2 == 0 ) { m_nelem += 4 * m_nx/(3*m_nh(i)/2); } - else { m_nelem += m_nx/ m_nh(i) ; } - // - store the starting node of the next row + // - add the nodes of this layer + if ( m_refine(i) == 0 ) { m_nnode += (5*m_nelx(i)+1); } + else { m_nnode += ( m_nelx(i)+1); } + // - add the elements of this layer + if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ); } + else { m_nelem += ( m_nelx(i) ); } + // - store the starting node of the next layer m_startNode(i+1) = m_nnode; } - - // add the top row of nodes to the number of rows (starting node already stored) - if ( m_nh(N-1)%2 == 0 ) { m_nnode += 1 + m_nx/(3*m_nh(N-1)/2); } - else { m_nnode += 1 + m_nx/ m_nh(N-1); } + // - add the top row of nodes + m_nnode += m_nelx(nely-1)+1; } // -------------------------------------- number of elements --------------------------------------- inline size_t FineLayer::nelem() { return m_nelem; } // ---------------------------------------- number of nodes ---------------------------------------- inline size_t FineLayer::nnode() { return m_nnode; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t FineLayer::nne() { return m_nne; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t FineLayer::ndim() { return m_ndim; } // ---------------------------- actual number of nodes in one direction ---------------------------- inline size_t FineLayer::shape(size_t i) { - // check the index - assert( i >= 0 and i < 2 ); - - // horizontal direction - if ( i == 0 ) - return m_nx; + assert( i >= 0 and i <= 1 ); - // vertical direction - // - zero-initialize - size_t n = 0; - // - cumulative sum of the size per row - for ( size_t i = 0 ; i < static_cast(m_nh.size()) ; ++i ) - n += m_nh(i); + if ( i == 0 ) return m_nelx.maxCoeff(); + else return m_nhy .sum(); - return n; } // --------------------------------- coordinates (nodal positions) --------------------------------- inline MatD FineLayer::coor() { // allocate output MatD out(m_nnode, m_ndim); - // initialize position in horizontal direction (of the finest nodes) - ColD x = ColD::LinSpaced(m_nx+1, 0.0, m_h*static_cast(m_nx)); - - // zero-initialize current node and height: loop from bottom to top; number of rows - double h = 0; + // current node, number of element layers size_t inode = 0; - size_t N = static_cast(m_nh.size()); + size_t nely = static_cast(m_nhy.size()); + + // y-position of each main node layer (i.e. excluding node layers for refinement/coarsening) + // - allocate + ColD y(nely+1); + // - initialize + y(0) = 0.0; + // - compute + for ( size_t iy = 1 ; iy < nely+1 ; ++iy ) + y(iy) = y(iy-1) + m_nhy(iy-1) * m_h; + + // loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of nodes + // -------------------------------------------------------------------------------------------- - // lower half - for ( size_t i = 0 ; i < (N-1)/2 ; ++i ) + for ( size_t iy = 0 ; ; ++iy ) { - // - refinement row - if ( m_nh(i)%2 == 0 ) - { - // -- bottom row: coarse nodes - for ( size_t j = 0 ; j < m_nx+1 ; j += 3*m_nh(i)/2 ) - { - out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; - } - h += static_cast(m_nh(i)/2); - // -- middle row: part the fine nodes - for ( size_t j = 0 ; j < m_nx ; j += 3*m_nh(i)/2 ) - { - out(inode,0) = x(j+m_nh(i)/2); out(inode,1) = h*m_h; inode++; - out(inode,0) = x(j+m_nh(i) ); out(inode,1) = h*m_h; inode++; - } - h += static_cast(m_nh(i)/2); + // get positions along the x- and z-axis + ColD x = ColD::LinSpaced(m_nelx(iy)+1, 0.0, m_Lx); + + // add nodes of the bottom layer of this element + for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { + out(inode,0) = x(ix); + out(inode,1) = y(iy); + ++inode; } - // - normal row - else + + // stop at middle layer + if ( iy == (nely-1)/2 ) + break; + + // add extra nodes of the intermediate layer, for refinement in x-direction + if ( m_refine(iy) == 0 ) { - for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(i) ) - { - out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; + // - get position offset in x- and y-direction + double dx = m_h * static_cast(m_nhx(iy)/3); + double dy = m_h * static_cast(m_nhy(iy)/2); + // - add nodes of the intermediate layer + for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { + for ( size_t j = 0 ; j < 2 ; ++j ) { + out(inode,0) = x(ix) + dx * static_cast(j+1); + out(inode,1) = y(iy) + dy; + ++inode; + } } - h += static_cast(m_nh(i)); } } - // top half - for ( size_t i = (N-1)/2 ; i < N ; ++i ) + // loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes + // -------------------------------------------------------------------------------------- + + for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { - // - coarsening row - if ( m_nh(i)%2 == 0 ) - { - // -- bottom row: fine nodes - for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(i)/2 ) - { - out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; - } - h += static_cast(m_nh(i)/2); - // -- middle row: part the fine nodes - for ( size_t j = 0 ; j < m_nx ; j += 3*m_nh(i)/2 ) - { - out(inode,0) = x(j+m_nh(i)/2); out(inode,1) = h*m_h; inode++; - out(inode,0) = x(j+m_nh(i) ); out(inode,1) = h*m_h; inode++; - } - h += static_cast(m_nh(i)/2); - } - // - normal row - else + // get positions along the x- and z-axis + ColD x = ColD::LinSpaced(m_nelx(iy)+1, 0.0, m_Lx); + + // add extra nodes of the intermediate layer, for refinement in x-direction + if ( m_refine(iy) == 0 ) { - for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(i) ) - { - out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; + // - get position offset in x- and y-direction + double dx = m_h * static_cast(m_nhx(iy)/3); + double dy = m_h * static_cast(m_nhy(iy)/2); + // - add nodes of the intermediate layer + for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { + for ( size_t j = 0 ; j < 2 ; ++j ) { + out(inode,0) = x(ix) + dx * static_cast(j+1); + out(inode,1) = y(iy) + dy; + ++inode; + } } - h += static_cast(m_nh(i)); } - } - // top row - if ( m_nh(N-1)%2 == 0 ) - { - for ( size_t j = 0 ; j < m_nx+1 ; j += 3*m_nh(N-1)/2 ) - { - out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; - } - } - else - { - for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(N-1) ) - { - out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; + // add nodes of the top layer of this element + for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { + out(inode,0) = x(ix ); + out(inode,1) = y(iy+1); + ++inode; } } return out; } // ---------------------------- connectivity (node-numbers per element) ---------------------------- inline MatS FineLayer::conn() { // allocate output MatS out(m_nelem, m_nne); - // current element, number of rows, beginning nodes + // current element, number of element layers, starting nodes of each node layer size_t ielem = 0; - size_t N = static_cast(m_nh.size()); + size_t nely = static_cast(m_nhy.size()); size_t bot,mid,top; - // bottom half - for ( size_t i = 0 ; i < (N-1)/2 ; ++i ) + // loop over all element layers + for ( size_t iy = 0 ; iy < nely ; ++iy ) { - // - starting nodes of bottom and top layer (and middle layer, only relevant for even shape) - if ( true ) bot = m_startNode(i ); - if ( m_nh(i)%2 == 0 ) mid = m_startNode(i )+m_nx/(3*m_nh(i)/2)+1; - if ( true ) top = m_startNode(i+1); + // - get: starting nodes of bottom(, middle) and top layer + bot = m_startNode(iy ); + mid = m_startNode(iy ) + m_nnd(iy); + top = m_startNode(iy+1); - // - even sized shape: refinement layer - if ( m_nh(i)%2 == 0 ) + // - define connectivity: no coarsening/refinement + if ( m_refine(iy) == -1 ) { - for ( size_t j = 0 ; j < m_nx/(3*m_nh(i)/2) ; ++j ) - { - // -- lower - out(ielem,0) = j +bot; - out(ielem,1) = j +1+bot; - out(ielem,2) = j*2+1+mid; - out(ielem,3) = j*2 +mid; - ielem++; - // -- upper right - out(ielem,0) = j +1+bot; - out(ielem,1) = j*3+3+top; - out(ielem,2) = j*3+2+top; - out(ielem,3) = j*2+1+mid; - ielem++; - // -- upper middle - out(ielem,0) = j*2 +mid; - out(ielem,1) = j*2+1+mid; - out(ielem,2) = j*3+2+top; - out(ielem,3) = j*3+1+top; - ielem++; - // -- upper left - out(ielem,0) = j +bot; - out(ielem,1) = j*2 +mid; - out(ielem,2) = j*3+1+top; - out(ielem,3) = j*3 +top; + for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { + out(ielem,0) = bot + (ix ); + out(ielem,1) = bot + (ix+1); + out(ielem,2) = top + (ix+1); + out(ielem,3) = top + (ix ); ielem++; } } - // - odd sized shape: normal layer - else - { - for ( size_t j = 0 ; j < m_nx/m_nh(i) ; ++j ) - { - out(ielem,0) = j +bot; - out(ielem,1) = j+1+bot; - out(ielem,2) = j+1+top; - out(ielem,3) = j +top; - ielem++; - } - } - } - // top half - for ( size_t i = (N-1)/2 ; i < N ; ++i ) - { - // - starting nodes of bottom and top layer (and middle layer, only relevant for even shape) - if ( true ) bot = m_startNode(i ); - if ( m_nh(i)%2 == 0 ) mid = m_startNode(i )+m_nx/(m_nh(i)/2)+1; - if ( true ) top = m_startNode(i+1); - - // - even sized shape: refinement layer - if ( m_nh(i)%2 == 0 ) + // - define connectivity: refinement along the x-direction (below the middle layer) + else if ( m_refine(iy) == 0 and iy <= (nely-1)/2 ) { - for ( size_t j = 0 ; j < m_nx/(3*m_nh(i)/2) ; ++j ) - { - // -- lower left - out(ielem,0) = j*3 +bot; - out(ielem,1) = j*3+1+bot; - out(ielem,2) = j*2 +mid; - out(ielem,3) = j +top; + for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { + // -- bottom element + out(ielem,0) = bot + ( ix ); + out(ielem,1) = bot + ( ix+1); + out(ielem,2) = mid + (2*ix+1); + out(ielem,3) = mid + (2*ix ); ielem++; - // -- lower middle - out(ielem,0) = j*3+1+bot; - out(ielem,1) = j*3+2+bot; - out(ielem,2) = j*2+1+mid; - out(ielem,3) = j*2 +mid; + // -- top-right element + out(ielem,0) = bot + ( ix+1); + out(ielem,1) = top + (3*ix+3); + out(ielem,2) = top + (3*ix+2); + out(ielem,3) = mid + (2*ix+1); ielem++; - // -- lower right - out(ielem,0) = j*3+2+bot; - out(ielem,1) = j*3+3+bot; - out(ielem,2) = j +1+top; - out(ielem,3) = j*2+1+mid; + // -- top-center element + out(ielem,0) = mid + (2*ix ); + out(ielem,1) = mid + (2*ix+1); + out(ielem,2) = top + (3*ix+2); + out(ielem,3) = top + (3*ix+1); ielem++; - // -- upper - out(ielem,0) = j +top; - out(ielem,1) = j*2 +mid; - out(ielem,2) = j*2+1+mid; - out(ielem,3) = j +1+top; + // -- top-left element + out(ielem,0) = bot + ( ix ); + out(ielem,1) = mid + (2*ix ); + out(ielem,2) = top + (3*ix+1); + out(ielem,3) = top + (3*ix ); ielem++; } } - // - odd sized shape: normal layer - else + + // - define connectivity: coarsening along the x-direction (above the middle layer) + else if ( m_refine(iy) == 0 and iy > (nely-1)/2 ) { - for ( size_t j = 0 ; j < m_nx/m_nh(i) ; ++j ) - { - out(ielem,0) = j +bot; - out(ielem,1) = j+1+bot; - out(ielem,2) = j+1+top; - out(ielem,3) = j +top; + for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { + // -- lower-left element + out(ielem,0) = bot + (3*ix ); + out(ielem,1) = bot + (3*ix+1); + out(ielem,2) = mid + (2*ix ); + out(ielem,3) = top + ( ix ); + ielem++; + // -- lower-center element + out(ielem,0) = bot + (3*ix+1); + out(ielem,1) = bot + (3*ix+2); + out(ielem,2) = mid + (2*ix+1); + out(ielem,3) = mid + (2*ix ); + ielem++; + // -- lower-right element + out(ielem,0) = bot + (3*ix+2); + out(ielem,1) = bot + (3*ix+3); + out(ielem,2) = top + ( ix+1); + out(ielem,3) = mid + (2*ix+1); + ielem++; + // -- upper element + out(ielem,0) = mid + (2*ix ); + out(ielem,1) = mid + (2*ix+1); + out(ielem,2) = top + ( ix+1); + out(ielem,3) = top + ( ix ); ielem++; } } } return out; } // ------------------------------ element numbers of the middle layer ------------------------------ inline ColS FineLayer::elementsMiddleLayer() { - ColS out(m_nx); + // number of element layers in y-direction, the index of the middle layer + size_t nely = static_cast(m_nhy.size()); + size_t iy = (nely-1)/2; - size_t j = m_startElem((m_startElem.size()-1)/2); + ColS out(m_nelx(iy)); - for ( size_t i = 0 ; i < m_nx ; ++i ) out(i) = i+j; + for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) + out(ix) = m_startElem(iy) + ix; return out; } // ------------------------------ node-numbers along the bottom edge ------------------------------- inline ColS FineLayer::nodesBottomEdge() { - ColS nodes; - - if ( m_nh(0)%2 == 0 ) nodes.conservativeResize(m_nx/(3*m_nh(0)/2)+1); - else nodes.conservativeResize(m_nx/ m_nh(0) +1); + ColS out(m_nelx(0)+1); - for ( size_t i = 0 ; i < static_cast(nodes.size()) ; ++i ) nodes(i) = i; + for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix ) + out(ix) = m_startNode(0) + ix; - return nodes; + return out; } // -------------------------------- node-numbers along the top edge -------------------------------- inline ColS FineLayer::nodesTopEdge() { - ColS nodes; + size_t nely = static_cast(m_nhy.size()); - if ( m_nh(0)%2 == 0 ) nodes.conservativeResize(m_nx/(3*m_nh(0)/2)+1); - else nodes.conservativeResize(m_nx/ m_nh(0) +1); + ColS out(m_nelx(nely-1)+1); - size_t j = m_startNode(m_startNode.size()-1); + for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix ) + out(ix) = m_startNode(nely) + ix; - for ( size_t i = 0 ; i < static_cast(nodes.size()) ; ++i ) nodes(i) = i+j; - - return nodes; + return out; } -// ----------------------------- node-numbers along the node left edge ----------------------------- +// ------------------------------- node-numbers along the left edge -------------------------------- inline ColS FineLayer::nodesLeftEdge() { - return m_startNode; + size_t nely = static_cast(m_nhy.size()); + + ColS out(nely+1); + + for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) + out(iy) = m_startNode(iy); + + for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) + out(iy+1) = m_startNode(iy+1); + + return out; } // ------------------------------- node-numbers along the right edge ------------------------------- inline ColS FineLayer::nodesRightEdge() { - size_t N = static_cast(m_nh.size()); - ColS nodes(N+1); + size_t nely = static_cast(m_nhy.size()); - // bottom half - for ( size_t i = 0 ; i < (N-1)/2 ; ++i ) - { - if ( m_nh(i)%2 == 0 ) nodes(i) = m_startNode(i) + m_nx/(3*m_nh(i)/2); - else nodes(i) = m_startNode(i) + m_nx/ m_nh(i) ; - } + ColS out(nely+1); - // top half - for ( size_t i = (N-1)/2 ; i < N ; ++i ) - { - if ( m_nh(i)%2 == 0 ) nodes(i) = m_startNode(i) + m_nx/(m_nh(i)/2); - else nodes(i) = m_startNode(i) + m_nx/ m_nh(i) ; - } + for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) + out(iy) = m_startNode(iy) + m_nelx(iy); + + for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) + out(iy+1) = m_startNode(iy+1) + m_nelx(iy); + + return out; +} + +// ---------------------- node-numbers along the bottom edge, without corners ---------------------- - // last node - nodes(N) = m_nnode-1; +inline ColS FineLayer::nodesBottomOpenEdge() +{ + ColS out(m_nelx(0)-1); + + for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix ) + out(ix-1) = m_startNode(0) + ix; + + return out; +} + +// ----------------------- node-numbers along the top edge, without corners ------------------------ + +inline ColS FineLayer::nodesTopOpenEdge() +{ + size_t nely = static_cast(m_nhy.size()); + + ColS out(m_nelx(nely-1)-1); + + for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix ) + out(ix-1) = m_startNode(nely) + ix; + + return out; +} + +// ----------------------- node-numbers along the left edge, without corners ----------------------- + +inline ColS FineLayer::nodesLeftOpenEdge() +{ + size_t nely = static_cast(m_nhy.size()); + + ColS out(nely-1); + + for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) + out(iy-1) = m_startNode(iy); + + for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) + out(iy) = m_startNode(iy+1); - return nodes; + return out; +} + +// ---------------------- node-numbers along the right edge, without corners ----------------------- + +inline ColS FineLayer::nodesRightOpenEdge() +{ + size_t nely = static_cast(m_nhy.size()); + + ColS out(nely-1); + + for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) + out(iy-1) = m_startNode(iy) + m_nelx(iy); + + for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) + out(iy) = m_startNode(iy+1) + m_nelx(iy); + + return out; } // ----------------------------- node-number of the bottom-left corner ----------------------------- inline size_t FineLayer::nodesBottomLeftCorner() { - ColS nodes = nodesBottomEdge(); - - return nodes(0); + return m_startNode(0); } // ---------------------------- node-number of the bottom-right corner ----------------------------- inline size_t FineLayer::nodesBottomRightCorner() { - ColS nodes = nodesBottomEdge(); - - return nodes(nodes.size()-1); + return m_startNode(0) + m_nelx(0); } // ------------------------------ node-number of the top-left corner ------------------------------- inline size_t FineLayer::nodesTopLeftCorner() { - ColS nodes = nodesTopEdge(); + size_t nely = static_cast(m_nhy.size()); - return nodes(0); + return m_startNode(nely); } // ------------------------------ node-number of the top-right corner ------------------------------ inline size_t FineLayer::nodesTopRightCorner() { - ColS nodes = nodesTopEdge(); + size_t nely = static_cast(m_nhy.size()); - return nodes(nodes.size()-1); + return m_startNode(nely) + m_nelx(nely-1); } -// ----------------------------- node-number of the corners (aliases) ------------------------------ +// -------------------------------------------- aliases -------------------------------------------- inline size_t FineLayer::nodesLeftBottomCorner() { return nodesBottomLeftCorner(); } -inline size_t FineLayer::nodesLeftTopCorner() { return nodesTopLeftCorner(); } inline size_t FineLayer::nodesRightBottomCorner() { return nodesBottomRightCorner(); } +inline size_t FineLayer::nodesLeftTopCorner() { return nodesTopLeftCorner(); } inline size_t FineLayer::nodesRightTopCorner() { return nodesTopRightCorner(); } // ------------------------------ node-numbers of periodic node-pairs ------------------------------ inline MatS FineLayer::nodesPeriodic() { - // edges - ColS bot = nodesBottomEdge(); - ColS top = nodesTopEdge(); - ColS rgt = nodesRightEdge(); - ColS lft = nodesLeftEdge(); + // edges (without corners) + ColS bot = nodesBottomOpenEdge(); + ColS top = nodesTopOpenEdge(); + ColS lft = nodesLeftOpenEdge(); + ColS rgt = nodesRightOpenEdge(); // allocate nodal ties - MatS nodes(bot.size()-2 + lft.size()-2 + 3, 2); + // - number of tying per category + size_t tedge = bot.size() + lft.size(); + size_t tnode = 3; + // - allocate + MatS out(tedge+tnode, 2); // counter size_t i = 0; // tie all corners to one corner - nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesBottomRightCorner(); ++i; - nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopRightCorner(); ++i; - nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopLeftCorner(); ++i; + out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesBottomRightCorner(); ++i; + out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopRightCorner(); ++i; + out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopLeftCorner(); ++i; - // tie edges to each other (exclude corners) - for ( auto j = 1 ; j < bot.size()-1 ; ++j ) { nodes(i,0) = bot(j); nodes(i,1) = top(j); ++i; } - for ( auto j = 1 ; j < lft.size()-1 ; ++j ) { nodes(i,0) = lft(j); nodes(i,1) = rgt(j); ++i; } + // tie all corresponding edges to each other + for ( auto j = 0 ; j(nodePer.rows()); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs for ( size_t i = 0 ; i < nper ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) out(nodePer(i,1),j) = out(nodePer(i,0),j); // renumber "out" to be sequential return GooseFEM::Mesh::renumber(out); } // ================================================================================================= }}} // namespace ... // ================================================================================================= #endif - - diff --git a/src/GooseFEM/MeshQuad4.h b/src/GooseFEM/MeshQuad4.h index 35a776f..5960947 100644 --- a/src/GooseFEM/MeshQuad4.h +++ b/src/GooseFEM/MeshQuad4.h @@ -1,122 +1,136 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESHQUAD4_H #define GOOSEFEM_MESHQUAD4_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // ===================================== GooseFEM::Mesh::Quad4 ===================================== namespace GooseFEM { namespace Mesh { namespace Quad4 { // ========================================== REGULAR MESH ========================================= class Regular { private: - size_t m_nx; // number of 'pixels' horizontal direction (length == "m_nx * m_h") - size_t m_ny; // number of 'pixels' vertical direction (length == "m_ny * m_h") - double m_h; // size of the element edge (equal in both directions) + double m_h; // elementary element edge-size (in both directions) + size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h") + size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h") size_t m_nelem; // number of elements size_t m_nnode; // number of nodes size_t m_nne=4; // number of nodes-per-element size_t m_ndim=2; // number of dimensions public: - // mesh with "nx*ny" 'pixels' and edge size "h" - Regular(size_t nx, size_t ny, double h=1.); + // mesh with "nelx*nely" 'elements' of edge size "h" + Regular(size_t nelx, size_t nely, double h=1.); // sizes size_t nelem(); // number of elements size_t nnode(); // number of nodes size_t nne(); // number of nodes-per-element size_t ndim(); // number of dimensions // mesh - MatD coor(); // nodal positions [ nnode , ndim ] - MatS conn(); // connectivity [ nelem , nne ] + MatD coor(); // nodal positions [nnode ,ndim] + MatS conn(); // connectivity [nelem ,nne ] // boundary nodes: edges - ColS nodesBottomEdge(); // nodes along the bottom edge - ColS nodesTopEdge(); // nodes along the top edge - ColS nodesLeftEdge(); // nodes along the left edge - ColS nodesRightEdge(); // nodes along the right edge + ColS nodesBottomEdge(); // node-numbers along the bottom edge + ColS nodesTopEdge(); // node-numbers along the top edge + ColS nodesLeftEdge(); // node-numbers along the left edge + ColS nodesRightEdge(); // node-numbers along the right edge + // boundary nodes: edges, without corners + ColS nodesBottomOpenEdge(); // node-numbers along the bottom edge + ColS nodesTopOpenEdge(); // node-numbers along the top edge + ColS nodesLeftOpenEdge(); // node-numbers along the left edge + ColS nodesRightOpenEdge(); // node-numbers along the right edge // boundary nodes: corners - size_t nodesBottomLeftCorner(); // bottom - left corner node - size_t nodesBottomRightCorner(); // bottom - right corner node - size_t nodesTopLeftCorner(); // top - left corner node - size_t nodesTopRightCorner(); // top - right corner node + size_t nodesBottomLeftCorner(); // node-number of the bottom - left corner + size_t nodesBottomRightCorner(); // node-number of the bottom - right corner + size_t nodesTopLeftCorner(); // node-number of the top - left corner + size_t nodesTopRightCorner(); // node-number of the top - right corner // boundary nodes: corners (aliases) - size_t nodesLeftBottomCorner(); // bottom - left corner node - size_t nodesLeftTopCorner(); // top - left corner node - size_t nodesRightBottomCorner(); // bottom - right corner node - size_t nodesRightTopCorner(); // top - right corner node + size_t nodesLeftBottomCorner(); // alias, see above: nodesBottomLeftCorner + size_t nodesLeftTopCorner(); // alias, see above: nodesBottomRightCorner + size_t nodesRightBottomCorner(); // alias, see above: nodesTopLeftCorner + size_t nodesRightTopCorner(); // alias, see above: nodesTopRightCorner // periodicity MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent) - size_t nodesOrigin(); // bottom-left node, to be used as reference for periodicity + size_t nodesOrigin(); // bottom-left node, used as reference for periodicity MatS dofs(); // DOF-numbers for each component of each node (sequential) - MatS dofsPeriodic(); // DOF-numbers for each component of each node (sequential) + MatS dofsPeriodic(); // ,, for the case that the periodicity if fully eliminated }; // ====================== MESH WITH A FINE LAYER THAT EXPONENTIALLY COARSENS ======================= class FineLayer { private: - double m_h; // elementary element size (middle-layer in x-direction has "L = m_nx * m_h") - size_t m_nx; // number of elements in x-direction - ColS m_nh; // element size in y-direction of each layer - ColS m_startNode; // start node of each layer - ColS m_startElem; // start element of each layer + double m_h; // elementary element edge-size (in all directions) + double m_Lx; // mesh size in "x" + ColS m_nelx; // number of elements in "x" (per el.layer in "y") + ColS m_nnd; // total number of nodes in the main node layer (per nd.layer in "y") + ColS m_nhx, m_nhy; // element size in each direction (per el.layer in "y") + ColI m_refine; // refine direction (-1:no refine, 0:"x") (per el.layer in "y") + ColS m_startElem; // start element (per el.layer in "y") + ColS m_startNode; // start node (per nd.layer in "y") size_t m_nelem; // number of elements size_t m_nnode; // number of nodes size_t m_nne=4; // number of nodes-per-element size_t m_ndim=2; // number of dimensions public: - // mesh with "nx*ny" 'pixels' and edge size "h"; the elements in y-direction are coarsened - FineLayer(size_t nx, size_t ny, double h=1., size_t nfine=1, size_t nskip=0); + // mesh with "nelx*nely" elements of edge size "h"; elements are coarsened in "y"-direction + FineLayer(size_t nelx, size_t nely, double h=1., size_t nfine=1); // sizes size_t nelem(); // number of elements size_t nnode(); // number of nodes size_t nne(); // number of nodes-per-element size_t ndim(); // number of dimensions - size_t shape(size_t i); // actual shape in horizontal and vertical direction + size_t shape(size_t i); // actual shape in a certain direction // mesh - MatD coor(); // nodal positions [ nnode , ndim ] - MatS conn(); // connectivity [ nelem , nne ] - // boundary nodes: edges + MatD coor(); // nodal positions [nnode ,ndim] + MatS conn(); // connectivity [nelem ,nne ] + // element sets ColS elementsMiddleLayer(); // elements in the middle, fine, layer - ColS nodesBottomEdge(); // nodes along the bottom edge - ColS nodesTopEdge(); // nodes along the top edge - ColS nodesLeftEdge(); // nodes along the left edge - ColS nodesRightEdge(); // nodes along the right edge + // boundary nodes: edges + ColS nodesBottomEdge(); // node-numbers along the bottom edge + ColS nodesTopEdge(); // node-numbers along the top edge + ColS nodesLeftEdge(); // node-numbers along the left edge + ColS nodesRightEdge(); // node-numbers along the right edge + // boundary nodes: edges, without corners + ColS nodesBottomOpenEdge(); // node-numbers along the bottom edge + ColS nodesTopOpenEdge(); // node-numbers along the top edge + ColS nodesLeftOpenEdge(); // node-numbers along the left edge + ColS nodesRightOpenEdge(); // node-numbers along the right edge // boundary nodes: corners - size_t nodesBottomLeftCorner(); // bottom - left corner node - size_t nodesBottomRightCorner(); // bottom - right corner node - size_t nodesTopLeftCorner(); // top - left corner node - size_t nodesTopRightCorner(); // top - right corner node + size_t nodesBottomLeftCorner(); // node-number of the bottom - left corner + size_t nodesBottomRightCorner(); // node-number of the bottom - right corner + size_t nodesTopLeftCorner(); // node-number of the top - left corner + size_t nodesTopRightCorner(); // node-number of the top - right corner // boundary nodes: corners (aliases) - size_t nodesLeftBottomCorner(); // bottom - left corner node - size_t nodesLeftTopCorner(); // top - left corner node - size_t nodesRightBottomCorner(); // bottom - right corner node - size_t nodesRightTopCorner(); // top - right corner node + size_t nodesLeftBottomCorner(); // alias, see above: nodesBottomLeftCorner + size_t nodesLeftTopCorner(); // alias, see above: nodesBottomRightCorner + size_t nodesRightBottomCorner(); // alias, see above: nodesTopLeftCorner + size_t nodesRightTopCorner(); // alias, see above: nodesTopRightCorner // periodicity MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent) - size_t nodesOrigin(); // bottom-left node, to be used as reference for periodicity + size_t nodesOrigin(); // bottom-left node, used as reference for periodicity MatS dofs(); // DOF-numbers for each component of each node (sequential) - MatS dofsPeriodic(); // DOF-numbers for each component of each node (sequential) + MatS dofsPeriodic(); // ,, for the case that the periodicity if fully eliminated }; // ================================================================================================= }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/MeshTri3.cpp b/src/GooseFEM/MeshTri3.cpp index 4206e23..bb0840d 100644 --- a/src/GooseFEM/MeshTri3.cpp +++ b/src/GooseFEM/MeshTri3.cpp @@ -1,553 +1,606 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESHTRI3_CPP #define GOOSEFEM_MESHTRI3_CPP // ------------------------------------------------------------------------------------------------- #include "MeshTri3.h" // ===================================== GooseFEM::Mesh::Tri3 ===================================== namespace GooseFEM { namespace Mesh { namespace Tri3 { // ===================================== CLASS - REGULAR MESH ====================================== // ------------------------------------------ constructor ------------------------------------------ -inline Regular::Regular(size_t nx, size_t ny, double h): m_nx(nx), m_ny(ny), m_h(h) +inline Regular::Regular(size_t nelx, size_t nely, double h): +m_h(h), m_nelx(nelx), m_nely(nely) { - assert( m_nx >= 1 ); - assert( m_ny >= 1 ); + assert( m_nelx >= 1 ); + assert( m_nely >= 1 ); - m_nnode = (m_nx+1) * (m_ny+1) ; - m_nelem = m_nx * m_ny * 2; + m_nnode = (m_nelx+1) * (m_nely+1) ; + m_nelem = m_nelx * m_nely * 2; } // -------------------------------------- number of elements --------------------------------------- inline size_t Regular::nelem() { return m_nelem; } // ---------------------------------------- number of nodes ---------------------------------------- inline size_t Regular::nnode() { return m_nnode; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t Regular::nne() { return m_nne; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t Regular::ndim() { return m_ndim; } // --------------------------------- coordinates (nodal positions) --------------------------------- inline MatD Regular::coor() { - MatD coor( m_nnode , m_ndim ); + MatD out(m_nnode , m_ndim); - ColD x = ColD::LinSpaced( m_nx+1 , 0.0 , m_h * static_cast(m_nx) ); - ColD y = ColD::LinSpaced( m_ny+1 , 0.0 , m_h * static_cast(m_ny) ); + ColD x = ColD::LinSpaced( m_nelx+1 , 0.0 , m_h * static_cast(m_nelx) ); + ColD y = ColD::LinSpaced( m_nely+1 , 0.0 , m_h * static_cast(m_nely) ); size_t inode = 0; - for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy ) { - for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix ) { - coor(inode,0) = x(ix); - coor(inode,1) = y(iy); + for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) { + for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) { + out(inode,0) = x(ix); + out(inode,1) = y(iy); ++inode; } } - return coor; + return out; } // ---------------------------- connectivity (node-numbers per element) ---------------------------- inline MatS Regular::conn() { - MatS conn( m_nelem , m_nne ); + MatS out(m_nelem , m_nne); size_t ielem = 0; - for ( size_t iy = 0 ; iy < m_ny ; ++iy ) { - for ( size_t ix = 0 ; ix < m_nx ; ++ix ) { - conn(ielem,0) = (iy+0)*(m_nx+1) + (ix+0); - conn(ielem,1) = (iy+0)*(m_nx+1) + (ix+1); - conn(ielem,2) = (iy+1)*(m_nx+1) + (ix+0); + for ( size_t iy = 0 ; iy < m_nely ; ++iy ) { + for ( size_t ix = 0 ; ix < m_nelx ; ++ix ) { + out(ielem,0) = (iy )*(m_nelx+1) + (ix ); + out(ielem,1) = (iy )*(m_nelx+1) + (ix+1); + out(ielem,2) = (iy+1)*(m_nelx+1) + (ix ); ++ielem; - conn(ielem,0) = (iy+0)*(m_nx+1) + (ix+1); - conn(ielem,1) = (iy+1)*(m_nx+1) + (ix+1); - conn(ielem,2) = (iy+1)*(m_nx+1) + (ix+0); + out(ielem,0) = (iy )*(m_nelx+1) + (ix+1); + out(ielem,1) = (iy+1)*(m_nelx+1) + (ix+1); + out(ielem,2) = (iy+1)*(m_nelx+1) + (ix ); ++ielem; } } - return conn; + return out; } // ------------------------------ node-numbers along the bottom edge ------------------------------- inline ColS Regular::nodesBottomEdge() { - ColS nodes(m_nx+1); + ColS out(m_nelx+1); - for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix ) - nodes(ix) = ix; + for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) + out(ix) = ix; - return nodes; + return out; } // -------------------------------- node-numbers along the top edge -------------------------------- inline ColS Regular::nodesTopEdge() { - ColS nodes(m_nx+1); + ColS out(m_nelx+1); - for ( size_t ix = 0 ; ix < m_nx+1 ; ++ix ) - nodes(ix) = ix + m_ny*(m_nx+1); + for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) + out(ix) = ix + m_nely*(m_nelx+1); - return nodes; + return out; } -// ----------------------------- node-numbers along the node left edge ----------------------------- +// ------------------------------- node-numbers along the left edge -------------------------------- inline ColS Regular::nodesLeftEdge() { - ColS nodes(m_ny+1); + ColS out(m_nely+1); - for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy ) - nodes(iy) = iy*(m_nx+1); + for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) + out(iy) = iy*(m_nelx+1); - return nodes; + return out; } // ------------------------------- node-numbers along the right edge ------------------------------- inline ColS Regular::nodesRightEdge() { - ColS nodes(m_ny+1); + ColS out(m_nely+1); + + for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) + out(iy) = iy*(m_nelx+1) + m_nelx; + + return out; +} + +// ---------------------- node-numbers along the bottom edge, without corners ---------------------- + +inline ColS Regular::nodesBottomOpenEdge() +{ + ColS out(m_nelx-1); - for ( size_t iy = 0 ; iy < m_ny+1 ; ++iy ) - nodes(iy) = iy*(m_nx+1) + m_nx; + for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) + out(ix-1) = ix; - return nodes; + return out; +} + +// ----------------------- node-numbers along the top edge, without corners ------------------------ + +inline ColS Regular::nodesTopOpenEdge() +{ + ColS out(m_nelx-1); + + for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) + out(ix-1) = ix + m_nely*(m_nelx+1); + + return out; +} + +// ----------------------- node-numbers along the left edge, without corners ----------------------- + +inline ColS Regular::nodesLeftOpenEdge() +{ + ColS out(m_nely-1); + + for ( size_t iy = 1 ; iy < m_nely ; ++iy ) + out(iy-1) = iy*(m_nelx+1); + + return out; +} + +// ---------------------- node-numbers along the right edge, without corners ----------------------- + +inline ColS Regular::nodesRightOpenEdge() +{ + ColS out(m_nely-1); + + for ( size_t iy = 1 ; iy < m_nely ; ++iy ) + out(iy-1) = iy*(m_nelx+1) + m_nelx; + + return out; } // ----------------------------- node-number of the bottom-left corner ----------------------------- inline size_t Regular::nodesBottomLeftCorner() { return 0; } // ---------------------------- node-number of the bottom-right corner ----------------------------- inline size_t Regular::nodesBottomRightCorner() { - return m_nx; + return m_nelx; } // ------------------------------ node-number of the top-left corner ------------------------------- inline size_t Regular::nodesTopLeftCorner() { - return m_ny*(m_nx+1); + return m_nely*(m_nelx+1); } // ------------------------------ node-number of the top-right corner ------------------------------ inline size_t Regular::nodesTopRightCorner() { - return (m_ny+1)*(m_nx+1) - 1; + return m_nely*(m_nelx+1) + m_nelx; } // ----------------------------- node-number of the corners (aliases) ------------------------------ inline size_t Regular::nodesLeftBottomCorner() { return nodesBottomLeftCorner(); } inline size_t Regular::nodesLeftTopCorner() { return nodesTopLeftCorner(); } inline size_t Regular::nodesRightBottomCorner() { return nodesBottomRightCorner(); } inline size_t Regular::nodesRightTopCorner() { return nodesTopRightCorner(); } // ------------------------------ node-numbers of periodic node-pairs ------------------------------ inline MatS Regular::nodesPeriodic() { - // edges - ColS bot = nodesBottomEdge(); - ColS top = nodesTopEdge(); - ColS rgt = nodesRightEdge(); - ColS lft = nodesLeftEdge(); + // edges (without corners) + ColS bot = nodesBottomOpenEdge(); + ColS top = nodesTopOpenEdge(); + ColS lft = nodesLeftOpenEdge(); + ColS rgt = nodesRightOpenEdge(); // allocate nodal ties - MatS nodes(bot.size()-2 + lft.size()-2 + 3, 2); + // - number of tying per category + size_t tedge = bot.size() + lft.size(); + size_t tnode = 3; + // - allocate + MatS out(tedge+tnode, 2); // counter size_t i = 0; // tie all corners to one corner - nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesBottomRightCorner(); ++i; - nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopRightCorner(); ++i; - nodes(i,0) = nodesBottomLeftCorner(); nodes(i,1) = nodesTopLeftCorner(); ++i; + out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesBottomRightCorner(); ++i; + out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopRightCorner(); ++i; + out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopLeftCorner(); ++i; - // tie edges to each other (exclude corners) - for ( auto j = 1 ; j < bot.size()-1 ; ++j ) { nodes(i,0) = bot(j); nodes(i,1) = top(j); ++i; } - for ( auto j = 1 ; j < lft.size()-1 ; ++j ) { nodes(i,0) = lft(j); nodes(i,1) = rgt(j); ++i; } + // tie all corresponding edges to each other + for ( auto j = 0 ; j(nodePer.rows()); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs for ( size_t i = 0 ; i < nper ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) out(nodePer(i,1),j) = out(nodePer(i,0),j); // renumber "out" to be sequential return GooseFEM::Mesh::renumber(out); } // ========================================= MESH ANALYSIS ========================================= // ------------------------------ get the orientation of each element ------------------------------ inline ColI getOrientation(const MatD &coor, const MatS &conn) { assert( conn.cols() == 3 ); assert( coor.cols() == 2 ); Eigen::Vector2d v1,v2; double k; size_t nelem = static_cast( conn.rows() ); ColI out( nelem ); for ( size_t ielem = 0 ; ielem < nelem ; ++ielem ) { v1 = coor.row( conn(ielem,0) ) - coor.row( conn(ielem,1) ); v2 = coor.row( conn(ielem,2) ) - coor.row( conn(ielem,1) ); k = v1(0) * v2(1) - v2(0) * v1(1); if ( k < 0 ) out( ielem ) = -1; else out( ielem ) = +1; } return out; } // ------------------------------ set the orientation of each element ------------------------------ inline MatS setOrientation(const MatD &coor, const MatS &conn, int orientation) { assert( conn.cols() == 3 ); assert( coor.cols() == 2 ); assert( orientation == -1 || orientation == +1 ); ColI val = getOrientation(coor, conn); return setOrientation( coor, conn, val, orientation ); } // -------------------- set the orientation of each element to a certain value --------------------- inline MatS setOrientation(const MatD &coor, const MatS &conn, const ColI &val, int orientation) { assert( conn.cols() == 3 ); assert( coor.cols() == 2 ); assert( conn.rows() == val.size() ); assert( orientation == -1 || orientation == +1 ); // avoid compiler warning UNUSED(coor); size_t nelem = static_cast(conn.rows()); MatS out = conn; for ( size_t ielem = 0 ; ielem < nelem ; ++ielem ) if ( ( orientation == -1 and val(ielem) > 0 ) or ( orientation == +1 and val(ielem) < 0 ) ) std::swap( out(ielem,2) , out(ielem,1) ); return out; } // ===================================== RE-TRIANGULATION API ====================================== inline MatS retriangulate(const MatD &coor, const MatS &conn, int orientation) { // get the orientation of all elements ColI dir = getOrientation( coor, conn ); // check the orientation bool eq = std::abs( dir.sum() ) == conn.rows(); // new connectivity MatS out; // perform re-triangulation // - use "TriUpdate" if ( eq ) { TriUpdate obj(coor,conn); obj.eval(); out = obj.conn(); } // - using TriCompute else { throw std::runtime_error("Work-in-progress, has to be re-triangulated using 'TriCompute'"); } return setOrientation(coor,out,orientation); } // ================ RE-TRIANGULATION SUPPORT CLASS - UPDATE EXISTING TRIANGULATION ================= // ------------------------------------------ constructor ------------------------------------------ inline TriUpdate::TriUpdate(const MatD &coor, const MatS &conn): m_conn(conn), m_coor(coor) { assert( conn.cols() == 3 ); assert( coor.cols() == 2 ); // store shapes m_nnode = static_cast( coor.rows() ); m_ndim = static_cast( coor.cols() ); m_nelem = static_cast( conn.rows() ); m_nne = static_cast( conn.cols() ); // resize internal arrays m_elem.resize(2); m_node.resize(4); // set default to out-of-bounds, to make clear that nothing happened yet m_elem.setConstant(m_nelem); m_node.setConstant(m_nnode); edge(); } // -------------------------- compute neighbors per edge of all elements --------------------------- inline void TriUpdate::edge() { m_edge.resize( m_nelem , m_nne ); m_edge.setConstant( m_nelem ); // signal that nothing has been set std::vector idx = {0,1,2}; // lists to convert connectivity -> edges std::vector jdx = {1,2,0}; // lists to convert connectivity -> edges std::vector edge; edge.reserve( m_nelem*idx.size() ); for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t i = 0 ; i < m_nne ; ++i ) edge.push_back( Edge( m_conn(e,idx[i]), m_conn(e,jdx[i]) , e , i , true ) ); std::sort( edge.begin() , edge.end() , Edge_sort ); for ( size_t i = 0 ; i < edge.size()-1 ; ++i ) { if ( edge[i].n1 == edge[i+1].n1 and edge[i].n2 == edge[i+1].n2 ) { m_edge( edge[i ].elem , edge[i ].edge ) = edge[i+1].elem; m_edge( edge[i+1].elem , edge[i+1].edge ) = edge[i ].elem; } } } // ---------------------------- update edges around renumbered elements ---------------------------- inline void TriUpdate::chedge(size_t edge, size_t old_elem, size_t new_elem) { size_t m; size_t neigh = m_edge( old_elem , edge ); if ( neigh >= m_nelem ) return; for ( m = 0 ; m < m_nne ; ++m ) if ( m_edge( neigh , m ) == old_elem ) break; m_edge( neigh , m ) = new_elem; } // --------------------------------- re-triangulate the full mesh ---------------------------------- inline bool TriUpdate::eval() { bool change = false; while ( increment() ) { change = true; } return change; } // ----------------------------------- one re-triangulation step ----------------------------------- inline bool TriUpdate::increment() { size_t ielem,jelem,iedge,jedge; double phi1,phi2; ColS c(4); ColS n(4); // loop over all elements for ( ielem = 0 ; ielem < m_nelem ; ++ielem ) { // loop over all edges for ( iedge = 0 ; iedge < m_nne ; ++iedge ) { // only proceed if the edge is shared with another element if ( m_edge(ielem,iedge) >= m_nelem ) continue; // read "jelem" jelem = m_edge(ielem,iedge); // find the edge involved for "jelem" for ( jedge=0; jedge M_PI ) { // update connectivity m_conn(ielem,0) = c(0); m_conn(ielem,1) = c(3); m_conn(ielem,2) = c(2); m_conn(jelem,0) = c(0); m_conn(jelem,1) = c(1); m_conn(jelem,2) = c(3); // change list with neighbors for the elements around (only two neighbors change) if ( iedge==0 ) { chedge(2,ielem,jelem); } else if ( iedge==1 ) { chedge(0,ielem,jelem); } else if ( iedge==2 ) { chedge(1,ielem,jelem); } if ( jedge==0 ) { chedge(2,jelem,ielem); } else if ( jedge==1 ) { chedge(0,jelem,ielem); } else if ( jedge==2 ) { chedge(1,jelem,ielem); } // convert to four static nodes if ( iedge==0 ) { n(0)=m_edge(ielem,2); n(3)=m_edge(ielem,1); } else if ( iedge==1 ) { n(0)=m_edge(ielem,0); n(3)=m_edge(ielem,2); } else if ( iedge==2 ) { n(0)=m_edge(ielem,1); n(3)=m_edge(ielem,0); } if ( jedge==0 ) { n(1)=m_edge(jelem,1); n(2)=m_edge(jelem,2); } else if ( jedge==1 ) { n(1)=m_edge(jelem,2); n(2)=m_edge(jelem,0); } else if ( jedge==2 ) { n(1)=m_edge(jelem,0); n(2)=m_edge(jelem,1); } // store the neighbors for the changed elements m_edge(ielem,0) = jelem; m_edge(jelem,0) = n(0) ; m_edge(ielem,1) = n(2) ; m_edge(jelem,1) = n(1) ; m_edge(ielem,2) = n(3) ; m_edge(jelem,2) = ielem; // store information for transfer algorithm m_node = c; m_elem(0) = ielem; m_elem(1) = jelem; return true; } } } return false; } // ================================ SUPPORT CLASS - EDGE DEFINITION ================================ // ------------------------------------------ constructor ------------------------------------------ Edge::Edge(size_t i, size_t j, size_t el, size_t ed, bool sort): n1(i), n2(j), elem(el), edge(ed) { if ( sort && n1>n2 ) std::swap(n1,n2); } // --------------------------------------- compare two edges --------------------------------------- inline bool Edge_cmp(Edge a, Edge b) { if ( a.n1 == b.n1 and a.n2 == b.n2 ) return true; return false; } // ----------------------- sort edges by comparing the first and second node ----------------------- inline bool Edge_sort(Edge a, Edge b) { if ( a.n1 < b.n1 or a.n2 < b.n2 ) return true; return false; } // ================================================================================================= }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/MeshTri3.h b/src/GooseFEM/MeshTri3.h index 245d8e0..baffb16 100644 --- a/src/GooseFEM/MeshTri3.h +++ b/src/GooseFEM/MeshTri3.h @@ -1,141 +1,146 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESHTRI3_H #define GOOSEFEM_MESHTRI3_H // ------------------------------------------------------------------------------------------------- #include "Mesh.h" // ===================================== GooseFEM::Mesh::Tri3 ===================================== namespace GooseFEM { namespace Mesh { namespace Tri3 { // ========================================== REGULAR MESH ========================================= class Regular { private: - size_t m_nx; // number of 'pixels' horizontal direction (length == "m_nx * m_h") - size_t m_ny; // number of 'pixels' vertical direction (length == "m_ny * m_h") - double m_h; // size of the element edge (equal in both directions) - size_t m_nelem; // number of elements - size_t m_nnode; // number of nodes - size_t m_nne=3; // number of nodes-per-element - size_t m_ndim=2; // number of dimensions + double m_h; // elementary element edge-size (in both directions) + size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h") + size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h") + size_t m_nelem; // number of elements + size_t m_nnode; // number of nodes + size_t m_nne=3; // number of nodes-per-element + size_t m_ndim=2; // number of dimensions public: - // mesh with "nx" pixels in horizontal direction, "ny" in vertical direction and "h" the edge size - Regular(size_t nx, size_t ny, double h=1.); + // mesh with "2*nelx*nely" 'elements' of edge size "h" + Regular(size_t nelx, size_t nely, double h=1.); // sizes - size_t nelem(); // number of elements - size_t nnode(); // number of nodes - size_t nne(); // number of nodes-per-element - size_t ndim(); // number of dimensions + size_t nelem(); // number of elements + size_t nnode(); // number of nodes + size_t nne(); // number of nodes-per-element + size_t ndim(); // number of dimensions // mesh - MatD coor(); // nodal positions [ nnode , ndim ] - MatS conn(); // connectivity [ nelem , nne ] + MatD coor(); // nodal positions [nnode ,ndim] + MatS conn(); // connectivity [nelem ,nne ] // boundary nodes: edges - ColS nodesBottomEdge(); // nodes along the bottom edge - ColS nodesTopEdge(); // nodes along the top edge - ColS nodesLeftEdge(); // nodes along the left edge - ColS nodesRightEdge(); // nodes along the right edge + ColS nodesBottomEdge(); // node-numbers along the bottom edge + ColS nodesTopEdge(); // node-numbers along the top edge + ColS nodesLeftEdge(); // node-numbers along the left edge + ColS nodesRightEdge(); // node-numbers along the right edge + // boundary nodes: edges, without corners + ColS nodesBottomOpenEdge(); // node-numbers along the bottom edge + ColS nodesTopOpenEdge(); // node-numbers along the top edge + ColS nodesLeftOpenEdge(); // node-numbers along the left edge + ColS nodesRightOpenEdge(); // node-numbers along the right edge // boundary nodes: corners - size_t nodesBottomLeftCorner(); // bottom - left corner node - size_t nodesBottomRightCorner(); // bottom - right corner node - size_t nodesTopLeftCorner(); // top - left corner node - size_t nodesTopRightCorner(); // top - right corner node + size_t nodesBottomLeftCorner(); // node-number of the bottom - left corner + size_t nodesBottomRightCorner(); // node-number of the bottom - right corner + size_t nodesTopLeftCorner(); // node-number of the top - left corner + size_t nodesTopRightCorner(); // node-number of the top - right corner // boundary nodes: corners (aliases) - size_t nodesLeftBottomCorner(); // bottom - left corner node - size_t nodesLeftTopCorner(); // top - left corner node - size_t nodesRightBottomCorner(); // bottom - right corner node - size_t nodesRightTopCorner(); // top - right corner node + size_t nodesLeftBottomCorner(); // alias, see above: nodesBottomLeftCorner + size_t nodesLeftTopCorner(); // alias, see above: nodesBottomRightCorner + size_t nodesRightBottomCorner(); // alias, see above: nodesTopLeftCorner + size_t nodesRightTopCorner(); // alias, see above: nodesTopRightCorner // periodicity - MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent) - size_t nodesOrigin(); // bottom-left node, to be used as reference for periodicity - MatS dofs(); // DOF-numbers for each component of each node (sequential) - MatS dofsPeriodic(); // DOF-numbers for each component of each node (sequential) + MatS nodesPeriodic(); // periodic node pairs [:,2]: (independent, dependent) + size_t nodesOrigin(); // bottom-left node, used as reference for periodicity + MatS dofs(); // DOF-numbers for each component of each node (sequential) + MatS dofsPeriodic(); // ,, for the case that the periodicity if fully eliminated }; // ========================================= MESH ANALYSIS ========================================= // read / set the orientation (-1 / +1) of all triangles inline ColI getOrientation(const MatD &coor, const MatS &conn ); inline MatS setOrientation(const MatD &coor, const MatS &conn, int orientation=-1); inline MatS setOrientation(const MatD &coor, const MatS &conn, const ColI &val, int orientation=-1); // ======================================= RE-TRIANGULATION ======================================== // simple interface to compute the full re-triangulation; it uses, depending on the input mesh: // (1) the minimal evasive "TriUpdate" // (2) the more rigorous "TriCompute" inline MatS retriangulate(const MatD &coor, const MatS &conn, int orientation=-1); // ------------------------- support class - update existing triangulation ------------------------- // minimal evasive re-triangulation which only flips edges of the existing connectivity class TriUpdate { private: MatS m_edge; // the element that neighbors along each edge (m_nelem: no neighbor) MatS m_conn; // connectivity (updated) MatD m_coor; // nodal positions (does not change) size_t m_nelem; // #elements size_t m_nnode; // #nodes size_t m_nne; // #nodes-per-element size_t m_ndim; // #dimensions ColS m_elem; // the two elements involved in the last element change (see below) ColS m_node; // the four nodes involved in the last element change (see below) // old: m_elem(0) = [ m_node(0) , m_node(1) , m_node(2) ] // m_elem(1) = [ m_node(1) , m_node(3) , m_node(2) ] // new: m_elem(0) = [ m_node(0) , m_node(3) , m_node(2) ] // m_elem(1) = [ m_node(0) , m_node(1) , m_node(3) ] // compute neighbors per edge of all elements void edge(); // update edges around renumbered elements void chedge(size_t edge, size_t old_elem, size_t new_elem); public: TriUpdate(){}; TriUpdate(const MatD &coor, const MatS &conn); bool eval (); // re-triangulate the full mesh (returns "true" if changed) bool increment(); // one re-triangulation step (returns "true" if changed) MatS conn () { return m_conn; }; // return (new) connectivity MatS ch_elem () { return m_elem; }; // return element involved in last element change MatS ch_node () { return m_node; }; // return nodes involved in last element change }; // ================================ SUPPORT CLASS - EDGE DEFINITION ================================ // support class to allow the storage of a list of edges class Edge { public: size_t n1 ; // node 1 (edge from node 1-2) size_t n2 ; // node 2 (edge from node 1-2) size_t elem; // element to which the edge belong size_t edge; // edge index within the element (e.g. edge==1 -> n1=conn(0,elem), n2=conn(1,elem)) Edge(){}; Edge(size_t i, size_t j, size_t el, size_t ed, bool sort=false); }; // ------------------------------------------------------------------------------------------------- // compare edges inline bool Edge_cmp (Edge a, Edge b); // check equality inline bool Edge_sort(Edge a, Edge b); // check if "a" is smaller than "b" in terms of node-numbers // ================================================================================================= }}} // namespace ... #endif