Skip to content

Commit db045d0

Browse files
pletzerpletzer
authored
Merge pull request #46 from pletzer/vector_check
Vector check
2 parents 54ff548 + 16dbba2 commit db045d0

File tree

8 files changed

+345
-108
lines changed

8 files changed

+345
-108
lines changed

‎mint/tests/test_polyline_integral.py‎

Lines changed: 53 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -1,20 +1,32 @@
11
from mint import Grid, PolylineIntegral
22
import numpy
3-
import sys
4-
from pathlib import Path
3+
import pytest
54

65

76
def potentialFunc(p):
87
x, y = p[:2]
98
return x + 2*y
109

1110

12-
def test_simple():
11+
def singularPotentialFunc(p):
12+
x, y = p[:2]
13+
return numpy.arctan2(y, x)/(2.*numpy.pi)
14+
15+
16+
@pytest.mark.parametrize("nx", [3])
17+
@pytest.mark.parametrize("ny", [2])
18+
@pytest.mark.parametrize("potFunc", [potentialFunc, singularPotentialFunc])
19+
@pytest.mark.parametrize("xyz", [numpy.array([(1., 0., 0.),
20+
(0., 1., 0.)]),
21+
numpy.array([(0., 0., 0.),
22+
(1., 0., 0.),
23+
(1., 1., 0.),
24+
(0., 1., 0.)])])
25+
def test_simple(nx, ny, potFunc, xyz):
1326

1427
# create the grid and the edge data
1528
gr = Grid()
1629

17-
nx, ny = 3, 2
1830
points = numpy.zeros((nx*ny, 4, 3), numpy.float64)
1931
data = numpy.zeros((nx*ny, 4))
2032
dx = 1.0 / float(nx)
@@ -46,40 +58,36 @@ def test_simple():
4658
# | |
4759
# +-->--+
4860
# 0
49-
data[k, 0] = potentialFunc(points[k, 1, :]) - potentialFunc(points[k, 0, :])
50-
data[k, 1] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 1, :])
51-
data[k, 2] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 3, :])
52-
data[k, 3] = potentialFunc(points[k, 3, :]) - potentialFunc(points[k, 0, :])
61+
data[k, 0] = potFunc(points[k, 1, :]) - potFunc(points[k, 0, :])
62+
data[k, 1] = potFunc(points[k, 2, :]) - potFunc(points[k, 1, :])
63+
data[k, 2] = potFunc(points[k, 2, :]) - potFunc(points[k, 3, :])
64+
data[k, 3] = potFunc(points[k, 3, :]) - potFunc(points[k, 0, :])
5365

5466
# increment the cell counter
5567
k += 1
5668

5769
gr.setPoints(points)
5870
gr.dump('test_polyline_integral.vtk')
5971

60-
6172
pli = PolylineIntegral()
62-
63-
# create the polyline through which the flux will be integrated
64-
xyz = numpy.array([(0., 0., 0.),
65-
(1., 0., 0.),
66-
(1., 1., 0.),
67-
(0., 1., 0.)])
68-
6973
# no periodicity in x
7074
pli.build(gr, xyz, counterclock=False, periodX=0.0)
7175

7276
flux = pli.getIntegral(data)
73-
exactFlux = potentialFunc(xyz[-1, :]) - potentialFunc(xyz[0, :])
77+
exactFlux = potFunc(xyz[-1, :]) - potFunc(xyz[0, :])
7478
print(f'total flux: {flux:.3f} exact flux: {exactFlux:.3f}')
7579
assert abs(flux - exactFlux) < 1.e-10
7680

77-
def test_partially_outside():
7881

82+
@pytest.mark.parametrize("nx", [3])
83+
@pytest.mark.parametrize("ny", [2])
84+
@pytest.mark.parametrize("potFunc", [potentialFunc])
85+
def test_partially_outside(nx, ny, potFunc):
86+
87+
print('target line is partially outside the domain, expect a warning!')
7988
# create the grid and the edge data
8089
gr = Grid()
8190

82-
nx, ny = 3, 2
8391
points = numpy.zeros((nx*ny, 4, 3), numpy.float64)
8492
data = numpy.zeros((nx*ny, 4))
8593
dx = 1.0 / float(nx)
@@ -111,10 +119,10 @@ def test_partially_outside():
111119
# | |
112120
# +-->--+
113121
# 0
114-
data[k, 0] = potentialFunc(points[k, 1, :]) - potentialFunc(points[k, 0, :])
115-
data[k, 1] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 1, :])
116-
data[k, 2] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 3, :])
117-
data[k, 3] = potentialFunc(points[k, 3, :]) - potentialFunc(points[k, 0, :])
122+
data[k, 0] = potFunc(points[k, 1, :]) - potFunc(points[k, 0, :])
123+
data[k, 1] = potFunc(points[k, 2, :]) - potFunc(points[k, 1, :])
124+
data[k, 2] = potFunc(points[k, 2, :]) - potFunc(points[k, 3, :])
125+
data[k, 3] = potFunc(points[k, 3, :]) - potFunc(points[k, 0, :])
118126

119127
# increment the cell counter
120128
k += 1
@@ -135,19 +143,23 @@ def test_partially_outside():
135143
flux = pli.getIntegral(data)
136144

137145
# because the first point is outside the domain, only the contribution
138-
# stemming from the path inside the domain will be computed. Let's
146+
# stemming from the path inside the domain will be computed. Let's
139147
# correct for this by moving the first point inwards
140148
xyz[0, 0] = 0.
141149
exactFlux = potentialFunc(xyz[-1, :]) - potentialFunc(xyz[0, :])
142150
print(f'total flux: {flux:.3f} exact flux: {exactFlux:.3f}')
143151
assert abs(flux - exactFlux) < 1.e-10
144152

145-
def test_completely_outside():
146153

154+
@pytest.mark.parametrize("nx", [3])
155+
@pytest.mark.parametrize("ny", [2])
156+
@pytest.mark.parametrize("potFunc", [potentialFunc])
157+
def test_completely_outside(nx, ny, potFunc):
158+
159+
print('target line is outside the domain, expect warnings!')
147160
# create the grid and the edge data
148161
gr = Grid()
149162

150-
nx, ny = 3, 2
151163
points = numpy.zeros((nx*ny, 4, 3), numpy.float64)
152164
data = numpy.zeros((nx*ny, 4))
153165
dx = 1.0 / float(nx)
@@ -179,10 +191,10 @@ def test_completely_outside():
179191
# | |
180192
# +-->--+
181193
# 0
182-
data[k, 0] = potentialFunc(points[k, 1, :]) - potentialFunc(points[k, 0, :])
183-
data[k, 1] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 1, :])
184-
data[k, 2] = potentialFunc(points[k, 2, :]) - potentialFunc(points[k, 3, :])
185-
data[k, 3] = potentialFunc(points[k, 3, :]) - potentialFunc(points[k, 0, :])
194+
data[k, 0] = potFunc(points[k, 1, :]) - potFunc(points[k, 0, :])
195+
data[k, 1] = potFunc(points[k, 2, :]) - potFunc(points[k, 1, :])
196+
data[k, 2] = potFunc(points[k, 2, :]) - potFunc(points[k, 3, :])
197+
data[k, 3] = potFunc(points[k, 3, :]) - potFunc(points[k, 0, :])
186198

187199
# increment the cell counter
188200
k += 1
@@ -208,8 +220,16 @@ def test_completely_outside():
208220

209221
if __name__ == '__main__':
210222

211-
test_simple()
212-
test_partially_outside()
213-
test_completely_outside()
223+
# polyline through which the line integral will be computed
224+
xyz = numpy.array([(1., 0., 0.),
225+
(0., 1., 0.)])
226+
test_simple(3, 2, singularPotentialFunc, xyz)
214227

228+
xyz = numpy.array([(0., 0., 0.),
229+
(1., 0., 0.),
230+
(1., 1., 0.),
231+
(0., 1., 0.)])
232+
test_simple(3, 2, potentialFunc, xyz)
215233

234+
test_partially_outside(2, 3, potentialFunc)
235+
test_completely_outside(2, 3, potentialFunc)

‎mint/tests/test_vector_interp.py‎

Lines changed: 22 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -128,11 +128,17 @@ def test_rectilinear():
128128
# set one edge to 1, all other edges to zero
129129
data[cellId, edgeIndex] = 1.0
130130

131-
# get the interpolated vectors
132-
vectorData = vi.getVectors(numpy.array(data))
131+
# get the edge interpolated vectors
132+
vectorData = vi.getEdgeVectors(data)
133133
assert(abs(vectorData.max() - 1.) < 1.e-12)
134134
assert(abs(vectorData.min() - 0.) < 1.e-12)
135135

136+
# get the lateral flux interpolated vectors
137+
vectorData = vi.getFaceVectors(data)
138+
# face vectors take the sign of the area vector, negative if pointing down
139+
assert(abs(numpy.fabs(vectorData).max() - 1.) < 1.e-12)
140+
assert(abs(numpy.fabs(vectorData).min() - 0.) < 1.e-12)
141+
136142
# reset this edge's value back to its original
137143
data[cellId, edgeIndex] = 0.0
138144

@@ -174,15 +180,19 @@ def test_slanted():
174180
# set one edge to 1, all other edges to zero
175181
data[cellId, edgeIndex] = 1.0
176182

177-
# get the interpolated vectors
178-
vectorData = vi.getVectors(numpy.array(data))
183+
# get the edge interpolated vectors
184+
vectorData = vi.getEdgeVectors(data)
185+
fileName = f'slanted_edgeVectors_cellId{cellId}edgeIndex{edgeIndex}.vtk'
186+
saveVectorsVTKFile(targetPoints, vectorData, fileName)
187+
188+
# get the lateral face interpolated vectors
189+
vectorData = vi.getFaceVectors(data)
190+
fileName = f'slanted_faceVectors_cellId{cellId}edgeIndex{edgeIndex}.vtk'
191+
saveVectorsVTKFile(targetPoints, vectorData, fileName)
179192

180193
# reset this edge's value back to its original
181194
data[cellId, edgeIndex] = 0.0
182195

183-
fileName = f'slanted_cellId{cellId}edgeIndex{edgeIndex}.vtk'
184-
saveVectorsVTKFile(targetPoints, vectorData, fileName)
185-
186196

187197
def test_degenerate():
188198

@@ -221,8 +231,11 @@ def test_degenerate():
221231
# set one edge to 1, all other edges to zero
222232
data[cellId, edgeIndex] = 1.0
223233

224-
# get the interpolated vectors
225-
vectorData = vi.getVectors(numpy.array(data))
234+
# get the edge interpolated vectors
235+
vectorData = vi.getEdgeVectors(data)
236+
237+
# get the face interpolated vectors
238+
vectorData = vi.getFaceVectors(data)
226239

227240
# reset this edge's value back to its original
228241
data[cellId, edgeIndex] = 0.0

‎mint/vector_interp.py‎

Lines changed: 35 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -5,14 +5,20 @@
55
import logging
66

77

8+
FILE = 'vectorinterp.py'
9+
DOUBLE_ARRAY_PTR = numpy.ctypeslib.ndpointer(dtype=numpy.float64)
10+
11+
812
def error_handler(filename, methodname, ier):
913
msg = f'ier={ier} after calling {methodname} in {filename}!'
1014
logging.error(msg)
1115
raise RuntimeError(msg)
1216

1317

14-
FILE = 'vectorinterp.py'
15-
DOUBLE_ARRAY_PTR = numpy.ctypeslib.ndpointer(dtype=numpy.float64)
18+
def warning_handler(filename, methodname, ier, detailedMsg=''):
19+
msg = f'ier={ier} after calling {methodname} in {filename}!\n'
20+
msg += detailedMsg
21+
logging.warning(msg)
1622

1723

1824
class VectorInterp(object):
@@ -95,19 +101,38 @@ def findPoints(self, targetPoints, tol2):
95101
targetPoints, tol2)
96102
return numBad
97103

98-
def getVectors(self, data):
104+
def getEdgeVectors(self, data):
99105
"""
100-
Get the vectors at the target points.
106+
Get the edge vectors at given target points.
101107
102108
:param data: edge data array of size numCells times 4
103109
:returns vector array of size numTargetPoints times 3
104-
:note: call this after invoking findPoints
110+
:note: call this after invoking findPoints.
105111
"""
106-
LIB.mnt_vectorinterp_getVectors.argtypes = [POINTER(c_void_p),
107-
DOUBLE_ARRAY_PTR,
108-
DOUBLE_ARRAY_PTR]
112+
LIB.mnt_vectorinterp_getEdgeVectors.argtypes = [POINTER(c_void_p),
113+
DOUBLE_ARRAY_PTR,
114+
DOUBLE_ARRAY_PTR]
115+
res = numpy.zeros((self.numTargetPoints, 3), numpy.float64)
116+
ier = LIB.mnt_vectorinterp_getEdgeVectors(self.obj, data, res)
117+
if ier:
118+
msg = "Some target lines fall outside the grid."
119+
warning_handler(FILE, 'getEdgeVectors', ier, detailedMsg=msg)
120+
return res
121+
122+
def getFaceVectors(self, data):
123+
"""
124+
Get the lateral face vectors at given target points.
125+
126+
:param data: edge data array of size numCells times 4
127+
:returns vector array of size numTargetPoints times 3
128+
:note: call this after invoking findPoints.
129+
"""
130+
LIB.mnt_vectorinterp_getFaceVectors.argtypes = [POINTER(c_void_p),
131+
DOUBLE_ARRAY_PTR,
132+
DOUBLE_ARRAY_PTR]
109133
res = numpy.zeros((self.numTargetPoints, 3), numpy.float64)
110-
ier = LIB.mnt_vectorinterp_getVectors(self.obj, data, res)
134+
ier = LIB.mnt_vectorinterp_getFaceVectors(self.obj, data, res)
111135
if ier:
112-
error_handler(FILE, 'getVectors', ier)
136+
msg = "Some target lines fall outside the grid."
137+
warning_handler(FILE, 'getFaceVectors', ier, detailedMsg=msg)
113138
return res

0 commit comments

Comments
 (0)