## Making FEM mesh dependent only on CharacteristicLength

About the development of the FEM module/workbench.

Moderator: bernd

orxshi
Posts: 15
Joined: Thu Dec 19, 2019 4:50 pm

### Making FEM mesh dependent only on CharacteristicLength

I would like to make FEM mesh to be completely dependent on CharacteristicLength but it also depends on number of lines/points used in construction of geometry. In the script, I increase the last arguments of linspace

Code: Select all

``````xm = numpy.linspace(0.00001, 2, 10)
xp = numpy.linspace(0.40001, 1.018, 10)
p = numpy.linspace(0, 2*math.pi, 10)``````
in order to make the geometry look smoother. So making the geometry finer makes FEM mesh finer as well. I want to control mesh intensity only with CharacteristicLength which is 0.05 in this case:

Code: Select all

``mr_fus = ObjectsFem.makeMeshRegion(FreeCAD.ActiveDocument, FreeCAD.ActiveDocument.FEMMeshGmsh, 0.05, 'fus')``
Is this possible?

Code: Select all

``````import sys
import math
import numpy

import Draft
import Part
import BOPTools.JoinFeatures

ZERO = 1e-10

def U(c, xl):

u = c[5]

if c[6] != 0:

t = c[0] + c[1] * abs((xl + c[2]) / c[3]) ** c[4]
if abs(t) <= ZERO:
t = 0

u += c[6] * t ** (1./c[7])

return u;

class Fuselage:

def H(xl):

if xl < 0.4:
c = [1.0, -1.0, -0.4, 0.4, 1.8, 0.0, 0.25, 1.8]
elif xl < 0.8:
c = [0.0, 0.0, 0.0, 0.0, 0.0, 0.25, 0.0, 0.0]
elif xl < 1.9:
c = [1.0, -1.0, -0.8, 1.1, 1.5, 0.05, 0.2, 0.6]
elif xl < 2.0:
c = [1.0, -1.0, -1.9, 0.1, 2.0, 0.0, 0.05, 2.0]

return U(c, xl);

def W(xl):

if xl < 0.4:
c = [1.0, -1.0, -0.4, 0.4, 2.0, 0.0, 0.25, 2.0]
elif xl < 0.8:
c = [0.0, 0.0, 0.0, 0.0, 0.0, 0.25, 0.0, 0.0]
elif xl < 1.9:
c = [1.0, -1.0, -0.8, 1.1, 1.5, 0.05, 0.2, 0.6]
elif xl < 2.0:
c = [1.0, -1.0, -1.9, 0.1, 2.0, 0.0, 0.05, 2.0]

return U(c, xl);

def Z(xl):

if xl < 0.4:
c = [1.0, -1.0, -0.4, 0.4, 1.8, -0.08, 0.08, 1.8]
elif xl < 0.8:
c = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
elif xl < 1.9:
c = [1.0, -1.0, -0.8, 1.1, 1.5, 0.04, -0.04, 0.6]
elif xl < 2.0:
c = [0.0, 0.0, 0.0, 0.0, 0.0, 0.04, 0.0, 0.0]

return U(c, xl);

def N(xl):

if xl < 0.4:
c = [2.0, 3.0, 0.0, 0.4, 1.0, 0.0, 1.0, 1.0]
elif xl < 0.8:
c = [0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0]
elif xl < 1.9:
c = [5.0, -3.0, -0.8, 1.1, 1.0, 0.0, 1.0, 1.0]
elif xl < 2.0:
c = [2.0, 0.0, 0.0, 0.0, 0.0, 0.04, 1.0, 1.0]

return U(c, xl);

class Pylon:

def H(xl):

if xl < 0.8:
c = [1.0, -1.0, -0.8, 0.4, 3.0, 0.0, 0.145, 3.0]
elif xl < 1.018:
c = [1.0, -1.0, -0.8, 0.218, 2.0, 0.0, 0.145, 2.0]

return U(c, xl);

def W(xl):

if xl < 0.8:
c = [1.0, -1.0, -0.8, 0.4, 3.0, 0.0, 0.166, 3.0]
elif xl < 1.018:
c = [1.0, -1.0, -0.8, 0.218, 2.0, 0.0, 0.166, 2.0]

return U(c, xl);

def Z(xl):

if xl < 0.4:
c = [0.0, 0.0, 0.0, 0.0, 0.0, 0.125, 0.0, 0.0]
elif xl < 1.018:
c = [1.0, -1.0, -0.8, 1.1, 1.5, 0.065, 0.06, 0.6]

return U(c, xl);

def N(xl):

if xl < 0.4:
c = [0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0]
elif xl < 1.018:
c = [0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0]

return U(c, xl);

def yl(xl, phi, part):
return r(part.H(xl), part.W(xl), part.N(xl), phi) * math.sin(phi)

def zl(xl, phi, part):
return r(part.H(xl), part.W(xl), part.N(xl), phi) * math.cos(phi) + part.Z(xl)

def r(H, W, N, phi):
a = abs(0.5 * H * math.sin(phi)) ** N + abs(0.5 * W * math.cos(phi)) ** N
b = (0.25 * H * W) ** N
return (b / a) ** (1./N)

xm = numpy.linspace(0.00001, 2, 10)
xp = numpy.linspace(0.40001, 1.018, 10)
p = numpy.linspace(0, 2*math.pi, 10)

def makepart(part, x):
polygons = []
for i in range(len(x)-1):
points = []
for j in range(len(p)-1):
points.append(FreeCAD.Vector(x[i], yl(x[i], p[j], part), zl(x[i], p[j], part)))
points.append(points[0])
polygons.append(Part.makePolygon(points))

loft = Part.makeLoft(polygons)
cap1 = Part.Face(polygons[0])
cap2 = Part.Face(polygons[-1])
shell = Part.Shell(loft.Faces+[cap1, cap2])
return shell

fuselage = makepart(Fuselage, xm)
pylon = makepart(Pylon, xp)

# make compound
heli = Part.makeCompound([fuselage, pylon])

# make heli a solid
s = Part.Solid(heli)
Part.show(heli)

import ObjectsFem
mesh.ElementDimension = 3

temp = []
for i in range(1,len(App.ActiveDocument.Shape.Shape.Faces)):
temp.append((App.ActiveDocument.Shape, 'Face' + str(i)))

mr_fus.References = temp

import femmesh.gmshtools as gmshtools
gmsh_mesh = gmshtools.GmshTools(mesh)
gmsh_mesh.create_mesh()``````
Kunda1
Posts: 6440
Joined: Thu Jan 05, 2017 9:03 pm

### Re: Making FEM mesh dependent only on CharacteristicLength

Moved to FEM subforum
Want to contribute back to FC? Checkout:
#lowhangingfruit | Use the Source, Luke. | How to Help FreeCAD | How to report FC bugs and features
UR_
Posts: 1153
Joined: Tue Jan 03, 2017 8:42 pm

### Re: Making FEM mesh dependent only on CharacteristicLength

Your script failed to make a solid

Code: Select all

``````# make heli a solid
s = Part.Solid(heli)
Part.show(heli)
``````

newdoc.Shape:
VERTEX : 36
EDGE : 54
WIRE : 22
FACE : 22
SHELL : 2
SOLID : 0
COMPSOLID : 0
COMPOUND : 1
SHAPE : 137

because there are many BOPCheck errors listed by Part:CheckGeometry

newdoc.Shape : Face1 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face22 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Edge4 : Edge : BOPAlgo SelfIntersect
newdoc.Shape : Face22 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face1 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face12 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face9 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face18 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Edge29 : Edge : BOPAlgo SelfIntersect
newdoc.Shape : Face1 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Edge48 : Edge : BOPAlgo SelfIntersect
newdoc.Shape : Face9 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Edge49 : Edge : BOPAlgo SelfIntersect
newdoc.Shape : Face9 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face9 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face19 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face1 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face13 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Edge51 : Edge : BOPAlgo SelfIntersect
newdoc.Shape : Face9 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Edge33 : Edge : BOPAlgo SelfIntersect
newdoc.Shape : Face1 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face9 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face20 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face1 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face14 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Edge4 : Edge : BOPAlgo SelfIntersect
newdoc.Shape : Edge31 : Edge : BOPAlgo SelfIntersect
newdoc.Shape : Face9 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Face22 : Face : BOPAlgo SelfIntersect
newdoc.Shape : Edge37 : Edge : BOPAlgo SelfIntersect
newdoc.Shape : Face1 : Face : BOPAlgo SelfIntersect

Should be fixed before Gmsh being fed.

### Who is online

Users browsing this forum: Bing [Bot] and 4 guests