mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-26 11:51:31 +00:00
Jukka Aho
141 lines
3.6 KiB
Julia
141 lines
3.6 KiB
Julia
# This file is a part of JuliaFEM.
|
|
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
|
|
|
# There are here for now until I figure a better place for them.
|
|
|
|
|
|
"""
|
|
Fit field s.t. || ∫ (Nᵢ(ξ)αᵢ - f(el, ξ)) dS || -> min!
|
|
|
|
Parameters
|
|
----------
|
|
f::Function
|
|
Needs to take (el::Element, xi::Vector) as argument
|
|
fixed_coeffs::Int[]
|
|
These coefficients are not changed during fitting -> constrained optimizatio
|
|
"""
|
|
function fit_field!(el::Element, field, f, fixed_coeffs=Int[])
|
|
w = [
|
|
128/225,
|
|
(332+13*sqrt(70))/900,
|
|
(332+13*sqrt(70))/900,
|
|
(332-13*sqrt(70))/900,
|
|
(332-13*sqrt(70))/900]
|
|
xi = Vector[
|
|
[0.0],
|
|
[ 1/3*sqrt(5 - 2*sqrt(10/7))],
|
|
[-1/3*sqrt(5 - 2*sqrt(10/7))],
|
|
[ 1/3*sqrt(5 + 2*sqrt(10/7))],
|
|
[-1/3*sqrt(5 + 2*sqrt(10/7))]]
|
|
n = get_number_of_basis_functions(el)
|
|
fld = get_field(el, field)
|
|
nfld = length(fld[1])
|
|
#Logging.debug("dim of field $field: $nfld")
|
|
|
|
M = zeros(n, n)
|
|
b = zeros(n, nfld)
|
|
for i=1:length(w)
|
|
detJ = get_detJ(el, xi[i])
|
|
N = get_basis(el, xi[i])
|
|
M += w[i]*N*N'*detJ
|
|
fi = f(el, xi[i])
|
|
for j=1:nfld
|
|
b[:, j] += w[i]*N*fi[j]*detJ
|
|
end
|
|
end
|
|
|
|
coeffs = zeros(n)
|
|
for j=1:nfld
|
|
for k=1:n
|
|
coeffs[k] = fld[k][j]
|
|
end
|
|
if length(fixed_coeffs) != 0
|
|
# constrained problem, some coefficients are fixed
|
|
N = Int[] # rest of coeffs
|
|
S = Int[] # fixed coeffs
|
|
for i = 1:n
|
|
if i in fixed_coeffs
|
|
push!(S, i)
|
|
else
|
|
push!(N, i)
|
|
end
|
|
end
|
|
lhs = M[N,N]
|
|
rhs = b[N,j] - M[N,S]*coeffs[S]
|
|
coeffs[N] = lhs \ rhs
|
|
else
|
|
coeffs[:] = M \ b[:,j]
|
|
end
|
|
for k=1:n
|
|
fld[k][j] = coeffs[k]
|
|
end
|
|
end
|
|
set_field(el, field, fld)
|
|
return
|
|
end
|
|
|
|
|
|
"""
|
|
Fit field s.t. || ∫ ∂/∂ξ(∑Nᵢ(ξ)αᵢ)f(el, ξ) dS || -> min!
|
|
"""
|
|
function fit_derivative_field!(el::Element, field, f, fixed_coeffs=Int[])
|
|
w = [
|
|
128/225,
|
|
(332+13*sqrt(70))/900,
|
|
(332+13*sqrt(70))/900,
|
|
(332-13*sqrt(70))/900,
|
|
(332-13*sqrt(70))/900]
|
|
xi = Vector[
|
|
[0.0],
|
|
[ 1/3*sqrt(5 - 2*sqrt(10/7))],
|
|
[-1/3*sqrt(5 - 2*sqrt(10/7))],
|
|
[ 1/3*sqrt(5 + 2*sqrt(10/7))],
|
|
[-1/3*sqrt(5 + 2*sqrt(10/7))]]
|
|
n = get_number_of_basis_functions(el)
|
|
fld = get_field(el, field)
|
|
nfld = length(fld[1])
|
|
#Logging.debug("dim of field $field: $nfld")
|
|
|
|
M = zeros(n, n)
|
|
b = zeros(n, nfld)
|
|
for i=1:length(w)
|
|
detJ = get_detJ(el, xi[i])
|
|
dNdxi = get_dbasisdxi(el, xi[i])
|
|
dNdX = dNdxi / detJ
|
|
M += w[i]*dNdX*dNdX'*detJ
|
|
fi = f(el, xi[i])
|
|
for j=1:nfld
|
|
b[:, j] += w[i]*dNdX*fi[j]*detJ
|
|
end
|
|
end
|
|
|
|
coeffs = zeros(n)
|
|
for j=1:nfld
|
|
for k=1:n
|
|
coeffs[k] = fld[k][j]
|
|
end
|
|
if length(fixed_coeffs) != 0
|
|
#Logging.info("constrained problem, some coefficients are fixed")
|
|
N = Int[] # rest of coeffs
|
|
S = Int[] # fixed coeffs
|
|
for i = 1:n
|
|
if i in fixed_coeffs
|
|
push!(S, i)
|
|
else
|
|
push!(N, i)
|
|
end
|
|
end
|
|
lhs = M[N,N]
|
|
rhs = b[N,j] - M[N,S]*coeffs[S]
|
|
coeffs[N] = lhs \ rhs
|
|
else
|
|
coeffs[:] = M \ b[:,j]
|
|
end
|
|
for k=1:n
|
|
fld[k][j] = coeffs[k]
|
|
end
|
|
end
|
|
set_field(el, field, fld)
|
|
return
|
|
end
|