Skip to content

Commit

Permalink
Moved grid edge creation to ExtendableGrids
Browse files Browse the repository at this point in the history
  • Loading branch information
@j-fu
j-fu committed Oct 24, 2020
1 parent 96425b5 commit 1370572
Show file tree
Hide file tree
Showing 5 changed files with 17 additions and 184 deletions.
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "VoronoiFVM"
uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
authors = ["Juergen Fuhrmann <[email protected]>"]
version = "0.8.9"
version = "0.8.10"

[deps]
DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
Expand All @@ -21,7 +21,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[compat]
DiffResults = "^1.0.0"
DocStringExtensions = "^0.8.0"
ExtendableGrids = "^0.4.0"
ExtendableGrids = "^0.4.2"
ExtendableSparse = "^0.3.0"
ForwardDiff = "^0.10.0"
IterativeSolvers = "^0.8.0"
Expand Down
2 changes: 1 addition & 1 deletion src/vfvm_abstractsystem.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ function enable_species!(this::AbstractSystem,ispec::Integer, regions::AbstractV
this.region_species[ispec,ireg]=ispec
for icell=1:num_cells(this.grid)
if _cellregions[icell]==ireg
for iloc=1:size(_cellnodes,1)
for iloc=1:num_targets(_cellnodes,icell)
iglob=_cellnodes[iloc,icell]
this.node_dof[ispec,iglob]=ispec
end
Expand Down
43 changes: 14 additions & 29 deletions src/vfvm_formfactors.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ $(SIGNATURES)
Calculate node volume and voronoi surface contributions for cell.
"""

function cellfactors!(::Type{Edge1D},::Type{Cartesian1D},coord,cellnodes,icell::Int,nodefac,edgefac) where Tv
function cellfactors!(::Type{Edge1D},::Type{Cartesian1D},coord,cellnodes,icell::Int,nodefac,edgefac)
K=cellnodes[1,icell]
L=cellnodes[2,icell]
xK=coord[1,K]
Expand All @@ -78,7 +78,7 @@ function cellfactors!(::Type{Edge1D},::Type{Cartesian1D},coord,cellnodes,icell::
nothing
end

function cellfactors!(::Type{Edge1D},::Type{<:Polar1D}, coord,cellnodes,icell::Int,nodefac,edgefac) where Tv
function cellfactors!(::Type{Edge1D},::Type{<:Polar1D}, coord,cellnodes,icell::Int,nodefac,edgefac)
K=cellnodes[1,icell]
L=cellnodes[2,icell]
xK=coord[1,K]
Expand All @@ -90,11 +90,10 @@ function cellfactors!(::Type{Edge1D},::Type{<:Polar1D}, coord,cellnodes,icell::I
r1=xK[1]
end
rhalf=0.5*(r1+r0);
πv::Tv=π
# cpar[1]= πv*(r1*r1-r0*r0); # circular volume
nodefac[1]= πv*(rhalf*rhalf-r0*r0); # circular volume between midline and boundary
nodefac[2]= πv*(r1*r1-rhalf*rhalf); # circular volume between midline and boundary
edgefac[1]= 2.0*πv*rhalf/(r1-r0); # circular surface / width
# cpar[1]= π*(r1*r1-r0*r0); # circular volume
nodefac[1]= π*(rhalf*rhalf-r0*r0); # circular volume between midline and boundary
nodefac[2]= π*(r1*r1-rhalf*rhalf); # circular volume between midline and boundary
edgefac[1]= 2.0*π*rhalf/(r1-r0); # circular surface / width
nothing
end

Expand Down Expand Up @@ -218,20 +217,21 @@ $(SIGNATURES)
Calculate node volume contributions for boundary face.
"""

function bfacefactors!(::Type{Vertex0D},::Type{Cartesian1D},coord,bfacenodes,ibface::Int,nodefac) where Tv
function bfacefactors!(::Type{Vertex0D},::Type{Cartesian1D},coord,bfacenodes,ibface::Int,nodefac)
nodefac[1]=1.0
nothing
end

function bfacefactors!(::Type{Vertex0D},::Type{<:Polar1D},coord,bfacenodes,ibface::Int,nodefac) where Tv
function bfacefactors!(::Type{Vertex0D},::Type{<:Polar1D},coord,bfacenodes,ibface::Int,nodefac)
inode=bfacenodes[1,ibface]
r=coord[1,inode]
nodefac[1]=2*pi*r
nodefac[1]=2*π*r
@show nodefac[1]
nothing
end

# TODO: Test
function bfacefactors!(::Type{Edge1D},::Type{<:Cartesian2D},coord,bfacenodes,ibface::Int,nodefac) where Tv
function bfacefactors!(::Type{Edge1D},::Type{<:Cartesian2D},coord,bfacenodes,ibface::Int,nodefac)
i1=bfacenodes[1,ibface]
i2=bfacenodes[2,ibface]
dx=coord[1,i1]-coord[1,i2]
Expand All @@ -243,7 +243,7 @@ function bfacefactors!(::Type{Edge1D},::Type{<:Cartesian2D},coord,bfacenodes,ibf
end

# TODO: Test
function bfacefactors!(::Type{Edge1D},::Type{<:Cylindrical2D},coord,bfacenodes,ibface::Int,nodefac) where Tv
function bfacefactors!(::Type{Edge1D},::Type{<:Cylindrical2D},coord,bfacenodes,ibface::Int,nodefac)
i1=bfacenodes[1,ibface]
i2=bfacenodes[2,ibface]
r1=coord[1,i1]
Expand All @@ -254,27 +254,12 @@ function bfacefactors!(::Type{Edge1D},::Type{<:Cylindrical2D},coord,bfacenodes,i
rmid=(r1+r2)/2
dz=z1-z2
l=sqrt(dr*dr+dz*dz)
nodefac[1]=pi*(r1+rmid)*l/2
nodefac[2]=pi*(r2+rmid)*l/2
nodefac[1]=π*(r1+rmid)*l/2
nodefac[2]=π*(r2+rmid)*l/2
nothing
end


ExtendableGrids.num_nodes(::Type{Vertex0D})=1
num_edges(::Type{Vertex0D})=0

const cen_Edge1D=reshape([1 2],:,1)
local_celledgenodes(::Type{Edge1D})=cen_Edge1D
ExtendableGrids.num_nodes(::Type{Edge1D})=2
num_edges(::Type{Edge1D})=1

const cen_Triangle2D=[ 2 3 1; 3 1 2]
local_celledgenodes(::Type{Triangle2D})=cen_Triangle2D
ExtendableGrids.num_nodes(::Type{Triangle2D})=3
num_edges(::Type{Triangle2D})=3




##################################################################
"""
Expand Down
122 changes: 0 additions & 122 deletions src/vfvm_xgrid.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@ const Grid = ExtendableGrids.simplexgrid

num_cellregions(grid::ExtendableGrid)=grid[NumCellRegions]
num_bfaceregions(grid::ExtendableGrid)=grid[NumBFaceRegions]
num_edges(grid::ExtendableGrid)=haskey(grid,EdgeNodes) ? num_sources(grid[EdgeNodes]) : 0

cellregions(grid::ExtendableGrid)= grid[CellRegions]
bfaceregions(grid::ExtendableGrid)= grid[BFaceRegions]
Expand Down Expand Up @@ -53,124 +52,3 @@ function spherical_symmetric!(grid::ExtendableGrid)
end


"""
$(SIGNATURES)
Prepare edge adjacencies (celledges, edgecells, edgenodes)
"""
abstract type CellEdges <: AbstractGridAdjacency end
abstract type EdgeCells <: AbstractGridAdjacency end
abstract type EdgeNodes <: AbstractGridAdjacency end



function prepare_edges!(grid::ExtendableGrid)
Ti=eltype(grid[CellNodes])
cellnodes=grid[CellNodes]
geom=grid[CellGeometries][1]
# Create cell-node incidence matrix
ext_cellnode_adj=ExtendableSparseMatrix{Ti,Ti}(num_nodes(grid),num_cells(grid))
for icell=1:num_cells(grid)
for inode=1:num_nodes(geom)
ext_cellnode_adj[cellnodes[inode,icell],icell]=1
end
end
flush!(ext_cellnode_adj)
# Get SparseMatrixCSC from the ExtendableMatrix
cellnode_adj=ext_cellnode_adj.cscmatrix

# Create node-node incidence matrix for neigboring
# nodes.
nodenode_adj=cellnode_adj*transpose(cellnode_adj)

# To get unique edges, we set the lower triangular part
# including the diagonal to 0
for icol=1:length(nodenode_adj.colptr)-1
for irow=nodenode_adj.colptr[icol]:nodenode_adj.colptr[icol+1]-1
if nodenode_adj.rowval[irow]>=icol
nodenode_adj.nzval[irow]=0
end
end
end
dropzeros!(nodenode_adj)


# Now we know the number of edges and
nedges=length(nodenode_adj.nzval)


if dim_space(grid)==2
# Let us do the Euler test (assuming no holes in the domain)
v=num_nodes(grid)
e=nedges
f=num_cells(grid)+1
@assert v-e+f==2
end
if dim_space(grid)==1
@assert nedges==num_cells(grid)
end

# Calculate edge nodes and celledges
edgenodes=zeros(Ti,2,nedges)
celledges=zeros(Ti,3,num_cells(grid))
cen=local_celledgenodes(Triangle2D)

for icell=1:num_cells(grid)
for iedge=1:num_edges(geom)
n1=cellnodes[cen[1,iedge],icell]
n2=cellnodes[cen[2,iedge],icell]

# We need to look in nodenod_adj for upper triangular part entries
# therefore, we need to swap accordingly before looking
if (n1<n2)
n0=n1
n1=n2
n2=n0;
end

for irow=nodenode_adj.colptr[n1]:nodenode_adj.colptr[n1+1]-1
if nodenode_adj.rowval[irow]==n2
# If the coresponding entry has been found, set its
# value. Note that this introduces a different edge orientation
# compared to the one found locally from cell data
celledges[iedge,icell]=irow
edgenodes[1,irow]=n1
edgenodes[2,irow]=n2
end
end
end
end


# Create sparse incidence matrix for the cell-edge adjacency
ext_celledge_adj=ExtendableSparseMatrix{Ti,Ti}(nedges,num_cells(grid))
for icell=1:num_cells(grid)
for iedge=1:num_edges(geom)
ext_celledge_adj[celledges[iedge,icell],icell]=1
end
end
flush!(ext_celledge_adj)
celledge_adj=ext_celledge_adj.cscmatrix

# The edge cell matrix is the transpose
edgecell_adj=SparseMatrixCSC(transpose(celledge_adj))

# Get the adjaency array from the matrix
edgecells=zeros(Ti,2,nedges)
for icol=1:length(edgecell_adj.colptr)-1
ii=1
for irow=edgecell_adj.colptr[icol]:edgecell_adj.colptr[icol+1]-1
edgecells[ii,icol]=edgecell_adj.rowval[irow]
ii+=1
end
end

grid[EdgeCells]=edgecells
grid[CellEdges]=celledges
grid[EdgeNodes]=edgenodes
true
end

ExtendableGrids.instantiate(grid, ::Type{CellEdges})=prepare_edges!(grid) && grid[CellEdges]
ExtendableGrids.instantiate(grid, ::Type{EdgeCells})=prepare_edges!(grid) && grid[EdgeCells]
ExtendableGrids.instantiate(grid, ::Type{EdgeNodes})=prepare_edges!(grid) && grid[EdgeNodes]
30 changes: 0 additions & 30 deletions test/test_prepare_edges.jl
View file

This file was deleted.

2 comments on commit 1370572

@j-fu
Copy link
Member Author

@j-fu j-fu commented on 1370572 Oct 24, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/23593

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.8.10 -m "<description of version>" 13705729b9c11e1b603ad84b188884c4bc0db467
git push origin v0.8.10

Please sign in to comment.