Hello,
I’m trying to link a node to face of a 3D geometry with 1D elements.
I succeed to creat “edges” with netgen as you can see on the picture.
However, I don’t understand how to convert these netgen “edges” into ngsolve element, and after use it in my FES.
Do you know how to do this ?
Thanks.
Thomas
My code :
from netgen.occ import *
from ngsolve import *
from ngsolve import *
from netgen.csg import Pnt
from netgen.meshing import Element1D, MeshPoint
import math
# Variable
E = 212000
nu = 0.3
Rho = 7849e-9
# Dirichlet Boundaries
Face_Dirichletx = 'Face12'
Face_Dirichlety = 'Face12'
Face_Dirichletz = 'Face12'
# Displacement boundary
FaceDisplacement = 'Face5'
# Degree of Element
EleDeg = 2
tolerance = 0.01
def arrondi(valeur):
if abs(valeur) < 0.1:
return 0.0
else:
return valeur
def vec_norm(r):
norm = math.sqrt(r[0] ** 2 + r[1] ** 2 + r[2] ** 2)
return norm
def vec_somme(rA, rB):
xsom = rA[0] + rB[0]
ysom = rA[1] + rB[1]
zsom = rA[2] + rB[2]
som = [xsom, ysom, zsom]
return som
def prod_vec(rA, rB):
xm = rA[1] * rB[2] - rA[2] * rB[1]
ym = rA[2] * rB[0] - rA[0] * rB[2]
zm = rA[0] * rB[1] - rA[1] * rB[0]
m = [xm, ym, zm]
return m
def dist_line(rp, m, e):
dist = vec_norm(vec_somme(m, prod_vec(e, rp))) / vec_norm(e)
return dist
def dist_hole(rA, rB, rp):
xe = rB[0] - rA[0]
ye = rB[1] - rA[1]
ze = rB[2] - rA[2]
e = [xe, ye, ze]
m = prod_vec(rA, rB)
d = dist_line(rp, m, e)
return d
shp_path = "TestVisM8.stp"
# center_points, ext_points = get_centers(shp)
center_points = [[0.0, 0.0, 0.0, 8.0]]
# print(center_points)
ext_points = [[0.0, 0.0, 10.0], [0.0, 0.0, -10.0]]
# print(ext_points)
geometry = OCCGeometry(shp_path)
mesh = Mesh(geometry.GenerateMesh(maxh=10)).Curve(3)
n_center = len(center_points)
N = len(mesh.vertices)
for i in range(n_center):
l = i + 1
indice = mesh.ngmesh.Add(MeshPoint(Pnt(center_points[i][0], center_points[i][1], center_points[i][2])))
idx = mesh.ngmesh.AddRegion("Center_point_" + str(l), dim=1)
j = 2 * i
k = 2 * i + 1
print(center_points[i])
print(ext_points[j], ext_points[k])
for v in mesh.vertices:
rp = list(v.point)
rA = ext_points[j]
rB = ext_points[k]
d = dist_hole(rA, rB, rp)
if d <= (0.5 * center_points[i][-1]) + tolerance:
mesh.ngmesh.Add(Element1D([indice, v.nr + 1], index=idx))
# print (v, v.nr, v.point)
# print (v, v.nr, v.point)
mesh = Mesh(mesh.ngmesh)
# FE Space ---------------------------------------------------------------------------------
fes = VectorH1(mesh, order=EleDeg, dirichletx=Face_Dirichletx, dirichlety=Face_Dirichlety, dirichletz=Face_Dirichletz)
u, v = fes.TnT()
gfu = GridFunction(fes)
The attached STEP file :
TestVisM8.stp (16.4 KB)