#
# Copyright (C) 2010-2026 The ESPResSo project
#
# This file is part of ESPResSo.
#
# ESPResSo is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# ESPResSo is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
import espressomd
from espressomd.interactions import OifLocalForces, OifGlobalForces
from .oif_utils import (
large_number, small_epsilon, discard_epsilon, custom_str, norm,
vec_distance, get_triangle_normal, area_triangle, angle_btw_triangles,
oif_calc_stretching_force, oif_calc_bending_force,
oif_calc_local_area_force, oif_calc_global_area_force, oif_calc_volume_force
)
[docs]
class FixedPoint:
"""
Represents mesh points, not connected to any ESPResSo particle.
"""
def __init__(self, pos, id):
if not isinstance(id, int):
raise TypeError("Id must be integer.")
if not ((len(pos) == 3) and isinstance(pos[0], float) and isinstance(
pos[1], float) and isinstance(pos[2], float)):
raise TypeError("Pos must be a list of three floats.")
self.x = pos[0]
self.y = pos[1]
self.z = pos[2]
self.id = id
[docs]
def get_pos(self):
return [self.x, self.y, self.z]
[docs]
def get_id(self):
return self.id
[docs]
class PartPoint:
"""
Represents mesh points, connected to ESPResSo particle.
"""
# part is physical ESPResSo particle corresponding to that particular point
def __init__(self, part, id, part_id):
if not (isinstance(part, espressomd.particle_data.ParticleHandle)
and isinstance(id, int) and isinstance(part_id, int)):
raise TypeError("Arguments to PartPoint are incorrect.")
self.part = part
self.part_id = part_id # because in adding bonds to the particles in OifCell
# one needs to know the global id of the particle.
self.id = id
[docs]
def get_pos(self):
return self.part.pos
[docs]
def get_vel(self):
return self.part.v
[docs]
def get_mass(self):
return self.part.mass
[docs]
def get_type(self):
return self.part.type
[docs]
def get_force(self):
return self.part.f
[docs]
def set_pos(self, pos):
self.part.pos = pos
[docs]
def set_vel(self, vel):
self.part.v = vel
[docs]
def set_force(self, force):
self.part.ext_force = force
[docs]
def kill_motion(self):
self.part.fix = [1, 1, 1]
[docs]
def unkill_motion(self):
self.part.unfix()
[docs]
class Edge:
"""
Represents edges in a mesh.
"""
def __init__(self, A, B):
if not all(isinstance(x, (PartPoint, FixedPoint)) for x in [A, B]):
raise TypeError(
"Arguments to Edge must be FixedPoint or PartPoint.")
self.A = A
self.B = B
[docs]
def length(self):
return vec_distance(self.A.get_pos(), self.B.get_pos())
[docs]
class Triangle:
"""
Represents triangles in a mesh.
"""
def __init__(self, A, B, C):
if not all(isinstance(x, (PartPoint, FixedPoint)) for x in [A, B, C]):
raise TypeError(
"Arguments to Triangle must be FixedPoint or PartPoint.")
self.A = A
self.B = B
self.C = C
[docs]
def area(self):
area = area_triangle(
self.A.get_pos(), self.B.get_pos(), self.C.get_pos())
return area
[docs]
class Angle:
"""
Represents angles in a mesh.
"""
def __init__(self, A, B, C, D):
if not all(isinstance(x, (PartPoint, FixedPoint))
for x in [A, B, C, D]):
raise TypeError(
"Arguments to Angle must be FixedPoint or PartPoint.")
self.A = A
self.B = B
self.C = C
self.D = D
[docs]
def size(self):
angle_size = angle_btw_triangles(
self.A.get_pos(), self.B.get_pos(), self.C.get_pos(), self.D.get_pos())
return angle_size
[docs]
class ThreeNeighbors:
"""
Represents three best spatially distributed neighbors of a point in a mesh.
"""
def __init__(self, A, B, C):
if not all(isinstance(x, (PartPoint, FixedPoint)) for x in [A, B, C]):
raise TypeError(
"Arguments to ThreeNeighbors must be FixedPoint or PartPoint.")
self.A = A
self.B = B
self.C = C
[docs]
def outer_normal(self):
outer_normal = get_triangle_normal(
self.A.get_pos(), self.B.get_pos(), self.C.get_pos())
return outer_normal
[docs]
class Mesh:
"""
Represents a triangular mesh.
"""
def __init__(
self, nodes_file=None, triangles_file=None, system=None, resize=(1.0, 1.0, 1.0),
particle_type=-1, particle_mass=1.0, normal=False, check_orientation=True):
if (system is None) or (not isinstance(system, espressomd.System)):
raise Exception(
"Mesh: No system provided or wrong type given.")
self.system = system
self.normal = normal
self.nodes_file = nodes_file
self.triangles_file = triangles_file
self.points = []
self.edges = []
self.triangles = []
self.angles = []
self.neighbors = []
self.ids_extremal_points = [0, 0, 0, 0, 0, 0, 0]
if not ((nodes_file is None) or (triangles_file is None)):
if not (isinstance(nodes_file, str)
and isinstance(triangles_file, str)):
raise TypeError("Mesh: Filenames must be strings.")
if not ((len(resize) == 3) and isinstance(resize[0], float) and isinstance(
resize[1], float) and isinstance(resize[2], float)):
raise TypeError("Mesh: Pos must be a list of three floats.")
if not isinstance(particle_type, int):
raise TypeError("Mesh: particle_type must be integer.")
if not isinstance(particle_mass, float):
raise TypeError("Mesh: particle_mass must be float.")
if not isinstance(normal, bool):
raise TypeError("Mesh: normal must be bool.")
if not isinstance(check_orientation, bool):
raise TypeError("Mesh: check_orientation must be bool.")
# reading the mesh point positions from file
in_file = open(nodes_file, "r")
nodes_coord = in_file.read().split("\n")
in_file.close()
# removes a blank line at the end of the file if there is any:
nodes_coord = filter(None, nodes_coord)
# here we have list of lines with triplets of strings
for line in nodes_coord: # extracts coordinates from the string line
line = np.array([float(x) for x in line.split()])
coords = np.array(resize) * line
tmp_fixed_point = FixedPoint(coords, len(self.points))
self.points.append(tmp_fixed_point)
# searching for extremal points IDs
x_min = large_number
x_max = -large_number
y_min = large_number
y_max = -large_number
z_min = large_number
z_max = -large_number
for tmp_fixed_point in self.points:
coords = tmp_fixed_point.get_pos()
if coords[0] < x_min:
x_min = coords[0]
self.ids_extremal_points[0] = tmp_fixed_point.get_id()
if coords[0] > x_max:
x_max = coords[0]
self.ids_extremal_points[1] = tmp_fixed_point.get_id()
if coords[1] < y_min:
y_min = coords[1]
self.ids_extremal_points[2] = tmp_fixed_point.get_id()
if coords[1] > y_max:
y_max = coords[1]
self.ids_extremal_points[3] = tmp_fixed_point.get_id()
if coords[2] < z_min:
z_min = coords[2]
self.ids_extremal_points[4] = tmp_fixed_point.get_id()
if coords[2] > z_max:
z_max = coords[2]
self.ids_extremal_points[5] = tmp_fixed_point.get_id()
# reading the triangle incidences from file
in_file = open(triangles_file, "r")
triangles_incid = in_file.read().split("\n")
in_file.close()
# removes a blank line at the end of the file if there is any:
triangles_incid = filter(None, triangles_incid)
for line in triangles_incid: # extracts incidences from the string line
incid = np.array([int(x) for x in line.split()])
tmp_triangle = Triangle(
self.points[incid[0]], self.points[incid[1]], self.points[incid[2]])
self.triangles.append(tmp_triangle)
if check_orientation is True:
# check whether all triangles in file had the same orientation;
# if not, correct the orientation
self.check_orientation()
# creating list of edge incidences from triangle incidences
# using temporary list of edge incidences
tmp_edge_incidences = []
for triangle in self.triangles:
pa = triangle.A.id
pb = triangle.B.id
pc = triangle.C.id
if ([pa, pb] not in tmp_edge_incidences) and (
[pb, pa] not in tmp_edge_incidences):
tmp_edge_incidences.append([pa, pb])
if ([pb, pc] not in tmp_edge_incidences) and (
[pc, pb] not in tmp_edge_incidences):
tmp_edge_incidences.append([pb, pc])
if ([pa, pc] not in tmp_edge_incidences) and (
[pc, pa] not in tmp_edge_incidences):
tmp_edge_incidences.append([pa, pc])
for tmp_incid in tmp_edge_incidences:
tmp_edge = Edge(
self.points[tmp_incid[0]], self.points[tmp_incid[1]])
self.edges.append(tmp_edge)
# creating list angles (former bending incidences) from triangle
# incidences
for edge in self.edges:
pa = edge.A.id
pb = edge.B.id
pc = -1
pd = -1
detected = 0
# detected = number of detected triangles with current edge common
# Algorithm is as follows: we run over all triangles and check
# whether two vertices are those from current edge. If we find such triangle,
# we put the ID of the third vertex to pc and we check if the orientation pa, pb, pc is the same as
# was in the triangle list (meaning, that we found one of the following three triples
# in the triangle list: pa, pb, pc or pb, pc, pa or pc, pa, pb).
# If we have the same orientation, we set orient = 1, otherwise orient = -1.
# Then we go further looking for the second triangle.
# The second triangle should have the opposite orientation.
# The normal of the first triangle will be P1P2 x P1P3, of the
# second triangle will be P2P4 x P2P3
orient = 0
for triangle in self.triangles:
# Run over all triangles and determine the two triangles
# with the common current edge
if (pa == triangle.A.id) and (pb == triangle.B.id):
if detected == 0:
# if no triangle with such edge was detected before
pc = triangle.C.id
detected = 1
orient = 1
else:
# if this is the second triangle with this edge,
# then also quit the for-loop over triangles
pd = triangle.C.id
break
if (pa == triangle.B.id) and (pb == triangle.C.id):
if detected == 0:
pc = triangle.A.id
detected = 1
orient = 1
else:
pd = triangle.A.id
break
if (pa == triangle.C.id) and (pb == triangle.A.id):
if detected == 0:
pc = triangle.B.id
detected = 1
orient = 1
else:
pd = triangle.B.id
break
if (pa == triangle.B.id) and (pb == triangle.A.id):
if detected == 0:
pc = triangle.C.id
detected = 1
orient = -1
else:
pd = triangle.C.id
break
if (pa == triangle.C.id) and (pb == triangle.B.id):
if detected == 0:
pc = triangle.A.id
detected = 1
orient = -1
else:
pd = triangle.A.id
break
if (pa == triangle.A.id) and (pb == triangle.C.id):
if detected == 0:
pc = triangle.B.id
detected = 1
orient = -1
else:
pd = triangle.B.id
break
if orient == 1:
tmp = pd
pd = pc
pc = tmp
tmp_angle = Angle(
self.points[pc], self.points[pa], self.points[pb], self.points[pd])
self.angles.append(tmp_angle)
# creating list of three neighbors for membrane collision
if normal is True:
for point in self.points:
tmp_neighbors = []
# cycle through edges and select those that contain point
for edge in self.edges:
# take an edge and copy the nodes of the edge to pa, pb
if edge.A.id == point.id:
tmp_neighbors.append(edge.B)
if edge.B.id == point.id:
tmp_neighbors.append(edge.A)
# create vectors to all neighbors and normalize them
tmp_vectors_to_neighbors = []
p_coords = np.array(point.get_pos())
for neighbor in tmp_neighbors:
tmp_vector = neighbor.get_pos() - p_coords
tmp_length = norm(tmp_vector)
if tmp_length < small_epsilon:
raise Exception("Mesh: Degenerate edge.")
tmp_vector /= tmp_length
tmp_vectors_to_neighbors.append(tmp_vector)
# check all triplets of neighbors and select the one that is best spatially distributed
# by adding the corresponding three normalized vectors
# and selecting the one with smallest resultant vector
n_neighbors = len(tmp_neighbors)
min_length = large_number
best_neighbors = [
tmp_neighbors[0], tmp_neighbors[1], tmp_neighbors[2]]
for i in range(0, n_neighbors):
for j in range(i + 1, n_neighbors):
for k in range(j + 1, n_neighbors):
tmp_result_vector = tmp_vectors_to_neighbors[i] + tmp_vectors_to_neighbors[j] + \
tmp_vectors_to_neighbors[k]
tmp_result_vector_length = norm(
tmp_result_vector)
if tmp_result_vector_length < min_length:
min_length = tmp_result_vector_length
best_neighbors = [
tmp_neighbors[i], tmp_neighbors[j], tmp_neighbors[k]]
# find one triangle that contains this point and compute
# its normal vector
for triangle in self.triangles:
if triangle.A.id == point.id or triangle.B.id == point.id or triangle.C.id == point.id:
tmp_normal_triangle = get_triangle_normal(
triangle.A.get_pos(), triangle.B.get_pos(),
triangle.C.get_pos())
break
# properly orient selected neighbors and save them to the
# list of neighbors
tmp_normal_neighbors = get_triangle_normal(
best_neighbors[
0].get_pos(), best_neighbors[1].get_pos(),
best_neighbors[2].get_pos())
tmp_length_normal_triangle = norm(tmp_normal_triangle)
tmp_length_normal_neighbors = norm(tmp_normal_neighbors)
tmp_product = np.dot(tmp_normal_triangle, tmp_normal_neighbors) / \
(tmp_length_normal_triangle *
tmp_length_normal_neighbors)
tmp_angle = np.arccos(tmp_product)
if tmp_angle > np.pi / 2.0:
selected_neighbors = ThreeNeighbors(
best_neighbors[0], best_neighbors[1], best_neighbors[2])
else:
selected_neighbors = ThreeNeighbors(
best_neighbors[0], best_neighbors[2], best_neighbors[1])
self.neighbors.append(selected_neighbors)
else:
for point in self.points:
selected_neighbors = ThreeNeighbors(point, point, point)
self.neighbors.append(selected_neighbors)
[docs]
def copy(self, origin=None, particle_type=-1, particle_mass=1.0,
rotate=None):
"""
Create particles in the system.
"""
mesh = Mesh(system=self.system)
mesh.ids_extremal_points = self.ids_extremal_points
if rotate is not None:
# variables for rotation
ca = np.cos(rotate[0])
sa = np.sin(rotate[0])
cb = np.cos(rotate[1])
sb = np.sin(rotate[1])
cc = np.cos(rotate[2])
sc = np.sin(rotate[2])
rotation = np.array(
[[cb * cc, sa * sb * cc - ca * sc, sc * sa + cc * sb * ca],
[cb * sc, ca * cc + sa * sb * sc, sc * sb * ca - cc * sa],
[-sb, cb * sa, ca * cb]])
for point in self.points:
# PartPoints are created
tmp_pos = point.get_pos()
tmp_rotate_pos = np.array(point.get_pos())
# rotation of nodes
if rotate is not None:
tmp_pos = rotation.dot(tmp_rotate_pos)
tmp_pos = [discard_epsilon(tmp_pos[0]), discard_epsilon(
tmp_pos[1]), discard_epsilon(tmp_pos[2])]
if origin is not None:
tmp_pos += np.array(origin)
# to remember the global id of the ESPResSo particle
new_part_id = len(self.system.part)
new_part = self.system.part.add(
pos=tmp_pos, type=particle_type, mass=particle_mass,
mol_id=particle_type, id=new_part_id)
new_part_point = PartPoint(new_part, len(mesh.points), new_part_id)
mesh.points.append(new_part_point)
for edge in self.edges:
new_edge = Edge(mesh.points[edge.A.id], mesh.points[edge.B.id])
mesh.edges.append(new_edge)
for triangle in self.triangles:
new_triangle = Triangle(
mesh.points[triangle.A.id], mesh.points[triangle.B.id], mesh.points[triangle.C.id])
mesh.triangles.append(new_triangle)
for angle in self.angles:
new_angle = Angle(
mesh.points[angle.A.id], mesh.points[
angle.B.id], mesh.points[angle.C.id],
mesh.points[angle.D.id])
mesh.angles.append(new_angle)
for neighbors in self.neighbors:
new_neighbors = ThreeNeighbors(
mesh.points[neighbors.A.id], mesh.points[neighbors.B.id],
mesh.points[neighbors.C.id])
mesh.neighbors.append(new_neighbors)
return mesh
[docs]
def check_orientation(self):
tmp_triangle_list = []
tmp_triangle_list_ok = []
t_ok = None
corrected_triangle = None
for triangle in self.triangles:
tmp_triangle_list.append(triangle)
# move the first triangle to the checked and corrected list
tmp_triangle_list_ok.append(tmp_triangle_list[0])
tmp_triangle_list.pop(0)
while tmp_triangle_list:
i = 0
while i < len(tmp_triangle_list):
tmp_triangle = tmp_triangle_list[i]
for correct_triangle in tmp_triangle_list_ok:
# check if triangles have a common edge, if so, check
# orientation
are_neighbors = True
if tmp_triangle.A.id == correct_triangle.A.id:
if tmp_triangle.B.id == correct_triangle.B.id:
t_ok = False # this is situation 123 and 124
corrected_triangle = Triangle(
tmp_triangle.A, tmp_triangle.C, tmp_triangle.B)
else:
if tmp_triangle.B.id == correct_triangle.C.id:
t_ok = True # this is situation 123 and 142
else:
if tmp_triangle.C.id == correct_triangle.B.id:
t_ok = True # this is situation 123 and 134
else:
if tmp_triangle.C.id == correct_triangle.C.id:
t_ok = False # this is situation 123 and 143
corrected_triangle = Triangle(
tmp_triangle.A, tmp_triangle.C, tmp_triangle.B)
else:
are_neighbors = False
else:
if tmp_triangle.A.id == correct_triangle.B.id:
if tmp_triangle.B.id == correct_triangle.C.id:
t_ok = False # this is situation 123 and 412
corrected_triangle = Triangle(
tmp_triangle.A, tmp_triangle.C, tmp_triangle.B)
else:
if tmp_triangle.B.id == correct_triangle.A.id:
t_ok = True # this is situation 123 and 214
else:
if tmp_triangle.C.id == correct_triangle.C.id:
t_ok = True # this is situation 123 and 413
else:
if tmp_triangle.C.id == correct_triangle.A.id:
t_ok = False # this is situation 123 and 314
corrected_triangle = Triangle(
tmp_triangle.A, tmp_triangle.C,
tmp_triangle.B)
else:
are_neighbors = False
else:
if tmp_triangle.A.id == correct_triangle.C.id:
if tmp_triangle.B.id == correct_triangle.A.id:
t_ok = False # this is situation 123 and 241
corrected_triangle = Triangle(
tmp_triangle.A, tmp_triangle.C, tmp_triangle.B)
else:
if tmp_triangle.B.id == correct_triangle.B.id:
t_ok = True # this is situation 123 and 421
else:
if tmp_triangle.C.id == correct_triangle.A.id:
t_ok = True # this is situation 123 and 341
else:
if tmp_triangle.C.id == correct_triangle.B.id:
t_ok = False # this is situation 123 and 431
corrected_triangle = Triangle(
tmp_triangle.A, tmp_triangle.C,
tmp_triangle.B)
else:
are_neighbors = False
else:
if tmp_triangle.B.id == correct_triangle.A.id:
if tmp_triangle.C.id == correct_triangle.B.id:
t_ok = False # this is situation 123 and 234
corrected_triangle = Triangle(
tmp_triangle.A, tmp_triangle.C, tmp_triangle.B)
else:
if tmp_triangle.C.id == correct_triangle.C.id:
t_ok = True # this is situation 123 and 243
else:
are_neighbors = False
else:
if tmp_triangle.B.id == correct_triangle.B.id:
if tmp_triangle.C.id == correct_triangle.C.id:
t_ok = False # this is situation 123 and 423
corrected_triangle = Triangle(
tmp_triangle.A, tmp_triangle.C,
tmp_triangle.B)
else:
if tmp_triangle.C.id == correct_triangle.A.id:
t_ok = True # this is situation 123 and 324
else:
are_neighbors = False
else:
if tmp_triangle.B.id == correct_triangle.C.id:
if tmp_triangle.C.id == correct_triangle.A.id:
t_ok = False # this is situation 123 and 342
corrected_triangle = Triangle(
tmp_triangle.A, tmp_triangle.C,
tmp_triangle.B)
else:
if tmp_triangle.C.id == correct_triangle.B.id:
t_ok = True # this is situation 123 and 432
else:
are_neighbors = False
else:
are_neighbors = False
if are_neighbors:
# move the tmp_triangle to the checked and corrected
# list
if t_ok:
tmp_triangle_list_ok.append(tmp_triangle)
else:
tmp_triangle_list_ok.append(corrected_triangle)
tmp_triangle_list.pop(i)
break
i += 1
# replace triangles with checked triangles
i = 0
for tmp_triangle in tmp_triangle_list_ok:
self.triangles[i] = Triangle(
tmp_triangle.A, tmp_triangle.C, tmp_triangle.B)
i += 1
# all triangles now have the same orientation, check if it is correct
tmp_volume = self.volume()
if tmp_volume < 0:
# opposite orientation, flip all triangles
i = 0
for tmp_triangle in self.triangles:
self.triangles[i] = Triangle(
tmp_triangle.A, tmp_triangle.C, tmp_triangle.B)
i += 1
return 0
[docs]
def surface(self):
surface = 0.0
for triangle in self.triangles:
surface += triangle.area()
return surface
[docs]
def volume(self):
volume = 0.0
for triangle in self.triangles:
tmp_normal = get_triangle_normal(
triangle.A.get_pos(), triangle.B.get_pos(), triangle.C.get_pos())
tmp_normal_length = norm(tmp_normal)
tmp_sum_z_coords = 1.0 / 3.0 * \
(triangle.A.get_pos()[2] + triangle.B.get_pos()[2] +
triangle.C.get_pos()[2])
volume -= (triangle.area() * tmp_normal[2] / tmp_normal_length *
tmp_sum_z_coords)
return volume
[docs]
def get_n_nodes(self):
"""
Get the number of mesh points.
"""
return len(self.points)
[docs]
def get_n_triangles(self):
"""
Get the number of triangles.
"""
return len(self.triangles)
[docs]
def get_n_edges(self):
"""
Get the number of edges.
"""
return len(self.edges)
[docs]
def output_mesh_triangles(self, triangles_file=None):
"""
Write the current triangles into a file that
can be used as input in the next simulations.
This is useful after checking orientation,
if any of the triangles were corrected.
"""
# this is useful after the mesh correction
# output of mesh nodes can be done from OifCell (this is because their
# position may change)
if triangles_file is None:
raise Exception(
"OifMesh: No file_name provided for triangles.")
output_file = open(triangles_file, "w")
for t in self.triangles:
output_file.write(
str(t.A.id) + " " + str(t.B.id) + " " + str(t.C.id) + "\n")
output_file.close()
return 0
[docs]
def mirror(self, mirror_x=0, mirror_y=0, mirror_z=0, out_file_name=""):
if out_file_name == "":
raise Exception(
"Cell.Mirror: output meshnodes file for new mesh is missing.")
if mirror_x not in (0, 1) or mirror_y not in (
0, 1) or mirror_z not in (0, 1):
raise Exception(
"Mesh.Mirror: for mirroring only values 0 or 1 are accepted. 1 indicates that the corresponding coordinate will be flipped.")
if mirror_x + mirror_y + mirror_z > 1:
raise Exception(
"Mesh.Mirror: flipping allowed only for one axis.")
if mirror_x + mirror_y + mirror_z == 1:
out_file = open(out_file_name, "w")
for p in self.points:
coor = p.get_pos()
if mirror_x == 1:
coor[0] *= -1.0
if mirror_y == 1:
coor[1] *= -1.0
if mirror_z == 1:
coor[2] *= -1.0
out_file.write(custom_str(coor[0]) + " " + custom_str(
coor[1]) + " " + custom_str(coor[2]) + "\n")
out_file.close()
return 0
[docs]
class OifCellType:
"""
Template to create elastic objects ("cells") from a mesh file.
The switches ``ks``, ``kb`` and ``kal`` set elastic parameters for
local interactions: ``ks`` for edge stiffness, ``kb`` for angle
preservation stiffness and ``kal`` for triangle area preservation
stiffness. Currently, the stiffness is implemented to be uniform over
the whole object, but with some tweaking, it is possible to have
non-uniform local interactions.
At least one of the elastic moduli should be set.
Note the difference between stretching (``ks``) and linear stretching
(``kslin``): these two options cannot be used simultaneously.
Linear stretching behaves like linear spring, where the stretching
force is calculated as :math:`\\mathbf{f_{ij}} = k_s \\mathbf{d_{ij}}`,
where :math:`\\mathbf{d_{ij}}` is the prolongation of the given edge.
By default, the stretching is non-linear (neo-Hookean).
The optional viscous damping (``kvisc``) is meant to reduce oscillations
in the stretching force. This corresponds to the case :math:`\\gamma^T = 0`
and :math:`\\gamma^C = k_{\\mathrm{visc}}` in equation 8 of :cite:`fedosov10a`
(which uses the :math:`\\times` symbol for the dot product).
Parameters
----------
nodes_file : :obj:`str`
Path to the mesh data file. Each line contains three real numbers.
These are the *x, y, z* coordinates of individual surface mesh nodes
of the objects centered at [0,0,0] and normalized
so that the "radius" of the object is 1.
triangles_file : :obj:`str`
Path to the triangles data file. Each line contains three integers
representing the 0-indexed line number of the mesh nodes
in ``nodes_file`` that form a triangle.
system : :obj:`espressomd.system.System`
System in which particles will be created.
ks : :obj:`float`
Elastic modulus for stretching forces.
kslin : :obj:`float`
Elastic modulus for linear stretching forces.
kb : :obj:`float`
Elastic modulus for bending forces.
kal : :obj:`float`
Elastic modulus for local area forces.
kvisc : :obj:`float`
Viscous damping for stretching forces.
kag : :obj:`float`
Elastic modulus for global area forces.
kv : :obj:`float`
Elastic modulus for volume forces.
resize : (3,) array_like of :obj:`float`
Scaling factor by which the coordinates stored in ``nodes_file``
will be stretched in the *x, y, z* directions.
normal : :obj:`bool`
If ``True``, create list of neighbors to allow membrane collisions,
and thus cell-cell interactions.
check_orientation : :obj:`bool`
If ``True``, check whether ``triangles_file`` contains triangles
with correct orientation. If not, it corrects the orientation
and created cells with corrected triangles.
It is useful for new or unknown meshes, but not necessary for meshes
that have already been tried out. Since it can take a few minutes for
larger meshes (with thousands of nodes), it can be set to ``False``.
In that case, the check is skipped when creating the ``CellType`` and
a warning is displayed.
The order of indices in ``triangles_file`` is important. Normally, each
triangle ABC should be oriented in such a way, that the normal vector
computed as vector product ABxAC must point inside the object.
For example, a sphere (or any other sufficiently convex object) contains
such triangles that the normals of these triangles point towards the
center of the sphere (almost).
The check runs over all triangles, makes sure that they have the
correct orientation and then calculates the volume of the object. If
the result is negative, it flips the orientation of all triangles.
Note, this method tells the user about the correction it makes.
If there is any, it might be useful to save the corrected triangulation
for future simulations using the method
:meth:`Mesh.output_mesh_triangles`, so that the
check does not have to be used repeatedly.
"""
def __init__(
self, nodes_file="", triangles_file="", system=None, resize=(1.0, 1.0, 1.0), ks=0.0, kslin=0.0,
kb=0.0, kal=0.0, kag=0.0, kv=0.0, kvisc=0.0, normal=False, check_orientation=True):
if (system is None) or (not isinstance(system, espressomd.System)):
raise Exception(
"OifCellType: No system provided or wrong type")
if (nodes_file == "") or (triangles_file == ""):
raise Exception(
"OifCellType: Either nodes_file or triangles_file is missing")
if not (isinstance(nodes_file, str)
and isinstance(triangles_file, str)):
raise TypeError("OifCellType: Filenames must be strings.")
if not ((len(resize) == 3) and isinstance(resize[0], float) and isinstance(
resize[1], float) and isinstance(resize[2], float)):
raise TypeError(
"OifCellType: Resize must be a list of three floats.")
if not (isinstance(ks, float) and isinstance(kslin, float) and isinstance(kb, float) and isinstance(
kal, float) and isinstance(kag, float) and isinstance(kv, float) and isinstance(kvisc, float)):
raise TypeError("OifCellType: Elastic parameters must be floats.")
if not isinstance(normal, bool):
raise TypeError("OifCellType: normal must be bool.")
if not isinstance(check_orientation, bool):
raise TypeError("OifCellType: check_orientation must be bool.")
if (ks != 0.0) and (kslin != 0.0):
raise Exception(
"OifCellType: Cannot use linear and nonlinear stretching at the same time.")
self.system = system
self.mesh = Mesh(
nodes_file=nodes_file, triangles_file=triangles_file, system=system, resize=resize,
normal=normal, check_orientation=check_orientation)
self.local_force_interactions = []
self.resize = resize
self.ks = ks
self.kslin = kslin
self.kb = kb
self.kal = kal
self.kag = kag
self.kv = kv
self.kvisc = kvisc
self.normal = normal
if (ks != 0.0) or (kslin != 0.0) or (kb != 0.0) or (kal != 0.0):
for angle in self.mesh.angles:
r0 = vec_distance(angle.B.get_pos(), angle.C.get_pos())
phi = angle_btw_triangles(
angle.A.get_pos(), angle.B.get_pos(), angle.C.get_pos(), angle.D.get_pos())
area1 = area_triangle(
angle.A.get_pos(), angle.B.get_pos(), angle.C.get_pos())
area2 = area_triangle(
angle.D.get_pos(), angle.B.get_pos(), angle.C.get_pos())
tmp_local_force_inter = OifLocalForces(
r0=r0, ks=ks, kslin=kslin, phi0=phi, kb=kb, A01=area1, A02=area2,
kal=kal, kvisc=kvisc)
self.local_force_interactions.append(
[tmp_local_force_inter, [angle.A, angle.B, angle.C, angle.D]])
self.system.bonded_inter.add(tmp_local_force_inter)
if (kag != 0.0) or (kv != 0.0):
surface = self.mesh.surface()
volume = self.mesh.volume()
self.global_force_interaction = OifGlobalForces(
A0_g=surface, ka_g=kag, V0=volume, kv=kv)
self.system.bonded_inter.add(self.global_force_interaction)
[docs]
def print_info(self):
"""
Print information about this template.
"""
print("\nThe following OifCellType was created: ")
print("\t nodes_file: " + self.mesh.nodes_file)
print("\t triangles_file: " + self.mesh.triangles_file)
print("\t n_nodes: " + str(self.mesh.get_n_nodes()))
print("\t n_triangles: " + str(self.mesh.get_n_triangles()))
print("\t n_edges: " + str(self.mesh.get_n_edges()))
print("\t ks: " + custom_str(self.ks))
print("\t kslin: " + custom_str(self.kslin))
print("\t kb: " + custom_str(self.kb))
print("\t kal: " + custom_str(self.kal))
print("\t kag: " + custom_str(self.kag))
print("\t kv: " + custom_str(self.kv))
print("\t kvisc: " + custom_str(self.kvisc))
print("\t normal: " + str(self.normal))
print("\t resize: " + str(self.resize))
print(" ")
[docs]
class OifCell:
"""
Represent a concrete elastic object ("cell").
Parameters
----------
cell_type : :obj:`OifCellType`
Template containing mesh nodes, triangle incidences,
elasticity parameters and bond parameters, and a system,
into which particles and bonds will be instantiated.
particle_type : :obj:`int`
Must start at 0 for the first cell and monotonically increase for
subsequent cells. Volume calculation of individual objects and
interactions between objects are set up using these types.
origin : (3,) array_like of :obj:`float`
Center of the cell object.
particle_mass : :obj:`float`
Mass of each particle forming the cell.
rotate : (3,) array_like of :obj:`float`
angles in radians, by which the cell will be rotated about the
*x, y, z* axis. Default value is (0.0, 0.0, 0.0).
Value (:math:`\\pi/2, 0.0, 0.0`) means that the object will
be rotated by :math:`\\pi/2` radians clockwise around the *x*
axis when looking in the positive direction of the axis.
"""
def __init__(self, cell_type=None, origin=None, particle_type=None,
particle_mass=1.0, rotate=None):
if (cell_type is None) or (not isinstance(cell_type, OifCellType)):
raise Exception(
"OifCell: No cellType provided or wrong type.")
if (origin is None) or \
(not ((len(origin) == 3) and isinstance(origin[0], float) and isinstance(origin[1], float) and isinstance(origin[2], float))):
raise TypeError("Origin must be tuple.")
if (particle_type is None) or (not isinstance(particle_type, int)):
raise Exception(
"OifCell: No particle_type specified or wrong type.")
if not isinstance(particle_mass, float):
raise Exception("OifCell: particle mass must be float.")
if (rotate is not None) and not ((len(rotate) == 3) and isinstance(
rotate[0], float) and isinstance(rotate[1], float) and isinstance(rotate[2], float)):
raise TypeError("Rotate must be list of three floats.")
self.cell_type = cell_type
self.cell_type.system.max_oif_objects = self.cell_type.system.max_oif_objects + 1
self.mesh = cell_type.mesh.copy(
origin=origin, particle_type=particle_type, particle_mass=particle_mass, rotate=rotate)
self.particle_mass = particle_mass
self.particle_type = particle_type
self.origin = origin
self.rotate = rotate
for inter in self.cell_type.local_force_interactions:
esp_inter = inter[0]
points = inter[1]
n_points = len(points)
if n_points == 2:
p0 = self.mesh.points[
points[0].id] # Getting PartPoints from id's of FixedPoints
p1 = self.mesh.points[points[1].id]
p0.part.add_bond((esp_inter, p1.part_id))
if n_points == 3:
p0 = self.mesh.points[points[0].id]
p1 = self.mesh.points[points[1].id]
p2 = self.mesh.points[points[2].id]
p0.part.add_bond((esp_inter, p1.part_id, p2.part_id))
if n_points == 4:
p0 = self.mesh.points[points[0].id]
p1 = self.mesh.points[points[1].id]
p2 = self.mesh.points[points[2].id]
p3 = self.mesh.points[points[3].id]
p1.part.add_bond(
(esp_inter, p0.part_id, p2.part_id, p3.part_id))
if (self.cell_type.kag != 0.0) or (self.cell_type.kv != 0.0):
for triangle in self.mesh.triangles:
triangle.A.part.add_bond(
(self.cell_type.global_force_interaction, triangle.B.part_id,
triangle.C.part_id))
[docs]
def get_origin(self):
"""
Get the cell center in unfolded coordinates.
"""
center = np.array([0.0, 0.0, 0.0])
for p in self.mesh.points:
center += p.get_pos()
return center / len(self.mesh.points)
[docs]
def set_origin(self, new_origin=(0.0, 0.0, 0.0)):
"""
Set the cell center.
"""
old_origin = self.get_origin()
for p in self.mesh.points:
new_position = p.get_pos() - old_origin + new_origin
p.set_pos(new_position)
[docs]
def get_approx_origin(self):
"""
Get the approximate location of the cell center.
It is computed as average of 6 mesh points that have extremal *x*,
*y* and *z* coordinates at the time of object loading.
"""
approx_center = np.array([0.0, 0.0, 0.0])
for pid in self.mesh.ids_extremal_points:
approx_center += self.mesh.points[pid].get_pos()
return approx_center / len(self.mesh.ids_extremal_points)
[docs]
def get_origin_folded(self):
"""
Get the cell center in folded coordinates.
"""
origin = self.get_origin()
return np.mod(origin, self.cell_type.system.box_l)
[docs]
def get_velocity(self):
"""
Get the average velocity of the cell mesh points.
"""
velocity = np.array([0.0, 0.0, 0.0])
for p in self.mesh.points:
velocity += p.get_vel()
return velocity / len(self.mesh.points)
[docs]
def set_velocity(self, new_velocity=(0.0, 0.0, 0.0)):
"""
Set the velocity of all cell mesh points.
"""
for p in self.mesh.points:
p.set_vel(new_velocity)
[docs]
def pos_bounds(self):
"""
Compute six extremal coordinates of the cell.
More precisely, run through all mesh points and return
the minimal and maximal :math:`x`-coordinate, :math:`y`-coordinate and
:math:`z`-coordinate in the order (:math:`x_{max}`, :math:`x_{min}`,
:math:`y_{max}`, :math:`y_{min}`, :math:`z_{max}`, :math:`z_{min}`).
"""
x_min = large_number
x_max = -large_number
y_min = large_number
y_max = -large_number
z_min = large_number
z_max = -large_number
for p in self.mesh.points:
coords = p.get_pos()
if coords[0] < x_min:
x_min = coords[0]
if coords[0] > x_max:
x_max = coords[0]
if coords[1] < y_min:
y_min = coords[1]
if coords[1] > y_max:
y_max = coords[1]
if coords[2] < z_min:
z_min = coords[2]
if coords[2] > z_max:
z_max = coords[2]
return [x_min, x_max, y_min, y_max, z_min, z_max]
[docs]
def surface(self):
"""
Get the cell surface.
"""
return self.mesh.surface()
[docs]
def volume(self):
"""
Get the cell volume.
"""
return self.mesh.volume()
[docs]
def diameter(self):
"""
Compute the maximal pairwise distance between all mesh points.
Computational complexity is :math:`\\mathcal{O}(N^2)`.
"""
n_points = len(self.mesh.points)
pos = np.array(list(map(lambda p: p.get_pos(), self.mesh.points)))
max_dist_sq = 0.0
for i in range(0, n_points - 1):
dist_sq = np.max((np.square(pos[i] - pos[i + 1:])).sum(axis=1))
if dist_sq > max_dist_sq:
max_dist_sq = dist_sq
return np.sqrt(max_dist_sq)
[docs]
def get_n_nodes(self):
"""
Get the number of mesh points.
"""
return self.mesh.get_n_nodes()
[docs]
def set_force(self, new_force=(0.0, 0.0, 0.0)):
"""
Set an external force to all mesh points.
Note, that this command sets the external force in each integration
step. So if you want to use the external force only in one iteration,
you need to set zero external force in the following integration step.
"""
for p in self.mesh.points:
p.set_force(new_force)
[docs]
def output_vtk_pos(self, file_name=None):
"""
Write the mesh information to a VTK file in unfolded coordinates.
ParaView can directly visualize this file.
"""
if file_name is None:
raise Exception("OifCell: No file_name provided for vtk output.")
n_points = len(self.mesh.points)
n_triangles = len(self.mesh.triangles)
output_file = open(file_name, "w")
output_file.write("# vtk DataFile Version 3.0\n")
output_file.write("Data\n")
output_file.write("ASCII\n")
output_file.write("DATASET POLYDATA\n")
output_file.write("POINTS " + str(n_points) + " float\n")
for p in self.mesh.points:
coords = p.get_pos()
output_file.write(custom_str(coords[0]) + " " + custom_str(
coords[1]) + " " + custom_str(coords[2]) + "\n")
output_file.write("TRIANGLE_STRIPS " + str(
n_triangles) + " " + str(4 * n_triangles) + "\n")
for t in self.mesh.triangles:
output_file.write(
"3 " + str(t.A.id) + " " + str(t.B.id) + " " + str(t.C.id) + "\n")
output_file.close()
[docs]
def output_vtk_pos_folded(self, file_name=None):
"""
Write the mesh information to a VTK file in folded coordinates.
ParaView can directly visualize this file.
"""
if file_name is None:
raise Exception(
"OifCell: No file_name provided for vtk output.")
n_points = len(self.mesh.points)
n_triangles = len(self.mesh.triangles)
# get coordinates of the origin
center = np.array([0.0, 0.0, 0.0])
for p in self.mesh.points:
center += p.get_pos()
center /= len(self.mesh.points)
center_folded = np.floor(center / self.cell_type.system.box_l)
# this gives how many times the origin is folded in all three
# directions
output_file = open(file_name, "w")
output_file.write("# vtk DataFile Version 3.0\n")
output_file.write("Data\n")
output_file.write("ASCII\n")
output_file.write("DATASET POLYDATA\n")
output_file.write("POINTS " + str(n_points) + " float\n")
for p in self.mesh.points:
coords = p.get_pos() - center_folded * self.cell_type.system.box_l
output_file.write(custom_str(coords[0]) + " " + custom_str(
coords[1]) + " " + custom_str(coords[2]) + "\n")
output_file.write("TRIANGLE_STRIPS " + str(
n_triangles) + " " + str(4 * n_triangles) + "\n")
for t in self.mesh.triangles:
output_file.write(
"3 " + str(t.A.id) + " " + str(t.B.id) + " " + str(t.C.id) + "\n")
output_file.close()
[docs]
def append_point_data_to_vtk(self, file_name=None, data_name=None,
data=None, first_append=None):
"""
Append the specified scalar ``data`` to a VTK file.
This is useful for ParaView visualisation of local velocity magnitudes,
magnitudes of forces, ..., in the meshnodes and can be shown in
ParaView by selecting the ``data_name`` in the *Properties* toolbar.
It is possible to consecutively write multiple datasets;
for the first one, the ``first_append`` parameter is set to ``True``,
for the following datasets, it needs to be set to ``False``.
This is to ensure the proper structure of the output file.
"""
if file_name is None:
raise Exception(
"OifCell: append_point_data_to_vtk: No file_name provided.")
if data is None:
raise Exception(
"OifCell: append_point_data_to_vtk: No data provided.")
if data_name is None:
raise Exception(
"OifCell: append_point_data_to_vtk: No data_name provided.")
if first_append is None:
raise Exception("OifCell: append_point_data_to_vtk: Need to know whether this is the first data list to be "
"appended for this file.")
n_points = self.get_n_nodes()
if len(data) != n_points:
raise Exception(
"OifCell: append_point_data_to_vtk: Number of data points does not match number of mesh points.")
output_file = open(file_name, "a")
if first_append is True:
output_file.write("POINT_DATA " + str(n_points) + "\n")
output_file.write("SCALARS " + data_name + " float 1\n")
output_file.write("LOOKUP_TABLE default\n")
for p in self.mesh.points:
output_file.write(str(data[p.id]) + "\n")
output_file.close()
[docs]
def output_raw_data(self, file_name=None, data=None):
"""
Write the vector or matrix ``rawdata`` to a white-space delimited
text file. The rows are ordered by the mesh points indices.
"""
if file_name is None:
raise Exception(
"OifCell: output_raw_data: No file_name provided.")
if data is None:
raise Exception(
"OifCell: output_raw_data: No data provided.")
n_points = self.get_n_nodes()
if len(data) != n_points:
raise Exception(
"OifCell: output_raw_data: Number of data points does not match number of mesh points.")
output_file = open(file_name, "w")
for p in self.mesh.points:
output_file.write(" ".join(map(str, data[p.id])) + "\n")
output_file.close()
[docs]
def output_mesh_points(self, file_name=None):
"""
Write the positions of the mesh points to a white-space delimited
text file. This file can later be read by :meth:`set_mesh_points`.
The center of the object is located at (0.0, 0.0, 0.0).
This method is meant to store a deformed shape in order
to be loaded later as a new shape template.
"""
if file_name is None:
raise Exception(
"OifCell: No file_name provided for mesh nodes output.")
output_file = open(file_name, "w")
center = self.get_origin()
for p in self.mesh.points:
coords = p.get_pos() - center
output_file.write(custom_str(coords[0]) + " " + custom_str(
coords[1]) + " " + custom_str(coords[2]) + "\n")
output_file.close()
[docs]
def set_mesh_points(self, file_name=None):
"""
Read a file containing a list of mesh point coordinates to deform
the current cell. The current origin stays unchanged. The file should
contain the coordinates of the mesh points with the origin location at
(0.0, 0.0, 0.0). It can be generated by :meth:`output_mesh_points`.
"""
if file_name is None:
raise Exception(
"OifCell: No file_name provided for set_mesh_points.")
center = self.get_origin()
n_points = self.get_n_nodes()
in_file = open(file_name, "r")
nodes_coord = in_file.read().split("\n")
in_file.close()
# removes a blank line at the end of the file if there is any:
nodes_coord = filter(None, nodes_coord)
# here we have list of lines with triplets of strings
if len(nodes_coord) != n_points:
raise Exception("OifCell: Mesh nodes not set to new positions: "
"number of lines in the file does not equal number of Cell nodes.")
else:
i = 0
for line in nodes_coord: # extracts coordinates from the string line
line = line.split()
new_position = np.array(line).astype(float) + center
self.mesh.points[i].set_pos(new_position)
i += 1
[docs]
def print_info(self):
"""
Print information about the cell.
"""
print("\nThe following OifCell was created: ")
print("\t particle_mass: " + custom_str(self.particle_mass))
print("\t particle_type: " + str(self.particle_type))
print("\t rotate: " + str(self.rotate))
print("\t origin: " + str(self.origin[0]) + " " + str(
self.origin[1]) + " " + str(self.origin[2]))
[docs]
def elastic_forces(
self, el_forces=(0, 0, 0, 0, 0, 0), f_metric=(0, 0, 0, 0, 0, 0), vtk_file=None,
raw_data_file=None):
"""
This method can be used in two different ways. One is to compute
the elastic forces locally for each mesh point, and the other is to
compute the f-metric, which is an approximation of elastic energy.
To compute the elastic forces, use the vector ``el_forces``.
It is a sextuple of zeros and ones, e.g. ``el_forces = (1,0,0,1,0,0)``,
where the ones denote the elastic forces to be computed. The order
is (stretching, bending, local area, global area, volume, total).
The output can be saved in two different ways: either by setting
``vtk_file = filename.vtk``, which saves a VTK file that can be
visualized using ParaView. If more than one elastic force was
selected, they can be chosen in the Properties window in ParaView.
The other type of output is ``raw_data_file*=filename.dat``, which will
save a datafile with the selected type of elastic force - one force
per row, where each row corresponds to a single mesh node. Note that
only one type of elastic force can be written this way at a time.
Thus, if you need output for several elastic forces, this method
should be called several times.
To compute the f-metric, use the vector ``f_metric``. It is also
a sextuple of zeros and ones, e.g. ``f_metric = (1,1,0,0,0,0)``,
where the ones denote the elastic forces to be computed. The order
is (stretching, bending, local area, global area, volume, total).
The output is also a vector with six elements, each corresponding
to the requested f-metric/"naive energy" computed as a
sum of magnitudes of respective elastic forces over all mesh points.
"""
# the order of parameters in elastic_forces and in f_metric is as follows (ks, kb, kal, kag, kv, total)
# vtk_file means that a vtk file for visualisation of elastic forces will be written
# raw_data_file means that just the elastic forces will be written into
# the output file
stretching_forces_list = []
bending_forces_list = []
local_area_forces_list = []
global_area_forces_list = []
volume_forces_list = []
elastic_forces_list = []
stretching_forces_norms_list = []
bending_forces_norms_list = []
local_area_forces_norms_list = []
global_area_forces_norms_list = []
volume_forces_norms_list = []
elastic_forces_norms_list = []
ks_f_metric = 0.0
kb_f_metric = 0.0
kal_f_metric = 0.0
kag_f_metric = 0.0
kv_f_metric = 0.0
total_f_metric = 0.0
for i in range(0, 6):
if (el_forces[i] != 0) and (el_forces[i] != 1):
raise Exception("OifCell: elastic_forces: Incorrect argument. el_forces has to be a sixtuple of 0s and 1s, "
"specifying which elastic forces will be calculated. The order in the sixtuple is (ks, kb, "
"kal, kag, kv, total).")
for i in range(0, 6):
if (f_metric[i] != 0) and (f_metric[i] != 1):
raise Exception("OifCell: elastic_forces: Incorrect argument. f_metric has to be a sixtuple of 0s and 1s, "
"specifying which f_metric will be calculated. The order in the sixtuple is (ks, kb, kal, "
"kag, kv, total)")
# calculation of stretching forces and f_metric
if (el_forces[0] == 1) or (el_forces[5] == 1) or (
f_metric[0] == 1) or (f_metric[5] == 1):
# initialize list
stretching_forces_list = np.zeros((len(self.mesh.points), 3))
# calculation uses edges, but results are stored for nodes
for e in self.mesh.edges:
a_current_pos = e.A.get_pos()
b_current_pos = e.B.get_pos()
a_orig_pos = self.cell_type.mesh.points[e.A.id].get_pos()
b_orig_pos = self.cell_type.mesh.points[e.B.id].get_pos()
current_dist = e.length()
orig_dist = vec_distance(a_orig_pos, b_orig_pos)
tmp_stretching_force = oif_calc_stretching_force(
self.cell_type.ks, a_current_pos, b_current_pos,
orig_dist, current_dist)
stretching_forces_list[e.A.id] += tmp_stretching_force
stretching_forces_list[e.B.id] -= tmp_stretching_force
# calculation of stretching f_metric, if needed
if f_metric[0] == 1:
ks_f_metric = 0.0
for p in self.mesh.points:
ks_f_metric += norm(stretching_forces_list[p.id])
# calculation of bending forces and f_metric
if (el_forces[1] == 1) or (el_forces[5] == 1) or (
f_metric[1] == 1) or (f_metric[5] == 1):
# initialize list
bending_forces_list = np.zeros((len(self.mesh.points), 3))
# calculation uses bending incidences, but results are stored for
# nodes
for angle in self.mesh.angles:
a_current_pos = angle.A.get_pos()
b_current_pos = angle.B.get_pos()
c_current_pos = angle.C.get_pos()
d_current_pos = angle.D.get_pos()
a_orig_pos = self.cell_type.mesh.points[angle.A.id].get_pos()
b_orig_pos = self.cell_type.mesh.points[angle.B.id].get_pos()
c_orig_pos = self.cell_type.mesh.points[angle.C.id].get_pos()
d_orig_pos = self.cell_type.mesh.points[angle.D.id].get_pos()
current_angle = angle.size()
orig_angle = angle_btw_triangles(
a_orig_pos, b_orig_pos, c_orig_pos, d_orig_pos)
tmp_bending_forces = oif_calc_bending_force(
self.cell_type.kb, a_current_pos, b_current_pos, c_current_pos,
d_current_pos, orig_angle, current_angle)
tmp_bending_force1 = np.array(
[tmp_bending_forces[0], tmp_bending_forces[1], tmp_bending_forces[2]])
tmp_bending_force2 = np.array(
[tmp_bending_forces[3], tmp_bending_forces[4], tmp_bending_forces[5]])
bending_forces_list[angle.A.id] += tmp_bending_force1
bending_forces_list[angle.B.id] -= 0.5 * \
tmp_bending_force1 + 0.5 * tmp_bending_force2
bending_forces_list[angle.C.id] -= 0.5 * \
tmp_bending_force1 + 0.5 * tmp_bending_force2
bending_forces_list[angle.D.id] += tmp_bending_force2
# calculation of bending f_metric, if needed
if f_metric[1] == 1:
kb_f_metric = 0.0
for p in self.mesh.points:
kb_f_metric += norm(bending_forces_list[p.id])
# calculation of local area forces and f_metric
if (el_forces[2] == 1) or (el_forces[5] == 1) or (
f_metric[2] == 1) or (f_metric[5] == 1):
# initialize list
local_area_forces_list = np.zeros((len(self.mesh.points), 3))
# calculation uses triangles, but results are stored for nodes
for t in self.mesh.triangles:
a_current_pos = t.A.get_pos()
b_current_pos = t.B.get_pos()
c_current_pos = t.C.get_pos()
a_orig_pos = self.cell_type.mesh.points[t.A.id].get_pos()
b_orig_pos = self.cell_type.mesh.points[t.B.id].get_pos()
c_orig_pos = self.cell_type.mesh.points[t.C.id].get_pos()
current_area = t.area()
orig_area = area_triangle(a_orig_pos, b_orig_pos, c_orig_pos)
tmp_local_area_forces = oif_calc_local_area_force(
self.cell_type.kal, a_current_pos, b_current_pos,
c_current_pos, orig_area, current_area)
local_area_forces_list[t.A.id] += np.array(
[tmp_local_area_forces[0], tmp_local_area_forces[1],
tmp_local_area_forces[2]])
local_area_forces_list[t.B.id] += np.array(
[tmp_local_area_forces[3], tmp_local_area_forces[4],
tmp_local_area_forces[5]])
local_area_forces_list[t.C.id] += np.array(
[tmp_local_area_forces[6], tmp_local_area_forces[7],
tmp_local_area_forces[8]])
# calculation of local area f_metric, if needed
if f_metric[2] == 1:
kal_f_metric = 0.0
for p in self.mesh.points:
kal_f_metric += norm(local_area_forces_list[p.id])
# calculation of global area forces and f_metric
if (el_forces[3] == 1) or (el_forces[5] == 1) or (
f_metric[3] == 1) or (f_metric[5] == 1):
# initialize list
global_area_forces_list = np.zeros((len(self.mesh.points), 3))
# calculation uses triangles, but results are stored for nodes
for t in self.mesh.triangles:
a_current_pos = t.A.get_pos()
b_current_pos = t.B.get_pos()
c_current_pos = t.C.get_pos()
current_surface = self.mesh.surface()
orig_surface = self.cell_type.mesh.surface()
tmp_global_area_forces = oif_calc_global_area_force(
self.cell_type.kag, a_current_pos, b_current_pos,
c_current_pos, orig_surface, current_surface)
global_area_forces_list[t.A.id] += np.array(
[tmp_global_area_forces[0], tmp_global_area_forces[1],
tmp_global_area_forces[2]])
global_area_forces_list[t.B.id] += np.array(
[tmp_global_area_forces[3], tmp_global_area_forces[4],
tmp_global_area_forces[5]])
global_area_forces_list[t.C.id] += np.array(
[tmp_global_area_forces[6], tmp_global_area_forces[7],
tmp_global_area_forces[8]])
# calculation of global area f_metric, if needed
if f_metric[3] == 1:
kag_f_metric = 0.0
for p in self.mesh.points:
kag_f_metric += norm(global_area_forces_list[p.id])
# calculation of volume forces and f_metric
if (el_forces[4] == 1) or (el_forces[5] == 1) or (
f_metric[4] == 1) or (f_metric[5] == 1):
# initialize list
volume_forces_list = np.zeros((len(self.mesh.points), 3))
# calculation uses triangles, but results are stored for nodes
for t in self.mesh.triangles:
a_current_pos = t.A.get_pos()
b_current_pos = t.B.get_pos()
c_current_pos = t.C.get_pos()
current_volume = self.mesh.volume()
orig_volume = self.cell_type.mesh.volume()
tmp_volume_force = oif_calc_volume_force(
self.cell_type.kv, a_current_pos, b_current_pos, c_current_pos,
orig_volume, current_volume)
volume_forces_list[t.A.id] += tmp_volume_force
volume_forces_list[t.B.id] += tmp_volume_force
volume_forces_list[t.C.id] += tmp_volume_force
# calculation of volume f_metric, if needed
if f_metric[4] == 1:
kv_f_metric = 0.0
for p in self.mesh.points:
kv_f_metric += norm(volume_forces_list[p.id])
# calculation of total elastic forces and f_metric
if (el_forces[5] == 1) or (f_metric[5] == 1):
elastic_forces_list = []
for p in self.mesh.points:
total_elastic_forces = stretching_forces_list[p.id] + bending_forces_list[p.id] + \
local_area_forces_list[p.id] + global_area_forces_list[p.id] + \
volume_forces_list[p.id]
elastic_forces_list.append(total_elastic_forces)
# calculation of total f_metric, if needed
if f_metric[5] == 1:
total_f_metric = 0.0
for p in self.mesh.points:
total_f_metric += norm(elastic_forces_list[p.id])
# calculate norms of resulting forces
if sum(el_forces) != 0:
if el_forces[0] == 1:
stretching_forces_norms_list = []
for p in self.mesh.points:
stretching_forces_norms_list.append(
norm(stretching_forces_list[p.id]))
if el_forces[1] == 1:
bending_forces_norms_list = []
for p in self.mesh.points:
bending_forces_norms_list.append(
norm(bending_forces_list[p.id]))
if el_forces[2] == 1:
local_area_forces_norms_list = []
for p in self.mesh.points:
local_area_forces_norms_list.append(
norm(local_area_forces_list[p.id]))
if el_forces[3] == 1:
global_area_forces_norms_list = []
for p in self.mesh.points:
global_area_forces_norms_list.append(
norm(global_area_forces_list[p.id]))
if el_forces[4] == 1:
volume_forces_norms_list = []
for p in self.mesh.points:
volume_forces_norms_list.append(
norm(volume_forces_list[p.id]))
if el_forces[5] == 1:
elastic_forces_norms_list = []
for p in self.mesh.points:
elastic_forces_norms_list.append(
norm(elastic_forces_list[p.id]))
# output vtk (folded)
if vtk_file is not None:
if sum(el_forces) == 0:
raise Exception("OifCell: elastic_forces: The option elastic_forces was not used. "
"Nothing to output to vtk file.")
self.output_vtk_pos_folded(vtk_file)
first = True
if el_forces[0] == 1:
self.append_point_data_to_vtk(
file_name=vtk_file, data_name="ks_f_metric",
data=stretching_forces_norms_list, first_append=first)
first = False
if el_forces[1] == 1:
self.append_point_data_to_vtk(
file_name=vtk_file, data_name="kb_f_metric",
data=bending_forces_norms_list, first_append=first)
first = False
if el_forces[2] == 1:
self.append_point_data_to_vtk(
file_name=vtk_file, data_name="kal_f_metric",
data=local_area_forces_norms_list, first_append=first)
first = False
if el_forces[3] == 1:
self.append_point_data_to_vtk(
file_name=vtk_file, data_name="kag_f_metric",
data=global_area_forces_norms_list, first_append=first)
first = False
if el_forces[4] == 1:
self.append_point_data_to_vtk(
file_name=vtk_file, data_name="kav_f_metric",
data=volume_forces_norms_list, first_append=first)
first = False
if el_forces[5] == 1:
self.append_point_data_to_vtk(
file_name=vtk_file, data_name="total_f_metric",
data=elastic_forces_norms_list, first_append=first)
# output raw data
if raw_data_file is not None:
if sum(el_forces) != 1:
raise Exception("OifCell: elastic_forces: Only one type of elastic forces can be written into one "
"raw_data_file. If you need several, please call OifCell.elastic_forces multiple times - "
"once per elastic force.")
if el_forces[0] == 1:
self.output_raw_data(
file_name=raw_data_file, data=stretching_forces_list)
if el_forces[1] == 1:
self.output_raw_data(
file_name=raw_data_file, data=bending_forces_list)
if el_forces[2] == 1:
self.output_raw_data(
file_name=raw_data_file, data=local_area_forces_list)
if el_forces[3] == 1:
self.output_raw_data(
file_name=raw_data_file, data=global_area_forces_list)
if el_forces[4] == 1:
self.output_raw_data(
file_name=raw_data_file, data=volume_forces_list)
if el_forces[5] == 1:
self.output_raw_data(
file_name=raw_data_file, data=elastic_forces_list)
# return f_metric
if sum(f_metric) > 0:
results = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
if f_metric[0] == 1:
results[0] = ks_f_metric
if f_metric[1] == 1:
results[1] = kb_f_metric
if f_metric[2] == 1:
results[2] = kal_f_metric
if f_metric[3] == 1:
results[3] = kag_f_metric
if f_metric[4] == 1:
results[4] = kv_f_metric
if f_metric[5] == 1:
results[5] = total_f_metric
return results
else:
return 0