occt sectionPerimeter

InlineOCCT.hs

{-# LANGUAGE TemplateHaskell #-}

module InlineOCCT where

import Data.Acquire
import qualified Data.Map as Map
import Foreign hiding (with)
import Language.C.Inline.Context
import Language.C.Inline.Cpp
import Language.C.Inline.HaskellIdentifier
import Language.C.Types as C
import Language.C.Types.Parse (CIdentifier (..))
import Language.Haskell.TH as TH
import OpenCascade.GP.Types
import OpenCascade.TopoDS.Types
import Waterfall.Internal.ToOpenCascade (v3ToDir, v3ToPnt, v3ToVertex)

getHsVariable :: String -> HaskellIdentifier -> TH.ExpQ
getHsVariable err s = do
  mbHsName <- TH.lookupValueName $ unHaskellIdentifier s
  case mbHsName of
    Nothing ->
      fail $
        "Cannot capture Haskell variable "
          ++ unHaskellIdentifier s
          ++ ", because it's not in scope. ("
          ++ err
          ++ ")"
    Just hsName -> TH.varE hsName

occtContext :: Context
occtContext = cppCtx {ctxTypesTable = ctxTypesTable cppCtx <> tt, ctxAntiQuoters = aq}

tt :: TypesTable
tt =
  Map.fromList
    [ (f "Vertex", [t|Ptr Vertex|]),
      (f "Pnt", [t|Ptr Pnt|]),
      (f "Dir", [t|Ptr Dir|]),
      (f "Solid", [t|Ptr Shape|])
    ]

f str = TypeName $ either (error "tt") id $ cIdentifierFromString True str

p str = Ptr [] (TypeSpecifier mempty (f str))

aq :: AntiQuoters
aq =
  Map.fromList
    [ ("dir", SomeAntiQuoter dirAntiQuoter),
      ("pnt", SomeAntiQuoter pntAntiQuoter),
      ("solid", SomeAntiQuoter solidAntiQuoter)
    ]

dirAntiQuoter :: AntiQuoter HaskellIdentifier
dirAntiQuoter =
  AntiQuoter
    { aqParser = do
        hId <- C.parseIdentifier
        useCpp <- C.parseEnableCpp
        let cId = mangleHaskellIdentifier useCpp hId
        return (cId, p "gp_Dir", hId),
      aqMarshaller = \_purity _cTypes _cTy cId -> do
        hsExp <- getHsVariable "occtContext" cId
        hsExp' <- [|with (v3ToDir $(return hsExp))|]
        hsTy <- [t|Ptr Dir|]
        return (hsTy, hsExp')
    }

pntAntiQuoter :: AntiQuoter HaskellIdentifier
pntAntiQuoter =
  AntiQuoter
    { aqParser = do
        hId <- C.parseIdentifier
        useCpp <- C.parseEnableCpp
        let cId = mangleHaskellIdentifier useCpp hId
        return (cId, p "gp_Pnt", hId),
      aqMarshaller = \_purity _cTypes _cTy cId -> do
        hsExp <- getHsVariable "occtContext" cId
        hsExp' <- [|with (v3ToPnt $(return hsExp))|]
        hsTy <- [t|Ptr Pnt|]
        return (hsTy, hsExp')
    }

solidAntiQuoter :: AntiQuoter HaskellIdentifier
solidAntiQuoter =
  AntiQuoter
    { aqParser = do
        hId <- C.parseIdentifier
        useCpp <- C.parseEnableCpp
        let cId = mangleHaskellIdentifier useCpp hId
        return (cId, p "TopoDS_Shape", hId),
      aqMarshaller = \_purity _cTypes _cTy cId -> do
        hsExp <- getHsVariable "occtContext" cId
        hsExp' <- [|with (acquireSolid $(return hsExp))|]
        hsTy <- [t|Ptr Shape|]
        return (hsTy, hsExp')
    }

main.hs

import SectionPerimeter
import Waterfall
import Linear
main = do
  print =<< sectionPerimeter unitSphere 1 0
  print =<< sectionPerimeter unitCube (V3 1 0 0) 0.5
  print =<< sectionPerimeter unitCube 1 (V3 1 0 0)
{-
6.283185307179585
4.0
4.242640687119286 -- 3 sqrt 2
-}

package.yaml

executable:
  main:
    main.hs
  source-dirs:
    - .
dependencies:
  - base
  - inline-c
  - inline-c-cpp
  - resourcet
  - containers
  - template-haskell
  - opencascade-hs
  - waterfall-cad
  - linear

SectionPerimeter.hs

{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE ViewPatterns #-}

module SectionPerimeter where

import Foreign.C.Types
import InlineOCCT
import qualified Language.C.Inline as C
import qualified Language.C.Inline.Cpp as Cpp
import Linear.V3 (V3 (..))
import qualified OpenCascade.GP.Vec as GPVec
import Waterfall
import Waterfall.Internal.Solid

C.context occtContext
Cpp.include "<BRepExtrema_DistShapeShape.hxx>"
Cpp.include "<gp_Pnt.hxx>"
Cpp.include "<gp_Pln.hxx>"
Cpp.include "<gp_Vec.hxx>"
Cpp.include "<TopoDS_Shape.hxx>"
Cpp.include "<TopoDS.hxx>"
Cpp.include "<BRepBuilderAPI_MakeFace.hxx>"
Cpp.include "<BRepGProp.hxx>"
Cpp.include "<BRepAlgoAPI_Section.hxx>"
Cpp.include "<BRep_Tool.hxx>"
Cpp.include "<TopExp_Explorer.hxx>"
Cpp.include "<BRepGProp_Cinert.hxx>"

-- | @p = sectionPerimeter1 s n x@
--
-- section a solid @s@ with the plane defined by normal @n@ and point @x@,
-- giving the perimeter @p@ (strictly speaking, the total length of all wires in that plane)
sectionPerimeter :: Solid -> V3 Double -> V3 Double -> IO CDouble
sectionPerimeter solid n p =
  [C.block| double { 
    gp_Pln pl = gp_Pln(* $pnt:p,* $dir:n);
    TopoDS_Face planeFace = BRepBuilderAPI_MakeFace(pl);
    BRepAlgoAPI_Section section(* $solid:solid,planeFace);
    section.Build();
    
    if (!section.IsDone()) {
        return 0.0;
    }
    
    TopoDS_Shape result = section.Shape();
    Standard_Real totalLength = 0.0;
    TopExp_Explorer edgeExplorer(result, TopAbs_EDGE);
    for (; edgeExplorer.More(); edgeExplorer.Next()) {
        TopoDS_Edge edge = TopoDS::Edge(edgeExplorer.Current());
        Standard_Real first, last;
        Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, first, last);
        gp_Pnt start, end;
        curve->D0(first, start);
        curve->D0(last, end);
        auto props = BRepGProp_Cinert();
        if (pl.Distance(start) < 1e-6 && pl.Distance(end) < 1e-6) {
            BRepGProp::LinearProperties(edge, props);
            totalLength += props.Mass();
        }
    }
  return totalLength;
 }
|]