How to assess CoefficientFunction at Mesh Vertices

Main Forum Old-Forum - do not post here
1

Hi Guys,

I am new here, and would like to ask a simple question.

I want to use MatPlotLib for visualization of displacement and stress in an elastic deformation simulation.
Although the GridFunction seems to have 6 times more vectors than the Mesh has vertices, I figured that the first 1/6th of GridFunction vec data corresponds to the displacement of Mesh vertices.

Next, I want to compute the Stress which is defined as a coefficient function. However, I cannot figure out how to input the MeshPoints in the correct format, can someone help me finish writing the last few lines of this code?

# Load CAD file as geometry
from netgen.occ import *
geo = OCCGeometry("Rail.stp")

# Generate mesh
from ngsolve import *
gmsh = geo.GenerateMesh()
msh = Mesh(gmsh)
print(msh.GetBoundaries())
print(msh.GetMaterials())

# Material parameters
E, nu = 2.1e11, 0.2

# Lame parameters
mu  = 1/2*(E / (1+nu))
lam = E * nu / ((1+nu)*(1-2*nu))

# Linear strain tensor
def epsilon(u):
     return 0.5 * (u.Deriv().trans + u.Deriv())

# linear stress tensor
def stress(u):
     return 2*mu*epsilon(u) + lam*Trace(epsilon(u))*Id(msh.dim)

# Set fixed boundary
fes = H1(msh, order=2, dirichlet="bc_1", dim=msh.dim)
u,v = fes.TrialFunction(), fes.TestFunction()
 
# Define problem
a = BilinearForm(fes)
a += SymbolicBFI( 2*mu*InnerProduct(epsilon(u),epsilon(v)) + lam*Trace(u.Deriv())*Trace(v.Deriv()))
a.Assemble()

# Add load
force = CoefficientFunction( (0,0,-1e6) )

# Solve
f = LinearForm(fes)
f += SymbolicLFI( force*v, definedon=msh.Boundaries("bc_0"))
f.Assemble()

# Get displacement
disp = GridFunction(fes,name="displacement")
disp.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec

# Export deformed mesh as lists of vertices/triangles
i = 0; scale = 1000;
x = []; y = []; z = [];
for v in msh.vertices:
    x.append(v.point[0] + scale*disp.vec[i][0])
    y.append(v.point[1] + scale*disp.vec[i][1])
    z.append(v.point[2] + scale*disp.vec[i][2])
    i += 1
vertices = [x,y,z]

triangles = []
for e in msh.Elements(BND):
    triangles.append([e.vertices[0].nr, e.vertices[1].nr, e.vertices[2].nr])

# How to compute stress at mesh nodes??? (strs is a coefficient function, not a grid function)
strs = stress(disp)

Attachment: Rail.stp

2

Ah, I think I might have figures it out!
Could someone kindly confirm that the following is correct way to compute Von Mises:

# Get Displacement
disp = GridFunction(fes,name="displacement")
disp.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec
strs = stress(disp)

# Export deformed mesh and Von Mises stress
i = 0; scale = 1000;
x = []; y = []; z = []; mises = [];
for v in msh.vertices:
    x.append(v.point[0] + scale*disp.vec[i][0])
    y.append(v.point[1] + scale*disp.vec[i][1])
    z.append(v.point[2] + scale*disp.vec[i][2])
    s = strs(msh(v.point[0],v.point[1],v.point[2]))
    mises.append(sqrt(s[0]**2 + s[4]**2 + s[8]**2 - s[0]*s[4] - s[4]*s[8] - s[8]*s[0] + 3*(s[1]**2 + s[5]**2 + s[6]**2)))
    i += 1
3

Hi,

as you figured out, you can evaluate CoefficientFunctions in MeshPoints, which know the element number, and the point on the reference element (i.e. the barycentric coordinates).
The function msh(x,y,z) searches for one element containing the point (which is not unique).

You can loop over elements, and map points from the reference element to the physical element, and evaluate the CoefficientFunction there. This avoids the search.

Points on the reference element are stored in an IntegrationRule object. The 4 corner points of a tet, and 4 irrelevant weights:

ir_tet = IntegrationRule( [(1,0,0), (0,1,0), (0,0,1), (0,0,0)], [0,0,0,0] )

We loop over the volume elements, and map the points to the physical element:

for el in msh.Elements(VOL):
    elpnts = msh.GetTrafo(el)(ir_tet)
    print (strs(elpnts))

Since this loop is in Python it is not the fastest. We have a vectorized version of the loop using numpy arrays:

pts = msh.MapToAllElements(ir_tet, VOL)
pointvalues = strs(pts)

Best, Joachim

4

Thanks for this extra tips Joachim, I will keep this in mind.

In the end, I decided to use Plotly for plotting, and final result is very pleasing (see attached screenshot).

I am including below the cleaned-up code, in case it ends up being useful to someone else…

# Load CAD file as geometry
from netgen.occ import *
geo = OCCGeometry("Drawings/Rail.stp")

# Generate mesh
from ngsolve import *
gmsh = geo.GenerateMesh()
msh = Mesh(gmsh)
print(msh.GetBoundaries())
print(msh.GetMaterials())

# Material parameters
E, nu = 2.1e11, 0.2

# Lame parameters
mu  = 1/2*(E / (1+nu))
lam = E * nu / ((1+nu)*(1-2*nu))

# Linear strain tensor
def epsilon(u):
     return 0.5 * (u.Deriv().trans + u.Deriv())

# linear stress tensor
def stress(u):
     return 2*mu*epsilon(u) + lam*Trace(epsilon(u))*Id(msh.dim)
 
# Von mises stress
def vonmises(s):
    return sqrt(s[0]**2 + s[4]**2 + s[8]**2 - 
                s[0]*s[4] - s[4]*s[8] - s[8]*s[0] + 
                3*(s[1]**2 + s[5]**2 + s[6]**2))

# Define fixed boundary
fes = H1(msh, order=2, dirichlet="bc_1", dim=msh.dim)
u,v = fes.TrialFunction(), fes.TestFunction()
 
# Define PDE
a = BilinearForm(fes)
a += SymbolicBFI( 2*mu*InnerProduct(epsilon(u),epsilon(v)) + lam*Trace(u.Deriv())*Trace(v.Deriv()))
a.Assemble()

# Define load
force = CoefficientFunction( (0,0,-1e6) )

# Solve
f = LinearForm(fes)
f += SymbolicLFI( force*v, definedon=msh.Boundaries("bc_0"))
f.Assemble()

# Get Displacement
disp = GridFunction(fes,name="displacement")
disp.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec
strs = stress(disp)

# Extract deformed mesh and Von Mises stress
scale = 1000;
x = []; y = []; z = []; mises = [];
for v in msh.vertices:
    # Find mesh location
    vMsh = msh(v.point[0],v.point[1],v.point[2])
    dispVec = disp(vMsh)
    stressTsr = strs(vMsh)
    x.append(v.point[0] + scale*dispVec[0])
    y.append(v.point[1] + scale*dispVec[1])
    z.append(v.point[2] + scale*dispVec[2])
    mises.append(vonmises(stressTsr))

v1 = []; v2 = []; v3 =[];
for e in msh.Elements(BND):
    v1.append(e.vertices[0].nr)
    v2.append(e.vertices[1].nr)
    v3.append(e.vertices[2].nr)

# Plot Data
from plotly.subplots import make_subplots
import plotly.graph_objects as go

fig = make_subplots(rows=1, cols=1,specs=[[{'type': 'surface'}]])

fig.add_trace(go.Mesh3d(x=x, y=y, z=z, i=v1, j=v2, k=v3, intensity = mises, colorbar_title='Von Mises',
                colorscale=[[0, 'mediumturquoise'], [0.5, 'gold'], [1, 'red']],
                intensitymode='vertex', name='y', showscale=True),row=1,col=1)

fig.show()

800×456 73.6 KB