Interpolation of symmetric tensors

[gslndlb]

I would like to use symmetric tensors in Gridap.

using Gridap, Gridap.TensorValues

𝒯 = CartesianDiscreteModel((0,2,0,1),(40,20))
Ω = Interior(𝒯)
dΩ = Measure(Ω,2)

Γ = Boundary(𝒯)
dΓ = Measure(Γ, 2)

refₙ = ReferenceFE(lagrangian, SymTensorValue{2, Float64, 3}, 1)
testₙ = TestFESpace(𝒯, refₙ)
N = TransientTrialFESpace(testₙ)

xx = 0:0.05:1.95
yy = 0:0.05:.95

qxx = SymTensorValue(1/2., 0., -1/2.)
qi = [qxx for i = 1:length(xx), j = 1:length(yy)]
q₀ = interpolate_everywhere(qi, N(0.0))

It fails at the definition of the FESpace, so I followed an advice found on the gitter and pasted all functions from here at the beginning of my code:

Now, the definition of the FESpace works but it fails at interpolation with
ERROR: dof ids either positive or negative, not zero

The Julia file to reproduce the bug is quite long so I attach it here (as .txt because Github doesn't accept .jl)
minimal_SymTensorValue.jl.txt

[JordiManyer]

@gslndlb The error is quite indicative. I believe you are setting some of your DoFs with id 0, which creates an error. In Gridap, free dofs have positive (1,2,3...) ids, while constrained dofs (by dirichlet bcs, for instance) have negative ids (-1,-2,...). You cannot have a DoF with zero id.

If you run your code, you get that

julia> get_cell_dof_ids(testₙ)

800-element Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}:
 [1, 4, 124, 127, 3, 6, 126, 129, 0, 0, 0, 0]
 [4, 7, 127, 130, 6, 9, 129, 132, 0, 0, 0, 0]
 [7, 10, 130, 133, 9, 12, 132, 135, 0, 0, 0, 0]
 [10, 13, 133, 136, 12, 15, 135, 138, 0, 0, 0, 0]
 [13, 16, 136, 139, 15, 18, 138, 141, 0, 0, 0, 0]
 [16, 19, 139, 142, 18, 21, 141, 144, 0, 0, 0, 0]
 [19, 22, 142, 145, 21, 24, 144, 147, 0, 0, 0, 0]
 ⋮
 [2440, 2443, 2563, 2566, 2442, 2445, 2565, 2568, 0, 0, 0, 0]
 [2443, 2446, 2566, 2569, 2445, 2448, 2568, 2571, 0, 0, 0, 0]
 [2446, 2449, 2569, 2572, 2448, 2451, 2571, 2574, 0, 0, 0, 0]
 [2449, 2452, 2572, 2575, 2451, 2454, 2574, 2577, 0, 0, 0, 0]
 [2452, 2455, 2575, 2578, 2454, 2457, 2577, 2580, 0, 0, 0, 0]
 [2455, 2458, 2578, 2581, 2457, 2460, 2580, 2583, 0, 0, 0, 0]

So you have an error while filling the dof ids in the ReferenceFE. You are leaving 4 unfilled.

If you create your ReferenceFE manually, you get:

using Gridap.ReferenceFEs

reffe = LagrangianRefFE(SymTensorValue{2, Float64, 3},QUAD,1)

julia> get_face_own_dofs(reffe)
9-element Vector{Vector{Int64}}:
 [1, 5, 5]
 [2, 6, 6]
 [3, 7, 7]
 [4, 8, 8]
 []
 []
 []
 []
 []

As you can see, although you have 12 dofs you are repeating some of them (don't worry about the empty lists).

Digging a little more, one can see the error comes from this function from ReferenceFEs/LagrangianRefFEs.jl:

function _generate_face_own_dofs(face_own_nodes, node_and_comp_to_dof)
  faces = 1:length(face_own_nodes)
  T = eltype(node_and_comp_to_dof)
  comps = 1:num_components(T)
  face_own_dofs = [Int[] for i in faces]
  for face in faces
    nodes = face_own_nodes[face]
    # Node major
    for comp in comps
      for node in nodes
        comp_to_dofs = node_and_comp_to_dof[node]
        dof = comp_to_dofs[comp]
        push!(face_own_dofs[face],dof)
      end
    end
  end

  face_own_dofs
end

In your case, num_components(T) = 3 and comp_to_dofs will be a SymmetricTensorvalue{2,Int,3}. Then when accessing dof = comp_to_dofs[comp] for comp=1,2,3 we will accessing twice the off-diagonal component (since the linear access order is t11,t12,t21,t22). What you want to do is access the components 1,2,4. So you'll have to overload that function as well.

This said, I don't really know why you would want to do this. Why not use a VectorValue{3} if you want 3 components?

[gslndlb]

@JordiManyer thanks for your reply.
The reason I want to do this is because I'm working on equations from nematic liquid crystals, in which the unknown is the nematic tensor. It is a field of a second-order tensor, symmetric and traceless. So I could use 3-components vectors and overwrite standard tensor functions (trace, contraction and divergence) but I thought symmetric tensors were implemented so it would have saved me some time. A workaround is to use non-symmetric tensors, at the cost of performance (not such a problem in 2D because it only adds one component, but in 3D it adds 3 components...).

So if I understood well, the problem is that I need comp_to_dof[3] to be t21 in order other standard tensor-related functions to work (Trace, contractions...), but when building the Finite element space, it is expected that every component corresponds to an actual dof, which t21 is not in the case of symmetric tensors.

I am discovering Gridap, Julia and Finite elements at the same time, so I'm not sure I have the ability to do it. How complicated do you think it is? In particular I'm more used to object-oriented languages (Python and C++) so I'm a bit lost in the structure of the code, to understand which function applies to what.

[amartinhuertas]

@gslndlb ... I have succesfully used symmetric tensor fields in Gridap in the past.

However, I faced a set of issues while working with the current implementation of symmetric tensors in Gridap's main branch; not sure if the issues I found are equivalent as the ones you found. In any case, I gathered all fixes in the following file which is yet to be merged with Gridap:

https://github.com/gridap/GridapHybrid.jl/blob/main/src/tmp_fixes_symmetric_valued_lagrangian_reffes.jl

just in case you might want to give these fixes a try ...

[gslndlb]

@amartinhuertas Yes I included your functions in my Julia code (attached .jl file). I didn't include them in my post because it creates a very long code.

Without your functions I cannot build the FESpace. When I add your functions I can build the space but not interpolate.

[amartinhuertas]

@amartinhuertas Yes I included your functions in my Julia code (attached .jl file). I didn't include them in my post because it creates a very long code.

Ok, I overlooked your first post, sorry ...

[amartinhuertas]

just in case it helps, the driver program that leverages symmetric tensors i available here:

https://github.com/gridap/GridapHybrid.jl/blob/main/test/StressAssistedHenckyHDGTests.jl

[JordiManyer]

So if I understood well, the problem is that I need comp_to_dof[3] to be t21 in order other standard tensor-related functions to work (Trace, contractions...), but when building the Finite element space, it is expected that every component corresponds to an actual dof, which t21 is not in the case of symmetric tensors.

Yes, I believe that is exactly the problem. And conceptually it makes sense that every component corresponds to a single DoF. That is why in most of the fixes we are simply replacing comp_to_dofs[comp] by comp_to_dofs.data[comp], which accesses the three components inside the tensor.

I think that for this to work perfectly and without quirks, we should somehow add an additional interface that is capable of distinguish between the TensorValue (4 components) and the real info stored inside (3 dofs). But that might require some time to think trough, make sure we are not missing anything.

In the meantime (and going back to my initial question) I think you might find it easier to do something like this:

function Π(f::CellField)
  return Π_op∘f
end

function Π_op(x::VectorValue{3})
  return SymTensorValue(x.data)
end

𝒯 = CartesianDiscreteModel((0,2,0,1),(40,20))
Ω = Interior(𝒯)
dΩ = Measure(Ω,2)

Γ = Boundary(𝒯)
dΓ = Measure(Γ, 2)

refₙ = ReferenceFE(lagrangian, VectorValue{3,Float64}, 1)
testₙ = TestFESpace(𝒯, refₙ)
N = TransientTrialFESpace(testₙ)

xx = 0:0.05:1.95
yy = 0:0.05:.95
qxx = VectorValue(1/2., 0., -1/2.)
qi = [qxx for i = 1:length(xx), j = 1:length(yy)]
q₀ = interpolate_everywhere(qi, N(0.0))


a(u,v) = ∫(Π(u) ⊙ Π(v))*dΩ
A = assemble_matrix(a, N, N)

Conceptually, this is what I'm doing:
You work with 3-component vector-valued spaces (which gives you the number of dofs that you want). This spaces are mature and will give you everything you want.
On top of that, you add a "projection" function Π that takes your 3-component vectors and produces SymTensorValues out of them. Conceptually, you mapping your Vector-valued space into a Symmetric-Tensor-Valued space (which has the same dimension, so it's mathematically sound). This way, you can work with Π(u) and do everything you wanted to do. Gridap will handle the conversion behind the scenes. Minimal code, all the functionality.

Some things might be a little tricky to setup at first if you don't know what you are doing, but you should be fine. I can give you support (like how to setup your weak form, etc...) if you need anything specific and are happy to share your code/needs (either here, slack, etc...). For instance, I believe the gradients are not 100% straightforward.

What do you think @amartinhuertas, am I missing anything important?

[JordiManyer]

Just for reference, the fix needed for Gridap.ReferenceFEs._generate_face_own_dofs was

function Gridap.ReferenceFEs._generate_face_own_dofs(
  face_own_nodes, 
  node_and_comp_to_dof::AbstractVector{<:Gridap.TensorValues.SymTensorValue}
  )
  faces = 1:length(face_own_nodes)
  T = eltype(node_and_comp_to_dof)
  comps = 1:num_components(T)
  face_own_dofs = [Int[] for i in faces]
  for face in faces
    nodes = face_own_nodes[face]
    # Node major
    for comp in comps
      for node in nodes
        comp_to_dofs = node_and_comp_to_dof[node]
        dof = comp_to_dofs.data[comp]
        push!(face_own_dofs[face],dof)
      end
    end
  end

  face_own_dofs
end

It does fix some things, but there are still more places inside the code where you would have to do a similar thing.... For instance, get_cell_dof_ids still gives the wrong answer...

[gslndlb]

Thank you for your answers! I followed the idea of working with Vectors instead of tensors and to use a projection.
I used a slightly different approach which is to decompose on a basis of 2x2 tensor that respect the constraints of my unknown tensor (symmetric and traceless). My new unkowns are the coefficients of the decomposition on this basis. This is more flexible than using SymTensorValue, and removes 2 components (one for the tensor because of the traceless constraint and removes the Lagrange multiplier that used to impose the traceless constraint).

If the basis is orthonormal, most terms in the weak form can be expressed even without projection. Some specific terms need a projection.

I think it's generally a good idea to use a relevant basis rather than implementing all possible constraints (symmetric, traceless and others that one could think of). It will be tedious to express the basis in dimension d >= 4 but this seems a niche situation.

[JordiManyer]

@gslndlb I'm glad this solution works for you. Shoot me a message on Gridap's slack/gitter is you need anything along the way. I will be closing this issue now.