Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

2D flux field visualization #51

Merged
merged 33 commits into from
Feb 14, 2023
Merged

2D flux field visualization #51

merged 33 commits into from
Feb 14, 2023

Conversation

termi-official
Copy link
Member

@termi-official termi-official commented Dec 13, 2022

using FerriteViz
import GLMakie

include("docs/src/ferrite-examples/heat-equation.jl");
dh,u=manufactured_heat_problem(Triangle, Lagrange{2,RefTetrahedron,2}(), 21)
(dh_grad, u_grad) = FerriteViz.interpolate_gradient_field(dh, u, :u);

FerriteViz.surface(dh_grad, u_grad)

L2-gradient-field

Stepped up the game and checked analytically the gradient field. It converges. :)

Copy link
Member

@koehlerson koehlerson left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nice, I think it would be nice to support vector valued fields in 2D for this.

Can you elaborate on https://github.com/Ferrite-FEM/FerriteViz.jl/blob/master/src/utils.jl#L292-L297 ?

src/utils.jl Show resolved Hide resolved
src/utils.jl Outdated Show resolved Hide resolved
src/utils.jl Outdated Show resolved Hide resolved
src/utils.jl Outdated Show resolved Hide resolved
@termi-official
Copy link
Member Author

The problem mentioned in https://github.com/Ferrite-FEM/FerriteViz.jl/blob/master/src/utils.jl#L292-L297 is that for volumetric entities we need an interpolation on entities with dimension 2 (for the most common problems). Now, L2 spaces have support on the local portion of a face, but not on the actual shared face. However, the current implementation just extract the faces from the provided interpolation to build the face triangles - which is an empty set for L2 spaces.

@termi-official termi-official marked this pull request as ready for review December 15, 2022 13:22
@termi-official
Copy link
Member Author

termi-official commented Dec 15, 2022

Apart from test coverage for interpolate_gradient_field this one should be fine.

Edit: 3D is future work when clip planes are implemented.

@koehlerson
Copy link
Member

Edit: 3D is future work when clip planes are implemented.

How does the docs example with 3D plasticity work then? :D

@termi-official
Copy link
Member Author

Derp.

@koehlerson
Copy link
Member

koehlerson commented Dec 15, 2022

take the incompressible elasticity, there we can visualize both, stresses and strains, or alternatively heat example with the heat flux (so negative gradient with some constant) to show off how to get from gradient to flux

Copy link
Member

@koehlerson koehlerson left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

few suggestions

src/utils.jl Outdated Show resolved Hide resolved
src/utils.jl Outdated Show resolved Hide resolved
src/utils.jl Show resolved Hide resolved
@termi-official
Copy link
Member Author

Okay, I think now only the next Ferrite release is missing (because of the L2 fix you pushed @koehlerson ).

@koehlerson
Copy link
Member

Okay, I think now only the next Ferrite release is missing (because of the L2 fix you pushed @koehlerson ).

builds locally with linear ansatz for the displacements 👍

@koehlerson
Copy link
Member

we can also think about having Ferrite#master as docs dependency

@termi-official
Copy link
Member Author

I would vote for sticking with the safe route, manually upgrading at each release to capture regressions.

@termi-official
Copy link
Member Author

termi-official commented Dec 15, 2022

Gradient field converges for the heat problem towards the analytical gradient in 2d

julia> maximum(u_grad) - π/2
0.0072785284134289086

and 3d

julia> maximum(u_grad) - π/2
0.01160556577741545

@termi-official termi-official changed the title Prototype for flux field visualization 2D flux field visualization Dec 15, 2022
@koehlerson
Copy link
Member

I think we could add a Union in the field argument such that either one gradient or multiple ones are computed. We should also provide the option to add a deformationfield to the dh_grad, the latter will be annoying for setting up the total DoF vector for dh_grad

@termi-official
Copy link
Member Author

Before we start introducing more and more hotfixes, I think we should focus on reorganizing the code.

@termi-official
Copy link
Member Author

I will skip the arrows plot in favor of #55 and #54 . Should be good to go (minus getting Ferrite master into test CI).

@termi-official termi-official added the enhancement New feature or request label Jan 13, 2023
Copy link
Member

@koehlerson koehlerson left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Otherwise LGTM

docs/src/atopics.md Outdated Show resolved Hide resolved
CHANGELOG.md Show resolved Hide resolved
src/utils.jl Outdated Show resolved Hide resolved
Dennis Ogiermann and others added 2 commits February 13, 2023 13:13
Co-authored-by: Maximilian Köhler <maximilian.koehler@ruhr-uni-bochum.de>
src/utils.jl Outdated Show resolved Hide resolved
@termi-official
Copy link
Member Author

Correctness of shear component computation has been roughly tested via some crappy simple shear

include("docs/src/ferrite-examples/incompressible-elasticity.jl")

function create_cook_grid(nx, ny)
    grid = generate_grid(Triangle, (nx, ny));
    # facesets for boundary conditions
    addfaceset!(grid, "traction", x -> x[1] ≈ 2.0);
    return grid;
end;

function create_bc(dh)
    dbc = ConstraintHandler(dh)
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "left"), (x,t) -> zero(Tensors.Vec{2}), [1,2]))
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "right"), (x,t) -> Tensors.Vec{2}((0.0, 0.1)), [1,2]))
    close!(dbc)
    t = 0.0
    update!(dbc, t)
    return dbc
end;

u_quadratic,dh_quadratic = solve(quadratic, linear, mp);

(dh_quadratic_grad, u_quadratic_grad) = FerriteViz.interpolate_gradient_field(dh_quadratic, u_quadratic, :u)
plotter_quadratic = FerriteViz.MakiePlotter(dh_quadratic_grad, u_quadratic_grad)
cmap = :jet

f = GLMakie.Figure()
axs = [GLMakie.Axis(f[1, 1], title="u_{x}"),GLMakie.Axis(f[1, 2], title="u_{y}"),
       GLMakie.Axis(f[3, 1], title="u_{x,y}"),GLMakie.Axis(f[3, 2], title="u_{y,x}")]
p1 = FerriteViz.solutionplot!(axs[1], dh_quadratic, u_quadratic, process=u->u[1], colormap=cmap)
p2 = FerriteViz.solutionplot!(axs[2], dh_quadratic, u_quadratic, process=u->u[2], colormap=cmap)
f[2,1] = GLMakie.Colorbar(f[1,1], p1, vertical=false)
f[2,2] = GLMakie.Colorbar(f[1,2], p2, vertical=false)

p4 = FerriteViz.solutionplot!(axs[3], plotter_quadratic, process=∇u->∇u[2], colormap=cmap, colorrange = (-1e-2, 1e-2))
p5 = FerriteViz.solutionplot!(axs[4], plotter_quadratic, process=∇u->∇u[3], colormap=cmap, colorrange = (-0e-2, 6e-2))
f[4,1] = GLMakie.Colorbar(f[3,1], p4, vertical=false)
f[4,2] = GLMakie.Colorbar(f[3,2], p5, vertical=false)

Screenshot from 2023-02-13 18-07-44

@termi-official
Copy link
Member Author

CI failed because it turns out that Tensors.jl gradients are the outermost indices, but access to a Vec (case: scalar problem) with two indices does not work. :) Should be fine now.

@termi-official termi-official merged commit 3fa406e into master Feb 14, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants