fverdugo
commented
3 years ago Hi @Omega-xyZac Can you share the code that generates this?
Closed zjwegert closed 3 years ago
fverdugo
commented
3 years ago Hi @Omega-xyZac Can you share the code that generates this?
fverdugo
commented
3 years ago It should be possible to have arbitrarily complex operations at gauss point level without blowing up the complexity of the type signature of the lazy array.
zjwegert
commented
3 years ago Hi @fverdugo, thank you for the fast reply. I have a very large code base so sharing isn't really an option. However, that example is generated using the below:
....
A = (ε((1/k_0)*data_1.uh) + δⁱʲ(1,1)) ⊗ δⁱʲ(1,1)+
(ε((1/k_0)*data_2.uh) + δⁱʲ(2,2)) ⊗ δⁱʲ(2,2)+
(ε((1/k_0)*data_3.uh) + δⁱʲ(3,3)) ⊗ δⁱʲ(3,3)+
(ε((1/k_0)*data_4.uh) + 1/2*δⁱʲ(1,2)) ⊗ δⁱʲ(1,2)+
(ε((1/k_0)*data_5.uh) + 1/2*δⁱʲ(1,3)) ⊗ δⁱʲ(1,3)+
(ε((1/k_0)*data_6.uh) + 1/2*δⁱʲ(2,3)) ⊗ δⁱʲ(2,3)
ℬ = ε((1/k_0)*data_7.uh) ⊗ eᵢ(1) +
ε((1/k_0)*data_8.uh) ⊗ eᵢ(2) +
ε((1/k_0)*data_9.uh) ⊗ eᵢ(3)
𝒢 = (-∇((1/α_0)*data_1.ϕh)) ⊗ δⁱʲ(1,1)+
(-∇((1/α_0)*data_2.ϕh)) ⊗ δⁱʲ(2,2)+
(-∇((1/α_0)*data_3.ϕh)) ⊗ δⁱʲ(3,3)+
(-∇((1/α_0)*data_4.ϕh)) ⊗ δⁱʲ(1,2)+
(-∇((1/α_0)*data_5.ϕh)) ⊗ δⁱʲ(1,3)+
(-∇((1/α_0)*data_6.ϕh)) ⊗ δⁱʲ(2,3)
H = (-∇((1/α_0)*data_7.ϕh) + eᵢ(1)) ⊗ eᵢ(1) +
(-∇((1/α_0)*data_8.ϕh) + eᵢ(2)) ⊗ eᵢ(2) +
(-∇((1/α_0)*data_9.ϕh) + eᵢ(3)) ⊗ eᵢ(3)
# Base Material
C_b = material.C; # (4-tensor)
e_b = material.e; # (3-tensor)
Κ_b = material.Κ; # (2-tensor)
nely,nelx,nelz = size(struc,1),size(struc,2),size(struc,3);
density = data_1.density_vector
Pᵤᵤ,Pᵤᵩ,Pᵩᵩ=penal; Pᵩᵤ = Pᵤᵩ;
C = density.^Pᵤᵤ.*fill(C_b,nely*nelx*nelz)
e = density.^Pᵤᵩ.*fill(e_b,nely*nelx*nelz)
Κ = density.^Pᵩᵩ.*fill(Κ_b,nely*nelx*nelz)
# Averaging
dΩ = data_1.measure
@GTensor CA[i,j,k,l] = C[i,j,m,p]*A[m,p,k,l]; @GTensor e𝒢[i,j,k,l] = e[m,i,j]*𝒢[m,k,l]
Cᵉ = get_array(∫(CA - e𝒢)dΩ);
....where data_1,data_2,... are of data type
struct SolverFEData
loading::Int8
vol::Float64
uh::SingleFieldFEFunction
ϕh::SingleFieldFEFunction
X::Array
Uϕ::FESpace
trian::CartesianGrid
measure::Measure
density_vector::Any
material_vals::Any
endand the @GTensor macro is the script that I've mentioned in a previous issue for doing general contractions by leveraging the TensorOperations package. Though, these particular operations could be replaced by something like .⋅² and .⋅.
Hopefully that is decipherable..
fverdugo
commented
3 years ago I cannot decipher the code, but it seems that you are implementing complex operations in a not so efficient way for the gridap backend.
You need to implement the complex operation in a function f that takes values at a generic integration point, and then use Operation(f)(...) or f\circ(....) just like in the linear elasticity or the hyperelasticity tutorial
zjwegert
commented
3 years ago Thank you for the suggestion! I think I've found a faster way to do this by integrating the quantities A, ℬ, etc. prior to contractions. I will also try your suggestion. I'll close this issue for now :).
zjwegert
commented
3 years ago Hi again @fverdugo, I tried your suggestion and implemented some functions that take values at generic integration points:
Af(εuh...) = (((1/k_0)*εuh[1]) + δⁱʲ(1,1)) ⊗ δⁱʲ(1,1)+
(((1/k_0)*εuh[2]) + δⁱʲ(2,2)) ⊗ δⁱʲ(2,2)+
(((1/k_0)*εuh[3]) + δⁱʲ(3,3)) ⊗ δⁱʲ(3,3)+
(((1/k_0)*εuh[4]) + 1/2*δⁱʲ(1,2)) ⊗ δⁱʲ(1,2)+
(((1/k_0)*εuh[5]) + 1/2*δⁱʲ(1,3)) ⊗ δⁱʲ(1,3)+
(((1/k_0)*εuh[6]) + 1/2*δⁱʲ(2,3)) ⊗ δⁱʲ(2,3);
ℬf(εuh...) = ((1/k_0)*εuh[1]) ⊗ eᵢ(1) +
((1/k_0)*εuh[2]) ⊗ eᵢ(2) +
((1/k_0)*εuh[3]) ⊗ eᵢ(3)
𝒢f(∇ϕh...) = (-((1/α_0)*∇ϕh[1])) ⊗ δⁱʲ(1,1)+
(-((1/α_0)*∇ϕh[2])) ⊗ δⁱʲ(2,2)+
(-((1/α_0)*∇ϕh[3])) ⊗ δⁱʲ(3,3)+
(-((1/α_0)*∇ϕh[4])) ⊗ δⁱʲ(1,2)+
(-((1/α_0)*∇ϕh[5])) ⊗ δⁱʲ(1,3)+
(-((1/α_0)*∇ϕh[6])) ⊗ δⁱʲ(2,3)
Hf(∇ϕh...) = (-((1/α_0)*∇ϕh[1]) + eᵢ(1)) ⊗ eᵢ(1) +
(-((1/α_0)*∇ϕh[2]) + eᵢ(2)) ⊗ eᵢ(2) +
(-((1/α_0)*∇ϕh[3]) + eᵢ(3)) ⊗ eᵢ(3)
A = Af ∘ (ε(data_1.uh),ε(data_2.uh),ε(data_3.uh),ε(data_4.uh),ε(data_5.uh),ε(data_6.uh))
ℬ = ℬf ∘ (ε(data_7.uh),ε(data_8.uh),ε(data_9.uh))
𝒢 = 𝒢f ∘ (∇(data_1.ϕh),∇(data_2.ϕh),∇(data_3.ϕh),∇(data_4.ϕh),∇(data_5.ϕh),∇(data_6.ϕh))
H = Hf ∘ (∇(data_7.ϕh),∇(data_8.ϕh),∇(data_9.ϕh))
# Base Material
C_b = material.C; e_b = material.e; Κ_b = material.Κ;
nely,nelx,nelz = size(struc,1),size(struc,2),size(struc,3);
density = data_1.density_vector
Pᵤᵤ,Pᵤᵩ,Pᵩᵩ=penal; Pᵩᵤ = Pᵤᵩ;
C = density.^Pᵤᵤ.*fill(C_b,nely*nelx*nelz)
e = density.^Pᵤᵩ.*fill(e_b,nely*nelx*nelz)
Κ = density.^Pᵩᵩ.*fill(Κ_b,nely*nelx*nelz)
# Averaging
dΩ = data_1.measure
@GTensor CA[i,j,k,l] = C[i,j,m,p]*A[m,p,k,l]; @GTensor e𝒢[i,j,k,l] = e[m,i,j]*𝒢[m,k,l]
Cᵉ = ∫(CA - e𝒢)dΩ; sum_C = sum(Cᵉ)
@GTensor eA[i,j,k] = e[i,l,m]*A[l,m,j,k]; @GTensor Κ𝒢[i,j,k] = Κ[i,l]*𝒢[l,j,k]
eᵉ = ∫(eA + Κ𝒢)dΩ; sum_e = sum(eᵉ)
@GTensor eℬ[i,j] = e[i,l,m]*ℬ[l,m,j]; @GTensor ΚH[i,j] = Κ[i,l]*H[l,j]
Κᵉ = ∫(eℬ + ΚH)dΩ; sum_Κ = sum(Κᵉ)Is this along the lines of what you meant? Do you have any other ideas for how to speed this up? Appreciate any suggestions you might have.
Suppose we are contracting, operating then integrating on tensors generated from an FE function. Then you can encounter very complex LazyArrays such as the example below. This is fine unless you want to sum up the array. Summing an array such as the one below can take almost 60 seconds - not a surprise given the typing.
Is there a faster way we could sum such an array? Maybe not given the nature of LazyArrays but I thought it worth the question.