Fenics as Solver

About the development of the FEM module/workbench.

Moderator: bernd

looo
Posts: 2655
Joined: Mon Nov 11, 2013 5:29 pm

Re: Fenics as Solver

Postby looo » Wed Dec 14, 2016 6:10 pm

JeffWitz wrote:I think that it is not the good way to proceed as FreeCAD seems to go to multi-physics simulations.
It's only a test, to understand the pipeline. For sure there has to be a genericResult object too.
Dirichlet : set the quantities we want to compute at given nodes (First type),
Neumann : Set the derivative of the quantity at the nodes (Second type),
Cauchy : set both of them (Third type),
Robin : set both of them with a linear law (Third type).
I think there will be use-cases which need other kind of bc. So better try to do it in a more generic way.

@bernd thanks for this example. But I still do not understand why it isn't working in my case.
Consider this example. I have an analysis with a Box_Mesh. Now I want to add a result and set the temperature:

Code: Select all

result = FreeCAD.ActiveDocument.addObject('Fem::FemResultObject', "result")
App.ActiveDocument.Analysis.Member += [result]
result.Mesh = App.ActiveDocument.Box_Mesh
result.NodeNumbers = App.ActiveDocument.Box_Mesh.FemMesh.NodeCount
result.Temperature = range(App.ActiveDocument.Box_Mesh.FemMesh.NodeCount)
But no result is shown for the temperature. The message is:

Code: Select all

FEM: Result node numbers are not equal to FEM Mesh NodeCount.
Do I have to set all the other values too?
User avatar
bernd
Posts: 8021
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland

Re: Fenics as Solver

Postby bernd » Wed Dec 14, 2016 7:07 pm

looo wrote:But I still do not understand why it isn't working in my case
try this

start FreeCAD --> startwb --> load FEM 3D example --> activate the analysis

Code: Select all

result = FreeCAD.ActiveDocument.addObject('Fem::FemResultObject', "result")
App.ActiveDocument.Analysis.Member += [result]
result.Mesh = App.ActiveDocument.Box_Mesh
result.NodeNumbers =  result.Mesh.FemMesh.Nodes.keys()
tmp = [0] * len(result.NodeNumbers)
tmp[0] = 10
result.Temperature = tmp

EDIT: update code to run on 0.18.14709

Code: Select all

import ObjectsFem
result = ObjectsFem.makeResultMechanical(App.ActiveDocument, 'MyResult')
App.ActiveDocument.Analysis.Group += [result]
App.ActiveDocument.recompute()

result.Mesh = App.ActiveDocument.Box_Mesh
result.NodeNumbers = [int(key) for key in result.Mesh.FemMesh.Nodes.keys()]  # keys is list of long
tmp = [0] * len(result.NodeNumbers)
tmp[0] = 10
result.Temperature = tmp
App.ActiveDocument.recompute()

double click on result --> activate temperature --> see attached screen
screen.jpg
screen.jpg (289.7 KiB) Viewed 1334 times
looo
Posts: 2655
Joined: Mon Nov 11, 2013 5:29 pm

Re: Fenics as Solver

Postby looo » Wed Dec 14, 2016 9:17 pm

thanks. This line was the important one:

Code: Select all

result.NodeNumbers = App.ActiveDocument.Box_Mesh.FemMesh.NodeCount
Now I have plugged in a solver I have written some time ago and set the displacement. Nearly worked straight forward. Visualizing the displacement gave me a little error.

Code: Select all

Traceback (most recent call last):
  File "/home/lo/anaconda/envs/fc_gen_sol/Mod/Fem/_TaskPanelShowResult.py", line 167, in z_displacement_selected
    self.select_displacement_type("U3")
  File "/home/lo/anaconda/envs/fc_gen_sol/Mod/Fem/_TaskPanelShowResult.py", line 259, in select_displacement_type
    (minm, avg, maxm) = self.get_result_stats(disp_type)
  File "/home/lo/anaconda/envs/fc_gen_sol/Mod/Fem/_TaskPanelShowResult.py", line 134, in get_result_stats
    match_table = {"U1": (Stats[0], Stats[1], Stats[2]),
IndexError: list index out of range
torsion.png
torsion.png (432.68 KiB) Viewed 1318 times
User avatar
bernd
Posts: 8021
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland

Re: Fenics as Solver

Postby bernd » Thu Dec 15, 2016 6:44 am

something seams wrong with the stats in your result object. Might be you do not have one and we do not assume this could gone happen. Could you post a simple example file and some code to reproduce?
looo
Posts: 2655
Joined: Mon Nov 11, 2013 5:29 pm

Re: Fenics as Solver

Postby looo » Thu Dec 15, 2016 7:51 pm

Yes stats are not set. Only displacement is filled...
Maybe we need another kind of result-object. I would like to have a result object which simple has a dict with all the result values. One node-based and one element based. (the element based isn't that important, but it is needed for some kind of analysis). The taskpanel for this result-object then has a dropdownlist (combobox) where the result for visualization can be selected. (Similar to the way how paraview works). But maybe this is more difficult as it involves c++ development...

I will try to post some examples soon.
HoWil
Posts: 831
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Postby HoWil » Thu Dec 15, 2016 8:48 pm

looo wrote:thanks. This line was the important one:

Code: Select all

result.NodeNumbers = App.ActiveDocument.Box_Mesh.FemMesh.NodeCount
Now I have plugged in a solver I have written some time ago and set the displacement. Nearly worked straight forward. Visualizing the displacement gave me a little error.

Code: Select all

Traceback (most recent call last):
  File "/home/lo/anaconda/envs/fc_gen_sol/Mod/Fem/_TaskPanelShowResult.py", line 167, in z_displacement_selected
    self.select_displacement_type("U3")
  File "/home/lo/anaconda/envs/fc_gen_sol/Mod/Fem/_TaskPanelShowResult.py", line 259, in select_displacement_type
    (minm, avg, maxm) = self.get_result_stats(disp_type)
  File "/home/lo/anaconda/envs/fc_gen_sol/Mod/Fem/_TaskPanelShowResult.py", line 134, in get_result_stats
    match_table = {"U1": (Stats[0], Stats[1], Stats[2]),
IndexError: list index out of range
torsion.png
Nice!
Hi looo,
Is it possible for you to provide a fenics example equivalent to the shipped calculix cantilever? This would help others testing the current available code.
Br
Howil
looo
Posts: 2655
Joined: Mon Nov 11, 2013 5:29 pm

Re: Fenics as Solver

Postby looo » Thu Dec 15, 2016 9:45 pm

yes, at some point we will try to do that. But I am not very deep into fenics. I think we will first try to do some simpler equations (really simple ones like shown in the picture). Getting the formulations for the full 3d linear continuum mechanics will involve some more work.
User avatar
bernd
Posts: 8021
Joined: Sun Sep 08, 2013 8:07 pm
Location: Zürich, Switzerland

Re: Fenics as Solver

Postby bernd » Mon Dec 19, 2016 3:15 pm

looo wrote:Yes stats are not set. Only displacement is filled...
Maybe we need another kind of result-object. I would like to have a result object which simple has a dict with all the result values. One node-based and one element based. (the element based isn't that important, but it is needed for some kind of analysis). The taskpanel for this result-object then has a dropdownlist (combobox) where the result for visualization can be selected. (Similar to the way how paraview works). But maybe this is more difficult as it involves c++ development...

I will try to post some examples soon.
sounds like an very interesting approach. Keep in mind in addition to the FreeCAD FEM Result TaskPanel the FreeCAD FEM VTK result plot tools are based on the FreeCAD FEM result object too. The Result Task Panel is very good to get a fast Result overview but if one really starts working with the results the vtk tools are indispensable.
joha2
Posts: 230
Joined: Tue Oct 11, 2016 9:48 pm

Re: Fenics as Solver

Postby joha2 » Fri Dec 30, 2016 12:40 pm

For the linear elasticity problem we may only need to adapt the code from the fenics tutorial given in
https://fenicsproject.org/pub/tutorial/ ... ftut:elast
It seems that it is not too different from the calculix linear elasticity test case.
HoWil
Posts: 831
Joined: Sun Jun 14, 2015 7:31 pm
Location: Austria

Re: Fenics as Solver

Postby HoWil » Fri Dec 30, 2016 8:20 pm

joha2 wrote:For the linear elasticity problem we may only need to adapt the code from the fenics tutorial given in
https://fenicsproject.org/pub/tutorial/ ... ftut:elast
It seems that it is not too different from the calculix linear elasticity test case.
Dear joha2,
I also found this example after my last post here and did adopt it with some parts from other examples from over here (https://github.com/FEniCS/dolfin/blob/m ... sticity.py for pvd export, ...).

In the provided example is 'gravity' induced force on the complete volume used and I tried to change this to the calculix-example used in FC.
But I struggled in the end with introducing a point load at the end of the beam.

For the 'gravity' or 'self weight' problem do match the results from calculix-FC and the Fenics script well (fine mesh with 'scalefactor = 1'), even if the solving time is much much higher for fenics.

Displacement from Fenics:

Code: Select all

Solving linear variational problem.
min/max u: -3.97926569654e-09 0.00207006247747
And FC gives:
2.24684 mm

Since I do not have time yet to further look into it, I post the last version of my file hoping it will help.

BR,
HoWil

Code: Select all

"""
FEniCS tutorial demo program: Linear elastic problem.

  -div(sigma(u)) = f

The model is used to simulate an elastic beam clamped at
its left end and deformed under its own weight.
"""

from __future__ import print_function
from fenics import *

# Scaled variables
Len = 8; W = 1; H = 1 # dimensions in m

#mu = 1
rho = 7900   # kg/m^3 has to be adopted


# Elasticity parameters
E = 210.0e9   # Pa
nu = 0.3   # dimensionless
mu = E/(2.0*(1.0 + nu))
lambda_ = E*nu/((1.0 + nu)*(1.0 - 2.0*nu))


#delta = W/Len
#gamma = 0.4*delta**2
#beta = 1.25
#lambda_ = beta
g = 9.81 # m/s^2 gamma

# Create mesh and define function space
scalefactor = 1
mesh = BoxMesh(Point(0, 0, 0), Point(Len, W, H), 30/scalefactor, 15/scalefactor, 15/scalefactor)
V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1E-14

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

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

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g)) 
T = Constant((0, 0, -10000000))   # force vector # does not work yet!!!!!
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds(0) # self load + load on face

Fz = 100000.
#ps = PointSource(V.sub(0),Point(Len, 0., 0.), Fz) # !!!!!!! does not work

#delta = PointSource(V, Point(0., 0., 0.), 0.)
#delta.apply(b)


# Compute solution
u = Function(V)
#ps.apply(u.vector()) # !!!!!!! does not work
solve(a == L, u, bc)

# Plot solution
plot(u, title='Displacement', mode='displacement')

# Plot stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  # deviatoric stress
von_Mises = sqrt(3./2*inner(s, s))
V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)
plot(von_Mises, title='Stress intensity')

# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
      u_magnitude.vector().array().min(),
      u_magnitude.vector().array().max())

# Save solution to file in VTK format
File('elasticity/displacement.pvd') << u
File('elasticity/von_mises.pvd') << von_Mises
File('elasticity/magnitude.pvd') << u_magnitude

# Hold plot
interactive()

Attachments
FC_fenics_demo.fcstd
(29.62 KiB) Downloaded 34 times