Making FEM mesh dependent only on CharacteristicLength

About the development of the FEM module/workbench.

Moderator: bernd

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

Making FEM mesh dependent only on CharacteristicLength

Post by orxshi » Sun Dec 29, 2019 11:38 am

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

sys.path.append('/usr/lib/freecad/lib/')

import FreeCAD
import Draft
import Part
import BOPTools.JoinFeatures

doc = FreeCAD.newDocument('newdoc')

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 = ObjectsFem.makeMeshGmsh(FreeCAD.ActiveDocument, 'FEMMeshGmsh')
mesh.ElementDimension = 3
FreeCAD.ActiveDocument.ActiveObject.Part = FreeCAD.ActiveDocument.Shape

mr_fus = ObjectsFem.makeMeshRegion(FreeCAD.ActiveDocument, FreeCAD.ActiveDocument.FEMMeshGmsh, 0.05, 'fus')

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()
User avatar
Kunda1
Posts: 7070
Joined: Thu Jan 05, 2017 9:03 pm

Re: Making FEM mesh dependent only on CharacteristicLength

Post by Kunda1 » Sun Dec 29, 2019 12:18 pm

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: 1175
Joined: Tue Jan 03, 2017 8:42 pm

Re: Making FEM mesh dependent only on CharacteristicLength

Post by UR_ » Mon Dec 30, 2019 10:07 am

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

Who is online

Users browsing this forum: No registered users and 4 guests