Extend quadrature type to include a storage type parameter #991

termi-official · 2024-06-26T12:24:05Z

Resolves #989

src/Quadrature/quadrature.jl

codecov · 2024-06-26T12:40:39Z

Codecov Report

Attention: Patch coverage is 93.87755% with 3 lines in your changes missing coverage. Please review.

Project coverage is 93.20%. Comparing base (3d12883) to head (785010f).

Files	Patch %	Lines
src/Quadrature/quadrature.jl	92.68%	3 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #991      +/-   ##
==========================================
+ Coverage   93.18%   93.20%   +0.01%     
==========================================
  Files          38       38              
  Lines        5766     5766              
==========================================
+ Hits         5373     5374       +1     
+ Misses        393      392       -1

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

KristofferC · 2024-06-26T13:16:33Z

src/Quadrature/quadrature.jl

-    function QuadratureRule{shape, T}(weights::Vector{T}, points::Vector{Vec{dim, T}}) where {dim, shape <: AbstractRefShape{dim}, T}
+# struct QuadratureRule{shape, T, rdim, WeightStorageType <: AbstractVector{T}, PointStorageType <: AbstractVector{<:Vec{rdim,T}}}
+struct QuadratureRule{shape, T, rdim, WeightStorageType <: AbstractVector{T}, PointStorageType <: AbstractVector} # The one above crashes the deprecation in deprecations.jl for some reason...
+    weights::WeightStorageType


Why do they need to be able to have different storage type and different array type? Couldn't the change here only allow for swapping out the Vector part while keeping the rest?

Fredrik just asked me the same :D I do not know how, because the vectors store different things. I may be able to store some compressed variant of the Vec (as a AbstractVector{T} with num weights * dim entries), but not sure if it is worth the effort.

An alternative is to have an AbstractQuadratureRule that needs to support some interface and then if you need really exotic ones you can always BYOQR (Bring Your Own Quadrature Rule). And then the FEValues would just take a QR <: AbstractQuadratureRule.

KristofferC · 2024-06-26T13:19:39Z

Just to see what people think, is removing the shape parameter a possibility? Like, it doesn't really define a quadrature rule. It is there now to try improve the "static" type safety so you don't accidentally forget to change the quadrature rule if you change the element type but is it worth spending a whole type parameter on that?

termi-official · 2024-06-26T13:24:05Z

Just to see what people think, is removing the shape parameter a possibility? Like, it doesn't really define a quadrature rule. It is there now to try improve the "static" type safety so you don't accidentally forget to change the quadrature rule if you change the element type but is it worth spending a whole type parameter on that?

I am really split on this one. On the one hand I find it annoying to move around with the rule, so I cannot plug my QuadratureRule in everywhere. On the other hand I am also happy to have this safety net when I accidentally plug a quadrature rule into a place where it should not be such that I get an early error message.

KnutAM · 2024-06-26T13:36:02Z

Just to see what people think, is removing the shape parameter a possibility

How should the constructor get information about the quadrature points then? QuadratureRule(shape::AbstractRefShape, order::Int)?

I don't really see the problem with having the refshape, but with this change, I think using

struct QuadratureRule{RefShape, WT, PT}
    weights::WT
    points::PT
    function QuadratureRule{RefShape}(weights::AbstractVector{T}, points::AbstractVector{Vec{rdim, T}}) where {rdim, RefShape <: AbstractRefShape{rdim}, T}
        if length(weights) != length(points)
            throw(ArgumentError("number of weights and number of points do not match"))
        end
        new{RefShape, typeof(weights), typeof(points)}(weights, points)
    end
end

should suffice. The inner constructor ensures that WT and PT always are AbstractArrays.

Edit: Removed T from "function QuadratureRule{RefShape, T}"

termi-official · 2024-06-26T13:37:10Z

How should the constructor get information about the quadrature points then? QuadratureRule(shape::AbstractRefShape, order::Int)?

We can leave the constructor as is.

KnutAM · 2024-06-26T13:38:28Z

We can leave the constructor as is.

We could, but as interface it isn't so nice that QuadratureRule{RefTriangle}(2) isa QuadratureRule{RefTriangle} is false

KnutAM · 2024-06-26T13:51:34Z

And with what I wrote above, we should adopt the same API for giving the number type as for FEValues, e.g.,

function QuadratureRule{shape}(::Type{T}, order::Int) where {shape <: AbstractRefShape, T <: Number}
    quad_type = (shape === RefPrism || shape === RefPyramid) ? (:polyquad) : (:legendre)
    return QuadratureRule{shape}(T, quad_type, order)
end

KristofferC · 2024-06-26T14:45:00Z

How should the constructor get information about the quadrature points then?

There would be a separation between the convenience function to create different quadrature rules (from the tables) and the type constructor itself (that would just take points and weights). So you would have something like create_gauss_quadrature(RefTriangle, 3) which creates a QuadratureRule.

KnutAM · 2024-06-26T19:18:52Z

improve the "static" type safety so you don't accidentally forget to change the quadrature rule if you change the element type but is it worth spending a whole type parameter on that?

After thinking about this today, my answer would be yes, I think that is worth spending a type parameter on.

I'm almost 100 % sure that without this parameter, bugs will occur in user codes (including my own) and silently give wrong results.

KnutAM

This would be my suggestion for how to do this

KnutAM · 2024-06-26T19:26:55Z

src/Quadrature/quadrature.jl

+struct FacetQuadratureRule{shape, T, dim, QRType <: QuadratureRule{shape, T, dim}, FacetRulesType <: Union{NTuple{<:Any, QRType}, AbstractVector{QRType}}}
+    face_rules::FacetRulesType


For the Array-based QuadratureRule this only contains references to two arrays, so I think we can simplify to always use a homogeneous tuple. nfaces(RefShape) should give a Const return value if needed to get N.

Suggested change

struct FacetQuadratureRule{shape, T, dim, QRType <: QuadratureRule{shape, T, dim}, FacetRulesType <: Union{NTuple{<:Any, QRType}, AbstractVector{QRType}}}

face_rules::FacetRulesType

struct FacetQuadratureRule{RefShape, QR<:QuadratureRule{RefShape}, N}

face_rules::NTuple{N, QR}

If everyone is okay with this then we can also follow this idea. I mean even tesseract only has 8 facets.

Might as well allow any container type then?

But does everything in Ferrite work if we use a non-indexable container? I mean, the array interface is the supertype for all indexable containers, right?

KnutAM · 2024-06-26T19:27:17Z

src/Quadrature/quadrature.jl

-    points::Vector{Vec{dim,T}}
-    function QuadratureRule{shape, T}(weights::Vector{T}, points::Vector{Vec{dim, T}}) where {dim, shape <: AbstractRefShape{dim}, T}
+# struct QuadratureRule{shape, T, rdim, WeightStorageType <: AbstractVector{T}, PointStorageType <: AbstractVector{<:Vec{rdim,T}}}
+struct QuadratureRule{shape, T, rdim, WeightStorageType <: AbstractVector{T}, PointStorageType <: AbstractVector} # The one above crashes the deprecation in deprecations.jl for some reason...


struct QuadratureRule{RefShape, WT, PT} weights::WT points::PT function QuadratureRule{RefShape}(weights::AbstractVector{T}, points::AbstractVector{Vec{rdim, T}}) where {rdim, RefShape <: AbstractRefShape{rdim}, T} if length(weights) != length(points) throw(ArgumentError("number of weights and number of points do not match")) end new{RefShape, typeof(weights), typeof(points)}(weights, points) end end

And change all (L60-L155)

- function QuadratureRule{ConcreteRefShape, T}(quad_type::Symbol, order::Int) where T + function QuadratureRule{ConcreteRefShape}(::Type{T}, quad_type::Symbol, order::Int) where {T<:Number}

So if I understand your suggestion correctly, then you would remove all type constraints in the type parametrization?

Yes, but these should be checked in the inner constructor.

(I left one in FacetQuadratureRule but that was just since it makes it more readable IMO, but that could also just be checked in the inner constructor of course)

KnutAM · 2024-06-26T19:34:42Z

src/Quadrature/quadrature.jl

+function QuadratureRule{shape}(weights::AbstractVector{T}, points::AbstractVector{Vec{rdim, T}}) where {rdim, shape <: AbstractRefShape{rdim}, T}
    QuadratureRule{shape, T}(weights, points)


- QuadratureRule{shape, T}(weights, points) + QuadratureRule{shape}(weights, points)

lijas · 2024-06-27T07:00:27Z

I think there should be some consideration on code simplicity here. Personally I think trying to have every functionality in one QuadratureRule struct is a bit too much, and would much more prefer to just create GPUQuadratureRule, and use some create_quadrule(...) interface, at least now in the beginning when GPU stuff are in the prototyping stage.

KnutAM · 2024-06-27T07:45:54Z

I think there should be some consideration on code simplicity here.

I agree, we should think carefully about this as FEValues are already quite intricately typed, and as @KristofferC mentions above, the interface is quite minimial for BYOQR. I guess the point is, that a GPU-compatible implementation allowing compatible arrays would also support the current functionality. So then it feels like we are duplicating the code, but perhaps that is worth it?

If we go for separate AbstractQuadratureRules, I would still vote for changing FacetQuadratureRule to contain a tuple, since then at least that can be used directly, and (at least I cannot see that) it introduces any significant difference in complexity.

Furthermore, I believe that the change to have the constructor as QuadratureRule{RefShape}([::Type{T},] [quad_type::Symbol], order::Int) makes sense, since this follows the way its done for CellValues. It will also make the GPUQuadratureRule and QuadratureRule have the same constructor.

termi-official · 2024-06-27T09:06:50Z

I think there should be some consideration on code simplicity here. Personally I think trying to have every functionality in one QuadratureRule struct is a bit too much, and would much more prefer to just create GPUQuadratureRule, and use some create_quadrule(...) interface, at least now in the beginning when GPU stuff are in the prototyping stage.

I agree, we should think carefully about this as FEValues are already quite intricately typed, and as @KristofferC mentions above, the interface is quite minimial for BYOQR. I guess the point is, that a GPU-compatible implementation allowing compatible arrays would also support the current functionality. So then it feels like we are duplicating the code, but perhaps that is worth it?

I agree on the argument that we should keep the code simple and this is exactly the reason why I propose the change. I think it is simpler to have the storage parametrized than maintaining additional >300 line of duplicated code for each possible storage type. I hope that this blog post here can explain the DRY principle better than I can https://scalastic.io/en/solid-dry-kiss/ .

Can we discuss the audit of the CellValues separately?

If we go for separate AbstractQuadratureRules, I would still vote for changing FacetQuadratureRule to contain a tuple, since then at least that can be used directly, and (at least I cannot see that) it introduces any significant difference in complexity.

If there is no disadvantage associated to this change, why not? Should we still provide a constructor which takes an abstract vector tho?

Furthermore, I believe that the change to have the constructor as QuadratureRule{RefShape}([::Type{T},] [quad_type::Symbol], order::Int) makes sense, since this follows the way its done for CellValues. It will also make the GPUQuadratureRule and QuadratureRule have the same constructor.

Since I have been repeatedly asked to make only small changes and one by one: Can we discuss this issue separately? This PR is really just about the storage type and not the quadrature custruction interface.

KnutAM · 2024-06-27T09:43:18Z

Can we discuss the audit of the CellValues separately?

That comment was simply regarding the cognitive load for users, not directly asking for changes to CellValues.

Should we still provide a constructor which takes an abstract vector tho?

Not sure, but this is (in practice) always internal right, so not required to provide such a convenience interface?

Since I have been repeatedly asked to make only small changes and one by one: Can we discuss this issue separately? This PR is really just about the storage type and not the quadrature custruction interface.

If we change the type-parameterization from QuadratureRule{RefShape, T, rdim} to QuadratureRule{RefShape, WT, PT}, then this PR is, IMO, related to the constructor, since it would result in QuadratureRule{RefShape, Float64}(2) isa QuadratureRule{RefShape, Float64} being false.
If the final descision here becomes to e.g. only change FacetQuadratureRules and introduce an AbstractQuadratureRule supertype, then of course the constructor change could be made in a separate PR.

Yet, I find the discussion of how it should look relevant here, since we are discussing how to have QuadratureRules that support GPU-compatible arrays. I think we should strive for an as homogeneous interface as possible, and removing the T parameter as discussed here, does question if having users provide that as a type-parameter instead of an argument is a good interface.

KnutAM · 2024-06-27T09:54:17Z

I think it is simpler to have the storage parametrized than maintaining additional >300 line of duplicated code for each possible storage type.

Perhaps I'm missing something, but a custom QuadratureRule should be quite small?

struct MyQR{RefShape, WT, PT} <: AbstractQuadratureRule{RefShape}
    weights::WT
    points::PT
    function MyQR{RefShape}(weights::AbstractVector{T}, points::AbstractVector{<:Vec{rdim, T}) where {RefShape, T, rdim}
        @assert length(weights) == length(points)
        @assert rdim == getrefdim(RefShape)
        return new{RefShape, typeof(weights), typeof(points)}(weights, points)
    end
end
function MyQR{RefShape}(::Type{StorageType}, args...) where {RefShape, StorageType<:AbstractVector}
    weights, points = get_weights_and_points(RefShape, args...) 
    # convert to StorageType
    return MyQR{RefShape}(converted_weights, converted_points)
end
getweights(qr::MyQR) = qr.weights
getpoints(qr::MyQR) = qr.points

where get_weights_and_points(RefShape, [T::Type{<:Number}], [rule_type::Symbol], order::Int) is factored out from current code bound to QuadratureRule.

(But of course, this implementation would also cover the current usage, so I agree that it is code duplication. I guess the "difficulty" with this general implementation is the less clear parameterization WT and PT over Vector{T} and Vector{Vec{rdim,T}}). Do I understand your concern correctly here @lijas ?

termi-official · 2024-06-27T11:32:48Z

Sorry, I cannot follow what exact your argument Knut and Elias. So let me try to summarize what I understand here: Instead of defining a single type with parametrized storage constrained to collections, we should have an interface, that tells we should define one type per collection type, which weights and points of concrete collections each and change the bulk rename QuadratureRule to AbstractQuadratureRule in most parts of the codebase?

Also

If we change the type-parameterization from QuadratureRule{RefShape, T, rdim} to QuadratureRule{RefShape, WT, PT}, then this PR is, IMO, related to the constructor, since it would result in QuadratureRule{RefShape, Float64}(2) isa QuadratureRule{RefShape, Float64} being false.

That is correct for your proposal. But please note that for the current state of the PR QuadratureRule{RefShape, Float64}(2) isa QuadratureRule{RefShape, Float64} is true.

lijas · 2024-06-27T11:48:10Z

https://scalastic.io/en/solid-dry-kiss/

Just skimming through some of the principles outlined there, it seems like the first three principles (srp, ocp, isp), could be used as arguments against this PR :P

I think maintaining a monolithic QuadratureRule is much more difficult that maintaining a AbstractQuadratureRule Interface. And new users who read the code for the first time, will not have to put a lot of mental effort in to what WeightType and PointType does.

But I dont want to block this PR, I just wanted to offer my view, so if you think it will simplify stuff later on, go for it :) .

KnutAM · 2024-06-27T12:03:25Z

I agree that we can argue for both options being easier or harder, depending on the point of view.

we should have an interface ... and change the bulk rename QuadratureRule to AbstractQuadratureRule in most parts of the codebase?

Yes, which would also allow users to "hack" in using their own quadrature rule (e.g. if I want a quadrature rule with tuple elements, or if you want to skip the weights (as we do for e.g. point values). If this would make things easier to understand: I'm not sure.

that tells we should define one type per collection type, which weights and points of concrete collections each

The interface wouldn't have to say that, just the QuadratureRule-type would fullfill that. At the end of the day, this is the same discussion we had for the CellValues in #764 if we should maintain both a simple variant and a complex variant (that can do the job of the simple one too), and there we decided against it.

That is correct for your proposal. But please note that for the current state of the PR QuadratureRule{RefShape, Float64}(2) isa QuadratureRule{RefShape, Float64} is true.

Yes, but since noone objected to that proposal, I assume that is one possible solution, and then I think we should discuss the implications of that here?

will not have to put a lot of mental effort in to what WeightType and PointType does.

An option could perhaps be just adding the comments

struct QuadratureRule{RefShape, WT, PT}
    weights::WT # E.g. Vector{Float64}
    points::PT  # E.g. Vector{Vec{2, Float64}}
    function QuadratureRule{RefShape}(weights::AbstractVector{T}, points::AbstractVector{Vec{rdim, T}}) 
...

?

I'm open for both solutions, but I do think that in the current state, the parameterization is too complex (hence my review suggestions).

fredrikekre · 2024-06-27T12:11:13Z

I like the definition from Knuts last comment. It is simple but flexible. We can of course also have AbstractQUadrature if the type definition there does not suffice in some case.

termi-official · 2024-06-27T12:53:03Z

https://scalastic.io/en/solid-dry-kiss/

Just skimming through some of the principles outlined there, it seems like the first three principles (srp, ocp, isp), could be used as arguments against this PR :P

No. Why? The QuadratureRule is responsible for managing the access to quadrature points and the weights, hence all proposals adhere to the single responsibility principle. The proposal for the adding a supertype and making the storage parameter in this PR also follow the open-closed principle, whereas the QuadratureRule on master does not, as we cannot extend it (i.e. "not open"). The proposal to leave the type parameters unconstrained violates the Liskov substituion principle, as we can use storage types which might not be ordered and hence the relation between the weights and points might be messed up. But for this one I am not sure if this will really be an issue in practice, as I think users who use such advanced data structures should be aware of this fact?

I think maintaining a monolithic QuadratureRule is much more difficult that maintaining a AbstractQuadratureRule Interface. And new users who read the code for the first time, will not have to put a lot of mental effort in to what WeightType and PointType does.

So, if I understand correctly, then your point is essentially that I should add more documentation and devdocs?

Yes, which would also allow users to "hack" in using their own quadrature rule (e.g. if I want a quadrature rule with tuple elements, or if you want to skip the weights (as we do for e.g. point values).

But if you do not have weights, then it would not be a quadrature rule in first place, right? (But something like QuadraturePointVector or similar).

The interface wouldn't have to say that, just the QuadratureRule-type would fullfill that.

But what would define the interface for an AbstractQuadratureRule then?

An option could perhaps be just adding the comments

struct QuadratureRule{RefShape, WT, PT}
    weights::WT # E.g. Vector{Float64}
    points::PT  # E.g. Vector{Vec{2, Float64}}
    function QuadratureRule{RefShape}(weights::AbstractVector{T}, points::AbstractVector{Vec{rdim, T}}) 
...

If no one else is concerned by the hypothetical correctness issue pointed out above, then I can also do this (as I also do not think that this will be really a practical issue).

fredrikekre · 2024-06-27T13:32:58Z

But what would define the interface for an AbstractQuadratureRule then?

Just getweights and getpoints which return something that iterates weights and points in the same order should be enough I think.

KristofferC · 2024-06-30T10:55:54Z

One things I want to say (which is independent of the parameterization) is that I think the current way of creating quadrature rules is not ideal. For example we write:

qr = QuadratureRule{RefTetrahedron}(2)

but there are many quadrature rules and it isn't clear which one this creates etc. This was changed by me in #58 and 8 years later I am suggesting we change it back to standard functions like:

create_legendre_quadrature(RefTetrahedron::RefShape, 1::Int)
create_lobotto_quadrature(RefCube::RefShape, 1::Int)

etc.

Extend quadrature type to include a storate type parameter

dbb32e7

termi-official commented Jun 26, 2024

View reviewed changes

src/Quadrature/quadrature.jl Outdated Show resolved Hide resolved

Derp

df40354

termi-official added the awaiting review PR is finished from the authors POV, waiting for feedback label Jun 26, 2024

KristofferC reviewed Jun 26, 2024

View reviewed changes

KnutAM reviewed Jun 26, 2024

View reviewed changes

termi-official and others added 3 commits June 28, 2024 13:21

Apply suggestions from code review.

dd4c7b6

Fix

5641411

Merge branch 'master' into do/quadrature-storage-parameter

638359e

fredrikekre changed the title ~~Extend quadrature type to include a storate type parameter~~ Extend quadrature type to include a storage type parameter Jun 30, 2024

fix FacetValues when fqr is a tuple

6307ccc

fredrikekre added 2 commits June 30, 2024 13:00

use vector by default

61cdee6

fix

785010f

fredrikekre merged commit 098d83d into master Jun 30, 2024
11 checks passed

fredrikekre deleted the do/quadrature-storage-parameter branch June 30, 2024 11:16

fredrikekre added this to the v1.0.0 milestone Jun 30, 2024

This was referenced Jun 30, 2024

Use separate functions for creating quadrature rule instead of the QuadratureRule constructor #1001

Open

Use separate functions for creating quadrature rule instead of the QuadratureRule constructor #1002

Closed

		struct FacetQuadratureRule{shape, T, dim, QRType <: QuadratureRule{shape, T, dim}, FacetRulesType <: Union{NTuple{<:Any, QRType}, AbstractVector{QRType}}}
		face_rules::FacetRulesType

		function QuadratureRule{shape}(weights::AbstractVector{T}, points::AbstractVector{Vec{rdim, T}}) where {rdim, shape <: AbstractRefShape{rdim}, T}
		QuadratureRule{shape, T}(weights, points)

Extend quadrature type to include a storage type parameter #991

Extend quadrature type to include a storage type parameter #991

Conversation

termi-official commented Jun 26, 2024

codecov bot commented Jun 26, 2024 • edited Loading

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KristofferC commented Jun 26, 2024

termi-official commented Jun 26, 2024

KnutAM commented Jun 26, 2024 • edited Loading

termi-official commented Jun 26, 2024

KnutAM commented Jun 26, 2024

KnutAM commented Jun 26, 2024 • edited Loading

KristofferC commented Jun 26, 2024 • edited Loading

KnutAM commented Jun 26, 2024

KnutAM left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lijas commented Jun 27, 2024

KnutAM commented Jun 27, 2024

termi-official commented Jun 27, 2024

KnutAM commented Jun 27, 2024

KnutAM commented Jun 27, 2024 • edited Loading

termi-official commented Jun 27, 2024

lijas commented Jun 27, 2024 • edited Loading

KnutAM commented Jun 27, 2024

fredrikekre commented Jun 27, 2024

termi-official commented Jun 27, 2024

fredrikekre commented Jun 27, 2024

KristofferC commented Jun 30, 2024

codecov bot commented Jun 26, 2024 •

edited

Loading

KnutAM commented Jun 26, 2024 •

edited

Loading

KnutAM commented Jun 26, 2024 •

edited

Loading

KristofferC commented Jun 26, 2024 •

edited

Loading

KnutAM commented Jun 27, 2024 •

edited

Loading

lijas commented Jun 27, 2024 •

edited

Loading