Fenics as Solver

About the development of the FEM module/workbench.

Moderator: bernd

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

Re: Fenics as Solver

Postby HoWil » Sat Apr 28, 2018 7:40 pm

joha2 wrote:
Sat Apr 28, 2018 3:04 pm
As you can see the field distribution is as expected.
Thank you very much. I will try to build some model myself and (hopefully soon) report back.
BR,
HoWil
HoWil
Posts: 838
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Postby HoWil » Sun May 06, 2018 6:41 pm

For all interested.... for each 'Mesh-Group' element like 'positive' or 'negative' you have to use 'Label' as identifier for the mesh export. Otherwise you get errors like the following when trying to read/import the mesh back into python:

Code: Select all

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/usr/lib/freecad/Mod/Fem/femsolver/fenics/fenics_tools.py", line 88, in __init__
    self.readXDMFFile(xdmffilename, group_value_dict)
  File "/usr/lib/freecad/Mod/Fem/femsolver/fenics/fenics_tools.py", line 103, in readXDMFFile
    xdmffile.read(self.markers[key], key)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to open MeshFunction for reading.
*** Reason:  Mesh Grid with data Attribute not found in XDMF.
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2017.2.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------
HoWil
Posts: 838
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Postby HoWil » Sun May 06, 2018 7:13 pm

joha2 wrote:
Sat Apr 28, 2018 3:04 pm
to have two distinct surfaces where to provide boundary conditions.
Hi Johannes,
Thanks again!
I had some time to play with your script and could reproduce your example.

Am I right that until now only the export of boundaries/face-groups is implemented?
If so, any chance that you can extend your scirpt for exporting domain groups as well?
BR,
HoWil
joha2
Posts: 236
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Postby joha2 » Sun May 06, 2018 7:25 pm

HoWil wrote:
Sun May 06, 2018 7:13 pm
I had some time to play with your script and could reproduce your example.
Nice! Hopefully the interface helps also for more complicated examples.
Btw. I made a further comment within the tutorial post to use the labels of the meshgroups in the script.
HoWil wrote:
Sun May 06, 2018 7:13 pm
Am I right that until now only the export of boundaries/face-groups is implemented?
If so, any chance that you can extend your scirpt for exporting domain groups as well?
There is also an interface class for cell functions to be able to use materials or exterior sources. In one of the earlier posts I used it for demo purposes. Is there anything within those interface classes which could be improved to be more useful?

Best wishes
Johannes
HoWil
Posts: 838
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Postby HoWil » Sun May 06, 2018 8:23 pm

joha2 wrote:
Sun May 06, 2018 7:25 pm
There is also an interface class for cell functions to be able to use materials or exterior sources. In one of the earlier posts I used it for demo purposes. Is there anything within those interface classes which could be improved to be more useful?
Best wishes
Johannes
Do you mean this post https://forum.freecadweb.org/viewtopic. ... 30#p183771?
I can see that you use volumes in your dialogue but cant reproduce it.
Will try it again tomorrow.
BR,
HoWil
joha2
Posts: 236
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Postby joha2 » Mon May 07, 2018 5:55 pm

HoWil wrote:
Sun May 06, 2018 8:23 pm
Do you mean this post https://forum.freecadweb.org/viewtopic. ... 30#p183771?
I can see that you use volumes in your dialogue but cant reproduce it.
Yes, I mean this post. Unfortunately I did not provide how to do it. You have to use BooleanFragments. Maybe this can be improved in the future, by creating mesh groups by tool solids. I provide some screenshots. Hopefully they are helpful :mrgreen:. Best wishes
Johannes
Attachments
fenics_volume_group3.png
fenics_volume_group3.png (35 KiB) Viewed 645 times
fenics_volume_group2.png
fenics_volume_group2.png (68.99 KiB) Viewed 645 times
fenics_volume_group.png
fenics_volume_group.png (49.33 KiB) Viewed 645 times
HoWil
Posts: 838
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Postby HoWil » Mon May 07, 2018 7:17 pm

joha2 wrote:
Mon May 07, 2018 5:55 pm
Yes, I mean this post. Unfortunately I did not provide how to do it. You have to use BooleanFragments. Maybe this can be improved in the future, by creating mesh groups by tool solids. I provide some screenshots. Hopefully they are helpful :mrgreen:. Best wishes
Johannes
I also use BooleanFragments but have no clue how to declare the created volume-groups/domains in FEniCS.
BTW I still struggle with selecting inner solids from e.g. a booleanfragments. Anyone any new hints (beside using Arch-Section, SelectionTools orElement-Selector-macro) @Bernd?
Inner faces are no longer a problem, but inner solids are still! Please see my attached model with an domain called 'metal' which is completely surrounded by air.
Screenshot from 2018-05-07 21-12-46.png
Screenshot from 2018-05-07 21-12-46.png (412.71 KiB) Viewed 634 times
BR,
HoWil
Attachments
FC_FEniCS_3D_selection_model.fcstd
(158.07 KiB) Downloaded 19 times
joha2
Posts: 236
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Postby joha2 » Tue May 08, 2018 6:45 am

HoWil wrote:
Mon May 07, 2018 7:17 pm
I also use BooleanFragments but have no clue how to declare the created volume-groups/domains in FEniCS.
Hey HoWil,
hopefully this inner solid thing gets solved, this is really necessary for field theoretic computations.
Here is some example, how to use the CellExpressionXDMF for defining a source. `rho` can be used like a normal expression in Fenics, because the necessary functions were overloaded. What is still missing: implementation of vector valued sources and some other things to improve the flexibility of the code.

Code: Select all

rho = CellExpressionXDMF("file.xdmf",
                         {"negative": lambda x: -1., "positive": lambda x: 1.},
                        default=0., degree=0)
The lambdas make the code quite flexible, but any other function is also permitted (e.g. for non equally distributed charges). Notice that also here the labels of the mesh groups are used. The CellExpressionXDMF structure could also be used for defining materials.

Best wishes
Johannes
HoWil
Posts: 838
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Postby HoWil » Thu May 10, 2018 7:23 pm

joha2 wrote:
Sat Apr 28, 2018 3:04 pm

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
"""

e1 = fe.Expression("1", degree=1)
e2 = fe.Expression("-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", 
	{"positive":{"type":"Dirichlet", "value":e1},
	"negative":{"type":"Dirichlet", "value":e2}})

"""
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)

u = fe.TrialFunction(V)
v = fe.TestFunction(V)
f = fe.Constant(0.0)
a = fe.inner(fe.grad(u), fe.grad(v))*fe.dx
L = f*v*fe.dx

u = fe.Function(V)
fe.solve(a == L, u, ff.getDirichletBCs(V))

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


fe.File("sol_fc_fenics_poisson.pvd") << u
Hi Johannes,
Is there a way to set custom markers to e.g. the faces called 'positive' and 'negative'?
Usually I do the follwoing ...

Code: Select all

# Create classes for defining parts of the boundaries and the interior
# of the domain
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], length)

# Initialize sub-domain instances
left = Left()
right = Right()

# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
but whant to combine this with your code without declaring 'left' and 'right' again.
I do need markes for integrating with e.g. the marker 2 with "ds(2)" along an boundary.
Thanks in advance,
HoWil
joha2
Posts: 236
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Postby joha2 » Fri May 11, 2018 8:11 am

HoWil wrote:
Thu May 10, 2018 7:23 pm
Hi Johannes,
Is there a way to set custom markers to e.g. the faces called 'positive' and 'negative'?
Usually I do the follwoing ...
<SCHNIPP>
<SCHNAPP>
but whant to combine this with your code without declaring 'left' and 'right' again.
I do need markes for integrating with e.g. the marker 2 with "ds(2)" along an boundary.
Thanks in advance,
HoWil
Heyhey HoWil,

you may use the export function with a further dict parameter,

Code: Select all

write_fenics_mesh_xdmf(fem_mesh_obj, fileString, group_values_dict)
The group_values_dict consists of pairs of (default_value, marked_value) for every key (corresponding to the label of that group).
(my default choice was 1 for marked and 0 for not marked.) Within the FacetFunctionXDMF you also have access to the mesh, the markers, the boundary conditions (only Dirichlet at the moment) and the integral measures via

Code: Select all

ff.mesh # mesh
ff.markers["group"] # markers
ff.ds["group"] # measure
ff.bcs["group"] # boundary conditions
Same story for CellExpressionXDMF:

Code: Select all

cf.mesh # mesh
cf.markers["group"] # markers
cf.dx["group"] # measure
Both classes need a new readXDMFfile call for changing those dicts. For further reference, please have a look at
https://github.com/FreeCAD/FreeCAD/blob ... s_tools.py
There are lots of lots of TODOs in there, which could be improved in the future to render the Fenics interface more useful, but I have not enough experience to make a useful decision at some points. I also have no useful default case for overlapping cell functions or facet functions.

For your left and right problem I would suggest to define groups in the mesh and use the default group_values_dict (1 for markes, 0 else). After that you may use

Code: Select all

ff.markers["left"]
ff.markers["right"]
ff.ds["left"]
ff.ds["right"]
to access the appropriate measures and markers.

Hope this helps.
Best wishes
Johannes