Let's discuss how to best integrate it in sfepy. The main code is in:
geom/meshutils.py
which converts to/from a lot of formats. It's quite long though and all in one class. My plan was to create a new class for each output/input format.
To have transparent read support of various formats, I have created the meshio.py file (sfe.fem.meshio) - there are curently two classes, for medit (.mesh) and for legacy VTK files. I suggest you split your class and put a class for each format in there. Or I may do it, but not today, I fear. I will atleast document the expected return value of the read() method of those classes.
This is all I need. Maybe even a new python file for each format? I.e. creating a new module "sfe.fem.io" and implement all the classes in there as separate files? If you are ok with that, I'll do it.
However, it's not that simple. No subclass in meshio.py implements the "write" method, this is all implemented in:
sfe/base/ioutils.py
There should imho be just one code (all in meshio.py?). What I did (see my patch), is, that the write methods accept an optinal keyword argument specifying the solution (scalar, vectors, or tensors, or both/all of it) and if the argument is present, the "write" method puts the solution together with the mesh into the file (if the format supports it).
My idea is to have one .py file for each format (for example in sfe/fem/io dir?) that implements everything there is to this format - reading mesh, writing mesh, writing a solution, so that it's not scattered all over sfepy. Do you agree? Are there any pitfalls regarding the sfe/base/ioutils.py file?
So let's create some basic mesh classes, for example based on the code in sfe/fem/mesh.py. Then I'll just subclass it, implement some output/input and that's it.
See above. What is needed is only a class with read() and write() (this one is not used anywhere yet) methods. The Mesh class in sfe/fem/mesh.py automatically selects the 'io' class then.
It's a little more complex, see above.
Another problem: how can I select some part of the geometry and assign a special boundary condition to it? In sfepy, as I understand it, it's handled by assigning a special number to each element (tetrahedron), but how can I generate such a mesh? It's not a business of sfepy. But generally I need to create a geometry (for example in gmsh), tell gmsh to mark some surface with a number. Then it is passed to tetgen, that is clever enough to assign this number to each element in the output. So all I need is to parse the output of tetgen and print the element numbers to a file for example.
Sometimes, I not only need the element IDs, but I also need the set of nodes sitting at a surface. All these features are implemented in the patch attached to the issue 17.
In my older FEM code, I simply read this file and assigned a correct BC to the elements in this file. How is this handled by sfepy?
Short answer: regions :-)
True, you can select elements by their id (the special number you mention), but there are many other ways
- select nodes by coordinates: region_03 = { 'name' : 'Gamma_Left', 'select' : 'nodes in (x < 0.00001)', }
- nodes of surface + subtracting/adding other regions: region_0 = { 'name' : 'Walls', 'select' : 'nodes of surface -n (r.Outlet +n r.Inlet)', 'canCells' : False, }
Excellent, regions is all I need. I am able to specify regions in gmesh, and it is preserved by tetgen in the mesh. So then it should easily work in sfepy even now. Awesome.
Ondrej