From b9df47247151260ce8941562b53c040d779e8d2f Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Thu, 6 Feb 2014 01:36:45 -0500
Subject: [PATCH] Add the (in progress) FEM.R1 module.

---
 numerical-analysis.cabal |   1 +
 src/FEM/R1.hs            | 169 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 170 insertions(+)
 create mode 100644 src/FEM/R1.hs

diff --git a/numerical-analysis.cabal b/numerical-analysis.cabal
index 92b716e..b666583 100644
--- a/numerical-analysis.cabal
+++ b/numerical-analysis.cabal
@@ -18,6 +18,7 @@ data-files: makefile
 
 library
   exposed-modules:
+    FEM.R1
     Integration.Gaussian,
     Integration.Simpson,
     Integration.Trapezoid,
diff --git a/src/FEM/R1.hs b/src/FEM/R1.hs
new file mode 100644
index 0000000..c6a76c8
--- /dev/null
+++ b/src/FEM/R1.hs
@@ -0,0 +1,169 @@
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- | Finite Element Method (FEM) to solve the problem,
+--
+--     -(A(x)*u'(x))' + c(x)*u*(x) = f(x) on (a,b)
+--
+--   with one of the two types of boundary conditions,
+--
+--     * Dirichlet: u(a) = u(b) = 0
+--
+--     * A(a)u'(a) = alpha, A(b)u'(b) = beta
+--
+--   if c(x) = 0, then it is assumed that the boundary conditions are
+--   Dirichlet.
+--
+module FEM.R1
+where
+
+import qualified Algebra.Algebraic as Algebraic ( C )
+import qualified Algebra.Field as Field ( C )
+import qualified Algebra.RealField as RealField ( C )
+import Data.Vector.Fixed ( Arity, S )
+import NumericPrelude
+import qualified Prelude as P
+
+import Linear.Matrix (
+  Col,
+  Mat(..),
+  (!!!),
+  construct,
+  nrows )
+import Polynomials.Orthogonal ( legendre )
+
+-- | Dirichlet boundary conditions. Since u(a)=u(b)=0 are fixed,
+--   there's no additional information conveyed by this type.
+data Dirichlet = Dirichlet
+
+-- | Neumann boundary conditions. @alpha@ specifies A(a)u'(b) and
+--   @beta@ specifies A(b)u'(b).
+data Neumann a = Neumann { alpha :: a, beta :: a }
+
+-- | Boundary conditions can be either Dirichlet or Neumann.
+type BoundaryConditions a = Either Dirichlet (Neumann a)
+
+type Interval a = (a,a)
+
+-- | The data for the PDE that we are attempting to solve.
+data PDE a =
+  PDE {
+    -- | A(x)
+    big_A   :: (a -> a),
+    -- | c(x)
+    c      :: (a -> a),
+    -- | f(x)
+    f      :: (a -> a),
+    -- | The domain in R^1 as an interval
+    domain :: Interval a,
+    bdy    :: BoundaryConditions a }
+
+
+
+-- | Non-PDE parameters for the finite element method. The additional
+--   type parameter @n@ should be a type-level representation of the
+--   largest element in @max_degrees@. It needs to be known statically
+--   for the dimensions of the pointer matrix.
+data Params m n a =
+  Params {
+    -- | A partition of the domain.
+    mesh :: Col m (Interval a),
+
+    -- | A list of the highest-degree polynomial that we will use over
+    --   each interval in our mesh.
+    max_degrees :: Col m Int }
+
+
+-- | The pointer matrix mapping local basis elements to global
+--   ones. The (i,j)th entry of the matrix contains the index of the
+--   global basis function corresponding to the jth local basis
+--   function over element i.
+--
+--   TODO: This works for Dirichlet boundary conditions only.
+--
+--   Note that the number of columns in the matrix depends on the
+--
+--   Examples:
+--
+--   >>> import Data.Vector.Fixed ( N5, N6 )
+--   >>> import Linear.Matrix ( Col5, fromList )
+--
+--   >>> let p = fromList [[3],[3],[5],[4],[5]] :: Col5 Int
+--   >>> let mesh = undefined :: Col5 (Int,Int)
+--   >>> let params = Params mesh p :: Params N5 N5 Int
+--   >>> let row1 = [0,1,5,6,0,0] :: [Int]
+--   >>> let row2 = [1,2,7,8,0,0] :: [Int]
+--   >>> let row3 = [2,3,9,10,11,12] :: [Int]
+--   >>> let row4 = [3,4,13,14,15,0] :: [Int]
+--   >>> let row5 = [4,0,16,17,18,19] :: [Int]
+--   >>> let expected = fromList [row1,row2,row3,row4,row5] :: Mat N5 N6 Int
+--   >>> pointer params == expected
+--   True
+--
+pointer :: (Arity m, Arity n) => Params m n a -> Mat m (S n) Int
+pointer params =
+  construct lambda
+  where
+    -- | The number of polynomials in the local basis for element i.
+    numpolys :: Int -> Int
+    numpolys i = ((max_degrees params) !!! (i,0)) + 1
+
+    lambda i j
+
+      -- The first two (linear) basis functions are easy to do.
+      | j < 2 = if i + j >= (nrows $ mesh params) then 0 else i + j
+
+      -- No local basis function exists for this j.
+      | j >= (numpolys i) = 0
+
+      -- The first row we handle as a special case.
+      | i == 0 = (nrows $ mesh params) + j - 2
+
+      -- The first entry for this row not corresponding to one of the
+      -- linear polynomials (which we did in the first guard). This
+      -- grabs the biggest number in the previous row and begins counting
+      -- from there.
+      | j == 2 = (lambda (i-1) (numpolys (i-1) - 1)) + 1
+
+      -- If none of the other cases apply, we can compute our value by
+      -- looking to the left in the matrix.
+      | otherwise = (lambda i (j-1)) + 1
+
+
+
+-- | An affine transform taking the interval [x1,x2] into [-1,1].
+--
+--   Examples:
+--
+--   >>> let phi = affine (-6,9)
+--   >>> phi (-6)
+--   -1.0
+--   >>>  phi (9)
+--   1.0
+--
+affine :: Field.C a => (a,a) -> (a -> a)
+affine (x1,x2) x = (fromInteger 2)*(x - x1)/(x2 - x1) - (fromInteger 1)
+
+
+-- | Orthonormal basis functions over [-1,1]. The test case below
+--   comes from Sage where the orthogonality of the polynomials'
+--   derivatives can easily be tested.
+--
+--   Examples:
+--
+--   >>> import qualified Algebra.Absolute as Absolute ( abs )
+--
+--   >>> let expected = 6.33910180790284
+--   >>> let actual = big_N 3 1.5 :: Double
+--   >>> Absolute.abs (actual - expected) < 1e-12
+--   True
+--
+big_N :: forall a. (Algebraic.C a, RealField.C a) => Integer -> a -> a
+big_N k x
+  | k < 0 = error "requested a negative basis function"
+  | otherwise =
+      coeff * ( legendre (k+2) x - legendre k x )
+      where
+        four = fromInteger 4
+        six = fromInteger 6
+        coeff = one / (sqrt (four*(fromInteger k) + six)) :: a
-- 
2.43.2