Finite Element Mesh#
|
In this example we will demonstrate how to:
|
from dlubal.api import rfem, common
from math import inf
import re
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
# -------------------------------------------------------
# This example demonstrates creating a simple slab with loading,
# retrieving mesh data via get_mesh_table() and visualising it.
# -------------------------------------------------------
def define_structure() -> list:
# Simple rectangular slab (6 x 4 m) with line support along one edge
# and nodal hinged supports at the remaining corners.
return [
# Materials
rfem.structure_core.Material(
no=1,
name="C25/30",
),
# Thickness
rfem.structure_core.Thickness(
no=1,
material=1,
uniform_thickness=0.250
),
# Nodes
rfem.structure_core.Node(no=1, coordinate_1=0, coordinate_2=0),
rfem.structure_core.Node(no=2, coordinate_1=0, coordinate_2=4),
rfem.structure_core.Node(no=3, coordinate_1=3, coordinate_2=4),
rfem.structure_core.Node(no=4, coordinate_1=6, coordinate_2=4),
rfem.structure_core.Node(no=5, coordinate_1=6, coordinate_2=0),
# Lines
rfem.structure_core.Line(no=1, definition_nodes=[1, 2]),
rfem.structure_core.Line(no=2, definition_nodes=[2, 3]),
rfem.structure_core.Line(no=3, definition_nodes=[3, 4]),
rfem.structure_core.Line(no=4, definition_nodes=[4, 5]),
rfem.structure_core.Line(no=5, definition_nodes=[5, 1]),
# Surface
rfem.structure_core.Surface(
no=1,
thickness=1,
material=1,
boundary_lines=[1, 2, 3, 4, 5]
),
# Nodal Supports
rfem.types_for_nodes.NodalSupport(
no=1,
user_defined_name_enabled=True,
name="Hinged",
nodes=[3, 4, 5],
spring=common.Vector3d(x=inf, y=inf, z=inf),
rotational_restraint=common.Vector3d(x=0, y=0, z=inf),
),
# Line Supports
rfem.types_for_lines.LineSupport(
no=1,
user_defined_name_enabled=True,
name="Hinged",
lines=[1],
spring=common.Vector3d(x=inf, y=inf, z=inf),
rotational_restraint=common.Vector3d(x=0, y=0, z=0),
),
]
def define_loading() -> list:
return [
# Static Analysis Settings
rfem.loading.StaticAnalysisSettings(
no=1,
),
# Load Cases
rfem.loading.LoadCase(
no=1,
name="Self-weight",
static_analysis_settings=1,
self_weight_active=False,
),
# Uniform surface load of 250 kN/m² in global Z direction
rfem.loads.SurfaceLoad(
no=1,
surfaces=[1],
load_case=1,
uniform_magnitude=250000,
load_type=rfem.loads.SurfaceLoad.LoadType.LOAD_TYPE_FORCE,
load_distribution=rfem.loads.SurfaceLoad.LoadDistribution.LOAD_DISTRIBUTION_UNIFORM,
load_direction=rfem.loads.SurfaceLoad.LoadDirection.LOAD_DIRECTION_GLOBAL_Z_OR_USER_DEFINED_W_TRUE_LENGTH
)
]
def _get_col(df, *names):
lower = {col.lower(): col for col in df.columns}
for name in names:
if name.lower() in lower:
return lower[name.lower()]
raise KeyError(f"None of columns {names} found in: {list(df.columns)}")
def _extract_node_cols(df):
node_cols = []
for col in df.columns:
match = re.search(r"node\D*(\d+)", col.lower())
if match:
node_cols.append((int(match.group(1)), col))
node_cols.sort(key=lambda item: item[0])
return [col for _, col in node_cols]
def _parse_node_coords(df) -> dict:
id_col = _get_col(df, "id", "no")
x_col = _get_col(df, "x")
y_col = _get_col(df, "y")
z_col = _get_col(df, "z")
coords = {}
for _, row in df.iterrows():
try:
node_id = int(row[id_col])
except (TypeError, ValueError):
continue
coords[node_id] = (float(row[x_col]), float(row[y_col]), float(row[z_col]))
return coords
def plot_deformed_mesh(nodes_deformed_df, nodes_original_df, elements_2d_df,
nodal_support_nodes=None, line_support_node_pairs=None,
deformation_scale=0.9):
node_positions = _parse_node_coords(nodes_deformed_df)
node_positions_org = _parse_node_coords(nodes_original_df)
# Actual deflections — used for labels and print output
z_deflections = {
node_id: node_positions[node_id][2] - node_positions_org[node_id][2]
for node_id in node_positions
if node_id in node_positions_org
}
# Visually scaled positions — deformation_scale < 1 reduces, > 1 exaggerates
node_positions_plot = {
node_id: (
x_org + (node_positions[node_id][0] - x_org) * deformation_scale,
y_org + (node_positions[node_id][1] - y_org) * deformation_scale,
z_org + (node_positions[node_id][2] - z_org) * deformation_scale,
)
for node_id, (x_org, y_org, z_org) in node_positions_org.items()
if node_id in node_positions
}
max_deflection_node = max(z_deflections, key=lambda nid: abs(z_deflections[nid])) if z_deflections else None
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
xs = [c[0] for c in node_positions_plot.values()]
ys = [c[1] for c in node_positions_plot.values()]
zs = [c[2] for c in node_positions_plot.values()]
ax.scatter(xs, ys, zs, color="black", marker="o", s=0.1)
node_cols_2d = _extract_node_cols(elements_2d_df)
for _, row in elements_2d_df.iterrows():
poly_ids = [node_positions_plot[int(row[col])] for col in node_cols_2d
if row[col] == row[col] and row[col] is not None # skip NaN
and int(row[col]) in node_positions_plot]
if len(poly_ids) >= 3:
ax.add_collection3d(
Poly3DCollection(
[poly_ids],
facecolors="lightgray",
linewidths=0.2,
edgecolors="gray",
alpha=0.5,
)
)
# Draw nodal supports as triangles at original node positions
if nodal_support_nodes:
for nid in nodal_support_nodes:
if nid in node_positions_org:
x, y, z = node_positions_org[nid]
ax.scatter([x], [y], [z], marker="^", color="black", s=40, zorder=50, depthshade=False)
# Draw line supports as a row of triangles along each supported edge
if line_support_node_pairs:
for n1, n2 in line_support_node_pairs:
if n1 in node_positions_org and n2 in node_positions_org:
x1, y1, z1 = node_positions_org[n1]
x2, y2, z2 = node_positions_org[n2]
for t in [i / 6 for i in range(7)]:
ax.scatter(
[x1 + t * (x2 - x1)],
[y1 + t * (y2 - y1)],
[z1 + t * (z2 - z1)],
marker="^", color="black", s=40, zorder=50, depthshade=False,
)
if max_deflection_node is not None:
x_max, y_max, z_max = node_positions_plot[max_deflection_node]
z_def_mm = abs(z_deflections[max_deflection_node]) * 1000.0
print(f"Max deflection at node {max_deflection_node}: {z_def_mm:.1f} mm")
ax.scatter([x_max], [y_max], [z_max], color="red", marker="o", s=30,
depthshade=False, zorder=100)
ax.text(x_max, y_max, z_max, f" {z_def_mm:.1f} mm", color="red", fontweight="bold")
# Z axis: -1x max deflection above slab (margin), +2x max deflection below (e.g. -77mm to +154mm)
max_def = max(abs(d) for d in z_deflections.values()) if z_deflections else 0.1
ax.set_zlim(-max_def, 2 * max_def)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z ↓")
ax.set_title("RFEM Deformed Mesh")
ax.set_box_aspect((1, 1, 1))
ax.grid(False)
ax.invert_zaxis()
ax.view_init(elev=25, azim=-60)
plt.show()
# -------------------------------------------------------
# MAIN SCRIPT
# -------------------------------------------------------
with rfem.Application() as rfem_app:
# Create model
rfem_app.close_all_models(save_changes=False)
rfem_app.create_model(name="mesh_table.rf6")
rfem_app.delete_all_objects()
rfem_app.create_object_list(
objs=define_structure() + define_loading()
)
# Update mesh settings
mesh_settings = rfem_app.get_mesh_settings()
mesh_settings.surfaces_shape_of_finite_elements = rfem.mesh.MeshSettings.SURFACES_SHAPE_OF_FINITE_ELEMENTS_FOR_SURFACES__QUADRANGLES_ONLY
mesh_settings.general_target_length_of_fe = 0.1
rfem_app.set_mesh_settings(mesh_settings=mesh_settings)
rfem_app.calculate_all(skip_warnings=True)
# --- Retrieve mesh data ---
# Deformed node positions (for LC1)
mesh_nodes_deformed = rfem_app.get_mesh_table(
mesh_entity=rfem.mesh.MeshEntity.NODE,
mesh_shape=rfem.mesh.MeshShape.DEFORMED,
loading=rfem.ObjectId(
no=1,
object_type=rfem.OBJECT_TYPE_LOAD_CASE,
),
load_increment=1
).data
print(f"Mesh Nodes - Deformed Shape:\n{mesh_nodes_deformed}")
# Original (undeformed) node positions
mesh_nodes_original = rfem_app.get_mesh_table(
mesh_entity=rfem.mesh.MeshEntity.NODE,
mesh_shape=rfem.mesh.MeshShape.ORIGINAL,
).data
# 2D mesh elements
mesh_2d_elements = rfem_app.get_mesh_table(
mesh_entity=rfem.mesh.MeshEntity.ELEMENT_2D,
mesh_shape=rfem.mesh.MeshShape.ORIGINAL,
).data
# --- Visualise deformed mesh ---
plot_deformed_mesh(
nodes_deformed_df=mesh_nodes_deformed,
nodes_original_df=mesh_nodes_original,
elements_2d_df=mesh_2d_elements,
nodal_support_nodes=[3, 4, 5], # hinged supports at corners
line_support_node_pairs=[(1, 2)], # hinged support along line 1
)