Skip to content

Commit

Permalink
cleanup^
Browse files Browse the repository at this point in the history
  • Loading branch information
piotrbartman committed Jul 23, 2024
1 parent 14b24e1 commit 21fd5f4
Show file tree
Hide file tree
Showing 6 changed files with 71 additions and 96 deletions.
129 changes: 50 additions & 79 deletions conmech/plotting/membrane.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
from matplotlib import tri

from conmech.state.products.intersection_contact_limit_points import (
IntersectionContactLimitPoints,
VerticalIntersectionContactLimitPoints,
)
from conmech.state.state import State

Expand Down Expand Up @@ -37,7 +37,7 @@ def plot_in_columns(states: List[State], *args, **kwargs):


def plot_limit_points(
prod: IntersectionContactLimitPoints,
prod: VerticalIntersectionContactLimitPoints,
color="black",
title=None,
label=None,
Expand Down Expand Up @@ -73,6 +73,7 @@ def plot(state: State, *args, **kwargs):
plt.show()


# pylint: disable=too-many-arguments, too-many-locals
def do_plot(
fig,
state: State,
Expand All @@ -89,43 +90,51 @@ def do_plot(
x=0.0,
):
assert state.body.mesh.dimension == 2 # membrane have to be 2D
X = state.body.mesh.nodes[:, 0]
Y = state.body.mesh.nodes[:, 1]
mesh_node_x = state.body.mesh.nodes[:, 0]
mesh_node_y = state.body.mesh.nodes[:, 1]

soltri = tri.Triangulation(X, Y, triangles=state.body.mesh.elements)
soltri = tri.Triangulation(
mesh_node_x, mesh_node_y, triangles=state.body.mesh.elements
)
interpol = tri.LinearTriInterpolator # if in3d else tri.CubicTriInterpolator
v = interpol(soltri, getattr(state, field)[:, 0])
u = interpol(soltri, state.displacement[:, 0])

# higher resolution
X = np.linspace(0, 1, 100)
Y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(X, Y)
X = X.ravel()
Y = Y.ravel()
U = np.zeros_like(X)
for i in range(len(U)):
U[i] = u(X[i], Y[i])

B = np.linspace(0, 1, 100)
RB = np.linspace(1, 0, 100)
BX = np.concatenate((B, np.ones_like(B), RB, np.zeros_like(RB), np.zeros(1)))
BY = np.concatenate((np.zeros_like(B), B, np.ones_like(RB), RB, np.zeros(1)))
BU = np.empty_like(BX)
for i in range(len(BX)):
BU[i] = u(BX[i], BY[i])
x_disc = np.linspace(0, 1, 100)
y_disc = np.linspace(0, 1, 100)
node_x, node_y = np.meshgrid(x_disc, y_disc)
node_x = node_x.ravel()
node_y = node_y.ravel()
node_val = np.zeros_like(node_x)
# pylint: disable=consider-using-enumerate
for i in range(len(node_val)):
node_val[i] = u(node_x[i], node_y[i])

bound = np.linspace(0, 1, 100)
rev_bound = np.linspace(1, 0, 100)
bound_x = np.concatenate(
(bound, np.ones_like(bound), rev_bound, np.zeros_like(rev_bound), np.zeros(1))
)
bound_y = np.concatenate(
(np.zeros_like(bound), bound, np.ones_like(rev_bound), rev_bound, np.zeros(1))
)
bound_val = np.empty_like(bound_x)
# pylint: disable=consider-using-enumerate
for i in range(len(bound_x)):
bound_val[i] = u(bound_x[i], bound_y[i])

if in3d:
ax.view_init(elev=elev, azim=azim)
plot_surface(
fig,
ax,
X,
Y,
U,
BX,
BY,
BU,
node_x,
node_y,
node_val,
bound_x,
bound_y,
bound_val,
lambda x, y, z: v(x, y),
vmin,
vmax,
Expand All @@ -134,49 +143,47 @@ def do_plot(
title,
)
else:
plot_intersection(fig, ax, u, x, ymin=zmin, ymax=zmax)
plot_intersection(ax, u, x, ymin=zmin, ymax=zmax)


def plot_surface(
fig,
ax,
X,
Y,
U,
BX,
BY,
BU,
v: Callable,
node_x,
node_y,
node_val,
bound_x,
bound_y,
bound_val,
v_func: Callable,
vmin=None,
vmax=None,
zmin=None,
zmax=None,
title=None,
):
p3dc = ax.plot_trisurf(X, Y, U, alpha=1)
p3dc = ax.plot_trisurf(node_x, node_y, node_val, alpha=1)
ax.set_zlim(zmin=zmin, zmax=zmax)

mappable = map_colors(p3dc, v, "coolwarm", vmin, vmax)
mappable = map_colors(p3dc, v_func, "coolwarm", vmin, vmax)
plt.title(title)

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("u")

ax.plot3D(BX, BY, BU, color="blue")
ax.plot3D(bound_x, bound_y, bound_val, color="blue")

cbar_ax = fig.add_axes([0.10, 0.957, 0.8, 0.02])
plt.colorbar(mappable, shrink=0.6, aspect=15, cax=cbar_ax, orientation="horizontal")

# plt.show()

# v2_mayavi(True, X, Y, U, np.ones_like(U))


def plot_intersection(fig, ax, u, x, ymin, ymax):
def plot_intersection(ax, valfunc, insec_x, ymin, ymax):
space = np.linspace(0, 1, 100)
ax.set_ylim([ymin, ymax])
ax.plot(space, u(np.ones_like(space) * x, space))
ax.plot(space, valfunc(np.ones_like(space) * insec_x, space))


def map_colors(p3dc, func, cmap="viridis", vmin=None, vmax=None):
Expand All @@ -191,6 +198,7 @@ def map_colors(p3dc, func, cmap="viridis", vmin=None, vmax=None):
"""

# reconstruct the triangles from internal data
# pylint: disable=protected-access
x, y, z, _ = p3dc._vec
slices = p3dc._segslices
triangles = np.array([np.array((x[s], y[s], z[s])).T for s in slices])
Expand All @@ -210,40 +218,3 @@ def map_colors(p3dc, func, cmap="viridis", vmin=None, vmax=None):

# if the caller wants a colorbar, they need this
return ScalarMappable(cmap=cmap, norm=norm)


def v2_mayavi(transparency, X, Y, Z, O):
X = X.reshape(100, 100)
Y = Y.reshape(100, 100)
Z = Z.reshape(100, 100)
O = O.reshape(100, 100)
from mayavi import mlab

# mlab.test_contour3d()
# mlab.show()
# fig = mlab.figure()

# ax_ranges = [-2, 2, -2, 2, 0, 8]
# ax_scale = [1.0, 1.0, 0.4]
# ax_extent = ax_ranges * np.repeat(ax_scale, 2)

X = np.linspace(0, 1, 100)
Y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(X, Y, indexing="ij")
surf3 = mlab.surf(X, Y, Z)
# surf4 = mlab.surf(X, Y, O)

# surf3.actor.actor.scale = ax_scale
# surf4.actor.actor.scale = ax_scale
# mlab.view(60, 74, 17, [-2.5, -4.6, -0.3])
# mlab.outline(surf3, color=(.7, .7, .7),)# extent=ax_extent)
# mlab.axes(surf3, color=(.7, .7, .7),)# extent=ax_extent,
# ranges=ax_ranges,
# label='x', ylabel='y', zlabel='z')

# if transparency:
# surf3.actor.property.opacity = 0.5
# surf4.actor.property.opacity = 0.5
# fig.scene.renderer.use_depth_peeling = 1

mlab.show()
4 changes: 2 additions & 2 deletions conmech/simulations/problem_solver.py
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ def __set_step_solver(self, value):
statement,
self.body,
self.time_step,
getattr(self.problem, 'contact_law', None),
getattr(self.problem, "contact_law", None),
driving_vector=self.driving_vector,
)

Expand All @@ -127,7 +127,7 @@ def __set_second_step_solver(self, value):
statement,
self.body,
self.time_step,
getattr(self.problem, 'contact_law_2', None),
getattr(self.problem, "contact_law_2", None),
driving_vector=self.driving_vector,
)

Expand Down
8 changes: 4 additions & 4 deletions conmech/state/products/intersection_contact_limit_points.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,10 +24,10 @@
from conmech.state.products.product import Product


class IntersectionContactLimitPoints(Product):
class VerticalIntersectionContactLimitPoints(Product):
def __init__(self, x, obstacle_level):
super().__init__(f"limit points at {x:.2f}")
self.x = x
self.level = x
self.obstacle_level = obstacle_level

def update(self, state):
Expand All @@ -36,8 +36,8 @@ def update(self, state):
step = len(state.body.mesh.nodes) ** -1 / 2
u = interpolate(state, "displacement")

def u_intsec(y):
return u(self.x, y) - self.obstacle_level
def u_intsec(y_coord):
return u(self.level, y_coord) - self.obstacle_level

self.data[state.time] = estimate_zeros(u_intsec, y_min, y_max, step)

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,10 @@
from conmech.state.products.product import Product


class Intersection(Product):
class VerticalIntersection(Product):
def __init__(self, x):
super().__init__(f"intersection at {x:.2f}")
self.x = x
self.level = x

def update(self, state):
y_min = np.min(state.body.mesh.nodes[:, 1])
Expand All @@ -39,4 +39,4 @@ def update(self, state):
soltri = tri.Triangulation(X, Y, triangles=state.body.mesh.elements)
u = tri.LinearTriInterpolator(soltri, state.displacement[:, 0])

self.data[state.time] = (space, u(np.ones_like(space) * self.x, space))
self.data[state.time] = (space, u(np.ones_like(space) * self.level, space))
10 changes: 6 additions & 4 deletions examples/BOSK_2024_example_1.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,9 @@
from conmech.scenarios.problems import ContactWaveProblem
from conmech.simulations.problem_solver import WaveSolver
from conmech.properties.mesh_description import CrossMeshDescription
from conmech.state.products.intersection import Intersection
from conmech.state.products.verticalintersection import VerticalIntersection
from conmech.state.products.intersection_contact_limit_points import (
IntersectionContactLimitPoints,
VerticalIntersectionContactLimitPoints,
)
from conmech.state.state import State

Expand Down Expand Up @@ -72,8 +72,10 @@ def main(config: Config, setup, name, steps):
n_steps=steps,
output_step=output_step,
products=[
IntersectionContactLimitPoints(obstacle_level=OBSTACLE_LEVEL, x=1.0),
Intersection(x=1.0),
VerticalIntersectionContactLimitPoints(
obstacle_level=OBSTACLE_LEVEL, x=1.0
),
VerticalIntersection(x=1.0),
],
initial_displacement=setup.initial_displacement,
initial_velocity=setup.initial_velocity,
Expand Down
10 changes: 6 additions & 4 deletions examples/BOSK_2024_example_2.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,9 @@
from conmech.scenarios.problems import InteriorContactWaveProblem
from conmech.simulations.problem_solver import WaveSolver
from conmech.properties.mesh_description import CrossMeshDescription
from conmech.state.products.intersection import Intersection
from conmech.state.products.verticalintersection import VerticalIntersection
from conmech.state.products.intersection_contact_limit_points import (
IntersectionContactLimitPoints,
VerticalIntersectionContactLimitPoints,
)
from conmech.state.state import State

Expand Down Expand Up @@ -96,8 +96,10 @@ def main(config: Config, setup, name, steps):
n_steps=steps,
output_step=output_step,
products=[
IntersectionContactLimitPoints(obstacle_level=OBSTACLE_LEVEL, x=0.50),
Intersection(x=0.50),
VerticalIntersectionContactLimitPoints(
obstacle_level=OBSTACLE_LEVEL, x=0.50
),
VerticalIntersection(x=0.50),
],
initial_displacement=setup.initial_displacement,
initial_velocity=setup.initial_velocity,
Expand Down

0 comments on commit 21fd5f4

Please sign in to comment.