From 13705729b9c11e1b603ad84b188884c4bc0db467 Mon Sep 17 00:00:00 2001 From: Juergen Fuhrmann Date: Sat, 24 Oct 2020 14:42:24 +0200 Subject: [PATCH] Moved grid edge creation to ExtendableGrids --- Project.toml | 4 +- src/vfvm_abstractsystem.jl | 2 +- src/vfvm_formfactors.jl | 43 +++++-------- src/vfvm_xgrid.jl | 122 ------------------------------------- test/test_prepare_edges.jl | 30 --------- 5 files changed, 17 insertions(+), 184 deletions(-) delete mode 100644 test/test_prepare_edges.jl diff --git a/Project.toml b/Project.toml index 87d5e630f..b86a53e86 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "VoronoiFVM" uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3" authors = ["Juergen Fuhrmann "] -version = "0.8.9" +version = "0.8.10" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" @@ -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" diff --git a/src/vfvm_abstractsystem.jl b/src/vfvm_abstractsystem.jl index 8dff0d6eb..6df0036f9 100644 --- a/src/vfvm_abstractsystem.jl +++ b/src/vfvm_abstractsystem.jl @@ -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 diff --git a/src/vfvm_formfactors.jl b/src/vfvm_formfactors.jl index fc84f6040..5369884d4 100644 --- a/src/vfvm_formfactors.jl +++ b/src/vfvm_formfactors.jl @@ -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] @@ -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] @@ -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 @@ -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] @@ -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] @@ -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 - - - ################################################################## """ diff --git a/src/vfvm_xgrid.jl b/src/vfvm_xgrid.jl index 2bfb42360..ba0324e93 100644 --- a/src/vfvm_xgrid.jl +++ b/src/vfvm_xgrid.jl @@ -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] @@ -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