Skip to content

Commit

Permalink
Added unit_normal
Browse files Browse the repository at this point in the history
  • Loading branch information
@fverdugo
fverdugo committed Dec 2, 2024
1 parent a24525f commit 6b072bc
Show file tree
Hide file tree
Showing 3 changed files with 596 additions and 374 deletions.
158 changes: 132 additions & 26 deletions src/domain.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1791,55 +1791,161 @@ end
# unit_normal(domain,codomain,glue)
#end
#
function map_unit_normal(φ_fun,ϕ_fun,n)
q -> begin
p = φ_fun(q)
J = ForwardDiff.jacobian(ϕ_fun,p)
Jt = transpose(J)
pinvJt = transpose(inv(Jt*J)*Jt)
v = pinvJt*n
m = sqrt(inner(v,v))
if m < eps()
return zero(v)
else
return v/m
end
end
end

#function map_unit_normal(φ_fun,ϕ_fun,n)
# q -> begin
# p = φ_fun(q)
# J = ForwardDiff.jacobian(ϕ_fun,p)
# Jt = transpose(J)
# pinvJt = transpose(inv(Jt*J)*Jt)
# v = pinvJt*n
# m = sqrt(inner(v,v))
# if m < eps()
# return zero(v)
# else
# return v/m
# end
# end
#end

#function unit_normal(mesh::AbstractMesh,d)
# D = num_dims(mesh)
# @assert d == D-1
# phiD = physical_map(mesh,D)
# phidinv = inverse_physical_map(mesh,d)
# phidD = reference_map(mesh,d,D)
# ctype_to_refface = GT.reference_faces(mesh,D)
# Drid_to_lface_to_n_data= map(ctype_to_refface) do refface
# boundary = refface |> GT.geometry |> GT.boundary
# boundary |> GT.outwards_normals # TODO also rename?
# end
# topo = topology(mesh)
# dface_to_Dfaces_data, dface_to_ldfaces_data = GT.face_incidence_ext(topo,d,D)
# n = quantity() do index
# dface = face_index(index,d)
# Dface = face_index(index,D)
# Dface_to_Drid = get_symbol!(index,face_reference_id(mesh,D),"Dface_to_Drid")
# dface_to_Dfaces = get_symbol!(index,dface_to_Dfaces_data,"dface_to_Dfaces")
# dface_to_ldfaces = get_symbol!(index,dface_to_ldfaces_data,"dface_to_ldfaces")
# Drid_to_lface_to_n = get_symbol!(index,Drid_to_lface_to_n_data,"Drid_to_lface_to_n")
# expr = @term begin
# Drid = $Dface_to_Drid[$Dface]
# Dfaces = $dface_to_Dfaces[$dface]
# Dface_around = GalerkinToolkit.find_face_adound($Dface,Dfaces)
# ldface = $dface_to_ldfaces[$dface][Dface_around]
# Drid_to_lface_to_n[Drif][ldface]
# end
# p = Drid_to_lface_to_n_data |> eltype |> eltype |> zero
# expr_term([d,D],expr,p,index)
# end
# #x -> call(map_unit_normal,(phidD∘phidinv)(x),phiD(x),n)
# x -> call(map_unit_normal,phiD(x),n)
#end


#function unit_normal(mesh::AbstractMesh,d)
# D = num_dims(mesh)
# @assert d == D-1
# phiD = physical_map(mesh,D)
# phidinv = inverse_physical_map(mesh,d)
# phidD = reference_map(mesh,d,D)
# ctype_to_refface = GT.reference_faces(mesh,D)
# Drid_to_lface_to_n_data= map(ctype_to_refface) do refface
# boundary = refface |> GT.geometry |> GT.boundary
# boundary |> GT.outwards_normals # TODO also rename?
# end
# topo = topology(mesh)
# dface_to_Dfaces_data, dface_to_ldfaces_data = GT.face_incidence_ext(topo,d,D)
# n = quantity() do index
# dface = face_index(index,d)
# Dface = face_index(index,D)
# Dface_to_Drid = get_symbol!(index,face_reference_id(mesh,D),"Dface_to_Drid")
# dface_to_Dfaces = get_symbol!(index,dface_to_Dfaces_data,"dface_to_Dfaces")
# dface_to_ldfaces = get_symbol!(index,dface_to_ldfaces_data,"dface_to_ldfaces")
# Drid_to_lface_to_n = get_symbol!(index,Drid_to_lface_to_n_data,"Drid_to_lface_to_n")
# expr = @term begin
# Drid = $Dface_to_Drid[$Dface]
# Dfaces = $dface_to_Dfaces[$dface]
# Dface_around = GalerkinToolkit.find_face_adound($Dface,Dfaces)
# ldface = $dface_to_ldfaces[$dface][Dface_around]
# Drid_to_lface_to_n[Drif][ldface]
# end
# p = Drid_to_lface_to_n_data |> eltype |> eltype |> zero
# expr_term([d,D],expr,p,index)
# end
# #x -> call(map_unit_normal,(phidD∘phidinv)(x),phiD(x),n)
# x -> call(map_unit_normal,phiD(x),n)
#end


function unit_normal(mesh::AbstractMesh,d)
D = num_dims(mesh)
@assert d == D-1
phiD = physical_map(mesh,D)
phidinv = inverse_physical_map(mesh,d)
phidD = reference_map(mesh,d,D)
phiDinv = inverse_physical_map(mesh,D)
ctype_to_refface = GT.reference_faces(mesh,D)
Drid_to_lface_to_n_data= map(ctype_to_refface) do refface
boundary = refface |> GT.geometry |> GT.boundary
boundary |> GT.outwards_normals # TODO also rename?
end
topo = topology(mesh)
dface_to_Dfaces_data, dface_to_ldfaces_data = GT.face_incidence_ext(topo,d,D)
n = quantity() do index
quantity() do index
dface = face_index(index,d)
Dface = face_index(index,D)
Dface_to_rid = get_symbol!(a.index,face_reference_id(mesh,D),"Dface_to_Drid")
dface_to_Dfaces = get_symbol!(index(a),dface_to_Dfaces_data,"dface_to_Dfaces")
dface_to_ldfaces = get_symbol!(index(a),dface_to_ldfaces_data,"dface_to_ldfaces")
Drid_to_lface_to_n = get_symbol!(index(a),Drid_to_lface_to_n_data,"Drid_to_lface_to_n")
Dface_to_Drid = get_symbol!(index,face_reference_id(mesh,D),"Dface_to_Drid")
dface_to_Dfaces = get_symbol!(index,dface_to_Dfaces_data,"dface_to_Dfaces")
dface_to_ldfaces = get_symbol!(index,dface_to_ldfaces_data,"dface_to_ldfaces")
Drid_to_lface_to_n = get_symbol!(index,Drid_to_lface_to_n_data,"Drid_to_lface_to_n")
expr = @term begin
Drid = $Dface_to_Drid[$Dface]
Dfaces = $dface_to_Dfaces[$dface]
Dface_around = GalerkinToolkit.find_face_adound($Dface,Dfaces)
ldface = $dface_to_ldfaces[$dface][Dface_around]
Drid_to_lface_to_n[Drif][ldface]
$Drid_to_lface_to_n[Drid][ldface]
end
p = Drid_to_lface_to_n_data |> eltype |> eltype |> zero
expr_term([d,D],expr,p,index)
n = expr_term([d,D],expr,p,index)
n_ref = unit_normal_term(n)
ϕinv_fun = term(phiDinv,index)
n_phys = binary_call_term(∘,n_ref,ϕinv_fun)
n_phys
end
call(map_unit_normal,phidD∘phidinv,phiD,n)
#phiD = physical_map(mesh,D)
#phiDinv = inverse_physical_map(mesh,D)
#ctype_to_refface = GT.reference_faces(mesh,D)
#Drid_to_lface_to_n_data= map(ctype_to_refface) do refface
# boundary = refface |> GT.geometry |> GT.boundary
# boundary |> GT.outwards_normals # TODO also rename?
#end
#topo = topology(mesh)
#dface_to_Dfaces_data, dface_to_ldfaces_data = GT.face_incidence_ext(topo,d,D)
#quantity() do index
# dface = face_index(index,d)
# Dface = face_index(index,D)
# Dface_to_Drid = get_symbol!(index,face_reference_id(mesh,D),"Dface_to_Drid")
# dface_to_Dfaces = get_symbol!(index,dface_to_Dfaces_data,"dface_to_Dfaces")
# dface_to_ldfaces = get_symbol!(index,dface_to_ldfaces_data,"dface_to_ldfaces")
# Drid_to_lface_to_n = get_symbol!(index,Drid_to_lface_to_n_data,"Drid_to_lface_to_n")
# expr = @term begin
# Drid = $Dface_to_Drid[$Dface]
# Dfaces = $dface_to_Dfaces[$dface]
# Dface_around = GalerkinToolkit.find_face_adound($Dface,Dfaces)
# ldface = $dface_to_ldfaces[$dface][Dface_around]
# Drid_to_lface_to_n[Drif][ldface]
# end
# p = Drid_to_lface_to_n_data |> eltype |> eltype |> zero
# n = expr_term([d,D],expr,p,index)
# ϕ_fun = term(phiD,index)
#

# ϕinv_fun = term(phiDinv,index)
# n_ref = binary_call_term(map_unit_normal,ϕ_fun,n)
# n_phys = binary_call_term(∘,n_ref,ϕinv_fun)
# n_phys
#end
end


#
#function unit_normal(domain::ReferenceDomain,codomain::PhysicalDomain,glue::BoundaryGlue)
# Γref = domain
Expand Down
65 changes: 60 additions & 5 deletions src/symbolics.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,7 +131,7 @@ function call(f,args::AbstractTerm...)
elseif length(args) == 2
binary_call_term(f,args...)
else
multivar_call_term(f,args...)
multivar_call_term(f,args)
end
end

Expand Down Expand Up @@ -242,11 +242,29 @@ free_dims(a::BinaryCallTerm) = union(free_dims(a.arg1),free_dims(a.arg2))

multivar_call_term(args...) = MultivarCallTerm(args...)

struct MultiVarCallTerm{A,B} <: AbstractTerm
struct MultivarCallTerm{A,B} <: AbstractTerm
callee::A
args::B
end

AbstractTrees.children(a::MultivarCallTerm) = (a.callee,a.args...)

index(a::MultivarCallTerm) = index(a.args[1])

function expression(a::MultivarCallTerm)
exprs = map(expression,a.args)
g = get_symbol!(index(a),a.callee,"callee")
:($g($(exprs...)))
end

function prototype(a::MultivarCallTerm)
ps = map(prototype,a.args)
return_prototype(a.callee,ps...)
end

free_dims(a::MultivarCallTerm) = union(map(free_dims,a.args)...)


boundary_term(args...) = BoundaryTerm(args...)

struct BoundaryTerm{A,B,C,D} <: AbstractTerm
Expand Down Expand Up @@ -593,7 +611,7 @@ index(a::DiscreteFunctionTerm) = index(a.functions)

const LinearOperators = Union{typeof(call),typeof(ForwardDiff.gradient),typeof(ForwardDiff.jacobian)}

function expression(c::BinaryCallTerm{<:LinearOperators,<:DiscreteFunctionTerm,<:ReferencePointTerm})
function expression(c::BinaryCallTerm{<:LinearOperators,<:DiscreteFunctionTerm,<:Any})
f = c.callee
a = c.arg1
b = c.arg2
Expand Down Expand Up @@ -678,15 +696,15 @@ end
# expr_term(dims,expr,p,index(a))
#end

function expression(c::BinaryCallTerm{typeof(call),<:PhysicalMapTerm,<:ReferencePointTerm})
function expression(c::BinaryCallTerm{typeof(call),<:PhysicalMapTerm,<:Any})
f = c.callee
a = discrete_function_term(c.arg1)
b = c.arg2
fab = call(f,a,b)
expression(fab)
end

function expression(c::BinaryCallTerm{typeof(ForwardDiff.jacobian),<:PhysicalMapTerm,<:ReferencePointTerm})
function expression(c::BinaryCallTerm{typeof(ForwardDiff.jacobian),<:PhysicalMapTerm,<:Any})
f = c.callee
a = discrete_function_term(c.arg1)
b = c.arg2
Expand Down Expand Up @@ -893,6 +911,43 @@ function find_face_adound(face,faces)
findfirst(i->i==face,faces)
end

unit_normal_term(args...) = UnitNormalTerm(args...)

struct UnitNormalTerm{A} <: AbstractTerm
outwards_normals::A
end

free_dims(a::UnitNormalTerm) = free_dims(a.outwards_normals)
prototype(a::UnitNormalTerm) = prototype(a.outwards_normals)
index(a::UnitNormalTerm) = index(a.outwards_normals)

function expression(c::BinaryCallTerm{typeof(call),<:UnitNormalTerm,<:Any})
f = c.callee
D = num_dims(mesh(index(c)))
a = c.arg1
x = c.arg2
phiD = physical_map_term(Val{D}(),index(c))
J = call(ForwardDiff.jacobian,phiD,x)
n = call(map_unit_normal,J,a.outwards_normals)
expression(n)
end

function map_unit_normal(J,n)
Jt = transpose(J)
pinvJt = transpose(inv(Jt*J)*Jt)
v = pinvJt*n
m = sqrt(inner(v,v))
if m < eps()
return zero(v)
else
return v/m
end
end

function get_symbol!(index,val::typeof(map_unit_normal),name="";prefix=index.data.prefix)
:(GalerkinToolkit.map_unit_normal)
end

#discrete_function_term(args...) = DiscreteFunctionTerm(args...)

#struct DiscreteFunctionTerm{A,B} <: AbstractTerm
Expand Down
Loading

0 comments on commit 6b072bc

Please sign in to comment.