fix exemple and implementation of DI extension

Author: yolhan83

Project.toml

@@ -6,7 +6,6 @@ version = "0.1.4"
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
@@ -21,7 +20,6 @@ BloodFlowTrixiDataInterpolationsExt = "DataInterpolations"
 DataInterpolations = "8.0"
 ForwardDiff = "0.10, 1"
 LinearAlgebra = "1.10, 1.11"
-QuadGK = "2.10, 2.11"
 StaticArrays = "1.9"
 Trixi = "0.10, 0.11"
 WriteVTK = "1.21"

exemples/Model2D/diexemple.jl

@@ -4,15 +4,22 @@ using DataInterpolations
 using BloodFlowTrixi
 using StaticArrays, LinearAlgebra
 using QuadGK
+using LinearAlgebra
 
 eq = BloodFlowEquations2D(; h = 0.1)
+xyz_data = [SA[cos(0.2*si),sin(0.2*si),si] for si in range(0,40,100)]
+
+curve = interpolate_curve(xyz_data)
+L = curve.t[end-1]
+println("curve length : $L")
+BloodFlowTrixi.curvature(s) = norm(DataInterpolations.derivative(curve,s,2))
 
 mesh = P4estMesh(
     (1,2),
     polydeg = 1,
     periodicity = (true,false),
     coordinates_min = (0.0,0.0),
-    coordinates_max = (2*pi,40.0),
+    coordinates_max = (2*pi,L),
     initial_refinement_level = 4
 )
 
@@ -47,8 +54,8 @@ cb = CallbackSet(
     dt_adapt,analyse
 )
 
-sol = solve(ode, SSPRK33(),dt = dt_adapt(ode),callback= cb)
+sol = solve(ode, SSPRK33(),dt = dt_adapt(ode),callback= cb,saveat = 0.03,save_everystep = false)
 
 # artery center-line
-xyz_data = [[cos(0.2*si),sin(0.2*si),si] for si in range(0,40,100)]
+
 res = get3DData(eq,xyz_data,semi,sol,1)

ext/BloodFlowTrixiDataInterpolationsExt.jl

@@ -7,14 +7,15 @@ module BloodFlowTrixiDataInterpolationsExt
         using ..DataInterpolations
     end
     using StaticArrays, LinearAlgebra
-    using ForwardDiff, QuadGK
-    function BloodFlowTrixi.get3DData(eq::BloodFlowEquations2D,curve_data::AbstractArray,semi,sol,time_index ::Int = 1;vtk ::Bool=false,out ::T="./datas") where {T<:AbstractString}
+    using ForwardDiff
+    function BloodFlowTrixi.interpolate_curve(curve_data::AbstractArray)
         N = length(curve_data)
         quadinterp = QuadraticSpline(curve_data,range(0,1,N))
-        L = quadgk(s->norm(ForwardDiff.derivative(quadinterp,s)),0,1)[1]
-        s_data = range(0,L,N)
-        newinterp = CubicSpline(curve_data,s_data)
-        curve = SmoothArcLengthInterpolation(newinterp;m=length(s_data),in_place=false)
+        curve = SmoothArcLengthInterpolation(quadinterp;m=N,in_place=false)
+        return curve
+    end
+    function BloodFlowTrixi.get3DData(eq::BloodFlowEquations2D,curve_data::AbstractArray,semi,sol,time_index ::Int = 1;vtk ::Bool=false,out ::T="./datas") where {T<:AbstractString}
+        curve = interpolate_curve(curve_data)
         tanj(s) = ForwardDiff.derivative(curve,s)
         function nor(s) 
             res= ForwardDiff.derivative(tanj,s)
@@ -23,11 +24,11 @@ module BloodFlowTrixiDataInterpolationsExt
                 a,b,c = tanj(s)
                 # return a any normal vector
                 if a != 0
-                    return [-b,a,0]/sqrt(a^2+b^2)
+                    return SA[-b,a,0]/sqrt(a^2+b^2)
                 elseif b != 0
-                    return [-b,a,0]/sqrt(a^2+b^2)
+                    return SA[-b,a,0]/sqrt(a^2+b^2)
                 else
-                    return [-sign(c),0,0]
+                    return SA[-sign(c),0,0]
                 end
             end
             return res/n
@@ -36,5 +37,5 @@ module BloodFlowTrixiDataInterpolationsExt
         er(theta,s) = cos(theta).*nor(s) .+ sin(theta).*∧(tanj(s),nor(s))
         return BloodFlowTrixi.get3DData(eq,s->curve(s),er,semi,sol,time_index;vtk=vtk,out=out)
     end
-    export get3DData
+    export get3DData,interpolate_curve
 end

src/2DModel/viz.jl

@@ -148,4 +148,8 @@ function get3DData(eq::BloodFlowEquations2D,semi,sol,time_index ::Int = 1;vtk ::
     ∧(v,w) = SA[v[2]*w[3]-v[3]*w[2],v[3]*w[1]-v[1]*w[3],v[1]*w[2]-v[2]*w[1]]
     er(theta,s) = cos(theta).*nor(s) .+ sin(theta).*∧(tanj(s),nor(s))
     return get3DData(eq,curve,er,semi,sol,time_index;vtk=vtk,out=out)
+end
+
+function interpolate_curve(xyz_data)
+    @error "add DataInterpolations.jl to use this function"
 end

src/BloodFlowTrixi.jl

@@ -12,11 +12,11 @@ Docs under https://yolhan83.github.io/BloodFlowTrixi.jl
 $(isnothing(get(ENV, "CI", nothing)) ? ("\n" * "Package local path: " * pathof(BloodFlowTrixi)) : "") 
 """
 module BloodFlowTrixi
-    using Trixi,WriteVTK,StaticArrays,QuadGK
+    using Trixi,WriteVTK,StaticArrays
     # Write your package code here.
     abstract type AbstractBloodFlowEquations{NDIMS, NVARS} <:Trixi.AbstractEquations{NDIMS, NVARS} end
     include("1DModel/1dmodel.jl")
     include("2DModel/2dmodel.jl")
 
-    export BloodFlowEquations1D,BloodFlowEquations1DOrd2, BloodFlowEquations2D,flux_nonconservative,source_term_simple,source_term_simple_ord2,boundary_condition_pressure_in,initial_condition_simple,friction,pressure,radius,boundary_condition_outflow,boundary_condition_slip_wall,get3DData
+    export BloodFlowEquations1D,BloodFlowEquations1DOrd2, BloodFlowEquations2D,flux_nonconservative,source_term_simple,source_term_simple_ord2,boundary_condition_pressure_in,initial_condition_simple,friction,pressure,radius,boundary_condition_outflow,boundary_condition_slip_wall,get3DData,interpolate_curve
 end