Some comments about my mesh situation/implementation.

[KristofferC]

I wrote on the mailing list and was asked to post here for some discussion regarding a unified mesh representation.

I will write my use case and my implementation and maybe it can create some discussion / comments.

Currently, I am writing some code for Finite element analysis in Julia and obviously the concept of a mesh comes up. In my code I have two types of meshes, GeoMesh and FEMesh.

GeoMesh is a mesh in the geometrical sense. It contains nodes, elements and currently two types of sets, element sets and node sets. Later, I later also want to include edge and surface sets. At this moment I am not sure that having the sets a part of the actual mesh type is the best way.

My current implementation can be found here:
https://gist.github.com/KristofferC/a1a884a0af1818653c3d

By the way, I noticed today that I seem to have duplicted some stuff there that exists in GeometricalPredicates.jl

One problem with this implementation is that GeoMesh.elements is type unstable. The reason for this is that a priori I do not know what type of elements will be used. However, for my purpose that is ok since I just don't do any computations with GeoMesh, I just use it to hold some information. I realize for other use cases this is unacceptable.

FEMesh has a similar structure to GeoMesh except the type of nodes and elements are different. Here the elements are actual finite elements. These includes much more information than the simple geometrical elements. For example, a finite element include what type of interpolator is used to interpolate fields inside the element and what type of degrees of freedom exist in the nodes etc.

The FEMesh is created from the GeoMesh by a function that also takes a mapping of what finite element each geometric element should become.

To visualize the mesh and the results I currently use Python's VTK bindings and call them with PyCall similarly to this (except I also give the values of the strain/stress tensors at the nodes): https://gist.github.com/KristofferC/9a989a2c0518bb70009c

It would be nice if I could just use some of the stuff under the JuliaGeometry organization instead of rolling my own so if you have any idea how to do that, feel free to share.

Best regards, Kristoffer Carlsson.

[sjkelly]

Thanks for sharing! One quick note: it looks like your type stability issues could be solved by parameterizing GeoMesh2 like so:

immutable GeoMesh2{T<:AbstractGeoElem2} <: AbstractGeoMesh
nodes::Vector{GeoNode2}
elements::Vector{T}
elementsets::Dict{ASCIIString, ElementSet}
nodesets::Dict{ASCIIString, NodeSet}
end 

We currently have some work in Meshes.jl to this end. I have started making Face primitives more generic. @SimonDanisch has a good proposal for using a Dict type for this (I see you did something similar), but I think the worry was access time compared to indexing a vector or regular field access. The amount of time I've been able to dedicate to Meshes has not been consistent lately, but I am hoping to work on closing JuliaGeometry/Meshes.jl#21 based on ideas in JuliaGeometry/Meshes.jl#36 in the near future.

Do you have any more literature I could look at for FEM meshing? It would be really nice to be able to leverage some of the Implicit CSG in Meshes.jl to be able to go to regular mesh export, visuals, and FEA.

[KristofferC]

I don't think your solution will work because I don't want to make the restriction that the GeoMesh2.elements have to contain elements of only one type.

I had some discussion on the mailing list about it here: https://groups.google.com/forum/?fromgroups=#!topic/julia-users/x-1gM_zJsRc

In the end I ended up using the solution in one of the last responses which do not remove the instability but makes it insignificant if the number of distinct types of elements are small (which it normally is).

Regarding meshing, I am not sure how much of help I can be. I have mostly used existing meshers and imported the resulting mesh.

[KristofferC]

I think the only way to completely remove the type instability is to create some sort of Type Factory that generates a Mesh type at run time with one Vector of each type of element or someting along those lines.

[KristofferC]

I see now that: JuliaGeometry/Meshes.jl#21 (comment) might actually completely solve the type instability problem for me.

[PetrKryslUCSD]

Interesting discussion.

I approached this slightly differently:

The mesh consists of sets. This allows for the mesh to be combined of
various types. It also allows for representation of element types needed
for evaluation of the interior and boundary integrals.

All elements in JFinEALE descend from an abstract type. The methods working
on the elements are parameterized by the concrete type.

Petr

I see now that: JuliaGeometry/Meshes.jl#21 (comment)
JuliaGeometry/Meshes.jl#21 (comment)
might actually completely solve the type instability problem.

—
Reply to this email directly or view it on GitHub
#3 (comment).

[KristofferC]

So you have a mesh type that contains of a vector of sets and each set is a vector of different element type? It is sort of like what I have now except I have a dictionary with element types mapped to a vector of elements of that type.

When I think about it I run into a second type instability now. I want to dispatch on two things, the type of element to calculate shape functions etc and also the material to calculate the constitutive stiffness matrix.

My (unoptimized) code for the element stiffness now look like:

function stiffness(elem::AbstractFElement, nodes::Vector{FENode}, material::AbstractMaterial)
    Ke = zeros(elem.n_dofs, elem.n_dofs)

    for gp in elem.gps
        Be = Bmatrix(elem, gp, nodes)
        De = material_stiffness(material, gp)
        dV = weight(elem, gp, nodes)
        Ke += Be' * De * Be * dV
    end

    return Ke
end

A priori, there is no knowledge of what element, material combination will be used so I need a good structure where the compiler doesn't have to check what concrete material and element happened to occur in every single iteration.

One way to do it could be maybe be a Dict{ (DataType, DataType), Vector} where I map the material type and element type to a vector of a concrete element type.

I then loop over all of the vectors in the dict and I then know both the material and the concrete element type. Since the combination of material types and element types are small it should be a low performance hit.

[KristofferC]

One thing to think about is if the performance hit is that bad actually from the type instability. Right now I am running a linear elastic material with linear triangles which mean that neither the constitutive stiffness nor the shape functions take any time at all to calculate.

For more complicated materials that has plasticity for example the time to lookup the concrete type might be insignificant.

[PetrKryslUCSD]

Actually one needs to parameterize also by the model reduction (1, 2, 3
dim, when <3 is it plane strain, stress, or axisymm). For instance

function stiffness{MR<:DeformationModelReduction,
S<:FESet,A<:SysmatAssemblerBase,T<:Number}(::Type{MR},
self::FEMMDeformationLinear{S},
assembler::A,
geom::NodalField{JFFlt},
u::NodalField{T})

On Tue, Mar 17, 2015 at 4:27 AM, Kristoffer Carlsson <
notifications@github.com> wrote:

So you have a mesh type that contains of a vector of sets and each set is
a vector of different element type? It is sort of like what I have now
except I have a dictionary with element types mapped to a vector of
elements of that type.

When I think about it I run into a second type instability now. I want to
dispatch on two things, the type of element to calculate shape functions
etc and also the material to calculate the constitutive stiffness matrix.

My (unoptimized) code for the element stiffness now look like:

function stiffness(elem::AbstractFElement, nodes::Vector{FENode}, material::AbstractMaterial)
Ke = zeros(elem.n_dofs, elem.n_dofs)

for gp in elem.gps
    Be = Bmatrix(elem, gp, nodes)
    De = material_stiffness(material, gp)
    dV = weight(elem, gp, nodes)
    Ke += Be' * De * Be * dV
end

return Keend

A priori, there is no knowledge of what element, material combination will
be used so I need a good structure where the compiler doesn't have to check
what concrete material and element happened to occur in every single
iteration.

One way to do it could be maybe be a Dict{ (DataType, DataType), Vector}
where I map the material type and element type to a vector of a concrete
element type.

I then loop over all of the vectors in the dict and I then know both the
material and the concrete element type. Since the combination of material
types and element types are small it should be a low performance hit.

—
Reply to this email directly or view it on GitHub
#3 (comment).

Petr Krysl
Professor
University of California, San Diego
Skype: Petr.Krysl.UCSD.EDU, Office: 858-822-4787
http://hogwarts.ucsd.edu/~pkrysl/

[KristofferC]

Yes, it is true, I just haven't added that yet. However, I think that the type of model assumption (plane stress, plane strain) is easier to deal with than with the element and materials because the model assumption is constant in the model so its concrete type will be known to the compiler.

Do you have any examples in JFinEALE where you are using multiple element types and materials in an analysis so maybe I can get some inspiration from your method.

Regarding the overhead of the instability, an analysis of a mesh with 80 000 elements seems to take ~1.3 seconds extra to assemble the stiffness matrix compared to the type stable version. For linear elasticity this more than doubles the analysis time.

[KristofferC]

I realize that some of my problems might come because I previously have mostly been programming OOP. In a OOP language I would typically have an Element class that has a material, gauss points, etc. However, since you don't know what material will be used, in Julia you have to put an AbstractMaterial in the element and then you have already lost the type stability.

In Julia, it seems better to take an approach where the types are simple and uncoupled and functions take more arguments instead of a big type that contains all information.

[PetrKryslUCSD]

On Tue, Mar 17, 2015 at 7:05 AM, Kristoffer Carlsson <
notifications@github.com> wrote:

Yes, it is true, I just haven't added that yet. However, I think that the
type of model assumption (plane stress, plane strain) is easier to deal
with than with the element and materials because the model assumption is
constant in the model so its concrete type will be known to the compiler.

Good point.

Do you have any examples in JFinEALE where you are using multiple element
types and materials in an analysis so maybe I can get some inspiration from
your method.

Yes, there is an example (3-d static stress analysis). Check out
JFinEALEexamples.

Regarding the overhead of the instability, an analysis of a mesh with 80
000 elements seems to take ~1.3 seconds extra to assemble the stiffness
matrix compared to the type stable version. For linear elasticity this more
than doubles the analysis time.

Yes, the penalty can be significant.

Petr Krysl
Professor
University of California, San Diego
Skype: Petr.Krysl.UCSD.EDU, Office: 858-822-4787
http://hogwarts.ucsd.edu/~pkrysl/

[KristofferC]

Thank you for the example.

Your "regions" contain one element set and one material. All the elements in that set are of the same type.

Let us say you have a mesh with different type of elements and you want to create one big region for the mesh that has the same material. Would you then have to split up the region into smaller regions that contain only one element type. For example, a mesh with quads and triangles and the same linear elastic material would have two regions in your code?

[PetrKryslUCSD]

Yes, one region - one element set - one element type.

I don't see any reason why a heterogeneous mesh entity based on elements
should exist. Mesh entity based on sets, yes.

For a mesh entity based on types, the operations need to dispach
dynamically based on type. In JFinEALE the operations work with a stable
type, because each operation runs on a set of a homogeneous type.

P

On Sun, Mar 22, 2015 at 7:31 AM, Kristoffer Carlsson <
notifications@github.com> wrote:

Thank you for the example.

Your "regions" contain one element set and one material. All the elements
in that set are of the same type.

Let us say you have a mesh with different type of elements and you want to
create one big region for the mesh that has the same material. Would you
then have to split up the region into smaller regions that contain only one
element type. For example, a mesh with quads and triangles and the same
linear elastic material would have two regions in your code?

—
Reply to this email directly or view it on GitHub
#3 (comment).

Petr Krysl
Professor
University of California, San Diego
Skype: Petr.Krysl.UCSD.EDU, Office: 858-822-4787
http://hogwarts.ucsd.edu/~pkrysl/

[KristofferC]

I mostly thought about it from a user friendliness perspective. When you apply a material to a region you don't care about what type of finite elements happen to be there. You just want to say that this geometric region has these material properties.

Of course, it is possible to just generate the type stable regions by looking at the set intersections between sets defining material parameters and sets defining element types.

[SimonDanisch]

did we arrive at a conclusion, or did you accidentally close this?

[KristofferC]

I could keep it open, it just felt like the discussion was a bit too specific to be useful for most people. But I guess the best way to represent a mesh where the elements have different properties can be useful not only in FEM.

[SimonDanisch]

Oh no, I think it's a great discussion.
It gives me some ideas about how the mesh type must look in order to be usable for FEM, which will be a large use case for it!

[PetrKryslUCSD]

Actually, the FE operations are only done on the entire mesh rarely. The
interior, sometimes, the boundary, very rarely.

In the interior integrals, having a single operation applied to the entire
mesh is possible, but in my opinion if there are multiple element types,
one would likely want to apply different VERSIONS of an operation: as an
example, consider elasticity, where for good stress representation on a
mixed eight-node serendipity quadrilateral and 10 node triangle mesh one
could resort to u-p interpolation, which would be done in different ways
for the quadrilaterals and triangles. Thinks like that are more conducive
to operations written for particular element types and executed on a batch
of such elements rather than these operations being dispatched dynamically
as that would have to be presumably based not only on the element type but
also on the input parameters (constant versus linear versus bilinear
pressure, for instance).

In the boundary integrals, one would always want to operate on subsets of
the boundary – different boundary conditions on different subsets.

Petr

On Sun, Mar 22, 2015 at 8:50 AM, Kristoffer Carlsson <
notifications@github.com> wrote:

I mostly thought about it from a user friendliness perspective. When you
apply a material to a region you don't care about what type of finite
elements happen to be there. You just want to say that this geometric
region has these material properties.

Of course, it is possible to just generate the type stable regions by
looking at the set intersections between sets defining material parameters
and sets defining element types.

—
Reply to this email directly or view it on GitHub
#3 (comment).

Petr Krysl
Professor
University of California, San Diego
Skype: Petr.Krysl.UCSD.EDU, Office: 858-822-4787
http://hogwarts.ucsd.edu/~pkrysl/

[KristofferC]

I think the data structures for a mesh will have to be different depending on the application. In FEM where you do a significant amount of work on each element it is not very important with dense data structures. Right now, I am thinking of using Dict{Int, Element} for example instead of Vector{Element}. When you have an adaptive mesh and you are creating and removing elements it is easier with the indexing if everything is in a dictionary.

However, if you want to do very little work over every element and instead have a large mesh the dictionary lookups might become significant and a dense structure like the one in Meshes.jl is probably to prefer.

[SimonDanisch]

I have updated the readme in MeshIO.jl to explain my new mesh type a little. Hope you like what I have come up with, and we can continue working on other mesh types and the conversions between them ;)
Please open issues about anything you find problematic or which doesn't make sense!

Best,
Simon