added docstrings; adjusted example elements to P2 Bubble

docs/make.jl

@@ -1,6 +1,8 @@
 using Documenter, ExampleJuggler, PlutoStaticHTML, VoronoiFVM, DocumenterCitations
 using ExtendableGrids, GridVisualize, LinearAlgebra, RecursiveArrayTools, SciMLBase
 
+using ExtendableFEMBase, ExtendableFEM
+
 using OrdinaryDiffEqBDF, OrdinaryDiffEqLowOrderRK, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqSDIRK, OrdinaryDiffEqTsit5 
 
 function make(; with_examples = true,
@@ -33,7 +35,8 @@ function make(; with_examples = true,
             "misc.md",
             "internal.md",
             "allindex.md",
-            "devel.md",]
+            "devel.md",
+            "extensions.md"]
     ]
 
 

docs/src/extensions.md

@@ -0,0 +1,8 @@
+# [`ExtendableFEMBase`](https://github.com/chmerdon/ExtendableFEMBase.jl/) Extension
+
+sometextthatexplainssomething
+
+```@docs
+VoronoiFVM.edgevelocities(grid, vel::FEVectorBlock; kwargs...)
+VoronoiFVM.bfacevelocities(::Any,::FEVectorBlock)
+```

examples/Example240_FiniteElementVelocities.jl

@@ -145,6 +145,7 @@ function runtests()
     end
 end
 
+# test analytical distance and calc_divergences
 function full_fem_demo(; coord_system = Cartesian2D, nref = 1, Plotter = nothing,
         μ = 1.0e-02, D = 0.01, cin = 1.0, assembly = :edgewise, interpolation_eps = 1.0e-09)
     H = 1.0
@@ -197,8 +198,8 @@ function compute_velocity(flowgrid, coord_system, μ = 1.0e-02; interpolation_ep
     axisymmetric = coord_system == Cylindrical2D ? true : false
 
     # define finite element spaces
-    FE_v, FE_p = H1BR{2}, L2P0{1}
-    reconst_FEType = HDIVBDM1{2}
+    FE_v, FE_p = H1P2B{2,2}, L2P1{1}
+    reconst_FEType = HDIVBDM2{2}
     FES = [FESpace{FE_v}(flowgrid), FESpace{FE_p}(flowgrid; broken = true)]
 
     # describe problem

ext/VoronoiFVMExtendableFEMBaseExt.jl

@@ -13,6 +13,9 @@ using LinearAlgebra: dot, norm
 
 using Base: fill!
 
+using DocStringExtensions: DocStringExtensions, SIGNATURES, TYPEDEF,
+                           TYPEDFIELDS, TYPEDSIGNATURES
+
 id(u) = (u, Identity)
 
 function iscloser(pint, p1, p2, eps)
@@ -66,7 +69,15 @@ function prepare_segment_integration(vel; axisymmetric=false, reconst=false, kwa
     return seg_integrator, point_evaluator, cellfinder, flowgrid
 end
 
+"""
+$(SIGNATURES)
+
+Compute [`VoronoiFVM.edgevelocities`](@ref) for a finite element flow field computed 
+by [`ExtendableFEM`](https://github.com/chmerdon/ExtendableFEM.jl).
+"""
 function VoronoiFVM.edgevelocities(grid, vel::FEVectorBlock; kwargs...)
+    # construct an augmented type to gather precomputed information 
+    # in order to pass it to the repeated integrate call in VoronoiFVM.edgevelocities
     axisymmetric = grid[CoordinateSystem] <: Cylindrical2D ? true : false
     seg_integrator, point_evaluator, cf, flowgrid = prepare_segment_integration(vel; axisymmetric, kwargs...)
     aug_fevec_block = AugmentedFEVectorBlock(vel, seg_integrator, point_evaluator, cf, flowgrid)
@@ -79,6 +90,12 @@ function VoronoiFVM.edgevelocities(grid, vel::FEVectorBlock; kwargs...)
     return velovec
 end
 
+"""
+$(SIGNATURES)
+
+Compute [`VoronoiFVM.bfacevelocities`](@ref) for a finite element flow field computed 
+by [`ExtendableFEM`](https://github.com/chmerdon/ExtendableFEM.jl).
+"""
 function VoronoiFVM.bfacevelocities(grid, vel::FEVectorBlock; kwargs...)
     axisymmetric = grid[CoordinateSystem] <: Cylindrical2D ? true : false
     seg_integrator, point_evaluator, cf, flowgrid = prepare_segment_integration(vel; axisymmetric, kwargs...)
@@ -93,14 +110,20 @@ function VoronoiFVM.bfacevelocities(grid, vel::FEVectorBlock; kwargs...)
     return velovec
 end
 
+# We need two explicitly type-annotated methods for a working method specialization.
+# This one...
 function VoronoiFVM.integrate(::Type{<:Cartesian2D}, p1, p2, hnormal, aug_vec_block::AugmentedFEVectorBlock; kwargs...)
     _integrate_along_segments(p1, p2, hnormal, aug_vec_block; kwargs...)
 end
 
+# ... and that one.
 function VoronoiFVM.integrate(::Type{<:Cylindrical2D}, p1, p2, hnormal, aug_vec_block::AugmentedFEVectorBlock; kwargs...)
     _integrate_along_segments(p1, p2, hnormal, aug_vec_block; kwargs...)
 end
 
+# compute the path integral for the velocity in aug_vec_block between p1 and p2 by 
+# incrementally walking through each cell in the grid between p1 and p2
+# and summing up each cell's contribution
 function _integrate_along_segments(p1, p2, hnormal, aug_vec_block::AugmentedFEVectorBlock; interpolate_eps=1.0e-12, axisymmetric=false, kwargs...)
 
     edge_length = norm(p1 - p2, 2)

src/vfvm_formfactors.jl

@@ -326,6 +326,14 @@ end
 #
 # TODO: this should be generalized for more quadrules
 #
+"""
+$(SIGNATURES)
+
+Integrate `velofunc` along the edge ``\\sigma=\\overline{\\mathtt{coordl}\\,\\mathtt{coordr}}`` 
+between the ``x_K`` and ``x_L`` where ``\\mathtt{hnormal}=x_K-x_L`` using [Simpson's Rule](https://en.wikipedia.org/wiki/Simpson%27s_rule).
+To be precise, compute for a cartesian coordinate system:
+``\\int_{\\sigma} \\mathbf{v} \\cdot \\mathbf{n} \\,\\mathrm{ds} \\lvert x_K - x_L \\rvert / \\lvert\\sigma\\rvert``.
+"""
 function integrate(::Type{<:Cartesian2D}, coordl, coordr, hnormal, velofunc; kwargs...)
     wl = 1.0 / 6.0
     wm = 2.0 / 3.0
@@ -337,6 +345,13 @@ function integrate(::Type{<:Cartesian2D}, coordl, coordr, hnormal, velofunc; kwa
     return (wl * vxl + wm * vxm + wr * vxr) * hnormal[1] + (wl * vyl + wm * vym + wr * vyr) * hnormal[2]
 end
 
+"""
+$(SIGNATURES)
+
+Similar to `integrate(::Type{<:Cartesian2D},...)`, but computes instead
+``\\int_{\\sigma} r \\, \\mathbf{v} \\cdot \\mathbf{n} \\,\\mathrm{ds} \\lvert x_K - x_L \\rvert / \\left ( \\lvert\\sigma\\rvert r(\\mathrm{mid}(\\sigma)) \\right )``
+where ``r(\\mathrm{mid}(\\sigma))`` is the ``r``-coordinate of the mid-point of ``\\sigma``.
+"""
 function integrate(::Type{<:Cylindrical2D}, coordl, coordr, hnormal, velofunc; kwargs...)
     wl = 1.0 / 6.0
     wm = 2.0 / 3.0
@@ -365,7 +380,9 @@ end
 """
 $(SIGNATURES)
 
-Project velocity onto grid edges,
+Project velocity onto grid edges.
+That is, we compute the path integrals of the given `velofunc` along the 
+Voronoi cell edges as provided by [`integrate`](@ref).
 """
 function edgevelocities(grid, velofunc; kwargs...)
     @assert dim_space(grid) < 3
@@ -413,7 +430,7 @@ end
 """
 $(SIGNATURES)
 
-Project velocity onto boundary face normals
+Similar to [`edgevelocities`](@ref), but for boundary faces.
 """
 function bfacevelocities(grid, velofunc; kwargs...)
     @assert dim_space(grid) < 3