Add apply_analytical! function

[KnutAM]

Seems like more users may find this convenient, e.g. #434
In contrast to the solution in #434, this solution also supports different function and geometric interpolations.

Remaining

• Tests
• Add to transient heat flow example

[codecov-commenter]

Codecov Report

Base: 92.81% // Head: 93.04% // Increases project coverage by +0.23% 🎉

Coverage data is based on head (da7e532) compared to base (b1a8a30).
Patch coverage: 100.00% of modified lines in pull request are covered.

❗ Current head da7e532 differs from pull request most recent head 4f85093. Consider uploading reports for the commit 4f85093 to get more accurate results

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #532      +/-   ##
==========================================
+ Coverage   92.81%   93.04%   +0.23%     
==========================================
  Files          24       25       +1     
  Lines        3978     4039      +61     
==========================================
+ Hits         3692     3758      +66     
+ Misses        286      281       -5     

Impacted Files	Coverage Δ
src/FEValues/cell_values.jl	100.00% <ø> (ø)
src/Dofs/MixedDofHandler.jl	87.44% <100.00%> (+0.11%)	⬆️
src/Dofs/apply_analytical.jl	100.00% <100.00%> (ø)
src/Quadrature/quadrature.jl	94.54% <100.00%> (+0.31%)	⬆️
src/interpolations.jl	91.48% <0.00%> (+1.48%)	⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[fredrikekre]

I think a more general name would be better since this is not only useful for initial conditions. For example in the homogenization example:

Ferrite.jl/docs/src/literate/computational_homogenization.jl

Lines 519 to 522 in 8e447e2

	chM = ConstraintHandler(dh)
	add!(chM, Dirichlet(:u, Set(1:getnnodes(grid)), (x, t) -> εᴹ[Int(t)] ⋅ x, [1, 2]))
	close!(chM)
	uM = zeros(ndofs(dh))

This seems to be simply to transfer an analytical function to FE space (without doing the L2 projection or solving the pertinent t = 0 equation). I don't have a good suggestion though.

[KnutAM]

Good point! apply_analytical! then perhaps, follows in line with the apply! functions already existing for modifying the degree of freedom vector?

[termi-official]

This seems to only work correctly for pure ODE systems. Algebraic constraints in e.g. DAEs may be a problem, especially enforcing the correct compatibilities between the fields. We should at least document this if you want to merge this.

[KnutAM]

Do you mean that it is difficult for the user to define the correct functions f(x) in the case of a complicated DAE? In that case, isn't it clear that the user must supply a valid analytical solution? If I misunderstand you here, do you have a concrete example?

I can write a short section about using this for initial conditions, and add a warning there about the compatibility. Do you have a suggestion for an appropriate section?

[termi-official]

Indeed, that is my point - writing down these equations for the algebraic portion can be super painful. Usually it is simpler to just provide the solution to the ODE portion of the problem and let the initializer figure out what the compatible solution to the algebraic equations are. At least that is how it is usually handled by DAE solvers. Not sure what the best approach would be here, though.
Also as a side note: It happens quite often that users are not aware that they are solving DAEs and not ODEs.

[koehlerson]

Maybe worth it to compress the gif a little if you haven't already since it doubled in size

[KnutAM]

Maybe worth it to compress the gif a little if you haven't already since it doubled in size

I tried autocompression (from 30MB->7MB), now I did manually - doesn't look so nice, but at least much smaller.
@fredrikekre : Should I add it (and the svg) to gh-pages/assets/ instead?

[termi-official]

Ah sorry forgot to mention. We usually utilize gifsky for compressing the gifs. From our experience the file size is small and visual is still great.

[fredrikekre]

deal.ii: https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#a761f008bdeb7d94a69205ae824deefad

[fredrikekre]

Left some comments :)

[fredrikekre]

Feels like this exist somewhere, maybe BCValues?

[KnutAM]

Yes,

Ferrite.jl/src/FEValues/face_values.jl

Line 264 in 2f86559

M[i, qp, face] = value(geom_interpol, i, ξ)

and

Ferrite.jl/src/FEValues/face_values.jl

Line 278 in 2f86559

x += bcv.M[i,q_point,face] * xh[i] # geometric_value(fe_v, q_point, i) * xh[i]

but this seems specialized for faces? In that case, we would call f(x) even more times?
(I can create a similar construct pre-caching the geometric shape values at the reference points if desired, but I feel the code will become complicated then. )

[fredrikekre]

I see. Perhaps one could create a CellValues with a quadrature rule that has the reference nodes as quadrature points then.

[KnutAM]

Introduced CellDofValues now.
I'm not sure about the name (as usual), but it is internal.

[fredrikekre]

LGTM, but looks like CellDofValues could be removed?

[fredrikekre]

Seems like this could simply be CellScalarValues with a custom quadrature rule with points in the reference coordinates?

[KnutAM]

OK, I'll look into that when I have time.
I guess this should be a form of 1st order lobatto quadrature then. But since this is not implemented for RefTet I have to think a bit how to fix that.
Will ping once fixed!

[fredrikekre]

You can just plug whatever coordinates you want in a "fake" quadrature rule, and think of them as "points for evaluation" instead of quadrature points, for example

julia> using Ferrite

julia> ip = Lagrange{2,RefTetrahedron,1}();

julia> points_for_evaluation = Ferrite.reference_coordinates(ip);

julia> fake_quadrature_rule = QuadratureRule{2,RefTetrahedron,Float64}(fill(NaN, 3), points_for_evaluation);

julia> cellvalues = CellScalarValues(fake_quadrature_rule, ip)
CellScalarValues{2, Float64, RefTetrahedron} with 3 shape functions and 3 quadrature points

You can add a constructor like

function QuadratureRule{D,R}(weights::AbstractVector{T}, points::AbstractVector{Vec{D,T}}) where {D, R, T}
    QuadratureRule{D,R,T}(weights, points)
end

to make it easier (no need to specify T) (see also #513).

[KnutAM]

Thanks for those tips - now I also learned how to do custom quad rules!
I reduced away the dim input in the constructor as well - hope that is ok!

Ferrite.jl/src/Quadrature/quadrature.jl

Lines 43 to 46 in 4f85093

	function QuadratureRule{shape}(weights::AbstractVector{Tw}, points::AbstractVector{Vec{dim,Tp}}) where {dim, shape, Tw, Tp}
	T = promote_type(Tw, Tp)
	QuadratureRule{dim,shape,T}(weights, points)
	end