diff --git a/Main.hs b/Main.hs deleted file mode 100644 --- a/Main.hs +++ /dev/null @@ -1,8 +0,0 @@ -{-# OPTIONS_GHC -Wall #-} - -module Main where - -import Numeric.HFoil.Repl - -main :: IO () -main = run diff --git a/Numeric/HFoil.hs b/Numeric/HFoil.hs deleted file mode 100644 --- a/Numeric/HFoil.hs +++ /dev/null @@ -1,26 +0,0 @@ -{- | - Module : Numeric.Hfoil - Description : Top level module - - This is the top level module which exports the API - -} - -{-# OPTIONS_GHC -Wall #-} - -module Numeric.HFoil( -- * Airfoils - module Numeric.HFoil.Foil - -- * Flow solution - , module Numeric.HFoil.Flow - -- * Naca 4 utilities - , module Numeric.HFoil.Naca4 - -- * Drawing utilities (gloss backend) - , module Numeric.HFoil.Drawing - -- * Interactive command line application - , module Numeric.HFoil.Repl - ) where - -import Numeric.HFoil.Naca4 -import Numeric.HFoil.Foil -import Numeric.HFoil.Drawing -import Numeric.HFoil.Flow -import Numeric.HFoil.Repl diff --git a/Numeric/HFoil/Drawing.hs b/Numeric/HFoil/Drawing.hs deleted file mode 100644 --- a/Numeric/HFoil/Drawing.hs +++ /dev/null @@ -1,168 +0,0 @@ -{-# OPTIONS_GHC -Wall #-} -{-# LANGUAGE FlexibleContexts #-} - -module Numeric.HFoil.Drawing( drawLine - , drawLineV - , drawSolution - , drawFoil - , drawOnce - , drawNormals - , drawForces - , drawKuttas - ) where - -import Graphics.Gloss hiding(Vector,dim) -import Numeric.LinearAlgebra hiding(Element, scale,i) -import Foreign.Storable(Storable) -import qualified Numeric.LinearAlgebra as LA -import Text.Printf - -import Numeric.HFoil.Flow -import Numeric.HFoil.Foil - -xSize, ySize :: Int -xSize = 800 -ySize = 500 - -cpScale :: Double -cpScale = -0.25 - -border :: Float -border = 0.7 - -normalLengths :: Double -normalLengths = 0.01 - -drawText :: Color -> (Float, Float) -> Float -> String -> Picture -drawText col (x,y) size str = translate (0.5*x*(fromIntegral xSize)) (0.5*y*(fromIntegral ySize)) - $ scale size size - $ color col - $ Text str - -drawLine :: Real a => Color -> [(a,a)] -> Picture -drawLine col coords = scale (border*(fromIntegral xSize)) (border*(fromIntegral xSize)) - $ translate (-0.5) 0 - $ color col - $ line $ map (\(x,y) -> (realToFrac x, realToFrac y)) coords - -drawCircle :: Real a => Color -> (a, a) -> Float -> Picture -drawCircle col (x,y) size = scale (border*(fromIntegral xSize)) (border*(fromIntegral xSize)) - $ translate (-0.5 + realToFrac x) (realToFrac y) - $ color col - $ Circle size - -drawLineV :: (Real a, Storable a) => Color -> (Vector a, Vector a) -> Picture -drawLineV col (vx, vy) = drawLine col $ zip (toList vx) (toList vy) - -drawFoil :: (Real a, Storable a, Show a) => Foil a -> Picture -drawFoil (Foil elements _) = pictures $ map drawElement elements - -drawElement :: (Real a, Storable a) => Element a -> Picture -drawElement element = drawLineV white (fNodes element) - -drawNormals :: Foil Double -> Picture -drawNormals (Foil elements _) = pictures $ map (\(xy0, xy1) -> drawLine green [xy0, xy1]) (zip xy0s xy1s) - where - xy0s = zip (toList xm) (toList ym) - xy1s = zip (toList (xm + (LA.scale normalLengths xUnitNormal))) (toList (ym + (LA.scale normalLengths yUnitNormal))) - - (xUnitNormal, yUnitNormal) = (\(x,y) -> (join x, join y)) $ unzip $ map fUnitNormals elements - (xm, ym) = (\(x,y) -> (join x, join y)) $ unzip $ map fMidpoints elements - -colorFun :: (Fractional a, Real a) => a -> a -> a -> Color -colorFun min' max' x' = makeColor (1-x) (1-x) x 1 - where - x = realToFrac $ (x' - min')/(max'-min') - -drawKuttas :: (Real a, Storable a) => FlowSol a -> Picture -drawKuttas flow = pictures $ concatMap (\(k0,k1) -> [circ k0, circ k1]) kis - where - kis = solKuttaIndices flow - (xs',ys') = unzip $ map fMidpoints $ (\(Foil els _) -> els) (solFoil flow) - xs = join xs' - ys = join ys' - circ k = drawCircle yellow (xs @> k, ys @> k) 0.006 - -drawForces :: FlowSol Double -> Picture -drawForces flow = pictures $ map (\(xy0, xy1, cp) -> drawLine (colorFun minCp maxCp cp) [xy0, xy1]) - $ zip3 xy0s xy1s (toList (solCps flow)) - where - xy0s = zip (toList xm) (toList ym) - xy1s = zip (toList (xm + xPressures)) (toList (ym + yPressures)) - (xPressures, yPressures) = (\(x,y) -> (LA.scale c x/lengths, LA.scale c y/lengths)) (solForces flow) - lengths = join $ map fLengths $ (\(Foil x _) -> x) $ solFoil flow - (xm, ym) = (\(x,y) -> (join x, join y)) $ unzip $ map fMidpoints $ (\(Foil x _) -> x) $ solFoil flow - - c = 0.1 - - maxCp = maxElement (solCps flow) - minCp = minElement (solCps flow) - -drawColoredFoil :: [Color] -> Foil Double -> Picture -drawColoredFoil colors foil@(Foil elements _) = pictures $ zipWith drawColoredElement colors' elements - where - colors' = groupSomethingByFoil foil colors - -drawColoredElement :: [Color] -> Element Double -> Picture -drawColoredElement colors element = pictures $ map (\(xy0, xy1, col) -> drawLine col [xy0, xy1]) (zip3 xy0s xy1s colors) - where - xys = (\(x,y) -> zip (toList x) (toList y)) $ fNodes element - xy0s = tail xys - xy1s = init xys - -groupSomethingByFoil :: Storable a => Foil a -> [b] -> [[b]] -groupSomethingByFoil (Foil elements _) somethings = f somethings (map (dim . fAngles) elements) - where - f xs (n:ns) = (take n xs):(f (drop n xs) ns) - f [] []= [] - f _ _ = error "uh oh (groupSomethingByFoil)" - -drawSolution :: FlowSol Double -> Picture -drawSolution flow = pictures $ onscreenText ++ - [ drawColoredFoil colors foil - , drawCircle white (fst $ solCenterPressure flow, snd $ solCenterPressure flow) 0.006 - , drawCircle white (fst $ solCenterPressure flow, 0) 0.006 - , drawCircle green (0.25,0::Double) 0.006 - ] ++ zipWith (\x y -> drawLineV red (x, y)) xs - (takesV (map dim xs) (LA.scale cpScale cps)) -- cp graph - where - foil@(Foil elements name) = solFoil flow - cps = solCps flow - - xs = map (fst . fMidpoints) elements - - onscreenText = zipWith (\m y -> drawText white (0.45,y) 0.15 m) msgs - $ take (length msgs) [0.8,0.65..] - msgs = [ name - , printf ("alpha: %.6f deg") ((solAlpha flow)*180/pi) - , printf ("Cl: %.6f") (solCl flow) - , printf ("Cd: %.6f") (solCd flow) - , printf ("Cm: %.6f (c/4, 0)") (solCm flow) - ] - - colors = map (colorFun (minElement cps) (maxElement cps)) (toList cps) - - - -drawOnce :: [Picture] -> IO () -drawOnce pics = do - display - (InWindow - "hfoil" -- window title - (xSize, ySize) -- window size - (10, 650)) -- window position - black -- background color - (pictures pics) -- picture to display - --- let line = plot_lines_values ^= [[ (xc, yt (naca4 "0012") xc) --- | xc <- [0,0.01..0.99::Double]]] --- $ plot_lines_title ^= "naca 0012" --- $ defaultPlotLines --- --- chart = layout1_title ^= "naca yo" --- $ layout1_plots ^= [Left (toPlot line)] --- $ defaultLayout1 --- --- renderableToWindow (toRenderable chart) 640 480 --- _ <- renderableToPNGFile (toRenderable chart) 640 480 "mDiv_vs_tc.png" --- return () diff --git a/Numeric/HFoil/Flow.hs b/Numeric/HFoil/Flow.hs deleted file mode 100644 --- a/Numeric/HFoil/Flow.hs +++ /dev/null @@ -1,206 +0,0 @@ -{-# OPTIONS_GHC -Wall #-} -{-# LANGUAGE FlexibleContexts #-} - -module Numeric.HFoil.Flow( FlowSol(..) - , solveFlow - ) where - -import Data.Packed.ST(newMatrix, writeMatrix, freezeMatrix) -import Control.Monad(forM_) -import Control.Monad.ST(runST) -import Numeric.LinearAlgebra hiding(i) -import Foreign.Storable - -import Numeric.HFoil.Foil - -data FlowSol a = FlowSol { solFoil :: Foil a - , solVs :: Vector a - , solCps :: Vector a - , solVorticities :: [a] - , solAlpha :: a - , solForces :: (Vector a, Vector a) - , solForce :: (a,a) - , solCl :: a - , solCd :: a - , solCm :: a -- moment about (0,0.25) - , solCenterPressure :: (a,a) - , solKuttaIndices :: [(Int,Int)] - } - -solveFlow :: (Num (Vector a), RealFloat a, Field a) => Foil a -> a -> FlowSol a -solveFlow foil@(Foil elements _) alpha = - FlowSol { solFoil = foil - , solVs = vs - , solCps = cps - , solVorticities = vorticities - , solAlpha = alpha - , solForces = (xForces, yForces) - , solForce = (xf,yf) - , solCl = cl - , solCd = cd - , solCm = cm - , solCenterPressure = (xCp, yCp) - , solKuttaIndices = kuttaIndices - } - where - -- surface velocities - (vs,vorticities) = ( (mapVector (\q -> cos(q - alpha)) angles) + (mV <> qsGamma) - , map (qsGamma @>) [n-(length elements),n-1] - ) - where - (mA, mV, b) = getAVb geometries angles kuttaIndices alpha - qsGamma = flatten $ linearSolve mA b - n = dim angles - -- pressure coefficients - cps = 1 - vs*vs - - -- forces and force coefficients - xForces = -cps*(join $ map (fst . fNormals) elements) - yForces = -cps*(join $ map (snd . fNormals) elements) - xf = sumElements xForces - yf = sumElements yForces - - cd = xf*(cos alpha) + yf*(sin alpha) - cl = -xf*(sin alpha) + yf*(cos alpha) - - -- midpoints - (xs,ys) = (\(x,y) -> (join x, join y)) $ unzip $ map fMidpoints elements - - -- centers of pressure - (xCp, yCp) = ( (sumElements (xs*yForces)) / (sumElements yForces) - , (sumElements (ys*xForces)) / (sumElements xForces) - ) - - -- moment about cp and quarter chord - cm = sumElements $ (mapVector (\z -> z - 0.25) xs)*yForces - ys*xForces - - kuttaIndices = ki 0 (map dim angles') - where - ki _ [] = [] - ki n0 (n:ns) = (n0,n0+n-1):(ki (n0+n) ns) - - angles = join angles' - angles' = map fAngles elements - - -- (mSL,mCL,mSB,mCB) - geometries = setupGeometries angles (xms, yms) (xns, yns) (xnps, ynps) - where - joinTuples (x,y) = (join x, join y) - (xms,yms) = joinTuples $ unzip $ map fMidpoints elements - (xns,yns) = joinTuples $ unzip $ map fInits elements - (xnps,ynps) = joinTuples $ unzip $ map fTails elements - - --- make influence matrix A and also the matrix V where V*[sources; vortex] == tangential speeds ---getAVb :: (Num (Vector a), RealFloat a, Container Vector a) => --- Foil a -> a -> (Matrix a, Matrix a) -getAVb :: (Floating a, Num (Vector a), Container Vector a, Storable a) => - (Matrix a, Matrix a, Matrix a, Matrix a) - -> Vector a - -> [(Int, Int)] - -> a - -> (Matrix a, Matrix a, Matrix a) -getAVb (mSL,mCL,mSB,mCB) angles kuttaIndices alpha = - ( scale (1/(2*pi)) $ fromBlocks [[ mAij, fromColumns vsAin] - ,[fromRows vsAnj, fromLists ann]] - , scale (1/(2*pi)) $ fromBlocks [[ mSB - mCL - , fromColumns vsTin]] - , asColumn $ join $ (mapVector (\q -> sin $ q - alpha) angles): - (map (\(n0,nf) -> fromList [-cos((angles @> n0) - alpha) - cos((angles @> nf) - alpha)]) kuttaIndices) - ) -- (for middle entry mV, mSL + mCB == mAij) - where - n = rows mSL - - -- sources influence matrix Aij - mAij = mSL + mCB - - -- vortices influence vectors Ains - vsAin = map (\(n0,nf) -> fromList $ map sumElements $ toRows $ subMatrix (0,n0) (n, nf-n0+1) m) kuttaIndices - where - m = (mCL - mSB) - - -- vortices tangential velocity outputs - vsTin = map (\(n0,nf) -> fromList $ map sumElements $ toRows $ subMatrix (0,n0) (n, nf-n0+1) m) kuttaIndices - where - m = (mSL + mCB) - - -- sources kutta condition influence vector Anj - vsAnj = map (\(n0,nf) -> (getRow n0 mSB) - (getRow n0 mCL) + (getRow nf mSB) - (getRow nf mCL)) kuttaIndices - - -- vortices kutta condition influence scalars ann - ann = (flip map) kuttaIndices $ \(ni0,nif) -> (flip map) kuttaIndices $ \(nj0,njf) -> - sumElements $ (getSubRow ni0 (nj0,njf) mSL) + (getSubRow ni0 (nj0,njf) mCB) + (getSubRow nif (nj0,njf) mSL) + (getSubRow nif (nj0,njf) mCB) - - getRow i mat = flatten $ subMatrix (i,0) (1,n) mat - getSubRow i (j0,jf) mat = flatten $ subMatrix (i,j0) (1,jf-j0+1) mat - - -{- -calcuate 4 matrices which will be useful - -mSL(i,j) = sin(qi - ij) * ln(rij+1/rij) -mCL(i,j) = cos(qi - ij) * ln(rij+1/rij) -mSB(i,j) = sin(qi - ij) * beta(i,j) -mCB(i,j) = cos(qi - ij) * beta(i,j) - -where for i==j: ln(r/r) = 0 - beta = pi -are explicitly set --} -setupGeometries :: (RealFloat t, Storable t) => - Vector t - -> (Vector t, Vector t) - -> (Vector t, Vector t) - -> (Vector t, Vector t) - -> (Matrix t, Matrix t, Matrix t, Matrix t) -setupGeometries angles (xms, yms) (xns, yns) (xnps, ynps) = runST $ do - let n = dim angles - mSL' <- newMatrix 0 n n - mCL' <- newMatrix 0 n n - mSB' <- newMatrix 0 n n - mCB' <- newMatrix 0 n n - - forM_ [0..n-1] $ \i -> forM_ [0..n-1] $ \j -> do - let qi = angles @> i - qj = angles @> j - - xmi = xms @> i - ymi = yms @> i - xnj = xns @> j - ynj = yns @> j - xnjp = xnps @> j - ynjp = ynps @> j - - s = sin (qi - qj) - c = cos (qi - qj) - - l - | i == j = 0 - | otherwise = log(r1/r0) - where - r1 = distance (xmi, ymi) (xnjp, ynjp) - r0 = distance (xmi, ymi) (xnj, ynj ) - distance (x1,y1) (x0,y0) = sqrt $ dx*dx + dy*dy - where - dx = x1 - x0 - dy = y1 - y0 - - b - | i == j = pi - | otherwise = atan2 (dyjp*dxj - dxjp*dyj) - (dxjp*dxj + dyjp*dyj) - where - dyj = ymi - ynj - dxj = xmi - xnj - dyjp = ymi - ynjp - dxjp = xmi - xnjp - - _ <- writeMatrix mSL' i j (s*l) - _ <- writeMatrix mCL' i j (c*l) - _ <- writeMatrix mSB' i j (s*b) - writeMatrix mCB' i j (c*b) - sl' <- freezeMatrix mSL' - cl' <- freezeMatrix mCL' - sb' <- freezeMatrix mSB' - cb' <- freezeMatrix mCB' - return (sl',cl',sb',cb') diff --git a/Numeric/HFoil/Foil.hs b/Numeric/HFoil/Foil.hs deleted file mode 100644 --- a/Numeric/HFoil/Foil.hs +++ /dev/null @@ -1,157 +0,0 @@ -{-# OPTIONS_GHC -Wall #-} -{-# LANGUAGE FlexibleContexts #-} - -module Numeric.HFoil.Foil( Foil(..) - , Element(..) - , panelizeNaca4 - , loadFoil - , getUIUCFoil - ) where - -import System.Directory(doesFileExist) -import Numeric.LinearAlgebra hiding (Element) -import Data.Tuple.Utils(fst3) -import Network.HTTP(simpleHTTP, getRequest, getResponseBody) -import Foreign.Storable(Storable) - -import qualified Numeric.HFoil.Naca4 as Naca4 - -data Foil a = Foil [Element a] String - -instance (Storable a) => Show (Foil a) where - show (Foil [el] name) = "{"++name++": " ++ show (1 + dim (fLengths el)) ++ " nodes}" - show (Foil els name) = "{"++name++": " ++ show (length els) ++ " elements, " ++ - show nodesPerEl ++ " nodes == "++show (sum nodesPerEl)++" total nodes}" - where - nodesPerEl = map (\x -> 1 + dim (fLengths x)) els - -data Element a = Element { fNodes :: (Vector a, Vector a) - , fLengths :: Vector a - , fAngles :: Vector a - , fMidpoints :: (Vector a, Vector a) - , fTangents :: (Vector a, Vector a) - , fNormals :: (Vector a, Vector a) - , fUnitNormals :: (Vector a, Vector a) - , fInits :: (Vector a, Vector a) - , fTails :: (Vector a, Vector a) - } - -toElement :: (Num (Vector a), RealFloat a, Container Vector a) => - [(a, a)] -> Element a -toElement xynodes = Element { fNodes = (xNodes, yNodes) - , fLengths = lengths - , fAngles = zipVectorWith atan2 yTangents xTangents - , fMidpoints = (xMids, yMids) - , fTangents = (xTangents, yTangents) - , fNormals = (xNormals, yNormals) - , fUnitNormals = (xUnitNormals, yUnitNormals) - , fInits = (xInits, yInits) - , fTails = (xTails, yTails) - } - where - n = (dim xNodes) - 1 - (xNodes, yNodes) = (\(xs,ys) -> (fromList xs, fromList ys)) $ unzip xynodes - (xInits, yInits) = (subVector 0 n xNodes, subVector 0 n yNodes) - (xTails, yTails) = (subVector 1 n xNodes, subVector 1 n yNodes) - (xTangents, yTangents) = (xTails - xInits, yTails - yInits) - (xMids, yMids) = (0.5*(xInits + xTails), 0.5*(yInits + yTails)) - (xNormals, yNormals) = (-yTangents, xTangents) - lengths = mapVector sqrt $ xTangents*xTangents + yTangents*yTangents - (xUnitNormals, yUnitNormals) = (xNormals/lengths, yNormals/lengths) - - --- why isn't this standard??? -poorMansStrip :: String -> String -poorMansStrip str = reverse $ dropWhile (== ' ') $ reverse $ dropWhile (== ' ') str - -getUIUCFoil :: (Num (Vector a), Read a, RealFloat a, Container Vector a) => - String -> IO (Either String (Foil a)) -getUIUCFoil name' = do - let name = poorMansStrip name' - let file = "http://www.ae.illinois.edu/m-selig/ads/coord/" ++ name ++ ".dat" - dl <- simpleHTTP (getRequest file) >>= getResponseBody - return (parseRawFoil dl name) - -parseRawFoil :: (Num (Vector a), Read a, RealFloat a, Container Vector a) => - String -> String -> Either String (Foil a) -parseRawFoil raw name - -- bad data - | any (\x -> 2 /= length x) (concat groupsOfLines) = Left (raw ++ "\nError parsing the above data") - -- single element - | length elements > 0 = Right (Foil (map toElement elements) name) - | otherwise = Left (raw ++ "\nError parsing the above data") - where - elements = map (map (\(x:y:[]) -> (read x, read y))) groupsOfLines - - -- group lines by splitting at empty lines - groupsOfLines :: [[[String]]] - groupsOfLines = f (lines raw) - where - f [] = [] - f ([]:xs) = f xs - f xs = (map words fst'):(f snd') - where - (fst', snd') = break (\x -> 0 == length x) xs - -loadFoil :: FilePath -> IO (Either String (Foil Double)) -loadFoil filename' = do - let filename = poorMansStrip filename' - exists <- doesFileExist filename - if exists - then do rawData <- readFile filename - return (parseRawFoil rawData filename) -- use filename as name - else return (Left ("file \"" ++ filename ++ "\" couldn't be found")) - -panelizeNaca4 :: (Enum a, Floating (Vector a), RealFloat a, Field a) => - Naca4.Naca4 a -> Int -> Foil a -panelizeNaca4 foil nPanels = Foil [toElement $ [(1,0)]++reverse lower++[(0,0)]++upper++[(1,0)]] - (Naca4.naca4_name foil) - where - (upper, lower) = unzip $ map (Naca4.coords foil) xcs - xcs = toList $ fst3 $ bunchPanels (Naca4.yt foil) (Naca4.dyt foil) xcs0 0 0 - xcs0 = fromList $ init $ tail $ toList $ linspace nXcs (0,1) - nXcs = (nPanels + (nPanels `mod` 2)) `div` 2 + 1 - -bunchPanels :: (Enum a, Floating (Vector a), Floating a, Ord a, Field a) => - (a -> a) -> (a -> a) -> Vector a -> Int -> Int -> (Vector a, Int, Int) -bunchPanels yt dyt xcs nIter nBadSteps - | nIter > 300 = error "panel buncher exceeded 300 iterations" - | nBadSteps > 1000 = error "panel buncher exceeded 1000 bad steps" - | sum (toList (abs deltaXcs)) < 1e-12 = (xcs, nIter, nBadSteps) - | otherwise = bunchPanels yt dyt goodXcs (nIter+1) (nBadSteps + length badOnes) - where - (badOnes, goodXcs:_) = break (\xs -> all (>0) (toList xs)) nextXcs - nextXcs = map (\alpha -> xcs + (scale alpha deltaXcs)) (map (2.0**) [0,-1..]) - deltaXcs = xcsStep yt dyt xcs - -xcsStep :: (Enum a, Floating (Vector a), Field a) => (a -> a) -> (a -> a) -> Vector a -> Vector a -xcsStep yt dyt xcs = flatten $ -(linearSolveLS mat2 (asColumn rs)) - where - n = dim xcs - - xs = join [fromList [0], xcs, fromList[1]] - ys = join [fromList [0], mapVector yt xcs, fromList[0]] - dxs = (subVector 1 (n+1) xs) - (subVector 0 (n+1) xs) - dys = (subVector 1 (n+1) ys) - (subVector 0 (n+1) ys) - deltas = sqrt (dxs*dxs + dys*dys) - dydxs = mapVector dyt xcs - - mat = r1 - r0 - where - r0 = fromBlocks [[(1><1)[0,0..] - eye = ident n - frontBunchingParam = 2.0 - diff = (fromBlocks [[zeros, eye]]) - (fromBlocks [[eye*(1+frontBunchingParam/(fromIntegral n)), zeros]]) - - mat2 = diff <> mat - rs = diff <> deltas diff --git a/Numeric/HFoil/Naca4.hs b/Numeric/HFoil/Naca4.hs deleted file mode 100644 --- a/Numeric/HFoil/Naca4.hs +++ /dev/null @@ -1,68 +0,0 @@ -{-# OPTIONS_GHC -Wall #-} - -module Numeric.HFoil.Naca4( Naca4(..) - , coords - , yt - , dyt - , naca4 - ) where - -naca4 :: (Read a, Fractional a) => String -> Naca4 a -naca4 name@(m_:p_:t0:t1:[]) = Naca4 m p t ("NACA "++name) - where - m = 0.01 * read [m_] - p = 0.1 * read [p_] - t = 0.01 * read (t0:[t1]) -naca4 _ = error "not a 4 digit airfoil" - --- m: max camber in hundredths of chord --- p: position of max camber in tenths of chord --- t: max thickness in hundredths of chord -data Naca4 a = Naca4 { naca4_m :: a - , naca4_p :: a - , naca4_t :: a - , naca4_name :: String - } deriving Show - --- xc: x/chord -yc :: (Ord a, Fractional a) => Naca4 a -> a -> a -yc (Naca4 {naca4_p = p, naca4_m = m}) xc - | xc < 0 = error "xc < 0" - | xc <= p = m/(p*p)*(2*p*xc - xc*xc) - | xc <= 1 = m/((1-p)*(1-p))*((1-2*p) + 2*p*xc - xc*xc) - | otherwise = error "xc > 1" - -dyc :: (Ord a, Fractional a) => Naca4 a -> a -> a -dyc (Naca4 {naca4_p = p, naca4_m = m}) xc - | xc < 0 = error "xc < 0" - | xc <= p = m/(p*p)*(2*p - 2*xc) - | xc <= 1 = m/((1-p)*(1-p))*(2*p - 2*xc) - | otherwise = error "xc > 1" - -yt :: (Ord a, Floating a) => Naca4 a -> a -> a -yt (Naca4 {naca4_t = t}) xc - | xc < 0 = error "xc < 0" - | xc > 1 = error $ "xc > 1" - | otherwise = 5*t*(0.2969*sqrt(xc) - 0.1260*(xc) - 0.3537*(xc)**2 + 0.2843*(xc)**3 - 0.1015*(xc)**4) - -dyt :: (Ord a, Floating a) => Naca4 a -> a -> a -dyt (Naca4 {naca4_t = t}) xc - | xc < 0 = error "xc < 0" - | xc > 1 = error "xc > 1" - | otherwise = 5*t*(0.5*0.2969/sqrt(xc) - 0.1260 - 2*0.3537*(xc) + 3*0.2843*(xc)**2 - 4*0.1015*(xc)**3) - -coords :: (Ord a, Floating a) => Naca4 a -> a -> ((a,a),(a,a)) -coords foil xc - | naca4_m foil == 0 = ((xc,yt_), (xc,-yt_)) - | otherwise = ((xu,yu ), (xl, yl )) - where - yt_ = yt foil xc - yc_ = yc foil xc - - yu = yc_ + yt_ * (cos theta) - yl = yc_ - yt_ * (cos theta) - - xu = xc - yt_ * (sin theta) - xl = xc + yt_ * (sin theta) - - theta = atan $ dyc foil xc diff --git a/Numeric/HFoil/Repl.hs b/Numeric/HFoil/Repl.hs deleted file mode 100644 --- a/Numeric/HFoil/Repl.hs +++ /dev/null @@ -1,125 +0,0 @@ -{-# OPTIONS_GHC -Wall #-} - -module Numeric.HFoil.Repl( run - ) where - -import System.Console.Haskeline hiding(display) -import Graphics.Gloss.Interface.IO.Animate hiding(scale, Vector) -import Control.Monad.IO.Class -import Control.Concurrent(forkIO) -import Control.Concurrent.MVar(newMVar, readMVar, swapMVar) - -import Numeric.HFoil.Foil -import Numeric.HFoil.Naca4 -import Numeric.HFoil.Drawing -import Numeric.HFoil.Flow - ----- configuration -nPanels :: Int -nPanels = 200 - -xSize, ySize :: Int -xSize = 800 -ySize = 500 - -data Config = Config { confForces :: Bool - , confKuttas :: Bool - , confNormals :: Bool - } - -defaultConfig :: Config -defaultConfig = Config { confForces = False - , confKuttas = False - , confNormals = False - } - -run :: IO () -run = do - let naca0 = "2412" - alfaDeg0 = 4 - flow0 = solveFlow (panelizeNaca4 (naca4 naca0) nPanels) (pi/180*alfaDeg0) - mpics <- newMVar $ [drawSolution flow0] - - putStrLn "Welcome to hfoil\n" - - _ <- forkIO $ runInputT defaultSettings - $ topLoop (\pics -> swapMVar mpics pics >>= (\_ -> return ())) defaultConfig - - animateIO - (InWindow - "hfoil" -- window title - (xSize, ySize) -- window size - (10, 650)) -- window position - black -- background color - (\_ -> readMVar mpics >>= return . pictures) -- draw function - -foilLoop :: ([Picture] -> IO ()) -> Config -> Foil Double -> InputT IO () -foilLoop draw conf foil@(Foil _ name) = do - minput <- getInputLine $ "\ESC[1;32m\STXhfoil."++name++">> \ESC[0m\STX" - case minput of - Nothing -> return () - Just "quit" -> do outputStrLn "gloss won't let you quit :(\ntry ctrl-c or hit ESC in drawing window" - foilLoop draw conf foil - Just ('a':'l':'f':'a':' ':[]) -> do outputStrLn $ "unrecognized command" - foilLoop draw conf foil - Just ('a':'l':'f':'a':' ':alphaDeg) -> do let flow = solveFlow foil (pi/180*(read alphaDeg)) - forces = case (confForces conf) of - True -> [drawForces flow] - False -> [] - kuttas = case (confKuttas conf) of - True -> [drawKuttas flow] - False -> [] - normals = case (confNormals conf) of - True -> [drawNormals (solFoil flow)] - False -> [] - liftIO $ draw $ forces++kuttas++normals++[drawSolution flow] - foilLoop draw conf foil - Just ('f':'o':'r':'c':'e':'s':[]) -> do - let newConf = conf {confForces = not (confForces conf)} - outputStrLn $ "force drawing set to "++ show (not (confForces conf)) - foilLoop draw newConf foil - Just ('k':'u':'t':'t':'a':'s':[]) -> do - let newConf = conf {confKuttas = not (confKuttas conf)} - outputStrLn $ "kutta drawing set to "++ show (not (confKuttas conf)) - foilLoop draw newConf foil - Just ('n':'o':'r':'m':'a':'l':'s':[]) -> do - let newConf = conf {confNormals = not (confNormals conf)} - outputStrLn $ "normals drawing set to "++ show (not (confNormals conf)) - foilLoop draw newConf foil - Just "" -> return () - Just input -> do outputStrLn $ "unrecognized command \"" ++ input ++ "\"" - foilLoop draw conf foil - -topLoop :: ([Picture] -> IO ()) -> Config -> InputT IO () -topLoop draw conf = do - minput <- getInputLine "\ESC[1;32m\STXhfoil>> \ESC[0m\STX" - case minput of - Nothing -> return () - Just "quit" -> do outputStrLn "gloss won't let you quit :(\ntry ctrl-c or hit ESC in drawing window" - topLoop draw conf - Just ('n':'a':'c':'a':' ':spec) -> do parseNaca draw conf spec - topLoop draw conf - Just ('l':'o':'a':'d':' ':name) -> do - foil <- liftIO (loadFoil name) - case foil of Left errMsg -> outputStrLn errMsg - Right foil' -> do liftIO $ draw [drawFoil foil', drawNormals foil'] - foilLoop draw conf foil' - topLoop draw conf - Just ('u':'i':'u':'c':' ':name) -> do - efoil <- liftIO (getUIUCFoil name) - case efoil of Left errMsg -> outputStrLn errMsg - Right foil -> do liftIO $ draw [drawFoil foil, drawNormals foil] - foilLoop draw conf foil - topLoop draw conf - - Just "" -> topLoop draw conf - Just input -> do outputStrLn $ "unrecognized command \"" ++ input ++ "\"" - topLoop draw conf - -parseNaca :: ([Picture] -> IO ()) -> Config -> String -> InputT IO () -parseNaca draw conf str - | length str == 4 = do let foil = panelizeNaca4 (naca4 str :: Naca4 Double) nPanels - liftIO $ draw [drawFoil foil, drawNormals foil] - foilLoop draw conf foil - | otherwise = do outputStrLn $ "Not 4 digits" - return () diff --git a/apps/Main.hs b/apps/Main.hs new file mode 100644 --- /dev/null +++ b/apps/Main.hs @@ -0,0 +1,8 @@ +{-# OPTIONS_GHC -Wall #-} + +module Main where + +import HFoil.Repl + +main :: IO () +main = run diff --git a/hfoil.cabal b/hfoil.cabal --- a/hfoil.cabal +++ b/hfoil.cabal @@ -1,31 +1,31 @@ Name: hfoil -Version: 0.1.2 +Version: 0.2.0 Stability: Experimental Synopsis: Hess-Smith panel code for inviscid 2-d airfoil analysis License: BSD3 License-file: LICENSE Author: Greg Horn Maintainer: gregmainland@gmail.com -Copyright: Greg Horn, 2012 +Copyright: (c) Greg Horn 2012, 2015 Category: Numerical, Science Build-type: Simple -Cabal-version: >=1.6 +Cabal-version: >= 1.8 Description: { Library and command line REPL with plotting to do simple inviscid hess-smith panel code. . Features include: . -* Single and multi-element airfoils +* Cheap and shameless xfoil ripoff for relp/plotting interface . * Naca 4-series support with Gauss-Newton paneling . * Broken UIUC database integration (type \"uiuc [foilname]\") . -* Shameless xfoil ripoff for relp/plotting interface +* Haskeline interface with tab-completion (oooh) . -* Only works with development version of Gloss that's not on Hackage yet +* Single and multi-element airfoils . -* Haskeline interface with tab-completion +* Inviscid, incompressible, 2-dimensional flow only . . To get started, do cabal install or whatever, then run the \"hfoil\" binary. @@ -38,22 +38,24 @@ Default: True Library - Exposed-modules: Numeric.HFoil - Numeric.HFoil.Naca4 - Numeric.HFoil.Foil - Numeric.HFoil.Drawing - Numeric.HFoil.Flow - Numeric.HFoil.Repl + hs-source-dirs: src + Exposed-modules: HFoil + HFoil.Naca4 + HFoil.Foil + HFoil.Drawing + HFoil.Flow + HFoil.Repl -- Other-modules: - Build-depends: gloss >= 1.7 && <= 1.8, - base >= 4 && < 5, - MissingH >= 1.1, - hmatrix >= 0.12, - haskeline >= 0.6, - transformers >= 0.2, - directory >= 1.1, - HTTP >= 4000.2 + Build-depends: base >= 4 && < 5 + , not-gloss >= 0.7.2.0 + , linear + , hmatrix >= 0.12 + , haskeline >= 0.6 + , transformers >= 0.2 + , directory >= 1.1 + , HTTP >= 4000.2 + , parsec >= 3 Ghc-Options: -Wall Ghc-Prof-Options: -prof -auto-all @@ -65,13 +67,14 @@ else Buildable: False + Build-depends: base >= 4 && < 5 + , hfoil >= 0.1.2 + + hs-source-dirs: apps Main-Is: Main.hs - -- (gloss needs -O2 -threaded) + -- (not-gloss needs -O2 -threaded) -- does it really? Ghc-Options: -O2 -threaded -Wall GHC-Prof-Options: -auto-all -caf-all -rtsopts - - if os(OSX) - extra-lib-dirs: /usr/lib source-repository head type: git diff --git a/src/HFoil.hs b/src/HFoil.hs new file mode 100644 --- /dev/null +++ b/src/HFoil.hs @@ -0,0 +1,27 @@ +{- | + Module : Hfoil + Description : Top level module + + This is the top level module which exports the API + -} + +{-# OPTIONS_GHC -Wall #-} + +module HFoil + ( -- * Airfoils + module HFoil.Foil + -- * Flow solution + , module HFoil.Flow + -- * Naca 4 utilities + , module HFoil.Naca4 + -- * Drawing utilities (gloss backend) + , module HFoil.Drawing + -- * Interactive command line application + , module HFoil.Repl + ) where + +import HFoil.Naca4 +import HFoil.Foil +import HFoil.Drawing +import HFoil.Flow +import HFoil.Repl diff --git a/src/HFoil/Drawing.hs b/src/HFoil/Drawing.hs new file mode 100644 --- /dev/null +++ b/src/HFoil/Drawing.hs @@ -0,0 +1,149 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE FlexibleContexts #-} + +module HFoil.Drawing + ( drawSolution + , drawFoil + , drawNormals + , drawForces + , drawKuttas + , drawOnce + ) where + +import Numeric.LinearAlgebra hiding( Element, scale, i ) +import Foreign.Storable ( Storable ) +import qualified Numeric.LinearAlgebra as LA +import Text.Printf ( printf ) +import Linear ( V3(..) ) + +import Vis + +import HFoil.Flow +import HFoil.Foil + +cpScale :: Fractional a => a +cpScale = -0.25 + +normalLengths :: Fractional a => a +normalLengths = 0.01 + +drawLine :: Num a => Color -> [(a,a)] -> VisObject a +drawLine col coords = Line (map (\(x,y) -> V3 x y 0) coords) col + +drawCircle :: Num a => Color -> (a, a) -> a -> VisObject a +drawCircle col (x,y) size = Trans (V3 x y 0) $ Sphere size Solid col +--drawCircle col (x,y) size = scale (border*(fromIntegral xSize)) (border*(fromIntegral xSize)) +-- $ translate (-0.5 + realToFrac x) (realToFrac y) +-- $ color col +-- $ Circle size + +drawLineV :: (Num a, Storable a) => Color -> (Vector a, Vector a) -> VisObject a +drawLineV col (vx, vy) = Line (zipWith (\x y -> V3 x y 0) (toList vx) (toList vy)) col + +drawFoil :: (Num a, Storable a) => Foil a -> VisObject a +drawFoil (Foil elements _) = VisObjects $ map drawElement elements + +drawElement :: (Num a, Storable a) => Element a -> VisObject a +drawElement element = drawLineV white (fNodes element) + +drawNormals :: Foil Double -> VisObject Double +drawNormals (Foil elements _) = VisObjects $ map (\(xy0, xy1) -> drawLine green [xy0, xy1]) (zip xy0s xy1s) + where + xy0s = zip (toList xm) (toList ym) + xy1s = zip (toList (xm + (LA.scale normalLengths xUnitNormal))) (toList (ym + (LA.scale normalLengths yUnitNormal))) + + (xUnitNormal, yUnitNormal) = (\(x,y) -> (vjoin x, vjoin y)) $ unzip $ map fUnitNormals elements + (xm, ym) = (\(x,y) -> (vjoin x, vjoin y)) $ unzip $ map fMidpoints elements + +colorFun :: (Fractional a, Real a) => a -> a -> a -> Color +colorFun min' max' x' = makeColor (1-x) (1-x) x 1 + where + x = realToFrac $ (x' - min') / (max'-min') + +drawKuttas :: (Real a, Fractional a, Storable a) => FlowSol a -> VisObject a +drawKuttas flow = VisObjects $ concatMap (\(k0,k1) -> [circ k0, circ k1]) kis + where + kis = solKuttaIndices flow + (xs',ys') = unzip $ map fMidpoints $ (\(Foil els _) -> els) (solFoil flow) + xs = vjoin xs' + ys = vjoin ys' + circ k = drawCircle yellow (xs @> k, ys @> k) 0.006 + +drawForces :: FlowSol Double -> VisObject Double +drawForces flow = VisObjects $ map (\(xy0, xy1, cp) -> drawLine (colorFun minCp maxCp cp) [xy0, xy1]) + $ zip3 xy0s xy1s (toList (solCps flow)) + where + xy0s = zip (toList xm) (toList ym) + xy1s = zip (toList (xm + xPressures)) (toList (ym + yPressures)) + (xPressures, yPressures) = (\(x,y) -> (LA.scale c x/lengths, LA.scale c y/lengths)) (solForces flow) + lengths = vjoin $ map fLengths $ (\(Foil x _) -> x) $ solFoil flow + (xm, ym) = (\(x,y) -> (vjoin x, vjoin y)) $ unzip $ map fMidpoints $ (\(Foil x _) -> x) $ solFoil flow + + c = 0.1 + + maxCp = maxElement (solCps flow) + minCp = minElement (solCps flow) + + +drawColoredFoil :: (Num a, Storable a) => [Color] -> Foil a -> VisObject a +drawColoredFoil colors foil@(Foil elements _) = VisObjects $ zipWith drawColoredElement colors' elements + where + colors' = groupSomethingByFoil foil colors + +drawColoredElement :: (Num a, Storable a) => [Color] -> Element a -> VisObject a +drawColoredElement colors element = VisObjects $ map (\(xy0, xy1, col) -> drawLine col [xy0, xy1]) (zip3 xy0s xy1s colors) + where + xys = (\(x,y) -> zip (toList x) (toList y)) $ fNodes element + xy0s = tail xys + xy1s = init xys + +groupSomethingByFoil :: Storable a => Foil a -> [b] -> [[b]] +groupSomethingByFoil (Foil elements _) somethings = f somethings (map (LA.dim . fAngles) elements) + where + f xs (n:ns) = (take n xs):(f (drop n xs) ns) + f [] []= [] + f _ _ = error "uh oh (groupSomethingByFoil)" + +drawSolution :: FlowSol Double -> VisObject Double +drawSolution flow = VisObjects $ onscreenText ++ + [ drawColoredFoil colors foil + , drawCircle white (fst $ solCenterPressure flow, snd $ solCenterPressure flow) 0.006 + , drawCircle white (fst $ solCenterPressure flow, 0) 0.006 + , drawCircle green (0.25,0) 0.006 + ] ++ zipWith (\x y -> drawLineV red (x, y)) xs + (takesV (map LA.dim xs) (LA.scale cpScale cps)) -- cp graph + where + foil@(Foil elements name) = solFoil flow + cps = solCps flow + + xs = map (fst . fMidpoints) elements + + onscreenText = + zipWith (\s k -> Text2d s (30,fromIntegral $ 30*k) Fixed9By15 (makeColor 1 1 1 1)) + msgs (reverse [1..length msgs]) + + msgs = [ name + , printf ("alpha: %.6f deg") ((solAlpha flow)*180/pi) + , printf ("Cl: %.6f") (solCl flow) + , printf ("Cd: %.6f") (solCd flow) + , printf ("Cm: %.6f (c/4, 0)") (solCm flow) + ] + + colors = map (colorFun (minElement cps) (maxElement cps)) (toList cps) + + +drawOnce :: Real a => [VisObject a] -> IO () +drawOnce pics = display (defaultOpts {optWindowName = "hfoil"}) (VisObjects pics) + +-- let line = plot_lines_values ^= [[ (xc, yt (naca4 "0012") xc) +-- | xc <- [0,0.01..0.99::Double]]] +-- $ plot_lines_title ^= "naca 0012" +-- $ defaultPlotLines +-- +-- chart = layout1_title ^= "naca yo" +-- $ layout1_plots ^= [Left (toPlot line)] +-- $ defaultLayout1 +-- +-- renderableToWindow (toRenderable chart) 640 480 +-- _ <- renderableToPNGFile (toRenderable chart) 640 480 "mDiv_vs_tc.png" +-- return () diff --git a/src/HFoil/Flow.hs b/src/HFoil/Flow.hs new file mode 100644 --- /dev/null +++ b/src/HFoil/Flow.hs @@ -0,0 +1,205 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE FlexibleContexts #-} + +module HFoil.Flow + ( FlowSol(..) + , solveFlow + ) where + +import Data.Packed.ST (newMatrix, writeMatrix, freezeMatrix) +import Control.Monad ( forM_ ) +import Control.Monad.ST ( runST ) +import Numeric.LinearAlgebra hiding(i) +import Foreign.Storable ( Storable ) + +import HFoil.Foil + +data FlowSol a = FlowSol { solFoil :: Foil a + , solVs :: Vector a + , solCps :: Vector a + , solVorticities :: [a] + , solAlpha :: a + , solForces :: (Vector a, Vector a) + , solForce :: (a,a) + , solCl :: a + , solCd :: a + , solCm :: a -- moment about (0,0.25) + , solCenterPressure :: (a,a) + , solKuttaIndices :: [(Int,Int)] + } + +solveFlow :: (Num (Vector a), RealFloat a, Field a) => Foil a -> a -> FlowSol a +solveFlow foil@(Foil elements _) alpha = + FlowSol { solFoil = foil + , solVs = vs + , solCps = cps + , solVorticities = vorticities + , solAlpha = alpha + , solForces = (xForces, yForces) + , solForce = (xf,yf) + , solCl = cl + , solCd = cd + , solCm = cm + , solCenterPressure = (xCp, yCp) + , solKuttaIndices = kuttaIndices + } + where + -- surface velocities + (vs,vorticities) = ( (mapVector (\q -> cos(q - alpha)) angles) + (mV <> qsGamma) + , map (qsGamma @>) [n-(length elements),n-1] + ) + where + (mA, mV, b) = getAVb geometries angles kuttaIndices alpha + qsGamma = flatten $ linearSolve mA b + n = dim angles + -- pressure coefficients + cps = 1 - vs*vs + + -- forces and force coefficients + xForces = -cps*(vjoin $ map (fst . fNormals) elements) + yForces = -cps*(vjoin $ map (snd . fNormals) elements) + xf = sumElements xForces + yf = sumElements yForces + + cd = xf*(cos alpha) + yf*(sin alpha) + cl = -xf*(sin alpha) + yf*(cos alpha) + + -- midpoints + (xs,ys) = (\(x,y) -> (vjoin x, vjoin y)) $ unzip $ map fMidpoints elements + + -- centers of pressure + (xCp, yCp) = ( (sumElements (xs*yForces)) / (sumElements yForces) + , (sumElements (ys*xForces)) / (sumElements xForces) + ) + + -- moment about cp and quarter chord + cm = sumElements $ (mapVector (\z -> z - 0.25) xs)*yForces - ys*xForces + + kuttaIndices = ki 0 (map dim angles') + where + ki _ [] = [] + ki n0 (n:ns) = (n0,n0+n-1):(ki (n0+n) ns) + + angles = vjoin angles' + angles' = map fAngles elements + + -- (mSL,mCL,mSB,mCB) + geometries = setupGeometries angles (xms, yms) (xns, yns) (xnps, ynps) + where + vjoinTuples (x,y) = (vjoin x, vjoin y) + (xms,yms) = vjoinTuples $ unzip $ map fMidpoints elements + (xns,yns) = vjoinTuples $ unzip $ map fInits elements + (xnps,ynps) = vjoinTuples $ unzip $ map fTails elements + + +-- make influence matrix A and also the matrix V where V*[sources; vortex] == tangential speeds +getAVb :: (Floating a, Num (Vector a), Container Vector a, Storable a) => + (Matrix a, Matrix a, Matrix a, Matrix a) + -> Vector a + -> [(Int, Int)] + -> a + -> (Matrix a, Matrix a, Matrix a) +getAVb (mSL,mCL,mSB,mCB) angles kuttaIndices alpha = + ( scale (1/(2*pi)) $ fromBlocks [[ mAij, fromColumns vsAin] + ,[fromRows vsAnj, fromLists ann]] + , scale (1/(2*pi)) $ fromBlocks [[ mSB - mCL + , fromColumns vsTin]] + , asColumn $ vjoin $ (mapVector (\q -> sin $ q - alpha) angles): + (map (\(n0,nf) -> fromList [-cos((angles @> n0) - alpha) - cos((angles @> nf) - alpha)]) kuttaIndices) + ) -- (for middle entry mV, mSL + mCB == mAij) + where + n = rows mSL + + -- sources influence matrix Aij + mAij = mSL + mCB + + -- vortices influence vectors Ains + vsAin = map (\(n0,nf) -> fromList $ map sumElements $ toRows $ subMatrix (0,n0) (n, nf-n0+1) m) kuttaIndices + where + m = (mCL - mSB) + + -- vortices tangential velocity outputs + vsTin = map (\(n0,nf) -> fromList $ map sumElements $ toRows $ subMatrix (0,n0) (n, nf-n0+1) m) kuttaIndices + where + m = (mSL + mCB) + + -- sources kutta condition influence vector Anj + vsAnj = map (\(n0,nf) -> (getRow n0 mSB) - (getRow n0 mCL) + (getRow nf mSB) - (getRow nf mCL)) kuttaIndices + + -- vortices kutta condition influence scalars ann + ann = (flip map) kuttaIndices $ \(ni0,nif) -> (flip map) kuttaIndices $ \(nj0,njf) -> + sumElements $ (getSubRow ni0 (nj0,njf) mSL) + (getSubRow ni0 (nj0,njf) mCB) + (getSubRow nif (nj0,njf) mSL) + (getSubRow nif (nj0,njf) mCB) + + getRow i mat = flatten $ subMatrix (i,0) (1,n) mat + getSubRow i (j0,jf) mat = flatten $ subMatrix (i,j0) (1,jf-j0+1) mat + + +{- +calcuate 4 matrices which will be useful + +mSL(i,j) = sin(qi - ij) * ln(rij+1/rij) +mCL(i,j) = cos(qi - ij) * ln(rij+1/rij) +mSB(i,j) = sin(qi - ij) * beta(i,j) +mCB(i,j) = cos(qi - ij) * beta(i,j) + +where for i==j: ln(r/r) = 0 + beta = pi +are explicitly set +-} +setupGeometries :: (RealFloat t, Storable t) => + Vector t + -> (Vector t, Vector t) + -> (Vector t, Vector t) + -> (Vector t, Vector t) + -> (Matrix t, Matrix t, Matrix t, Matrix t) +setupGeometries angles (xms, yms) (xns, yns) (xnps, ynps) = runST $ do + let n = dim angles + mSL' <- newMatrix 0 n n + mCL' <- newMatrix 0 n n + mSB' <- newMatrix 0 n n + mCB' <- newMatrix 0 n n + + forM_ [0..n-1] $ \i -> forM_ [0..n-1] $ \j -> do + let qi = angles @> i + qj = angles @> j + + xmi = xms @> i + ymi = yms @> i + xnj = xns @> j + ynj = yns @> j + xnjp = xnps @> j + ynjp = ynps @> j + + s = sin (qi - qj) + c = cos (qi - qj) + + l + | i == j = 0 + | otherwise = log(r1/r0) + where + r1 = distance (xmi, ymi) (xnjp, ynjp) + r0 = distance (xmi, ymi) (xnj, ynj ) + distance (x1,y1) (x0,y0) = sqrt $ dx*dx + dy*dy + where + dx = x1 - x0 + dy = y1 - y0 + + b + | i == j = pi + | otherwise = atan2 (dyjp*dxj - dxjp*dyj) + (dxjp*dxj + dyjp*dyj) + where + dyj = ymi - ynj + dxj = xmi - xnj + dyjp = ymi - ynjp + dxjp = xmi - xnjp + + _ <- writeMatrix mSL' i j (s*l) + _ <- writeMatrix mCL' i j (c*l) + _ <- writeMatrix mSB' i j (s*b) + writeMatrix mCB' i j (c*b) + sl' <- freezeMatrix mSL' + cl' <- freezeMatrix mCL' + sb' <- freezeMatrix mSB' + cb' <- freezeMatrix mCB' + return (sl',cl',sb',cb') diff --git a/src/HFoil/Foil.hs b/src/HFoil/Foil.hs new file mode 100644 --- /dev/null +++ b/src/HFoil/Foil.hs @@ -0,0 +1,246 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE FlexibleContexts #-} + +module HFoil.Foil + ( Foil(..) + , Element(..) + , panelizeNaca4 + , loadFoil + , getUIUCFoil + ) where + +import System.Directory ( doesFileExist ) +import Numeric.LinearAlgebra hiding ( Element ) +import Network.HTTP ( simpleHTTP, getRequest, getResponseBody ) +import Foreign.Storable ( Storable ) +import Text.Read ( readMaybe ) +import Text.Parsec ( Parsec, ParsecT, Stream, (), parse + , option, char, try, lookAhead, manyTill, digit, endOfLine + , many1, optional, anyChar, sepBy1, many, oneOf, skipMany ) + +import qualified HFoil.Naca4 as Naca4 + +data Foil a = Foil [Element a] String + +instance (Storable a) => Show (Foil a) where + show (Foil [el] name) = "{"++name++": " ++ show (1 + dim (fLengths el)) ++ " nodes}" + show (Foil els name) = "{"++name++": " ++ show (length els) ++ " elements, " ++ + show nodesPerEl ++ " nodes == "++show (sum nodesPerEl)++" total nodes}" + where + nodesPerEl = map (\x -> 1 + dim (fLengths x)) els + +data Element a = Element { fNodes :: (Vector a, Vector a) + , fLengths :: Vector a + , fAngles :: Vector a + , fMidpoints :: (Vector a, Vector a) + , fTangents :: (Vector a, Vector a) + , fNormals :: (Vector a, Vector a) + , fUnitNormals :: (Vector a, Vector a) + , fInits :: (Vector a, Vector a) + , fTails :: (Vector a, Vector a) + } + + +-- make sure the nodes aren't reversed +toElement :: (Num (Vector a), RealFloat a, Container Vector a) + => [(a, a)] -> Element a +toElement nodes + | diff < 0 = el + | otherwise = toElement' (reverse nodes) + where + el = toElement' nodes + angles = toList (fAngles el) + diff = sum $ zipWith (-) angles (drop 1 angles) + +toElement' :: (Num (Vector a), RealFloat a, Container Vector a) + => [(a, a)] -> Element a +toElement' xynodes = + Element { fNodes = (xNodes, yNodes) + , fLengths = lengths + , fAngles = zipVectorWith atan2 yTangents xTangents + , fMidpoints = (xMids, yMids) + , fTangents = (xTangents, yTangents) + , fNormals = (xNormals, yNormals) + , fUnitNormals = (xUnitNormals, yUnitNormals) + , fInits = (xInits, yInits) + , fTails = (xTails, yTails) + } + where + n = (dim xNodes) - 1 + (xNodes, yNodes) = (\(xs,ys) -> (fromList xs, fromList ys)) $ unzip xynodes + (xInits, yInits) = (subVector 0 n xNodes, subVector 0 n yNodes) + (xTails, yTails) = (subVector 1 n xNodes, subVector 1 n yNodes) + (xTangents, yTangents) = (xTails - xInits, yTails - yInits) + (xMids, yMids) = (0.5*(xInits + xTails), 0.5*(yInits + yTails)) + (xNormals, yNormals) = (-yTangents, xTangents) + lengths = mapVector sqrt $ xTangents*xTangents + yTangents*yTangents + (xUnitNormals, yUnitNormals) = (xNormals/lengths, yNormals/lengths) + + +-- why isn't this standard??? +poorMansStrip :: String -> String +poorMansStrip str = reverse $ dropWhile (== ' ') $ reverse $ dropWhile (== ' ') str + +getUIUCFoil :: String -> IO (Either String (Foil Double)) +getUIUCFoil name' = do + let name = poorMansStrip name' + file = "http://m-selig.ae.illinois.edu/ads/coord/" ++ name ++ ".dat" + dl <- simpleHTTP (getRequest file) >>= getResponseBody + return (parseRawFoil dl name) + +loadFoil :: FilePath -> IO (Either String (Foil Double)) +loadFoil filename' = do + let filename = poorMansStrip filename' + exists <- doesFileExist filename + if exists + then do rawData <- readFile filename + return (parseRawFoil rawData filename) -- use filename as name + else return (Left ("file \"" ++ filename ++ "\" couldn't be found")) + +fst3 :: (a,b,c) -> a +fst3 (x,_,_) = x + +panelizeNaca4 :: (Enum a, Floating (Vector a), RealFloat a, Field a) => + Naca4.Naca4 a -> Int -> Foil a +panelizeNaca4 foil nPanels = Foil [toElement $ [(1,0)]++reverse lower++[c0]++upper++[(1,0)]] + (Naca4.naca4_name foil) + where + c0 = fst (Naca4.coords foil 0) + (upper, lower) = unzip $ map (Naca4.coords foil) xcs + xcs = toList $ fst3 $ bunchPanels (Naca4.yt foil) (Naca4.dyt foil) xcs0 0 0 + xcs0 = fromList $ init $ tail $ toList $ linspace nXcs (0,1) + nXcs = (nPanels + (nPanels `mod` 2)) `div` 2 + 1 + +bunchPanels :: forall a . (Enum a, Floating (Vector a), Floating a, Ord a, Field a) => + (a -> a) -> (a -> a) -> Vector a -> Int -> Int -> (Vector a, Int, Int) +bunchPanels yt dyt xcs nIter nBadSteps + | nIter > 300 = error "panel buncher exceeded 300 iterations" + | nBadSteps > 1000 = error "panel buncher exceeded 1000 bad steps" + | sum (toList (abs deltaXcs)) < 1e-12 = (xcs, nIter, nBadSteps) + | otherwise = bunchPanels yt dyt goodXcs (nIter+1) (nBadSteps + length badOnes) + where + (badOnes, goodXcs:_) = break (\xs -> all (>0) (toList xs)) nextXcs + nextXcs :: [Vector a] + nextXcs = map (\alpha -> xcs + (scale alpha deltaXcs)) (map (2.0**) [0,-1..]) + deltaXcs :: Vector a + deltaXcs = xcsStep yt dyt xcs + +xcsStep :: (Enum a, Floating (Vector a), Field a) => (a -> a) -> (a -> a) -> Vector a -> Vector a +xcsStep yt dyt xcs = flatten $ -(linearSolveLS mat2 (asColumn rs)) + where + n = dim xcs + + -- x coords + xs = vjoin [fromList [0], xcs, fromList [1]] + -- y coords + ys = vjoin [fromList [0], mapVector yt xcs, fromList [0]] + -- x deltas + dxs = (subVector 1 (n+1) xs) - (subVector 0 (n+1) xs) + -- y deltas + dys = (subVector 1 (n+1) ys) - (subVector 0 (n+1) ys) + -- delta magnitures + deltas = sqrt (dxs*dxs + dys*dys) + -- slopes + dydxs = mapVector dyt xcs + + mat = r1 - r0 + where + r0 = fromBlocks [[(1><1)[0,0..] + eye = ident n + frontBunchingParam = 2.0 + diff = (fromBlocks [[zeros, eye]]) - (fromBlocks [[eye*(1+frontBunchingParam/(fromIntegral n)), zeros]]) + + mat2 = diff <> mat + rs = diff <> deltas + +discardIdenticalNeighbors :: Eq a => [a] -> [a] +discardIdenticalNeighbors (x0:x1:xs) + | x0 == x1 = discardIdenticalNeighbors (x1:xs) + | otherwise = x0 : discardIdenticalNeighbors (x1:xs) +discardIdenticalNeighbors xs = xs + +stitchIncompleteElements :: [[(Double, Double)]] -> [[(Double,Double)]] +stitchIncompleteElements (as:bs:others) + | and [ x0a < 0.5 + , x0b < 0.5 + , 0.5 < xFa + , 0.5 < xFb + ] = (discardIdenticalNeighbors ((reverse as) ++ bs)) : stitchIncompleteElements others + | and [ 0.5 < x0a + , 0.5 < x0b + , xFa < 0.5 + , xFb < 0.5 + ] = (discardIdenticalNeighbors (as ++ reverse bs)) : stitchIncompleteElements others + | otherwise = as : bs : stitchIncompleteElements others + where + x0a = fst $ head as + xFa = fst $ last as + + x0b = fst $ head bs + xFb = fst $ last bs +stitchIncompleteElements [] = [] +stitchIncompleteElements xs = xs + +parseRawFoil :: String -> String -> Either String (Foil Double) +parseRawFoil raw name = foil + where + foil = case parse foilP name raw of + Left pe -> Left (raw ++ "\nError parsing the above data: " ++ (show pe)) + Right (_, els) -> Right (Foil (map toElement (stitchIncompleteElements els)) name) + +headerP :: Parsec String () String +headerP = manyTill anyChar (try (lookAhead elementsP)) + +foilP :: Parsec String () (String, [[(Double, Double)]]) +foilP = do + s <- headerP "header" + els <- elementsP "elements" + return (s, els) + +elementP :: Parsec String () [(Double, Double)] +elementP = many2 coordP + where + many2 :: (Stream s m t) => ParsecT s u m a -> ParsecT s u m [a] + many2 p = do + x0 <- p + x1 <- p + xs <- many p + return (x0:x1:xs) + + +elementsP :: Parsec String () [[(Double, Double)]] +elementsP = sepBy1 elementP endOfLine + +doubleP :: Parsec String () Double +doubleP = do + neg <- option ' ' (char '-') + leading <- option "0" (many1 digit) + dot' <- char '.' + trailing <- option "0" (many1 digit) + let doublish = neg : leading ++ dot' : trailing + case readMaybe doublish of + Just x -> return x + Nothing -> error $ "failed to read this supposed double: " ++ show doublish + +coordP :: Parsec String () (Double, Double) +coordP = do + let space' = oneOf [' ', '\t'] + spaces' = skipMany space' "space-like" + spaces' "leading whitespace" + x <- doubleP "first coordinate" + spaces' "middle whitespace" + y <- doubleP "second coordinate" + spaces' "trailing whitespace" + _ <- optional endOfLine "end of line" + return (x,y) diff --git a/src/HFoil/Naca4.hs b/src/HFoil/Naca4.hs new file mode 100644 --- /dev/null +++ b/src/HFoil/Naca4.hs @@ -0,0 +1,71 @@ +{-# OPTIONS_GHC -Wall #-} + +module HFoil.Naca4 + ( Naca4(..) + , coords + , yt + , dyt + , naca4 + ) where + +import Text.Read ( readMaybe ) + +naca4 :: (Read a, Fractional a) => String -> Maybe (Naca4 a) +naca4 name@(m_:p_:t0:t1:[]) = do + m <- fmap (0.01 *) $ readMaybe [m_] + p <- fmap (0.1 *) $ readMaybe [p_] + t <- fmap (0.01 *) $ readMaybe [t0,t1] + return $ Naca4 m p t ("NACA " ++ name) +naca4 _ = Nothing + +-- m: max camber in hundredths of chord +-- p: position of max camber in tenths of chord +-- t: max thickness in hundredths of chord +data Naca4 a = Naca4 { naca4_m :: a + , naca4_p :: a + , naca4_t :: a + , naca4_name :: String + } deriving Show + +-- xc: x/chord +yc :: (Ord a, Fractional a) => Naca4 a -> a -> a +yc (Naca4 {naca4_p = p, naca4_m = m}) xc + | xc < 0 = error "xc < 0" + | xc < p = m/(p*p)*(2*p*xc - xc*xc) + | xc <= 1 = m/((1-p)*(1-p))*((1-2*p) + 2*p*xc - xc*xc) + | otherwise = error "xc > 1" + +dyc :: (Ord a, Fractional a) => Naca4 a -> a -> a +dyc (Naca4 {naca4_p = p, naca4_m = m}) xc + | xc < 0 = error "xc < 0" + | xc < p = m/(p*p)*(2*p - 2*xc) + | xc <= 1 = m/((1-p)*(1-p))*(2*p - 2*xc) + | otherwise = error "xc > 1" + +yt :: (Ord a, Floating a) => Naca4 a -> a -> a +yt (Naca4 {naca4_t = t}) xc + | xc < 0 = error "xc < 0" + | xc > 1 = error $ "xc > 1" + | otherwise = 5*t*(0.2969*sqrt(xc) - 0.1260*(xc) - 0.3537*(xc)**2 + 0.2843*(xc)**3 - 0.1015*(xc)**4) + +dyt :: (Ord a, Floating a) => Naca4 a -> a -> a +dyt (Naca4 {naca4_t = t}) xc + | xc < 0 = error "xc < 0" + | xc > 1 = error "xc > 1" + | otherwise = 5*t*(0.5*0.2969/sqrt(xc) - 0.1260 - 2*0.3537*(xc) + 3*0.2843*(xc)**2 - 4*0.1015*(xc)**3) + +coords :: (Ord a, Floating a) => Naca4 a -> a -> ((a,a),(a,a)) +coords foil xc + | naca4_m foil == 0 = ((xc,yt_), (xc,-yt_)) + | otherwise = ((xu,yu ), (xl, yl )) + where + yt_ = yt foil xc + yc_ = yc foil xc + + yu = yc_ + yt_ * (cos theta) + yl = yc_ - yt_ * (cos theta) + + xu = xc - yt_ * (sin theta) + xl = xc + yt_ * (sin theta) + + theta = atan $ dyc foil xc diff --git a/src/HFoil/Repl.hs b/src/HFoil/Repl.hs new file mode 100644 --- /dev/null +++ b/src/HFoil/Repl.hs @@ -0,0 +1,226 @@ +{-# OPTIONS_GHC -Wall #-} +-- {-# Language FlexibleContexts #-} + +module HFoil.Repl + ( run + ) where + +import Control.Concurrent ( forkIO ) +import Control.Concurrent.MVar ( newMVar, readMVar, swapMVar ) +import Control.Monad.IO.Class ( MonadIO, liftIO ) +import Control.Monad.Trans.Class ( lift ) +import Control.Monad.Trans.State.Strict ( StateT, evalStateT, get, modify ) +import Data.List ( isPrefixOf ) +import Linear ( Quaternion(..), V3(..) ) +import System.Console.Haskeline ( InputT, runInputT, defaultSettings, getInputLine, outputStrLn, setComplete ) +import System.Console.Haskeline.Completion ( CompletionFunc, Completion, completeWord, simpleCompletion ) +import Text.Read ( readMaybe ) + +import Vis + +import HFoil.Foil +import HFoil.Naca4 +import HFoil.Drawing +import HFoil.Flow + +nPanels :: Int +nPanels = 200 + +-- configuration +data Config = Config { confForces :: Bool + , confKuttas :: Bool + , confNormals :: Bool + } + +defaultConfig :: Config +defaultConfig = Config { confForces = False + , confKuttas = False + , confNormals = False + } + +data Mode = TopMode | FoilMode + +foilCommands :: [(String, String)] +foilCommands = + [ ("alfa", "alfa [#]") + , ("forces", "forces") + , ("kuttas", "kuttas") + , ("normals", "normals") + , ("help", "help") + ] + +topCommands :: [(String, String)] +topCommands = + [ ("naca", "naca xxxx") + , ("load", "load [filename]") + , ("uiuc", "uiuc [foil name]") + ] + +topHelp :: InputT (StateT Mode IO) () +topHelp = mapM_ (outputStrLn . snd) topCommands + +foilHelp :: InputT (StateT Mode IO) () +foilHelp = mapM_ (outputStrLn . snd) foilCommands + +comp :: CompletionFunc (StateT Mode IO) +comp = completeWord Nothing " \t" searchFunc + where + searchFunc :: String -> (StateT Mode IO) [Completion] + searchFunc str = do + mode <- get + let wordList = case mode of + TopMode -> map fst topCommands + FoilMode -> map fst foilCommands + return $ map simpleCompletion $ filter (str `isPrefixOf`) wordList + +run :: IO () +run = do + mpics <- newMVar $ [] + + putStrLn "Welcome to hfoil\n" + + let go :: InputT (StateT Mode IO) () + go = topLoop (\pics -> swapMVar mpics pics >>= (\_ -> return ())) defaultConfig + + settings = setComplete comp defaultSettings + _ <- forkIO $ flip evalStateT TopMode $ runInputT settings go + + let toScreen xs = + RotQuat (Quaternion 0 (V3 1 0 0)) + $ Trans (V3 (-0.5) 0 0) + $ VisObjects xs + -- $ VisObjects (Axes (0.3,15) : xs) + cam0 = + Camera0 + { phi0 = 90 + , theta0 = 90 + , rho0 = 2 + } + + animateIO + (defaultOpts {optWindowName = "hfoil", optInitialCamera = Just cam0}) + (\_ -> fmap toScreen (readMVar mpics)) + +data FoilState = + FoilState + { fsFlowSol :: Maybe (FlowSol Double) + , fsConf :: Config + , fsFoil :: Foil Double + } + +drawPicture :: MonadIO m => ([VisObject Double] -> IO ()) -> StateT FoilState m () +drawPicture draw = do + fs <- get + let conf = fsConf fs + foil = fsFoil fs + normals = case confNormals conf of + True -> [drawNormals foil] + False -> [] + + case fsFlowSol fs of + Nothing -> liftIO $ draw (drawFoil foil: normals) + Just flow -> do + let forces = case (confForces conf) of + True -> [drawForces flow] + False -> [] + kuttas = case (confKuttas conf) of + True -> [drawKuttas flow] + False -> [] + liftIO $ draw $ forces++kuttas++normals++[drawSolution flow] + +strip :: String -> String +strip = rstrip . lstrip + where + lstrip (' ':xs) = lstrip xs + lstrip x = x + + rstrip = reverse . lstrip . reverse + +foilLoop :: ([VisObject Double] -> IO ()) -> StateT FoilState (InputT (StateT Mode IO)) () +foilLoop draw = do + lift (lift (modify (const FoilMode))) + fs <- get + let foil@(Foil _ name) = fsFoil fs + conf = fsConf fs + drawPicture draw + minput <- lift $ getInputLine $ "\ESC[1;32m\STXhfoil."++name++">> \ESC[0m\STX" + + case fmap strip minput of + Nothing -> return () + Just "quit" -> do lift $ outputStrLn "not-gloss won't let you quit :(\ntry ctrl-c or hit ESC in drawing window" + foilLoop draw + Just ('a':'l':'f':'a':' ':alphaDeg') -> do + case readMaybe alphaDeg' of + Nothing -> do lift $ outputStrLn $ "parse fail on " ++ show alphaDeg' + foilLoop draw + Just alphaDeg -> do let flow :: FlowSol Double + flow = solveFlow foil (pi/180*alphaDeg) + modify (\fs' -> fs' {fsFlowSol = Just flow}) + foilLoop draw + Just "forces" -> do + let newConf = conf {confForces = not (confForces conf)} + lift $ outputStrLn $ "force drawing set to "++ show (not (confForces conf)) + modify (\fs' -> fs' {fsConf = newConf}) + foilLoop draw + Just "kuttas" -> do + let newConf = conf {confKuttas = not (confKuttas conf)} + lift $ outputStrLn $ "kutta drawing set to "++ show (not (confKuttas conf)) + modify (\fs' -> fs' {fsConf = newConf}) + foilLoop draw + Just "normals" -> do + let newConf = conf {confNormals = not (confNormals conf)} + lift $ outputStrLn $ "normals drawing set to "++ show (not (confNormals conf)) + modify (\fs' -> fs' {fsConf = newConf}) + foilLoop draw + Just "help" -> lift foilHelp >> foilLoop draw + Just "h" -> lift foilHelp >> foilLoop draw + Just "?" -> lift foilHelp >> foilLoop draw + Just "" -> return () + Just input -> do lift $ outputStrLn $ "unrecognized command \"" ++ input ++ "\"" + foilLoop draw + + +topLoop :: ([VisObject Double] -> IO ()) -> Config -> InputT (StateT Mode IO) () +topLoop draw conf = do + lift (modify (const TopMode)) + minput <- getInputLine "\ESC[1;32m\STXhfoil>> \ESC[0m\STX" + case minput of + Nothing -> return () + Just msg -> do runTop draw conf msg + topLoop draw conf + +runTop :: ([VisObject Double] -> IO ()) -> Config -> String -> InputT (StateT Mode IO) () +runTop draw conf msg = case strip msg of + "quit" -> outputStrLn "not-gloss won't let you quit :(\ntry ctrl-c or hit ESC in drawing window" + ('n':'a':'c':'a':' ':spec) -> do + case naca4 spec :: Maybe (Naca4 Double) of + Nothing -> outputStrLn "not a valid naca4" + Just n4 -> runFoil draw conf (panelizeNaca4 n4 nPanels) + ('l':'o':'a':'d':' ':name) -> do + mfoil <- liftIO (loadFoil name) + case mfoil of Left errMsg -> outputStrLn errMsg + Right foil -> do runFoil draw conf foil + ('u':'i':'u':'c':' ':name) -> do + efoil <- liftIO (getUIUCFoil name) + case efoil of Left errMsg -> outputStrLn errMsg + Right foil -> do let Foil els _ = foil + outputStrLn $ "got " ++ show (length els) ++ " elements" + runFoil draw conf foil + "help" -> topHelp + "h" -> topHelp + "?" -> topHelp + "" -> return () + other -> outputStrLn $ "unrecognized command \"" ++ other ++ "\"" + + +runFoil :: ([VisObject Double] -> IO ()) -> Config -> Foil Double -> InputT (StateT Mode IO) () +runFoil draw conf foil = do + let state0 = + FoilState + { fsFlowSol = Nothing + , fsConf = conf + , fsFoil = foil + } + let go :: StateT FoilState (InputT (StateT Mode IO)) () + go = foilLoop draw + flip evalStateT state0 go