diff --git a/data/in_hypar_mesh.json b/data/in_hypar_mesh.json new file mode 100644 index 00000000..f57e8656 --- /dev/null +++ b/data/in_hypar_mesh.json @@ -0,0 +1 @@ +{"compas": "1.7.1", "datatype": "compas.datastructures/Mesh", "data": {"facedata": {"5": {"path": [1, 0]}, "6": {"path": [1, 1]}, "7": {"path": [1, 2]}, "8": {"path": [1, 3]}, "9": {"path": [2, 0]}, "10": {"path": [2, 1]}, "11": {"path": [2, 2]}, "12": {"path": [2, 3]}, "13": {"path": [3, 0]}, "14": {"path": [3, 1]}, "15": {"path": [3, 2]}, "16": {"path": [3, 3]}, "17": {"path": [4, 0]}, "18": {"path": [4, 1]}, "19": {"path": [4, 2]}, "20": {"path": [4, 3]}}, "dva": {"y": 0.0, "x": 0.0, "z": 0.0}, "face": {"5": [9, 0, 10, 21], "6": [10, 2, 13, 21], "7": [13, 3, 15, 21], "8": [15, 7, 9, 21], "9": [12, 1, 11, 22], "10": [11, 8, 16, 22], "11": [16, 3, 17, 22], "12": [17, 4, 12, 22], "13": [17, 3, 13, 23], "14": [13, 2, 14, 23], "15": [14, 5, 18, 23], "16": [18, 4, 17, 23], "17": [15, 3, 16, 24], "18": [16, 8, 20, 24], "19": [20, 6, 19, 24], "20": [19, 7, 15, 24]}, "edgedata": {}, "dfa": {}, "max_vertex": 24, "dea": {}, "attributes": {"name": "Mesh"}, "max_face": 20, "vertex": {"0": {"z": 20.0, "y": -46.983559137412122, "x": -44.487094537112114}, "1": {"z": 20.0, "y": 43.016440862587885, "x": 45.512905462887865}, "2": {"z": 10.0, "y": -46.983559137412122, "x": 0.51290546288787808}, "3": {"z": 0.0, "y": -1.9835591374121182, "x": 0.51290546288788419}, "4": {"z": 10.0, "y": -1.9835591374121213, "x": 45.512905462887865}, "5": {"z": 0.0, "y": -46.983559137412122, "x": 45.512905462887865}, "6": {"z": 0.0, "y": 43.016440862587885, "x": -44.487094537112114}, "7": {"z": 10.0, "y": -1.9835591374121186, "x": -44.487094537112114}, "8": {"z": 10.0, "y": 43.016440862587885, "x": 0.51290546288789052}, "9": {"z": 15.0, "y": -24.483559137412122, "x": -44.487094537112114}, "10": {"z": 15.0, "y": -46.983559137412122, "x": -21.987094537112117}, "11": {"z": 15.0, "y": 43.016440862587885, "x": 23.012905462887879}, "12": {"z": 15.0, "y": 20.516440862587881, "x": 45.512905462887865}, "13": {"z": 5.0, "y": -24.483559137412119, "x": 0.51290546288788108}, "14": {"z": 5.0, "y": -46.983559137412122, "x": 23.012905462887872}, "15": {"z": 5.0, "y": -1.9835591374121184, "x": -21.987094537112117}, "16": {"z": 5.0, "y": 20.516440862587881, "x": 0.5129054628878873}, "17": {"z": 5.0, "y": -1.9835591374121198, "x": 23.012905462887872}, "18": {"z": 5.0, "y": -24.483559137412122, "x": 45.512905462887865}, "19": {"z": 5.0, "y": 20.516440862587885, "x": -44.487094537112114}, "20": {"z": 5.0, "y": 43.016440862587885, "x": -21.98709453711211}, "21": {"z": 10.0, "y": -24.483559137412119, "x": -21.987094537112114}, "22": {"z": 10.0, "y": 20.516440862587885, "x": 23.012905462887879}, "23": {"z": 5.0, "y": -24.483559137412119, "x": 23.012905462887872}, "24": {"z": 5.0, "y": 20.516440862587885, "x": -21.987094537112114}}}} \ No newline at end of file diff --git a/scripts/91_nfd_hypar.py b/scripts/91_nfd_hypar.py new file mode 100644 index 00000000..7b829f95 --- /dev/null +++ b/scripts/91_nfd_hypar.py @@ -0,0 +1,68 @@ +from os import path + +from compas.datastructures import Mesh, mesh_subdivide +from compas_view2 import app + +from compas_fd.nfd import nfd_ur_numpy +from _nfd_helpers import mesh_update + + +# ================================================= +# IO +# ================================================= + +HERE = path.dirname(__file__) +FILE_I = path.join(HERE, '..', 'data', 'in_hypar_mesh.json') +mesh = Mesh.from_json(FILE_I) + + +# ================================================= +# input mesh +# ================================================= + +mesh.vertices_attribute('is_anchor', True, + mesh.vertices_where({'vertex_degree': 2})) + +mesh = mesh_subdivide(mesh, scheme='quad', k=2) +mesh.quads_to_triangles() + +bounds = mesh.edges_on_boundaries()[0] +mesh.edges_attribute('q_pre', 20, bounds) + +dva = {'rx': .0, 'ry': .0, 'rz': .0, + 'px': .0, 'py': .0, 'pz': .0, + 'is_anchor': False} +dea = {'q_pre': .0} +dfa = {'s_pre': [1.0, 1.0, 0.0]} + +mesh.update_default_vertex_attributes(dva) +mesh.update_default_edge_attributes(dea) +mesh.update_default_face_attributes(dfa) + + +# ================================================= +# get mesh data +# ================================================= + +fixed = [mesh.key_index()[v] for v in mesh.vertices_where({'is_anchor': True})] +P = mesh.vertices_attributes(['px', 'py', 'pz']) +S = mesh.faces_attribute('s_pre') +Q = mesh.edges_attribute('q_pre') + + +# ================================================= +# solver +# ================================================= + +xyz, r, f, _, s = nfd_ur_numpy(mesh, fixed, S, force_density_goals=Q, + vertex_loads=P, kmax=10, stress_flag=1) +mesh_update(mesh, xyz, r, s, f) + + +# ================================================= +# viz +# ================================================= + +viewer = app.App() +viewer.add(mesh) +viewer.show() diff --git a/scripts/92_nfd_hypar_opening.py b/scripts/92_nfd_hypar_opening.py new file mode 100644 index 00000000..2d9e103e --- /dev/null +++ b/scripts/92_nfd_hypar_opening.py @@ -0,0 +1,72 @@ +from os import path + +from compas.datastructures import Mesh, mesh_subdivide +from compas_view2 import app + +from compas_fd.nfd import nfd_ur_numpy +from _nfd_helpers import mesh_update + + +# ================================================= +# IO +# ================================================= + +HERE = path.dirname(__file__) +FILE_I = path.join(HERE, '..', 'data', 'in_hypar_mesh.json') +mesh = Mesh.from_json(FILE_I) + + +# ================================================= +# input mesh +# ================================================= + +mesh.vertices_attribute('is_anchor', True, + mesh.vertices_where({'vertex_degree': 2})) + +central_vertex = 3 +mesh.delete_vertex(3) + +mesh = mesh_subdivide(mesh, scheme='quad', k=3) +mesh.quads_to_triangles() + +in_bounds, out_bounds = sorted(mesh.edges_on_boundaries(), key=len) +mesh.edges_attribute('q_pre', 11, in_bounds) +mesh.edges_attribute('q_pre', 30, out_bounds) + +dva = {'rx': .0, 'ry': .0, 'rz': .0, + 'px': .0, 'py': .0, 'pz': .0, + 'is_anchor': False} +dea = {'q_pre': .0} +dfa = {'s_pre': [1.0, 1.0, 0.0]} + +mesh.update_default_vertex_attributes(dva) +mesh.update_default_edge_attributes(dea) +mesh.update_default_face_attributes(dfa) + + +# ================================================= +# get mesh data +# ================================================= + +fixed = [mesh.key_index()[v] for v in mesh.vertices_where({'is_anchor': True})] +P = mesh.vertices_attributes(['px', 'py', 'pz']) +S = mesh.faces_attribute('s_pre') +Q = mesh.edges_attribute('q_pre') + + +# ================================================= +# solver +# ================================================= + +xyz, r, f, _, s = nfd_ur_numpy(mesh, fixed, S, force_density_goals=Q, vertex_loads=P, + kmax=10, stress_flag=1, stress_tol=.05, xyz_tol=.01) +mesh_update(mesh, xyz, r, s, f) + + +# ================================================= +# viz +# ================================================= + +viewer = app.App() +viewer.add(mesh) +viewer.show() diff --git a/scripts/93_nfd_hypar_quads.py b/scripts/93_nfd_hypar_quads.py new file mode 100644 index 00000000..dae62323 --- /dev/null +++ b/scripts/93_nfd_hypar_quads.py @@ -0,0 +1,67 @@ +from os import path + +from compas.datastructures import Mesh, mesh_subdivide +from compas_view2 import app + +from compas_fd.nfd import nfd_ur_numpy +from _nfd_helpers import mesh_update + + +# ================================================= +# IO +# ================================================= + +HERE = path.dirname(__file__) +FILE_I = path.join(HERE, '..', 'data', 'in_hypar_mesh.json') +mesh = Mesh.from_json(FILE_I) + + +# ================================================= +# input mesh +# ================================================= + +mesh.vertices_attribute('is_anchor', True, + mesh.vertices_where({'vertex_degree': 2})) + +mesh = mesh_subdivide(mesh, scheme='quad', k=3) + +bounds = mesh.edges_on_boundaries()[0] +mesh.edges_attribute('q_pre', 30, bounds) + +dva = {'rx': .0, 'ry': .0, 'rz': .0, + 'px': .0, 'py': .0, 'pz': .0, + 'is_anchor': False} +dea = {'q_pre': .0} +dfa = {'s_pre': [1.0, 1.0, 0.0]} + +mesh.update_default_vertex_attributes(dva) +mesh.update_default_edge_attributes(dea) +mesh.update_default_face_attributes(dfa) + + +# ================================================= +# get mesh data +# ================================================= + +fixed = [mesh.key_index()[v] for v in mesh.vertices_where({'is_anchor': True})] +P = mesh.vertices_attributes(['px', 'py', 'pz']) +S = mesh.faces_attribute('s_pre') +Q = mesh.edges_attribute('q_pre') + + +# ================================================= +# solver +# ================================================= + +xyz, r, f, _, s = nfd_ur_numpy(mesh, fixed, S, force_density_goals=Q, vertex_loads=P, + kmax=10, stress_flag=1, stress_tol=.05, xyz_tol=.01) +mesh_update(mesh, xyz, r, s, f) + + +# ================================================= +# viz +# ================================================= + +viewer = app.App() +viewer.add(mesh) +viewer.show() diff --git a/scripts/94_nfd_hypar_quads_tris.py b/scripts/94_nfd_hypar_quads_tris.py new file mode 100644 index 00000000..094e80e2 --- /dev/null +++ b/scripts/94_nfd_hypar_quads_tris.py @@ -0,0 +1,71 @@ +from os import path + +from compas.datastructures import Mesh, mesh_subdivide +from compas_view2 import app + +from compas_fd.nfd import nfd_ur_numpy +from _nfd_helpers import mesh_update + + +# ================================================= +# IO +# ================================================= + +HERE = path.dirname(__file__) +FILE_I = path.join(HERE, '..', 'data', 'in_hypar_mesh.json') +mesh = Mesh.from_json(FILE_I) + + +# ================================================= +# input mesh +# ================================================= + +mesh.vertices_attribute('is_anchor', True, + mesh.vertices_where({'vertex_degree': 2})) + +mesh = mesh_subdivide(mesh, scheme='quad', k=2) +faces_to_split = (mesh.get_any_face() for _ in range(200)) +for face in faces_to_split: + v = mesh.face_vertices(face) + mesh.split_face(face, v[0], v[2]) + +bounds = mesh.edges_on_boundaries()[0] +mesh.edges_attribute('q_pre', 15, bounds) + +dva = {'rx': .0, 'ry': .0, 'rz': .0, + 'px': .0, 'py': .0, 'pz': .0, + 'is_anchor': False} +dea = {'q_pre': .0} +dfa = {'s_pre': [1.0, 1.0, 0.0]} + +mesh.update_default_vertex_attributes(dva) +mesh.update_default_edge_attributes(dea) +mesh.update_default_face_attributes(dfa) + + +# ================================================= +# get mesh data +# ================================================= + +fixed = [mesh.key_index()[v] for v in mesh.vertices_where({'is_anchor': True})] +P = mesh.vertices_attributes(['px', 'py', 'pz']) +S = mesh.faces_attribute('s_pre') +Q = mesh.edges_attribute('q_pre') + + +# ================================================= +# solver +# ================================================= + +xyz, r, f, _, s = nfd_ur_numpy(mesh, fixed, S, force_density_goals=Q, vertex_loads=P, + kmax=10, stress_flag=1, stress_tol=.05, xyz_tol=.01) +mesh_update(mesh, xyz, r, s, f) + + +# ================================================= +# viz +# ================================================= + +viewer = app.App() +viewer.add(mesh) +viewer.show() diff --git a/scripts/95_nfd_hypar_support.py b/scripts/95_nfd_hypar_support.py new file mode 100644 index 00000000..d5b61920 --- /dev/null +++ b/scripts/95_nfd_hypar_support.py @@ -0,0 +1,95 @@ +from os import path + +from compas.datastructures import Mesh, mesh_subdivide +from compas_view2 import app + +from compas_fd.nfd import nfd_ur_numpy +from _nfd_helpers import mesh_update + + +# ================================================= +# IO +# ================================================= + +HERE = path.dirname(__file__) +FILE_I = path.join(HERE, '..', 'data', 'in_hypar_mesh.json') +mesh = Mesh.from_json(FILE_I) + + +# ================================================= +# input mesh +# ================================================= + +mesh.vertices_attribute('is_anchor', True, + mesh.vertices_where({'vertex_degree': 2})) + +mesh = mesh_subdivide(mesh, scheme='quad', k=3) +mesh.quads_to_triangles() + +# get ring of vertices around central vertex +central_vertex = 3 +ring_size = 3 +ring_all = mesh.vertex_neighborhood(central_vertex, ring=ring_size) +ring_in = mesh.vertex_neighborhood(central_vertex, ring=ring_size - 1) +ring = set(ring_all) - set(ring_in) +mesh.vertices_attributes(['z', 'is_anchor'], [25, True], ring) + +# triangulate top faces at supports +fva = mesh.face_vertex_after +for v in ring: + faces = mesh.vertex_faces(v) + + for face in faces: + vts = mesh.face_vertices(face) + v_count = len(vts) + if v_count != 4: + continue + w = fva(face, fva(face, v)) + if not mesh.vertex_attribute(w, 'is_anchor'): + continue + mesh.split_face(face, v, w) + +# delete vertices inside support ring +for v in ring_in: + mesh.delete_vertex(v) + +bounds = mesh.edges_on_boundaries()[0] +mesh.edges_attribute('q_pre', 30, bounds) + +dva = {'rx': .0, 'ry': .0, 'rz': .0, + 'px': .0, 'py': .0, 'pz': .0, + 'is_anchor': False} +dea = {'q_pre': .0} +dfa = {'s_pre': [1.0, 1.0, 0.0]} + +mesh.update_default_vertex_attributes(dva) +mesh.update_default_edge_attributes(dea) +mesh.update_default_face_attributes(dfa) + + +# ================================================= +# get mesh data +# ================================================= + +fixed = [mesh.key_index()[v] for v in mesh.vertices_where({'is_anchor': True})] +P = mesh.vertices_attributes(['px', 'py', 'pz']) +S = mesh.faces_attribute('s_pre') +Q = mesh.edges_attribute('q_pre') + + +# ================================================= +# solver +# ================================================= + +xyz, r, f, _, s = nfd_ur_numpy(mesh, fixed, S, force_density_goals=Q, vertex_loads=P, + kmax=10, stress_flag=1, stress_tol=.01, xyz_tol=.01) +mesh_update(mesh, xyz, r, s, f) + + +# ================================================= +# viz +# ================================================= + +viewer = app.App() +viewer.add(mesh) +viewer.show() diff --git a/scripts/96_nfd_hypar_anisotropic_stress_A.py b/scripts/96_nfd_hypar_anisotropic_stress_A.py new file mode 100644 index 00000000..b898ec98 --- /dev/null +++ b/scripts/96_nfd_hypar_anisotropic_stress_A.py @@ -0,0 +1,68 @@ +from os import path + +from compas.datastructures import Mesh, mesh_subdivide +from compas_view2 import app + +from compas_fd.nfd import nfd_ur_numpy +from _nfd_helpers import mesh_update + + +# ================================================= +# IO +# ================================================= + +HERE = path.dirname(__file__) +FILE_I = path.join(HERE, '..', 'data', 'in_hypar_mesh.json') +mesh = Mesh.from_json(FILE_I) + + +# ================================================= +# input mesh +# ================================================= + +mesh.vertices_attribute('is_anchor', True, + mesh.vertices_where({'vertex_degree': 2})) + +mesh = mesh_subdivide(mesh, scheme='quad', k=2) +mesh.quads_to_triangles() + +bounds = mesh.edges_on_boundaries()[0] +mesh.edges_attribute('q_pre', 30, bounds) + +dva = {'rx': .0, 'ry': .0, 'rz': .0, + 'px': .0, 'py': .0, 'pz': .0, + 'is_anchor': False} +dea = {'q_pre': .0} +dfa = {'s_pre': [2.5, 1.0, 0.0]} + +mesh.update_default_vertex_attributes(dva) +mesh.update_default_edge_attributes(dea) +mesh.update_default_face_attributes(dfa) + + +# ================================================= +# get mesh data +# ================================================= + +fixed = [mesh.key_index()[v] for v in mesh.vertices_where({'is_anchor': True})] +P = mesh.vertices_attributes(['px', 'py', 'pz']) +S = mesh.faces_attribute('s_pre') +Q = mesh.edges_attribute('q_pre') + + +# ================================================= +# solver +# ================================================= + +xyz, r, f, _, s = nfd_ur_numpy(mesh, fixed, S, force_density_goals=Q, vertex_loads=P, + kmax=10, stress_flag=3, stress_ref=(1, 1, 0)) +mesh_update(mesh, xyz, r, s, f) + + +# ================================================= +# viz +# ================================================= + +viewer = app.App() +viewer.add(mesh) +viewer.show() diff --git a/scripts/97_nfd_hypar_anisotropic_stress_B.py b/scripts/97_nfd_hypar_anisotropic_stress_B.py new file mode 100644 index 00000000..d9658608 --- /dev/null +++ b/scripts/97_nfd_hypar_anisotropic_stress_B.py @@ -0,0 +1,68 @@ +from os import path + +from compas.datastructures import Mesh, mesh_subdivide +from compas_view2 import app + +from compas_fd.nfd import nfd_ur_numpy +from _nfd_helpers import mesh_update + + +# ================================================= +# IO +# ================================================= + +HERE = path.dirname(__file__) +FILE_I = path.join(HERE, '..', 'data', 'in_hypar_mesh.json') +mesh = Mesh.from_json(FILE_I) + + +# ================================================= +# input mesh +# ================================================= + +mesh.vertices_attribute('is_anchor', True, + mesh.vertices_where({'vertex_degree': 2})) + +mesh = mesh_subdivide(mesh, scheme='quad', k=2) +mesh.quads_to_triangles() + +bounds = mesh.edges_on_boundaries()[0] +mesh.edges_attribute('q_pre', 30, bounds) + +dva = {'rx': .0, 'ry': .0, 'rz': .0, + 'px': .0, 'py': .0, 'pz': .0, + 'is_anchor': False} +dea = {'q_pre': .0} +dfa = {'s_pre': [1.0, 2.5, 0.0]} + +mesh.update_default_vertex_attributes(dva) +mesh.update_default_edge_attributes(dea) +mesh.update_default_face_attributes(dfa) + + +# ================================================= +# get mesh data +# ================================================= + +fixed = [mesh.key_index()[v] for v in mesh.vertices_where({'is_anchor': True})] +P = mesh.vertices_attributes(['px', 'py', 'pz']) +S = mesh.faces_attribute('s_pre') +Q = mesh.edges_attribute('q_pre') + + +# ================================================= +# solver +# ================================================= + +xyz, r, f, _, s = nfd_ur_numpy(mesh, fixed, S, force_density_goals=Q, vertex_loads=P, + kmax=10, stress_flag=3, stress_ref=(1, 0.5, 0)) +mesh_update(mesh, xyz, r, s, f) + + +# ================================================= +# viz +# ================================================= + +viewer = app.App() +viewer.add(mesh) +viewer.show() diff --git a/scripts/98_nfd_hypar_shell.py b/scripts/98_nfd_hypar_shell.py new file mode 100644 index 00000000..2ff11e61 --- /dev/null +++ b/scripts/98_nfd_hypar_shell.py @@ -0,0 +1,68 @@ +from os import path + +from compas.datastructures import Mesh, mesh_subdivide +from compas_view2 import app + +from compas_fd.nfd import nfd_ur_numpy +from _nfd_helpers import mesh_update + + +# ================================================= +# IO +# ================================================= + +HERE = path.dirname(__file__) +FILE_I = path.join(HERE, '..', 'data', 'in_hypar_mesh.json') +mesh = Mesh.from_json(FILE_I) + + +# ================================================= +# input mesh +# ================================================= + +mesh.vertices_attribute('is_anchor', True, + mesh.vertices_where({'vertex_degree': 2})) +mesh.vertices_attribute('z', 0, mesh.vertices_where({'vertex_degree': 2})) + +mesh = mesh_subdivide(mesh, scheme='quad', k=2) +# mesh.quads_to_triangles() + +bounds = mesh.edges_on_boundaries()[0] +mesh.edges_attribute('q_pre', -15, bounds) + +dva = {'rx': .0, 'ry': .0, 'rz': .0, + 'is_anchor': False} +dea = {'q_pre': .0} +dfa = {'s_pre': [-1.0, -1.0, 0.0], + 'px': .0, 'py': .0, 'pz': -.035} + +mesh.update_default_vertex_attributes(dva) +mesh.update_default_edge_attributes(dea) +mesh.update_default_face_attributes(dfa) + + +# ================================================= +# get mesh data +# ================================================= + +fixed = [mesh.key_index()[v] for v in mesh.vertices_where({'is_anchor': True})] +P = mesh.faces_attributes(['px', 'py', 'pz']) +S = mesh.faces_attribute('s_pre') +Q = mesh.edges_attribute('q_pre') + + +# ================================================= +# solver +# ================================================= + +xyz, r, f, _, s = nfd_ur_numpy(mesh, fixed, S, force_density_goals=Q, global_face_loads=P, kmax=5, stress_flag=1) +mesh_update(mesh, xyz, r, s, f) + + +# ================================================= +# viz +# ================================================= + +viewer = app.App() +viewer.add(mesh) +viewer.show() diff --git a/scripts/99_nfd_inflated_shell.py b/scripts/99_nfd_inflated_shell.py new file mode 100644 index 00000000..2e45fe1a --- /dev/null +++ b/scripts/99_nfd_inflated_shell.py @@ -0,0 +1,70 @@ +from os import path + +from compas.datastructures import Mesh, mesh_subdivide +from compas_view2 import app + +from compas_fd.nfd import nfd_ur_numpy +from _nfd_helpers import mesh_update + + +# ================================================= +# IO +# ================================================= + +HERE = path.dirname(__file__) +FILE_I = path.join(HERE, '..', 'data', 'in_hypar_mesh.json') +mesh = Mesh.from_json(FILE_I) + + +# ================================================= +# input mesh +# ================================================= + +mesh.vertices_attribute('is_anchor', True, + mesh.vertices_where({'vertex_degree': 2})) +mesh.vertices_attribute('z', 0, mesh.vertices_where({'vertex_degree': 2})) + +mesh = mesh_subdivide(mesh, scheme='quad', k=1) +# mesh.quads_to_triangles() + +bounds = mesh.edges_on_boundaries()[0] +mesh.edges_attribute('q_pre', -10, bounds) + +dva = {'rx': .0, 'ry': .0, 'rz': .0, + 'is_anchor': False} +dea = {'q_pre': .0} +dfa = {'s_pre': [-1.0, -1.0, 0.0], + 'px': .0, 'py': .0, 'pz': -.05} + +mesh.update_default_vertex_attributes(dva) +mesh.update_default_edge_attributes(dea) +mesh.update_default_face_attributes(dfa) + + +# ================================================= +# get mesh data +# ================================================= + +fixed = [mesh.key_index()[v] for v in mesh.vertices_where({'is_anchor': True})] +P1 = mesh.faces_attributes(['px', 'py', 'pz']) +P2 = [-0.2 * p for load_vec in P1 for p in load_vec] +S = mesh.faces_attribute('s_pre') +Q = mesh.edges_attribute('q_pre') + + +# ================================================= +# solver +# ================================================= + +xyz, r, f, _, s = nfd_ur_numpy(mesh, fixed, S, force_density_goals=Q, local_face_loads=P1, + global_face_loads=P2, kmax=5, stress_flag=1) +mesh_update(mesh, xyz, r, s, f) + + +# ================================================= +# viz +# ================================================= + +viewer = app.App() +viewer.add(mesh) +viewer.show() diff --git a/scripts/_nfd_helpers.py b/scripts/_nfd_helpers.py new file mode 100644 index 00000000..145751e2 --- /dev/null +++ b/scripts/_nfd_helpers.py @@ -0,0 +1,47 @@ + +__all__ = [ + 'mesh_update', + 'mesh_update_xyz' +] + + +def mesh_update(mesh, xyz, residuals, stresses, forces, + res_name='r', stress_name='s', force_name='f'): + """Update mesh with new coordinates, stresses, forces.""" + mesh_update_xyz(mesh, xyz) + mesh_update_residuals(mesh, residuals, res_name) + mesh_update_stresses(mesh, stresses, stress_name) + mesh_update_forces(mesh, forces, force_name) + + +def mesh_update_xyz(mesh, xyz): + """Update mesh with new vertex coordinates.""" + for vertex, coo in zip(mesh.vertices(), xyz): + mesh.vertex_attributes(vertex, 'xyz', coo) + + +def mesh_update_residuals(mesh, residuals, name='r'): + """Update mesh with new vertex coordinates.""" + for vertex, res in zip(mesh.vertices(), residuals): + mesh.vertex_attribute(vertex, name, res) + + +def mesh_update_stresses(mesh, stress_data, name='s'): + """Update mesh with face stresses.""" + stresses, vecs = stress_data + if stresses is None: + return + for face, stress in zip(mesh.faces(), stresses): + mesh.face_attribute(face, name, stress) + if vecs is None: + return + for face, vec in zip(mesh.faces(), vecs): + mesh.face_attribute(face, name + '_vecs', vec) + + +def mesh_update_forces(mesh, forces, name='f'): + """Update mesh with new edge forces.""" + if forces is None: + return + for edge, force in zip(mesh.edges(), forces): + mesh.edge_attribute(edge, name, force) diff --git a/src/compas_fd/__init__.py b/src/compas_fd/__init__.py index 66a2aa6b..7171402b 100644 --- a/src/compas_fd/__init__.py +++ b/src/compas_fd/__init__.py @@ -11,6 +11,7 @@ compas_fd.datastructures compas_fd.fd + compas_fd.nfd compas_fd.loads """ diff --git a/src/compas_fd/fd/result.py b/src/compas_fd/fd/result.py index a2058918..94e08055 100644 --- a/src/compas_fd/fd/result.py +++ b/src/compas_fd/fd/result.py @@ -1,6 +1,8 @@ from typing import Any from typing import List +from typing import Tuple from typing import NamedTuple +from typing import Union from nptyping import NDArray import numpy as np @@ -10,4 +12,5 @@ class Result(NamedTuple): vertices: NDArray[(Any, 3), np.float64] residuals: NDArray[(Any, 3), np.float64] forces: List[float] - lenghts: List[float] + lengths: List[float] + stresses: Union[Tuple, None] = None diff --git a/src/compas_fd/nfd/__init__.py b/src/compas_fd/nfd/__init__.py new file mode 100644 index 00000000..0932ef4e --- /dev/null +++ b/src/compas_fd/nfd/__init__.py @@ -0,0 +1,35 @@ +""" +****************** +compas_fd.nfd +****************** + +.. currentmodule:: compas_fd.nfd + + +Functions +========= + +.. autosummary:: + :toctree: generated/ + :nosignatures: + + nfd_ur_numpy + nfd_numpy + +""" +from __future__ import print_function +from __future__ import absolute_import +from __future__ import division + +import compas + +if not compas.IPY: + from .solvers import nfd_ur_numpy, nfd_numpy + +__all__ = [] + +if not compas.IPY: + __all__ += [ + 'nfd_ur_numpy', + 'nfd_numpy' + ] diff --git a/src/compas_fd/nfd/geometry.py b/src/compas_fd/nfd/geometry.py new file mode 100644 index 00000000..83757b61 --- /dev/null +++ b/src/compas_fd/nfd/geometry.py @@ -0,0 +1,588 @@ +from __future__ import annotations + +from typing import Any, Callable, Generator, List, Sequence, Tuple, NamedTuple +from typing_extensions import Annotated +from nptyping import NDArray + +from collections import namedtuple +from operator import itemgetter +from math import sin, cos, atan2, pi +from functools import wraps + +from numpy import append, argsort, asarray, cross, float64 +from numpy.linalg import eig, inv + +from compas.datastructures import Mesh +from compas.geometry import Frame, centroid_points, Rotation +from compas.geometry import Vector, angle_vectors_signed +from compas.utilities import pairwise + +from .math_utilities import euclidean_distance, arc_cos, is_isotropic, \ + transform_stress_angle, stress_vec_to_tensor, planar_rotation, transform_stress +from .goals import Goals + + +Vec = Annotated[Sequence[float], 3] +NpVec = NDArray[(3,), float64] +VertexArray = NDArray[(Any, 3), float64] + + +__all__ = [ + 'NaturalFace', + 'TriFace', + 'QuadFace', + 'NaturalEdge', + 'mesh_preprocess' +] + + +# ================================================= +# aliases +# ================================================= + +s, c = sin, cos + + +def s2(x): + return s(x) * s(x) + + +def c2(x): + return c(x) * c(x) + + +# ================================================= +# decorators +# ================================================= + +def cache_iter(func: Callable) -> Callable: + """Cache a method return value of the current iteration in the cache dictionary + of its class instance. After being cached, subsequent calls to the method will + return the stored result, until the cache dictionary is explicitely cleared.""" + cached_func = '_c_' + func.__name__ + + @wraps(func) + def wrapper(self, *args, **kwargs): + res = self.cache.get(cached_func, None) + if res is None: + res = func(self, *args, **kwargs) + self.cache[cached_func] = res + return res + return wrapper + + +def check_singularity(func: Callable) -> Callable: + """Verify if a numerical error occurred in the method call due to zero division.""" + def wrapper(*args, **kwargs): + try: + return func(*args, **kwargs) + except ZeroDivisionError: + raise Exception('Singularity in transformation matrix.') + return wrapper + + +# ================================================= +# natural geometry +# ================================================= + +class NaturalEdge: + """Represents geometry and force data of a single edge. + + Parameters + ---------- + mesh_xyz : ndarray of float + XYZ vertex coordinates of the whole mesh as a (n, 3) array. + vertices_ids : sequence of int + Sequence of contiguous vertex indices. + force_goal : float + Goal force for the edge. + force_density_goal : float + Goal force density for the edge, overwrites force goals. + + Attributes + ---------- + edge_count : int + The number of instantiated NaturalEdges. + xyz : ndarray of float + XYZ vertex coordinates of the edge as a (2, 3) array. + length: float + Length of edge. + force_density: float + Force density in edge for current geometry. + force : float + Force in edge for current geometry. + """ + + edge_count = 0 + + def __new__(cls, mesh_xyz, vertices_ids, force_goal, force_density_goal): + cls.edge_count += 1 + if not (force_goal or force_density_goal): + return + return super().__new__(cls) + + def __init__(self, mesh_xyz: VertexArray, vertices_ids: Sequence[int], + force_goal: float = None, force_density_goal: float = None) -> None: + self.cache = {} + self.edge_id = __class__.edge_count - 1 + self.vertices_ids = tuple(vertices_ids) + self.xyz = mesh_xyz + self.force_goal = force_goal + self.force_density_goal = force_density_goal + + def update_xyz(self, mesh_xyz: VertexArray) -> None: + """Update XYZ edge coordinates and clear geometry cache.""" + self.cache.clear() + self.xyz = mesh_xyz + + @property + @cache_iter + def length(self) -> float: + """Edge length from euclidean distance of vertces.""" + return euclidean_distance(*self.xyz) + + @property + def xyz(self) -> NDArray[(2, 3), float64]: + """XYZ vertex coordinates of the edge.""" + return self._xyz + + @xyz.setter + def xyz(self, mesh_xyz: VertexArray) -> None: + u, v = self.vertices_ids + self._xyz = (mesh_xyz[u], mesh_xyz[v]) + + @property + def force_density(self) -> float: + """Force density in edge for current geometry.""" + self._force_density = self.force_density_goal or (self.force_goal / self.length) + return self._force_density + + @property + @cache_iter + def force(self) -> float: + """Force in edge for current geometry.""" + return self._force_density * self.length + + @classmethod + def get_forces(cls, edges: Sequence[NaturalEdge]) -> List[float]: + """Get forces for sequence of edges.""" + forces = [.0] * cls.edge_count + for edge in edges: + forces[edge.edge_id] = edge.force + return forces + + @classmethod + def get_lengths(cls, edges: Sequence[NaturalEdge]) -> List[float]: + """Get lengths for sequence of edges.""" + return [edge.length for edge in edges] + + +class NaturalFace: + """Represents geometry and natural stress data of a single face. + + Parameters + ---------- + mesh_xyz : ndarray of float + XYZ vertex coordinates of the whole mesh as a (n, 3) array. + vertices_ids : sequence of int + Sequence of contiguous vertex indices. + stress_goal : sequence of float + Goal 2nd Piola-Kirchhoff (σx, σy, τxy) stress field for the face, + input over global directions and normalized over thickness. + (default is None, setting a uniform stress field (1, 1, 0)). + stress_ref : sequence of float + Normal of reference plane for non-isotropic stress field orientation. + is_subface : boolean + Whether the face is a tri subface of a quad face. + + Attributes + ---------- + face_count : int + The number of instantiated NaturalFaces. + xyz : ndarray of float + XYZ vertex coordinates of the face as a (n, 3) array. + """ + + face_count = 0 + + def __new__(cls, mesh_xyz: VertexArray, vertices_ids: Sequence[int], + stress_goal: Vec = None, stress_ref: Vec = None, + is_subface: bool = False) -> None: + v = len(vertices_ids) + if v == 3: + return super().__new__(TriFace) + elif v == 4: + return super().__new__(QuadFace) + else: + raise ValueError('Only tri and quad faces can be processed, ' + f'A face with {v} vertices was input.') + + def __init__(self, mesh_xyz: VertexArray, vertices_ids: Sequence[int], + stress_goal: Vec = None, stress_ref: Vec = None, + is_subface: bool = False) -> None: + self.cache = {} + self._set_ids(vertices_ids, is_subface) + self.xyz = mesh_xyz + self.stress_ref = stress_ref + self.stress_goal = stress_goal + + def update_xyz(self, mesh_xyz: VertexArray) -> None: + """Update XYZ face coordinates and clear geometry cache.""" + self.cache.clear() + self.xyz = mesh_xyz + + def _set_ids(self, vertices_ids: Sequence[int], is_subface: bool) -> None: + self.vertices_ids = vertices_ids + if is_subface: + self.face_id = -1 + else: + self.face_id = __class__.face_count + __class__.face_count += 1 + + def __len__(self) -> int: + return len(self.vertices_ids) + + @property + def xyz(self) -> VertexArray: + """XYZ vertex coordinates of the face.""" + return self._xyz + + @xyz.setter + def xyz(self, mesh_xyz: VertexArray) -> None: + self._xyz = [mesh_xyz[v] for v in self.vertices_ids] + + @property + @cache_iter + def stress_goal(self) -> NpVec: + """Goal 2nd Piola-Kirchhoff (σx, σy, τxy) stress field for the face.""" + if self.stress_ref is None or is_isotropic(self._stress_goal): + return self._stress_goal + return self._calc_stress_goal() + + @stress_goal.setter + def stress_goal(self, goal: Sequence[float]): + self._stress_goal = asarray([1, 1, 0]) if goal is None else asarray(goal) + if (not is_isotropic(self._stress_goal)) and (not self.stress_ref): + raise ValueError("Non-isotropic stress state requires reference normal.") + + def _calc_stress_goal(self) -> Vec: + """Calculate stress goal by intersection of the face with + the plane defined by stress reference vector as its normal.""" + f = self.frame + g = self._stress_goal + x, z, ref = f.xaxis, f.zaxis, asarray(self.stress_ref) + x_r = cross(z, ref) + θ = atan2(cross(x, x_r).dot(z), x.dot(x_r)) + return transform_stress_angle(g, θ) + + def transform_stress_to_global(self, vec: NpVec) -> NpVec: + """Transform stress pseudo-vector from local face frame to global frame.""" + R = asarray(Rotation.from_frame(self.frame))[:3, :3] + return R.dot(append(vec, .0)) + + @cache_iter + def principal_stress(self, to_global: bool = False) -> Tuple[NDArray[(2,), float64], NDArray[(2, 2), float64]]: + """Principal stress values and vectors from local stresses.""" + s = stress_vec_to_tensor(self.stress) + eigenvals, eigenvecs = eig(s) + # sort eigenvalues and corresponding vectors + ids = argsort(-eigenvals) + eigenvals, eigenvecs = eigenvals[ids], eigenvecs[ids, :] + if to_global: + eigenvecs = [self.transform_stress_to_global(v) for v in eigenvecs] + return eigenvals, eigenvecs + + face_stresses = namedtuple('stresses', ['amplitudes', 'vectors']) + + @classmethod + def get_stresses(cls, faces: Sequence[NaturalFace], stress_flag: int) -> NamedTuple: + """Strategy for calculating stresses for a sequence of faces.""" + s = cls.face_stresses + if stress_flag == -1: # only tri (sub)face stresses + return s(asarray(cls._get_tris_attr(faces, 'stress')), None) + if stress_flag == 0: # no stresses + return s(None, None) + if stress_flag == 1: # stresses in local frames + return s(asarray([f.stress for f in faces]), None) + if (stress_flag in (2, 3)): # principal stresses in local or global frame + eigs = (f.principal_stress(to_global=(stress_flag == 3)) for f in faces) + return s(*map(list, zip(*eigs))) + raise ValueError("stress_flag parameter must be set to {0, 1, 2, 3}.") + + @classmethod + def get_stress_goals(cls, faces: Sequence[NaturalFace]) -> NpVec: + """Get stress goals for multiple tri (sub)faces. For quad faces, + transformation into tri subface frames is hence not required.""" + return asarray(cls._get_tris_attr(faces, 'stress_goal')) + + @staticmethod + def _get_tris_attr(faces: Sequence[NaturalFace], attr: Any) -> List: + """Get attribute for multiple tri (sub)faces.""" + attrs = [] + for f in faces: + if len(f) == 3: + attrs.append(getattr(f, attr, None)) + else: + attrs.extend([getattr(t, attr, None) + for t in f.tri_faces]) + return attrs + + +class QuadFace(NaturalFace): + """Represents geometry and natural stress data of a single quad face. + A quad face is composed of 4 pairwise overlapping tri faces. + + Parameters + ---------- + mesh_xyz : ndarray of float + XYZ vertex coordinates of the whole mesh as a (n, 3) array. + vertices_ids : sequence of int + Sequence of contiguous vertex indices. + stress_goal : sequence of float + Goal 2nd Piola-Kirchhoff (σx, σy, τxy) stress field on the face, + input over global directions and normalized over thickness. + (default is None, setting a uniform stress field (1, 1, 0)). + stress_ref : sequence of float + Normal of reference plane for non-isotropic stress field orientation. + + Attributes + ---------- + frame : :class:`Frame` + Local coordinate frame. + rotations : generator of (2, 2) ndarrays of float + Planarized rotation matrices from quad to tri subfaces. + force_densities : list of float + 6 natural force densities of the edges. + stress : ndarray of float + Pseudo-vector of 2nd Piola-Kirchhoff stresses. + """ + + # vertex sequence of each the 4 tri subfaces of the quad face + __TRI_VERTICES = ((1, 3, 0), (3, 1, 2), (2, 0, 1), (0, 2, 3)) + + def __init__(self, mesh_xyz: VertexArray, vertices_ids: Sequence[int], + stress_goal: Vec = None, stress_ref: Vec = None) -> None: + super().__init__(mesh_xyz, vertices_ids, stress_goal, stress_ref) + self._set_tri_faces(mesh_xyz) + + def update_xyz(self, mesh_xyz: VertexArray) -> None: + """Update XYZ face coordinates and clear geometry cache, + both of quad face itself as of its composing tri subfaces.""" + super().update_xyz(mesh_xyz) + for t in self.tri_faces: + t.update_xyz(mesh_xyz) + + def _set_tri_faces(self, mesh_xyz: VertexArray) -> None: + """Construct the 4 tri subfaces of the quad face.""" + tri_vertices_ids = (itemgetter(*self.__TRI_VERTICES[i]) + (self.vertices_ids) for i in range(4)) + self.tri_faces = [TriFace(mesh_xyz, v_ids, self._stress_goal, + self.stress_ref, True) for v_ids in tri_vertices_ids] + + @property + @cache_iter + def area(self) -> float: + """Face area, as half sum of composing tri subface areas.""" + return sum(t.area for t in self.tri_faces) / 2 + + @property + @cache_iter + def frame(self) -> Frame: + """Face frame, with origin at centroid and + x-axis along midpoint of quad edge (1, 2).""" + v0, v1, v2, v3 = self.xyz + c0 = centroid_points(self.xyz) + c1 = centroid_points((v1, v2)) + c2 = centroid_points((v2, v3)) + return Frame.from_points(c0, c1, c2) + + @property + @cache_iter + def rotations(self) -> Generator[NDArray[(2, 2), float64], None, None]: + """Planar rotation matrices from quad face frame to + each of the 4 tri subface frames.""" + # angle of quad edge (0, 1) to quad x-axis + vec_mid = self.frame.xaxis + vec_edge = Vector.from_start_end(self.xyz[0], self.xyz[1]) + dθ = angle_vectors_signed(vec_edge, vec_mid, (0, 0, 1)) + + # in-plane angles of tri subfaces axes to quad axes + βa, γc = self.tri_faces[0].angles[1], self.tri_faces[2].angles[2] + θa = pi - βa - dθ + θb = -βa - dθ + θc = γc - dθ + θd = pi + γc - dθ + return (planar_rotation(θ) for θ in (θa, θb, θc, θd)) + + @property + def force_densities(self) -> NDArray[(6,), float64]: + """6 natural force densities of the quad face for current + geometry and goal 2nd Piola-Kirchhoff stresses.""" + na, nb, nc, nd = (t.force_densities for t in self.tri_faces) + sum_force_densities = (na[1] + nc[0], nb[0] + nc[1], nb[1] + nd[0], + na[0] + nd[1], na[2] + nb[2], nc[2] + nd[2]) + return asarray([n / 2 for n in sum_force_densities]) + + @property + @cache_iter + def stress(self) -> NpVec: + """Pseudo-vector of 2nd Piola-Kirchhoff stresses at current geometry.""" + tot_a = self.area * 2 + tris_a = (t.area for t in self.tri_faces) + tris_s = (t.stress for t in self.tri_faces) + quad_s = self._transform_stress_tris_to_quad(tris_s) + return sum(a * s for a, s in zip(tris_a, quad_s)) / tot_a + + def _transform_stress_tris_to_quad(self, tris_stress: NpVec + ) -> Generator[NpVec, None, None]: + """Transform stresses from tri subfaces frames to quad frame.""" + return (transform_stress(s, R) for s, R in zip(tris_stress, self.rotations)) + + +class TriFace(NaturalFace): + """Represents geometry and natural stress data of a single tri face. + + Parameters + ---------- + mesh_xyz : ndarray of float + XYZ vertex coordinates of the whole mesh as a (n, 3) array. + vertices_ids : sequence of int + Sequence of contiguous vertex indices. + stress_goal : sequence of float + Goal 2nd Piola-Kirchhoff (σx, σy, τxy) stress field for the face, + input over reference directions and normalized over thickness. + (default is None, setting a uniform stress field (1, 1, 0)). + stress_ref : sequence of float + Normal of reference plane for non-isotropic stress field orientation. + + Attributes + ---------- + frame : :class:`Frame` + Local coordinate frame. + edge_lengths : tuple of float + Ordered edge lengths. + angles : tuple of float + Angles in radians. + area: float + Face area. + transformation : ndarray of float + (3, 3) transformation matrix for tri face edges into stresses. + transformation_inv : ndarray of float + (3, 3) inverted transformation matrix for tri face stresses into edges. + force_densities : list of float + 3 natural force densities of the edges. + stress : ndarray of float + Pseudo-vector of 2nd Piola-Kirchhoff stresses. + """ + + @property + @cache_iter + def frame(self) -> Frame: + """Tri face frame, with origin at vertex 0 and x-axis along edge (0, 1).""" + return Frame.from_points(*self.xyz) + + @property + @cache_iter + def edge_lengths(self) -> Tuple[float, float, float]: + """Ordered edge lengths from euclidean distance between vertices.""" + ordered_xyz = self.xyz[1:] + self.xyz[:2] + return tuple(euclidean_distance(u, v) for u, v in pairwise(ordered_xyz)) + + @property + @cache_iter + def angles(self) -> Tuple[float, float, float]: + """Angles in radians, calculated by law of cosines from edge lengths.""" + l0, l1, l2 = self.edge_lengths + α = arc_cos((l0 ** 2 + l1 ** 2 - l2 ** 2) / (2 * l0 * l1)) + β = arc_cos((l1 ** 2 + l2 ** 2 - l0 ** 2) / (2 * l1 * l2)) + γ = pi - α - β + return (α, β, γ) + + @property + @cache_iter + def area(self) -> float: + """Face area.""" + l0, l1, l2 = self.edge_lengths + return l1 * l2 * sin(self.angles[1]) / 2 + + @property + @cache_iter + def transformation(self) -> NDArray[(3, 3), float64]: + """3x3 transformation matrix for tri face edge directions + into stresses in local frame. Its column vectors are unit + edge stresses described in the local face frame.""" + α, β, γ = self.angles + T = [[c2(γ), c2(β), 1], # noqa + [s2(γ), s2(β), 0], # noqa + [s(γ) * c(γ), -s(β) * c(β), 0]] # noqa + return asarray(T) + + @property + @check_singularity + def transformation_inv(self) -> NpVec: + """3x3 transformation matrix for tri face stresses + in local frame into edge directions.""" + return inv(self.transformation) + + @property + def force_densities(self) -> NpVec: + """3 natural force densities of the tri face for current + geometry and goal (2nd Piola-Kirchhoff) stresses.""" + a = self.area + g = self.stress_goal + Li2 = asarray([1 / (ln ** 2) for ln in self.edge_lengths]) + Ti = self.transformation_inv + self._force_densities = a * Li2 * Ti.dot(g) + return self._force_densities + + @property + @cache_iter + def stress(self) -> NpVec: + """Pseudo-vector of 2nd Piola-Kirchhoff stresses at current geometry.""" + L2 = asarray([ln ** 2 for ln in self.edge_lengths]) + a = self.area + T = self.transformation + # force densities of previous geometry + n = self._force_densities + stress = (T * L2).dot(n) / a + return stress + + +# ================================================= +# mesh data processing +# ================================================= + +def mesh_preprocess(mesh: Mesh, goals: Goals) -> Tuple[ + VertexArray, List[NaturalEdge], List[NaturalFace]]: + """Pre-process mesh to lists of elements. + + Parameters + ---------- + mesh : :class:`Mesh` + Instance of Mesh datastructure. + goals: :class:'Goals' + Instance of Goals class. + + Returns + ---------- + list of NaturalFace + Processed mesh faces. + list of NaturalEdge + Processed mesh edges. + ndarray of float + XYZ vertex coordinates of the whole mesh as a (n, 3) array. + """ + # vertex indices mapping + v_index = mesh.key_index() + mesh_xyz = [mesh.vertex_coordinates(v) for v in mesh.vertices()] + edges_vts = ([v_index[u], v_index[v]] for u, v in mesh.edges()) + faces_vts = ([v_index[w] for w in mesh.face_vertices(f)] + for f in mesh.faces()) + + # set natural faces and edges + faces = [NaturalFace(mesh_xyz, face_v, sp, goals.stress_ref) + for face_v, sp in zip(faces_vts, goals.stress)] + edges = [e for e in [NaturalEdge(mesh_xyz, edge_v, fp, qp) for edge_v, fp, qp + in zip(edges_vts, goals.forces, goals.force_densities)] if e] + + return asarray(mesh_xyz), edges, faces diff --git a/src/compas_fd/nfd/goals.py b/src/compas_fd/nfd/goals.py new file mode 100644 index 00000000..850209f3 --- /dev/null +++ b/src/compas_fd/nfd/goals.py @@ -0,0 +1,40 @@ +from __future__ import annotations + +from typing import Sequence +from typing_extensions import Annotated + +from compas.datastructures import Mesh + + +Vec = Annotated[Sequence[float], 3] + + +class Goals: + """Represents collection of geometric and stress goals for mesh elements. + + Parameters + ---------- + mesh : :class:`Mesh` + Instance of Mesh datastructure. + stress : sequence of sequence of float + Goal 2nd Piola-Kirchhoff (σx, σy, τxy) stress field per face. + forces : sequence of float + Goal force per edge (default is None). + force_densities : sequence of float + Goal force density per edge, overwrites force goals. + stress_ref : sequence of float + Normal of reference plane for non-isotropic stress field orientation + (Default is None). + """ + + def __init__(self, mesh: Mesh, stress: Sequence[Vec], forces: Sequence[float], + force_densities: Sequence[float], stress_ref: Vec = None): + self.f_count = mesh.number_of_faces() + self.stress = self._process_goal(stress, default=(1.0, 1.0, .0)) + self.force_densities = self._process_goal(force_densities) + self.forces = self._process_goal(forces) + self.stress_ref = stress_ref + + def _process_goal(self, goal, default=.0): + """Construct goals as per input or to predefined defaults.""" + return goal or [default for _ in range(self.f_count)] \ No newline at end of file diff --git a/src/compas_fd/nfd/math_utilities.py b/src/compas_fd/nfd/math_utilities.py new file mode 100644 index 00000000..9a58b6bc --- /dev/null +++ b/src/compas_fd/nfd/math_utilities.py @@ -0,0 +1,99 @@ +from __future__ import annotations + +from math import sin, cos, asin, acos, hypot +from numpy import asarray + +from typing import Sequence +from typing_extensions import Annotated +from nptyping import NDArray + +Vec = Annotated[Sequence[float], 3] +NpVec = NDArray[(3,), float] + + +__all__ = [ + 'arc_sin', + 'arc_cos', + 'euclidean_distance', + 'planar_rotation', + 'is_isotropic', + 'transform_stress', + 'transform_stress_angle', + 'stress_vec_to_tensor', + 'stress_tensor_to_vec' +] + + +s, c = sin, cos + + +def s2(angle: float) -> float: + """Calculate sine squared of angle in radians.""" + return s(angle) ** 2 + + +def c2(angle: float) -> float: + """Calculate cosine squared of angle in radians.""" + return c(angle) ** 2 + + +def arc_sin(num: float) -> float: + """Calculate inverse sine. Input values are bounded + between -.9999 and .9999 for numerical stability.""" + return asin(max(min(num, .9999), -.9999)) + + +def arc_cos(num: float) -> float: + """Calculate inverse cosine. Input values are bounded + between -.9999 and .9999 for numerical stability.""" + return acos(max(min(num, .9999), -.9999)) + + +def euclidean_distance(pt_a: Vec, pt_b: Vec) -> float: + """Calculate distance between two 3D points using the hypotenuse.""" + return hypot(pt_a[0] - pt_b[0], pt_a[1] - pt_b[1], pt_a[2] - pt_b[2]) + + +def planar_rotation(angle: float) -> NDArray[(2, 2), float]: + """Get the planar rotation matrix from an angle in radians.""" + return asarray([[c(angle), -s(angle)], + [s(angle), c(angle)]]) + + +def is_isotropic(vec: NpVec) -> bool: + """Check whether a planar stress pseudo-vector is isotropic.""" + return (vec[0] == vec[1]) and (vec[2] == 0) + + +def transform_stress(stress: NpVec, rotation: NDArray[(2, 2), float], + invert: bool = False) -> NpVec: + """Transform a planar stress pseudo-vector by a 2x2 rotation matrix. + Internally, the stress pseudo-vector is converted to a 2x2 tensor.""" + s = stress_vec_to_tensor(stress) + R = rotation + r = (R.dot(s).dot(R.T) if invert + else R.T.dot(s).dot(R)) + return stress_tensor_to_vec(r) + + +def transform_stress_angle(stress: NpVec, angle: float, invert: bool = False) -> NpVec: + """Transform a planar stress pseudo-vector by an angle in radians.""" + a = -angle if invert else angle + s2a = s2(a) + c2a = c2(a) + sca = s(a) * c(a) + T = [[c2a, s2a, 2 * sca], + [s2a, c2a, -2 * sca], + [-sca, sca, c2a - s2a]] + return asarray(T).dot(stress) + + +def stress_vec_to_tensor(vec: NpVec) -> NDArray[(2, 2), float64]: + """Convert planar stresses from pseudo-vector to tensor form.""" + return asarray([[vec[0], vec[2]], + [vec[2], vec[1]]]) + + +def stress_tensor_to_vec(tens: NDArray[(2, 2), float64]) -> NpVec: + """Convert planar stresses from tensor to pseudo-vector form.""" + return asarray([tens[0, 0], tens[1, 1], tens[0, 1]]) diff --git a/src/compas_fd/nfd/matrices.py b/src/compas_fd/nfd/matrices.py new file mode 100644 index 00000000..03e4a224 --- /dev/null +++ b/src/compas_fd/nfd/matrices.py @@ -0,0 +1,195 @@ +from __future__ import annotations + +from typing import Any, Sequence +from typing_extensions import Annotated +from nptyping import NDArray + +from numpy import asarray, copy, zeros, ix_, float64 +from scipy.sparse import coo_matrix + +from compas.geometry import Rotation +from compas_fd.nfd.geometry import NaturalFace, NaturalEdge + + +Vec = Annotated[Sequence[float], 3] + + +__all__ = [ + 'StiffnessMatrixAssembler', + 'LoadMatrixAssembler' +] + + +class StiffnessMatrixAssembler: + """Represents assembly of a global force density stiffness matrix. + A new instance of StiffnessMatrixAssembler should be generated at each + solver iteration, so that the stiffness matrix is rebuilt completely. + + Parameters + ---------- + free : sequence of int + Indices of free mesh vertices. + fixed : sequence of int + Indices of fixed mesh vertices. + edges : iterable of :class:`NaturalEdge` + Processed edges of the mesh. + faces : iterable of :class:`NaturalFace` + Processed faces of the mesh. + + Attributes + ---------- + matrix : ndarray + Full force density stiffness matrix. + free_matrix : ndarray + Force density stiffness matrix of free vertices. + fixed_matrix : ndarray + Force density stiffness matrix of fixed vertices. + """ + + def __init__(self, free: Sequence[int], fixed: Sequence[int], + edges: Sequence[NaturalEdge], faces: Sequence[NaturalFace]) -> None: + self.free = free + self.fixed = fixed + + self.data = [] + self.rows = [] + self.cols = [] + + self._add_faces(faces) + self._add_edges(edges) + + v_count = len(free + fixed) + self._mat = coo_matrix((self.data, (self.rows, self.cols)), + (v_count, v_count)).tocsr() + + def _add_faces(self, faces: Sequence[NaturalFace]) -> None: + """Call NaturalFaces to calculate their edge force densities, + and enter them in the intermediate coo matrix lists.""" + for face in faces: + fv_count = len(face) + if fv_count == 3: + self._add_tri_face(face) + elif fv_count == 4: + self._add_quad_face(face) + + def _add_tri_face(self, face: NaturalFace) -> None: + """Call TriFace to calculate its edge force densities, + and enter them in the intermediate coo matrix lists.""" + v0, v1, v2 = face.vertices_ids + n0, n1, n2 = face.force_densities + self.data += [n1 + n2, -n2, -n1, + -n2, n0 + n2, -n0, + -n1, -n0, n0 + n1] + self.rows += [v0] * 3 + [v1] * 3 + [v2] * 3 + self.cols += [v0, v1, v2] * 3 + + def _add_quad_face(self, face: NaturalFace) -> None: + """Call QuadFace to calculate its edge force densities, + and enter them in the intermediate coo matrix lists.""" + v0, v1, v2, v3 = face.vertices_ids + n0, n1, n2, n3, n4, n5 = face.force_densities + self.data += [n0 + n3 + n5, -n0, -n5, -n3, + -n0, n0 + n1 + n4, -n1, -n4, + -n5, -n1, n1 + n2 + n5, -n2, + -n3, -n4, -n2, n2 + n3 + n4] + self.rows += [v0] * 4 + [v1] * 4 + [v2] * 4 + [v3] * 4 + self.cols += [v0, v1, v2, v3] * 4 + + def _add_edges(self, edges: Sequence[NaturalEdge]) -> None: + """Call NaturalEdges to calculate their force densities, + and enter them in the intermediate coo matrix lists.""" + for edge in edges: + v0, v1 = edge.vertices_ids + n = edge.force_density + self.data += [n, n, -n, -n] + self.rows += [v0, v1, v0, v1] + self.cols += [v0, v1, v1, v0] + + @property + def matrix(self) -> NDArray[(Any, Any), float64]: + """Full force density stiffness matrix.""" + return self._mat + + @property + def free_matrix(self) -> NDArray[(Any, Any), float64]: + """Force Density stiffness matrix of free vertices.""" + return self._mat[ix_(self.free, self.free)] + + @property + def fixed_matrix(self) -> NDArray[(Any, Any), float64]: + """Force density stiffness matrix of fixed vertices.""" + return self._mat[ix_(self.free, self.fixed)] + + +class LoadMatrixAssembler: + """Represents assembly of a global load matrix. + + Parameters + ---------- + size : int + Row count of load matrix as number of vertices in the mesh. + faces : iterable of :class:`NaturalFace` + The processed faces of the mesh. + vertex_loads : sequence of sequence of float + The loads on the mesh vertices in global XYZ directions. + global_face_loads : sequence of sequence of float + The loads on the mesh faces in global XYZ directions. + local_face_loads : sequence of sequence of float + The loads on the mesh faces in local face XYZ directions. + Local loads get converted to global loads by rotation + from the local face frames. + + Attributes + ---------- + matrix : ndarray + The full load matrix as a 2D array. + + Notes + ----- + A LoadMatrixAssembler is to be instantiated once per session and to persist + over iterations. By calling instance method update(), the internal matrix + is updated corresponding to the geometry of the latest state. + """ + + def __init__(self, size: int, faces: Sequence[NaturalFace], vertex_loads: Sequence[Vec] = None, + global_face_loads: Sequence[Vec] = None, local_face_loads: Sequence[Vec] = None) -> None: + self.faces = faces + self._vertices_mat = self._load_mat(vertex_loads, zeros((size, 3))) + self._gfl = self._load_mat(global_face_loads) + self._lfl = self._load_mat(local_face_loads) + self._has_face_loads = ((self._gfl is not None) or (self._lfl is not None)) + self._mat = self._vertices_mat + self.update() + + @property + def matrix(self) -> NDArray[(Any, Any), float64]: + """Full vertex load matrix.""" + return self._mat + + def _load_mat(self, loads: Sequence[Vec], default: bool = None): + """Pre-process load matrices of each element type.""" + return asarray(loads).reshape(-1, 3) if loads is not None else default + + def update(self) -> None: + """Assemble all loads corresponding to current geometry into updated load matrix. + To be called at each iteration where geometry has been updated.""" + if not self._has_face_loads: + return self._mat + self._mat = copy(self._vertices_mat) + if self._gfl is not None: + self._add_face_loads(local=False) + if self._lfl is not None: + self._add_face_loads(local=True) + + def _add_face_loads(self, local: bool = False) -> None: + """Add all face loads with either global or local reference frame.""" + face_loads = self._lfl if local else self._gfl + for face, face_load in zip(self.faces, face_loads): + vertex_load = face_load * (face.area / len(face)) + if local: + R = asarray(Rotation.from_frame(face.frame))[:3, :3] + vertex_load = R.dot(vertex_load) + self._mat[face.vertices_ids, :] += vertex_load + + def _add_edge_loads(self) -> None: + raise NotImplementedError diff --git a/src/compas_fd/nfd/solvers.py b/src/compas_fd/nfd/solvers.py new file mode 100644 index 00000000..6648f0f6 --- /dev/null +++ b/src/compas_fd/nfd/solvers.py @@ -0,0 +1,176 @@ +from __future__ import annotations + +from typing import Sequence +from typing_extensions import Annotated + +from functools import partial + +from numpy import array, average +from numpy.linalg import norm +from scipy.sparse.linalg import spsolve + +from compas.datastructures import Mesh +from compas_fd.fd.result import Result + +from .geometry import NaturalFace, NaturalEdge, Goals, mesh_preprocess +from .matrices import StiffnessMatrixAssembler, LoadMatrixAssembler + + +Vec = Annotated[Sequence[float], 3] + + +__all__ = [ + 'nfd_ur_numpy', + 'nfd_numpy' +] + + +# ================================================= +# solver iterators +# ================================================= + +def nfd_ur_numpy(mesh: Mesh, fixed: Sequence[int], stress_goals: Sequence[float] = None, + force_goals: Sequence[float] = None, force_density_goals: Sequence[float] = None, + vertex_loads: Sequence[Vec] = None, global_face_loads: Sequence[Vec] = None, + local_face_loads: Sequence[Vec] = None, stress_flag: int = 1, stress_ref: Vec = None, + stress_tol: float = 1e-2, xyz_tol: float = 1e-2, kmax: int = 10) -> Result: + """Natural force density method with updated reference strategy. + Input stress fields are taken for the reference geometry + that is updated at each iteration by the natural force densities. + + Parameters + ---------- + mesh : :class:`Mesh` + Instance of Mesh datastructure. + fixed: sequence of int + Incices of fixed (anchored) vertices. + stress_goals : sequence of sequence of float + Goal 2nd Piola-Kirchhoff (σx, σy, τxy) stress field per face, + as input over reference directions and normalized over thickness. + (default is None, setting a uniform stress field (1, 1, 0)). + force_goals : sequence of float + Goal force per edge (default is None). + force_density_goals : sequence of float + Goal force density per edge, overwrites force goals (default is None). + vertex_loads : sequence of sequence of float + Global XYZ components of loads per vertex (default is None). + global_face_loads : sequence of sequence of float + Global XYZ components of loads per face area (default is None). + local_face_loads : sequence of sequence of float + Local face frame XYZ components of loads per face area (default is None). + stress_flag : {0, 1, 2, 3} + Flag for stress calculation at final solver iteration (default is 1). + 0: Do not calculate stresses. + 1: Cauchy stresses per face. + 2: Principal stress values and vectors per face. + 3: Principal stress values and vectors in global frame. + stress_ref : sequence of float + Normal of reference plane for non-isotropic stress field orientation. + stress_tol : float + Tolerance for averaged sum of squared errors + from stress vectors to goal stress vector (default is 1e-2). + xyz_tol : float + Tolerance for difference in coordinate displacements + between two consecutive iterations (default is 1e-2). + kmax : int + Maximum number of iterations (default is 10). + + Returns + ---------- + ndarray + XYZ vertex coordinates of the equilibrium geometry. + ndarray + Residual and reaction forces per vertex. + tuple + Stresses as (stress values, eigenvectors). + Content depends on the value passed to the stress_flag parameter. + list + Forces per edge. + + Notes + ----- + For more info, see [1]_, [2]_ + + References + ---------- + .. [1] Pauletti, R.M.O. and Pimenta, P.M., 2008. The natural force + density method for the shape finding of taut structures. + Computer Methods in Applied Mechanics and Engineering, + 197(49-50), pp.4419-4428. + + .. [2] Pauletti, R.M.O. and Fernandes, F.L., 2020. An outline of the natural + force density method and its extension to quadrilateral elements. + International Journal of Solids and Structures, 185, pp.423-438. + """ + # pre-process mesh data + goals = Goals(mesh, stress_goals, force_goals, force_density_goals, stress_ref) + xyz, edges, faces = mesh_preprocess(mesh, goals) + loads = LoadMatrixAssembler(xyz.shape[0], faces, vertex_loads, + global_face_loads, local_face_loads) + + if kmax == 1: + return _nfd_solve(xyz, fixed, faces, edges, loads, stress_flag) + + for k in range(kmax): + _xyz, r, f, lg, s = _nfd_solve(xyz, fixed, faces, edges, loads, -1) + stress_goals = NaturalFace.get_stress_goals(faces) + stress_res = average(norm(stress_goals - s.amplitudes, axis=1)) + xyz_Δ = max(norm(xyz - _xyz, axis=1)) + converged = (stress_res < stress_tol) or (xyz_Δ < xyz_tol) + xyz = _xyz + if converged: + break + + _output_message(converged, k, stress_res, xyz_Δ) + s = NaturalFace.get_stresses(faces, stress_flag) + + return Result(xyz, r, f, lg, s) + + +def _output_message(converged: bool, k: int, stress_res: float, xyz_Δ: float) -> None: + if converged: + print(f'Convergence reached after {k+1} iterations.') + else: + print(f'No convergence reached after {k+1} iterations.', + '\nAverage stress residual: ', round(stress_res, 5), + '\nMax displacement residual: ', round(xyz_Δ, 5)) + + +nfd_numpy = partial(nfd_ur_numpy, kmax=1) + + +# ================================================= +# solver +# ================================================= + +def _nfd_solve(xyz, fixed: Sequence[int], faces: Sequence[NaturalFace], + edges: Sequence[NaturalEdge], loads: Sequence[Vec], stress_flag: int) -> Result: + """Solve system for coordinates and dependent variables + using the one-shot natural force density method.""" + # pre-process vertex data + v = xyz.shape[0] + free = list(set(range(v)) - set(fixed)) + _xyz = array(xyz, copy=True) + + # assemble new stiffness matrix + sma = StiffnessMatrixAssembler(free, fixed, edges, faces) + D, Di, Df = sma.matrix, sma.free_matrix, sma.fixed_matrix + + # get updated load matrix + loads.update() + p = loads.matrix + + # solve for coordinates and update elements + _xyz[free] = spsolve(Di, p[free] - Df.dot(xyz[fixed])) + for face in faces: + face.update_xyz(_xyz) + for edge in edges: + edge.update_xyz(_xyz) + + # solve for dependent variables + r = p - D.dot(xyz) + f = NaturalEdge.get_forces(edges) + lg = NaturalEdge.get_lengths(edges) + s = NaturalFace.get_stresses(faces, stress_flag) + + return Result(_xyz, r, f, lg, s)