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

Bad Dof-Ordering when using FerriteMeshParser in Ferrite#master #49

Open
Joroks opened this issue Jul 15, 2024 · 3 comments
Open

Bad Dof-Ordering when using FerriteMeshParser in Ferrite#master #49

Joroks opened this issue Jul 15, 2024 · 3 comments

Comments

@Joroks
Copy link

Joroks commented Jul 15, 2024

Currently when using FerriteMeshParser to parse the .inp file I get the following stiffnes matrix:
grafik

If I want to get a somewhat good looking matrix I have to first sort the cellsets in the Grid:

values(grid.cellsets) .|> sort!

I then get
grafik

My guess is that this is because we don't use OrderedCollections like Ferrite#master does, so this should probably be a pretty simple fix?

@Joroks
Copy link
Author

Joroks commented Jul 15, 2024

It seems like this only happens when using SubDofHandlers. I'll quickly set up a minimal example

@Joroks
Copy link
Author

Joroks commented Jul 15, 2024

Here is the file I'm using: mesh.txt

using Ferrite # on current master branch
using FerriteMeshParser
grid = get_ferrite_grid("mesh.txt", meshformat=FerriteMeshParser.AbaqusMeshFormat(), generate_facetsets=false)
ip = Lagrange{RefHexahedron, 1}()

dh1 = DofHandler(grid)
add!(dh1, :u, ip)
close!(dh1)

dh2 = DofHandler(grid)
sdh = SubDofHandler(dh2, getcellset(grid, "FE"))
add!(sdh, :u, ip)
close!(dh2)

display(allocate_matrix(dh1))
display(allocate_matrix(dh2))

@KnutAM
Copy link
Member

KnutAM commented Jul 15, 2024

Thanks for the report and investigation!

My guess is that this is because we don't use OrderedCollections like Ferrite#master does, so this should probably be a pretty simple fix?

Adding my quick notes on how to make the change here (I don't have time to do it and test it right now), so for future me or others who wants to do this, the fix should be to change

facetset = sizehint!(Set{FacetIndex}(), length(nodeset))

and

cellsets = Dict(key => Set(getnumbers(rawelements)) for (key, rawelements) in rawelementsdict)
merge!(cellsets, Dict(key => Set(nums) for (key, nums) in rawelementsets))

to OrderedSet. The only challenge is to support backwards compatibility with Ferrite v0.3. This could be done by changing

@static if !isdefined(Ferrite, :SerendipityQuadraticHexahedron)
const SerendipityQuadraticHexahedron = Ferrite.Cell{3,20,6}
const SerendipityQuadraticQuadrilateral = Ferrite.Cell{2,8,4}
else
const SerendipityQuadraticHexahedron = Ferrite.SerendipityQuadraticHexahedron
const SerendipityQuadraticQuadrilateral = Ferrite.SerendipityQuadraticQuadrilateral
end
const FacetsDefined = isdefined(Ferrite, :FacetIndex) # Ferrite after v1.0 (Ferrite#914)
const FacetIndex = FacetsDefined ? Ferrite.FacetIndex : Ferrite.FaceIndex
const facets = FacetsDefined ? Ferrite.facets : Ferrite.faces
const addfacetset! = FacetsDefined ? Ferrite.addfacetset! : Ferrite.addfaceset!

to

const FerritePreV1 = !isdefined(Ferrite, :FacetIndex)

const SerendipityQuadraticHexahedron    = FerritePreV1 ? Ferrite.Cell{3,20,6} : Ferrite.SerendipityQuadraticHexahedron
const SerendipityQuadraticQuadrilateral = FerritePreV1 ? Ferrite.Cell{2,8,4}  : Ferrite.SerendipityQuadraticQuadrilateral

const FacetIndex   = FerritePreV1 ? Ferrite.FaceIndex   : Ferrite.FacetIndex
const facets       = FerritePreV1 ? Ferrite.faces       : Ferrite.facets
const addfacetset! = FerritePreV1 ? Ferrite.addfaceset! : Ferrite.addfacetset!

and use FerritePreV1 to convert to Set. For code-simplicity, I suggest we do the extra conversion on the return part even if that is slower. I suspect most people will uppgrade to Ferrite v1 anyways, and then it is just a matter of supporting older versions during the transition phase.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants