forked from finsberg/pulse
-
Notifications
You must be signed in to change notification settings - Fork 0
/
test_geometry.py
156 lines (109 loc) · 4.62 KB
/
test_geometry.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
import pytest
import dolfin
import numpy as np
from pulse.example_meshes import mesh_paths
from pulse.geometry import (Geometry, Marker, CRLBasis, HeartGeometry,
Microstructure, MarkerFunctions)
from pulse.dolfin_utils import QuadratureSpace
class Free(dolfin.SubDomain):
def inside(self, x, on_boundary):
return x[0] > (1.0 - dolfin.DOLFIN_EPS) and on_boundary
class Fixed(dolfin.SubDomain):
def inside(self, x, on_boundary):
return x[0] < dolfin.DOLFIN_EPS and on_boundary
fixed = Fixed()
fixed_marker = 1
free = Free()
free_marker = 2
def strain_markers_3d(mesh, nregions):
strain_markers = dolfin.MeshFunction("size_t", mesh, 3)
strain_markers.set_all(0)
xs = np.linspace(0, 1, nregions+1)
region = 0
for it_x in range(nregions):
for it_y in range(nregions):
for it_z in range(nregions):
region += 1
domain_str = ""
domain_str += "x[0] >= {}".format(xs[it_x])
domain_str += " && x[1] >= {}".format(xs[it_y])
domain_str += " && x[2] >= {}".format(xs[it_z])
domain_str += " && x[0] <= {}".format(xs[it_x+1])
domain_str += " && x[1] <= {}".format(xs[it_y+1])
domain_str += " && x[2] <= {}".format(xs[it_z+1])
len_sub = dolfin.CompiledSubDomain(domain_str)
len_sub.mark(strain_markers, region)
return strain_markers
@pytest.fixture
def unitcube_geometry():
N = 2
mesh = dolfin.UnitCubeMesh(N, N, N)
V_f = QuadratureSpace(mesh, 4)
l0 = dolfin.interpolate(
dolfin.Expression(("1.0", "0.0", "0.0"), degree=1), V_f)
r0 = dolfin.interpolate(
dolfin.Expression(("0.0", "1.0", "0.0"), degree=1), V_f)
c0 = dolfin.interpolate(
dolfin.Expression(("0.0", "0.0", "1.0"), degree=1), V_f)
crl_basis = CRLBasis(l0=l0, r0=r0, c0=c0)
cfun = strain_markers_3d(mesh, 2)
ffun = dolfin.MeshFunction("size_t", mesh, 2)
ffun.set_all(0)
fixed.mark(ffun, fixed_marker)
free.mark(ffun, free_marker)
marker_functions = MarkerFunctions(ffun=ffun, cfun=cfun)
fixed_marker_ = Marker(name='fixed', value=fixed_marker, dimension=2)
free_marker_ = Marker(name='free', value=free_marker, dimension=2)
markers = (fixed_marker_, free_marker_)
# Fibers
f0 = dolfin.interpolate(
dolfin.Expression(("1.0", "0.0", "0.0"), degree=1), V_f)
s0 = dolfin.interpolate(
dolfin.Expression(("0.0", "1.0", "0.0"), degree=1), V_f)
n0 = dolfin.interpolate(
dolfin.Expression(("0.0", "0.0", "1.0"), degree=1), V_f)
microstructure = Microstructure(f0=f0, s0=s0, n0=n0)
geometry = Geometry(mesh=mesh, markers=markers,
marker_functions=marker_functions,
microstructure=microstructure,
crl_basis=crl_basis)
return geometry
def test_create_geometry(unitcube_geometry):
# Check that attributes are set correctly
assert unitcube_geometry.f0 is not None
assert unitcube_geometry.s0 is not None
assert unitcube_geometry.n0 is not None
assert unitcube_geometry.microstructure is not None
assert unitcube_geometry.l0 is not None
assert unitcube_geometry.c0 is not None
assert unitcube_geometry.r0 is not None
assert unitcube_geometry.crl_basis is not None
assert unitcube_geometry.ffun is not None
assert unitcube_geometry.cfun is not None
assert unitcube_geometry.mesh is not None
assert unitcube_geometry.markers is not None
def test_mesvolumes(unitcube_geometry):
assert (float(unitcube_geometry.meshvol) - 1.0) < dolfin.DOLFIN_EPS_LARGE
assert all([(float(f) - 0.125) < dolfin.DOLFIN_EPS_LARGE
for f in unitcube_geometry.meshvols])
assert unitcube_geometry.nregions == 8
def test_crl_basis(unitcube_geometry):
assert unitcube_geometry.crl_basis_list
assert unitcube_geometry.nbasis == 3
def test_simple_ellipsoid():
geometry = HeartGeometry.from_file(mesh_paths['simple_ellipsoid'])
assert abs(geometry.cavity_volume() - 2.511304019359619) < 1e-14
assert geometry.is_biv is False
def test_biv_ellipsoid():
geometry = HeartGeometry.from_file(mesh_paths['biv_ellipsoid'])
assert abs(geometry.cavity_volume(chamber='lv')
- 0.7422747965172004) < 1e-14
assert abs(geometry.cavity_volume(chamber='rv')
- 0.8171389194959642) < 1e-14
assert geometry.is_biv is True
if __name__ == "__main__":
# geo = unitcube_geometry()
# test_create_geometry(geo)
# test_mesvolumes(geo)
test_simple_ellipsoid()
test_biv_ellipsoid()