Looking to contribute to Sfepy

[Alex Eftimiades]

I was working on a discrete exterior calculus framework about a year ago. When it came time to start benchmarking it, I was disappointed with my options for mesh generation. I wanted to test my framework on different geometries, topologies, and dimensions in a consistent way. Since I could not find a good way of doing this, I wrote my own in Cython. I apologize for the poor documentation. It turns out I ended up reinventing "Kuhn triangulation" as documented in section 2.2.1 and 2.2.2 in this thesis. I am hoping to get more involved collaborating with open source projects as opposed to working on my own projects. I noticed Sfepy distinctly noted lacking BSD licensed mesh generation code in a 2008 power point but seemed to have made some progress since then as documented here.

My implementation allows you to specify an arbitrary number of subdivisions along an arbitrary number of axes, and contains two algorithms for generating a corresponding simplicial mesh. Pictures of the result of each algorithm for a 2D discretization are attached. The process looks like this:

1. Create a dict that maps tuples of integers (number of "steps" along each axis) to a global index of the corresponding vertex.
>> shape = [nx, ny, ....]
>> g = grid_indices(shape)

2. Connect the vertexes according to one of two algorithms for generating a uniform simplicial mesh (the difference in the resulting meshes is readily seen in the attached pictures). This is stored as an NxD dimensional numpy array of unsigned integers for a mesh with N D-dimensional simplexes. Each integer in the array is the global index of a vertex in the mesh.
>> simplices = grid(g)

3. Customize topology of mesh by stitching pairs of vertexes together. This is accomplished by replacing each occurrence of an index in the above array with the index of the vertex you want to "stitch" it to. So in the case of a one dimensional chain of 10 vertexes, I could create a "ring"-like topology by replacing each occurrence of the integer "9"  (the rightmost vertex) with a "0" (the leftmost vertex) in my 9x1 dimensional numpy array of integers (as described in step 2). Since periodic boundary conditions are very common, I wrote a routine specifically for toroidal topologies that are periodic along arbitrary combinations of axes. The routine returns a dict that maps vertex indexes to other vertex indexes such that the key index is to be replaced with the value index. So in this example, it would return {9: 0} to indicate vertex 9 is to be replaced with vertex 0.
>> stitches = pbc_stitches(g, shape, [0,3]) #imposing periodic boundary conditions along the first and fourth axes.

4. I found it most useful to keep the original "unstitched" simplexes too so geometry can be specified separately by embedding the vertexes within a specific coordinate system. For example, I might wish to impose periodic boundary conditions along the second axis in a 3 dimensional system and leave the system in cartesian coordinates. I might alternatively choose to embed the points in a cylindrical coordinate system and have the angular coordinate periodic.

I don't know if you would find this useful. If not, is there anything else I can do to contribute to Sfepy? I read the how to contribute guide, but your development section did say to let you know on your mailing list that I was interested in contributing to Sfepy. I am generally looking to gain experience contributing to open source projects. Since I spend a lot of time writing Python programs that numerically solve PDEs (typically eigenvalue problems at the moment), I thought sfepy might be a good place to start.

Thanks,
Alex Eftimiades

[Robert Cimrman]

Hello Alex,

yes, contributions are welcome!

I am not sure if your code is directly applicable - am I right that it
uniformly meshes a rectangle or cube by a regular simplex mesh? But anyway, let
us try to find something that would be interesting for you and useful for
sfepy. Did you manage to install it without problems?

Cheers,
r.

On 03/27/2015 02:10 AM, Alex Eftimiades wrote:
>
>
> I was working on a discrete exterior calculus framework about a year ago.
> When it came time to start benchmarking it, I was disappointed with my
> options for mesh generation. I wanted to test my framework on different
> geometries, topologies, and dimensions in a consistent way. Since I could
> not find a good way of doing this, I wrote my own

> <https://github.com/aeftimia/pyfeec/blob/master/pyfeec/grid_utils.pyx> in

> Cython. I apologize for the poor documentation. It turns out I ended up
> reinventing "Kuhn triangulation" as documented in section 2.2.1 and 2.2.2
> in this thesis

> <https://dspace.library.cornell.edu/bitstream/1813/6140/1/92-1322.pdf>. I

> am hoping to get more involved collaborating with open source projects as
> opposed to working on my own projects. I noticed Sfepy distinctly noted
> lacking BSD licensed mesh generation code in a 2008 power point

> <http://www.ondrejcertik.com/media/euroscipy2008.pdf> but seemed to have

> made some progress since then as documented here

> <http://sfepy.org/doc-devel/src/sfepy/mesh/mesh_generators.html>.

> <http://sfepy.org/doc-devel/developer_guide.html#how-to-contribute>, but

[Alex Eftimiades]

On Friday, March 27, 2015 at 10:45:19 AM UTC-4, Robert Cimrman wrote:

Hello Alex,

yes, contributions are welcome!

I am not sure if your code is directly applicable - am I right that it
uniformly meshes a rectangle or cube by a regular simplex mesh?

The idea was to uniformly mesh manifolds of arbitrary topology and dimension by breaking it into partitions that can be mapped K-cells, then gluing arbitrary subsets of the cells together.
 

But anyway, let
us try to find something that would be interesting for you and useful for
sfepy. Did you manage to install it without problems?

I just ran the unit tests on the version on github, and everything seems to be working, though test_input_navier_stokes2d_iga.py complained about not being able to find "igakit".

I'm not exactly sure what direction(s) sfepy is looking to expand in, but I'd love to help in whatever way I can.
I've spent a fair amount of time discretizing wave equations (usually electromagnetic and acoustic) if that helps.

 

Cheers,
r.

Thanks,
Alex

[Robert Cimrman]

On 03/27/2015 05:31 PM, Alex Eftimiades wrote:
> On Friday, March 27, 2015 at 10:45:19 AM UTC-4, Robert Cimrman wrote:
>>
>> Hello Alex,
>>
>> yes, contributions are welcome!
>>
>> I am not sure if your code is directly applicable - am I right that it
>> uniformly meshes a rectangle or cube by a regular simplex mesh?
>
> The idea was to uniformly mesh manifolds of arbitrary topology and
> dimension by breaking it into partitions that can be mapped K-cells, then
> gluing arbitrary subsets of the cells together.

OK I see (sort of - I know nothing about discrete exterior calculus etc -
forgive me my stupid questions :)). Can it be used to triangulate other shapes
than K-cells? More examples would be nice, if you have some.

>> But anyway, let
>> us try to find something that would be interesting for you and useful for
>> sfepy. Did you manage to install it without problems?
>>
>
> I just ran the unit tests on the version on github, and everything seems to
> be working, though test_input_navier_stokes2d_iga.py complained about not
> being able to find "igakit".

Good. Yes, igakit is needed for that one. We should add it as a dependency into
setup.py etc.

> I'm not exactly sure what direction(s) sfepy is looking to expand in, but
> I'd love to help in whatever way I can.
> I've spent a fair amount of time discretizing wave equations (usually
> electromagnetic and acoustic) if that helps.

Yes, the more of your background I know, the better I could steer you to a
suitable topic. Electromagnetic problems are among the things that are missing
in sfepy, and would be nice to have. Do you have some experience with the
vector elements (e.g. of Nedelec or Raviart-Thomas type)? What is your
experience with FEM? Do you know other discretization methods, such as
discontinuous Galerkin?

As a start, I suggest you read the tutorial and the primer, so that you get a
feel of how the code can be used. Also check the examples.

Thanks for your interest!
r.

[Alex Eftimiades]

On 03/27/2015 04:48 PM, Robert Cimrman wrote:

On 03/27/2015 05:31 PM, Alex Eftimiades wrote:

On Friday, March 27, 2015 at 10:45:19 AM UTC-4, Robert Cimrman wrote:

Hello Alex,

yes, contributions are welcome!

I am not sure if your code is directly applicable - am I right that it
uniformly meshes a rectangle or cube by a regular simplex mesh?

The idea was to uniformly mesh manifolds of arbitrary topology and
dimension by breaking it into partitions that can be mapped K-cells, then
gluing arbitrary subsets of the cells together.

OK I see (sort of - I know nothing about discrete exterior calculus etc - forgive me my stupid questions :)). Can it be used to triangulate other shapes than K-cells? More examples would be nice, if you have some.

Discrete exterior calculus is not really mainstream--at least in its most general form. It seems to have been reinvented for different purposes a few times. On such example is the Yee scheme commonly used for electromagnetics simulations.

Anyway, about the the triangulation code: The idea was to be able to triangulate just about anything you'd come across in physics or engineering. I attached some pictures of eigenmodes of triangulated structures. The structures triangulated are a torus, a Möbius strip, and a Klein bottle. I've also triangulated circles, spheres and cylinders. I also attached the corresponding meshes (though with less elements so you could see the pattern). In principle, you could triangulate 8 dimensional regions embedded in 12 dimensional spaces if you had the inclination and computational power. The idea was to keep the framework as general as possible.

But anyway, let
us try to find something that would be interesting for you and useful for
sfepy. Did you manage to install it without problems?

I just ran the unit tests on the version on github, and everything seems to
be working, though test_input_navier_stokes2d_iga.py complained about not
being able to find "igakit".

Good. Yes, igakit is needed for that one. We should add it as a dependency into setup.py etc.

I'm not exactly sure what direction(s) sfepy is looking to expand in, but
I'd love to help in whatever way I can.
I've spent a fair amount of time discretizing wave equations (usually
electromagnetic and acoustic) if that helps.

Yes, the more of your background I know, the better I could steer you to a suitable topic. Electromagnetic problems are among the things that are missing in sfepy, and would be nice to have. Do you have some experience with the vector elements (e.g. of Nedelec or Raviart-Thomas type)? What is your experience with FEM? Do you know other discretization methods, such as discontinuous Galerkin?

My background has been rewriting an open source implementation of discrete exterior calculus, PyDEC so that it would work under more general conditions, in parallel, and in a generally faster manner. I admittedly lack industry experience with FEM, and I am hoping to gain some more practical experience from collaborating on open source projects.

As you can see from the attached pictures, I have some experience with eigenvalue problems and wave propagation. I am looking to expand my skill set working other problems.

As a start, I suggest you read the tutorial and the primer, so that you get a feel of how the code can be used. Also check the examples.

Will do.

[Robert Cimrman]

On 03/27/2015 11:24 PM, Alex Eftimiades wrote:
> On 03/27/2015 04:48 PM, Robert Cimrman wrote:
> Anyway, about the the triangulation code: The idea was to be able to
> triangulate just about anything you'd come across in physics or engineering. I
> attached some pictures of eigenmodes of triangulated structures. The structures
> triangulated are a torus, a Möbius strip, and a Klein bottle. I've also
> triangulated circles, spheres and cylinders. I also attached the corresponding
> meshes (though with less elements so you could see the pattern). In principle,
> you could triangulate 8 dimensional regions embedded in 12 dimensional spaces
> if you had the inclination and computational power. The idea was to keep the
> framework as general as possible.

Nice examples!

Now to have that capability in sfepy in some well-scriptable way would be
handy. The code supporst only 2D and 3D elements, so the generality you mention
is not really needed, but I guess it does not hurt. The examples show 2D
surfaces embedded in 3D, am I right? Can your code do also 3D volumes in 3D?

>>
>>> I'm not exactly sure what direction(s) sfepy is looking to expand in, but
>>> I'd love to help in whatever way I can.
>>> I've spent a fair amount of time discretizing wave equations (usually
>>> electromagnetic and acoustic) if that helps.
>>
>> Yes, the more of your background I know, the better I could steer you to a
>> suitable topic. Electromagnetic problems are among the things that are
>> missing in sfepy, and would be nice to have. Do you have some experience with
>> the vector elements (e.g. of Nedelec or Raviart-Thomas type)? What is your
>> experience with FEM? Do you know other discretization methods, such as
>> discontinuous Galerkin?
>
> My background has been rewriting an open source implementation of discrete
> exterior calculus, PyDEC

> <http://www.math.uiuc.edu/%7Ehirani/papers/BeHi2012_TOMS.pdf> so that it would
> work under more general conditions <http://arxiv.org/abs/1405.7879>, in

> parallel, and in a generally faster manner. I admittedly lack industry
> experience with FEM, and I am hoping to gain some more practical experience
> from collaborating on open source projects.

The paper mentions other codes (e.g. dolfin) - interestingly, I have
implemented our CMesh class according to the article by Anders Logg (Efficient
Representation of Computational Meshes. 2012) which uses the topological ideas
that might be related to what you do. So going in this direction - using
general mathematical ideas to unify many concepts related to
meshes/discretization - seems very useful. Now in petsc they have DMPlex grid...

> As you can see from the attached pictures, I have some experience with
> eigenvalue problems and wave propagation. I am looking to expand my skill set
> working other problems.

That is useful. What about eigenvalue problem solvers - just curious which ones
you have used.

r.

[Alex Eftimiades]

Yes, those were 2D objects embedded in 3D, and yes it can do 3D objects
in 3D space (in general, N-D objects in M-D space). I spent a few hours
today trying to get 3D visualization working on my computer again (I had
it working before). Fortunately, I got it. I attached a triangulated
sphere and (solid) torus.

I skimmed through the article you mentioned, "Efficient Representation
of Computational Meshes" and it looked very interesting--and indeed very
relevant to the types of approaches I take. I have only begun to look at
the paper since I spent so much time nearly crashing my computer to get
my graphics drivers working again for the 3D visualization.

As for eigenvalue problem solvers--I've just discretized the laplacian
and used standard linear algebra packages to diagonalize the resulting
matrices. I've used scipy and pycula for that.

At least in finite element exterior calculus, the discrete laplacian
looks like d.T.dot(mass_matrix).dot(d), where "d" is an incidence matrix
between edges and vertices for scalar fields, triangles and edges for
vector fields, etc (for higher dimensional arbitrary rank antisymmetric
tensor fields), and "mass_matrix" is some symmetric matrix that is
vaguely analogous to an inner product. You then solve a generalized
eigenvalue problem of the form

>> eigh(laplacian, other_mass_matrix)

I think this is similar to how finite elements works in practice--but
I'm more familiar the problem from this angle. I noticed the "Efficient
Representation of Computational Meshes" seems to store meshes in a way
in which this perspective would be very natural--which I why I thought
it looked very relevant to my usual approach.

So far, I've solved for normal modes of scalar and vector fields in
simple domains (such as rectangular and circular) with Dirichlet,
Neumann, and periodic boundary conditions just to benchmark my software
against known results. But the process looks exactly the same for any
domain once you set up the mesh.

Hopefully some of this might be useful for sfepy.

Thanks,
Alex

[Robert Cimrman]

On 03/28/2015 10:04 PM, Alex Eftimiades wrote:
> Yes, those were 2D objects embedded in 3D, and yes it can do 3D objects in 3D
> space (in general, N-D objects in M-D space). I spent a few hours today trying
> to get 3D visualization working on my computer again (I had it working before).
> Fortunately, I got it. I attached a triangulated sphere and (solid) torus.

Nice! So maybe you could start by documenting your code, so that we can try it
as well easily? That would be useful for you even in case it is not all needed
for our purposes.

> I skimmed through the article you mentioned, "Efficient Representation of
> Computational Meshes" and it looked very interesting--and indeed very relevant
> to the types of approaches I take. I have only begun to look at the paper since
> I spent so much time nearly crashing my computer to get my graphics drivers
> working again for the 3D visualization.
>
> As for eigenvalue problem solvers--I've just discretized the laplacian and used
> standard linear algebra packages to diagonalize the resulting matrices. I've
> used scipy and pycula for that.
>
> At least in finite element exterior calculus, the discrete laplacian looks like
> d.T.dot(mass_matrix).dot(d), where "d" is an incidence matrix between edges and
> vertices for scalar fields, triangles and edges for vector fields, etc (for
> higher dimensional arbitrary rank antisymmetric tensor fields), and
> "mass_matrix" is some symmetric matrix that is vaguely analogous to an inner
> product. You then solve a generalized eigenvalue problem of the form
>
> >> eigh(laplacian, other_mass_matrix)
>
> I think this is similar to how finite elements works in practice--but I'm more
> familiar the problem from this angle. I noticed the "Efficient Representation
> of Computational Meshes" seems to store meshes in a way in which this
> perspective would be very natural--which I why I thought it looked very
> relevant to my usual approach.

This the the way the meshes are stored now in sfepy as well, in C a cython. I
am planning just now to use the CMesh more and diminish the use of Mesh.

> So far, I've solved for normal modes of scalar and vector fields in simple
> domains (such as rectangular and circular) with Dirichlet, Neumann, and
> periodic boundary conditions just to benchmark my software against known
> results. But the process looks exactly the same for any domain once you set up
> the mesh.
>
> Hopefully some of this might be useful for sfepy.

I think it will. We will see more after you get more acquainted with our code
(pointing out poorly documented parts is a useful contribution too) and vice
versa.

r.