Fenics as Solver

About the development of the FEM module/workbench.

Moderator: bernd

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

Re: Fenics as Solver

Post by HoWil »

joha2 wrote: Sun May 27, 2018 10:03 am
HoWil wrote: Sat May 26, 2018 2:57 pm Can you please test if inner boundaries are already accessible jet? :?:
Hey HoWil,

as far as I know, inner boundaries should not be a problem. Since for my "heizung" example there is a room with an inner "heizung" volume which could be used as volume heat source via the cell function interface and as inner boundary condition via the facet function interface. I will have a look at your example during the next few days.

Best wishes
Johannes
Hi Johannes,
Thanks in advance for your help!
Is it maybe possible to share your "heizung" example, so I can compare it with my work?
BR,
HoWil
joha2
Posts: 303
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Post by joha2 »

HoWil wrote: Sun May 27, 2018 2:28 pm Is it maybe possible to share your "heizung" example, so I can compare it with my work?
Hey HoWil,

sorry for the long delay. The code for the "heizung" example is as follows (I did not remove the comments for using the cell function):

Code: Select all

# import Fenics main module

import fenics as fe

# import FreeCAD modules

import FreeCADGui
import FreeCAD

# import Fenics interface classes

from feminout.importFenicsMesh import writeFenicsXDMF
from femsolver.fenics.fenics_tools import CellExpressionXDMF, FacetFunctionXDMF

"""
set Python executable to an appropriate value to maintain instant functionality
during the make process
"""

import sys
sys.executable = '/usr/bin/python'

#fe.set_debug_level(fe.DBG)

"""
Select mesh in active document and export it to a temporary file
"""

sel = FreeCADGui.Selection.getSelection()
if len(sel) == 1 and FreeCAD.ActiveDocument is not None:
	sel = sel[0]	
	writeFenicsXDMF.write_fenics_mesh_xdmf(sel, "mesh.xdmf")

"""
Choose charge expressions via Fenics' expressions classes
"""

eaussen = fe.Expression("-10.0", degree=1)
einnen = fe.Expression("-5.0", degree=1)
eheizung = fe.Expression("10.0", degree=1)
eheizung_volume = fe.Expression("0.1", degree=1)

"""
Use Fenics interface class to import temporary mesh.
Use labels from mesh groups and use expressions for Dirichlet boundary conditions.
"""

ff = FacetFunctionXDMF("mesh.xdmf", 
	{"innen":{"type":"Dirichlet", "value":einnen},
	"aussen":{"type":"Dirichlet", "value":eaussen}, #}), 
    "heizung":{"type":"Dirichlet", "value":eheizung}})

#cf = CellExpressionXDMF("mesh.xdmf", {"heizung_volume":eheizung_volume}, degree=1)

lambda_air = fe.Constant(26.2)

"""
In the following use simple Fenics example for Poisson equation to calculate
an electrostatic type potential on the mesh and with the mentioned boundary conditions.
"""

V = fe.FunctionSpace(ff.mesh, "Lagrange", 1)
#W = fe.FunctionSpace(cf.mesh, "Lagrange", 1)

u = fe.TrialFunction(V)
v1 = fe.TestFunction(V)
#v2 = fe.TestFunction(W)

#dx_heizung = cf.dx["heizung_volume"](1)

f = fe.Constant(0.0)
my_form_lhs = lambda_air*fe.inner(fe.grad(u), fe.grad(v1))*fe.dx 
#my_form_rhs = cf*v2*dx_heizung
my_form_rhs = fe.Constant(0.0)*v1*fe.dx

#print(fe.assemble(1*cf.dx["heizung_volume"]))
#print(fe.assemble(1*cf.dx["heizung_volume"](0)))
#print(fe.assemble(1*cf.dx["heizung_volume"](1)))
print(fe.dx)
#print(cf.dx["heizung_volume"](1))

u = fe.Function(V)
fe.solve(my_form_lhs == my_form_rhs, u, ff.getDirichletBCs(V))

"""
PVD file can be inspected by using Paraview
"""


fe.File("sol_fc_fenics_raum_mit_heizung.pvd") << u


Further note that in the cut the inner boundaries have a temperature of 10 degrees like expected for the "heizung".


cut_with_temp_distribution.png
cut_with_temp_distribution.png (47.86 KiB) Viewed 2470 times

I hope that the macro code and the example geometry helps.

Best wishes
Johannes
Attachments
raum_mit_heizung.zip
(19.99 KiB) Downloaded 76 times
HoWil
Veteran
Posts: 1279
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Post by HoWil »

joha2 wrote: Sat Jun 02, 2018 8:14 am The code for the "heizung" example is as follows
Hi Johannes,

Nice model. I could run it and examined the results in paraview.
Nevertheless, I can reproduce my problem with accessing the inner boundaries :!: :shock:
Please try the following code with your example and look for the result of the integral over the face 'heizung' :

Code: Select all

print(fe.assemble(1*ff.ds['aussen'](1)))  # gives correct 2610000.0 mm²; area of 'outer' window and door
print(fe.assemble(1*ff.ds['innen'](1)))  # gives correct 1710000.0 mm²; area of 'inner' door
print(fe.assemble(1*ff.ds['heizung'](1)))  # gives 0.0 !! ; should result 2680000,00 mm²
BR,
HoWil

BTW: usually is mixing of different references a problem when creating a mesh. You have used the boolean-fragment 'raum_heizkoerper' as main Part for the gmsh-mesh 'raummesh' but than you mixed faces of 'raum_heizkoerper' and some of its child 'heizkoerper' (for 'heizung').

Tested with:
OS: Ubuntu 18.04 LTS
Word size of OS: 64-bit
Word size of FreeCAD: 64-bit
Version: 0.18.13780 (Git)
Build type: None
Branch: master
Hash: 9fb122008b99ea2d30ed3e6f7cc93a3b7717cce5
Python version: 2.7.15rc1
Qt version: 4.8.7
Coin version: 4.0.0a
OCC version: 7.2.0
Locale: English/UnitedStates (en_US)

EDIT and on
OS: Ubuntu 18.04 LTS
Version: 0.18.13886 (Git)
User avatar
looo
Veteran
Posts: 3941
Joined: Mon Nov 11, 2013 5:29 pm

Re: Fenics as Solver

Post by looo »

nice progress!
I tried running it with py3 and there was a problem somewhere in the writeFenicsXDMF function:

Code: Select all

Traceback (most recent call last):
  File "/home/lo/.FreeCAD/Macro/fenics1.py", line 33, in <module>
    writeFenicsXDMF.write_fenics_mesh_xdmf(sel, "mesh.xdmf")
  File "/home/lo/conda/conda-bld/freecad-debug_1522865186960/work/Mod/Fem/feminout/writeFenicsXDMF.py", line 266, in write_fenics_mesh_xdmf
    mesh_function_codim = dim_cell - FreeCAD_Group_Dimensions[mesh_function_type]
<class 'KeyError'>: (b'Face',)
version info:
OS: Ubuntu 16.04.4 LTS
Word size of OS: 64-bit
Word size of FreeCAD: 64-bit
Version: 0.18.13827 (Git)
Build type: Release
Branch: osx-conda-config
Hash: a285204fd9a509499410a747943765b8b545dfb9
Python version: 3.6.5
Qt version: 5.6.2
Coin version: 4.0.0a
OCC version: 7.2.0
Locale: German/Germany (de_DE)
joha2
Posts: 303
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Post by joha2 »

HoWil wrote: Sun Jun 03, 2018 6:57 pm Please try the following code with your example and look for the result of the integral over the face 'heizung' :

Code: Select all

print(fe.assemble(1*ff.ds['aussen'](1)))  # gives correct 2610000.0 mm²; area of 'outer' window and door
print(fe.assemble(1*ff.ds['innen'](1)))  # gives correct 1710000.0 mm²; area of 'inner' door
print(fe.assemble(1*ff.ds['heizung'](1)))  # gives 0.0 !! ; should result 2680000,00 mm²
I will test it. My conjecture is that due to inner and outer surfaces at inner boundaries and when you multiply each area element with some inner and the corresponding outer normal vector, their sum is identically zero. Maybe to correct this behaviour we have to make those boundaries real outer boundaries, i.e. we have to remove the inner part of "heizung" and make it hollow. :idea: Would that be an option?
HoWil wrote: Sun Jun 03, 2018 6:57 pm BTW: usually is mixing of different references a problem when creating a mesh. You have used the boolean-fragment 'raum_heizkoerper' as main Part for the gmsh-mesh 'raummesh' but than you mixed faces of 'raum_heizkoerper' and some of its child 'heizkoerper' (for 'heizung').
Thanks for the hint. The usage of boolean-fragments and references is still not very clear to me, so I did it in a way which works. So what would the best practise way to do it? :mrgreen:

Best wishes
Johannes
joha2
Posts: 303
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Post by joha2 »

looo wrote: Sun Jun 03, 2018 9:03 pm nice progress!
I tried running it with py3 and there was a problem somewhere in the writeFenicsXDMF function:

Code: Select all

Traceback (most recent call last):
  File "/home/lo/.FreeCAD/Macro/fenics1.py", line 33, in <module>
    writeFenicsXDMF.write_fenics_mesh_xdmf(sel, "mesh.xdmf")
  File "/home/lo/conda/conda-bld/freecad-debug_1522865186960/work/Mod/Fem/feminout/writeFenicsXDMF.py", line 266, in write_fenics_mesh_xdmf
    mesh_function_codim = dim_cell - FreeCAD_Group_Dimensions[mesh_function_type]
<class 'KeyError'>: (b'Face',)
So now is the question: Does this come from my definition of the mesh_function_type or from the return value of the meshgroup within FreeCAD? :mrgreen: :D I will check that. Thanks for the hint @looo.

Best wishes
Johannes
User avatar
looo
Veteran
Posts: 3941
Joined: Mon Nov 11, 2013 5:29 pm

Re: Fenics as Solver

Post by looo »

you are right. So fem_mesh.getGroupElementType(g) is expected to return default string? If this is true we have to change from PyBytes_FromString to PyUnicode_FromString.
User avatar
looo
Veteran
Posts: 3941
Joined: Mon Nov 11, 2013 5:29 pm

Re: Fenics as Solver

Post by looo »

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

Re: Fenics as Solver

Post by HoWil »

joha2 wrote: Mon Jun 04, 2018 6:25 am I will test it. My conjecture is that due to inner and outer surfaces at inner boundaries and when you multiply each area element with some inner and the corresponding outer normal vector, their sum is identically zero. Maybe to correct this behaviour we have to make those boundaries real outer boundaries, i.e. we have to remove the inner part of "heizung" and make it hollow. :idea: Would that be an option?
Yes, this is an option.
EDIT: I am not entirely sure if this is necessary. For boundary definition directly in fenics one does not need any knowledge of normal vectors and so on. I will try to prepare an example.

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

Re: Fenics as Solver

Post by HoWil »

joha2 wrote: Mon Jun 04, 2018 6:25 am I will test it. My conjecture is that due to inner and outer surfaces at inner boundaries and when you multiply each area element with some inner and the corresponding outer normal vector, their sum is identically zero. Maybe to correct this behaviour we have to make those boundaries real outer boundaries, i.e. we have to remove the inner part of "heizung" and make it hollow. :idea: Would that be an option?
Hi Johannes,
I tested the integration over an inner boundary in a minimal fenics example where I created the geometry and mesh in fenics and the result is the same. I get zero when simply integrating over such an inner boundary.
So there has to be a solution inside fenics.
I am thinking of using the normal vector in combination with an inner product but did not succeed until now.
Do you have any new ideas about that?
Br,
Howil
Post Reply