Merge pull request #62 from gridap/hex_to_tet_conversions

Hex to tet conversions

.travis.yml

@@ -20,6 +20,8 @@ jobs:
       julia: 1.1
       os: linux
       script:
+        - julia --project=docs/ -e 'using Pkg;
+                                    Pkg.add(PackageSpec(name="UnstructuredGrids", rev="master"))'
         - julia --project=docs/ -e 'using Pkg;
                                     Pkg.develop(PackageSpec(path=pwd()));
                                     Pkg.instantiate()'

Manifest.toml

@@ -197,7 +197,9 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
 [[UnstructuredGrids]]
 deps = ["WriteVTK"]
-git-tree-sha1 = "8d383ac05dcbbaa3b1222b441370e485e2abc831"
+git-tree-sha1 = "e780e3e9a2895f689cd68465221b866a6820112e"
+repo-rev = "master"
+repo-url = "https://github.com/gridap/UnstructuredGrids.jl.git"
 uuid = "9cd2e232-b78e-5341-8154-33bf4eafe965"
 version = "0.2.0"
 

src/CellValues/CellValuesGallery.jl

@@ -20,6 +20,7 @@ import Base: eltype
 import Base: size
 import Base: getindex
 import Base: IndexStyle
+import Gridap: compress
 
 struct CellValueFromArray{T,N,V<:AbstractArray{T,N}} <: IndexCellValue{T,N}
   v::V
@@ -58,6 +59,10 @@ size(self::CellVectorFromDataAndPtrs) = (length(self.ptrs)-1,)
 
 IndexStyle(::Type{CellVectorFromDataAndPtrs{T,V,P}}) where {T,V,P} = IndexLinear()
 
+function compress(cv::CellVectorFromDataAndPtrs)
+  (cv.data, cv.ptrs)
+end
+
 struct CellVectorFromDataAndStride{T,V} <: IndexCellArray{T,1,CachedSubVector{T,V},1}
   data::V
   stride::Int

src/CellValues/CellValuesOperations.jl

@@ -1,13 +1,15 @@
 module CellValuesOperations
 
 using Gridap
+using Gridap.Helpers
 
 using Gridap.Kernels: CellSumKernel
 using Gridap.Kernels: CellNewAxisKernel
 
 export cellsum
 export cellmean
 export cellnewaxis
+export compress
 import Base: +,-,*,/,\, ==, ≈
 import LinearAlgebra: inv, det
 import TensorValues: inner, outer, meas
@@ -111,4 +113,8 @@ function cellnewaxis(ca::CellMap;dim::Int)
   apply(k,ca)
 end
 
+function compress(cv::CellVector)
+  @notimplemented
+end
+
 end # module

src/CellValues/files.jl

@@ -29,12 +29,12 @@ include("MapApply.jl")
 include("CellMapApply.jl")
 @reexport using Gridap.CellMapApply
 
-include("CellValuesGallery.jl")
-#@reexport using Gridap.CellValuesGallery
-
 include("CellValuesOperations.jl")
 @reexport using Gridap.CellValuesOperations
 
+include("CellValuesGallery.jl")
+#@reexport using Gridap.CellValuesGallery
+
 include("CellValuesAppend.jl")
 @reexport using Gridap.CellValuesAppend
 

src/Geometry/GeometryWrappers.jl

@@ -23,9 +23,41 @@ RefCell(polytope::Polytope) = _ref_cell_from_polytope(polytope)
 Construct a Grid from a polytope
 """
 function Grid(polytope::Polytope{D},dim::Int) where D
+  if dim < D
+    _grid_d(polytope,dim)
+  elseif dim==D
+    _grid_D(polytope)
+  else
+    @unreachable
+  end
+end
+
+# Helpers
+
+function _grid_D(polytope::Polytope{D}) where D
+
+  order = 1
+  reffe = LagrangianRefFE{D,Float64}(polytope,order)
+  points = reffe.dofbasis.nodes
+
+  cell_to_nodes = [ [i for i in 1:length(points)], ]
+
+  cells_data,cells_ptrs = generate_data_and_ptrs(cell_to_nodes)
+  extrusion = polytope.extrusion.array.data
+
+  l = 1
+  co = ConstantCellValue(order,l)
+  ct = ConstantCellValue(extrusion,l)
+
+  UnstructuredGrid(points,cells_data,cells_ptrs,ct,co)
+
+end
+
+function _grid_d(polytope::Polytope{D},dim::Int) where D
   @assert dim < D
 
-  reffe = LagrangianRefFE{D,Float64}(polytope, 1)
+  order = 1
+  reffe = LagrangianRefFE{D,Float64}(polytope,order)
   points = reffe.dofbasis.nodes
 
   dim_to_jface_to_vertices, dim_to_jface_to_code = _faces(polytope)
@@ -49,7 +81,7 @@ function Grid(polytope::Polytope{D},dim::Int) where D
 
 end
 
-# Helpers
+
 
 function _ref_cell_from_polytope(polytope::Polytope{D}) where D
 

src/Geometry/Simplexify.jl

@@ -0,0 +1,254 @@
+module Simplexify
+
+using Gridap
+using Gridap.Helpers
+using UnstructuredGrids.Kernels: refine_grid_connectivity
+using UnstructuredGrids.Kernels: generate_tface_to_face
+using Gridap.DiscreteModels: DiscreteModelFromData
+
+export simplexify
+
+function simplexify(grid::Grid{D,Z}) where {D,Z}
+  ugrid = UnstructuredGrid(grid)
+
+  ltcell_to_lpoints = _generate_ltcell_to_lpoints(celltypes(ugrid))
+
+  tcells_data, tcells_ptrs = refine_grid_connectivity(
+    ugrid.cells_data, ugrid.cells_ptrs, ltcell_to_lpoints)
+
+  order = 1
+  _check_order(order,cellorders(ugrid))
+
+  ntcells = length(tcells_ptrs) - 1
+  ex = tuple(fill(TET_AXIS,Z)...)
+  _ct = ConstantCellValue(ex,ntcells)
+  _co = ConstantCellValue(order,ntcells)
+
+  UnstructuredGrid(
+    points(ugrid),
+    tcells_data,
+    tcells_ptrs,
+    _ct,
+    _co)
+
+end
+
+function simplexify(model::DiscreteModel{D}) where D
+
+  grid = Grid(model)
+
+  graph = FullGridGraph(model)
+
+  facelabels = FaceLabels(model)
+
+  tgrid = simplexify(grid)
+
+  tgraph = FullGridGraph(tgrid)
+
+  tgrids = _generate_tgrids(tgrid,tgraph)
+
+  ltcell_to_lnodes = _generate_ltcell_to_lpoints(celltypes(grid))
+
+  tfacelabels = _generate_tfacelabels(
+    tgrids,grid,facelabels,tgraph,graph,ltcell_to_lnodes)
+
+  DiscreteModelFromData(tgrids,tgraph,tfacelabels)
+
+end
+
+# Helpers
+
+const UNSET_ID = 0
+
+function _generate_tfacelabels(
+  tgrids,grid,facelabels,tgraph,graph,ltcell_to_lnodes)
+
+  dim_to_tface_to_label = [ fill(UNSET_ID,ncells(g)) for g in tgrids ]
+
+  dim_to_face_to_label = facelabels.dim_to_nface_to_label
+
+  tgrid = tgrids[end]
+  tpolytope = CellPolytopes(tgrid).value
+  polytope = CellPolytopes(grid).value
+
+  _fill_dim_to_tface_to_label!(
+    dim_to_tface_to_label,
+    dim_to_face_to_label,
+    tpolytope,
+    polytope,
+    tgraph,
+    graph,
+    ltcell_to_lnodes)
+
+  FaceLabels(
+    dim_to_tface_to_label, facelabels.tag_to_labels, facelabels.tag_to_name)
+
+end
+
+function _fill_dim_to_tface_to_label!(
+  dim_to_tface_to_label,
+  dim_to_face_to_label,
+  tpolytope,
+  polytope,
+  tgraph,
+  graph,
+  ltcell_to_lnodes)
+
+  D = length(dim_to_face_to_label)-1
+  d = 0
+  dim_to_tface_to_label[d+1] = dim_to_face_to_label[d+1]
+
+  trefcell = RefCell(tpolytope)
+  refcell = RefCell(polytope)
+
+  for d in 1:(D-1)
+
+    cell_to_faces = connections(graph,D,d)
+    cell_to_faces_data, cell_to_faces_ptrs = compress(cell_to_faces)
+
+    tcell_to_tfaces = connections(tgraph,D,d)
+
+    ntfaces = maximum(tcell_to_tfaces.data)
+
+    ltface_to_ltnodes = trefcell.faces[d+1]
+    lface_to_lnodes = refcell.faces[d+1]
+
+    tface_to_face = generate_tface_to_face(
+      cell_to_faces_data,
+      cell_to_faces_ptrs,
+      tcell_to_tfaces.data,
+      tcell_to_tfaces.ptrs,
+      ltcell_to_lnodes,
+      ltface_to_ltnodes,
+      lface_to_lnodes,
+      ntfaces)
+
+    _update_labels!(
+      dim_to_tface_to_label[d+1],dim_to_face_to_label[d+1],tface_to_face)
+
+  end
+
+  d = D
+  ncells = length(dim_to_face_to_label[d+1])
+  nltcells = length(ltcell_to_lnodes)
+  tcell_to_cell = _generate_tcell_to_cell(ncells,nltcells)
+  _update_labels!(
+    dim_to_tface_to_label[d+1],dim_to_face_to_label[d+1],tcell_to_cell)
+
+  for d = 1:(D-1)
+
+    for j in (d+1):D
+
+      dface_to_jfaces = connections(tgraph,d,j)
+      dface_to_label = dim_to_tface_to_label[d+1]
+      jface_to_label = dim_to_tface_to_label[j+1]
+      _fix_dface_to_label!(dface_to_label,jface_to_label,dface_to_jfaces)
+
+    end
+
+  end
+
+end
+
+function _fix_dface_to_label!(dface_to_label,jface_to_label,dface_to_jfaces)
+
+  ndfaces = length(dface_to_label)
+  @assert ndfaces == length(dface_to_jfaces)
+
+  for dface in 1:ndfaces
+
+    dlabel = dface_to_label[dface]
+    if dlabel != UNSET_ID
+      continue
+    end
+
+    jfaces = dface_to_jfaces[dface]
+    for jface in jfaces
+      jlabel = jface_to_label[jface]
+      if jlabel != UNSET_ID
+        dface_to_label[dface] = jlabel
+        break
+      end
+    end
+
+  end
+
+end
+
+function _update_labels!(tface_to_label,face_to_label,tface_to_face)
+  for tface in 1:length(tface_to_label)
+    face = tface_to_face[tface]
+    if face != 0
+      tface_to_label[tface] = face_to_label[face]
+    end
+  end
+end
+
+function _generate_tcell_to_cell(ncells,nltcells)
+  ntcells = ncells * nltcells
+  tcell_to_cell = Vector{Int}(undef,ntcells)
+  tcell = 1
+  for cell in 1:ncells
+    for ltcell in 1:nltcells
+      tcell_to_cell[tcell] = cell
+      tcell += 1
+    end
+  end
+  tcell_to_cell
+end
+
+function _generate_tgrids(tgrid::Grid{D},tgraph) where D
+  tgrids = Grid[]
+  for d in 0:(D-1)
+    tface_to_vertices = connections(tgraph,d,0)
+    tgrid_d = _grid_of_dim(points(tgrid),tface_to_vertices,d)
+    push!(tgrids,tgrid_d)
+  end
+  push!(tgrids,tgrid)
+  tgrids
+end
+
+function _grid_of_dim(coords,face_to_vertices,Z)
+
+  nfaces = length(face_to_vertices)
+  fcode = tuple([TET_AXIS for i in 1:Z]...)
+
+  order = 1
+
+  _points = coords
+  _cells_data = face_to_vertices.data
+  _cells_ptrs = face_to_vertices.ptrs
+  _ctypes = ConstantCellValue(fcode,nfaces)
+  _corders = ConstantCellValue(order,nfaces)
+
+  UnstructuredGrid(
+    _points, _cells_data, _cells_ptrs, _ctypes, _corders)
+
+end
+
+function _check_order(order,co)
+  @notimplemented
+end
+
+function _check_order(order,co::ConstantCellValue)
+  @notimplementedif co.value != order
+end
+
+function _generate_ltcell_to_lpoints(ct)
+  @notimplemented
+end
+
+function _generate_ltcell_to_lpoints(ct::ConstantCellValue)
+  extrusion = ct.value
+  if extrusion == (HEX_AXIS, HEX_AXIS)
+    ltcell_to_lpoints = [[1,2,3],[4,3,2]]
+    return ltcell_to_lpoints
+  elseif extrusion == (HEX_AXIS, HEX_AXIS, HEX_AXIS)
+    ltcell_to_lpoints = [
+      [7,3,2,1], [7,5,2,1], [7,4,3,2], [7,4,8,2], [7,6,5,2], [7,6,8,2]]
+  else
+    @notimplemented
+  end
+end
+
+end # module

src/Geometry/UnstructuredGeometry.jl

@@ -55,6 +55,12 @@ function UnstructuredGrid(
     corders)
 end
 
+function UnstructuredGrid(grid::Grid)
+  @notimplemented
+end
+
+UnstructuredGrid(grid::UnstructuredGrid) = grid
+
 points(self::UnstructuredGrid) = CellValueFromArray(self.points)
 
 cells(self::UnstructuredGrid) = CellVectorFromDataAndPtrs(self.cells_data,self.cells_ptrs)