Commit 90b6d221 authored by Daniel Boeckenhoff's avatar Daniel Boeckenhoff
Browse files

mesh3D cutting works fine with doctest

parent 5518e2cd
...@@ -751,7 +751,7 @@ class Tensors(AbstractNdarray): ...@@ -751,7 +751,7 @@ class Tensors(AbstractNdarray):
rtol=None, atol=None, equal_nan=False, rtol=None, atol=None, equal_nan=False,
return_bool=True): return_bool=True):
""" """
Test, whether the instance has the same content as other. Evaluate, whether the instance has the same content as other.
Args: Args:
optional: optional:
rtol (float) rtol (float)
...@@ -762,15 +762,18 @@ class Tensors(AbstractNdarray): ...@@ -762,15 +762,18 @@ class Tensors(AbstractNdarray):
if issubclass(type(other), Tensors) and self.coordSys != other.coordSys: if issubclass(type(other), Tensors) and self.coordSys != other.coordSys:
other = other.copy() other = other.copy()
other.transform(self.coordSys) other.transform(self.coordSys)
x, y = np.asarray(self), np.asarray(other)
if rtol is None and atol is None: if rtol is None and atol is None:
if return_bool: mask = (x == y)
return np.array_equal(self, other) if equal_nan:
return self == other both_nan = np.isnan(x) & np.isnan(y)
elif rtol is None: mask[both_nan] = both_nan[both_nan]
rtol = 0. else:
elif atol is None: if rtol is None:
atol = 0. rtol = 0.
mask = np.isclose(self, other, rtol=rtol, atol=atol, equal_nan=equal_nan) if atol is None:
atol = 0.
mask = np.isclose(x, y, rtol=rtol, atol=atol, equal_nan=equal_nan)
if return_bool: if return_bool:
return bool(np.all(mask)) return bool(np.all(mask))
return mask return mask
......
...@@ -789,13 +789,10 @@ class Mesh3D(tfields.TensorMaps): ...@@ -789,13 +789,10 @@ class Mesh3D(tfields.TensorMaps):
if not any(~mask) or all(~mask): if not any(~mask) or all(~mask):
inst = inst[mask] inst = inst[mask]
elif at_intersection == 'split' or at_intersection == 'splitRough': elif at_intersection == 'split' or at_intersection == 'splitRough':
# add points and faces intersecting with the plane
expression_parts = tfields.lib.symbolics.split_expression(expression)
''' '''
define a new mesh that will be merged with the existing one add vertices and faces that are at the border of the cuts
the new mesh will describe the faces that are at the border of the
cuts
''' '''
expression_parts = tfields.lib.symbolics.split_expression(expression)
if len(expression_parts) > 1: if len(expression_parts) > 1:
new_mesh = inst.copy() new_mesh = inst.copy()
if at_intersection == 'splitRough': if at_intersection == 'splitRough':
...@@ -863,23 +860,21 @@ class Mesh3D(tfields.TensorMaps): ...@@ -863,23 +860,21 @@ class Mesh3D(tfields.TensorMaps):
""" """
Add the intersection points and faces Add the intersection points and faces
""" """
print "______________"
newP = _intersect(triangle_points, plane, vertices_rejected) newP = _intersect(triangle_points, plane, vertices_rejected)
last_idx = len(vertices) - 1 last_idx = len(vertices) - 1
for tri_list in newP: for tri_list in newP:
print vertices new_face = []
face = []
for item in tri_list: for item in tri_list:
if isinstance(item, int): if isinstance(item, int):
# reference to old vertex # reference to old vertex
face.append(item) new_face.append(face[item])
elif isinstance(item, complex): elif isinstance(item, complex):
# reference to new vertex that has been # reference to new vertex that has been
# concatenated already # concatenated already
face.append(last_idx + int(item.imag)) new_face.append(last_idx + int(item.imag))
else: else:
# new vertex # new vertex
face.append(len(vertices)) new_face.append(len(vertices))
vertices = np.append(vertices, vertices = np.append(vertices,
[map(float, item)], [map(float, item)],
axis=0) axis=0)
...@@ -887,13 +882,9 @@ class Mesh3D(tfields.TensorMaps): ...@@ -887,13 +882,9 @@ class Mesh3D(tfields.TensorMaps):
np.full((1,) + field.shape[1:], np.nan), np.full((1,) + field.shape[1:], np.nan),
axis=0) axis=0)
for field in fields] for field in fields]
fields[-1][-1] = len(vertices) faces = np.append(faces, [new_face], axis=0)
print vertices
print newP
print face
faces = np.append(faces, [face], axis=0)
faces_fields = [np.append(field, faces_fields = [np.append(field,
np.full((1,) + field.shape[1:], np.nan), [field[i]],
axis=0) axis=0)
for field in faces_fields] for field in faces_fields]
faces_fields[-1][-1] = i faces_fields[-1][-1] = i
...@@ -912,9 +903,9 @@ class Mesh3D(tfields.TensorMaps): ...@@ -912,9 +903,9 @@ class Mesh3D(tfields.TensorMaps):
inst.maps[0] = inst.maps[0][mask] inst.maps[0] = inst.maps[0][mask]
else: else:
raise ValueError("Sympy expression is not splitable.") raise ValueError("Sympy expression is not splitable.")
inst = inst.cleaned()
elif at_intersection == 'remove': elif at_intersection == 'remove':
pass inst = inst[mask]
else: else:
raise AttributeError("No at_intersection method called {at_intersection} " raise AttributeError("No at_intersection method called {at_intersection} "
"implemented".format(**locals())) "implemented".format(**locals()))
...@@ -922,7 +913,6 @@ class Mesh3D(tfields.TensorMaps): ...@@ -922,7 +913,6 @@ class Mesh3D(tfields.TensorMaps):
if _in_recursion: if _in_recursion:
template = None template = None
else: else:
print inst.fields
template_field = inst.fields.pop(-1) template_field = inst.fields.pop(-1)
template_maps = [] template_maps = []
for mp in inst.maps: for mp in inst.maps:
...@@ -932,7 +922,6 @@ class Mesh3D(tfields.TensorMaps): ...@@ -932,7 +922,6 @@ class Mesh3D(tfields.TensorMaps):
template = tfields.Mesh3D(tfields.Tensors(inst), template = tfields.Mesh3D(tfields.Tensors(inst),
template_field, template_field,
maps=template_maps) maps=template_maps)
print template.fields
return inst, template return inst, template
def _cut_template(self, template): def _cut_template(self, template):
...@@ -956,7 +945,7 @@ class Mesh3D(tfields.TensorMaps): ...@@ -956,7 +945,7 @@ class Mesh3D(tfields.TensorMaps):
>>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]], >>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]],
... [1, 0]) ... [1, 0])
>>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3], >>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3],
... [1, 0, 3, 2, 4] ... [1, 0, 3, 2, 4],
... maps=[tmap]) ... maps=[tmap])
Use template as instruction to make a fast cut Use template as instruction to make a fast cut
...@@ -970,11 +959,24 @@ class Mesh3D(tfields.TensorMaps): ...@@ -970,11 +959,24 @@ class Mesh3D(tfields.TensorMaps):
""" """
# Redirect fields # Redirect fields
fields = []
if template.fields: if template.fields:
fields = [field[template.fields[0].astype(int)] template_field = np.array(template.fields[0])
for field in self.fields] if len(self) > 0:
else: '''
fields = [] if new vertices have been created in the template, it is
in principle unclear what fields we have to refer to.
Thus in creating the template, we gave np.nan.
To make it fast, we replace nan with 0 as a dummy and correct
the field entries afterwards with np.nan.
'''
nan_mask = np.isnan(template_field)
template_field[nan_mask] = 0 # dummy reference to index 0.
template_field = template_field.astype(int)
for field in self.fields:
projected_field = field[template_field]
projected_field[nan_mask] = np.nan # correction for nan
fields.append(projected_field)
# Redirect maps fields # Redirect maps fields
maps = [] maps = []
...@@ -1042,7 +1044,7 @@ class Mesh3D(tfields.TensorMaps): ...@@ -1042,7 +1044,7 @@ class Mesh3D(tfields.TensorMaps):
>>> float(m_split[:, 0].min()) >>> float(m_split[:, 0].min())
1.5 1.5
>>> len(m_split) >>> len(m_split)
29 15
>>> m_split.nfaces() >>> m_split.nfaces()
15 15
...@@ -1050,10 +1052,16 @@ class Mesh3D(tfields.TensorMaps): ...@@ -1050,10 +1052,16 @@ class Mesh3D(tfields.TensorMaps):
additionally an instruction to conduct the exact same cut fast (template) additionally an instruction to conduct the exact same cut fast (template)
>>> m_split_2, template = m.cut(cutExpr, at_intersection='split', >>> m_split_2, template = m.cut(cutExpr, at_intersection='split',
... return_template=True) ... return_template=True)
>>> assert m_split.equal(m_split_2) >>> m_split_template = m.cut(template)
>>> assert m_split.equal(m.cut(template)) >>> assert m_split.equal(m_split_2, equal_nan=True)
>>> m_split.fields >>> assert m_split.equal(m_split_template, equal_nan=True)
>>> m.cut(template).fields >>> assert len(template.fields) == 1
>>> assert len(m_split.fields) == 1
>>> assert len(m_split_template.fields) == 1
>>> assert m_split.fields[0].equal(
... list(range(8, 16)) + [np.nan] * 7, equal_nan=True)
>>> assert m_split_template.fields[0].equal(
... list(range(8, 16)) + [np.nan] * 7, equal_nan=True)
This seems irrelevant at first but Consider, the map field or the This seems irrelevant at first but Consider, the map field or the
tensor field changes: tensor field changes:
...@@ -1121,7 +1129,7 @@ class Mesh3D(tfields.TensorMaps): ...@@ -1121,7 +1129,7 @@ class Mesh3D(tfields.TensorMaps):
if __name__ == '__main__': if __name__ == '__main__':
import doctest import doctest
doctest.run_docstring_examples(Mesh3D.cut, globals()) # doctest.run_docstring_examples(Mesh3D.cut, globals())
# doctest.run_docstring_examples(Mesh3D._cut_template, globals()) # doctest.run_docstring_examples(Mesh3D._cut_template, globals())
quit() # quit()
doctest.testmod() doctest.testmod()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment