use Jacobi.jl instead of C++ program for Gegenbauer nodes and weights

.gitignore

-misc/gegenbauer/*.txt
-misc/gegenbauer/gegenbauer_rule
 outputs
 pars/*/
 */**/*.so*

.gitlab-ci.yml

@@ -5,8 +5,6 @@ image: julia:1.7
 before_script:
   # workaround for https://github.com/JuliaDocs/Documenter.jl/issues/686
   - apt-get -qq update; apt-get -y install git
-  # required to build gegenbauer rule
-  - apt-get -y install g++
   - julia --project=@. -e "import Pkg; Pkg.instantiate()"
 
 unittests:

misc/gegenbauer/README

-# Gauss-Gegenbauer Rule
-
-Source: https://people.sc.fsu.edu/~jburkardt/m_src/gegenbauer_rule/gegenbauer_rule.html

src/ScalarEq/setup.jl

@@ -24,6 +24,9 @@ function Project(env::Environment, prms)
   !isnothing(av_rsolver) && register_variables!(cache, av_rsolver, dont_register=true)
   display(cache)
 
+  # setup initial data
+  initialdata!(env, P, prms["ScalarEq"])
+
   # setup boundary conditions
   bdryconds = make_BoundaryConditions(env, equation, rsolver, prms["ScalarEq"])
   ldg_bdryconds, av_bdryconds = if hrsc isa AbstractArtificialViscosity
@@ -39,9 +42,6 @@ function Project(env::Environment, prms)
   # somewhere. For periodic BCs we don't need extra storage at all.
   # But what other types of boundary conditions will be needed in the future?
 
-  # setup initial data
-  initialdata!(env, P, prms["ScalarEq"])
-
   # setup callbacks
   projectcb = make_callback(env, P, bdryconds)
 

src/dg1d.jl

@@ -56,7 +56,6 @@ export Bisection, NewtonSolver, find_root
 include("rootfinders.jl")
 include("fd.jl")
 include("lagrange.jl")
-include("gb.jl")
 include("glgl.jl")
 include("lgl.jl")
 include("evolve.jl")

src/gb.jl

-module GB
-
-
-function fname_rule(n_pts)
-  thisdir = dirname(@__FILE__)
-  fname = "$thisdir/../misc/gegenbauer/rule_$n_pts"
-  fname = normpath(fname)
-  return fname
-end # function
-
-
-
-"""
-Compute nodes and weights for Gegenbauer rule using external C++ program.
-The C++ program dg1d/misc/gegenbauer/gegenbauer_rule should be build with the commands
-```
-  \$ cd dg1d/misc/gegenbauer
-  \$ g++ gegenbauer_rule.cpp -o gegenbauer_rule
-```
-
-n_pts ... = order + 1
-"""
-function generate_nodes_weights(n_pts)
-
-  # only interested in standard interval x ϵ [-1,1] with weight (1-t^2)^2
-  α = 2.0
-  A = -1.0
-  B = 1.0
-
-  thisdir = dirname(@__FILE__)
-  exe = "$thisdir/../misc/gegenbauer/gegenbauer_rule"
-  exe = normpath(exe)
-
-  if ! isfile(exe)
-    src = normpath("$thisdir/../misc/gegenbauer/gegenbauer_rule.cpp")
-    run(`g++ $src -o $exe`)
-  end
-
-  fname = fname_rule(n_pts)
-
-  pipe = pipeline(`$exe $n_pts $α $A $B $fname`, stdout=devnull, stderr=devnull)
-  run(pipe)
-end # function
-
-
-
-"""
-Read Gauss-Gegenbauer quadrature nodes and weights from file located at
-evm/gegenbauer. If no file is found for requested number of points n_pts,
-we generate it by calling rule(n_pts).
-n_pts must be greater than 0.
-"""
-function rule(n_pts)
-
-  fname = fname_rule(n_pts)
-  fname_w = "$(fname)_w.txt"
-  fname_r = "$(fname)_r.txt"
-  fname_x = "$(fname)_x.txt"
-
-  rule_exists = isfile(fname_w) && isfile(fname_r) && isfile(fname_x)
-
-  if ! rule_exists 
-    generate_nodes_weights(n_pts)
-  end
-
-  x = zeros(n_pts)
-  w = zeros(n_pts)
-
-  file = open(fname_x, "r")
-  idx = 1
-  while !eof(file)
-    l = readline(file)
-    l = strip(l)
-    x[idx] = parse(Float64, l)
-    idx += 1
-  end
-  close(file)
-
-  file = open(fname_w, "r")
-  idx = 1
-  while !eof(file)
-    l = readline(file)
-    l = strip(l)
-    w[idx] = parse(Float64, l)
-    idx += 1
-  end
-  close(file)
-
-  return x, w
-end # function
-
-
-end # GB

src/glgl.jl

@@ -2,11 +2,11 @@ module GLGL
 
 
 using LinearAlgebra, Jacobi
-import dg1d: GB, purge_zeros!, lagrange_poly
+import dg1d: purge_zeros!, lagrange_poly
 
 
 @doc """
-Computes quadrature nodes, weights and derivative weights 
+Computes quadrature nodes, weights and derivative weights
 for the Generalized-Gauss-Lobatto quadrature rule.
 n_pts ... number of quadrature nodes
 Return values
@@ -19,7 +19,10 @@ function rule(n_pts)
   N = n_pts - 2
   gb_x, gb_w = [], []
   if N > 0
-    gb_x, gb_w = GB.rule(N)
+    # Integration nodes and weights for Gauss-Gegenbauer quadrature
+    α = 2; β = 2
+    gb_x = zgj(N,    α, β)
+    gb_w = wgj(gb_x, α, β)
   end
 
   x = zeros(n_pts)