Fenics as Solver

About the development of the FEM module/workbench.

Moderator: bernd

Post Reply
joha2
Posts: 303
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Post by joha2 »

HoWil, first of all big thanks for your help! I still have difficulties preparing a simulation with a gravity loadcase in FreeCAD 0.16 (are there no volume forces implemented in this version?). For the (exact) reproducibility of the results, it might be useful to import the fenics mesh into FC or vice versa and use the same elements. Is this possible?

Anyway, I will have a look at your fenics code, although I am still a beginner. Another question coming to my mind in this context is: Which type of weak formulation do we have to use in fenics to provide u boundary conditions (i.e. fixed boundary)? The standard formulation int dV div sigma(u) v = -int dV sigma(u) epsilon(v) + int dA sigma n v seems not useful to provide direct u boundary conditions.

Best wishes and happy new year
Johannes
HoWil
Veteran
Posts: 1279
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Post by HoWil »

joha2 wrote:HoWil, first of all big thanks for your help! I still have difficulties preparing a simulation with a gravity loadcase in FreeCAD 0.16 (are there no volume forces implemented in this version?).
It is only available in 0.17. If you are on ubuntu simply install the nightly version in parallel to 0.16.
joha2 wrote: For the (exact) reproducibility of the results, it might be useful to import the fenics mesh into FC or vice versa and use the same elements. Is this possible?
You are completely right. I do not have a clue how to export/import the mesh. It seams like fenics can use gmsh meshes (see http://people.sc.fsu.edu/~jburkardt/exa ... /gmsh.html) maybe bernd can help here ; or .off files which are directly accessible from the export menu (https://fenicsproject.org/olddocs/dolfi ... ation.html).
joha2 wrote: Anyway, I will have a look at your fenics code, although I am still a beginner.
joha2 wrote: Me too....
quote]
Another question coming to my mind in this context is: Which type of weak formulation do we have to use in fenics to provide u boundary conditions (i.e. fixed boundary)? The standard formulation int dV div sigma(u) v = -int dV sigma(u) epsilon(v) + int dA sigma n v seems not useful to provide direct u boundary conditions.

Best wishes and happy new year
Johannes
I am not sure if I understand your question right. Please see my example where

Code: Select all

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

was also used here....
https://github.com/hplgit/fenics-tutori ... sticity.py
and in a similar manner here...
https://github.com/FEniCS/dolfin/blob/m ... sticity.py

Again, I am just starting learning to use FeniCS maybe others here are more experienced in this things.
User avatar
bernd
Veteran
Posts: 12849
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland
Contact:

Re: Fenics as Solver

Post by bernd »

HoWil wrote:
joha2 wrote: For the (exact) reproducibility of the results, it might be useful to import the fenics mesh into FC or vice versa and use the same elements. Is this possible?
You are completely right. I do not have a clue how to export/import the mesh. It seams like fenics can use gmsh meshes (see http://people.sc.fsu.edu/~jburkardt/exa ... /gmsh.html) maybe bernd can help here ; or .off files which are directly accessible from the export menu (https://fenicsproject.org/olddocs/dolfi ... ation.html).
What mesh format does fenics uses?
joha2
Posts: 303
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Post by joha2 »

Usually they use some kind of simple xml format. Although as far as I know there are some simple converters to other mesh formats available. The xml format looks like the following:

Code: Select all

<?xml version="1.0" encoding="UTF-8"?>

<dolfin xmlns:dolfin="http://www.fenics.org/dolfin/">
  <mesh celltype="tetrahedron" dim="3">
    <vertices size="2129">
      <vertex index="0" x="1" y="0" z="0" />
      <vertex index="1" x="0" y="1" z="0" />
      <vertex index="2" x="-1" y="0" z="0" />
      <vertex index="3" x="0" y="-1" z="0" />
      <vertex index="4" x="0" y="0" z="-1" />
      <vertex index="5" x="0" y="0" z="1" />
      <vertex index="6" x="0.994522" y="0.104528" z="0" />
      <vertex index="7" x="0.978148" y="0.207912" z="0" />
.
.
.
      <vertex index="2128" x="0.318039" y="0.64495" z="-0.325478" />
    </vertices>
    <cells size="9312">
      <tetrahedron index="0" v0="336" v1="429" v2="437" v3="1554"/>
      <tetrahedron index="1" v0="336" v1="386" v2="1530" v3="1613"/>
.
.
.
      <tetrahedron index="9311" v0="2125" v1="2126" v2="2127" v3="2128"/>
    </cells>
  </mesh>
</dolfin>

HoWil
Veteran
Posts: 1279
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Post by HoWil »

bernd wrote:What mesh format does fenics uses?
joha2 wrote:Usually they use some kind of simple xml format. Although as far as I know there are some simple converters to other mesh formats available. The xml format looks like the following:
Bernd can you test "dolfin-convert.py, a Python script for converting mesh files" on at least a 2D gmsh-mesh.
Please see http://people.sc.fsu.edu/~jburkardt/py_ ... nvert.html where you also find the .py files.

There is also written that there are problems with 3D-gmsh-meshes:
Note that the examples of converting a 3D GMSH file to DOLFIN XML file have failed. This is because the GMSH 3D tetrahedral elements are listed as type 29, whereas DOLFIN-CONVERT is expecting an element type of 4. Thanks to Sean Breckling for reporting this issue, 13 May 2016.
EDIT:
I also found this
http://scicomp.stackexchange.com/questi ... .0699#7527
but could not use 'meshconvert.convert2xml(ifilename, ofilename, iformat=iformat)' as explained ... do not know what I am doing wrong...

Code: Select all

In [1]: from dolfin_utils import meshconvert

In [2]: meshconvert.convert2xml('shape2mesh.msh', 'shape2mesh.xml', iformat='gmsh')
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-2-b7193789cfe8> in <module>()
----> 1 meshconvert.convert2xml('shape2mesh.msh', 'shape2mesh.xml', iformat='gmsh')

AttributeError: 'module' object has no attribute 'convert2xml'
Here is also a gmsh mesh for testing:
shape2mesh.msh.zip
3D cantilever mesh as gmsh input example
(7.59 KiB) Downloaded 70 times
User avatar
looo
Veteran
Posts: 3941
Joined: Mon Nov 11, 2013 5:29 pm

Re: Fenics as Solver

Post by looo »

Short question: whats the reason for having interface files instead of python data. It makes sense if the solver uses these files directly, but with fenics it's possible to create the mesh with python, so why not use python dicts and lists to create the interface?
HoWil
Veteran
Posts: 1279
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Post by HoWil »

The main point for using fenics as an 'external' tool is the license issue see viewtopic.php?f=18&t=4677&start=10#p146173.
But you are right one should test it while directly handing over the data to fenics, or at least saving the 'python'-mesh data (points) as python pickle file *) or. hdf for lager meshes (more complicate structure).
This is more universal in my opinion.
BR
Howil

*) it is really quite simple and standard in python ; see examples over here https://docs.python.org/3.5/library/pickle.html
User avatar
looo
Veteran
Posts: 3941
Joined: Mon Nov 11, 2013 5:29 pm

Re: Fenics as Solver

Post by looo »

pickle is also an option.

The implementation I am working on treats solvers as external "plugins" which can be shipped by the FreeCAD-plugin-installer, pip, conda... . I think there is no license issue this way. As long as the solver-plugin isn't installed there is no use of gpl software.
This way it's possible to create the solver in the main thread and run some methods in a separated QThread. So passing simple python structure is no problem. Creating the mesh is then part of the external solver.

The fenics-mesh could be setup this way, but I haven't yet found resources how to setup the boundary conditions.
https://github.com/looooo/FreeCAD/blob/ ... .py#L14L27

currently I am trying to run some functions in a separate thread (prepare, solve, writeMesh). But I am not really satisfied with this. Somehow it would be nice to have FreeCAD and solver separated, but then the interaction is quite limited.
https://github.com/looooo/FreeCAD/blob/ ... .py#L19L29
joha2
Posts: 303
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Post by joha2 »

Finally I was able to export the mesh from FreeCAD to fenics to reproduce the linear elasticity example. The results are OK, but not really identical to each other. Moreover they differ from the former results of HoWil by a factor of 2. But this can be compensated by increasing the number of elements. (first result is automatic meshing in FreeCAD, second is max element size 200mm)

The fenics output is given by:

Code: Select all

(1) min/max u: -3.78642684121e-07 0.00104275547029

(2) min/max u: -5.16023990138e-08 0.00199829391954
The value shown in FreeCAD FEM is:

Code: Select all

(1) umax = 1.04304 mm

(2) umax = 1.99835 mm
First of all, I will describe the procedure to export the mesh from FreeCAD into fenics:
  1. Create FEM Analysis in FC and use gmsh meshing with 3D first order elements (higher order elements are not recognized by dolfin-convert)
  2. Perform FEM Analysis and look at the results
  3. Use mesh.FemMesh.write("./boxmesh.unv") to export the mesh from FreeCAD into an .unv file
  4. Use gmsh to read the .unv file and save the mesh as .msh file
  5. Use dolfin-convert.py to transform the .msh file into .xml readable by fenics
  6. Replace the mesh generation command in the demo file by mesh = Mesh('./boxmesh.xml')
  7. If you used mm units in FreeCAD and meter units in fenics you have to rescale the mesh by using x = mesh.coordinates(); x[:,:] *= 0.001
  8. Perform fenics simulation
The transfer seems quite complicated, but it's just for testing reasons. The difference from the HoWil's demo file seems
to come from using first order elements which are maybe too stiff. Does anybody have analytical results for comparison? :-)

Best wishes
Johannes
Last edited by joha2 on Sat Feb 11, 2017 9:01 am, edited 1 time in total.
HoWil
Veteran
Posts: 1279
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Post by HoWil »

Hi Johannes,
This sounds promising. I remember I had trouble to run dolfin-convert.py.
I will into it within the next days.
Can you please provide some project files? Maybe a zip out the complete directory including dolfin-convert.py and all the scripts you used ?
A long as the mesh is not to fine everything should fit into the zip.
Br
Howil
Post Reply