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 »

Hello Johannes,
joha2 wrote: Fri May 11, 2018 8:11 am 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.)
do you mean something like

Code: Select all

group_values_dict = {"positive":(0,4)}  # (default_value, marked_value)
BTW:Why do I have to save a default value for each domain/boundary?
I tested the above but it did not work for me.
I used the following code to export the mesh form FC:

Code: Select all

fileString = 'mesh_simplest_air_in_box_model.xdmf'
group_values_dict = {"positive":(0,4)}  # (default_value, marked_value)

# 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'

# 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")
    writeFenicsXDMF.write_fenics_mesh_xdmf(sel, fileString, group_values_dict)
As you can see from the following screenshot is "positive" the boundary on the positive x-axis which matches the boundary-selection with class "Right" order boundary-marker "2" (please see code below).
Screenshot from 2018-05-11 16-46-27.png
Screenshot from 2018-05-11 16-46-27.png (347.88 KiB) Viewed 2765 times

Furthermore, I do get different values when selecting directly with e.g. ds(2) and ff.ds["positive"] for integration (both should address the same boundary) EDIT:but only ds(2) is correct\EDIT. This is shown in the minimal example below (it should be ready for copy paste and run ;) EDITplease look for "assemble( (u) *ff.ds["positive"])"\EDIT):

Code: Select all

from dolfin import *

import sys
sys.path.append(r'/usr/lib/freecad/lib')
import FreeCAD as FC

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

# # Setting up and solving the model
length = 10.0  # x-dir in m, distance between plates
depth = 10.0  # y-dir in m, depth of capacitor plate
height = 10.0  # z-dir in m, heigth of capacitor plate

V1 = -3. # V, potential on first plate
V2 = 10. # V, potential on second plate
e1 = Expression(str(V1), degree=1)
e2 = Expression(str(V2), degree=1)

boundarydict = {"negative":{"type":"Dirichlet", "value":e1}, #, "marked":1
                "positive":{"type":"Dirichlet", "value":e2}, #, "marked":2
               }
ff = FacetFunctionXDMF("mesh_simplest_air_in_box_model.xdmf", boundarydict)
#ff = FacetFunctionXDMF("mesh_simplest_air_in_box_model.orig.xdmf", boundarydict)
mesh = ff.mesh

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)

# 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
left.mark(boundaries, 1)
right.mark(boundaries, 2)

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

# Define input data
a0 = Constant(1.0)  # epsilon_r 
g_L = 0  
g_R = 0  
f = 0 

# Define function space and basis functions
V = FunctionSpace(mesh, "Lagrange", 1)  # V = FunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)

bcs = ff.getDirichletBCs(V)

# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Define variational form
F = (inner(a0*grad(u), grad(v))*dx(0) 
     - g_L*v*ds(1) - g_R*v*ds(2)
     - f*v*dx(0) )

# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)

# Solve problem
u = Function(V)
solve(a == L, u,  bcs)

print("\n\n#############################################")
print("The following two values should be identical")
print("assemble( (u) *ds(2)) = "+str(assemble( u*ds(2))))
print("assemble( (u) *ff.ds['positive']) = "+str(assemble( (u) *ff.ds["positive"])))

Do you have any idea whats going wrong??
Thanks in advance,
BR,
HoWil
Attachments
mesh_simplest_air_in_box_model.xdmf.zip
(60.02 KiB) Downloaded 102 times
FC_simplest_air_in_box_model.fcstd
(123.27 KiB) Downloaded 93 times
joha2
Posts: 303
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Post by joha2 »

HoWil wrote: Fri May 11, 2018 2:59 pm Do you have any idea whats going wrong??
Hey HoWil,

I will answer in the forum, since maybe the answer is also useful for others. I analyzed my code and saw that I indeed saved the whole Measure object into the ds and dx respectively. This means since the default marked and unmarked values are 1 and 0 you have to use

Code: Select all

ds["positive"](1)
for the marked part of the boundary and

Code: Select all

ds["positive"](0)
for the unmarked part of the boundary. If you use other values for marking and unmarking the numbers are changed. So now I need a design decision. Is it more useful for the user to have the marked ds returned by ds["positive"] (i.e. ds["positive"](1)) or is it more useful to have the full measure returned such that the user might also have access to the parts which are not marked?

Looking forward for your suggestions.

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

Re: Fenics as Solver

Post by HoWil »

joha2 wrote: Sat May 19, 2018 9:25 pm I analyzed my code and saw that I indeed saved the whole Measure object into the ds and dx respectively. This means since the default marked and unmarked values are 1 and 0 you have to use

Code: Select all

ds["positive"](1)
for the marked part of the boundary and

Code: Select all

ds["positive"](0)
for the unmarked part of the boundary.
Hi Johannes,
Thank you very much :!:
I will test it this evening and will report back.
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: Sat May 19, 2018 9:25 pm This means since the default marked and unmarked values are 1 and 0 you have to use

Code: Select all

ds["positive"](1)
for the marked part of the boundary and

Code: Select all

ds["positive"](0)
for the unmarked part of the boundary.
Hi Johannes,
I tested your suggestion above and it worked for me (for boundaries). What I am still not able to do is to address domains, here the trick with (0) does not work?!
I have set-up an minimal example where I read a mesh of a asymmetric box from FC. The model basically looks like:
Screenshot from 2018-05-20 20-53-48.png
Screenshot from 2018-05-20 20-53-48.png (116.3 KiB) Viewed 2670 times
I can now integrate over the unequal faces "positive" and "negative" after I import this mesh into fenics in jupyter-notebook. Please see:
Screenshot from 2018-05-20 21-19-31.png
Screenshot from 2018-05-20 21-19-31.png (170.24 KiB) Viewed 2670 times
But addressing the domain 'air' with "assemble(1.0*cf.dx["air"](1))" does still not work.
The full given error is:

Code: Select all

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-15-0814d3517c36> in <module>()
----> 1 assemble(1.0*cf.dx["air"](1))

/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.pyc in assemble(form, tensor, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, backend)
    191 
    192     # Create dolfin Form object referencing all data needed by assembler
--> 193     dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
    194 
    195     # Create tensor

/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.pyc in _create_dolfin_form(form, form_compiler_parameters, function_spaces)
     65         return Form(form,
     66                     form_compiler_parameters=form_compiler_parameters,
---> 67                     function_spaces=function_spaces)
     68     else:
     69         raise TypeError("Invalid form type %s" % (type(form),))

/usr/lib/python2.7/dist-packages/dolfin/fem/form.pyc in __init__(self, form, form_compiler_parameters, function_spaces)
    168         subdomains = self.subdomains.get("cell")
    169         if subdomains is not None:
--> 170             self.set_cell_domains(subdomains)
    171         subdomains = self.subdomains.get("exterior_facet")
    172         if subdomains is not None:

TypeError: in method 'Form_set_cell_domains', argument 2 of type 'std::shared_ptr< dolfin::MeshFunction< std::size_t > const >'
Thank you in advance,
BR,
HoWil
joha2
Posts: 303
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Post by joha2 »

HoWil wrote: Sun May 20, 2018 7:25 pm Thank you in advance,
Hey HoWil,

I think I isolated the problem. I submitted a PR to bernds femdev branch https://github.com/berndhahnebach/FreeCAD_bhb/pull/38 which fixed the type of the CellFunction into 'size_t' instead of 'int'. Therefore for my macro code the following lines seem to work:

Code: Select all

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))
They give the following results which seem to be okay.

Code: Select all

62456800000.0
62216800000.0
240000000.0
dx(subdomain_id=everywhere)
dx(subdomain_id=1, domain=<Mesh #11>, subdomain_data=<MeshFunction of topological dimension 3 containing 279285 values>)
The problem with my implementation is that it is very difficult to use FacetFunctions and CellFunctions (and their measures) in the same functional. So I think I have to change the design slightly such that there is an underlying structure which reads the XDMFFile and saves the mesh and afterwards creates FacetFunctions and CellFunctions separately. (The mesh types in cf.mesh and ff.mesh seem to be incompatible, maybe due to topological issues). My workaround is:

Code: Select all

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

u = fe.Function(V)
fe.solve(my_form_lhs == my_form_rhs, u, ff.getDirichletBCs(V))
But this seems to be very unconvenient :-( Do you have another idea?

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

Re: Fenics as Solver

Post by HoWil »

joha2 wrote: Mon May 21, 2018 3:17 pm I think I isolated the problem. I submitted a PR to bernds femdev branch https://github.com/berndhahnebach/FreeCAD_bhb/pull/38 which fixed the type of the CellFunction into 'size_t' instead of 'int'.
Thanks for your help :!: So I will hopefully be able to to test it tomorrow.
joha2 wrote: Mon May 21, 2018 3:17 pm

Code: Select all

heizung = cf.dx["heizung_volume"](1)
But this seems to be very unconvenient :-( Do you have another idea?
Do you mean the fact that we do have to know/insert the name of the domain (here "heizung_volume") as well as the making number (here "1")?
Yes this is a bit inconvenient. I personally do not understand why I have to insert in FC a label, a marker number and a default value for each domain and boundary.
BR
HoWil
joha2
Posts: 303
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Post by joha2 »

HoWil wrote: Mon May 21, 2018 7:04 pm Yes this is a bit inconvenient. I personally do not understand why I have to insert in FC a label, a marker number and a default value for each domain and boundary.
I meant the usage of two different function spaces to implement the different contributions from the cell function and the facet function :mrgreen: :)
But yes, this is a bit unusable. I think the FC label is good for remembering the sub domains and the marker number is good for using also the complement of the marked domain. My problem is that I have too less experience with Fenics to be sure which is the most general functionality the user needs to be not restricted in performing their calculation.

Most people with Fenics experience did not yet answer to my questions concerning these usability design questions :(

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

Re: Fenics as Solver

Post by HoWil »

joha2 wrote: Mon May 21, 2018 8:44 pm I meant the usage of two different function spaces to implement the different contributions from the cell function and the facet function :mrgreen: :)
I see, yes this looks a bit inconvenient and confusing. Will test your code tomorrow since I am not sure that your commit is already in the ubuntu-daily repo (or it simply did not work for me).
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 May 21, 2018 3:17 pm My workaround is....
Dear Johannes,
I finally was able to use your example code. It took me a while until it worked and I am sure I still do not fully understand your workaround :D .
Nevertheless, I think that I found another problem :? . I think it is not possible to use inner boundaries :!:
Attached, you find a minimal example script where I use you really great! library to export a FC-model called 'FC_air_in_asym_box_model_with_si_v1p2.fcstd' to fenics. The example is written in a jupyter notebook; if this is a problem please let me know. Inside the .zip you find also a printed .pdf version.

The model looks like in the following screenshot where I need to access (/integrate over) the inner boundary called 'integration'.
Screenshot from 2018-05-26 16-18-18.png
Screenshot from 2018-05-26 16-18-18.png (225.42 KiB) Viewed 2562 times
If I select an outer face like 'negative' everything works fine! Please see also:
Screenshot from 2018-05-26 16-46-20.png
Screenshot from 2018-05-26 16-46-20.png (90.94 KiB) Viewed 2562 times
I simply integrate over the surface to get the area as validation. I basically use something like:

Code: Select all

boundary_neg = ff.ds["negative"](1)
assemble(1.0*boundary_neg)
what works fine for all other outer boundaries like 'positive' or 'negative'.
When using the same method for an inner face, the integration over 'integration' results zero without any error or warning (see In19 in the screenshot above).

Can you please test if inner boundaries are already accessible jet? :?:
Thanks in advance,
HoWil
Attachments
mesh_notebook_and_pdf.zip
(255.56 KiB) Downloaded 65 times
FC_air_in_asym_box_model_with_si_v1p2.fcstd
(1021.46 KiB) Downloaded 67 times
joha2
Posts: 303
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Post by joha2 »

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
Post Reply