gridap / Gridap.jl

Grid-based approximation of partial differential equations in Julia
MIT License
711 stars 97 forks source link

Integration of ThirdOrderTensorValue Bug #416

Closed zjwegert closed 4 years ago

zjwegert commented 4 years ago

Hi all, Please help! There seems to be a a bug with integrating a 3-tensor. See below:

domain = (0.0, 1, 0.0, 1, 0.0, 1);
model = CartesianDiscreteModel(domain,(20,20,20));
trian = Triangulation(model)
quad = CellQuadrature(trian,2)

ThirdOrder = ThirdOrderTensorValue(1:27 ...)
FourthOrder = SymFourthOrderTensorValue(1:36 ...)

integrate(FourthOrder,trian,quad) ## Working!
integrate(ThirdOrder,trian,quad) ## Error

Throws:

AssertionError: L == D1 * D2 * D3
in include_string at base\loading.jl:1088
in top-level scope at test.jl:9
in integrate at Gridap\y0NoW\src\Geometry\Triangulations.jl:273
in integrate at Gridap\y0NoW\src\CellData\CellQuadratures.jl:108
in integrate at Gridap\y0NoW\src\Fields\Integrate.jl:17
in evaluate_field_array at Gridap\y0NoW\src\Fields\FieldArrays.jl:13
in _evaluate_field_array at Gridap\y0NoW\src\Fields\FieldArrays.jl:18
in apply at Gridap\y0NoW\src\Arrays\Apply.jl:61
in apply at Gridap\y0NoW\src\Arrays\Apply.jl:306
in apply_kernel at Gridap\y0NoW\src\Arrays\Kernels.jl:101
in kernel_cache at Gridap\y0NoW\src\Fields\FieldArrays.jl:24
in field_cache at Gridap\y0NoW\src\Fields\ConstantFields.jl:4
in zeros at base\array.jl:521
in zeros at base\array.jl:526
in zero at base\number.jl:242
in convert at base\number.jl:7
in ThirdOrderTensorValue{3,3,3,Int64,27} at Gridap\y0NoW\src\TensorValues\ThirdOrderTensorValueTypes.jl:48
in ThirdOrderTensorValue{3,3,3,Int64,L} where L at Gridap\y0NoW\src\TensorValues\ThirdOrderTensorValueTypes.jl:8
zjwegert commented 4 years ago

I added some display statements

struct ThirdOrderTensorValue{D1,D2,D3,T,L} <: MultiValue{Tuple{D1,D2,D3},T,3,L}
    data::NTuple{L,T}
    function ThirdOrderTensorValue{D1,D2,D3,T}(data::NTuple{L,T}) where {D1,D2,D3,T,L}
        display(data)
        display([L,D1,D2,D3])
        @assert L == D1*D2*D3
        new{D1,D2,D3,T,L}(data)
    end
end

and received the output for writing ThirdOrder

(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27)
4-element Array{Int64,1}:
 27
  3
  3

and when running integrate:

(0,)
4-element Array{Int64,1}:
 1
 3
 3
 3
zjwegert commented 4 years ago

I think I've found the issue, ThirdOrderTensorValue did not yet have a zero method. I've implemented the below. Will make a pull soon.

@generated function zero(::Type{<:ThirdOrderTensorValue{D1,D2,D3,T}}) where {D1,D2,D3,T}
  L=D1*D2*D3
  quote
    ThirdOrderTensorValue{D1,D2,D3,T}(tfill(Base.zero(T),Val{$L}()))
  end
end
zero(::Type{<:ThirdOrderTensorValue{D1,D2,D3,T,L}}) where {D1,D2,D3,T,L} = ThirdOrderTensorValue{D1,D2,D3,T}(tfill(Base.zero(T),Val{L}()))
zero(::ThirdOrderTensorValue{D1,D2,D3,T,L}) where {D1,D2,D3,T,L} = zero(ThirdOrderTensorValue{D1,D2,D3,T,L})

Let me know what you think.