[yt-svn] commit/yt: 9 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Mar 20 14:02:54 PDT 2017
9 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/ea796b6687a9/
Changeset: ea796b6687a9
Branch: yt
User: al007
Date: 2017-03-16 20:44:20+00:00
Summary: Preliminary work on QUAD9 support.
Affected #: 7 files
diff -r 95401107ce0793c6b74d10c336f0633aaee990d5 -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 yt/utilities/lib/autogenerated_element_samplers.pxd
--- a/yt/utilities/lib/autogenerated_element_samplers.pxd
+++ b/yt/utilities/lib/autogenerated_element_samplers.pxd
@@ -25,6 +25,19 @@
double* phys_x) nogil
+cdef void Q2Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void Q2Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
cdef void Tet2Function3D(double* fx,
double* x,
double* vertices,
diff -r 95401107ce0793c6b74d10c336f0633aaee990d5 -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -68,6 +68,31 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+cdef void Q2Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = -0.5*(1 + x[0])*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[8] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(1 + x[1])*vertices[12] + (1 + x[0])*(-1 + x[1])*(1 + x[1])*(-1 + x[1])*vertices[16] - phys_x[0] - 0.5*x[0]*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[14] + 0.25*x[0]*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[0] + 0.25*x[0]*(-1 + x[0])*x[1]*(1 + x[1])*vertices[6] - 0.5*x[0]*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[10] + 0.25*x[0]*(1 + x[0])*x[1]*(-1 + x[1])*vertices[2] + 0.25*x[0]*(1 + x[0])*x[1]*(1 + x[1])*vertices[4];
+ fx[1] = -0.5*(1 + x[0])*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[9] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(1 + x[1])*vertices[13] + (1 + x[0])*(-1 + x[1])*(1 + x[1])*(-1 + x[1])*vertices[17] - phys_x[1] - 0.5*x[0]*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[15] + 0.25*x[0]*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[1] + 0.25*x[0]*(-1 + x[0])*x[1]*(1 + x[1])*vertices[7] - 0.5*x[0]*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[11] + 0.25*x[0]*(1 + x[0])*x[1]*(-1 + x[1])*vertices[3] + 0.25*x[0]*(1 + x[0])*x[1]*(1 + x[1])*vertices[5];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Q2Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -0.5*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[14] + 0.25*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[0] - 0.5*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[8] + 0.25*(-1 + x[0])*x[1]*(1 + x[1])*vertices[6] - 0.5*(-1 + x[0])*x[1]*(1 + x[1])*vertices[12] + (-1 + x[1])*(1 + x[1])*(-1 + x[1])*vertices[16] - 0.5*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[10] + 0.25*(1 + x[0])*x[1]*(-1 + x[1])*vertices[2] - 0.5*(1 + x[0])*x[1]*(-1 + x[1])*vertices[8] + 0.25*(1 + x[0])*x[1]*(1 + x[1])*vertices[4] - 0.5*(1 + x[0])*x[1]*(1 + x[1])*vertices[12] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[10] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[14] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[0] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[2] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[4] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[6];
+ scol[0] = -0.5*(-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[8] + 0.25*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[0] - 0.5*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[14] + 0.25*(-1 + x[1])*x[0]*(1 + x[0])*vertices[2] - 0.5*(-1 + x[1])*x[0]*(1 + x[0])*vertices[10] + pow(-1 + x[1], 2)*(1 + x[0])*vertices[16] + (1 + x[1])*(-2 + 2*x[1])*(1 + x[0])*vertices[16] - 0.5*(1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[12] + 0.25*(1 + x[1])*x[0]*(-1 + x[0])*vertices[6] - 0.5*(1 + x[1])*x[0]*(-1 + x[0])*vertices[14] + 0.25*(1 + x[1])*x[0]*(1 + x[0])*vertices[4] - 0.5*(1 + x[1])*x[0]*(1 + x[0])*vertices[10] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[8] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[12] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[0] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[6] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[2] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[4];
+ rcol[1] = -0.5*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[15] + 0.25*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[1] - 0.5*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[9] + 0.25*(-1 + x[0])*x[1]*(1 + x[1])*vertices[7] - 0.5*(-1 + x[0])*x[1]*(1 + x[1])*vertices[13] + (-1 + x[1])*(1 + x[1])*(-1 + x[1])*vertices[17] - 0.5*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[11] + 0.25*(1 + x[0])*x[1]*(-1 + x[1])*vertices[3] - 0.5*(1 + x[0])*x[1]*(-1 + x[1])*vertices[9] + 0.25*(1 + x[0])*x[1]*(1 + x[1])*vertices[5] - 0.5*(1 + x[0])*x[1]*(1 + x[1])*vertices[13] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[11] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[15] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[1] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[3] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[5] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[7];
+ scol[1] = -0.5*(-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[9] + 0.25*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[1] - 0.5*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[15] + 0.25*(-1 + x[1])*x[0]*(1 + x[0])*vertices[3] - 0.5*(-1 + x[1])*x[0]*(1 + x[0])*vertices[11] + pow(-1 + x[1], 2)*(1 + x[0])*vertices[17] + (1 + x[1])*(-2 + 2*x[1])*(1 + x[0])*vertices[17] - 0.5*(1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[13] + 0.25*(1 + x[1])*x[0]*(-1 + x[0])*vertices[7] - 0.5*(1 + x[1])*x[0]*(-1 + x[0])*vertices[15] + 0.25*(1 + x[1])*x[0]*(1 + x[0])*vertices[5] - 0.5*(1 + x[1])*x[0]*(1 + x[0])*vertices[11] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[9] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[13] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[1] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[7] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[3] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[5];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
cdef void Tet2Function3D(double* fx,
double* x,
double* vertices,
diff -r 95401107ce0793c6b74d10c336f0633aaee990d5 -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -237,6 +237,22 @@
cdef int check_inside(self, double* mapped_coord) nogil
+
+cdef class Q2Sampler2D(NonlinearSolveSampler2D):
+
+ cdef void map_real_to_unit(self,
+ double* mapped_x,
+ double* vertices,
+ double* physical_x) nogil
+
+
+ cdef double sample_at_unit_point(self,
+ double* coord,
+ double* vals) nogil
+
+ cdef int check_inside(self, double* mapped_coord) nogil
+
+
cdef class T2Sampler2D(NonlinearSolveSampler2D):
cdef void map_real_to_unit(self,
diff -r 95401107ce0793c6b74d10c336f0633aaee990d5 -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -24,6 +24,8 @@
Q1Jacobian3D, \
Q1Function2D, \
Q1Jacobian2D, \
+ Q2Function2D, \
+ Q2Jacobian2D, \
W1Function3D, \
W1Jacobian3D, \
T2Function2D, \
@@ -788,6 +790,8 @@
else:
for i in range(2):
mapped_x[i] = x[i]
+ with gil:
+ print(mapped_x[i])
cdef class Q1Sampler2D(NonlinearSolveSampler2D):
@@ -830,6 +834,63 @@
return 0
return 1
+cdef class Q2Sampler2D(NonlinearSolveSampler2D):
+
+ '''
+
+ This implements sampling inside a 2D, quadratic, quadrilateral mesh element.
+
+ '''
+
+ def __init__(self):
+ super(Q2Sampler2D, self).__init__()
+ self.num_mapped_coords = 2
+ self.dim = 2
+ self.func = Q2Function2D
+ self.jac = Q2Jacobian2D
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
+ cdef double F, rm, rp, sm, sp
+ cdef double[9] phi
+ cdef double rv = 0
+
+ zet = coord[0]
+ eta = coord[1]
+ zetm = coord[0] - 1.
+ zetp = coord[0] + 1.
+ etam = coord[1] - 1.
+ etap = coord[1] + 1.
+
+ phi[0] = zet * zetm * eta * etam / 4.
+ phi[1] = zet * zetp * eta * etam / 4.
+ phi[2] = zet * zetp * eta * etap / 4.
+ phi[3] = zet * zetm * eta * etap / 4.
+ phi[4] = zetp * zetm * eta * etam / -2.
+ phi[5] = zet * zetp * etap * etam / -2.
+ phi[6] = zetp * zetm * eta * etap / -2.
+ phi[7] = zet * zetm * etap * etam / -2.
+ phi[8] = zetp * zetm * etap * etam
+
+ for i in range(9):
+ rv += vals[i] * phi[i]
+
+ return rv
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_inside(self, double* mapped_coord) nogil:
+ # for quads, we check whether the mapped_coord is between
+ # -1 and 1 in both directions.
+ if (fabs(mapped_coord[0]) - 1.0 > self.inclusion_tol or
+ fabs(mapped_coord[1]) - 1.0 > self.inclusion_tol):
+ return 0
+ return 1
+
+
cdef class T2Sampler2D(NonlinearSolveSampler2D):
'''
@@ -1093,6 +1154,23 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+def test_quad2_sampler(np.ndarray[np.float64_t, ndim=2] vertices,
+ np.ndarray[np.float64_t, ndim=1] field_values,
+ np.ndarray[np.float64_t, ndim=1] physical_x):
+
+ cdef double val
+
+ sampler = Q2Sampler2D()
+
+ val = sampler.sample_at_real_point(<double*> vertices.data,
+ <double*> field_values.data,
+ <double*> physical_x.data)
+
+ return val
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
def test_tri2_sampler(np.ndarray[np.float64_t, ndim=2] vertices,
np.ndarray[np.float64_t, ndim=1] field_values,
np.ndarray[np.float64_t, ndim=1] physical_x):
diff -r 95401107ce0793c6b74d10c336f0633aaee990d5 -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -30,6 +30,7 @@
P1Sampler3D, \
Q1Sampler3D, \
Q1Sampler2D, \
+ Q2Sampler2D, \
S2Sampler3D, \
W1Sampler3D, \
T2Sampler2D, \
@@ -609,6 +610,8 @@
sampler = P1Sampler1D()
elif ndim == 2 and nvertices == 4:
sampler = Q1Sampler2D()
+ elif ndim == 2 and nvertices == 9:
+ sampler = Q2Sampler2D()
elif ndim == 2 and nvertices == 6:
sampler = T2Sampler2D()
elif ndim == 3 and nvertices == 10:
diff -r 95401107ce0793c6b74d10c336f0633aaee990d5 -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 yt/utilities/lib/tests/test_element_mappings.py
--- a/yt/utilities/lib/tests/test_element_mappings.py
+++ b/yt/utilities/lib/tests/test_element_mappings.py
@@ -21,6 +21,7 @@
test_hex_sampler, \
test_tri_sampler, \
test_quad_sampler, \
+ test_quad2_sampler, \
test_hex20_sampler, \
test_wedge_sampler, \
test_tri2_sampler, \
@@ -32,10 +33,12 @@
NV = vertices.shape[0]
NDIM = vertices.shape[1]
x = np.empty(NDIM)
+ # import pdb; pdb.set_trace()
for i in range(NV):
x = vertices[i]
val = sampler(vertices, field_values, x)
- assert_almost_equal(val, field_values[i])
+ print("Value is %s, field_value is %s" % (val, field_values[i]))
+ # assert_almost_equal(val, field_values[i])
def test_P1Sampler1D():
@@ -77,6 +80,24 @@
check_all_vertices(test_quad_sampler, vertices, field_values)
+def test_Q2Sampler2D():
+
+ vertices = np.array([[2., 3.],
+ [7., 4.],
+ [10., 15.],
+ [4., 12.],
+ [3., 3.],
+ [9., 8.],
+ [6., 14.],
+ [3., 6.],
+ [4., 8.]])
+
+ field_values = np.array([7., 27., 40., 12., 13., 30., 22., 9., 16.])
+
+
+ check_all_vertices(test_quad2_sampler, vertices, field_values)
+
+
def test_Q1Sampler3D():
vertices = np.array([[2.00657905, 0.6888599, 1.4375],
[1.8658198, 1.00973171, 1.4375],
diff -r 95401107ce0793c6b74d10c336f0633aaee990d5 -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 yt/utilities/mesh_types.yaml
--- a/yt/utilities/mesh_types.yaml
+++ b/yt/utilities/mesh_types.yaml
@@ -24,6 +24,22 @@
(1 + x[0])*(1 + x[1])/4.,
(1 - x[0])*(1 + x[1])/4.]
+Quad9:
+ mesh_type: Q2
+ num_dim: 2
+ num_vertices: 9
+ num_mapped_coords: 2
+ shape_functions: |
+ [x[0] * (x[0] - 1) * x[1] * (x[1] - 1) / 4.,
+ x[0] * (x[0] + 1) * x[1] * (x[1] - 1) / 4.,
+ x[0] * (x[0] + 1) * x[1] * (x[1] + 1) / 4.,
+ x[0] * (x[0] - 1) * x[1] * (x[1] + 1) / 4.,
+ (x[0] + 1) * (x[0] - 1) * x[1] * (x[1] - 1) / -2.,
+ x[0] * (x[0] + 1) * (x[1] + 1) * (x[1] - 1) / -2.,
+ (x[0] + 1) * (x[0] - 1) * x[1] * (x[1] + 1) / -2.,
+ x[0] * (x[0] - 1) * (x[1] + 1) * (x[1] - 1) / -2.,
+ (x[0] + 1) * (x[1] - 1) * (x[1] + 1) * (x[1] - 1)]
+
Wedge6:
mesh_type: W1
num_dim: 3
https://bitbucket.org/yt_analysis/yt/commits/701f9eaeac29/
Changeset: 701f9eaeac29
Branch: yt
User: al007
Date: 2017-03-17 17:18:37+00:00
Summary: Element mapping now working.
Affected #: 4 files
diff -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 -r 701f9eaeac29503f04cc4122e37f6b8bfd370a04 yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -72,8 +72,8 @@
double* x,
double* vertices,
double* phys_x) nogil:
- fx[0] = -0.5*(1 + x[0])*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[8] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(1 + x[1])*vertices[12] + (1 + x[0])*(-1 + x[1])*(1 + x[1])*(-1 + x[1])*vertices[16] - phys_x[0] - 0.5*x[0]*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[14] + 0.25*x[0]*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[0] + 0.25*x[0]*(-1 + x[0])*x[1]*(1 + x[1])*vertices[6] - 0.5*x[0]*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[10] + 0.25*x[0]*(1 + x[0])*x[1]*(-1 + x[1])*vertices[2] + 0.25*x[0]*(1 + x[0])*x[1]*(1 + x[1])*vertices[4];
- fx[1] = -0.5*(1 + x[0])*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[9] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(1 + x[1])*vertices[13] + (1 + x[0])*(-1 + x[1])*(1 + x[1])*(-1 + x[1])*vertices[17] - phys_x[1] - 0.5*x[0]*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[15] + 0.25*x[0]*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[1] + 0.25*x[0]*(-1 + x[0])*x[1]*(1 + x[1])*vertices[7] - 0.5*x[0]*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[11] + 0.25*x[0]*(1 + x[0])*x[1]*(-1 + x[1])*vertices[3] + 0.25*x[0]*(1 + x[0])*x[1]*(1 + x[1])*vertices[5];
+ fx[0] = (1 + x[0])*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[16] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[8] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(1 + x[1])*vertices[12] - phys_x[0] - 0.5*x[0]*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[14] + 0.25*x[0]*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[0] + 0.25*x[0]*(-1 + x[0])*x[1]*(1 + x[1])*vertices[6] - 0.5*x[0]*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[10] + 0.25*x[0]*(1 + x[0])*x[1]*(-1 + x[1])*vertices[2] + 0.25*x[0]*(1 + x[0])*x[1]*(1 + x[1])*vertices[4];
+ fx[1] = (1 + x[0])*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[17] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[9] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(1 + x[1])*vertices[13] - phys_x[1] - 0.5*x[0]*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[15] + 0.25*x[0]*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[1] + 0.25*x[0]*(-1 + x[0])*x[1]*(1 + x[1])*vertices[7] - 0.5*x[0]*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[11] + 0.25*x[0]*(1 + x[0])*x[1]*(-1 + x[1])*vertices[3] + 0.25*x[0]*(1 + x[0])*x[1]*(1 + x[1])*vertices[5];
@cython.boundscheck(False)
@@ -84,10 +84,10 @@
double* x,
double* vertices,
double* phys_x) nogil:
- rcol[0] = -0.5*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[14] + 0.25*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[0] - 0.5*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[8] + 0.25*(-1 + x[0])*x[1]*(1 + x[1])*vertices[6] - 0.5*(-1 + x[0])*x[1]*(1 + x[1])*vertices[12] + (-1 + x[1])*(1 + x[1])*(-1 + x[1])*vertices[16] - 0.5*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[10] + 0.25*(1 + x[0])*x[1]*(-1 + x[1])*vertices[2] - 0.5*(1 + x[0])*x[1]*(-1 + x[1])*vertices[8] + 0.25*(1 + x[0])*x[1]*(1 + x[1])*vertices[4] - 0.5*(1 + x[0])*x[1]*(1 + x[1])*vertices[12] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[10] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[14] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[0] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[2] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[4] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[6];
- scol[0] = -0.5*(-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[8] + 0.25*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[0] - 0.5*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[14] + 0.25*(-1 + x[1])*x[0]*(1 + x[0])*vertices[2] - 0.5*(-1 + x[1])*x[0]*(1 + x[0])*vertices[10] + pow(-1 + x[1], 2)*(1 + x[0])*vertices[16] + (1 + x[1])*(-2 + 2*x[1])*(1 + x[0])*vertices[16] - 0.5*(1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[12] + 0.25*(1 + x[1])*x[0]*(-1 + x[0])*vertices[6] - 0.5*(1 + x[1])*x[0]*(-1 + x[0])*vertices[14] + 0.25*(1 + x[1])*x[0]*(1 + x[0])*vertices[4] - 0.5*(1 + x[1])*x[0]*(1 + x[0])*vertices[10] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[8] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[12] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[0] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[6] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[2] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[4];
- rcol[1] = -0.5*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[15] + 0.25*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[1] - 0.5*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[9] + 0.25*(-1 + x[0])*x[1]*(1 + x[1])*vertices[7] - 0.5*(-1 + x[0])*x[1]*(1 + x[1])*vertices[13] + (-1 + x[1])*(1 + x[1])*(-1 + x[1])*vertices[17] - 0.5*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[11] + 0.25*(1 + x[0])*x[1]*(-1 + x[1])*vertices[3] - 0.5*(1 + x[0])*x[1]*(-1 + x[1])*vertices[9] + 0.25*(1 + x[0])*x[1]*(1 + x[1])*vertices[5] - 0.5*(1 + x[0])*x[1]*(1 + x[1])*vertices[13] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[11] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[15] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[1] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[3] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[5] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[7];
- scol[1] = -0.5*(-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[9] + 0.25*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[1] - 0.5*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[15] + 0.25*(-1 + x[1])*x[0]*(1 + x[0])*vertices[3] - 0.5*(-1 + x[1])*x[0]*(1 + x[0])*vertices[11] + pow(-1 + x[1], 2)*(1 + x[0])*vertices[17] + (1 + x[1])*(-2 + 2*x[1])*(1 + x[0])*vertices[17] - 0.5*(1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[13] + 0.25*(1 + x[1])*x[0]*(-1 + x[0])*vertices[7] - 0.5*(1 + x[1])*x[0]*(-1 + x[0])*vertices[15] + 0.25*(1 + x[1])*x[0]*(1 + x[0])*vertices[5] - 0.5*(1 + x[1])*x[0]*(1 + x[0])*vertices[11] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[9] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[13] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[1] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[7] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[3] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[5];
+ rcol[0] = -0.5*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[14] + (-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[16] + 0.25*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[0] - 0.5*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[8] + 0.25*(-1 + x[0])*x[1]*(1 + x[1])*vertices[6] - 0.5*(-1 + x[0])*x[1]*(1 + x[1])*vertices[12] - 0.5*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[10] + (1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[16] + 0.25*(1 + x[0])*x[1]*(-1 + x[1])*vertices[2] - 0.5*(1 + x[0])*x[1]*(-1 + x[1])*vertices[8] + 0.25*(1 + x[0])*x[1]*(1 + x[1])*vertices[4] - 0.5*(1 + x[0])*x[1]*(1 + x[1])*vertices[12] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[10] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[14] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[0] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[2] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[4] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[6];
+ scol[0] = -0.5*(-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[8] + (-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[16] + 0.25*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[0] - 0.5*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[14] + 0.25*(-1 + x[1])*x[0]*(1 + x[0])*vertices[2] - 0.5*(-1 + x[1])*x[0]*(1 + x[0])*vertices[10] - 0.5*(1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[12] + (1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[16] + 0.25*(1 + x[1])*x[0]*(-1 + x[0])*vertices[6] - 0.5*(1 + x[1])*x[0]*(-1 + x[0])*vertices[14] + 0.25*(1 + x[1])*x[0]*(1 + x[0])*vertices[4] - 0.5*(1 + x[1])*x[0]*(1 + x[0])*vertices[10] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[8] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[12] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[0] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[6] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[2] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[4];
+ rcol[1] = -0.5*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[15] + (-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[17] + 0.25*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[1] - 0.5*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[9] + 0.25*(-1 + x[0])*x[1]*(1 + x[1])*vertices[7] - 0.5*(-1 + x[0])*x[1]*(1 + x[1])*vertices[13] - 0.5*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[11] + (1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[17] + 0.25*(1 + x[0])*x[1]*(-1 + x[1])*vertices[3] - 0.5*(1 + x[0])*x[1]*(-1 + x[1])*vertices[9] + 0.25*(1 + x[0])*x[1]*(1 + x[1])*vertices[5] - 0.5*(1 + x[0])*x[1]*(1 + x[1])*vertices[13] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[11] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[15] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[1] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[3] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[5] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[7];
+ scol[1] = -0.5*(-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[9] + (-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[17] + 0.25*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[1] - 0.5*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[15] + 0.25*(-1 + x[1])*x[0]*(1 + x[0])*vertices[3] - 0.5*(-1 + x[1])*x[0]*(1 + x[0])*vertices[11] - 0.5*(1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[13] + (1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[17] + 0.25*(1 + x[1])*x[0]*(-1 + x[0])*vertices[7] - 0.5*(1 + x[1])*x[0]*(-1 + x[0])*vertices[15] + 0.25*(1 + x[1])*x[0]*(1 + x[0])*vertices[5] - 0.5*(1 + x[1])*x[0]*(1 + x[0])*vertices[11] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[9] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[13] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[1] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[7] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[3] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[5];
@cython.boundscheck(False)
diff -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 -r 701f9eaeac29503f04cc4122e37f6b8bfd370a04 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -746,7 +746,7 @@
def __init__(self):
super(NonlinearSolveSampler2D, self).__init__()
self.tolerance = 1.0e-9
- self.max_iter = 10
+ self.max_iter = 100
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -756,12 +756,13 @@
double* vertices,
double* physical_x) nogil:
cdef int i
- cdef double d
+ cdef double d, lam
cdef double[2] f
- cdef double[2] x
+ cdef double[2] x, xk, s
cdef double[4] A
cdef int iterations = 0
- cdef double err
+ cdef double err_c, err_plus
+ cdef double alpha = 1e-4
# initial guess
for i in range(2):
@@ -769,29 +770,41 @@
# initial error norm
self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 2)
+ err_c = maxnorm(f, 2)
# begin Newton iteration
- while (err > self.tolerance and iterations < self.max_iter):
+ while (err_c > self.tolerance and iterations < self.max_iter):
self.jac(&A[0], &A[2], x, vertices, physical_x)
d = (A[0]*A[3] - A[1]*A[2])
- x[0] -= ( A[3]*f[0] - A[2]*f[1]) / d
- x[1] -= (-A[1]*f[0] + A[0]*f[1]) / d
+ s[0] = -( A[3]*f[0] - A[2]*f[1]) / d
+ s[1] = -(-A[1]*f[0] + A[0]*f[1]) / d
+ xk[0] = x[0] + s[0]
+ xk[1] = x[1] + s[1]
- self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 2)
+ self.func(f, xk, vertices, physical_x)
+ err_plus = maxnorm(f, 2)
+ lam = 1
+ while err_plus > err_c * (1. - alpha * lam) and lam > 1e-6:
+ lam = lam / 2
+ xk[0] = x[0] + lam * s[0]
+ xk[1] = x[1] + lam * s[1]
+ self.func(f, xk, vertices, physical_x)
+ err_plus = maxnorm(f, 2)
+
+ x[0] = xk[0]
+ x[1] = xk[1]
+ err_c = err_plus
+
iterations += 1
- if (err > self.tolerance):
+ if (err_c > self.tolerance):
# we did not converge, set bogus value
for i in range(2):
mapped_x[i] = -99.0
else:
for i in range(2):
mapped_x[i] = x[i]
- with gil:
- print(mapped_x[i])
cdef class Q1Sampler2D(NonlinearSolveSampler2D):
diff -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 -r 701f9eaeac29503f04cc4122e37f6b8bfd370a04 yt/utilities/lib/tests/test_element_mappings.py
--- a/yt/utilities/lib/tests/test_element_mappings.py
+++ b/yt/utilities/lib/tests/test_element_mappings.py
@@ -37,8 +37,6 @@
for i in range(NV):
x = vertices[i]
val = sampler(vertices, field_values, x)
- print("Value is %s, field_value is %s" % (val, field_values[i]))
- # assert_almost_equal(val, field_values[i])
def test_P1Sampler1D():
@@ -94,7 +92,6 @@
field_values = np.array([7., 27., 40., 12., 13., 30., 22., 9., 16.])
-
check_all_vertices(test_quad2_sampler, vertices, field_values)
diff -r ea796b6687a9d6b9fccd30a77b5d71ba9a542ca1 -r 701f9eaeac29503f04cc4122e37f6b8bfd370a04 yt/utilities/mesh_types.yaml
--- a/yt/utilities/mesh_types.yaml
+++ b/yt/utilities/mesh_types.yaml
@@ -38,7 +38,7 @@
x[0] * (x[0] + 1) * (x[1] + 1) * (x[1] - 1) / -2.,
(x[0] + 1) * (x[0] - 1) * x[1] * (x[1] + 1) / -2.,
x[0] * (x[0] - 1) * (x[1] + 1) * (x[1] - 1) / -2.,
- (x[0] + 1) * (x[1] - 1) * (x[1] + 1) * (x[1] - 1)]
+ (x[0] + 1) * (x[0] - 1) * (x[1] + 1) * (x[1] - 1)]
Wedge6:
mesh_type: W1
https://bitbucket.org/yt_analysis/yt/commits/a9c997840e0c/
Changeset: a9c997840e0c
Branch: yt
User: al007
Date: 2017-03-17 17:42:10+00:00
Summary: Add documentation of Newton routine. Add line search to 3D sampler.
Affected #: 1 file
diff -r 701f9eaeac29503f04cc4122e37f6b8bfd370a04 -r a9c997840e0c95dae72547fd2bfaea3be34df5ed yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -342,15 +342,31 @@
double* mapped_x,
double* vertices,
double* physical_x) nogil:
+ '''
+
+ x: solution vector; holds unit/mapped coordinates
+ xk: temporary vector for holding solution of current iteration
+ f: residual vector
+ r, s, t: three columns of Jacobian matrix corresponding to unit/mapped
+ coordinates r, s, and t
+ d: Jacobian determinant
+ s_n: Newton step vector
+ lam: fraction of Newton step by which to change x
+ alpha: constant proportional to how much residual required to decrease
+ err_c: Error of current iteration
+ err_plus: Error of next iteration
+
+ '''
cdef int i
- cdef double d
+ cdef double d, lam
cdef double[3] f
cdef double[3] r
cdef double[3] s
cdef double[3] t
- cdef double[3] x
+ cdef double[3] x, xk, s_n
cdef int iterations = 0
- cdef double err
+ cdef double err_c, err_plus
+ cdef double alpha = 1e-4
# initial guess
for i in range(3):
@@ -358,20 +374,38 @@
# initial error norm
self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 3)
+ err_c = maxnorm(f, 3)
# begin Newton iteration
- while (err > self.tolerance and iterations < self.max_iter):
+ while (err_c > self.tolerance and iterations < self.max_iter):
self.jac(r, s, t, x, vertices, physical_x)
d = determinant_3x3(r, s, t)
- x[0] = x[0] - (determinant_3x3(f, s, t)/d)
- x[1] = x[1] - (determinant_3x3(r, f, t)/d)
- x[2] = x[2] - (determinant_3x3(r, s, f)/d)
- self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 3)
+
+ s_n[0] = - (determinant_3x3(f, s, t)/d)
+ s_n[1] = - (determinant_3x3(r, f, t)/d)
+ s_n[2] = - (determinant_3x3(r, s, f)/d)
+ xk[0] = x[0] + s_n[0]
+ xk[1] = x[1] + s_n[1]
+ xk[2] = x[2] + s_n[2]
+ self.func(f, xk, vertices, physical_x)
+ err_plus = maxnorm(f, 3)
+
+ lam = 1
+ while err_plus > err_c * (1. - alpha * lam) and lam > 1e-6:
+ lam = lam / 2
+ xk[0] = x[0] + lam * s_n[0]
+ xk[1] = x[1] + lam * s_n[1]
+ xk[2] = x[2] + lam * s_n[2]
+ self.func(f, xk, vertices, physical_x)
+ err_plus = maxnorm(f, 3)
+
+ x[0] = xk[0]
+ x[1] = xk[1]
+ x[2] = xk[2]
+ err_c = err_plus
iterations += 1
- if (err > self.tolerance):
+ if (err_c > self.tolerance):
# we did not converge, set bogus value
for i in range(3):
mapped_x[i] = -99.0
@@ -746,7 +780,7 @@
def __init__(self):
super(NonlinearSolveSampler2D, self).__init__()
self.tolerance = 1.0e-9
- self.max_iter = 100
+ self.max_iter = 10
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -755,10 +789,24 @@
double* mapped_x,
double* vertices,
double* physical_x) nogil:
+ '''
+
+ x: solution vector; holds unit/mapped coordinates
+ xk: temporary vector for holding solution of current iteration
+ f: residual vector
+ A: Jacobian matrix (derivative of residual vector wrt x)
+ d: Jacobian determinant
+ s_n: Newton step vector
+ lam: fraction of Newton step by which to change x
+ alpha: constant proportional to how much residual required to decrease
+ err_c: Error of current iteration
+ err_plus: Error of next iteration
+
+ '''
cdef int i
cdef double d, lam
cdef double[2] f
- cdef double[2] x, xk, s
+ cdef double[2] x, xk, s_n
cdef double[4] A
cdef int iterations = 0
cdef double err_c, err_plus
@@ -777,25 +825,24 @@
self.jac(&A[0], &A[2], x, vertices, physical_x)
d = (A[0]*A[3] - A[1]*A[2])
- s[0] = -( A[3]*f[0] - A[2]*f[1]) / d
- s[1] = -(-A[1]*f[0] + A[0]*f[1]) / d
- xk[0] = x[0] + s[0]
- xk[1] = x[1] + s[1]
-
+ s_n[0] = -( A[3]*f[0] - A[2]*f[1]) / d
+ s_n[1] = -(-A[1]*f[0] + A[0]*f[1]) / d
+ xk[0] = x[0] + s_n[0]
+ xk[1] = x[1] + s_n[1]
self.func(f, xk, vertices, physical_x)
err_plus = maxnorm(f, 2)
+
lam = 1
while err_plus > err_c * (1. - alpha * lam) and lam > 1e-6:
lam = lam / 2
- xk[0] = x[0] + lam * s[0]
- xk[1] = x[1] + lam * s[1]
+ xk[0] = x[0] + lam * s_n[0]
+ xk[1] = x[1] + lam * s_n[1]
self.func(f, xk, vertices, physical_x)
err_plus = maxnorm(f, 2)
x[0] = xk[0]
x[1] = xk[1]
err_c = err_plus
-
iterations += 1
if (err_c > self.tolerance):
https://bitbucket.org/yt_analysis/yt/commits/727c3552f8f0/
Changeset: 727c3552f8f0
Branch: yt
User: al007
Date: 2017-03-17 18:11:34+00:00
Summary: Fix tests.
Affected #: 1 file
diff -r a9c997840e0c95dae72547fd2bfaea3be34df5ed -r 727c3552f8f08c541241e4e308407d66aaac92dd yt/utilities/lib/tests/test_element_mappings.py
--- a/yt/utilities/lib/tests/test_element_mappings.py
+++ b/yt/utilities/lib/tests/test_element_mappings.py
@@ -33,10 +33,10 @@
NV = vertices.shape[0]
NDIM = vertices.shape[1]
x = np.empty(NDIM)
- # import pdb; pdb.set_trace()
for i in range(NV):
x = vertices[i]
val = sampler(vertices, field_values, x)
+ assert_almost_equal(val, field_values[i])
def test_P1Sampler1D():
@@ -84,11 +84,11 @@
[7., 4.],
[10., 15.],
[4., 12.],
- [3., 3.],
- [9., 8.],
- [6., 14.],
- [3., 6.],
- [4., 8.]])
+ [4.5, 3.5],
+ [8.5, 9.5],
+ [7., 13.5],
+ [3., 7.5],
+ [5.75, 8.5]])
field_values = np.array([7., 27., 40., 12., 13., 30., 22., 9., 16.])
https://bitbucket.org/yt_analysis/yt/commits/002cbe00d500/
Changeset: 002cbe00d500
Branch: yt
User: al007
Date: 2017-03-17 20:41:45+00:00
Summary: Make answer test actually show something meaningful.
Affected #: 1 file
diff -r 727c3552f8f08c541241e4e308407d66aaac92dd -r 002cbe00d50029a91c201f1233c1feeaebabe363 yt/visualization/volume_rendering/tests/test_mesh_render.py
--- a/yt/visualization/volume_rendering/tests/test_mesh_render.py
+++ b/yt/visualization/volume_rendering/tests/test_mesh_render.py
@@ -164,6 +164,8 @@
ytcfg["yt", "ray_tracing_engine"] = engine
ds = data_dir_load(tet10, kwargs={'step':-1})
sc = create_scene(ds, field)
+ ms = sc.get_source(0)
+ ms.color_bounds = (-.01, .2)
im = sc.render()
return compare(ds, im, "%s_render_answers_tet10_%s_%s" %
(engine, field[0], field[1]))
https://bitbucket.org/yt_analysis/yt/commits/57f6ebc27d0b/
Changeset: 57f6ebc27d0b
Branch: yt
User: al007
Date: 2017-03-20 14:11:43+00:00
Summary: Add answer test for new quad2 sampler.
Affected #: 2 files
diff -r 002cbe00d50029a91c201f1233c1feeaebabe363 -r 57f6ebc27d0b05318e077a26ddb6168f67c305ba tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -63,9 +63,10 @@
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
- yt/analysis_modules/photon_simulator/tests/test_sloshing.py
- local_unstructured_002:
+ local_unstructured_003:
- yt/visualization/volume_rendering/tests/test_mesh_render.py
- yt/visualization/tests/test_mesh_slices.py:test_tri2
+ - yt/visualization/tests/test_mesh_slices.py:test_quad2
- yt/visualization/tests/test_mesh_slices.py:test_multi_region
local_boxlib_003:
diff -r 002cbe00d50029a91c201f1233c1feeaebabe363 -r 57f6ebc27d0b05318e077a26ddb6168f67c305ba yt/visualization/tests/test_mesh_slices.py
--- a/yt/visualization/tests/test_mesh_slices.py
+++ b/yt/visualization/tests/test_mesh_slices.py
@@ -49,6 +49,15 @@
for field in ds.field_list:
yield compare(ds, field, "answers_tri2_%s_%s" % (field[0], field[1]))
+quad2 = "SecondOrderQuads/lid_driven_out.e"
+
+ at requires_ds(quad2)
+def test_quad2():
+ ds = data_dir_load(quad2, kwargs={'step':-1})
+ field_list = [('all', 'T'), ('all', 'vel_x'), ('all', 'vel_y')]
+ for field in field_list:
+ yield compare(ds, field, "answers_quad2_%s_%s" % (field[0], field[1]))
+
multi_region = "MultiRegion/two_region_example_out.e"
@requires_ds(multi_region)
https://bitbucket.org/yt_analysis/yt/commits/7b25ca463d3e/
Changeset: 7b25ca463d3e
Branch: yt
User: al007
Date: 2017-03-20 18:01:20+00:00
Summary: Introduce min_lam variable and add a little more doco.
Affected #: 1 file
diff -r 57f6ebc27d0b05318e077a26ddb6168f67c305ba -r 7b25ca463d3e96d24c67211f7b0e7a2bfd511c7e yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -344,6 +344,10 @@
double* physical_x) nogil:
'''
+ A thorough description of Newton's method and modifications for global
+ convergence can be found in Dennis's text "Numerical Methods for
+ Unconstrained Optimization and Nonlinear Equations."
+
x: solution vector; holds unit/mapped coordinates
xk: temporary vector for holding solution of current iteration
f: residual vector
@@ -352,11 +356,16 @@
d: Jacobian determinant
s_n: Newton step vector
lam: fraction of Newton step by which to change x
- alpha: constant proportional to how much residual required to decrease
+ alpha: constant proportional to how much residual required to decrease.
+ 1e-4 is value of alpha recommended by Dennis
err_c: Error of current iteration
err_plus: Error of next iteration
+ min_lam: minimum fraction of Newton step that the line search is allowed
+ to take. General experience suggests that lambda values smaller
+ than 1e-3 will not significantly reduce the residual, but we
+ set to 1e-6 just to be safe
+ '''
- '''
cdef int i
cdef double d, lam
cdef double[3] f
@@ -367,6 +376,7 @@
cdef int iterations = 0
cdef double err_c, err_plus
cdef double alpha = 1e-4
+ cdef double min_lam = 1e-6
# initial guess
for i in range(3):
@@ -391,7 +401,7 @@
err_plus = maxnorm(f, 3)
lam = 1
- while err_plus > err_c * (1. - alpha * lam) and lam > 1e-6:
+ while err_plus > err_c * (1. - alpha * lam) and lam > min_lam:
lam = lam / 2
xk[0] = x[0] + lam * s_n[0]
xk[1] = x[1] + lam * s_n[1]
@@ -791,6 +801,10 @@
double* physical_x) nogil:
'''
+ A thorough description of Newton's method and modifications for global
+ convergence can be found in Dennis's text "Numerical Methods for
+ Unconstrained Optimization and Nonlinear Equations."
+
x: solution vector; holds unit/mapped coordinates
xk: temporary vector for holding solution of current iteration
f: residual vector
@@ -798,11 +812,16 @@
d: Jacobian determinant
s_n: Newton step vector
lam: fraction of Newton step by which to change x
- alpha: constant proportional to how much residual required to decrease
+ alpha: constant proportional to how much residual required to decrease.
+ 1e-4 is value of alpha recommended by Dennis
err_c: Error of current iteration
err_plus: Error of next iteration
+ min_lam: minimum fraction of Newton step that the line search is allowed
+ to take. General experience suggests that lambda values smaller
+ than 1e-3 will not significantly reduce the residual, but we
+ set to 1e-6 just to be safe
+ '''
- '''
cdef int i
cdef double d, lam
cdef double[2] f
@@ -811,6 +830,7 @@
cdef int iterations = 0
cdef double err_c, err_plus
cdef double alpha = 1e-4
+ cdef double min_lam = 1e-6
# initial guess
for i in range(2):
@@ -833,7 +853,7 @@
err_plus = maxnorm(f, 2)
lam = 1
- while err_plus > err_c * (1. - alpha * lam) and lam > 1e-6:
+ while err_plus > err_c * (1. - alpha * lam) and lam > min_lam:
lam = lam / 2
xk[0] = x[0] + lam * s_n[0]
xk[1] = x[1] + lam * s_n[1]
https://bitbucket.org/yt_analysis/yt/commits/f085d71d2e84/
Changeset: f085d71d2e84
Branch: yt
User: al007
Date: 2017-03-20 19:52:49+00:00
Summary: Increment answer test number again.
Affected #: 1 file
diff -r 7b25ca463d3e96d24c67211f7b0e7a2bfd511c7e -r f085d71d2e84f29e23711c5870166321ccde9adb tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -63,7 +63,7 @@
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
- yt/analysis_modules/photon_simulator/tests/test_sloshing.py
- local_unstructured_003:
+ local_unstructured_004:
- yt/visualization/volume_rendering/tests/test_mesh_render.py
- yt/visualization/tests/test_mesh_slices.py:test_tri2
- yt/visualization/tests/test_mesh_slices.py:test_quad2
https://bitbucket.org/yt_analysis/yt/commits/eaf22dd8f960/
Changeset: eaf22dd8f960
Branch: yt
User: ngoldbaum
Date: 2017-03-20 21:02:44+00:00
Summary: Merged in al007/yt (pull request #2549)
Add visualization support for QUAD9 elements
Approved-by: Andrew Myers <atmyers2 at gmail.com>
Approved-by: Nathan Goldbaum <ngoldbau at illinois.edu>
Approved-by: yt-fido <yt.fido at gmail.com>
Affected #: 10 files
diff -r 398bded61777ec01e1bdf868edb966b7e52ece57 -r eaf22dd8f960f1767b66025127fb4413d0654853 tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -63,9 +63,10 @@
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
- yt/analysis_modules/photon_simulator/tests/test_sloshing.py
- local_unstructured_002:
+ local_unstructured_004:
- yt/visualization/volume_rendering/tests/test_mesh_render.py
- yt/visualization/tests/test_mesh_slices.py:test_tri2
+ - yt/visualization/tests/test_mesh_slices.py:test_quad2
- yt/visualization/tests/test_mesh_slices.py:test_multi_region
local_boxlib_003:
diff -r 398bded61777ec01e1bdf868edb966b7e52ece57 -r eaf22dd8f960f1767b66025127fb4413d0654853 yt/utilities/lib/autogenerated_element_samplers.pxd
--- a/yt/utilities/lib/autogenerated_element_samplers.pxd
+++ b/yt/utilities/lib/autogenerated_element_samplers.pxd
@@ -25,6 +25,19 @@
double* phys_x) nogil
+cdef void Q2Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void Q2Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
cdef void Tet2Function3D(double* fx,
double* x,
double* vertices,
diff -r 398bded61777ec01e1bdf868edb966b7e52ece57 -r eaf22dd8f960f1767b66025127fb4413d0654853 yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -68,6 +68,31 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+cdef void Q2Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = (1 + x[0])*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[16] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[8] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(1 + x[1])*vertices[12] - phys_x[0] - 0.5*x[0]*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[14] + 0.25*x[0]*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[0] + 0.25*x[0]*(-1 + x[0])*x[1]*(1 + x[1])*vertices[6] - 0.5*x[0]*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[10] + 0.25*x[0]*(1 + x[0])*x[1]*(-1 + x[1])*vertices[2] + 0.25*x[0]*(1 + x[0])*x[1]*(1 + x[1])*vertices[4];
+ fx[1] = (1 + x[0])*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[17] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[9] - 0.5*(1 + x[0])*(-1 + x[0])*x[1]*(1 + x[1])*vertices[13] - phys_x[1] - 0.5*x[0]*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[15] + 0.25*x[0]*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[1] + 0.25*x[0]*(-1 + x[0])*x[1]*(1 + x[1])*vertices[7] - 0.5*x[0]*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[11] + 0.25*x[0]*(1 + x[0])*x[1]*(-1 + x[1])*vertices[3] + 0.25*x[0]*(1 + x[0])*x[1]*(1 + x[1])*vertices[5];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Q2Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -0.5*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[14] + (-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[16] + 0.25*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[0] - 0.5*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[8] + 0.25*(-1 + x[0])*x[1]*(1 + x[1])*vertices[6] - 0.5*(-1 + x[0])*x[1]*(1 + x[1])*vertices[12] - 0.5*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[10] + (1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[16] + 0.25*(1 + x[0])*x[1]*(-1 + x[1])*vertices[2] - 0.5*(1 + x[0])*x[1]*(-1 + x[1])*vertices[8] + 0.25*(1 + x[0])*x[1]*(1 + x[1])*vertices[4] - 0.5*(1 + x[0])*x[1]*(1 + x[1])*vertices[12] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[10] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[14] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[0] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[2] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[4] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[6];
+ scol[0] = -0.5*(-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[8] + (-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[16] + 0.25*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[0] - 0.5*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[14] + 0.25*(-1 + x[1])*x[0]*(1 + x[0])*vertices[2] - 0.5*(-1 + x[1])*x[0]*(1 + x[0])*vertices[10] - 0.5*(1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[12] + (1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[16] + 0.25*(1 + x[1])*x[0]*(-1 + x[0])*vertices[6] - 0.5*(1 + x[1])*x[0]*(-1 + x[0])*vertices[14] + 0.25*(1 + x[1])*x[0]*(1 + x[0])*vertices[4] - 0.5*(1 + x[1])*x[0]*(1 + x[0])*vertices[10] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[8] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[12] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[0] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[6] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[2] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[4];
+ rcol[1] = -0.5*(-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[15] + (-1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[17] + 0.25*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[1] - 0.5*(-1 + x[0])*x[1]*(-1 + x[1])*vertices[9] + 0.25*(-1 + x[0])*x[1]*(1 + x[1])*vertices[7] - 0.5*(-1 + x[0])*x[1]*(1 + x[1])*vertices[13] - 0.5*(1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[11] + (1 + x[0])*(1 + x[1])*(-1 + x[1])*vertices[17] + 0.25*(1 + x[0])*x[1]*(-1 + x[1])*vertices[3] - 0.5*(1 + x[0])*x[1]*(-1 + x[1])*vertices[9] + 0.25*(1 + x[0])*x[1]*(1 + x[1])*vertices[5] - 0.5*(1 + x[0])*x[1]*(1 + x[1])*vertices[13] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[11] - 0.5*x[0]*(1 + x[1])*(-1 + x[1])*vertices[15] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[1] + 0.25*x[0]*x[1]*(-1 + x[1])*vertices[3] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[5] + 0.25*x[0]*x[1]*(1 + x[1])*vertices[7];
+ scol[1] = -0.5*(-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[9] + (-1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[17] + 0.25*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[1] - 0.5*(-1 + x[1])*x[0]*(-1 + x[0])*vertices[15] + 0.25*(-1 + x[1])*x[0]*(1 + x[0])*vertices[3] - 0.5*(-1 + x[1])*x[0]*(1 + x[0])*vertices[11] - 0.5*(1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[13] + (1 + x[1])*(1 + x[0])*(-1 + x[0])*vertices[17] + 0.25*(1 + x[1])*x[0]*(-1 + x[0])*vertices[7] - 0.5*(1 + x[1])*x[0]*(-1 + x[0])*vertices[15] + 0.25*(1 + x[1])*x[0]*(1 + x[0])*vertices[5] - 0.5*(1 + x[1])*x[0]*(1 + x[0])*vertices[11] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[9] - 0.5*x[1]*(1 + x[0])*(-1 + x[0])*vertices[13] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[1] + 0.25*x[1]*x[0]*(-1 + x[0])*vertices[7] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[3] + 0.25*x[1]*x[0]*(1 + x[0])*vertices[5];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
cdef void Tet2Function3D(double* fx,
double* x,
double* vertices,
diff -r 398bded61777ec01e1bdf868edb966b7e52ece57 -r eaf22dd8f960f1767b66025127fb4413d0654853 yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -237,6 +237,22 @@
cdef int check_inside(self, double* mapped_coord) nogil
+
+cdef class Q2Sampler2D(NonlinearSolveSampler2D):
+
+ cdef void map_real_to_unit(self,
+ double* mapped_x,
+ double* vertices,
+ double* physical_x) nogil
+
+
+ cdef double sample_at_unit_point(self,
+ double* coord,
+ double* vals) nogil
+
+ cdef int check_inside(self, double* mapped_coord) nogil
+
+
cdef class T2Sampler2D(NonlinearSolveSampler2D):
cdef void map_real_to_unit(self,
diff -r 398bded61777ec01e1bdf868edb966b7e52ece57 -r eaf22dd8f960f1767b66025127fb4413d0654853 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -24,6 +24,8 @@
Q1Jacobian3D, \
Q1Function2D, \
Q1Jacobian2D, \
+ Q2Function2D, \
+ Q2Jacobian2D, \
W1Function3D, \
W1Jacobian3D, \
T2Function2D, \
@@ -340,15 +342,41 @@
double* mapped_x,
double* vertices,
double* physical_x) nogil:
+ '''
+
+ A thorough description of Newton's method and modifications for global
+ convergence can be found in Dennis's text "Numerical Methods for
+ Unconstrained Optimization and Nonlinear Equations."
+
+ x: solution vector; holds unit/mapped coordinates
+ xk: temporary vector for holding solution of current iteration
+ f: residual vector
+ r, s, t: three columns of Jacobian matrix corresponding to unit/mapped
+ coordinates r, s, and t
+ d: Jacobian determinant
+ s_n: Newton step vector
+ lam: fraction of Newton step by which to change x
+ alpha: constant proportional to how much residual required to decrease.
+ 1e-4 is value of alpha recommended by Dennis
+ err_c: Error of current iteration
+ err_plus: Error of next iteration
+ min_lam: minimum fraction of Newton step that the line search is allowed
+ to take. General experience suggests that lambda values smaller
+ than 1e-3 will not significantly reduce the residual, but we
+ set to 1e-6 just to be safe
+ '''
+
cdef int i
- cdef double d
+ cdef double d, lam
cdef double[3] f
cdef double[3] r
cdef double[3] s
cdef double[3] t
- cdef double[3] x
+ cdef double[3] x, xk, s_n
cdef int iterations = 0
- cdef double err
+ cdef double err_c, err_plus
+ cdef double alpha = 1e-4
+ cdef double min_lam = 1e-6
# initial guess
for i in range(3):
@@ -356,20 +384,38 @@
# initial error norm
self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 3)
+ err_c = maxnorm(f, 3)
# begin Newton iteration
- while (err > self.tolerance and iterations < self.max_iter):
+ while (err_c > self.tolerance and iterations < self.max_iter):
self.jac(r, s, t, x, vertices, physical_x)
d = determinant_3x3(r, s, t)
- x[0] = x[0] - (determinant_3x3(f, s, t)/d)
- x[1] = x[1] - (determinant_3x3(r, f, t)/d)
- x[2] = x[2] - (determinant_3x3(r, s, f)/d)
- self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 3)
+
+ s_n[0] = - (determinant_3x3(f, s, t)/d)
+ s_n[1] = - (determinant_3x3(r, f, t)/d)
+ s_n[2] = - (determinant_3x3(r, s, f)/d)
+ xk[0] = x[0] + s_n[0]
+ xk[1] = x[1] + s_n[1]
+ xk[2] = x[2] + s_n[2]
+ self.func(f, xk, vertices, physical_x)
+ err_plus = maxnorm(f, 3)
+
+ lam = 1
+ while err_plus > err_c * (1. - alpha * lam) and lam > min_lam:
+ lam = lam / 2
+ xk[0] = x[0] + lam * s_n[0]
+ xk[1] = x[1] + lam * s_n[1]
+ xk[2] = x[2] + lam * s_n[2]
+ self.func(f, xk, vertices, physical_x)
+ err_plus = maxnorm(f, 3)
+
+ x[0] = xk[0]
+ x[1] = xk[1]
+ x[2] = xk[2]
+ err_c = err_plus
iterations += 1
- if (err > self.tolerance):
+ if (err_c > self.tolerance):
# we did not converge, set bogus value
for i in range(3):
mapped_x[i] = -99.0
@@ -753,13 +799,38 @@
double* mapped_x,
double* vertices,
double* physical_x) nogil:
+ '''
+
+ A thorough description of Newton's method and modifications for global
+ convergence can be found in Dennis's text "Numerical Methods for
+ Unconstrained Optimization and Nonlinear Equations."
+
+ x: solution vector; holds unit/mapped coordinates
+ xk: temporary vector for holding solution of current iteration
+ f: residual vector
+ A: Jacobian matrix (derivative of residual vector wrt x)
+ d: Jacobian determinant
+ s_n: Newton step vector
+ lam: fraction of Newton step by which to change x
+ alpha: constant proportional to how much residual required to decrease.
+ 1e-4 is value of alpha recommended by Dennis
+ err_c: Error of current iteration
+ err_plus: Error of next iteration
+ min_lam: minimum fraction of Newton step that the line search is allowed
+ to take. General experience suggests that lambda values smaller
+ than 1e-3 will not significantly reduce the residual, but we
+ set to 1e-6 just to be safe
+ '''
+
cdef int i
- cdef double d
+ cdef double d, lam
cdef double[2] f
- cdef double[2] x
+ cdef double[2] x, xk, s_n
cdef double[4] A
cdef int iterations = 0
- cdef double err
+ cdef double err_c, err_plus
+ cdef double alpha = 1e-4
+ cdef double min_lam = 1e-6
# initial guess
for i in range(2):
@@ -767,21 +838,34 @@
# initial error norm
self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 2)
+ err_c = maxnorm(f, 2)
# begin Newton iteration
- while (err > self.tolerance and iterations < self.max_iter):
+ while (err_c > self.tolerance and iterations < self.max_iter):
self.jac(&A[0], &A[2], x, vertices, physical_x)
d = (A[0]*A[3] - A[1]*A[2])
- x[0] -= ( A[3]*f[0] - A[2]*f[1]) / d
- x[1] -= (-A[1]*f[0] + A[0]*f[1]) / d
+ s_n[0] = -( A[3]*f[0] - A[2]*f[1]) / d
+ s_n[1] = -(-A[1]*f[0] + A[0]*f[1]) / d
+ xk[0] = x[0] + s_n[0]
+ xk[1] = x[1] + s_n[1]
+ self.func(f, xk, vertices, physical_x)
+ err_plus = maxnorm(f, 2)
+
+ lam = 1
+ while err_plus > err_c * (1. - alpha * lam) and lam > min_lam:
+ lam = lam / 2
+ xk[0] = x[0] + lam * s_n[0]
+ xk[1] = x[1] + lam * s_n[1]
+ self.func(f, xk, vertices, physical_x)
+ err_plus = maxnorm(f, 2)
- self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 2)
+ x[0] = xk[0]
+ x[1] = xk[1]
+ err_c = err_plus
iterations += 1
- if (err > self.tolerance):
+ if (err_c > self.tolerance):
# we did not converge, set bogus value
for i in range(2):
mapped_x[i] = -99.0
@@ -830,6 +914,63 @@
return 0
return 1
+cdef class Q2Sampler2D(NonlinearSolveSampler2D):
+
+ '''
+
+ This implements sampling inside a 2D, quadratic, quadrilateral mesh element.
+
+ '''
+
+ def __init__(self):
+ super(Q2Sampler2D, self).__init__()
+ self.num_mapped_coords = 2
+ self.dim = 2
+ self.func = Q2Function2D
+ self.jac = Q2Jacobian2D
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
+ cdef double F, rm, rp, sm, sp
+ cdef double[9] phi
+ cdef double rv = 0
+
+ zet = coord[0]
+ eta = coord[1]
+ zetm = coord[0] - 1.
+ zetp = coord[0] + 1.
+ etam = coord[1] - 1.
+ etap = coord[1] + 1.
+
+ phi[0] = zet * zetm * eta * etam / 4.
+ phi[1] = zet * zetp * eta * etam / 4.
+ phi[2] = zet * zetp * eta * etap / 4.
+ phi[3] = zet * zetm * eta * etap / 4.
+ phi[4] = zetp * zetm * eta * etam / -2.
+ phi[5] = zet * zetp * etap * etam / -2.
+ phi[6] = zetp * zetm * eta * etap / -2.
+ phi[7] = zet * zetm * etap * etam / -2.
+ phi[8] = zetp * zetm * etap * etam
+
+ for i in range(9):
+ rv += vals[i] * phi[i]
+
+ return rv
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_inside(self, double* mapped_coord) nogil:
+ # for quads, we check whether the mapped_coord is between
+ # -1 and 1 in both directions.
+ if (fabs(mapped_coord[0]) - 1.0 > self.inclusion_tol or
+ fabs(mapped_coord[1]) - 1.0 > self.inclusion_tol):
+ return 0
+ return 1
+
+
cdef class T2Sampler2D(NonlinearSolveSampler2D):
'''
@@ -1093,6 +1234,23 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+def test_quad2_sampler(np.ndarray[np.float64_t, ndim=2] vertices,
+ np.ndarray[np.float64_t, ndim=1] field_values,
+ np.ndarray[np.float64_t, ndim=1] physical_x):
+
+ cdef double val
+
+ sampler = Q2Sampler2D()
+
+ val = sampler.sample_at_real_point(<double*> vertices.data,
+ <double*> field_values.data,
+ <double*> physical_x.data)
+
+ return val
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
def test_tri2_sampler(np.ndarray[np.float64_t, ndim=2] vertices,
np.ndarray[np.float64_t, ndim=1] field_values,
np.ndarray[np.float64_t, ndim=1] physical_x):
diff -r 398bded61777ec01e1bdf868edb966b7e52ece57 -r eaf22dd8f960f1767b66025127fb4413d0654853 yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -30,6 +30,7 @@
P1Sampler3D, \
Q1Sampler3D, \
Q1Sampler2D, \
+ Q2Sampler2D, \
S2Sampler3D, \
W1Sampler3D, \
T2Sampler2D, \
@@ -609,6 +610,8 @@
sampler = P1Sampler1D()
elif ndim == 2 and nvertices == 4:
sampler = Q1Sampler2D()
+ elif ndim == 2 and nvertices == 9:
+ sampler = Q2Sampler2D()
elif ndim == 2 and nvertices == 6:
sampler = T2Sampler2D()
elif ndim == 3 and nvertices == 10:
diff -r 398bded61777ec01e1bdf868edb966b7e52ece57 -r eaf22dd8f960f1767b66025127fb4413d0654853 yt/utilities/lib/tests/test_element_mappings.py
--- a/yt/utilities/lib/tests/test_element_mappings.py
+++ b/yt/utilities/lib/tests/test_element_mappings.py
@@ -21,6 +21,7 @@
test_hex_sampler, \
test_tri_sampler, \
test_quad_sampler, \
+ test_quad2_sampler, \
test_hex20_sampler, \
test_wedge_sampler, \
test_tri2_sampler, \
@@ -77,6 +78,23 @@
check_all_vertices(test_quad_sampler, vertices, field_values)
+def test_Q2Sampler2D():
+
+ vertices = np.array([[2., 3.],
+ [7., 4.],
+ [10., 15.],
+ [4., 12.],
+ [4.5, 3.5],
+ [8.5, 9.5],
+ [7., 13.5],
+ [3., 7.5],
+ [5.75, 8.5]])
+
+ field_values = np.array([7., 27., 40., 12., 13., 30., 22., 9., 16.])
+
+ check_all_vertices(test_quad2_sampler, vertices, field_values)
+
+
def test_Q1Sampler3D():
vertices = np.array([[2.00657905, 0.6888599, 1.4375],
[1.8658198, 1.00973171, 1.4375],
diff -r 398bded61777ec01e1bdf868edb966b7e52ece57 -r eaf22dd8f960f1767b66025127fb4413d0654853 yt/utilities/mesh_types.yaml
--- a/yt/utilities/mesh_types.yaml
+++ b/yt/utilities/mesh_types.yaml
@@ -24,6 +24,22 @@
(1 + x[0])*(1 + x[1])/4.,
(1 - x[0])*(1 + x[1])/4.]
+Quad9:
+ mesh_type: Q2
+ num_dim: 2
+ num_vertices: 9
+ num_mapped_coords: 2
+ shape_functions: |
+ [x[0] * (x[0] - 1) * x[1] * (x[1] - 1) / 4.,
+ x[0] * (x[0] + 1) * x[1] * (x[1] - 1) / 4.,
+ x[0] * (x[0] + 1) * x[1] * (x[1] + 1) / 4.,
+ x[0] * (x[0] - 1) * x[1] * (x[1] + 1) / 4.,
+ (x[0] + 1) * (x[0] - 1) * x[1] * (x[1] - 1) / -2.,
+ x[0] * (x[0] + 1) * (x[1] + 1) * (x[1] - 1) / -2.,
+ (x[0] + 1) * (x[0] - 1) * x[1] * (x[1] + 1) / -2.,
+ x[0] * (x[0] - 1) * (x[1] + 1) * (x[1] - 1) / -2.,
+ (x[0] + 1) * (x[0] - 1) * (x[1] + 1) * (x[1] - 1)]
+
Wedge6:
mesh_type: W1
num_dim: 3
diff -r 398bded61777ec01e1bdf868edb966b7e52ece57 -r eaf22dd8f960f1767b66025127fb4413d0654853 yt/visualization/tests/test_mesh_slices.py
--- a/yt/visualization/tests/test_mesh_slices.py
+++ b/yt/visualization/tests/test_mesh_slices.py
@@ -49,6 +49,15 @@
for field in ds.field_list:
yield compare(ds, field, "answers_tri2_%s_%s" % (field[0], field[1]))
+quad2 = "SecondOrderQuads/lid_driven_out.e"
+
+ at requires_ds(quad2)
+def test_quad2():
+ ds = data_dir_load(quad2, kwargs={'step':-1})
+ field_list = [('all', 'T'), ('all', 'vel_x'), ('all', 'vel_y')]
+ for field in field_list:
+ yield compare(ds, field, "answers_quad2_%s_%s" % (field[0], field[1]))
+
multi_region = "MultiRegion/two_region_example_out.e"
@requires_ds(multi_region)
diff -r 398bded61777ec01e1bdf868edb966b7e52ece57 -r eaf22dd8f960f1767b66025127fb4413d0654853 yt/visualization/volume_rendering/tests/test_mesh_render.py
--- a/yt/visualization/volume_rendering/tests/test_mesh_render.py
+++ b/yt/visualization/volume_rendering/tests/test_mesh_render.py
@@ -164,6 +164,8 @@
ytcfg["yt", "ray_tracing_engine"] = engine
ds = data_dir_load(tet10, kwargs={'step':-1})
sc = create_scene(ds, field)
+ ms = sc.get_source(0)
+ ms.color_bounds = (-.01, .2)
im = sc.render()
return compare(ds, im, "%s_render_answers_tet10_%s_%s" %
(engine, field[0], field[1]))
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list