support precision other than Float64 in mesh and quad

src/mesh.jl

@@ -195,9 +195,9 @@ end
 Base.getindex(msh::Mesh, ent::EntityKey) = getindex(msh, Domain(ent))
 
 """
-    meshgen(Ω, n)
-    meshgen(Ω, n_dict)
-    meshgen(Ω; meshsize)
+    meshgen(Ω, n; T = Float64)
+    meshgen(Ω, n_dict; T = Float64)
+    meshgen(Ω; meshsize, T = Float64)
 
 Generate a `Mesh` for the domain `Ω` where each curve is meshed using
 `n` elements. Passing a dictionary allows for a finer control; in such cases,
@@ -209,6 +209,8 @@ is computed as so as to obtain an *average* mesh size of `meshsize`. Note that
 the actual mesh size may vary significantly for each element if the
 parametrization is far from uniform.
 
+The mesh is created with primitive data of type `T`.
+
 This function requires the entities forming `Ω` to have an explicit
 parametrization.
 
@@ -217,7 +219,7 @@ parametrization.
     the quality of the underlying parametrization. For complex surfaces, you are
     better off using a proper mesher such as `gmsh`.
 """
-function meshgen(Ω::Domain, T = Float64, args...; kwargs...)
+function meshgen(Ω::Domain, args...; T = Float64, kwargs...)
     # extract the ambient dimension for these entities (i.e. are we in 2d or
     # 3d). Only makes sense if all entities have the same ambient dimension.
     N = ambient_dimension(first(Ω))

src/quadrature.jl

@@ -141,12 +141,12 @@ function Quadrature(msh::AbstractMesh; qorder)
 end
 
 @noinline function _build_quadrature!(
-    quad,
+    quad::Quadrature{N,T},
     els::AbstractVector{E},
     qrule::ReferenceQuadrature,
-) where {E}
-    N = ambient_dimension(quad)
-    x̂, ŵ = qrule() # nodes and weights on reference element
+) where {N,T,E}
+    x̂ = map(x̂ -> T.(x̂), qcoords(qrule))
+    ŵ = map(ŵ -> T.(ŵ), qweights(qrule))
     num_nodes = length(ŵ)
     M = geometric_dimension(domain(E))
     codim = N - M
@@ -158,8 +158,8 @@ end
             jac = jacobian(el, x̂i)
             μ = _integration_measure(jac)
             w = μ * ŵi
-            ν = codim == 1 ? _normal(jac) : nothing
-            qnode = QuadratureNode(x, w, ν)
+            ν = codim == 1 ? T.(_normal(jac)) : nothing
+            qnode = QuadratureNode(T.(x), T.(w), ν)
             push!(quad.qnodes, qnode)
         end
     end

test/meshgen_test.jl

@@ -13,7 +13,13 @@ using StaticArrays
     arc1 = Inti.parametric_curve(f, 0, 0.5)
     arc2 = Inti.parametric_curve(f, 0.5, 1)
     Γ = Inti.Domain(arc1, arc2)
-    msh = Inti.meshgen(Γ, 100)
+    msh = Inti.meshgen(Γ, 100; T = Float32)
+    quad = Inti.Quadrature(msh; qorder = 2)
+    @test Inti.integrate(x -> 1, quad) ≈ 2 * π * r
+    @test Inti.measure(arc1) ≈ π * r
+    @test Inti.measure(arc2) ≈ π * r
+
+    msh = Inti.meshgen(Γ, 100; T = Float64)
     quad = Inti.Quadrature(msh; qorder = 2)
     @test Inti.integrate(x -> 1, quad) ≈ 2 * π * r
     @test Inti.measure(arc1) ≈ π * r